将遵循方向的点集分组在一起的算法

2024-03-30

注意：我将这个问题放在 MATLAB 和 Python 标签中，因为我最精通这些语言。但是，我欢迎任何语言的解决方案。

问题序言

我用鱼眼镜头拍摄了一张图像。该图像由带有一堆方形物体的图案组成。我想要对该图像执行的操作是检测每个正方形的质心，然后使用这些点对图像进行去畸变 - 具体来说，我正在寻找正确的畸变模型参数。应该注意的是，并不是所有的方块都需要被检测到。只要他们中的大多数人都是这样，那就完全没问题......但这不是这篇文章的重点。我已经编写了参数估计算法，但问题是它需要图像中出现共线的点。

我想问的基本问题是给定这些点，将它们分组在一起以便每个组由水平线或垂直线组成的最佳方法是什么？

我的问题的背景

对于我要问的问题来说，这并不重要，但如果您想知道我从哪里获得数据并进一步了解我要问的问题，请阅读。如果你不感兴趣，那么你可以直接跳到问题设置下面的部分。

我正在处理的图像示例如下所示：

这是一张 960 x 960 的图像。该图像最初的分辨率较高，但我对图像进行了二次采样，以加快处理时间。正如您所看到的，图像中分散着一堆方形图案。另外，我计算的质心是相对于上面的二次采样图像的。

我设置的用于检索质心的管道如下：

执行 Canny 边缘检测
专注于最大程度减少误报的感兴趣区域。该感兴趣区域基本上是没有任何黑色胶带覆盖其一侧的正方形。
找到所有不同的闭合轮廓
对于每个不同的闭合轮廓...

A。执行 Harris 角点检测

b.判断结果是否有4个角点

C。如果是这样，那么这个轮廓属于一个正方形并找到这个形状的质心

d.如果没有，则跳过此形状
将步骤 #4 中检测到的所有质心放入矩阵中以供进一步检查。

这是上图的示例结果。每个检测到的方块都有四个点，根据其相对于方块本身的位置进行颜色编码。对于我检测到的每个质心，我会在图像本身中该质心所在的位置写入一个 ID。

在上图中，检测到了 37 个正方形。

问题设置

假设我有一些图像像素点存储在N x 3矩阵。前两列是x（水平）和y（垂直）坐标，在图像坐标空间中，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 and 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()

We get:

回到问题

正如您所看到的，这些点或多或少呈网格图案，您可以看到我们可以在点之间形成线条。具体来说，你可以看到有可以水平和垂直形成的线条。

例如，如果您引用我的问题的背景部分中的图像，我们可以看到有 5 组点可以以水平方式分组。例如，点 23、26、29、31、33、34、35、37 和 36 形成一组。点 19、21、24、28、30 和 32 形成另一组，依此类推。同样，在垂直方向上，我们可以看到点 26、19、12 和 3 形成一组，点 29、21、13 和 5 形成另一组，依此类推。

要问的问题

我的问题是这样的：考虑到点可以处于任何方向，有什么方法可以成功地分别将水平分组和垂直分组中的点分组？

状况

There 必须至少为三个每行点数。如果有任何低于此的内容，则这不符合分段的资格。因此，点 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], ....};

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], ...]

我尝试过的

从算法上讲，我尝试的是撤消在这些点上经历的旋转。我表演过主成分分析 http://en.wikipedia.org/wiki/Principal_component_analysis我尝试根据计算的正交基向量投影这些点，以便这些点或多或少位于直矩形网格上。

一旦我有了它，这只是进行一些扫描线处理的简单问题，您可以根据水平或垂直坐标上的差异变化对点进行分组。您可以通过以下任一方式对坐标进行排序x or 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()

上面的代码不仅重新投影数据，而且还将原始点和投影点一起绘制在一张图中。然而，当我尝试重新投影我的数据时，这是我得到的图：

红色的点是原始图像坐标，而蓝色的点被重新投影到基础向量上以尝试消除旋转。它仍然没有完全完成这项工作。这些点仍然存在一些方向，因此如果我尝试执行扫描线算法，则水平跟踪下方的线或垂直跟踪侧面的点将无意中分组，这是不正确的。

也许我对这个问题想得太多了，但是如果您对此有任何见解，我们将不胜感激。如果答案确实很棒，我会倾向于给予高额赏金，因为我已经在这个问题上停留了很长一段时间了。

我希望这个问题不是冗长的。如果您不知道如何解决这个问题，那么无论如何，我感谢您花时间阅读我的问题。

期待您的任何见解。非常感谢！

注 1：它有许多设置 -> 对于其他图像可能需要更改这些设置才能获得您想要的结果请参阅 % 设置 - 尝试使用这些值

注 2：它没有找到您想要的所有行 -> 但它是一个起点......

要调用此函数，请在命令提示符中调用：

>> [h, v] = testLines;

We get:

>> 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

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)