将遵循一个方向的点集组合在一起的算法
Algorithm to group sets of points together that follow a direction
注意：我把这个问题放在 MATLAB 和 Python 标签中，因为我最精通这些语言。但是，我欢迎任何语言的解决方案。
问题序言
我用鱼眼镜头拍了张照片。此图像由带有一堆方形对象的图案组成。我想对这张图像做的是检测每个正方形的质心，然后使用这些点对图像进行去畸变——具体来说，我正在寻找正确的畸变模型参数。需要注意的是，并不是所有的方块都需要检测。只要他们中的大多数是，那就完全没问题了……但这不是本文的重点 post。参数估计算法我已经写过了，但问题是它需要在图像中出现共线的点。
我想问的基本问题是给定这些点，将它们组合在一起的最佳方法是什么，以便每组由水平线或垂直线组成？
我的问题背景
这对于我要问的问题来说并不重要，但如果您想知道我从哪里获得数据并进一步了解我要问的问题，请阅读。如果您不感兴趣，那么您可以直接跳到下面的问题设置部分。
我正在处理的图像示例如下所示：
这是一张 960 x 960 的图片。该图像原本分辨率较高，但我对图像进行二次采样以加快处理时间。如您所见，图像中散布着一堆正方形图案。另外，我计算的质心是关于上述二次采样图像的。
我设置的检索质心的管道如下：
- 执行 Canny 边缘检测
- 关注可最大程度减少误报的感兴趣区域。这个感兴趣的区域基本上是没有任何覆盖其一侧的黑色胶带的正方形。
- 找到所有不同的封闭轮廓
对于每个不同的封闭轮廓...
一个。执行哈里斯角点检测
b。判断结果是否有4个角点
c。如果是这样，那么这个轮廓属于一个正方形并找到这个形状的质心
d。如果没有，则跳过此形状
将步骤 #4 中检测到的所有质心放入矩阵中以供进一步检查。
这是上图的示例结果。每个检测到的方块都有四个点，根据它相对于方块本身的位置进行颜色编码。对于我检测到的每个质心，我都会在图像本身的质心所在的位置写一个 ID。
通过上图，检测到37个方块。
问题设置
假设我有一些图像像素点存储在 N x 3 矩阵中。前两列是 x（水平）和 y（垂直）坐标，其中在图像坐标 space 中，y 坐标是 inverted，表示正数y向下移动。第三列是与该点关联的 ID。
这里是一些用 MATLAB 编写的代码，它获取这些点，将它们绘制到二维网格上，并用矩阵的第三列标记每个点。如果您阅读以上背景，这些是我上面概述的算法检测到的点。
data = [ 475. , 605.75, 1.;
571. , 586.5 , 2.;
233. , 558.5 , 3.;
669.5 , 562.75, 4.;
291.25, 546.25, 5.;
759. , 536.25, 6.;
362.5 , 531.5 , 7.;
448. , 513.5 , 8.;
834.5 , 510. , 9.;
897.25, 486. , 10.;
545.5 , 491.25, 11.;
214.5 , 481.25, 12.;
271.25, 463. , 13.;
646.5 , 466.75, 14.;
739. , 442.75, 15.;
340.5 , 441.5 , 16.;
817.75, 421.5 , 17.;
423.75, 417.75, 18.;
202.5 , 406. , 19.;
519.25, 392.25, 20.;
257.5 , 382. , 21.;
619.25, 368.5 , 22.;
148. , 359.75, 23.;
324.5 , 356. , 24.;
713. , 347.75, 25.;
195. , 335. , 26.;
793.5 , 332.5 , 27.;
403.75, 328. , 28.;
249.25, 308. , 29.;
495.5 , 300.75, 30.;
314. , 279. , 31.;
764.25, 249.5 , 32.;
389.5 , 249.5 , 33.;
475. , 221.5 , 34.;
565.75, 199. , 35.;
802.75, 173.75, 36.;
733. , 176.25, 37.];
figure; hold on;
axis ij;
scatter(data(:,1), data(:,2),40, 'r.');
text(data(:,1)+10, data(:,2)+10, num2str(data(:,3)));
与 Python 类似，使用 numpy 和 matplotlib，我们有：
import numpy as np
import matplotlib.pyplot as plt
data = np.array([[ 475. , 605.75, 1. ],
[ 571. , 586.5 , 2. ],
[ 233. , 558.5 , 3. ],
[ 669.5 , 562.75, 4. ],
[ 291.25, 546.25, 5. ],
[ 759. , 536.25, 6. ],
[ 362.5 , 531.5 , 7. ],
[ 448. , 513.5 , 8. ],
[ 834.5 , 510. , 9. ],
[ 897.25, 486. , 10. ],
[ 545.5 , 491.25, 11. ],
[ 214.5 , 481.25, 12. ],
[ 271.25, 463. , 13. ],
[ 646.5 , 466.75, 14. ],
[ 739. , 442.75, 15. ],
[ 340.5 , 441.5 , 16. ],
[ 817.75, 421.5 , 17. ],
[ 423.75, 417.75, 18. ],
[ 202.5 , 406. , 19. ],
[ 519.25, 392.25, 20. ],
[ 257.5 , 382. , 21. ],
[ 619.25, 368.5 , 22. ],
[ 148. , 359.75, 23. ],
[ 324.5 , 356. , 24. ],
[ 713. , 347.75, 25. ],
[ 195. , 335. , 26. ],
[ 793.5 , 332.5 , 27. ],
[ 403.75, 328. , 28. ],
[ 249.25, 308. , 29. ],
[ 495.5 , 300.75, 30. ],
[ 314. , 279. , 31. ],
[ 764.25, 249.5 , 32. ],
[ 389.5 , 249.5 , 33. ],
[ 475. , 221.5 , 34. ],
[ 565.75, 199. , 35. ],
[ 802.75, 173.75, 36. ],
[ 733. , 176.25, 37. ]])
plt.figure()
plt.gca().invert_yaxis()
plt.plot(data[:,0], data[:,1], 'r.', markersize=14)
for idx in np.arange(data.shape[0]):
plt.text(data[idx,0]+10, data[idx,1]+10, str(int(data[idx,2])), size='large')
plt.show()
我们得到：
回到问题
如您所见，这些点或多或少呈网格状，您可以看到我们可以在这些点之间形成线条。具体可以看到有横向和纵向可以组成的线
比如你参考我问题背景部分的图片，我们可以看到有5组点可以水平分组。例如，点23、26、29、31、33、34、35、37、36为一组。第 19、21、24、28、30 和 32 点构成另一组，依此类推。同样在垂直方向上，我们可以看到26、19、12、3点为一组，29、21、13、5点为一组，依此类推。
要问的问题
我的问题是：假设点可以在任何方向，有什么方法可以分别成功地将水平分组和垂直分组中的点分组？
条件
每行必须至少三个点。如果小于此，则这不符合细分市场的条件。因此，点 36 和 10 不符合垂直线的条件，同样孤立点 23 不应该符合垂直线的条件，但它是第一个水平分组的一部分。
以上校准图案可以任意方向。然而，对于我正在处理的情况，你能得到的最糟糕的方向是你在上面的背景部分看到的。
预期输出
输出将是一对列表，其中第一个列表包含元素，其中每个元素为您提供一系列点 ID，这些点 ID 形成一条水平线。同样，第二个列表包含元素，其中每个元素为您提供一系列点 ID，这些点 ID 形成一条垂直线。
因此，水平序列的预期输出如下所示：
MATLAB
horiz_list = {[23, 26, 29, 31, 33, 34, 35, 37, 36], [19, 21, 24, 28, 30, 32], ...};
vert_list = {[26, 19, 12, 3], [29, 21, 13, 5], ....};
Python
horiz_list = [[23, 26, 29, 31, 33, 34, 35, 37, 36], [19, 21, 24, 28, 30, 32], ....]
vert_list = [[26, 19, 12, 3], [29, 21, 13, 5], ...]
我试过的
从算法上讲，我尝试的是撤消在这些点上经历的旋转。我已经执行 Principal Components Analysis 并尝试相对于计算的正交基向量投影点，以便这些点或多或少位于直矩形网格上。
一旦我有了它，做一些扫描线处理就是一件简单的事情，你可以根据水平或垂直坐标上的差异变化对点进行分组。您可以按 x 或 y 值对坐标进行排序，然后检查这些排序后的坐标并寻找较大的变化。一旦遇到这种变化，就可以将变化之间的点组合在一起以形成线条。对每个维度执行此操作将为您提供水平或垂直分组。
关于 PCA，这是我在 MATLAB 和 Python:
中所做的
MATLAB
%# Step #1 - Get just the data - no IDs
data_raw = data(:,1:2);
%# Decentralize mean
data_nomean = bsxfun(@minus, data_raw, mean(data_raw,1));
%# Step #2 - Determine covariance matrix
%# This already decentralizes the mean
cov_data = cov(data_raw);
%# Step #3 - Determine right singular vectors
[~,~,V] = svd(cov_data);
%# Step #4 - Transform data with respect to basis
F = V.'*data_nomean.';
%# Visualize both the original data points and transformed data
figure;
plot(F(1,:), F(2,:), 'b.', 'MarkerSize', 14);
axis ij;
hold on;
plot(data(:,1), data(:,2), 'r.', 'MarkerSize', 14);
Python
import numpy as np
import numpy.linalg as la
# Step #1 and Step #2 - Decentralize mean
centroids_raw = data[:,:2]
mean_data = np.mean(centroids_raw, axis=0)
# Transpose for covariance calculation
data_nomean = (centroids_raw - mean_data).T
# Step #3 - Determine covariance matrix
# Doesn't matter if you do this on the decentralized result
# or the normal result - cov subtracts the mean off anyway
cov_data = np.cov(data_nomean)
# Step #4 - Determine right singular vectors via SVD
# Note - This is already V^T, so there's no need to transpose
_,_,V = la.svd(cov_data)
# Step #5 - Transform data with respect to basis
data_transform = np.dot(V, data_nomean).T
plt.figure()
plt.gca().invert_yaxis()
plt.plot(data[:,0], data[:,1], 'b.', markersize=14)
plt.plot(data_transform[:,0], data_transform[:,1], 'r.', markersize=14)
plt.show()
以上代码不仅重新投影数据，而且还将原始点和投影点一起绘制在一个图中。然而，当我尝试重新投影我的数据时，这是我得到的图：
红色的点是原始图像坐标，而蓝色的点被重新投影到基向量上以尝试消除旋转。它仍然不能完全完成工作。关于这些点仍然有一些方向，所以如果我尝试执行我的扫描线算法，来自下面的线的点用于水平跟踪或到侧面的点用于垂直跟踪将无意中分组，这是不正确的。
也许我想多了这个问题，但如果您对此有任何见解，我们将不胜感激。如果答案确实很棒，我会倾向于奖励高额赏金，因为我已经在这个问题上停留了很长一段时间。
我希望这个问题不是冗长的。如果您不知道如何解决这个问题，那么我感谢您花时间阅读我的问题。
期待您的任何见解。非常感谢！
虽然我无法提出比您已经尝试过的更好的方法来对任何给定的质心点列表进行分组，但我希望以下想法可以帮助您：
由于您对图像的内容非常具体（包含一个正方形区域），我想知道您是否真的需要根据 problem setup 中给出的数据对质心点进行分组，或者是否您也可以使用 Background to the problem 中描述的数据。由于您已经确定了每个检测到的正方形的角以及它们在给定正方形中的位置，因此在我看来，通过比较角坐标来确定给定正方形的邻居是非常准确的。
因此，为了找到任何正方形右邻的候选者，我建议您将该正方形的右上角和右下角与任何其他正方形（或其中的任何正方形）的左上角和左下角进行比较一定距离）。仅允许较小的垂直差异和稍大的水平差异，您可以 "match" 两个正方形，如果 两者 它们对应的角点足够靠近。
通过使用角之间允许的 vertical/horizontal 差异的上限，您甚至可以将这些边界内的最佳匹配方块指定为邻居
一个问题可能是你没有检测到所有的方块，所以在方块 30 和 32 之间有一个相当大的 space。因为你说你需要 'at least' 每行 3 个方块，您可以简单地忽略该水平线上的方块 32。如果这不适合您，您可以尝试匹配尽可能多的方块，然后使用先前计算的数据将 "missing" 方块分配给网格中的一个点：
在关于方块 32 的示例中，您会检测到它有上下邻居 27 和 37。此外，您应该已经能够确定方块 27 位于第 1 行内，而 37 位于第 3 行内，所以你可以将方块 32 分配给中间的 "best matching" 行，在这种情况下显然是 2。
这种通用方法基本上是您已经尝试过的方法，但应该更准确，因为您现在比较的是两条线的方向和距离，而不是简单地比较网格中两点的位置。
也作为您之前尝试的旁路节点 - 您可以使用黑色角线稍微校正图像的初始旋转吗？这可能会使进一步的失真算法（比如您在评论中与 knedlsepp 讨论的那些算法）更加准确。 （编辑：我刚才确实阅读了Parag的评论-通过线的角度比较点当然与事先旋转图像基本相同）
注意 1：它有许多设置 -> 对于其他图像可能需要更改这些设置才能获得您想要的结果请参阅 % 设置 - 使用这些值
注意 2：它不会找到您想要的所有行 -> 但它是一个起点....
要调用此函数，请在命令提示符中调用此函数：
>> [h, v] = testLines;
我们得到：
>> celldisp(h)
h{1} =
1 2 4 6 9 10
h{2} =
3 5 7 8 11 14 15 17
h{3} =
1 2 4 6 9 10
h{4} =
3 5 7 8 11 14 15 17
h{5} =
1 2 4 6 9 10
h{6} =
3 5 7 8 11 14 15 17
h{7} =
3 5 7 8 11 14 15 17
h{8} =
1 2 4 6 9 10
h{9} =
1 2 4 6 9 10
h{10} =
12 13 16 18 20 22 25 27
h{11} =
13 16 18 20 22 25 27
h{12} =
3 5 7 8 11 14 15 17
h{13} =
3 5 7 8 11 14 15
h{14} =
12 13 16 18 20 22 25 27
h{15} =
3 5 7 8 11 14 15 17
h{16} =
12 13 16 18 20 22 25 27
h{17} =
19 21 24 28 30
h{18} =
21 24 28 30
h{19} =
12 13 16 18 20 22 25 27
h{20} =
19 21 24 28 30
h{21} =
12 13 16 18 20 22 24 25
h{22} =
12 13 16 18 20 22 24 25 27
h{23} =
23 26 29 31 33 34 35
h{24} =
23 26 29 31 33 34 35 37
h{25} =
23 26 29 31 33 34 35 36 37
h{26} =
33 34 35 37 36
h{27} =
31 33 34 35 37
>> celldisp(v)
v{1} =
33 28 18 8 1
v{2} =
34 30 20 11 2
v{3} =
26 19 12 3
v{4} =
35 22 14 4
v{5} =
29 21 13 5
v{6} =
25 15 6
v{7} =
31 24 16 7
v{8} =
37 32 27 17 9
还会生成一个图形，通过每个适当的点集绘制线：
function [horiz_list, vert_list] = testLines
global counter;
global colours;
close all;
data = [ 475. , 605.75, 1.;
571. , 586.5 , 2.;
233. , 558.5 , 3.;
669.5 , 562.75, 4.;
291.25, 546.25, 5.;
759. , 536.25, 6.;
362.5 , 531.5 , 7.;
448. , 513.5 , 8.;
834.5 , 510. , 9.;
897.25, 486. , 10.;
545.5 , 491.25, 11.;
214.5 , 481.25, 12.;
271.25, 463. , 13.;
646.5 , 466.75, 14.;
739. , 442.75, 15.;
340.5 , 441.5 , 16.;
817.75, 421.5 , 17.;
423.75, 417.75, 18.;
202.5 , 406. , 19.;
519.25, 392.25, 20.;
257.5 , 382. , 21.;
619.25, 368.5 , 22.;
148. , 359.75, 23.;
324.5 , 356. , 24.;
713. , 347.75, 25.;
195. , 335. , 26.;
793.5 , 332.5 , 27.;
403.75, 328. , 28.;
249.25, 308. , 29.;
495.5 , 300.75, 30.;
314. , 279. , 31.;
764.25, 249.5 , 32.;
389.5 , 249.5 , 33.;
475. , 221.5 , 34.;
565.75, 199. , 35.;
802.75, 173.75, 36.;
733. , 176.25, 37.];
figure; hold on;
axis ij;
% Change due to Benoit_11
scatter(data(:,1), data(:,2),40, 'r.'); text(data(:,1)+10, data(:,2)+10, num2str(data(:,3)));
text(data(:,1)+10, data(:,2)+10, num2str(data(:,3)));
% Process your data as above then run the function below(note it has sub functions)
counter = 0;
colours = 'bgrcmy';
[horiz_list, vert_list] = findClosestPoints ( data(:,1), data(:,2) );
function [horiz_list, vert_list] = findClosestPoints ( x, y )
% calc length of points
nX = length(x);
% set up place holder flags
modelledH = false(nX,1);
modelledV = false(nX,1);
horiz_list = {};
vert_list = {};
% loop for all points
for p=1:nX
% have we already modelled a horizontal line through these?
% second last param - true - horizontal, false - vertical
if modelledH(p)==false
[modelledH, index] = ModelPoints ( p, x, y, modelledH, true, true );
horiz_list = [horiz_list index];
else
[~, index] = ModelPoints ( p, x, y, modelledH, true, false );
horiz_list = [horiz_list index];
end
% make a temp copy of the x and y and remove any of the points modelled
% from the horizontal -> this is to avoid them being found in the
% second call.
tempX = x;
tempY = y;
tempX(index) = NaN;
tempY(index) = NaN;
tempX(p) = x(p);
tempY(p) = y(p);
% Have we found a vertial line?
if modelledV(p)==false
[modelledV, index] = ModelPoints ( p, tempX, tempY, modelledV, false, true );
vert_list = [vert_list index];
end
end
end
function [modelled, index] = ModelPoints ( p, x, y, modelled, method, fullRun )
% p - row in your original data matrix
% x - data(:,1)
% y - data(:,2)
% modelled - array of flags to whether rows have been modelled
% method - horizontal or vertical (used to calc graadients)
% fullRun - full calc or just to get indexes
% this could be made better by storing the indexes of each horizontal in the method above
% Settings - play around with these values
gradDelta = 0.2; % find points where gradient is less than this value
gradLimit = 0.45; % if mean gradient of line is above this ignore
numberOfPointsToCheck = 7; % number of points to check when look along the line
% to find other points (this reduces chance of it
% finding other points far away
% I optimised this for your example to be 7
% Try varying it and you will see how it effect the result.
% Find the index of points which are inline.
[index, grad] = CalcIndex ( x, y, p, gradDelta, method );
% check gradient of line
if abs(mean(grad))>gradLimit
index = [];
return
end
% add point of interest to index
index = [p index];
% loop through all points found above to find any other points which are in
% line with these points (this allows for slight curvature
combineIndex = [];
for ii=2:length(index)
% Find inedex of the points found above (find points on curve)
[index2] = CalcIndex ( x, y, index(ii), gradDelta, method, numberOfPointsToCheck, grad(ii-1) );
% Check that the point on this line are on the original (i.e. inline -> not at large angle
if any(ismember(index,index2))
% store points found
combineIndex = unique([index2 combineIndex]);
end
end
% copy to index
index = combineIndex;
if fullRun
% do some plotting
% TODO: here you would need to calculate your arrays to output.
xx = x(index);
[sX,sOrder] = sort(xx);
% Check its found at least 3 points
if length ( index(sOrder) ) > 2
% flag the modelled on the points found
modelled(index(sOrder)) = true;
% plot the data
plot ( x(index(sOrder)), y(index(sOrder)), colours(mod(counter,numel(colours)) + 1));
counter = counter + 1;
end
index = index(sOrder);
end
end
function [index, gradCheck] = CalcIndex ( x, y, p, gradLimit, method, nPoints2Consider, refGrad )
% x - data(:,1)
% y - data(:,2)
% p - point of interest
% method (x/y) or (y\x)
% nPoints2Consider - only look at N points (options)
% refgrad - rather than looking for gradient of closest point -> use this
% - reference gradient to find similar points (finds points on curve)
nX = length(x);
% calculate gradient
for g=1:nX
if method
grad(g) = (x(g)-x(p))\(y(g)-y(p));
else
grad(g) = (y(g)-y(p))\(x(g)-x(p));
end
end
% find distance to all other points
delta = sqrt ( (x-x(p)).^2 + (y-y(p)).^2 );
% set its self = NaN
delta(delta==min(delta)) = NaN;
% find the closest points
[m,order] = sort(delta);
if nargin == 7
% for finding along curve
% set any far away points to be NaN
grad(order(nPoints2Consider+1:end)) = NaN;
% find the closest points to the reference gradient within the allowable limit
index = find(abs(grad-refGrad)<gradLimit==1);
% store output
gradCheck = grad(index);
else
% find the points which are closes to the gradient of the closest point
index = find(abs(grad-grad(order(1)))<gradLimit==1);
% store gradients to output
gradCheck = grad(index);
end
end
end
我正在使用已发布图像的裁剪版本作为输入。这里我只讨论网格的方向可以被认为接近 horizontal/vertical 的情况。这可能无法完全解决您的范围，但我认为它可能会给您一些指导。
将图像二值化，以填充扭曲的方块。这里我使用了一个简单的 Otsu 阈值。然后对这个二值图像进行距离变换。
在距离变换图像中，我们看到正方形之间的间隙为峰值。
要获得水平方向的线，取距离图像的每个 列 的局部最大值，然后找到连通分量。
要获得垂直方向的线，取距离图像的每个 行 的局部最大值，然后找到连通分量。
下图显示了这样找到的水平线和垂直线，角点为圆圈。
对于相当长的连通分量，可以拟合一条曲线（直线或多项式），然后对角点进行分类，比如根据到曲线的距离、该点在曲线的哪一侧等.
我是用 Matlab 做的。我没有尝试曲线拟合和分类部分。
clear all;
close all;
im = imread('0RqUd-1.jpg');
gr = rgb2gray(im);
% any preprocessing to get a binary image that fills the distorted squares
bw = ~im2bw(gr, graythresh(gr));
di = bwdist(bw); % distance transform
di2 = imdilate(di, ones(3)); % propagate max
corners = corner(gr); % simple corners
% find regional max for each column of dist image
regmxh = zeros(size(di2));
for c = 1:size(di2, 2)
regmxh(:, c) = imregionalmax(di2(:, c));
end
% label connected components
ccomph = bwlabel(regmxh, 8);
% find regional max for each row of dist image
regmxv = zeros(size(di2));
for r = 1:size(di2, 1)
regmxv(r, :) = imregionalmax(di2(r, :));
end
% label connected components
ccompv = bwlabel(regmxv, 8);
figure, imshow(gr, [])
hold on
plot(corners(:, 1), corners(:, 2), 'ro')
figure, imshow(di, [])
figure, imshow(label2rgb(ccomph), [])
hold on
plot(corners(:, 1), corners(:, 2), 'ro')
figure, imshow(label2rgb(ccompv), [])
hold on
plot(corners(:, 1), corners(:, 2), 'ro')
要获得任意方向网格的这些线，您可以将距离图像视为图形并找到最佳路径。请参阅 here 了解基于图形的不错方法。
注意：我把这个问题放在 MATLAB 和 Python 标签中，因为我最精通这些语言。但是，我欢迎任何语言的解决方案。
问题序言
我用鱼眼镜头拍了张照片。此图像由带有一堆方形对象的图案组成。我想对这张图像做的是检测每个正方形的质心，然后使用这些点对图像进行去畸变——具体来说，我正在寻找正确的畸变模型参数。需要注意的是，并不是所有的方块都需要检测。只要他们中的大多数是，那就完全没问题了……但这不是本文的重点 post。参数估计算法我已经写过了，但问题是它需要在图像中出现共线的点。
我想问的基本问题是给定这些点，将它们组合在一起的最佳方法是什么，以便每组由水平线或垂直线组成？
我的问题背景
这对于我要问的问题来说并不重要，但如果您想知道我从哪里获得数据并进一步了解我要问的问题，请阅读。如果您不感兴趣，那么您可以直接跳到下面的问题设置部分。
我正在处理的图像示例如下所示：
这是一张 960 x 960 的图片。该图像原本分辨率较高，但我对图像进行二次采样以加快处理时间。如您所见，图像中散布着一堆正方形图案。另外，我计算的质心是关于上述二次采样图像的。
我设置的检索质心的管道如下：
- 执行 Canny 边缘检测
- 关注可最大程度减少误报的感兴趣区域。这个感兴趣的区域基本上是没有任何覆盖其一侧的黑色胶带的正方形。
- 找到所有不同的封闭轮廓
对于每个不同的封闭轮廓...
一个。执行哈里斯角点检测
b。判断结果是否有4个角点
c。如果是这样，那么这个轮廓属于一个正方形并找到这个形状的质心
d。如果没有，则跳过此形状
将步骤 #4 中检测到的所有质心放入矩阵中以供进一步检查。
这是上图的示例结果。每个检测到的方块都有四个点，根据它相对于方块本身的位置进行颜色编码。对于我检测到的每个质心，我都会在图像本身的质心所在的位置写一个 ID。
通过上图，检测到37个方块。
问题设置
假设我有一些图像像素点存储在 N x 3 矩阵中。前两列是 x（水平）和 y（垂直）坐标，其中在图像坐标 space 中，y 坐标是 inverted，表示正数y向下移动。第三列是与该点关联的 ID。
这里是一些用 MATLAB 编写的代码，它获取这些点，将它们绘制到二维网格上，并用矩阵的第三列标记每个点。如果您阅读以上背景，这些是我上面概述的算法检测到的点。
data = [ 475. , 605.75, 1.;
571. , 586.5 , 2.;
233. , 558.5 , 3.;
669.5 , 562.75, 4.;
291.25, 546.25, 5.;
759. , 536.25, 6.;
362.5 , 531.5 , 7.;
448. , 513.5 , 8.;
834.5 , 510. , 9.;
897.25, 486. , 10.;
545.5 , 491.25, 11.;
214.5 , 481.25, 12.;
271.25, 463. , 13.;
646.5 , 466.75, 14.;
739. , 442.75, 15.;
340.5 , 441.5 , 16.;
817.75, 421.5 , 17.;
423.75, 417.75, 18.;
202.5 , 406. , 19.;
519.25, 392.25, 20.;
257.5 , 382. , 21.;
619.25, 368.5 , 22.;
148. , 359.75, 23.;
324.5 , 356. , 24.;
713. , 347.75, 25.;
195. , 335. , 26.;
793.5 , 332.5 , 27.;
403.75, 328. , 28.;
249.25, 308. , 29.;
495.5 , 300.75, 30.;
314. , 279. , 31.;
764.25, 249.5 , 32.;
389.5 , 249.5 , 33.;
475. , 221.5 , 34.;
565.75, 199. , 35.;
802.75, 173.75, 36.;
733. , 176.25, 37.];
figure; hold on;
axis ij;
scatter(data(:,1), data(:,2),40, 'r.');
text(data(:,1)+10, data(:,2)+10, num2str(data(:,3)));
与 Python 类似，使用 numpy 和 matplotlib，我们有：
import numpy as np
import matplotlib.pyplot as plt
data = np.array([[ 475. , 605.75, 1. ],
[ 571. , 586.5 , 2. ],
[ 233. , 558.5 , 3. ],
[ 669.5 , 562.75, 4. ],
[ 291.25, 546.25, 5. ],
[ 759. , 536.25, 6. ],
[ 362.5 , 531.5 , 7. ],
[ 448. , 513.5 , 8. ],
[ 834.5 , 510. , 9. ],
[ 897.25, 486. , 10. ],
[ 545.5 , 491.25, 11. ],
[ 214.5 , 481.25, 12. ],
[ 271.25, 463. , 13. ],
[ 646.5 , 466.75, 14. ],
[ 739. , 442.75, 15. ],
[ 340.5 , 441.5 , 16. ],
[ 817.75, 421.5 , 17. ],
[ 423.75, 417.75, 18. ],
[ 202.5 , 406. , 19. ],
[ 519.25, 392.25, 20. ],
[ 257.5 , 382. , 21. ],
[ 619.25, 368.5 , 22. ],
[ 148. , 359.75, 23. ],
[ 324.5 , 356. , 24. ],
[ 713. , 347.75, 25. ],
[ 195. , 335. , 26. ],
[ 793.5 , 332.5 , 27. ],
[ 403.75, 328. , 28. ],
[ 249.25, 308. , 29. ],
[ 495.5 , 300.75, 30. ],
[ 314. , 279. , 31. ],
[ 764.25, 249.5 , 32. ],
[ 389.5 , 249.5 , 33. ],
[ 475. , 221.5 , 34. ],
[ 565.75, 199. , 35. ],
[ 802.75, 173.75, 36. ],
[ 733. , 176.25, 37. ]])
plt.figure()
plt.gca().invert_yaxis()
plt.plot(data[:,0], data[:,1], 'r.', markersize=14)
for idx in np.arange(data.shape[0]):
plt.text(data[idx,0]+10, data[idx,1]+10, str(int(data[idx,2])), size='large')
plt.show()
我们得到：
回到问题
如您所见，这些点或多或少呈网格状，您可以看到我们可以在这些点之间形成线条。具体可以看到有横向和纵向可以组成的线
比如你参考我问题背景部分的图片，我们可以看到有5组点可以水平分组。例如，点23、26、29、31、33、34、35、37、36为一组。第 19、21、24、28、30 和 32 点构成另一组，依此类推。同样在垂直方向上，我们可以看到26、19、12、3点为一组，29、21、13、5点为一组，依此类推。
要问的问题
我的问题是：假设点可以在任何方向，有什么方法可以分别成功地将水平分组和垂直分组中的点分组？
条件
每行必须至少三个点。如果小于此，则这不符合细分市场的条件。因此，点 36 和 10 不符合垂直线的条件，同样孤立点 23 不应该符合垂直线的条件，但它是第一个水平分组的一部分。
以上校准图案可以任意方向。然而，对于我正在处理的情况，你能得到的最糟糕的方向是你在上面的背景部分看到的。
预期输出
输出将是一对列表，其中第一个列表包含元素，其中每个元素为您提供一系列点 ID，这些点 ID 形成一条水平线。同样，第二个列表包含元素，其中每个元素为您提供一系列点 ID，这些点 ID 形成一条垂直线。
因此，水平序列的预期输出如下所示：
MATLAB
horiz_list = {[23, 26, 29, 31, 33, 34, 35, 37, 36], [19, 21, 24, 28, 30, 32], ...};
vert_list = {[26, 19, 12, 3], [29, 21, 13, 5], ....};
Python
horiz_list = [[23, 26, 29, 31, 33, 34, 35, 37, 36], [19, 21, 24, 28, 30, 32], ....]
vert_list = [[26, 19, 12, 3], [29, 21, 13, 5], ...]
我试过的
从算法上讲，我尝试的是撤消在这些点上经历的旋转。我已经执行 Principal Components Analysis 并尝试相对于计算的正交基向量投影点，以便这些点或多或少位于直矩形网格上。
一旦我有了它，做一些扫描线处理就是一件简单的事情，你可以根据水平或垂直坐标上的差异变化对点进行分组。您可以按 x 或 y 值对坐标进行排序，然后检查这些排序后的坐标并寻找较大的变化。一旦遇到这种变化，就可以将变化之间的点组合在一起以形成线条。对每个维度执行此操作将为您提供水平或垂直分组。
关于 PCA，这是我在 MATLAB 和 Python:
中所做的MATLAB
%# Step #1 - Get just the data - no IDs
data_raw = data(:,1:2);
%# Decentralize mean
data_nomean = bsxfun(@minus, data_raw, mean(data_raw,1));
%# Step #2 - Determine covariance matrix
%# This already decentralizes the mean
cov_data = cov(data_raw);
%# Step #3 - Determine right singular vectors
[~,~,V] = svd(cov_data);
%# Step #4 - Transform data with respect to basis
F = V.'*data_nomean.';
%# Visualize both the original data points and transformed data
figure;
plot(F(1,:), F(2,:), 'b.', 'MarkerSize', 14);
axis ij;
hold on;
plot(data(:,1), data(:,2), 'r.', 'MarkerSize', 14);
Python
import numpy as np
import numpy.linalg as la
# Step #1 and Step #2 - Decentralize mean
centroids_raw = data[:,:2]
mean_data = np.mean(centroids_raw, axis=0)
# Transpose for covariance calculation
data_nomean = (centroids_raw - mean_data).T
# Step #3 - Determine covariance matrix
# Doesn't matter if you do this on the decentralized result
# or the normal result - cov subtracts the mean off anyway
cov_data = np.cov(data_nomean)
# Step #4 - Determine right singular vectors via SVD
# Note - This is already V^T, so there's no need to transpose
_,_,V = la.svd(cov_data)
# Step #5 - Transform data with respect to basis
data_transform = np.dot(V, data_nomean).T
plt.figure()
plt.gca().invert_yaxis()
plt.plot(data[:,0], data[:,1], 'b.', markersize=14)
plt.plot(data_transform[:,0], data_transform[:,1], 'r.', markersize=14)
plt.show()
以上代码不仅重新投影数据，而且还将原始点和投影点一起绘制在一个图中。然而，当我尝试重新投影我的数据时，这是我得到的图：
红色的点是原始图像坐标，而蓝色的点被重新投影到基向量上以尝试消除旋转。它仍然不能完全完成工作。关于这些点仍然有一些方向，所以如果我尝试执行我的扫描线算法，来自下面的线的点用于水平跟踪或到侧面的点用于垂直跟踪将无意中分组，这是不正确的。
也许我想多了这个问题，但如果您对此有任何见解，我们将不胜感激。如果答案确实很棒，我会倾向于奖励高额赏金，因为我已经在这个问题上停留了很长一段时间。
我希望这个问题不是冗长的。如果您不知道如何解决这个问题，那么我感谢您花时间阅读我的问题。
期待您的任何见解。非常感谢！
虽然我无法提出比您已经尝试过的更好的方法来对任何给定的质心点列表进行分组，但我希望以下想法可以帮助您：
由于您对图像的内容非常具体（包含一个正方形区域），我想知道您是否真的需要根据 problem setup 中给出的数据对质心点进行分组，或者是否您也可以使用 Background to the problem 中描述的数据。由于您已经确定了每个检测到的正方形的角以及它们在给定正方形中的位置，因此在我看来，通过比较角坐标来确定给定正方形的邻居是非常准确的。
因此，为了找到任何正方形右邻的候选者，我建议您将该正方形的右上角和右下角与任何其他正方形（或其中的任何正方形）的左上角和左下角进行比较一定距离）。仅允许较小的垂直差异和稍大的水平差异，您可以 "match" 两个正方形，如果 两者 它们对应的角点足够靠近。
通过使用角之间允许的 vertical/horizontal 差异的上限，您甚至可以将这些边界内的最佳匹配方块指定为邻居
一个问题可能是你没有检测到所有的方块，所以在方块 30 和 32 之间有一个相当大的 space。因为你说你需要 'at least' 每行 3 个方块，您可以简单地忽略该水平线上的方块 32。如果这不适合您，您可以尝试匹配尽可能多的方块，然后使用先前计算的数据将 "missing" 方块分配给网格中的一个点：
在关于方块 32 的示例中，您会检测到它有上下邻居 27 和 37。此外，您应该已经能够确定方块 27 位于第 1 行内，而 37 位于第 3 行内，所以你可以将方块 32 分配给中间的 "best matching" 行，在这种情况下显然是 2。
这种通用方法基本上是您已经尝试过的方法，但应该更准确，因为您现在比较的是两条线的方向和距离，而不是简单地比较网格中两点的位置。
也作为您之前尝试的旁路节点 - 您可以使用黑色角线稍微校正图像的初始旋转吗？这可能会使进一步的失真算法（比如您在评论中与 knedlsepp 讨论的那些算法）更加准确。 （编辑：我刚才确实阅读了Parag的评论-通过线的角度比较点当然与事先旋转图像基本相同）
注意 1：它有许多设置 -> 对于其他图像可能需要更改这些设置才能获得您想要的结果请参阅 % 设置 - 使用这些值
注意 2：它不会找到您想要的所有行 -> 但它是一个起点....
要调用此函数，请在命令提示符中调用此函数：
>> [h, v] = testLines;
我们得到：
>> celldisp(h)
h{1} =
1 2 4 6 9 10
h{2} =
3 5 7 8 11 14 15 17
h{3} =
1 2 4 6 9 10
h{4} =
3 5 7 8 11 14 15 17
h{5} =
1 2 4 6 9 10
h{6} =
3 5 7 8 11 14 15 17
h{7} =
3 5 7 8 11 14 15 17
h{8} =
1 2 4 6 9 10
h{9} =
1 2 4 6 9 10
h{10} =
12 13 16 18 20 22 25 27
h{11} =
13 16 18 20 22 25 27
h{12} =
3 5 7 8 11 14 15 17
h{13} =
3 5 7 8 11 14 15
h{14} =
12 13 16 18 20 22 25 27
h{15} =
3 5 7 8 11 14 15 17
h{16} =
12 13 16 18 20 22 25 27
h{17} =
19 21 24 28 30
h{18} =
21 24 28 30
h{19} =
12 13 16 18 20 22 25 27
h{20} =
19 21 24 28 30
h{21} =
12 13 16 18 20 22 24 25
h{22} =
12 13 16 18 20 22 24 25 27
h{23} =
23 26 29 31 33 34 35
h{24} =
23 26 29 31 33 34 35 37
h{25} =
23 26 29 31 33 34 35 36 37
h{26} =
33 34 35 37 36
h{27} =
31 33 34 35 37
>> celldisp(v)
v{1} =
33 28 18 8 1
v{2} =
34 30 20 11 2
v{3} =
26 19 12 3
v{4} =
35 22 14 4
v{5} =
29 21 13 5
v{6} =
25 15 6
v{7} =
31 24 16 7
v{8} =
37 32 27 17 9
还会生成一个图形，通过每个适当的点集绘制线：
function [horiz_list, vert_list] = testLines
global counter;
global colours;
close all;
data = [ 475. , 605.75, 1.;
571. , 586.5 , 2.;
233. , 558.5 , 3.;
669.5 , 562.75, 4.;
291.25, 546.25, 5.;
759. , 536.25, 6.;
362.5 , 531.5 , 7.;
448. , 513.5 , 8.;
834.5 , 510. , 9.;
897.25, 486. , 10.;
545.5 , 491.25, 11.;
214.5 , 481.25, 12.;
271.25, 463. , 13.;
646.5 , 466.75, 14.;
739. , 442.75, 15.;
340.5 , 441.5 , 16.;
817.75, 421.5 , 17.;
423.75, 417.75, 18.;
202.5 , 406. , 19.;
519.25, 392.25, 20.;
257.5 , 382. , 21.;
619.25, 368.5 , 22.;
148. , 359.75, 23.;
324.5 , 356. , 24.;
713. , 347.75, 25.;
195. , 335. , 26.;
793.5 , 332.5 , 27.;
403.75, 328. , 28.;
249.25, 308. , 29.;
495.5 , 300.75, 30.;
314. , 279. , 31.;
764.25, 249.5 , 32.;
389.5 , 249.5 , 33.;
475. , 221.5 , 34.;
565.75, 199. , 35.;
802.75, 173.75, 36.;
733. , 176.25, 37.];
figure; hold on;
axis ij;
% Change due to Benoit_11
scatter(data(:,1), data(:,2),40, 'r.'); text(data(:,1)+10, data(:,2)+10, num2str(data(:,3)));
text(data(:,1)+10, data(:,2)+10, num2str(data(:,3)));
% Process your data as above then run the function below(note it has sub functions)
counter = 0;
colours = 'bgrcmy';
[horiz_list, vert_list] = findClosestPoints ( data(:,1), data(:,2) );
function [horiz_list, vert_list] = findClosestPoints ( x, y )
% calc length of points
nX = length(x);
% set up place holder flags
modelledH = false(nX,1);
modelledV = false(nX,1);
horiz_list = {};
vert_list = {};
% loop for all points
for p=1:nX
% have we already modelled a horizontal line through these?
% second last param - true - horizontal, false - vertical
if modelledH(p)==false
[modelledH, index] = ModelPoints ( p, x, y, modelledH, true, true );
horiz_list = [horiz_list index];
else
[~, index] = ModelPoints ( p, x, y, modelledH, true, false );
horiz_list = [horiz_list index];
end
% make a temp copy of the x and y and remove any of the points modelled
% from the horizontal -> this is to avoid them being found in the
% second call.
tempX = x;
tempY = y;
tempX(index) = NaN;
tempY(index) = NaN;
tempX(p) = x(p);
tempY(p) = y(p);
% Have we found a vertial line?
if modelledV(p)==false
[modelledV, index] = ModelPoints ( p, tempX, tempY, modelledV, false, true );
vert_list = [vert_list index];
end
end
end
function [modelled, index] = ModelPoints ( p, x, y, modelled, method, fullRun )
% p - row in your original data matrix
% x - data(:,1)
% y - data(:,2)
% modelled - array of flags to whether rows have been modelled
% method - horizontal or vertical (used to calc graadients)
% fullRun - full calc or just to get indexes
% this could be made better by storing the indexes of each horizontal in the method above
% Settings - play around with these values
gradDelta = 0.2; % find points where gradient is less than this value
gradLimit = 0.45; % if mean gradient of line is above this ignore
numberOfPointsToCheck = 7; % number of points to check when look along the line
% to find other points (this reduces chance of it
% finding other points far away
% I optimised this for your example to be 7
% Try varying it and you will see how it effect the result.
% Find the index of points which are inline.
[index, grad] = CalcIndex ( x, y, p, gradDelta, method );
% check gradient of line
if abs(mean(grad))>gradLimit
index = [];
return
end
% add point of interest to index
index = [p index];
% loop through all points found above to find any other points which are in
% line with these points (this allows for slight curvature
combineIndex = [];
for ii=2:length(index)
% Find inedex of the points found above (find points on curve)
[index2] = CalcIndex ( x, y, index(ii), gradDelta, method, numberOfPointsToCheck, grad(ii-1) );
% Check that the point on this line are on the original (i.e. inline -> not at large angle
if any(ismember(index,index2))
% store points found
combineIndex = unique([index2 combineIndex]);
end
end
% copy to index
index = combineIndex;
if fullRun
% do some plotting
% TODO: here you would need to calculate your arrays to output.
xx = x(index);
[sX,sOrder] = sort(xx);
% Check its found at least 3 points
if length ( index(sOrder) ) > 2
% flag the modelled on the points found
modelled(index(sOrder)) = true;
% plot the data
plot ( x(index(sOrder)), y(index(sOrder)), colours(mod(counter,numel(colours)) + 1));
counter = counter + 1;
end
index = index(sOrder);
end
end
function [index, gradCheck] = CalcIndex ( x, y, p, gradLimit, method, nPoints2Consider, refGrad )
% x - data(:,1)
% y - data(:,2)
% p - point of interest
% method (x/y) or (y\x)
% nPoints2Consider - only look at N points (options)
% refgrad - rather than looking for gradient of closest point -> use this
% - reference gradient to find similar points (finds points on curve)
nX = length(x);
% calculate gradient
for g=1:nX
if method
grad(g) = (x(g)-x(p))\(y(g)-y(p));
else
grad(g) = (y(g)-y(p))\(x(g)-x(p));
end
end
% find distance to all other points
delta = sqrt ( (x-x(p)).^2 + (y-y(p)).^2 );
% set its self = NaN
delta(delta==min(delta)) = NaN;
% find the closest points
[m,order] = sort(delta);
if nargin == 7
% for finding along curve
% set any far away points to be NaN
grad(order(nPoints2Consider+1:end)) = NaN;
% find the closest points to the reference gradient within the allowable limit
index = find(abs(grad-refGrad)<gradLimit==1);
% store output
gradCheck = grad(index);
else
% find the points which are closes to the gradient of the closest point
index = find(abs(grad-grad(order(1)))<gradLimit==1);
% store gradients to output
gradCheck = grad(index);
end
end
end
我正在使用已发布图像的裁剪版本作为输入。这里我只讨论网格的方向可以被认为接近 horizontal/vertical 的情况。这可能无法完全解决您的范围，但我认为它可能会给您一些指导。
将图像二值化，以填充扭曲的方块。这里我使用了一个简单的 Otsu 阈值。然后对这个二值图像进行距离变换。
在距离变换图像中，我们看到正方形之间的间隙为峰值。
要获得水平方向的线，取距离图像的每个 列 的局部最大值，然后找到连通分量。
要获得垂直方向的线，取距离图像的每个 行 的局部最大值，然后找到连通分量。
下图显示了这样找到的水平线和垂直线，角点为圆圈。
对于相当长的连通分量，可以拟合一条曲线（直线或多项式），然后对角点进行分类，比如根据到曲线的距离、该点在曲线的哪一侧等.
我是用 Matlab 做的。我没有尝试曲线拟合和分类部分。
clear all;
close all;
im = imread('0RqUd-1.jpg');
gr = rgb2gray(im);
% any preprocessing to get a binary image that fills the distorted squares
bw = ~im2bw(gr, graythresh(gr));
di = bwdist(bw); % distance transform
di2 = imdilate(di, ones(3)); % propagate max
corners = corner(gr); % simple corners
% find regional max for each column of dist image
regmxh = zeros(size(di2));
for c = 1:size(di2, 2)
regmxh(:, c) = imregionalmax(di2(:, c));
end
% label connected components
ccomph = bwlabel(regmxh, 8);
% find regional max for each row of dist image
regmxv = zeros(size(di2));
for r = 1:size(di2, 1)
regmxv(r, :) = imregionalmax(di2(r, :));
end
% label connected components
ccompv = bwlabel(regmxv, 8);
figure, imshow(gr, [])
hold on
plot(corners(:, 1), corners(:, 2), 'ro')
figure, imshow(di, [])
figure, imshow(label2rgb(ccomph), [])
hold on
plot(corners(:, 1), corners(:, 2), 'ro')
figure, imshow(label2rgb(ccompv), [])
hold on
plot(corners(:, 1), corners(:, 2), 'ro')
要获得任意方向网格的这些线，您可以将距离图像视为图形并找到最佳路径。请参阅 here 了解基于图形的不错方法。