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

2023-09-10 23:46:59 作者:酒醒梦一场

注意:我把这个问题同时在MATLAB和Python的标签,我最精通这些语言。不过,我欢迎任何语言解决方案。

我已经拍摄的图像与鱼眼镜头。此图像由一堆方形物体的图案。我想要做的这个形象是每个检测这些广场的心,然后用这些点来进行图像的undistortion - 确切地说,我寻找合适的失真模型参数。应当指出的是,并非所有的平方需要检测。只要一个很好的大部分是,那么这是完全罚款....但不是这个职位的地步。参数估计算法我已写的,但问题是,它要求出现共线图像中的点。

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

一个。执行Harris角点检测

a. Perform a Harris Corner Detection

乙。确定是否结果有4个角点

b. Determine if the result has 4 corner points

℃。如果这样做,那么这个轮廓属于一个正方形,并​​发现这种形状的心

c. If this does, then this contour belonged to a square and find the centroid of this shape

Ð。如果没有,则跳过这种形状

d. If it doesn't, then skip this shape

将步骤#4所有检测到的重心变成一个矩阵作进一步检查。

Place all detected centroids from Step #4 into a matrix for further examination.

下面是与上面的图像的例子的结果。根据那里它是相对于该方本身位置的每一个检测的正方形有四个点色$ C $光盘。对于我所检测到的每个质心,我写的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×3 矩阵的一些图像的像素点。前两列是 X (水平)和(垂直)坐标,其中在图像坐标空间,在坐标倒,这意味着正向下移动。第三列是与点相关联的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.

下面是一些code,MATLAB写的是将这些点,他们绘制到2D网格和标签的每个点与基体的第三列。如果你看了上面的背景下,这些都是我的上述算法检测到的点。

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 ,我们有:

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.

我的问题是这样的:?什么是一个方法,可以成功地在水平分组和垂直分组分开,因为该点可以在任何方向群分的

My question is this: What is a method that can successfully group points in horizontal groupings and vertical groupings separately, given that the points could be in any orientation?

有必须至少有三个每行分。如果说有什么不足的是,那么这个没有资格作为一个段。因此,在点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

回归 分类与聚类 三大方向剖解机器学习算法的优缺点 人工智能 算法与数学之美 CSDN博客

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 的值,然后检查这些排序坐标,寻找一个大的更改。当你遇到这样的变化,那么你就可以在变化的组分在一起,形成你的线条。对于这样做是为了每一个维度会给你无论是水平或垂直的分组。

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

上面code不仅进行投影变换的数据,但它也在单个图中绘出无论原点和投影点一起。但是,当我试图重新投影我的数据,这是我得到的情节:

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.

也许我该得太多问题,但您对此有何看法将大大AP preciated。如果答案确实是一流的,我会倾向于授予高奖金,因为我已经被困在这个问题上相当长的一段时间。

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:它有许多设置 - >这对于其他图像可能需要改变,以得到你想要的结果的参见%设置 - 打这些值,的

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