[發明專利]基于Canny算子和Hillbert?Huang變換的SIFT圖像特征點提取方法

在審

• 申請號: 201710118809.8
• 申請日: 2017-03-02
• 公開（公告）號: CN106709499A
• 公開（公告）日: 2017-05-24
• 發明（設計）人: 王曉楠;黃登山
• 申請（專利權）人: 西北工業大學
• 主分類號: G06K9/46
• 分類號: G06K9/46;G06T7/13;G06T5/00
• 代理公司: 西北工業大學專利中心61204
• 代理人: 金鳳

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()方法在圖片中標記出匹配的特征點。