CN106709499A - 基于Canny算子和Hillbert‑Huang变换的SIFT图像特征点提取方法 - Google Patents

Abstract

本发明提供了一种基于Canny算子和Hillbert‑Huang变换的SIFT图像特征点提取方法，涉及图像处理领域，是基于Canny算子和Hillbert‑Huang变换，改进的SIFT图像特征点提取的方法，由于采用了Canny算子和Hillbert‑Huang变换融合的技术方案，改进了原有SIFT算法实时性不高，特征点有时较少以及对边缘光滑的目标无法准确提取特征点的缺陷，能够有效地增强了SIFT特征点的抗噪声能力，提高了算法的稳定性，为三维重建和目标定位等提供了良好的基础。

Description

基于Canny算子和Hillbert-Huang变换的SIFT图像特征点提 取方法

技术领域

本发明涉及图像处理领域，尤其是涉及一种图像特征点匹配的方法。

背景技术

图像匹配作为工业、农业、商业以及军事领域中不可缺少的一项重要研究手段，在近年来一直被人们不断进行深入探索，在国内外主要有以下几个手段：颜色或灰度特征提取、纹理边缘特征提取、图像代数特征提取、图像变换系数特征提取，其中，SIFT(Scale-invariant feature transform)算法以其受光照强度影响小，对目标遮挡处理结果好，良好的鲁棒性以及实时高速的特点，在特征点检测中被许多学者广泛应用。

SIFT算法是由Lowe在2004年提出的一种图像特征点匹配算法，通过确定尺度空间、检测关键点位置、确定关键点位置以及对选取关键点描述几个步骤，最终实现图像匹配，然而SIFT算法也存在缺陷，SIFT算法中所利用DOG算子进行边缘检测后所得到的关键点，由于边缘响应的原因，有很多并不是真正的关键点。

在边缘检测中，Canny算子具有滤波、增强以及检测等多方面优化，受到噪声干扰性小，但是Canny算子处理步骤中，利用高斯算子对图像平滑处理时，会出现过度平滑和边缘移位的情况，二维的Hillbert-Huang变换用于图像分解中，可以将图像分解成对应的内蕴模函数，使边缘轮廓更加清晰。

发明内容

为了克服现有技术的不足，本发明提供一种基于Canny算子和Hillbert-Huang变换的改进的SIFT图像特征点提取的方法，对SIFT图像匹配算法进行改进，通过实验结果验证，图像特征点提取准确率有较大提高。

本发明解决其技术问题所采用的技术方案包括以下步骤：

步骤1：对输入的原图像I(x,y)定义尺度空间L(x，y，σ)，构建高斯金字塔：

L(x，y，σ)＝G(x，y，σ)*I(x，y) (1)

其中G(x，y，σ)为高斯函数，m，n为高斯模板的维度，σ为尺度空间因子；

步骤2：构建高斯差分金字塔

(1)将图像金字塔分o组，一组称为一个Octave，每组分为多层，组内总层数为S，在每一组图像的顶层用高斯模糊生成3幅图像，高斯金字塔每组则有S+3层图像；

(2)对原图像采用隔点采样法进行降采样，具体实施为：比例因子为2对原图像每隔一行一列取一点；

(3)确定高斯金字塔的层数：原图像为金字塔的底部第一层，每次降采样所得到的新图像为金字塔的上一层的一张图像，每个金字塔共n层，金字塔的层数n根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：

n＝log₂{min(M，N)}-t (3)

其中，M，N为原图像的像素值，t为塔顶图像的最小维数的对数值；

(4)计算高斯金字塔组内某一层图像的尺度

其中σ₀是基准层尺度，根据SIFT算法取值为1.6，o为组octave的索引，s为组内层的索引，根据SIFT算法取值为3，关键点的尺度坐标σ即关键点所在的组和组内的层，不同组相同层的组内尺度坐标相同；

(5)计算DOG算子

D(x，y，σ)＝(G(x，y，kσ)-G(x，y，σ))*I(x，y)＝L(x，y，kσ)-L(x，y，σ) (5)

其中k为组内总层数的倒数，即：

步骤3：对图像进行空间极值点的检测，具体步骤如下：

将步骤2中计算得到的DOG算子D(x，y，σ)泰勒展开，得到

其中X＝(x，y，σ)^T，使公式(5)的泰勒展开式求导结果值为零，得到中心偏移量

其中，代表相对插值中心的偏移量，根据SIFT算法在任何一个维度上，即x，y，σ三者中的任何一个数值超过0.5，则视为插值中心已经偏移；

步骤4：将Canny算子和Hillbert-Huang变换融合，具体步骤如下：

(1)求解固有模态函数(Intrinsic Mode Function，IMF)，其中，I(x,y)为原图像，mI(x，y)表示对图像求均值曲面，m^kI(x，y)表示第k次对图像mI(x，y)求均值曲面，原图像与mI(x，y)的差值为第一个IMF分量c₁，得到第一个IMF分量后，将原图像与c₁的差值mI(x，y)作为待处理部分继续分解，则可得到第二个IMF分量c₂，以此类推可得到第n个IMF分量的分解表达式，具体公式如下：

c₁＝I(x，y)-mI(x，y)

c₂＝I(x，y)-c₁-m[I(x，y)-c₁]＝mI(x，y)-m²I(x，y)

c₃＝I(x，y)-c₁-c₂-m[I(x，y)-c₁-c₂]＝m²I(x，y)-m³I(x，y)

c_n＝m^n-1I(x，y)-mⁿI(x，y) (9)

将式(9)中c₁至c_n相加，得到二维IMF：

(2)计算图像的梯度

1)确定图像I(x，y)的x方向偏微分E_x：

2)确定图像I(x，y)的y方向的偏微分E_y：

3)确定图像I(x，y)上点(i，j)的边缘梯度强度A(i，j)：

4)确定图像I(x，y)上点(i，j)的梯度方向α(i，j)：

(3)利用非极大抑制消除错误点并得到单像素边缘点

(4)利用双阈值二值化得到边界

具体实施步骤为：

1)预先设定检测需要的高低阈值，定义低阈值为T₁，高阈值为T₂，且T₂＝2T₁，T₁＝12；

2)对图像进行双阈值处理，对于任意边缘像素值在T₁与T₂之间的，如果能通过边缘连接到一个像素大于T₂而且边缘所有像素大于最小阈值T₁的则保留，否则丢弃，边缘梯度强度A(i，j)大于高阈值则是边缘，边缘梯度强度A(i，j)低于低阈值则不是边缘，边缘梯度强度A(i，j)介于高阈值和低阈值之间的，判断8邻域内是否存在高于高阈值的梯度幅值，存在则是边缘，否则不是边缘；

(5)边界跟踪，得到边缘图像；

1)将边缘梯度强度A(i，j)小于T₁的像素的灰度值设为0，得到图像1；

2)将边缘梯度强度A(i，j)小于T₂的像素的灰度值设为0，得到图像2；

3)对图像2进行扫描，当遇到第一个非零灰度的像素p(x,y)时，以p(x,y)为开始点跟踪轮廓线，直到图像2扫描结束，轮廓线的终点为q(x,y)，此时跟踪结束；

4)找出图像2中q(x,y)点在图像1中相同位置对应的点s(x,y)，当s(x,y)点的8邻近区域中有非零像素s(x,y)存在，则将其包括到图像2中作为r(x,y)点，从r(x,y)开始，重复步骤1)，直到图像1和图像2中均没有非零像素为止，此时完成了对包含p(x,y)的轮廓线的连结；

5)返回到步骤1)，寻找下一条轮廓线，重复步骤1)、2)和3)，直到图像2中找不到新轮廓线，即无非零灰度像素为止，对于Canny算子边缘检测消除边缘响应后得到关键点X＝(x，y，σ)^T，关键点X＝(x，y，σ)^T组成特征点点集R₂；

步骤5：

1)使用SIFT方法得到特征点集R₁，将特征点点集R₁中的点和步骤4中得到的特征点点集R₂所有点两两比较，判断坐标是否相等，相同则舍去R₁中的点，不相同则将R₁中的点与R₂中对应点所在的3×3邻域内的点集R₃比较，相同则舍去R₁中的点，否则，将R₁中的点与R₂中的其他没有比较的边缘点比较，如果任何一次比较有和R₁中参与比较的点坐标相同的点，则舍去R₁中的该点，否则保留；

2)计算标志点f₁和f₂，其中f₁＝R₁-R₂，f₂＝size((R₃-R₁)，1)，f₁和f₂为标志点，f₁是点集R₁与R₂的差集，当f₁为0时，去除点集R₁中的对应点，f₁为1时则保留点集R₁中的对应点作为候选点，当f₂为8时，保留点集R₁中的对应点作为候选点，f₂为7时，则舍去R₁中的该点；

3)将2)中确定的候选点计算其方向，具体步骤如下：

求解关键点梯度的幅值：

求解关键点梯度的方向：

θ(x，y)＝tan^-1((L(x，y+1)-L(x，y)-1))/(L(x+1，y)-L(x-1，y))) (16)

步骤6：

将步骤1至步骤5利用Matlab进行仿真，并利用Matlab中plot()方法在图片中标记出匹配的特征点。

本发明的有益效果是由于采用了Canny算子和Hillbert-Huang变换融合的技术方案，改进了原有SIFT算法实时性不高，特征点有时较少以及对边缘光滑的目标无法准确提取特征点的缺陷，能够有效地增强了SIFT特征点的抗噪声能力，提高了算法的稳定性，为三维重建和目标定位等提供了良好的基础。

附图说明

图1是本发明的融合算法技术路线示意图。

图2是原SIFT算法特征点、本发明特征点以及除去特征点的柱状图。

图3是测试图片一。

图4是测试图片二。

图5是测试图片三。

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。

针对SIFT算法存在的缺陷，SIFT算法中所利用DOG算子进行边缘检测后所得到的关键点，由于边缘响应的原因，有很多并不是真正的关键点。在边缘检测中，Canny算子具有滤波、增强以及检测等多方面优化，受到噪声干扰性小，但是Canny算子处理步骤中，利用高斯算子对图像平滑处理时，会出现过度平滑和边缘移位的情况。二维的Hillbert-Huang变换用于图像分解中，可以将图像分解成对应的内蕴模函数，使边缘轮廓更加清晰。所以将Canny和Hillbert-Huang变换进行融合，进行边缘点检测，再将得到的边缘点与SIFT算法中得到的边缘点比较，两者同时检测出的关键点，即为应该选取的特征点。

本发明的主要技术方案包括以下步骤：

步骤1：对输入的原图像I(x,y)定义尺度空间L(x，y，σ)，构建高斯金字塔：

L(x，y，σ)＝G(x，y，σ)*I(x，y) (1)

其中G(x，y，σ)为高斯函数，m，n为高斯模板的维度，σ为尺度空间因子，值越小表示图像被平滑的越少，相应的尺度越小；

步骤2：构建高斯差分金字塔

(1)将图像金字塔分o组，一组称为一个Octave，每组分为多层，组内总层数为S，每一组图像的首末两层无法进行极值比较，为了满足尺度变化的连续性，在每一组图像的顶层用高斯模糊生成3幅图像，高斯金字塔每组则有S+3层图像；

(2)对原图像进行降采样，采用隔点采样法，具体实施为：比例因子为2对原图像每隔一行一列取一点；

(3)确定高斯金字塔的层数：原图像为金字塔的底部第一层，每次降采样所得到的新图像为金字塔的上一层的一张图像，每个金字塔共n层，金字塔的层数n根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：

n＝log₂{min(M，N)}-t (3)

其中，M，N为原图像的像素值，t为塔顶图像的最小维数的对数值；

(4)计算高斯金字塔组内某一层图像的尺度

其中σ₀是基准层尺度，根据SIFT算法取值为1.6，o为组octave的索引，s为组内层的索引，根据SIFT算法取值为3，关键点的尺度坐标σ即按关键点所在的组和组内的层，不同组相同层的组内尺度坐标相同；

(5)计算DOG算子

D(x，y，σ)＝(G(x，y，kσ)-G(x，y，σ))*I(x，y)＝L(x，y，kσ)-L(x，y，σ) (5)

其中k为组内总层数的倒数，即：

步骤3：对图像进行空间极值点的检测，具体步骤如下：

(1)将步骤2中计算得到的DOG算子D(x，y，σ)泰勒展开，得到

其中X＝(x，y，σ)^T，使公式(5)的泰勒展开式求导结果值为零，得到中心偏移量

其中，代表相对插值中心的偏移量，根据SIFT算法在任何一个维度上，即x，y，σ三者中的任何一个数值超过0.5，则视为插值中心已经偏移；

步骤4：将Canny算子和Hillbert-Huang变换融合，具体步骤如下：

(1)求解固有模态函数(Intrinsic Mode Function，IMF)，其中，I(x,y)为原图像，ml(x，y)表示对图像求均值曲面，m^kI(x，y)表示第k次对图像mI(x，y)求均值曲面，原图像与mI(x，y)的差值为第一个IMF分量c₁,得到第一个IMF分量后，将原图像与c₁的差值mI(x，y)作为待处理部分继续分解，则可得到第二个IMF分量c₂，以此类推可得到第n个IMF分量的分解表达式，具体公式如下：

c₁＝I(x，y)-mI(x，y)

c₂＝I(x，y)-c₁-m[I(x，y)-c₁]＝mI(x，y)-m²I(x，y)

c₃＝I(x，y)-c₁-c₂-m[I(x，y)-c₁-c₂]＝m²I(x，y)-m³I(x，y)

c_n＝m^n-1I(x，y)-mⁿI(x，y) (9)

将式(9)中c₁至c_n相加，得到二维IMF：

(2)计算图像的梯度

1)确定图像I(x，y)的x方向偏微分E_x：

2)确定图像I(x，y)的y方向的偏微分E_y：

3)确定图像I(x，y)上点(i，j)的边缘梯度强度A(i，j)：

4)确定图像I(x，y)上点(i，j)的梯度方向α(i，j)：

(3)利用非极大抑制来消除错误点并得到单像素边缘点

(4)利用双阈值二值化得到边界

具体实施步骤为：

1)预先设定检测需要的高低阈值，定义低阈值为T₁，高阈值为T₂，且T₂＝2T₁，T₁＝12；

2)对图像进行双阈值处理，对于任意边缘像素值在T₁与T₂之间的，如果能通过边缘连接到一个像素大于T₂而且边缘所有像素大于最小阈值T₁的则保留，否则丢弃，边缘梯度强度A(i，j)大于高阈值则是边缘，边缘梯度强度A(i，j)低于低阈值则不是边缘，边缘梯度强度A(i，j)介于高阈值和低阈值之间的，判断8邻域内是否存在高于高阈值的梯度幅值，存在则是边缘，否则不是边缘；

(5)边界跟踪，得到边缘图像；

1)将边缘梯度强度A(i，j)小于T₁的像素的灰度值设为0，得到图像1；

2)将边缘梯度强度A(i，j)小于T₂的像素的灰度值设为0，得到图像2，由于图像2的阈值较高，去除大部分噪音，但同时也损失了有用的边缘信息，而图像1的阈值较低，保留了较多的信息，以图像2为基础，以图像1为补充来连结图像的边缘；

3)对图像2进行扫描，当遇到第一个非零灰度的像素p(x,y)时，以p(x,y)为开始点跟踪轮廓线，直到图像2扫描结束，轮廓线的终点为q(x,y)，此时跟踪结束；

4)找出图像2中q(x,y)点在图像1中相同位置对应的点s(x,y)，当s(x,y)点的8邻近区域中有非零像素s(x,y)存在，则将其包括到图像2中作为r(x,y)点，从r(x,y)开始，重复步骤1)，直到图像1和图像2中均没有非零像素为止，此时完成了对包含p(x,y)的轮廓线的连结；

5)返回到步骤1)，寻找下一条轮廓线，重复步骤1)、2)和3)，直到图像2中找不到新轮廓线，即无非零灰度像素为止，对于Canny算子边缘检测消除边缘响应后得到关键点X＝(x，y，σ)^T，关键点X＝(x，y，σ)^T组成特征点点集R₂；

步骤5：

1)使用SIFT方法得到特征点集R₁，将特征点点集R₁中的点和步骤4中得到的特征点点集R₂所有点两两比较，判断坐标是否相等，相同则舍去R₁中的点，不相同则将R₁中的点与R₂中对应点所在的3×3邻域内的点集R₃比较，相同则舍去R₁中的点，否则，将R₁中的点与R₂中的其他没有比较的边缘点比较，如果任何一次比较有和R₁中参与比较的点坐标相同的点，则舍去R₁中的该点，否则保留；

2)计算标志点f₁和f₂，其中f₁＝R₁-R₂，f₂＝size((R₃-R₁)，1)，f₁和f₂为标志点，f₁是点集R₁与R₂的差集，当f₁为0时，去除点集R₁中的对应点，f₁为1时则保留点集R₁中的对应点作为候选点，当f₂为8时，保留点集R₁中的对应点作为候选点，f₂为7时，则舍R₁中的该点；

3)将2)中确定的候选点计算其方向，具体步骤如下：

求解关键点梯度的幅值：

求解关键点梯度的方向：

θ(x，y)＝tan^-1((L(x，y+1)-L(x，y-1))/(L(x+1，y)-L(x-1，y))) (16)

步骤6：

将步骤1至步骤5利用Matlab进行仿真，并利用Matlab中plot()方法在图片中标记出匹配的特征点，选取图3、图4和图5的三幅图片进行实验，计算得到特征点的数量。

图2是使用三个测试图片得到的匹配特征点的三组柱状图，其中测试一为使用图3的图片的实验结果，测试二为使用图4的图片的实验结果，测试三为使用图5的图片的实验结果，每一组柱状图内的三组数据依次分别为原SIFT算法检测的特征点个数、本发明检测的特征点个数和除去的特征点的个数，所述的除去的特征点为使用本发明与原SIFT算法比较，所除去的特征点的个数，即原SIFT算法与本发明所检测特征点的差值，图2的纵坐标为特征点的个数，从图2能够明显看出本发明比原SIFT算法结果更为精确。

本发明提取的特征点的数目，其中图5的图像较为复杂，如果用原有的SIFT算法进行特征点提取，将会出现很多由于边缘响应产生的错点，本发明不仅鲁棒性增加，更提高了精确度。

Claims (1)

1.一种基于Canny算子和Hillbert-Huang变换的SIFT图像特征点提取的方法，其特征在于包括下述步骤：

步骤1：对输入的原图像I(x,y)定义尺度空间L(x，y，σ)，构建高斯金字塔：

L(x，y，σ)＝G(x，y，σ)*I(x，y) (1)

$G(x,y,\mathrm{\σ})=\frac{1}{2{\mathrm{\π\σ}}^2}{e}^{-\frac{{(x-m/2)}^2+{(y-n/2)}^2}{2{\mathrm{\σ}}^2}}---\left(2\right)$

其中G(x，y，σ)为高斯函数，m，n为高斯模板的维度，σ为尺度空间因子；

步骤2：构建高斯差分金字塔

(1)将图像金字塔分o组，一组称为一个Octave，每组分为多层，组内总层数为S，在每一组图像的顶层用高斯模糊生成3幅图像，高斯金字塔每组则有S+3层图像；

(2)对原图像采用隔点采样法进行降采样，具体实施为：比例因子为2对原图像每隔一行一列取一点；

(3)确定高斯金字塔的层数：原图像为金字塔的底部第一层，每次降采样所得到的新图像为金字塔的上一层的一张图像，每个金字塔共n层，金字塔的层数n根据图像的原始大小和塔顶图像的大小共同决定，其计算公式如下：

n＝log₂{min(M，N)}-t (3)

其中，M，N为原图像的像素值，t为塔顶图像的最小维数的对数值；

(4)计算高斯金字塔组内某一层图像的尺度

$\mathrm{\σ}(o,s)={\mathrm{\σ}}_0{2}^{o+\frac{1}{s}}---\left(4\right)$

其中σ₀是基准层尺度，根据SIFT算法取值为1.6，o为组octave的索引，s为组内层的索引，根据SIFT算法取值为3，关键点的尺度坐标σ即关键点所在的组和组内的层，不同组相同层的组内尺度坐标相同；

(5)计算DOG算子

D(x，y，σ)＝(G(x，y，kσ)-G(x，y，σ))*I(x，y)＝L(x，y，kσ)-L(x，y，σ) (5)

其中k为组内总层数的倒数，即：

$k=2^{\frac{1}{s}}---\left(6\right)$

步骤3：对图像进行空间极值点的检测，具体步骤如下：

将步骤2中计算得到的DOG算子D(x，y，σ)泰勒展开，得到

$D\left(X\right)=D+\frac{\∂D^T}{\∂X}X+\frac{1}{2}{X}^T\frac{{\∂}^{2}D}{\∂X^2}X---\left(7\right)$

其中X＝(x，y，σ)^T，使公式(5)的泰勒展开式求导结果值为零，得到中心偏移量

$\hat{X}=-\frac{{\∂}^2{D}^{-1}}{\∂X^2}\frac{\∂D}{\∂X}---\left(8\right)$

其中，代表相对插值中心的偏移量，根据SIFT算法在任何一个维度上，即x，y，σ三者中的任何一个数值超过0.5，则视为插值中心已经偏移；

步骤4：将Canny算子和Hillbert-Huang变换融合，具体步骤如下：

(1)求解固有模态函数(Intrinsic Mode Function，IMF)，其中，I(x,y)为原图像，mI(x，y)表示对图像求均值曲面，m^kI(x，y)表示第k次对图像mI(x，y)求均值曲面，原图像与mI(x，y)的差值为第一个IMF分量c₁，得到第一个IMF分量后，将原图像与c₁的差值mI(x，y)作为待处理部分继续分解，则可得到第二个IMF分量c₂，以此类推可得到第n个IMF分量的分解表达式，具体公式如下：

c₁＝I(x，y)-mI(x，y)

c₂＝I(x，y)-c₁-m[I(x，y)-c₁]＝mI(x，y)-m²I(x，y)

c₃＝I(x，y)-c₁-c₂-m[I(x，y)-c₁-c₂]＝m²I(x，y)-m³I(x，y)

c_n＝m^n-1I(x，y)-mⁿI(x，y) (9)

将式(9)中c₁至c_n相加，得到二维IMF：

$I(x,y)={\mathrm{\Σ}}_{i=1}^n{c}_n+m^{n}I(x,y)---\left(10\right)$

(2)计算图像的梯度

1)确定图像I(x，y)的x方向偏微分E_x：

$E_x=\frac{\∂G(x,y,\mathrm{\σ})}{\∂x}*I(x,y)---\left(11\right)$

2)确定图像I(x，y)的y方向的偏微分E_y：

$E_y=\frac{\∂D(x,y,\mathrm{\σ})}{\∂x}*I(x,y)---\left(12\right)$

3)确定图像I(x，y)上点(i，j)的边缘梯度强度A(i，j)：

$A(i,j)=\sqrt{{E_x}^2(i,j)+{E_y}^2(i,j)}---\left(13\right)$

4)确定图像I(x，y)上点(i，j)的梯度方向α(i，j)：

$\mathrm{\α}(i,j)=\arctan \frac{E_x(i,j)}{E_y(i,j)}---\left(14\right)$

(3)利用非极大抑制消除错误点并得到单像素边缘点

(4)利用双阈值二值化得到边界

具体实施步骤为：

1)预先设定检测需要的高低阈值，定义低阈值为T₁，高阈值为T₂，且T₂＝2T₁，T₁＝12；

2)对图像进行双阈值处理，对于任意边缘像素值在T₁与T₂之间的，如果能通过边缘连接到一个像素大于T₂而且边缘所有像素大于最小阈值T₁的则保留，否则丢弃，边缘梯度强度A(i，j)大于高阈值则是边缘，边缘梯度强度A(i，j)低于低阈值则不是边缘，边缘梯度强度A(i，j)介于高阈值和低阈值之间的，判断8邻域内是否存在高于高阈值的梯度幅值，存在则是边缘，否则不是边缘；

(5)边界跟踪，得到边缘图像；

1)将边缘梯度强度A(i，j)小于T₁的像素的灰度值设为0，得到图像1；

2)将边缘梯度强度A(i，j)小于T₂的像素的灰度值设为0，得到图像2；

3)对图像2进行扫描，当遇到第一个非零灰度的像素p(x,y)时，以p(x,y)为开始点跟踪轮廓线，直到图像2扫描结束，轮廓线的终点为q(x,y)，此时跟踪结束；

4)找出图像2中q(x,y)点在图像1中相同位置对应的点s(x,y)，当s(x,y)点的8邻近区域中有非零像素s(x,y)存在，则将其包括到图像2中作为r(x,y)点，从r(x,y)开始，重复步骤1)，直到图像1和图像2中均没有非零像素为止，此时完成了对包含p(x,y)的轮廓线的连结；

5)返回到步骤1)，寻找下一条轮廓线，重复步骤1)、2)和3)，直到图像2中找不到新轮廓线，即无非零灰度像素为止，对于Canny算子边缘检测消除边缘响应后得到关键点X＝(x，y，σ)^T，关键点X＝(x，y，σ)^T组成特征点点集R₂；

步骤5：

1)使用SIFT方法得到特征点集R₁，将特征点点集R₁中的点和步骤4中得到的特征点点集R₂所有点两两比较，判断坐标是否相等，相同则舍去R₁中的点，不相同则将R₁中的点与R₂中对应点所在的3×3邻域内的点集R₃比较，相同则舍去R₁中的点，否则，将R₁中的点与R₂中的其他没有比较的边缘点比较，如果任何一次比较有和R₁中参与比较的点坐标相同的点，则舍去R₁中的该点，否则保留；

2)计算标志点f₁和f₂，其中f₁＝R₁-R₂，f₂＝size((R₃-R₁)，1)，f₁和f₂为标志点，f₁是点集R₁与R₂的差集，当f₁为0时，去除点集R₁中的对应点，f₁为1时则保留点集R₁中的对应点作为候选点，当f₂为8时，保留点集R₁中的对应点作为候选点，f₂为7时，则舍去R₁中的该点；

3)将2)中确定的候选点计算其方向，具体步骤如下：

求解关键点梯度的幅值：

$m(x,y)=\sqrt{{\left(t_1\right(x+1,y)-L(x-1,y\left)\right)}^2+{\left(L\right(x,y+1)-L(x,y-1\left)\right)}^2}---\left(15\right)$

求解关键点梯度的方向：

θ(x，y)＝tan^-1((L(x，y+1)-L(x，y-1))/(L(x+1，y)-L(x-1，y))) (16)

步骤6：

将步骤1至步骤5利用Matlab进行仿真，并利用Matlab中plot()方法在图片中标记出匹配的特征点。