将遵循一个方向的一组点组合在一起的算法 [英] Algorithm to group sets of points together that follow a direction
问题描述
注意:我将这个问题放在 MATLAB 和 Python 标签中，因为我最精通这些语言.但是，我欢迎任何语言的解决方案.
我用鱼眼镜头拍了一张照片.该图像由带有一堆方形对象的图案组成.我想要对这个图像做的是检测每个正方形的质心，然后使用这些点来执行图像的不失真 - 具体来说，我正在寻找正确的失真模型参数.应该注意的是，并不是所有的方块都需要被检测到.只要他们中的大多数是，那完全没问题……但这不是这篇文章的重点.参数估计算法我已经写好了，但问题是它需要在图像中出现共线的点.
I have taken an image with a fisheye lens. This image consists of a pattern with a bunch of square objects. What I want to do with this image is detect the centroid of each of these squares, then use these points to perform an undistortion of the image - specifically, I am seeking the right distortion model parameters. It should be noted that not all of the squares need to be detected. As long as a good majority of them are, then that's totally fine.... but that isn't the point of this post. The parameter estimation algorithm I have already written, but the problem is that it requires points that appear collinear in the image.
我想问的基本问题是给定这些点，将它们组合在一起以使每组由一条水平线或一条垂直线组成的最佳方法是什么?
The base question I want to ask is given these points, what is the best way to group them together so that each group consists of a horizontal line or vertical line?
这对于我提出的问题并不重要，但如果您想知道我从哪里获得数据并进一步了解我提出的问题，请阅读.如果您不感兴趣，可以直接跳到下面的问题设置部分.
This isn't really important with regards to the question I'm asking, but if you'd like to know where I got my data from and to further understand the question I'm asking, please read. If you're not interested, then you can skip right to the Problem setup section below.
我正在处理的图像示例如下所示:
An example of an image I am dealing with is shown below:
这是一张 960 x 960 的图像.该图像最初具有更高的分辨率，但我对图像进行了二次采样以加快处理时间.如您所见，图像中散布着一堆方形图案.此外，我计算的质心是关于上述二次采样图像的.
It is a 960 x 960 image. The image was originally higher resolution, but I subsample the image to facilitate faster processing time. As you can see, there are a bunch of square patterns that are dispersed in the image. Also, the centroids I have calculated are with respect to the above subsampled image.
我为检索质心而设置的管道如下:
The pipeline I have set up to retrieve the centroids is the following:
- 执行 Canny 边缘检测
- 专注于可最大程度减少误报的感兴趣区域.这个感兴趣的区域基本上是正方形，没有任何覆盖其一侧的黑色胶带.
- 找出所有不同的闭合轮廓
对于每个不同的闭合轮廓...
- Perform a Canny Edge Detection
- Focus on a region of interest that minimizes false positives. This region of interest is basically the squares without any of the black tape that covers one of their sides.
- Find all distinct closed contours
For each distinct closed contour...
一个.执行哈里斯角检测
a. Perform a Harris Corner Detection
B.判断结果是否有4个角点
b. Determine if the result has 4 corner points
c.如果是这样，那么这个轮廓属于一个正方形并找到这个形状的质心
c. If this does, then this contour belonged to a square and find the centroid of this shape
d.如果没有，则跳过此形状
d. If it doesn't, then skip this shape
将第 4 步中所有检测到的质心放入矩阵中以供进一步检查.
Place all detected centroids from Step #4 into a matrix for further examination.
这是上图的示例结果.每个检测到的正方形都有四个点，根据它相对于正方形本身的位置进行颜色编码.对于我检测到的每个质心，我会在该质心在图像本身的位置写一个 ID.
Here's an example result with the above image. Each detected square has the four points colour coded according to the location of where it is with respect to the square itself. For each centroid that I have detected, I write an ID right where that centroid is in the image itself.
通过上图，检测到 37 个方块.
With the above image, there are 37 detected squares.
假设我有一些图像像素点存储在 N x 3 矩阵中.前两列是 x(水平)和 y(垂直)坐标，其中在图像坐标空间中，y 坐标是 倒置，表示正的 y 向下移动.第三列是与点关联的 ID.
Suppose I have some image pixel points stored in a N x 3 matrix. The first two columns are the x (horizontal) and y (vertical) coordinates where in image coordinate space, the y coordinate is inverted, which means that positive y moves downwards. The third column is an ID associated with the point.
以下是一些用 MATLAB 编写的代码，它们采用这些点，将它们绘制到二维网格上，并用矩阵的第三列标记每个点.如果您阅读了上述背景，这些就是我上面概述的算法检测到的点.
Here is some code written in MATLAB that takes these points, plots them onto a 2D grid and labels each point with the third column of the matrix. If you read the above background, these are the points that were detected by my algorithm outlined above.
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，我们有:
Similarly in Python, using numpy and matplotlib, we have:
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()
我们得到:
如您所见，这些点或多或少呈网格状，您可以看到我们可以在这些点之间形成连线.具体来说，可以看到有横竖都可以形成的线条.
As you can see, these points are more or less in a grid pattern and you can see that we can form lines between the points. Specifically, you can see that there are lines that can be formed horizontally and vertically.
例如，如果你参考我问题的背景部分中的图像，我们可以看到有5组点可以按水平方式分组.例如，点 23、26、29、31、33、34、35、37 和 36 形成一组.点 19、21、24、28、30 和 32 形成另一个组，依此类推.同理，在垂直方向上，我们可以看到点26、19、12和3组成一组，点29、21、13和5组成另一组，以此类推.
For example, if you reference the image in the background section of my problem, we can see that there are 5 groups of points that can be grouped in a horizontal manner. For example, points 23, 26, 29, 31, 33, 34, 35, 37 and 36 form one group. Points 19, 21, 24, 28, 30 and 32 form another group and so on and so forth. Similarly in a vertical sense, we can see that points 26, 19, 12 and 3 form one group, points 29, 21, 13 and 5 form another group and so on.
我的问题是:考虑到点可以处于任何方向，有什么方法可以成功地将点分别分为水平分组和垂直分组?
每行必须至少有三个点.如果低于此值，则这不符合细分的条件.因此，点 36 和 10 不符合垂直线的条件，类似地，孤立点 23 不应该是垂直线，但它是第一个水平分组的一部分.
There must be at least three points per line. If there is anything less than that, then this does not qualify as a segment. Therefore, the points 36 and 10 don't qualify as a vertical line, and similarly the isolated point 23 shouldn't quality as a vertical line, but it is part of the first horizontal grouping.
上述校准图案可以在任何方向.然而，对于我正在处理的问题，你能得到的最糟糕的方向是你在背景部分看到的.
The above calibration pattern can be in any orientation. However, for what I'm dealing with, the worst kind of orientation you can get is what you see above in the background section.
<小时>
预期产出
输出将是一对列表，其中第一个列表包含元素，其中每个元素为您提供形成水平线的点 ID 序列.类似地，第二个列表包含元素，其中每个元素为您提供形成垂直线的点 ID 序列.
Expected Output
The output would be a pair of lists where the first list has elements where each element gives you a sequence of point IDs that form a horizontal line. Similarly, the second list has elements where each element gives you a sequence of point IDs that form a vertical line.
因此，水平序列的预期输出如下所示:
Therefore, the expected output for the horizontal sequences would look something like this:
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], ...]
我尝试了什么
从算法上讲，我尝试的是撤消在这些点上经历的旋转.我已经执行了 主成分分析，并尝试根据计算出的正交基向量来投影点以便这些点或多或少位于直的矩形网格上.
What I have tried
Algorithmically, what I have tried is to undo the rotation that is experienced at these points. I've performed Principal Components Analysis and I tried projecting the points with respect to the computed orthogonal basis vectors so that the points would more or less be on a straight rectangular grid.
一旦我有了它，这只是做一些扫描线处理的简单问题，您可以根据水平或垂直坐标上的差异变化对点进行分组.您可以按 x 或 y 值对坐标进行排序，然后检查这些排序后的坐标并寻找较大的变化.遇到此更改后，您可以将更改之间的点组合在一起以形成您的线条.针对每个维度执行此操作将为您提供水平或垂直分组.
Once I have that, it's just a simple matter of doing some scanline processing where you could group points based on a differential change on either the horizontal or vertical coordinates. You'd sort the coordinates by either the x or y values, then examine these sorted coordinates and look for a large change. Once you encounter this change, then you can group points in between the changes together to form your lines. Doing this with respect to each dimension would give you either the horizontal or vertical groupings.
关于 PCA，这是我在 MATLAB 和 Python 中所做的:
With regards to PCA, here's what I did in MATLAB and Python:
%# 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()
<小时>
上面的代码不仅对数据进行了重新投影，而且还将原始点和投影点绘制在一个图形中.然而，当我尝试重新投影我的数据时，这是我得到的图:
The above code not only reprojects the data, but it also plots both the original points and the projected points together in a single figure. However, when I tried reprojecting my data, this is the plot I get:
红色的点是原始图像坐标，而蓝色的点被重新投影到基向量上以尝试去除旋转.它仍然不能很好地完成工作.仍然有一些关于点的方向，所以如果我尝试执行我的扫描线算法，则来自下方用于水平跟踪的线或用于垂直跟踪的侧面的点将被无意中分组，这是不正确的.
The points in red are the original image coordinates while the points in blue are reprojected onto the basis vectors to try and remove the rotation. It still doesn't quite do the job. There is still some orientation with respect to the points so if I tried to do my scanline algorithm, points from the lines below for horizontal tracing or to the side for vertical tracing would be inadvertently grouped and this isn't correct.
也许我想得太多了，但您对此的任何见解都将不胜感激.如果答案确实很棒，我会倾向于给予高额赏金，因为我已经在这个问题上困扰了很长时间.
Perhaps I'm overthinking the problem, but any insights you have regarding this would be greatly appreciated. If the answer is indeed superb, I would be inclined to award a high bounty as I've been stuck on this problem for quite some time.
我希望这个问题不是冗长的.如果您不知道如何解决这个问题，那么我感谢您花时间阅读我的问题.
I hope this question wasn't long winded. If you don't have an idea of how to solve this, then I thank you for your time in reading my question regardless.
期待您的任何见解.非常感谢！
Looking forward to any insights that you may have. Thanks very much!
推荐答案
注意 1:它有许多设置 -> 对于其他图像可能需要更改以获得您想要的结果参见 % Settings - play围绕这些价值观
Note 1: It has a number of settings -> which for other images may need to altered to get the result you want see % Settings - play around with these values
注意 2:它没有找到您想要的所有行 -> 但它是一个起点....
Note 2: It doesn't find all of the lines you want -> but its a starting point....
要调用此函数，请在命令提示符中调用它:
To call this function, invoke this in the command prompt:
>> [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
还会生成一个图形，通过每个适当的点集绘制线条:
A figure is also generated that draws the lines through each proper set of points:
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 (yx)
% 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
这篇关于将遵循一个方向的一组点组合在一起的算法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
