[发明专利]基于Canny算子和Hillbert‑Huang变换的SIFT图像特征点提取方法

权利要求书

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()方法在图片中标记出匹配的特征点。