教你一步一步用c語(yǔ)言實(shí)現(xiàn)sift算法、上

作者：July、二零一一年三月十二日

出處：http://blog.csdn.net/v_JULY_v

參考：Rob Hess維護(hù)的sift 庫(kù)

環(huán)境：windows xp+vc6.0

條件：c語(yǔ)言實(shí)現(xiàn)。

說(shuō)明：本BLOG內(nèi)會(huì)陸續(xù)一一實(shí)現(xiàn)所有經(jīng)典算法。

------------------------

引言：??? 在我寫(xiě)的關(guān)于sift算法的前倆篇文章里頭，已經(jīng)對(duì)sift算法有了初步的介紹：九、圖像特征提取與匹配之SIFT算法，而后在：九(續(xù))、sift算法的編譯與實(shí)現(xiàn)里，我也簡(jiǎn)單記錄下了如何利用opencv，gsl等庫(kù)編譯運(yùn)行sift程序。

但據(jù)一朋友表示，是否能用c語(yǔ)言實(shí)現(xiàn)sift算法，同時(shí)，盡量不用到opencv，gsl等第三方庫(kù)之類的東西。而且，Rob Hess維護(hù)的sift 庫(kù)，也不好懂，有的人根本搞不懂是怎么一回事。

那么本文，就教你如何利用c語(yǔ)言一步一步實(shí)現(xiàn)sift算法，同時(shí)，你也就能真正明白sift算法到底是怎么一回事了。

ok，先看一下，本程序最終運(yùn)行的效果圖，sift 算法分為五個(gè)步驟(下文詳述)，對(duì)應(yīng)以下第二--第六幅圖：

sift算法的步驟??? 要實(shí)現(xiàn)一個(gè)算法，首先要完全理解這個(gè)算法的原理或思想。咱們先來(lái)簡(jiǎn)單了解下，什么叫sift算法：

sift，尺度不變特征轉(zhuǎn)換，是一種電腦視覺(jué)的算法用來(lái)偵測(cè)與描述影像中的局部性特征，它在空間尺度中尋找極值點(diǎn)，并提取出其位置、尺度、旋轉(zhuǎn)不變量，此算法由 David Lowe 在1999年所發(fā)表，2004年完善總結(jié)。

所謂，Sift算法就是用不同尺度(標(biāo)準(zhǔn)差)的高斯函數(shù)對(duì)圖像進(jìn)行平滑，然后比較平滑后圖像的差別，

差別大的像素就是特征明顯的點(diǎn)。

以下是sift算法的五個(gè)步驟：

一、建立圖像尺度空間(或高斯金字塔)，并檢測(cè)極值點(diǎn)

首先建立尺度空間，要使得圖像具有尺度空間不變形，就要建立尺度空間，sift算法采用了高斯函數(shù)來(lái)建立尺度空間，高斯函數(shù)公式為：

G(x,y,e) = [1/2*pi*e^2] * exp[ -(x^2 + y^2)/2e^2]

上述公式G(x,y,e)，即為尺度可變高斯函數(shù)。

而，一個(gè)圖像的尺度空間L(x,y,e) ,定義為原始圖像I(x,y)與上述的一個(gè)可變尺度的2維高斯函數(shù)G(x,y,e) 卷積運(yùn)算。

即，原始影像I(x，y)在不同的尺度e下，與高斯函數(shù)G(x,y,e)進(jìn)行卷積，得到L(x,y,e)，如下：

L(x,y,e) = G(x,y,e)*I(x,y)

以上的(x，y)是空間坐標(biāo)， e，是尺度坐標(biāo)，或尺度空間因子，e的大小決定平滑程度，大尺度對(duì)應(yīng)圖像的概貌特征，小尺度對(duì)應(yīng)圖像的細(xì)節(jié)特征。大的e值對(duì)應(yīng)粗糙尺度(低分辨率)，反之，對(duì)應(yīng)精細(xì)尺度(高分辨率)。

尺度，受e這個(gè)參數(shù)控制的表示。而不同的L(x,y,e)就構(gòu)成了尺度空間，具體計(jì)算的時(shí)候，即使連續(xù)的高斯函數(shù)，都被離散為(一般為奇數(shù)大小)(2*k+1) *(2*k+1)矩陣，來(lái)和數(shù)字圖像進(jìn)行卷積運(yùn)算。

隨著e的變化，建立起不同的尺度空間，或稱之為建立起圖像的高斯金字塔。

但，像上述L(x,y,e) = G(x,y,e)*I(x,y)的操作，在進(jìn)行高斯卷積時(shí)，整個(gè)圖像就要遍歷所有的像素進(jìn)行卷積(邊界點(diǎn)除外)，于此，就造成了時(shí)間和空間上的很大浪費(fèi)。

為了更有效的在尺度空間檢測(cè)到穩(wěn)定的關(guān)鍵點(diǎn)，也為了縮小時(shí)間和空間復(fù)雜度，對(duì)上述的操作作了一個(gè)改建：即，提出了高斯差分尺度空間(DOG scale-space)。利用不同尺度的高斯差分與原始圖像I(x,y)相乘 ，卷積生成。

D(x,y,e) = ((G(x,y,ke) - G(x,y,e)) * I(x,y)

= L(x,y,ke) - L(x,y,e)??? DOG算子計(jì)算簡(jiǎn)單，是尺度歸一化的LOG算子的近似。

ok，耐心點(diǎn)，咱們?cè)賮?lái)總結(jié)一下上述內(nèi)容：

1、高斯卷積

在組建一組尺度空間后，再組建下一組尺度空間，對(duì)上一組尺度空間的最后一幅圖像進(jìn)行二分之一采樣，得到下一組尺度空間的第一幅圖像，然后進(jìn)行像建立第一組尺度空間那樣的操作，得到第二組尺度空間，公式定義為

L(x,y,e) = G(x,y,e)*I(x,y)

圖像金字塔的構(gòu)建：圖像金字塔共O組，每組有S層，下一組的圖像由上一組圖像降采樣得到，效果圖，圖A如下(左為上一組，右為下一組)：

2、高斯差分

在尺度空間建立完畢后，為了能夠找到穩(wěn)定的關(guān)鍵點(diǎn)，采用高斯差分的方法來(lái)檢測(cè)那些在局部位置的極值點(diǎn)，即采用倆個(gè)相鄰的尺度中的圖像相減，即公式定義為：

D(x,y,e) = ((G(x,y,ke) - G(x,y,e)) * I(x,y)

= L(x,y,ke) - L(x,y,e)

效果圖，圖B：

SIFT的精妙之處在于采用圖像金字塔的方法解決這一問(wèn)題，我們可以把兩幅圖像想象成是連續(xù)的，分別以它們作為底面作四棱錐，就像金字塔，那么每一個(gè) 截面與原圖像相似，那么兩個(gè)金字塔中必然會(huì)有包含大小一致的物體的無(wú)窮個(gè)截面，但應(yīng)用只能是離散的，所以我們只能構(gòu)造有限層，層數(shù)越多當(dāng)然越好，但處理時(shí) 間會(huì)相應(yīng)增加，層數(shù)太少不行，因?yàn)橄蛳虏蓸拥慕孛嬷锌赡苷也坏匠叽绱笮∫恢碌膬蓚€(gè)物體的圖像。

咱們?cè)賮?lái)具體闡述下構(gòu)造D(x,y,e)的詳細(xì)步驟：??? 1、首先采用不同尺度因子的高斯核對(duì)圖像進(jìn)行卷積以得到圖像的不同尺度空間，將這一組圖像作為金子塔圖像的第一層。

2、接著對(duì)第一層圖像中的2倍尺度圖像(相對(duì)于該層第一幅圖像的2倍尺度)以2倍像素距離進(jìn)行下采樣來(lái)得到金子塔圖像的第二層中的第一幅圖像，對(duì)該圖像采用不同尺度因子的高斯核進(jìn)行卷積，以獲得金字塔圖像中第二層的一組圖像。

3、再以金字塔圖像中第二層中的2倍尺度圖像(相對(duì)于該層第一幅圖像的2倍尺度)以2倍像素距離進(jìn)行下采樣來(lái)得到金字塔圖像的第三層中的第一幅圖像，對(duì)該圖像采用不同尺度因子的高斯核進(jìn)行卷積，以獲得金字塔圖像中第三層的一組圖像。這樣依次類推，從而獲得了金字塔圖像的每一層中的一組圖像，如下圖所示：

4、對(duì)上圖得到的每一層相鄰的高斯圖像相減，就得到了高斯差分圖像，如下述第一幅圖所示。下述第二幅圖中的右列顯示了將每組中相鄰圖像相減所生成的高斯差分圖像的結(jié)果，限于篇幅，圖中只給出了第一層和第二層高斯差分圖像的計(jì)算(下述倆幅圖統(tǒng)稱為圖2)：

5、因?yàn)楦咚共罘趾瘮?shù)是歸一化的高斯拉普拉斯函數(shù)的近似，所以可以從高斯差分金字塔分層結(jié)構(gòu)提取出圖像中的極值點(diǎn)作為候選的特征點(diǎn)。對(duì)DOG 尺度空間每個(gè)點(diǎn)與相鄰尺度和相鄰位置的點(diǎn)逐個(gè)進(jìn)行比較，得到的局部極值位置即為特征點(diǎn)所處的位置和對(duì)應(yīng)的尺度。

二、檢測(cè)關(guān)鍵點(diǎn)??? 為了尋找尺度空間的極值點(diǎn)，每一個(gè)采樣點(diǎn)要和它所有的相鄰點(diǎn)比較，看其是否比它的圖像域和尺度域的相鄰點(diǎn)大或者小。如下圖，圖3所示，中間的檢測(cè)點(diǎn)和它同尺度的8個(gè)相鄰點(diǎn)和上下相鄰尺度對(duì)應(yīng)的9×2個(gè)點(diǎn)共26個(gè)點(diǎn)比較，以確保在尺度空間和二維圖像空間都檢測(cè)到極值點(diǎn)。

因?yàn)樾枰噜彸叨冗M(jìn)行比較，所以在一組高斯差分圖像中只能檢測(cè)到兩個(gè)尺度的極值點(diǎn)(如上述第二幅圖中右圖的五角星標(biāo)識(shí))，而其它尺度的極值點(diǎn)檢測(cè)則需要在圖像金字塔的上一層高斯差分圖像中進(jìn)行。依次類推，最終在圖像金字塔中不同層的高斯差分圖像中完成不同尺度極值的檢測(cè)。

當(dāng)然這樣產(chǎn)生的極值點(diǎn)并不都是穩(wěn)定的特征點(diǎn)，因?yàn)槟承O值點(diǎn)響應(yīng)較弱，而且DOG算子會(huì)產(chǎn)生較強(qiáng)的邊緣響應(yīng)。

三、關(guān)鍵點(diǎn)方向的分配??? 為了使描述符具有旋轉(zhuǎn)不變性，需要利用圖像的局部特征為給每一個(gè)關(guān)鍵點(diǎn)分配一個(gè)方向。利用關(guān)鍵點(diǎn)鄰域像素的梯度及方向分布的特性，可以得到梯度模值和方向如下：

其中，尺度為每個(gè)關(guān)鍵點(diǎn)各自所在的尺度。

在以關(guān)鍵點(diǎn)為中心的鄰域窗口內(nèi)采樣，并用直方圖統(tǒng)計(jì)鄰域像素的梯度方向。梯度直方圖的范圍是0～360度，其中每10度一個(gè)方向，總共36個(gè)方向。

直方圖的峰值則代表了該關(guān)鍵點(diǎn)處鄰域梯度的主方向，即作為該關(guān)鍵點(diǎn)的方向。

在計(jì)算方向直方圖時(shí)，需要用一個(gè)參數(shù)等于關(guān)鍵點(diǎn)所在尺度1.5倍的高斯權(quán)重窗對(duì)方向直方圖進(jìn)行加權(quán)，上圖中用藍(lán)色的圓形表示，中心處的藍(lán)色較重，表示權(quán)值最大，邊緣處顏色潛，表示權(quán)值小。如下圖所示，該示例中為了簡(jiǎn)化給出了8方向的方向直方圖計(jì)算結(jié)果，實(shí)際sift創(chuàng)始人David Lowe的原論文中采用36方向的直方圖。

方向直方圖的峰值則代表了該特征點(diǎn)處鄰域梯度的方向，以直方圖中最大值作為該關(guān)鍵點(diǎn)的主方向。為了增強(qiáng)匹配的魯棒性，只保留峰值大于主方向峰值80％的方向作為該關(guān)鍵點(diǎn)的輔方向。因此，對(duì)于同一梯度值的多個(gè)峰值的關(guān)鍵點(diǎn)位置，在相同位置和尺度將會(huì)有多個(gè)關(guān)鍵點(diǎn)被創(chuàng)建但方向不同。僅有15％的關(guān)鍵點(diǎn)被賦予多個(gè)方向，但可以明顯的提高關(guān)鍵點(diǎn)匹配的穩(wěn)定性。

至此，圖像的關(guān)鍵點(diǎn)已檢測(cè)完畢，每個(gè)關(guān)鍵點(diǎn)有三個(gè)信息：位置、所處尺度、方向。由此可以確定一個(gè)SIFT特征區(qū)域。

四、特征點(diǎn)描述符??? 通過(guò)以上步驟，對(duì)于每一個(gè)關(guān)鍵點(diǎn)，擁有三個(gè)信息：位置、尺度以及方向。接下來(lái)就是為每個(gè)關(guān)鍵點(diǎn)建立一個(gè)描述符，使其不隨各種變化而改變，比如光照變化、視角變化等等。并且描述符應(yīng)該有較高的獨(dú)特性，以便于提高特征點(diǎn)正確匹配的概率。

首先將坐標(biāo)軸旋轉(zhuǎn)為關(guān)鍵點(diǎn)的方向，以確保旋轉(zhuǎn)不變性。

接下來(lái)以關(guān)鍵點(diǎn)為中心取8×8的窗口。

上圖，圖5中左部分的中央黑點(diǎn)為當(dāng)前關(guān)鍵點(diǎn)的位置，每個(gè)小格代表關(guān)鍵點(diǎn)鄰域所在尺度空間的一個(gè)像素，箭頭方向代表該像素的梯度方向，箭頭長(zhǎng)度代表梯度模值，圖中藍(lán)色的圈代表高斯加權(quán)的范圍(越靠近關(guān)鍵點(diǎn)的像素梯度方向信息貢獻(xiàn)越大)。

然后在每4×4的小塊上計(jì)算8個(gè)方向的梯度方向直方圖，繪制每個(gè)梯度方向的累加值，即可形成一個(gè)種子點(diǎn)，如圖5右部分所示。此圖中一個(gè)關(guān)鍵點(diǎn)由2×2共4個(gè)種子點(diǎn)組成，每個(gè)種子點(diǎn)有8個(gè)方向向量信息。這種鄰域方向性信息聯(lián)合的思想增強(qiáng)了算法抗噪聲的能力，同時(shí)對(duì)于含有定位誤差的特征匹配也提供了較好的容錯(cuò)性。

實(shí)際計(jì)算過(guò)程中，為了增強(qiáng)匹配的穩(wěn)健性，Lowe建議對(duì)每個(gè)關(guān)鍵點(diǎn)使用4×4共16個(gè)種子點(diǎn)來(lái)描述，這樣對(duì)于一個(gè)關(guān)鍵點(diǎn)就可以產(chǎn)生128個(gè)數(shù)據(jù)，即最終形成128維的SIFT特征向量。此時(shí)SIFT特征向量已經(jīng)去除了尺度變化、旋轉(zhuǎn)等幾何變形因素的影響，再繼續(xù)將特征向量的長(zhǎng)度歸一化，則可以進(jìn)一步去除光照變化的影響。

五、最后一步：當(dāng)兩幅圖像的SIFT特征向量生成后，下一步我們采用關(guān)鍵點(diǎn)特征向量的歐式距離來(lái)作為兩幅圖像中關(guān)鍵點(diǎn)的相似性判定度量。取上圖中，圖像A中的某個(gè)關(guān)鍵點(diǎn)，并找出其與圖像B中歐式距離最近的前兩個(gè)關(guān)鍵點(diǎn)，在這兩個(gè)關(guān)鍵點(diǎn)中，如果最近的距離除以次近的距離少于某個(gè)比例閾值，則接受這一對(duì)匹配點(diǎn)。降低這個(gè)比例閾值，SIFT匹配點(diǎn)數(shù)目會(huì)減少，但更加穩(wěn)定。

sift算法的逐步c實(shí)現(xiàn)??? ok，上文攪了那么多的理論，如果你沒(méi)有看懂它，咋辦列?沒(méi)關(guān)系，下面，咱們來(lái)一步一步實(shí)現(xiàn)此sift算法，即使你沒(méi)有看到上述的理論，慢慢的，你也會(huì)明白sift算法到底是怎么一回事，sift算法到底是怎么實(shí)現(xiàn)的...。

yeah，請(qǐng)看：

前期工作：在具體編寫(xiě)核心函數(shù)之前，得先做幾個(gè)前期的準(zhǔn)備工作：

0、頭文件：

#ifdef?_CH_

#pragma?package?

#endif

#ifndef?_EiC

#include?

#include?"stdlib.h"

#include?"string.h"

#include?"malloc.h"

#include?"math.h"

#include?

#include?

#include?

#include?

#include?

#include?

#include?

#endif

#ifdef?_EiC

#define?WIN32

#endif

1、定義幾個(gè)宏，及變量，以免下文函數(shù)中，突然冒出一個(gè)變量，而您卻不知道怎么一回事：

#define?NUMSIZE?2

#define?GAUSSKERN?3.5

#define?PI?3.14159265358979323846

//Sigma?of?base?image?--?See?D.L.'s?paper.

#define?INITSIGMA?0.5

//Sigma?of?each?octave?--?See?D.L.'s?paper.

#define?SIGMA?sqrt(3)//1.6//

//Number?of?scales?per?octave.??See?D.L.'s?paper.

#define?SCALESPEROCTAVE?2

#define?MAXOCTAVES?4

intnumoctaves;

#define?CONTRAST_THRESHOLD???0.02

#define?CURVATURE_THRESHOLD??10.0

#define?DOUBLE_BASE_IMAGE_SIZE?1

#define?peakRelThresh?0.8

#define?LEN?128

//?temporary?storage

CvMemStorage*?storage?=?0;

2、然后，咱們還得，聲明幾個(gè)變量，以及建幾個(gè)數(shù)據(jù)結(jié)構(gòu)(數(shù)據(jù)結(jié)構(gòu)是一切程序事物的基礎(chǔ)麻，:D。)：

//Data?structure?for?a?float?image.

typedefstructImageSt?{/*金字塔每一層*/

floatlevelsigma;

intlevelsigmalength;

floatabsolute_sigma;

CvMat?*Level;//CvMat是OPENCV的矩陣類，其元素可以是圖像的象素值

}?ImageLevels;

typedefstructImageSt1?{/*金字塔每一階梯*/

introw,?col;//Dimensions?of?image.

floatsubsample;

ImageLevels?*Octave;

}?ImageOctaves;

ImageOctaves?*DOGoctaves;

//DOG?pyr，DOG算子計(jì)算簡(jiǎn)單，是尺度歸一化的LoG算子的近似。

ImageOctaves?*mag_thresh?;

ImageOctaves?*mag_pyr?;

ImageOctaves?*grad_pyr?;

//keypoint數(shù)據(jù)結(jié)構(gòu)，Lists?of?keypoints?are?linked?by?the?"next"?field.

typedefstructKeypointSt

{

floatrow,?col;/*?反饋回原圖像大小，特征點(diǎn)的位置?*/

floatsx,sy;/*?金字塔中特征點(diǎn)的位置*/

intoctave,level;/*金字塔中，特征點(diǎn)所在的階梯、層次*/

floatscale,?ori,mag;/*所在層的尺度sigma,主方向orientation?(range?[-PI,PI])，以及幅值*/

float*descrip;/*特征描述字指針：128維或32維等*/

structKeypointSt?*next;/*?Pointer?to?next?keypoint?in?list.?*/

}?*Keypoint;

//定義特征點(diǎn)具體變量

Keypoint?keypoints=NULL;//用于臨時(shí)存儲(chǔ)特征點(diǎn)的位置等

Keypoint?keyDescriptors=NULL;//用于最后的確定特征點(diǎn)以及特征描述字

3、聲明幾個(gè)圖像的基本處理函數(shù)：

圖像處理基本函數(shù)，其實(shí)也可以用OPENCV的函數(shù)代替，但本文，咱們選擇了用c語(yǔ)言實(shí)現(xiàn)，盡量不用第三方庫(kù)的東西，所以，還得自己編寫(xiě)這些函數(shù)：

CvMat?*?halfSizeImage(CvMat?*?im);//縮小圖像：下采樣

CvMat?*?doubleSizeImage(CvMat?*?im);//擴(kuò)大圖像：最近臨方法

CvMat?*?doubleSizeImage2(CvMat?*?im);//擴(kuò)大圖像：線性插值

floatgetPixelBI(CvMat?*?im,floatcol,floatrow);//雙線性插值函數(shù)

voidnormalizeVec(float*?vec,intdim);//向量歸一化

CvMat*?GaussianKernel2D(floatsigma);//得到2維高斯核

voidnormalizeMat(CvMat*?mat)?;//矩陣歸一化

float*?GaussianKernel1D(floatsigma,intdim)?;//得到1維高斯核

//在具體像素處寬度方向進(jìn)行高斯卷積

floatConvolveLocWidth(float*?kernel,intdim,?CvMat?*?src,intx,inty)?;

//在整個(gè)圖像寬度方向進(jìn)行1D高斯卷積

voidConvolve1DWidth(float*?kern,intdim,?CvMat?*?src,?CvMat?*?dst)?;

//在具體像素處高度方向進(jìn)行高斯卷積

floatConvolveLocHeight(float*?kernel,intdim,?CvMat?*?src,intx,inty)?;

//在整個(gè)圖像高度方向進(jìn)行1D高斯卷積

voidConvolve1DHeight(float*?kern,intdim,?CvMat?*?src,?CvMat?*?dst);

//用高斯函數(shù)模糊圖像

intBlurImage(CvMat?*?src,?CvMat?*?dst,floatsigma)?;

算法核心?? 本程序中，sift算法被分為以下五個(gè)步驟及其相對(duì)應(yīng)的函數(shù)(可能表述與上，或與前倆篇文章有所偏差，但都一個(gè)意思)：

//SIFT算法第一步：圖像預(yù)處理

CvMat?*ScaleInitImage(CvMat?*?im)?;//金字塔初始化

//SIFT算法第二步：建立高斯金字塔函數(shù)

ImageOctaves*?BuildGaussianOctaves(CvMat?*?image)?;//建立高斯金字塔

//SIFT算法第三步：特征點(diǎn)位置檢測(cè)，最后確定特征點(diǎn)的位置

intDetectKeypoint(intnumoctaves,?ImageOctaves?*GaussianPyr);

voidDisplayKeypointLocation(IplImage*?image,?ImageOctaves?*GaussianPyr);

//SIFT算法第四步：計(jì)算高斯圖像的梯度方向和幅值，計(jì)算各個(gè)特征點(diǎn)的主方向

voidComputeGrad_DirecandMag(intnumoctaves,?ImageOctaves?*GaussianPyr);

intFindClosestRotationBin?(intbinCount,floatangle);//進(jìn)行方向直方圖統(tǒng)計(jì)

voidAverageWeakBins?(double*?bins,intbinCount);//對(duì)方向直方圖濾波

//確定真正的主方向

boolInterpolateOrientation?(doubleleft,doublemiddle,doubleright,double*degreeCorrection,double*peakValue);

//確定各個(gè)特征點(diǎn)處的主方向函數(shù)

voidAssignTheMainOrientation(intnumoctaves,?ImageOctaves?*GaussianPyr,ImageOctaves?*mag_pyr,ImageOctaves?*grad_pyr);

//顯示主方向

voidDisplayOrientation?(IplImage*?image,?ImageOctaves?*GaussianPyr);

//SIFT算法第五步：抽取各個(gè)特征點(diǎn)處的特征描述字

voidExtractFeatureDescriptors(intnumoctaves,?ImageOctaves?*GaussianPyr);

//為了顯示圖象金字塔，而作的圖像水平、垂直拼接

CvMat*?MosaicHorizen(?CvMat*?im1,?CvMat*?im2?);

CvMat*?MosaicVertical(?CvMat*?im1,?CvMat*?im2?);

//特征描述點(diǎn)，網(wǎng)格

#define?GridSpacing?4

主體實(shí)現(xiàn)??? ok，以上所有的工作都就緒以后，那么接下來(lái)，咱們就先來(lái)編寫(xiě)main函數(shù)，因?yàn)槟阋豢粗骱瘮?shù)之后，你就立馬能發(fā)現(xiàn)sift算法的工作流程及其原理了。

(主函數(shù)中涉及到的函數(shù)，下一篇文章：一、教你一步一步用c語(yǔ)言實(shí)現(xiàn)sift算法、下，咱們自會(huì)一個(gè)一個(gè)編寫(xiě))：

intmain(void)

{

//聲明當(dāng)前幀IplImage指針

IplImage*?src?=?NULL;

IplImage*?image1?=?NULL;

IplImage*?grey_im1?=?NULL;

IplImage*?DoubleSizeImage?=?NULL;

IplImage*?mosaic1?=?NULL;

IplImage*?mosaic2?=?NULL;

CvMat*?mosaicHorizen1?=?NULL;

CvMat*?mosaicHorizen2?=?NULL;

CvMat*?mosaicVertical1?=?NULL;

CvMat*?image1Mat?=?NULL;

CvMat*?tempMat=NULL;

ImageOctaves?*Gaussianpyr;

introws,cols;

#define?Im1Mat(ROW,COL)?((float?*)(image1Mat->data.fl?+?image1Mat->step/sizeof(float)?*(ROW)))[(COL)]

//灰度圖象像素的數(shù)據(jù)結(jié)構(gòu)

#define?Im1B(ROW,COL)?((uchar*)(image1->imageData?+?image1->widthStep*(ROW)))[(COL)*3]

#define?Im1G(ROW,COL)?((uchar*)(image1->imageData?+?image1->widthStep*(ROW)))[(COL)*3+1]

#define?Im1R(ROW,COL)?((uchar*)(image1->imageData?+?image1->widthStep*(ROW)))[(COL)*3+2]

storage?=?cvCreateMemStorage(0);

//讀取圖片

if(?(src?=?cvLoadImage("street1.jpg",?1))?==?0?)//?test1.jpg?einstein.pgm?back1.bmp

return-1;

//為圖像分配內(nèi)存

image1?=?cvCreateImage(cvSize(src->width,?src->height),??IPL_DEPTH_8U,3);

grey_im1?=?cvCreateImage(cvSize(src->width,?src->height),??IPL_DEPTH_8U,1);

DoubleSizeImage?=?cvCreateImage(cvSize(2*(src->width),?2*(src->height)),??IPL_DEPTH_8U,3);

//為圖像陣列分配內(nèi)存，假設(shè)兩幅圖像的大小相同，tempMat跟隨image1的大小

image1Mat?=?cvCreateMat(src->height,?src->width,?CV_32FC1);

//轉(zhuǎn)化成單通道圖像再處理

cvCvtColor(src,?grey_im1,?CV_BGR2GRAY);

//轉(zhuǎn)換進(jìn)入Mat數(shù)據(jù)結(jié)構(gòu),圖像操作使用的是浮點(diǎn)型操作

cvConvert(grey_im1,?image1Mat);

doublet?=?(double)cvGetTickCount();

//圖像歸一化

cvConvertScale(?image1Mat,?image1Mat,?1.0/255,?0?);

intdim?=?min(image1Mat->rows,?image1Mat->cols);

numoctaves?=?(int)?(log((double)?dim)?/?log(2.0))?-?2;//金字塔階數(shù)

numoctaves?=?min(numoctaves,?MAXOCTAVES);

//SIFT算法第一步，預(yù)濾波除噪聲，建立金字塔底層

tempMat?=?ScaleInitImage(image1Mat)?;

//SIFT算法第二步，建立Guassian金字塔和DOG金字塔

Gaussianpyr?=?BuildGaussianOctaves(tempMat)?;

t?=?(double)cvGetTickCount()?-?t;

printf("the?time?of?build?Gaussian?pyramid?and?DOG?pyramid?is?%.1f/n",?t/(cvGetTickFrequency()*1000.)?);

#define?ImLevels(OCTAVE,LEVEL,ROW,COL)?((float?*)(Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl?+?Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float)?*(ROW)))[(COL)]

//顯示高斯金字塔

for(inti=0;?i

{

if(i==0)

{

mosaicHorizen1=MosaicHorizen(?(Gaussianpyr[0].Octave)[0].Level,?(Gaussianpyr[0].Octave)[1].Level?);

for(intj=2;j

mosaicHorizen1=MosaicHorizen(?mosaicHorizen1,?(Gaussianpyr[0].Octave)[j].Level?);

for(?j=0;j

mosaicHorizen1=halfSizeImage(mosaicHorizen1);

}

elseif(i==1)

{

mosaicHorizen2=MosaicHorizen(?(Gaussianpyr[1].Octave)[0].Level,?(Gaussianpyr[1].Octave)[1].Level?);

for(intj=2;j

mosaicHorizen2=MosaicHorizen(?mosaicHorizen2,?(Gaussianpyr[1].Octave)[j].Level?);

for(?j=0;j

mosaicHorizen2=halfSizeImage(mosaicHorizen2);

mosaicVertical1=MosaicVertical(?mosaicHorizen1,?mosaicHorizen2?);

}

else

{

mosaicHorizen1=MosaicHorizen(?(Gaussianpyr[i].Octave)[0].Level,?(Gaussianpyr[i].Octave)[1].Level?);

for(intj=2;j

mosaicHorizen1=MosaicHorizen(?mosaicHorizen1,?(Gaussianpyr[i].Octave)[j].Level?);

for(?j=0;j

mosaicHorizen1=halfSizeImage(mosaicHorizen1);

mosaicVertical1=MosaicVertical(?mosaicVertical1,?mosaicHorizen1?);

}

}

mosaic1?=?cvCreateImage(cvSize(mosaicVertical1->width,?mosaicVertical1->height),??IPL_DEPTH_8U,1);

cvConvertScale(?mosaicVertical1,?mosaicVertical1,?255.0,?0?);

cvConvertScaleAbs(?mosaicVertical1,?mosaic1,?1,?0?);

//??cvSaveImage("GaussianPyramid?of?me.jpg",mosaic1);

cvNamedWindow("mosaic1",1);

cvShowImage("mosaic1",?mosaic1);

cvWaitKey(0);

cvDestroyWindow("mosaic1");

//顯示DOG金字塔

for(?i=0;?i

{

if(i==0)

{

mosaicHorizen1=MosaicHorizen(?(DOGoctaves[0].Octave)[0].Level,?(DOGoctaves[0].Octave)[1].Level?);

for(intj=2;j

mosaicHorizen1=MosaicHorizen(?mosaicHorizen1,?(DOGoctaves[0].Octave)[j].Level?);

for(?j=0;j

mosaicHorizen1=halfSizeImage(mosaicHorizen1);

}

elseif(i==1)

{

mosaicHorizen2=MosaicHorizen(?(DOGoctaves[1].Octave)[0].Level,?(DOGoctaves[1].Octave)[1].Level?);

for(intj=2;j

mosaicHorizen2=MosaicHorizen(?mosaicHorizen2,?(DOGoctaves[1].Octave)[j].Level?);

for(?j=0;j

mosaicHorizen2=halfSizeImage(mosaicHorizen2);

mosaicVertical1=MosaicVertical(?mosaicHorizen1,?mosaicHorizen2?);

}

else

{

mosaicHorizen1=MosaicHorizen(?(DOGoctaves[i].Octave)[0].Level,?(DOGoctaves[i].Octave)[1].Level?);

for(intj=2;j

mosaicHorizen1=MosaicHorizen(?mosaicHorizen1,?(DOGoctaves[i].Octave)[j].Level?);

for(?j=0;j

mosaicHorizen1=halfSizeImage(mosaicHorizen1);

mosaicVertical1=MosaicVertical(?mosaicVertical1,?mosaicHorizen1?);

}

}

//考慮到DOG金字塔各層圖像都會(huì)有正負(fù)，所以，必須尋找最負(fù)的，以將所有圖像抬高一個(gè)臺(tái)階去顯示

doublemin_val=0;

doublemax_val=0;

cvMinMaxLoc(?mosaicVertical1,?&min_val,?&max_val,NULL,?NULL,?NULL?);

if(?min_val<0.0?)

cvAddS(?mosaicVertical1,?cvScalarAll(?(-1.0)*min_val?),?mosaicVertical1,?NULL?);

mosaic2?=?cvCreateImage(cvSize(mosaicVertical1->width,?mosaicVertical1->height),??IPL_DEPTH_8U,1);

cvConvertScale(?mosaicVertical1,?mosaicVertical1,?255.0/(max_val-min_val),?0?);

cvConvertScaleAbs(?mosaicVertical1,?mosaic2,?1,?0?);

//??cvSaveImage("DOGPyramid?of?me.jpg",mosaic2);

cvNamedWindow("mosaic1",1);

cvShowImage("mosaic1",?mosaic2);

cvWaitKey(0);

//SIFT算法第三步：特征點(diǎn)位置檢測(cè)，最后確定特征點(diǎn)的位置

intkeycount=DetectKeypoint(numoctaves,?Gaussianpyr);

printf("the?keypoints?number?are?%d?;/n",?keycount);

cvCopy(src,image1,NULL);

DisplayKeypointLocation(?image1?,Gaussianpyr);

cvPyrUp(?image1,?DoubleSizeImage,?CV_GAUSSIAN_5x5?);

cvNamedWindow("image1",1);

cvShowImage("image1",?DoubleSizeImage);

cvWaitKey(0);

cvDestroyWindow("image1");

//SIFT算法第四步：計(jì)算高斯圖像的梯度方向和幅值，計(jì)算各個(gè)特征點(diǎn)的主方向

ComputeGrad_DirecandMag(numoctaves,?Gaussianpyr);

AssignTheMainOrientation(?numoctaves,?Gaussianpyr,mag_pyr,grad_pyr);

cvCopy(src,image1,NULL);

DisplayOrientation?(?image1,?Gaussianpyr);

//??cvPyrUp(?image1,?DoubleSizeImage,?CV_GAUSSIAN_5x5?);

cvNamedWindow("image1",1);

//??cvResizeWindow("image1",?2*(image1->width),?2*(image1->height)?);

cvShowImage("image1",?image1);

cvWaitKey(0);

//SIFT算法第五步：抽取各個(gè)特征點(diǎn)處的特征描述字

ExtractFeatureDescriptors(?numoctaves,?Gaussianpyr);

cvWaitKey(0);

//銷毀窗口

cvDestroyWindow("image1");

cvDestroyWindow("mosaic1");

//釋放圖像

cvReleaseImage(&image1);

cvReleaseImage(&grey_im1);

cvReleaseImage(&mosaic1);

cvReleaseImage(&mosaic2);

return0;

}

總結(jié)

以上是生活随笔為你收集整理的C语言实现寻找极值点,九之再续：教你一步一步用c语言实现sift算法、上的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。