将一组点集合在一起的算法
ray*_*ica 37 python algorithm matlab image image-processing
注意:我将这个问题放在MATLAB和Python标签中,因为我是这些语言中最精通的.但是,我欢迎任何语言的解决方案.
问题序言
我用鱼眼镜头拍了一张照片.此图像由带有一堆方形对象的图案组成.我想用这个图像做的是检测每个方块的质心,然后使用这些点来执行图像的不失真 - 具体来说,我正在寻找正确的失真模型参数.应该注意,并非所有的方块都需要被检测到.只要他们中的绝大多数都是,那就完全没问题......但这不是这篇文章的重点.我已经编写过参数估计算法,但问题是它需要在图像中看起来共线的点.
我想问的基本问题给出了这些要点,将它们组合在一起的最佳方法是什么,以便每个组由水平线或垂直线组成？
我的问题的背景
对于我提出的问题,这并不重要,但如果您想知道我从哪里得到我的数据并进一步理解我提出的问题,请阅读.如果您不感兴趣,那么您可以直接跳到下面的问题设置部分.
我正在处理的图像示例如下所示:
这是一张960 x 960的图像.图像最初的分辨率较高,但我对图像进行二次采样以加快处理时间.如您所见,有一堆方形图案分散在图像中.此外,我计算的质心是关于上面的二次采样图像.
我设置的用于检索质心的管道如下:
- 执行Canny边缘检测
- 专注于最小化误报的感兴趣区域.这个感兴趣的区域基本上是没有任何覆盖其中一侧的黑色胶带的正方形.
- 找到所有截然不同的轮廓
对于每个不同的闭合轮廓......
一个.执行Harris角点检测
湾 确定结果是否有4个角点
C.如果是这样,则该轮廓属于正方形并找到该形状的质心
d.如果没有,则跳过此形状
将步骤#4中检测到的所有质心放入矩阵中进行进一步检查.
以下是上图中的示例结果.每个检测到的正方形具有四个点,根据它相对于正方形本身的位置进行颜色编码.对于我检测到的每个质心,我会在图像本身的质心处写一个ID.
使用上面的图像,有37个检测到的正方形.
问题设置
假设我有一些图像像素点存储在N x 3矩阵中.前两列是x(水平)和y(垂直)坐标,在图像坐标空间中,y坐标反转,这意味着正向下y移动.第三列是与该点相关联的ID.
下面是一些用MATLAB编写的代码,它采用这些点,将它们绘制到2D网格上,并用矩阵的第三列标记每个点.如果您阅读上述背景,则这些是我上面概述的算法检测到的点.
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.
因此,水平序列的预期输出看起来像这样:
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], ....};
蟒蛇
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], ...]
我试过了什么
从算法上讲,我所尝试的是撤消在这些点上经历的旋转.我已经执行了主成分分析,并且我尝试相对于计算的正交基矢量投射点,使得这些点或多或少地在直的矩形网格上.
一旦我有了这个,只需要做一些扫描线处理,你可以根据水平或垂直坐标上的差异变化对点进行分组.您可以使用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);
蟒蛇
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()
上面的代码不仅重新投影数据,而且还将原始点和投影点一起绘制在一个图中.但是,当我尝试重新投影我的数据时,这是我得到的情节:
红色的点是原始图像坐标,而蓝色的点被重新投影到基础矢量上以尝试去除旋转.它仍然不能完成这项工作.关于这些点仍然存在一些方向,因此如果我尝试使用扫描线算法,则从下面的行进行水平跟踪或者用于垂直跟踪的一侧的点将被无意中分组,这是不正确的.
也许我正在推翻这个问题,但是你对此有任何见解都会非常感激.如果答案确实很棒,我会倾向于给予高额奖励,因为我已经坚持这个问题很长一段时间了.
我希望这个问题不会长篇大论.如果您不知道如何解决这个问题,那么我感谢您抽出时间阅读我的问题.
期待您可能拥有的任何见解.非常感谢!
mat*_*gui 13
注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
虽然我不能建议一个更好的方法来分组任何给定的质心点列表,而不是你已经尝试过的,我希望以下想法可以帮助你:
由于您非常关注图像的内容(包含正方形区域),我想知道您是否实际上需要根据您给出的数据对质心点进行分组problem setup,或者您是否也可以使用所描述的数据Background to the problem.由于您已经确定了每个检测到的正方形的角以及它们在该给定正方形中的位置,因此在我看来,通过比较角坐标来确定给定正方形的邻居是非常准确的.
因此,为了找到任何正方形的右邻居的任何候选者,我建议你将该正方形的右上角和右下角与任何其他正方形(或某一距离内的任何正方形)的左上角和左下角进行比较.只允许小的垂直差异和略微更大的水平差异,如果两个对应的角点足够接近,则可以"匹配"两个正方形.
通过使用拐角之间允许的垂直/水平差异的上限,您甚至可以将这些边界内的最佳匹配方块指定为邻居
一个问题可能是你没有检测到所有的正方形,所以在30和32之间有一个相当大的空间.因为你说你需要'每行至少'3个方格,你可能只需忽略它该水平线上的正方形32.如果这不是您的选项,您可以尝试匹配尽可能多的方块,然后使用以前计算的数据将"缺失"方块分配给网格中的某个点:
在关于方形32的示例中,您将检测到它具有上下邻居27和37.此外,您应该能够确定方形27位于第1行中,而37位于第3行中,因此您可以指定方形32到中间的"最佳匹配"行,在这种情况下显然是2.
这种通用的方法基本上就是你已经尝试过的方法,但是由于你现在正在比较两条线的方向和距离,而不是简单地比较网格中两点的位置,所以应该更准确.
此外,作为您之前尝试的副节点 - 您是否可以使用黑色边线来纠正图像的初始旋转？这可能使得进一步的失真算法(如您在评论中与knedlsepp讨论过的那样)更准确.(编辑:我刚才读过Parag的评论 - 比较线条的角度当然与预先旋转图像基本相同)
我正在使用已发布图片的裁剪版本作为输入.这里我只讨论网格方向可以被认为接近水平/垂直的情况.这可能无法完全解决您的范围,但我认为它可能会给您一些指示.
将图像二值化以便填充扭曲的正方形.在这里,我使用简单的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')
要为任意定向网格获取这些线,您可以将距离图像视为图形并找到最佳路径.请参阅此处了解基于图表的方法.
| 归档时间： |
|
| 查看次数： |
2540 次 |
| 最近记录： |