教你一步一步用c語言實現sift演算法

CopperDong發表於2017-11-05
C語言演算法

                     教你一步一步用c語言實現sift演算法、上

作者：July、二零一一年三月十二日
出處：http://blog.csdn.net/v_JULY_v
參考：Rob Hess維護的sift 庫
環境：windows xp+vc6.0
條件：c語言實現。
說明：本BLOG內會陸續一一實現所有經典演算法。
------------------------

引言：
    在我寫的關於sift演算法的前倆篇文章裡頭，已經對sift演算法有了初步的介紹：九、影象特徵提取與匹配之SIFT算法，而後在：九（續）、sift演算法的編譯與實現裡，我也簡單記錄下瞭如何利用opencv，gsl等庫編譯執行sift程式。
    但據一朋友表示，是否能用c語言實現sift演算法，同時，儘量不用到opencv，gsl等第三方庫之類的東西。而且，Rob Hess維護的sift 庫，也不好懂，有的人根本搞不懂是怎麼一回事。
    那麼本文，就教你如何利用c語言一步一步實現sift演算法，同時，你也就能真正明白sift演算法到底是怎麼一回事了。

    ok，先看一下，本程式最終執行的效果圖，sift 演算法分為五個步驟（下文詳述），對應以下第二--第六幅圖：

 

sift演算法的步驟
    要實現一個演算法，首先要完全理解這個演算法的原理或思想。我們們先來簡單瞭解下，什麼叫sift演算法：
    sift，尺度不變特徵轉換，是一種電腦視覺的演算法用來偵測與描述影像中的區域性性特徵，它在空間尺度中尋找極值點，並提取出其位置、尺度、旋轉不變數，此演算法由 David Lowe 在1999年所發表，2004年完善總結。

    所謂，Sift演算法就是用不同尺度（標準差）的高斯函式對影象進行平滑，然後比較平滑後影象的差別，
差別大的畫素就是特徵明顯的點。

    以下是sift演算法的五個步驟：
    一、建立影象尺度空間(或高斯金字塔)，並檢測極值點

    首先建立尺度空間，要使得影象具有尺度空間不變形，就要建立尺度空間，sift演算法採用了高斯函式來建立尺度空間，高斯函式公式為：
        
    上述公式G(x,y,e)，即為尺度可變高斯函式。

    而，一個影象的尺度空間L（x,y,e) ,定義為原始影象I(x,y)與上述的一個可變尺度的2維高斯函式G(x,y,e) 卷積運算。
    即，原始影像I(x，y)在不同的尺度e下，與高斯函式G(x,y,e)進行卷積，得到L(x,y,e)，如下：
        

    以上的（x，y）是空間座標， e，是尺度座標，或尺度空間因子，e的大小決定平滑程度，大尺度對應影象的概貌特徵，小尺度對應影象的細節特徵。大的e值對應粗糙尺度(低解析度)，反之，對應精細尺度(高解析度)。

    尺度，受e這個引數控制的表示。而不同的L（x,y,e）就構成了尺度空間，具體計算的時候，即使連續的高斯函式，都被離散為（一般為奇數大小）(2*k+1) *(2*k+1)矩陣，來和數字影象進行卷積運算。
    隨著e的變化，建立起不同的尺度空間，或稱之為建立起影象的高斯金字塔。

    但，像上述L(x,y,e) = G(x,y,e)*I(x,y)的操作，在進行高斯卷積時，整個影象就要遍歷所有的畫素進行卷積(邊界點除外)，於此，就造成了時間和空間上的很大浪費。
    為了更有效的在尺度空間檢測到穩定的關鍵點，也為了縮小時間和空間複雜度，對上述的操作作了一個改建：即，提出了高斯差分尺度空間（DOG scale-space）。利用不同尺度的高斯差分與原始影象I(x,y)相乘 ，卷積生成。
        
    DOG運算元計算簡單，是尺度歸一化的LOG運算元的近似。

    ok，耐心點，我們們再來總結一下上述內容：
    1、高斯卷積
    在組建一組尺度空間後，再組建下一組尺度空間，對上一組尺度空間的最後一幅影象進行二分之一取樣，得到下一組尺度空間的第一幅影象，然後進行像建立第一組尺度空間那樣的操作，得到第二組尺度空間，公式定義為
         L(x,y,e) = G(x,y,e)*I(x,y)

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

    
    2、高斯差分
    在尺度空間建立完畢後，為了能夠找到穩定的關鍵點，採用高斯差分的方法來檢測那些在區域性位置的極值點，即採用倆個相鄰的尺度中的影象相減，即公式定義為：
        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的精妙之處在於採用影象金字塔的方法解決這一問題，我們可以把兩幅影象想象成是連續的，分別以它們作為底面作四稜錐，就像金字塔，那麼每一個 截面與原影象相似，那麼兩個金字塔中必然會有包含大小一致的物體的無窮個截面，但應用只能是離散的，所以我們只能構造有限層，層數越多當然越好，但處理時 間會相應增加，層數太少不行，因為向下取樣的截面中可能找不到尺寸大小一致的兩個物體的影象。

    我們們再來具體闡述下構造D(x,y,e)的詳細步驟：
    1、首先採用不同尺度因子的高斯核對影象進行卷積以得到影象的不同尺度空間，將這一組影象作為金子塔影象的第一層。
    2、接著對第一層影象中的2倍尺度影象（相對於該層第一幅影象的2倍尺度）以2倍畫素距離進行下采樣來得到金子塔影象的第二層中的第一幅影象，對該影象採用不同尺度因子的高斯核進行卷積，以獲得金字塔影象中第二層的一組影象。
    3、再以金字塔影象中第二層中的2倍尺度影象（相對於該層第一幅影象的2倍尺度）以2倍畫素距離進行下采樣來得到金字塔影象的第三層中的第一幅影象，對該影象採用不同尺度因子的高斯核進行卷積，以獲得金字塔影象中第三層的一組影象。這樣依次類推，從而獲得了金字塔影象的每一層中的一組影象，如下圖所示：

    4、對上圖得到的每一層相鄰的高斯影象相減，就得到了高斯差分影象，如下述第一幅圖所示。下述第二幅圖中的右列顯示了將每組中相鄰影象相減所生成的高斯差分影象的結果，限於篇幅，圖中只給出了第一層和第二層高斯差分影象的計算（下述倆幅圖統稱為圖2）：

影象金字塔的建立：對於一幅影象I,建立其在不同尺度(scale)的影象，也成為子八度（octave），這是為了scale-invariant，也就是在任何尺度都能夠有對應的特徵點，第一個子八度的scale為原圖大小，後面每個octave為上一個octave降取樣的結果，即原圖的1/4（長寬分別減半），構成下一個子八度（高一層金字塔）。

    5、因為高斯差分函式是歸一化的高斯拉普拉斯函式的近似，所以可以從高斯差分金字塔分層結構提取出影象中的極值點作為候選的特徵點。對DOG 尺度空間每個點與相鄰尺度和相鄰位置的點逐個進行比較，得到的區域性極值位置即為特徵點所處的位置和對應的尺度。

    二、檢測關鍵點
    為了尋找尺度空間的極值點，每一個取樣點要和它所有的相鄰點比較，看其是否比它的影象域和尺度域的相鄰點大或者小。如下圖，圖3所示，中間的檢測點和它同尺度的8個相鄰點和上下相鄰尺度對應的9×2個點共26個點比較，以確保在尺度空間和二維影象空間都檢測到極值點。

    因為需要同相鄰尺度進行比較，所以在一組高斯差分影象中只能檢測到兩個尺度的極值點（如上述第二幅圖中右圖的五角星標識），而其它尺度的極值點檢測則需要在影象金字塔的上一層高斯差分影象中進行。依次類推，最終在影象金字塔中不同層的高斯差分影象中完成不同尺度極值的檢測。
    當然這樣產生的極值點並不都是穩定的特徵點，因為某些極值點響應較弱，而且DOG運算元會產生較強的邊緣響應。

    三、關鍵點方向的分配
    為了使描述符具有旋轉不變性，需要利用影象的區域性特徵為給每一個關鍵點分配一個方向。利用關鍵點鄰域畫素的梯度及方向分佈的特性，可以得到梯度模值和方向如下：

   
    其中，尺度為每個關鍵點各自所在的尺度。
    在以關鍵點為中心的鄰域視窗內取樣，並用直方圖統計鄰域畫素的梯度方向。梯度直方圖的範圍是0～360度，其中每10度一個方向，總共36個方向。
直方圖的峰值則代表了該關鍵點處鄰域梯度的主方向，即作為該關鍵點的方向。

    在計算方向直方圖時，需要用一個引數等於關鍵點所在尺度1.5倍的高斯權重窗對方向直方圖進行加權，上圖中用藍色的圓形表示，中心處的藍色較重，表示權值最大，邊緣處顏色潛，表示權值小。如下圖所示，該示例中為了簡化給出了8方向的方向直方圖計算結果，實際sift創始人David Lowe的原論文中採用36方向的直方圖。

    方向直方圖的峰值則代表了該特徵點處鄰域梯度的方向，以直方圖中最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，只保留峰值大於主方向峰值80％的方向作為該關鍵點的輔方向。因此，對於同一梯度值的多個峰值的關鍵點位置，在相同位置和尺度將會有多個關鍵點被建立但方向不同。僅有15％的關鍵點被賦予多個方向，但可以明顯的提高關鍵點匹配的穩定性。
        
    至此，影象的關鍵點已檢測完畢，每個關鍵點有三個資訊：位置、所處尺度、方向。由此可以確定一個SIFT特徵區域。

    四、特徵點描述符
    通過以上步驟，對於每一個關鍵點，擁有三個資訊：位置、尺度以及方向。接下來就是為每個關鍵點建立一個描述符，使其不隨各種變化而改變，比如光照變化、視角變化等等。並且描述符應該有較高的獨特性，以便於提高特徵點正確匹配的概率。  
首先將座標軸旋轉為關鍵點的方向，以確保旋轉不變性。

    接下來以關鍵點為中心取8×8的視窗。
    上圖，圖5中左部分的中央黑點為當前關鍵點的位置，每個小格代表關鍵點鄰域所在尺度空間的一個畫素，箭頭方向代表該畫素的梯度方向，箭頭長度代表梯度模值，圖中藍色的圈代表高斯加權的範圍（越靠近關鍵點的畫素梯度方向資訊貢獻越大）。
    然後在每4×4的小塊上計算8個方向的梯度方向直方圖，繪製每個梯度方向的累加值，即可形成一個種子點，如圖5右部分所示。此圖中一個關鍵點由2×2共4個種子點組成，每個種子點有8個方向向量資訊。這種鄰域方向性資訊聯合的思想增強了演算法抗噪聲的能力，同時對於含有定位誤差的特徵匹配也提供了較好的容錯性。

    實際計算過程中，為了增強匹配的穩健性，Lowe建議對每個關鍵點使用4×4共16個種子點來描述，這樣對於一個關鍵點就可以產生128個資料，即最終形成128維的SIFT特徵向量。此時SIFT特徵向量已經去除了尺度變化、旋轉等幾何變形因素的影響，再繼續將特徵向量的長度歸一化，則可以進一步去除光照變化的影響。


    五、最後一步：當兩幅影象的SIFT特徵向量生成後，下一步我們採用關鍵點特徵向量的歐式距離來作為兩幅影象中關鍵點的相似性判定度量。取上圖中，影象A中的某個關鍵點，並找出其與影象B中歐式距離最近的前兩個關鍵點，在這兩個關鍵點中，如果最近的距離除以次近的距離少於某個比例閾值，則接受這一對匹配點。降低這個比例閾值，SIFT匹配點數目會減少，但更加穩定。關於sift 演算法的更多理論介紹請參看此文：http://blog.csdn.net/abcjennifer/article/details/7639681。

 

sift演算法的逐步c實現
    ok，上文攪了那麼多的理論，如果你沒有看懂它，咋辦列?沒關係，下面，我們們來一步一步實現此sift演算法，即使你沒有看到上述的理論，慢慢的，你也會明白sift演算法到底是怎麼一回事，sift演算法到底是怎麼實現的...。
   yeah，請看：

前期工作：
在具體編寫核心函式之前，得先做幾個前期的準備工作：

0、標頭檔案：

  1. #ifdef _CH_  
  2. #pragma package <opencv>  
  3. #endif  
  4.   
  5. #ifndef _EiC  
  6. #include <stdio.h>  
  7.   
  8. #include "stdlib.h"  
  9. #include "string.h"   
  10. #include "malloc.h"   
  11. #include "math.h"   
  12. #include <assert.h>  
  13. #include <ctype.h>  
  14. #include <time.h>  
  15. #include <cv.h>  
  16. #include <cxcore.h>  
  17. #include <highgui.h>  
  18. #include <vector>  
  19. #endif  
  20.   
  21. #ifdef _EiC  
  22. #define WIN32  
  23. #endif  

1、定義幾個巨集，及變數，以免下文函式中，突然冒出一個變數，而您卻不知道怎麼一回事：

  1. #define NUMSIZE 2  
  2. #define GAUSSKERN 3.5  
  3. #define PI 3.14159265358979323846  
  4.   
  5. //Sigma of base image -- See D.L.'s paper.  
  6. #define INITSIGMA 0.5  
  7. //Sigma of each octave -- See D.L.'s paper.  
  8. #define SIGMA sqrt(3)//1.6//  
  9.   
  10. //Number of scales per octave.  See D.L.'s paper.  
  11. #define SCALESPEROCTAVE 2  
  12. #define MAXOCTAVES 4  
  13. int     numoctaves;  
  14.   
  15. #define CONTRAST_THRESHOLD   0.02  
  16. #define CURVATURE_THRESHOLD  10.0  
  17. #define DOUBLE_BASE_IMAGE_SIZE 1  
  18. #define peakRelThresh 0.8  
  19. #define LEN 128  
  20.   
  21. // temporary storage  
  22. CvMemStorage* storage = 0;   

2、然後，我們們還得，宣告幾個變數，以及建幾個資料結構（資料結構是一切程式事物的基礎麻，:D。）：
  1. //Data structure for a float image.  
  2. typedef struct ImageSt {        /*金字塔每一層*/  
  3.    
  4.  float levelsigma;  
  5.  int levelsigmalength;  
  6.  float absolute_sigma;  
  7.  CvMat *Level;       //CvMat是OPENCV的矩陣類，其元素可以是影象的象素值         
  8. } ImageLevels;  
  9.   
  10. typedef struct ImageSt1 {      /*金字塔每一階梯*/  
  11.  int row, col;          //Dimensions of image.   
  12.  float subsample;  
  13.  ImageLevels *Octave;                
  14. } ImageOctaves;  
  15.   
  16. ImageOctaves *DOGoctaves;        
  17. //DOG pyr，DOG運算元計算簡單，是尺度歸一化的LoG運算元的近似。  
  18.     
  19. ImageOctaves *mag_thresh ;  
  20. ImageOctaves *mag_pyr ;  
  21. ImageOctaves *grad_pyr ;  
  22.   
  23. //keypoint資料結構，Lists of keypoints are linked by the "next" field.  
  24. typedef struct KeypointSt   
  25. {  
  26.  float row, col; /* 反饋回原影象大小，特徵點的位置 */  
  27.  float sx,sy;    /* 金字塔中特徵點的位置*/  
  28.  int octave,level;/*金字塔中，特徵點所在的階梯、層次*/  
  29.    
  30.  float scale, ori,mag; /*所在層的尺度sigma,主方向orientation (range [-PI,PI])，以及幅值*/  
  31.  float *descrip;       /*特徵描述字指標：128維或32維等*/  
  32.  struct KeypointSt *next;/* Pointer to next keypoint in list. */  
  33. } *Keypoint;  
  34.   
  35. //定義特徵點具體變數  
  36. Keypoint keypoints=NULL;      //用於臨時儲存特徵點的位置等  
  37. Keypoint keyDescriptors=NULL; //用於最後的確定特徵點以及特徵描述字  

3、宣告幾個影象的基本處理函式：

  1. CvMat * halfSizeImage(CvMat * im);     //縮小影象：下采樣  
  2. CvMat * doubleSizeImage(CvMat * im);   //擴大影象：最近臨方法  
  3. CvMat * doubleSizeImage2(CvMat * im);  //擴大影象：線性插值  
  4. float getPixelBI(CvMat * im, float col, float row);//雙線性插值函式  
  5. void normalizeVec(float* vec, int dim);//向量歸一化    
  6. CvMat* GaussianKernel2D(float sigma);  //得到2維高斯核  
  7. void normalizeMat(CvMat* mat) ;        //矩陣歸一化  
  8. float* GaussianKernel1D(float sigma, int dim) ; //得到1維高斯核  
  9.   
  10. //在具體畫素處寬度方向進行高斯卷積  
  11. float ConvolveLocWidth(float* kernel, int dim, CvMat * src, int x, int y) ;    
  12. //在整個影象寬度方向進行1D高斯卷積  
  13. void Convolve1DWidth(float* kern, int dim, CvMat * src, CvMat * dst) ;         
  14. //在具體畫素處高度方向進行高斯卷積  
  15. float ConvolveLocHeight(float* kernel, int dim, CvMat * src, int x, int y) ;   
  16. //在整個影象高度方向進行1D高斯卷積  
  17. void Convolve1DHeight(float* kern, int dim, CvMat * src, CvMat * dst);       
  18. //用高斯函式模糊影象    
  19. int BlurImage(CvMat * src, CvMat * dst, float sigma) ;            
             
演算法核心
   本程式中，sift演算法被分為以下五個步驟及其相對應的函式（可能表述與上，或與前倆篇文章有所偏差，但都一個意思）：
  1. //SIFT演算法第一步：影象預處理  
  2. CvMat *ScaleInitImage(CvMat * im) ;                  //金字塔初始化  
  3.   
  4. //SIFT演算法第二步：建立高斯金字塔函式  
  5. ImageOctaves* BuildGaussianOctaves(CvMat * image) ;  //建立高斯金字塔  
  6.   
  7. //SIFT演算法第三步：特徵點位置檢測，最後確定特徵點的位置  
  8. int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr);  
  9. void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr);  
  10.   
  11. //SIFT演算法第四步：計算高斯影象的梯度方向和幅值，計算各個特徵點的主方向  
  12. void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr);  
  13.   
  14. int FindClosestRotationBin (int binCount, float angle);  //進行方向直方圖統計  
  15. void AverageWeakBins (double* bins, int binCount);       //對方向直方圖濾波  
  16. //確定真正的主方向  
  17. bool InterpolateOrientation (double left, double middle,double right, double *degreeCorrection, double *peakValue);  
  18. //確定各個特徵點處的主方向函式  
  19. void AssignTheMainOrientation(int numoctaves, ImageOctaves *GaussianPyr,ImageOctaves *mag_pyr,ImageOctaves *grad_pyr);  
  20. //顯示主方向  
  21. void DisplayOrientation (IplImage* image, ImageOctaves *GaussianPyr);  
  22.   
  23. //SIFT演算法第五步：抽取各個特徵點處的特徵描述字  
  24. void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr);  
  25.   
  26. //為了顯示圖象金字塔，而作的影象水平、垂直拼接  
  27. CvMat* MosaicHorizen( CvMat* im1, CvMat* im2 );  
  28. CvMat* MosaicVertical( CvMat* im1, CvMat* im2 );  
  29.   
  30. //特徵描述點，網格    
  31. #define GridSpacing 4  


主體實現
    ok，以上所有的工作都就緒以後，那麼接下來，我們們就先來編寫main函式，因為你一看主函式之後，你就立馬能發現sift演算法的工作流程及其原理了。
    （主函式中涉及到的函式，下一篇文章：一、教你一步一步用c語言實現sift演算法、下，我們們自會一個一個編寫）：

  1. int main( void )  
  2. {  
  3.  //宣告當前幀IplImage指標  
  4.  IplImage* src = NULL;   
  5.  IplImage* image1 = NULL;   
  6.  IplImage* grey_im1 = NULL;   
  7.  IplImage* DoubleSizeImage = NULL;  
  8.    
  9.  IplImage* mosaic1 = NULL;   
  10.  IplImage* mosaic2 = NULL;   
  11.    
  12.  CvMat* mosaicHorizen1 = NULL;  
  13.  CvMat* mosaicHorizen2 = NULL;  
  14.  CvMat* mosaicVertical1 = NULL;  
  15.    
  16.  CvMat* image1Mat = NULL;  
  17.  CvMat* tempMat=NULL;  
  18.    
  19.  ImageOctaves *Gaussianpyr;  
  20.  int rows,cols;  
  21.   
  22. #define Im1Mat(ROW,COL) ((float *)(image1Mat->data.fl + image1Mat->step/sizeof(float) *(ROW)))[(COL)]  
  23.    
  24.  //灰度圖象畫素的資料結構  
  25. #define Im1B(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3]  
  26. #define Im1G(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+1]  
  27. #define Im1R(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+2]  
  28.    
  29.  storage = cvCreateMemStorage(0);   
  30.    
  31.  //讀取圖片  
  32.  if( (src = cvLoadImage( "street1.jpg", 1)) == 0 )  // test1.jpg einstein.pgm back1.bmp  
  33.   return -1;  
  34.   
  35.  //為影象分配記憶體   
  36.  image1 = cvCreateImage(cvSize(src->width, src->height),  IPL_DEPTH_8U,3);  
  37.  grey_im1 = cvCreateImage(cvSize(src->width, src->height),  IPL_DEPTH_8U,1);  
  38.  DoubleSizeImage = cvCreateImage(cvSize(2*(src->width), 2*(src->height)),  IPL_DEPTH_8U,3);  
  39.   
  40.  //為影象陣列分配記憶體，假設兩幅影象的大小相同，tempMat跟隨image1的大小  
  41.  image1Mat = cvCreateMat(src->height, src->width, CV_32FC1);  
  42.  //轉化成單通道影象再處理  
  43.  cvCvtColor(src, grey_im1, CV_BGR2GRAY);  
  44.  //轉換進入Mat資料結構,影象操作使用的是浮點型操作  
  45.  cvConvert(grey_im1, image1Mat);  
  46.    
  47.  double t = (double)cvGetTickCount();  
  48.  //影象歸一化  
  49.  cvConvertScale( image1Mat, image1Mat, 1.0/255, 0 );  
  50.    
  51.  int dim = min(image1Mat->rows, image1Mat->cols);  
  52.  numoctaves = (int) (log((double) dim) / log(2.0)) - 2;    //金字塔階數  
  53.  numoctaves = min(numoctaves, MAXOCTAVES);  
  54.   
  55.  //SIFT演算法第一步，預濾波除噪聲，建立金字塔底層  
  56.  tempMat = ScaleInitImage(image1Mat) ;  
  57.  //SIFT演算法第二步，建立Guassian金字塔和DOG金字塔  
  58.  Gaussianpyr = BuildGaussianOctaves(tempMat) ;  
  59.    
  60.  t = (double)cvGetTickCount() - t;  
  61.  printf( "the time of build Gaussian pyramid and DOG pyramid is %.1f/n", t/(cvGetTickFrequency()*1000.) );  
  62.    
  63. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
  64.  //顯示高斯金字塔  
  65.  for (int i=0; i<numoctaves;i++)  
  66.  {  
  67.   if (i==0)  
  68.   {  
  69.    mosaicHorizen1=MosaicHorizen( (Gaussianpyr[0].Octave)[0].Level, (Gaussianpyr[0].Octave)[1].Level );  
  70.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
  71.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[0].Octave)[j].Level );  
  72.    for ( j=0;j<NUMSIZE;j++)  
  73.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
  74.   }  
  75.   else if (i==1)  
  76.   {  
  77.    mosaicHorizen2=MosaicHorizen( (Gaussianpyr[1].Octave)[0].Level, (Gaussianpyr[1].Octave)[1].Level );  
  78.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
  79.     mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (Gaussianpyr[1].Octave)[j].Level );  
  80.    for ( j=0;j<NUMSIZE;j++)  
  81.     mosaicHorizen2=halfSizeImage(mosaicHorizen2);  
  82.    mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );  
  83.   }  
  84.   else  
  85.   {  
  86.    mosaicHorizen1=MosaicHorizen( (Gaussianpyr[i].Octave)[0].Level, (Gaussianpyr[i].Octave)[1].Level );  
  87.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
  88.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[i].Octave)[j].Level );  
  89.    for ( j=0;j<NUMSIZE;j++)  
  90.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
  91.    mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );  
  92.   }  
  93.  }  
  94.  mosaic1 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height),  IPL_DEPTH_8U,1);  
  95.  cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0, 0 );  
  96.  cvConvertScaleAbs( mosaicVertical1, mosaic1, 1, 0 );  
  97.    
  98.  //  cvSaveImage("GaussianPyramid of me.jpg",mosaic1);  
  99.  cvNamedWindow("mosaic1",1);  
  100.  cvShowImage("mosaic1", mosaic1);  
  101.  cvWaitKey(0);  
  102.  cvDestroyWindow("mosaic1");  
  103.  //顯示DOG金字塔  
  104.  for ( i=0; i<numoctaves;i++)  
  105.  {  
  106.   if (i==0)  
  107.   {  
  108.    mosaicHorizen1=MosaicHorizen( (DOGoctaves[0].Octave)[0].Level, (DOGoctaves[0].Octave)[1].Level );  
  109.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
  110.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[0].Octave)[j].Level );  
  111.    for ( j=0;j<NUMSIZE;j++)  
  112.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
  113.   }  
  114.   else if (i==1)  
  115.   {  
  116.    mosaicHorizen2=MosaicHorizen( (DOGoctaves[1].Octave)[0].Level, (DOGoctaves[1].Octave)[1].Level );  
  117.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
  118.     mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (DOGoctaves[1].Octave)[j].Level );  
  119.    for ( j=0;j<NUMSIZE;j++)  
  120.     mosaicHorizen2=halfSizeImage(mosaicHorizen2);  
  121.    mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );  
  122.   }  
  123.   else  
  124.   {  
  125.    mosaicHorizen1=MosaicHorizen( (DOGoctaves[i].Octave)[0].Level, (DOGoctaves[i].Octave)[1].Level );  
  126.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
  127.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[i].Octave)[j].Level );  
  128.    for ( j=0;j<NUMSIZE;j++)  
  129.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
  130.    mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );  
  131.   }  
  132.  }  
  133.  //考慮到DOG金字塔各層影象都會有正負，所以，必須尋找最負的，以將所有影象抬高一個臺階去顯示  
  134.  double min_val=0;  
  135.  double max_val=0;  
  136.  cvMinMaxLoc( mosaicVertical1, &min_val, &max_val,NULL, NULL, NULL );  
  137.  if ( min_val<0.0 )  
  138.   cvAddS( mosaicVertical1, cvScalarAll( (-1.0)*min_val ), mosaicVertical1, NULL );  
  139.  mosaic2 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height),  IPL_DEPTH_8U,1);  
  140.  cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0/(max_val-min_val), 0 );  
  141.  cvConvertScaleAbs( mosaicVertical1, mosaic2, 1, 0 );  
  142.    
  143.  //  cvSaveImage("DOGPyramid of me.jpg",mosaic2);  
  144.  cvNamedWindow("mosaic1",1);  
  145.  cvShowImage("mosaic1", mosaic2);  
  146.  cvWaitKey(0);  
  147.    
  148.  //SIFT演算法第三步：特徵點位置檢測，最後確定特徵點的位置  
  149.  int keycount=DetectKeypoint(numoctaves, Gaussianpyr);  
  150.  printf("the keypoints number are %d ;/n", keycount);  
  151.  cvCopy(src,image1,NULL);  
  152.  DisplayKeypointLocation( image1 ,Gaussianpyr);  
  153.    
  154.  cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
  155.  cvNamedWindow("image1",1);  
  156.  cvShowImage("image1", DoubleSizeImage);  
  157.  cvWaitKey(0);    
  158.  cvDestroyWindow("image1");  
  159.    
  160.  //SIFT演算法第四步：計算高斯影象的梯度方向和幅值，計算各個特徵點的主方向  
  161.  ComputeGrad_DirecandMag(numoctaves, Gaussianpyr);  
  162.  AssignTheMainOrientation( numoctaves, Gaussianpyr,mag_pyr,grad_pyr);  
  163.  cvCopy(src,image1,NULL);  
  164.  DisplayOrientation ( image1, Gaussianpyr);  
  165.    
  166.  //  cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
  167.  cvNamedWindow("image1",1);  
  168.  //  cvResizeWindow("image1", 2*(image1->width), 2*(image1->height) );  
  169.  cvShowImage("image1", image1);  
  170.  cvWaitKey(0);  
  171.    
  172.  //SIFT演算法第五步：抽取各個特徵點處的特徵描述字  
  173.  ExtractFeatureDescriptors( numoctaves, Gaussianpyr);  
  174.  cvWaitKey(0);  
  175.    
  176.  //銷燬視窗  
  177.  cvDestroyWindow("image1");  
  178.  cvDestroyWindow("mosaic1");  
  179.  //釋放影象  
  180.  cvReleaseImage(&image1);  
  181.  cvReleaseImage(&grey_im1);  
  182.  cvReleaseImage(&mosaic1);  
  183.  cvReleaseImage(&mosaic2);  
  184.  return 0;  
  185. }  

更多見下文：一、教你一步一步用c語言實現sift演算法、下。本文完。

   教你一步一步用c語言實現sift演算法、下


作者：July、二零一一年三月十二日
出處：http://blog.csdn.net/v_JULY_v。
參考：Rob Hess維護的sift 庫
環境：windows xp+vc6.0
條件：c語言實現。
說明：本BLOG內會陸續一一實現所有經典演算法。
------------------------


本文接上，教你一步一步用c語言實現sift演算法、上，而來：
函式編寫
    ok，接上文，我們們一個一個的來編寫main函式中所涉及到所有函式，這也是本文的關鍵部分：

  1. //下采樣原來的影象，返回縮小2倍尺寸的影象  
  2. CvMat * halfSizeImage(CvMat * im)   
  3. {  
  4.  unsigned int i,j;  
  5.  int w = im->cols/2;  
  6.  int h = im->rows/2;   
  7.  CvMat *imnew = cvCreateMat(h, w, CV_32FC1);  
  8.    
  9. #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]  
  10. #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]  
  11.  for ( j = 0; j < h; j++)   
  12.   for ( i = 0; i < w; i++)   
  13.    Imnew(j,i)=Im(j*2, i*2);  
  14.   return imnew;  
  15. }  
  16.   
  17. //上取樣原來的影象，返回放大2倍尺寸的影象  
  18. CvMat * doubleSizeImage(CvMat * im)   
  19. {  
  20.  unsigned int i,j;  
  21.  int w = im->cols*2;  
  22.  int h = im->rows*2;   
  23.  CvMat *imnew = cvCreateMat(h, w, CV_32FC1);  
  24.    
  25. #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]  
  26. #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]  
  27.    
  28.  for ( j = 0; j < h; j++)   
  29.   for ( i = 0; i < w; i++)   
  30.    Imnew(j,i)=Im(j/2, i/2);  
  31.     
  32.   return imnew;  
  33. }  
  34.   
  35. //上取樣原來的影象，返回放大2倍尺寸的線性插值影象  
  36. CvMat * doubleSizeImage2(CvMat * im)   
  37. {  
  38.  unsigned int i,j;  
  39.  int w = im->cols*2;  
  40.  int h = im->rows*2;   
  41.  CvMat *imnew = cvCreateMat(h, w, CV_32FC1);  
  42.    
  43. #define Im(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]  
  44. #define Imnew(ROW,COL) ((float *)(imnew->data.fl + imnew->step/sizeof(float) *(ROW)))[(COL)]  
  45.    
  46.  // fill every pixel so we don't have to worry about skipping pixels later  
  47.  for ( j = 0; j < h; j++)   
  48.  {  
  49.   for ( i = 0; i < w; i++)   
  50.   {  
  51.    Imnew(j,i)=Im(j/2, i/2);  
  52.   }  
  53.  }  
  54.  /* 
  55.  A B C 
  56.  E F G 
  57.  H I J 
  58.  pixels A C H J are pixels from original image 
  59.  pixels B E G I F are interpolated pixels 
  60.  */  
  61.  // interpolate pixels B and I  
  62.  for ( j = 0; j < h; j += 2)  
  63.   for ( i = 1; i < w - 1; i += 2)  
  64.    Imnew(j,i)=0.5*(Im(j/2, i/2)+Im(j/2, i/2+1));  
  65.   // interpolate pixels E and G  
  66.   for ( j = 1; j < h - 1; j += 2)  
  67.    for ( i = 0; i < w; i += 2)  
  68.     Imnew(j,i)=0.5*(Im(j/2, i/2)+Im(j/2+1, i/2));  
  69.    // interpolate pixel F  
  70.    for ( j = 1; j < h - 1; j += 2)  
  71.     for ( i = 1; i < w - 1; i += 2)  
  72.      Imnew(j,i)=0.25*(Im(j/2, i/2)+Im(j/2+1, i/2)+Im(j/2, i/2+1)+Im(j/2+1, i/2+1));  
  73.     return imnew;  
  74. }  
  75.   
  76. //雙線性插值，返回畫素間的灰度值  
  77. float getPixelBI(CvMat * im, float col, float row)   
  78. {  
  79.  int irow, icol;  
  80.  float rfrac, cfrac;  
  81.  float row1 = 0, row2 = 0;  
  82.  int width=im->cols;  
  83.  int height=im->rows;  
  84. #define ImMat(ROW,COL) ((float *)(im->data.fl + im->step/sizeof(float) *(ROW)))[(COL)]  
  85.    
  86.  irow = (int) row;  
  87.  icol = (int) col;  
  88.    
  89.  if (irow < 0 || irow >= height  
  90.   || icol < 0 || icol >= width)  
  91.   return 0;  
  92.  if (row > height - 1)  
  93.   row = height - 1;  
  94.  if (col > width - 1)  
  95.   col = width - 1;  
  96.  rfrac = 1.0 - (row - (float) irow);  
  97.  cfrac = 1.0 - (col - (float) icol);  
  98.  if (cfrac < 1)   
  99.  {  
  100.   row1 = cfrac * ImMat(irow,icol) + (1.0 - cfrac) * ImMat(irow,icol+1);  
  101.  }   
  102.  else   
  103.  {  
  104.   row1 = ImMat(irow,icol);  
  105.  }  
  106.  if (rfrac < 1)   
  107.  {  
  108.   if (cfrac < 1)   
  109.   {  
  110.    row2 = cfrac * ImMat(irow+1,icol) + (1.0 - cfrac) * ImMat(irow+1,icol+1);  
  111.   } else   
  112.   {  
  113.    row2 = ImMat(irow+1,icol);  
  114.   }  
  115.  }  
  116.  return rfrac * row1 + (1.0 - rfrac) * row2;  
  117. }  
  118.   
  119. //矩陣歸一化  
  120. void normalizeMat(CvMat* mat)   
  121. {  
  122. #define Mat(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]  
  123.  float sum = 0;  
  124.    
  125.  for (unsigned int j = 0; j < mat->rows; j++)   
  126.   for (unsigned int i = 0; i < mat->cols; i++)   
  127.    sum += Mat(j,i);  
  128.   for ( j = 0; j < mat->rows; j++)   
  129.    for (unsigned int i = 0; i < mat->rows; i++)   
  130.     Mat(j,i) /= sum;  
  131. }  
  132.   
  133. //向量歸一化  
  134. void normalizeVec(float* vec, int dim)   
  135. {  
  136.     unsigned int i;  
  137.  float sum = 0;  
  138.  for ( i = 0; i < dim; i++)  
  139.   sum += vec[i];  
  140.  for ( i = 0; i < dim; i++)  
  141.   vec[i] /= sum;  
  142. }  
  143.   
  144. //得到向量的歐式長度，2-範數  
  145. float GetVecNorm( float* vec, int dim )  
  146. {  
  147.  float sum=0.0;  
  148.  for (unsigned int i=0;i<dim;i++)  
  149.   sum+=vec[i]*vec[i];  
  150.     return sqrt(sum);  
  151. }  
  152.   
  153. //產生1D高斯核  
  154. float* GaussianKernel1D(float sigma, int dim)   
  155. {  
  156.    
  157.  unsigned int i;  
  158.  //printf("GaussianKernel1D(): Creating 1x%d vector for sigma=%.3f gaussian kernel/n", dim, sigma);  
  159.    
  160.  float *kern=(float*)malloc( dim*sizeof(float) );  
  161.  float s2 = sigma * sigma;  
  162.  int c = dim / 2;  
  163.  float m= 1.0/(sqrt(2.0 * CV_PI) * sigma);  
  164.     double v;   
  165.  for ( i = 0; i < (dim + 1) / 2; i++)   
  166.  {  
  167.   v = m * exp(-(1.0*i*i)/(2.0 * s2)) ;  
  168.   kern[c+i] = v;  
  169.   kern[c-i] = v;  
  170.  }  
  171.  //   normalizeVec(kern, dim);  
  172.  // for ( i = 0; i < dim; i++)  
  173.  //  printf("%f  ", kern[i]);  
  174.  //  printf("/n");  
  175.  return kern;  
  176. }  
  177.   
  178. //產生2D高斯核矩陣  
  179. CvMat* GaussianKernel2D(float sigma)   
  180. {  
  181.  // int dim = (int) max(3.0f, GAUSSKERN * sigma);  
  182.     int dim = (int) max(3.0f, 2.0 * GAUSSKERN *sigma + 1.0f);  
  183.  // make dim odd  
  184.  if (dim % 2 == 0)  
  185.   dim++;  
  186.  //printf("GaussianKernel(): Creating %dx%d matrix for sigma=%.3f gaussian/n", dim, dim, sigma);  
  187.  CvMat* mat=cvCreateMat(dim, dim, CV_32FC1);  
  188. #define Mat(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]  
  189.  float s2 = sigma * sigma;  
  190.  int c = dim / 2;  
  191.  //printf("%d %d/n", mat.size(), mat[0].size());  
  192.     float m= 1.0/(sqrt(2.0 * CV_PI) * sigma);  
  193.  for (int i = 0; i < (dim + 1) / 2; i++)   
  194.  {  
  195.   for (int j = 0; j < (dim + 1) / 2; j++)   
  196.   {  
  197.    //printf("%d %d %d/n", c, i, j);  
  198.    float v = m * exp(-(1.0*i*i + 1.0*j*j) / (2.0 * s2));  
  199.    Mat(c+i,c+j) =v;  
  200.    Mat(c-i,c+j) =v;  
  201.    Mat(c+i,c-j) =v;  
  202.    Mat(c-i,c-j) =v;  
  203.   }  
  204.  }  
  205.  // normalizeMat(mat);  
  206.  return mat;  
  207. }  
  208.   
  209. //x方向畫素處作卷積  
  210. float ConvolveLocWidth(float* kernel, int dim, CvMat * src, int x, int y)   
  211. {  
  212. #define Src(ROW,COL) ((float *)(src->data.fl + src->step/sizeof(float) *(ROW)))[(COL)]  
  213.     unsigned int i;  
  214.  float pixel = 0;  
  215.     int col;  
  216.  int cen = dim / 2;  
  217.  //printf("ConvolveLoc(): Applying convoluation at location (%d, %d)/n", x, y);  
  218.  for ( i = 0; i < dim; i++)   
  219.  {  
  220.   col = x + (i - cen);  
  221.   if (col < 0)  
  222.    col = 0;  
  223.   if (col >= src->cols)  
  224.    col = src->cols - 1;  
  225.   pixel += kernel[i] * Src(y,col);  
  226.  }  
  227.  if (pixel > 1)  
  228.   pixel = 1;  
  229.  return pixel;  
  230. }  
  231.   
  232. //x方向作卷積  
  233. void Convolve1DWidth(float* kern, int dim, CvMat * src, CvMat * dst)   
  234. {  
  235. #define DST(ROW,COL) ((float *)(dst->data.fl + dst->step/sizeof(float) *(ROW)))[(COL)]  
  236.  unsigned int i,j;  
  237.    
  238.  for ( j = 0; j < src->rows; j++)   
  239.  {  
  240.   for ( i = 0; i < src->cols; i++)   
  241.   {  
  242.    //printf("%d, %d/n", i, j);  
  243.    DST(j,i) = ConvolveLocWidth(kern, dim, src, i, j);  
  244.   }  
  245.  }  
  246. }  
  247.   
  248. //y方向畫素處作卷積  
  249. float ConvolveLocHeight(float* kernel, int dim, CvMat * src, int x, int y)   
  250. {  
  251. #define Src(ROW,COL) ((float *)(src->data.fl + src->step/sizeof(float) *(ROW)))[(COL)]  
  252.     unsigned int j;  
  253.  float pixel = 0;  
  254.  int cen = dim / 2;  
  255.  //printf("ConvolveLoc(): Applying convoluation at location (%d, %d)/n", x, y);  
  256.  for ( j = 0; j < dim; j++)   
  257.  {  
  258.   int row = y + (j - cen);  
  259.   if (row < 0)  
  260.    row = 0;  
  261.   if (row >= src->rows)  
  262.    row = src->rows - 1;  
  263.   pixel += kernel[j] * Src(row,x);  
  264.  }  
  265.  if (pixel > 1)  
  266.   pixel = 1;  
  267.  return pixel;  
  268. }  
  269.   
  270. //y方向作卷積  
  271. void Convolve1DHeight(float* kern, int dim, CvMat * src, CvMat * dst)   
  272. {  
  273. #define Dst(ROW,COL) ((float *)(dst->data.fl + dst->step/sizeof(float) *(ROW)))[(COL)]  
  274.     unsigned int i,j;  
  275.  for ( j = 0; j < src->rows; j++)   
  276.  {  
  277.   for ( i = 0; i < src->cols; i++)   
  278.   {  
  279.    //printf("%d, %d/n", i, j);  
  280.    Dst(j,i) = ConvolveLocHeight(kern, dim, src, i, j);  
  281.   }  
  282.  }  
  283. }  
  284.   
  285. //卷積模糊影象  
  286. int BlurImage(CvMat * src, CvMat * dst, float sigma)   
  287. {  
  288.  float* convkernel;  
  289.  int dim = (int) max(3.0f, 2.0 * GAUSSKERN * sigma + 1.0f);  
  290.     CvMat *tempMat;  
  291.  // make dim odd  
  292.  if (dim % 2 == 0)  
  293.   dim++;  
  294.  tempMat = cvCreateMat(src->rows, src->cols, CV_32FC1);  
  295.  convkernel = GaussianKernel1D(sigma, dim);  
  296.    
  297.  Convolve1DWidth(convkernel, dim, src, tempMat);  
  298.  Convolve1DHeight(convkernel, dim, tempMat, dst);  
  299.  cvReleaseMat(&tempMat);  
  300.  return dim;  
  301. }  

五個步驟

    ok，接下來，進入重點部分，我們們依據上文介紹的sift演算法的幾個步驟，來一一實現這些函式。
    為了版述清晰，再貼一下，主函式，順便再加強下對sift 演算法的五個步驟的認識：

1、SIFT演算法第一步：影象預處理
CvMat *ScaleInitImage(CvMat * im) ;                  //金字塔初始化
2、SIFT演算法第二步：建立高斯金字塔函式
ImageOctaves* BuildGaussianOctaves(CvMat * image) ;  //建立高斯金字塔
3、SIFT演算法第三步：特徵點位置檢測，最後確定特徵點的位置
int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr);
4、SIFT演算法第四步：計算高斯影象的梯度方向和幅值，計算各個特徵點的主方向
void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr);
5、SIFT演算法第五步：抽取各個特徵點處的特徵描述字
void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr);

ok，接下來一一具體實現這幾個函式：
SIFT演算法第一步
    SIFT演算法第一步：擴大影象，預濾波剔除噪聲，得到金字塔的最底層-第一階的第一層：

  1. CvMat *ScaleInitImage(CvMat * im)   
  2. {  
  3.     double sigma,preblur_sigma;  
  4.  CvMat *imMat;  
  5.  CvMat * dst;  
  6.  CvMat *tempMat;  
  7.  //首先對影象進行平滑濾波，抑制噪聲  
  8.  imMat = cvCreateMat(im->rows, im->cols, CV_32FC1);  
  9.     BlurImage(im, imMat, INITSIGMA);  
  10.  //針對兩種情況分別進行處理：初始化放大原始影象或者在原影象基礎上進行後續操作  
  11.  //建立金字塔的最底層  
  12.  if (DOUBLE_BASE_IMAGE_SIZE)   
  13.  {  
  14.   tempMat = doubleSizeImage2(imMat);//對擴大兩倍的影象進行二次取樣，取樣率為0.5，採用線性插值  
  15. #define TEMPMAT(ROW,COL) ((float *)(tempMat->data.fl + tempMat->step/sizeof(float) * (ROW)))[(COL)]  
  16.     
  17.   dst = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);  
  18.   preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);  
  19.         BlurImage(tempMat, dst, preblur_sigma);   
  20.     
  21.   // The initial blurring for the first image of the first octave of the pyramid.  
  22.         sigma = sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );  
  23.   //  sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4);  
  24.   //printf("Init Sigma: %f/n", sigma);  
  25.   BlurImage(dst, tempMat, sigma);       //得到金字塔的最底層-放大2倍的影象  
  26.   cvReleaseMat( &dst );   
  27.   return tempMat;  
  28.  }   
  29.  else   
  30.  {  
  31.   dst = cvCreateMat(im->rows, im->cols, CV_32FC1);  
  32.   //sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA);  
  33.   preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);  
  34.   sigma = sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );  
  35.   //printf("Init Sigma: %f/n", sigma);  
  36.   BlurImage(imMat, dst, sigma);        //得到金字塔的最底層：原始影象大小  
  37.   return dst;  
  38.  }   
  39. }  
  

SIFT演算法第二步
    SIFT第二步，建立Gaussian金字塔，給定金字塔第一階第一層影象後，計算高斯金字塔其他尺度影象，
每一階的數目由變數SCALESPEROCTAVE決定，給定一個基本影象，計算它的高斯金字塔影象，返回外部向量是階梯指標，內部向量是每一個階梯內部的不同尺度影象。

  1. //SIFT演算法第二步  
  2. ImageOctaves* BuildGaussianOctaves(CvMat * image)   
  3. {  
  4.     ImageOctaves *octaves;  
  5.  CvMat *tempMat;  
  6.     CvMat *dst;  
  7.  CvMat *temp;  
  8.    
  9.  int i,j;  
  10.  double k = pow(2, 1.0/((float)SCALESPEROCTAVE));  //方差倍數  
  11.  float preblur_sigma, initial_sigma , sigma1,sigma2,sigma,absolute_sigma,sigma_f;  
  12.  //計算金字塔的階梯數目  
  13.  int dim = min(image->rows, image->cols);  
  14.     int numoctaves = (int) (log((double) dim) / log(2.0)) - 2;    //金字塔階數  
  15.  //限定金字塔的階梯數  
  16.  numoctaves = min(numoctaves, MAXOCTAVES);  
  17.  //為高斯金塔和DOG金字塔分配記憶體  
  18.  octaves=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
  19.     DOGoctaves=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
  20.    
  21.  printf("BuildGaussianOctaves(): Base image dimension is %dx%d/n", (int)(0.5*(image->cols)), (int)(0.5*(image->rows)) );  
  22.  printf("BuildGaussianOctaves(): Building %d octaves/n", numoctaves);  
  23.    
  24.     // start with initial source image  
  25.  tempMat=cvCloneMat( image );  
  26.  // preblur_sigma = 1.0;//sqrt(2 - 4*INITSIGMA*INITSIGMA);  
  27.     initial_sigma = sqrt(2);//sqrt( (4*INITSIGMA*INITSIGMA) + preblur_sigma * preblur_sigma );  
  28.  //   initial_sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4);  
  29.       
  30.  //在每一階金字塔影象中建立不同的尺度影象  
  31.  for ( i = 0; i < numoctaves; i++)   
  32.  {     
  33.   //首先建立金字塔每一階梯的最底層，其中0階梯的最底層已經建立好  
  34.   printf("Building octave %d of dimesion (%d, %d)/n", i, tempMat->cols,tempMat->rows);  
  35.         //為各個階梯分配記憶體  
  36.   octaves[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE + 3) * sizeof(ImageLevels) );  
  37.   DOGoctaves[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE + 2) * sizeof(ImageLevels) );  
  38.   //儲存各個階梯的最底層  
  39.   (octaves[i].Octave)[0].Level=tempMat;  
  40.     
  41.   octaves[i].col=tempMat->cols;  
  42.   octaves[i].row=tempMat->rows;  
  43.   DOGoctaves[i].col=tempMat->cols;  
  44.   DOGoctaves[i].row=tempMat->rows;  
  45.   if (DOUBLE_BASE_IMAGE_SIZE)  
  46.             octaves[i].subsample=pow(2,i)*0.5;  
  47.   else  
  48.             octaves[i].subsample=pow(2,i);  
  49.     
  50.   if(i==0)       
  51.   {  
  52.    (octaves[0].Octave)[0].levelsigma = initial_sigma;  
  53.    (octaves[0].Octave)[0].absolute_sigma = initial_sigma;  
  54.    printf("0 scale and blur sigma : %f /n", (octaves[0].subsample) * ((octaves[0].Octave)[0].absolute_sigma));  
  55.   }  
  56.   else  
  57.   {  
  58.    (octaves[i].Octave)[0].levelsigma = (octaves[i-1].Octave)[SCALESPEROCTAVE].levelsigma;  
  59.             (octaves[i].Octave)[0].absolute_sigma = (octaves[i-1].Octave)[SCALESPEROCTAVE].absolute_sigma;  
  60.    printf( "0 scale and blur sigma : %f /n", ((octaves[i].Octave)[0].absolute_sigma) );  
  61.   }  
  62.   sigma = initial_sigma;  
  63.         //建立本階梯其他層的影象  
  64.   for ( j =  1; j < SCALESPEROCTAVE + 3; j++)   
  65.   {  
  66.    dst = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);//用於儲存高斯層  
  67.    temp = cvCreateMat(tempMat->rows, tempMat->cols, CV_32FC1);//用於儲存DOG層  
  68.    // 2 passes of 1D on original  
  69.    //   if(i!=0)  
  70.    //   {  
  71.    //       sigma1 = pow(k, j - 1) * ((octaves[i-1].Octave)[j-1].levelsigma);  
  72.    //          sigma2 = pow(k, j) * ((octaves[i].Octave)[j-1].levelsigma);  
  73.    //       sigma = sqrt(sigma2*sigma2 - sigma1*sigma1);  
  74.    sigma_f= sqrt(k*k-1)*sigma;  
  75.    //   }  
  76.    //   else  
  77.    //   {  
  78.    //       sigma = sqrt(SIGMA * SIGMA - INITSIGMA * INITSIGMA * 4)*pow(k,j);  
  79.    //   }    
  80.             sigma = k*sigma;  
  81.    absolute_sigma = sigma * (octaves[i].subsample);  
  82.    printf("%d scale and Blur sigma: %f  /n", j, absolute_sigma);  
  83.      
  84.    (octaves[i].Octave)[j].levelsigma = sigma;  
  85.             (octaves[i].Octave)[j].absolute_sigma = absolute_sigma;  
  86.             //產生高斯層  
  87.    int length=BlurImage((octaves[i].Octave)[j-1].Level, dst, sigma_f);//相應尺度  
  88.             (octaves[i].Octave)[j].levelsigmalength = length;  
  89.    (octaves[i].Octave)[j].Level=dst;  
  90.             //產生DOG層  
  91.             cvSub( ((octaves[i].Octave)[j]).Level, ((octaves[i].Octave)[j-1]).Level, temp, 0 );  
  92.    //         cvAbsDiff( ((octaves[i].Octave)[j]).Level, ((octaves[i].Octave)[j-1]).Level, temp );  
  93.             ((DOGoctaves[i].Octave)[j-1]).Level=temp;  
  94.   }  
  95.   // halve the image size for next iteration  
  96.   tempMat  = halfSizeImage( ( (octaves[i].Octave)[SCALESPEROCTAVE].Level ) );  
  97.  }  
  98.  return octaves;  
  99. }  
  

SIFT演算法第三步
    SIFT演算法第三步，特徵點位置檢測，最後確定特徵點的位置檢測DOG金字塔中的區域性最大值，找到之後，還要經過兩個檢驗才能確認為特徵點：一是它必須有明顯的差異，二是他不應該是邊緣點，（也就是說，在極值點處的主曲率比應該小於某一個閾值）。

  1. //SIFT演算法第三步，特徵點位置檢測，  
  2. int DetectKeypoint(int numoctaves, ImageOctaves *GaussianPyr)  
  3. {  
  4.     //計算用於DOG極值點檢測的主曲率比的閾值  
  5.  double curvature_threshold;  
  6.  curvature_threshold= ((CURVATURE_THRESHOLD + 1)*(CURVATURE_THRESHOLD + 1))/CURVATURE_THRESHOLD;  
  7. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(DOGoctaves[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + DOGoctaves[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
  8.    
  9.  int   keypoint_count = 0;     
  10.  for (int i=0; i<numoctaves; i++)    
  11.  {          
  12.   for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中間的scaleperoctave個層  
  13.   {    
  14.    //在影象的有效區域內尋找具有顯著性特徵的區域性最大值  
  15.    //float sigma=(GaussianPyr[i].Octave)[j].levelsigma;  
  16.    //int dim = (int) (max(3.0f, 2.0*GAUSSKERN *sigma + 1.0f)*0.5);  
  17.    int dim = (int)(0.5*((GaussianPyr[i].Octave)[j].levelsigmalength)+0.5);  
  18.    for (int m=dim;m<((DOGoctaves[i].row)-dim);m++)   
  19.     for(int n=dim;n<((DOGoctaves[i].col)-dim);n++)  
  20.     {       
  21.      if ( fabs(ImLevels(i,j,m,n))>= CONTRAST_THRESHOLD )  
  22.      {  
  23.         
  24.       if ( ImLevels(i,j,m,n)!=0.0 )  //1、首先是非零  
  25.       {  
  26.        float inf_val=ImLevels(i,j,m,n);  
  27.        if(( (inf_val <= ImLevels(i,j-1,m-1,n-1))&&  
  28.         (inf_val <= ImLevels(i,j-1,m  ,n-1))&&  
  29.         (inf_val <= ImLevels(i,j-1,m+1,n-1))&&  
  30.         (inf_val <= ImLevels(i,j-1,m-1,n  ))&&  
  31.         (inf_val <= ImLevels(i,j-1,m  ,n  ))&&  
  32.         (inf_val <= ImLevels(i,j-1,m+1,n  ))&&  
  33.         (inf_val <= ImLevels(i,j-1,m-1,n+1))&&  
  34.         (inf_val <= ImLevels(i,j-1,m  ,n+1))&&  
  35.         (inf_val <= ImLevels(i,j-1,m+1,n+1))&&    //底層的小尺度9  
  36.           
  37.         (inf_val <= ImLevels(i,j,m-1,n-1))&&  
  38.         (inf_val <= ImLevels(i,j,m  ,n-1))&&  
  39.         (inf_val <= ImLevels(i,j,m+1,n-1))&&  
  40.         (inf_val <= ImLevels(i,j,m-1,n  ))&&  
  41.         (inf_val <= ImLevels(i,j,m+1,n  ))&&  
  42.         (inf_val <= ImLevels(i,j,m-1,n+1))&&  
  43.         (inf_val <= ImLevels(i,j,m  ,n+1))&&  
  44.         (inf_val <= ImLevels(i,j,m+1,n+1))&&     //當前層8  
  45.           
  46.         (inf_val <= ImLevels(i,j+1,m-1,n-1))&&  
  47.         (inf_val <= ImLevels(i,j+1,m  ,n-1))&&  
  48.         (inf_val <= ImLevels(i,j+1,m+1,n-1))&&  
  49.         (inf_val <= ImLevels(i,j+1,m-1,n  ))&&  
  50.         (inf_val <= ImLevels(i,j+1,m  ,n  ))&&  
  51.         (inf_val <= ImLevels(i,j+1,m+1,n  ))&&  
  52.         (inf_val <= ImLevels(i,j+1,m-1,n+1))&&  
  53.         (inf_val <= ImLevels(i,j+1,m  ,n+1))&&  
  54.         (inf_val <= ImLevels(i,j+1,m+1,n+1))     //下一層大尺度9          
  55.         ) ||   
  56.         ( (inf_val >= ImLevels(i,j-1,m-1,n-1))&&  
  57.         (inf_val >= ImLevels(i,j-1,m  ,n-1))&&  
  58.         (inf_val >= ImLevels(i,j-1,m+1,n-1))&&  
  59.         (inf_val >= ImLevels(i,j-1,m-1,n  ))&&  
  60.         (inf_val >= ImLevels(i,j-1,m  ,n  ))&&  
  61.         (inf_val >= ImLevels(i,j-1,m+1,n  ))&&  
  62.         (inf_val >= ImLevels(i,j-1,m-1,n+1))&&  
  63.         (inf_val >= ImLevels(i,j-1,m  ,n+1))&&  
  64.         (inf_val >= ImLevels(i,j-1,m+1,n+1))&&  
  65.           
  66.         (inf_val >= ImLevels(i,j,m-1,n-1))&&  
  67.         (inf_val >= ImLevels(i,j,m  ,n-1))&&  
  68.         (inf_val >= ImLevels(i,j,m+1,n-1))&&  
  69.         (inf_val >= ImLevels(i,j,m-1,n  ))&&  
  70.         (inf_val >= ImLevels(i,j,m+1,n  ))&&  
  71.         (inf_val >= ImLevels(i,j,m-1,n+1))&&  
  72.         (inf_val >= ImLevels(i,j,m  ,n+1))&&  
  73.         (inf_val >= ImLevels(i,j,m+1,n+1))&&   
  74.           
  75.         (inf_val >= ImLevels(i,j+1,m-1,n-1))&&  
  76.         (inf_val >= ImLevels(i,j+1,m  ,n-1))&&  
  77.         (inf_val >= ImLevels(i,j+1,m+1,n-1))&&  
  78.         (inf_val >= ImLevels(i,j+1,m-1,n  ))&&  
  79.         (inf_val >= ImLevels(i,j+1,m  ,n  ))&&  
  80.         (inf_val >= ImLevels(i,j+1,m+1,n  ))&&  
  81.         (inf_val >= ImLevels(i,j+1,m-1,n+1))&&  
  82.         (inf_val >= ImLevels(i,j+1,m  ,n+1))&&  
  83.         (inf_val >= ImLevels(i,j+1,m+1,n+1))   
  84.         ) )      //2、滿足26箇中極值點  
  85.        {     
  86.         //此處可儲存  
  87.         //然後必須具有明顯的顯著性，即必須大於CONTRAST_THRESHOLD=0.02  
  88.         if ( fabs(ImLevels(i,j,m,n))>= CONTRAST_THRESHOLD )  
  89.         {  
  90.          //最後顯著處的特徵點必須具有足夠的曲率比，CURVATURE_THRESHOLD=10.0，首先計算Hessian矩陣  
  91.          // Compute the entries of the Hessian matrix at the extrema location.  
  92.          /* 
  93.          1   0   -1 
  94.          0   0   0 
  95.          -1   0   1         *0.25 
  96.          */  
  97.          // Compute the trace and the determinant of the Hessian.  
  98.          //Tr_H = Dxx + Dyy;  
  99.          //Det_H = Dxx*Dyy - Dxy^2;  
  100.          float Dxx,Dyy,Dxy,Tr_H,Det_H,curvature_ratio;  
  101.          Dxx = ImLevels(i,j,m,n-1) + ImLevels(i,j,m,n+1)-2.0*ImLevels(i,j,m,n);  
  102.          Dyy = ImLevels(i,j,m-1,n) + ImLevels(i,j,m+1,n)-2.0*ImLevels(i,j,m,n);  
  103.          Dxy = ImLevels(i,j,m-1,n-1) + ImLevels(i,j,m+1,n+1) - ImLevels(i,j,m+1,n-1) - ImLevels(i,j,m-1,n+1);  
  104.          Tr_H = Dxx + Dyy;  
  105.          Det_H = Dxx*Dyy - Dxy*Dxy;  
  106.          // Compute the ratio of the principal curvatures.  
  107.          curvature_ratio = (1.0*Tr_H*Tr_H)/Det_H;  
  108.          if ( (Det_H>=0.0) && (curvature_ratio <= curvature_threshold) )  //最後得到最具有顯著性特徵的特徵點  
  109.          {  
  110.           //將其儲存起來，以計算後面的特徵描述字  
  111.           keypoint_count++;  
  112.           Keypoint k;  
  113.           /* Allocate memory for the keypoint. */  
  114.           k = (Keypoint) malloc(sizeof(struct KeypointSt));  
  115.           k->next = keypoints;  
  116.           keypoints = k;  
  117.           k->row = m*(GaussianPyr[i].subsample);  
  118.           k->col =n*(GaussianPyr[i].subsample);  
  119.           k->sy = m;    //行  
  120.           k->sx = n;    //列  
  121.           k->octave=i;  
  122.           k->level=j;  
  123.           k->scale = (GaussianPyr[i].Octave)[j].absolute_sigma;        
  124.          }//if >curvature_thresh  
  125.         }//if >contrast  
  126.        }//if inf value  
  127.      }//if non zero  
  128.      }//if >contrast  
  129.     }  //for concrete image level col  
  130.   }//for levels  
  131.  }//for octaves  
  132.  return keypoint_count;  
  133. }  
  134.   
  135. //在影象中，顯示SIFT特徵點的位置  
  136. void DisplayKeypointLocation(IplImage* image, ImageOctaves *GaussianPyr)  
  137. {  
  138.    
  139.     Keypoint p = keypoints; // p指向第一個結點  
  140.  while(p) // 沒到表尾  
  141.  {     
  142.   cvLine( image, cvPoint((int)((p->col)-3),(int)(p->row)),   
  143.    cvPoint((int)((p->col)+3),(int)(p->row)), CV_RGB(255,255,0),  
  144.    1, 8, 0 );  
  145.         cvLine( image, cvPoint((int)(p->col),(int)((p->row)-3)),   
  146.    cvPoint((int)(p->col),(int)((p->row)+3)), CV_RGB(255,255,0),  
  147.    1, 8, 0 );  
  148.   //  cvCircle(image,cvPoint((uchar)(p->col),(uchar)(p->row)),  
  149.   //   (int)((GaussianPyr[p->octave].Octave)[p->level].absolute_sigma),  
  150.   //   CV_RGB(255,0,0),1,8,0);  
  151.   p=p->next;  
  152.  }   
  153. }  
  154.   
  155. // Compute the gradient direction and magnitude of the gaussian pyramid images  
  156. void ComputeGrad_DirecandMag(int numoctaves, ImageOctaves *GaussianPyr)  
  157. {  
  158.  // ImageOctaves *mag_thresh ;  
  159.     mag_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
  160.  grad_pyr=(ImageOctaves*) malloc( numoctaves * sizeof(ImageOctaves) );  
  161.  // float sigma=( (GaussianPyr[0].Octave)[SCALESPEROCTAVE+2].absolute_sigma ) / GaussianPyr[0].subsample;  
  162.  // int dim = (int) (max(3.0f, 2 * GAUSSKERN *sigma + 1.0f)*0.5+0.5);  
  163. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + GaussianPyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
  164.  for (int i=0; i<numoctaves; i++)    
  165.  {          
  166.   mag_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );  
  167.         grad_pyr[i].Octave= (ImageLevels*) malloc( (SCALESPEROCTAVE) * sizeof(ImageLevels) );  
  168.   for(int j=1;j<SCALESPEROCTAVE+1;j++)//取中間的scaleperoctave個層  
  169.   {    
  170.             CvMat *Mag = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
  171.    CvMat *Ori = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
  172.    CvMat *tempMat1 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
  173.    CvMat *tempMat2 = cvCreateMat(GaussianPyr[i].row, GaussianPyr[i].col, CV_32FC1);  
  174.    cvZero(Mag);  
  175.    cvZero(Ori);  
  176.    cvZero(tempMat1);  
  177.    cvZero(tempMat2);   
  178. #define MAG(ROW,COL) ((float *)(Mag->data.fl + Mag->step/sizeof(float) *(ROW)))[(COL)]     
  179. #define ORI(ROW,COL) ((float *)(Ori->data.fl + Ori->step/sizeof(float) *(ROW)))[(COL)]    
  180. #define TEMPMAT1(ROW,COL) ((float *)(tempMat1->data.fl + tempMat1->step/sizeof(float) *(ROW)))[(COL)]  
  181. #define TEMPMAT2(ROW,COL) ((float *)(tempMat2->data.fl + tempMat2->step/sizeof(float) *(ROW)))[(COL)]  
  182.    for (int m=1;m<(GaussianPyr[i].row-1);m++)   
  183.     for(int n=1;n<(GaussianPyr[i].col-1);n++)  
  184.     {  
  185.      //計算幅值  
  186.      TEMPMAT1(m,n) = 0.5*( ImLevels(i,j,m,n+1)-ImLevels(i,j,m,n-1) );  //dx  
  187.                     TEMPMAT2(m,n) = 0.5*( ImLevels(i,j,m+1,n)-ImLevels(i,j,m-1,n) );  //dy  
  188.                     MAG(m,n) = sqrt(TEMPMAT1(m,n)*TEMPMAT1(m,n)+TEMPMAT2(m,n)*TEMPMAT2(m,n));  //mag  
  189.      //計算方向  
  190.      ORI(m,n) =atan( TEMPMAT2(m,n)/TEMPMAT1(m,n) );  
  191.                     if (ORI(m,n)==CV_PI)  
  192.       ORI(m,n)=-CV_PI;  
  193.     }  
  194.     ((mag_pyr[i].Octave)[j-1]).Level=Mag;  
  195.     ((grad_pyr[i].Octave)[j-1]).Level=Ori;  
  196.     cvReleaseMat(&tempMat1);  
  197.     cvReleaseMat(&tempMat2);  
  198.   }//for levels  
  199.  }//for octaves  
  200. }  

SIFT演算法第四步
  1. //SIFT演算法第四步：計算各個特徵點的主方向，確定主方向  
  2. void AssignTheMainOrientation(int numoctaves, ImageOctaves *GaussianPyr,ImageOctaves *mag_pyr,ImageOctaves *grad_pyr)  
  3. {  
  4.     // Set up the histogram bin centers for a 36 bin histogram.  
  5.     int num_bins = 36;  
  6.     float hist_step = 2.0*PI/num_bins;  
  7.     float hist_orient[36];  
  8.     for (int i=0;i<36;i++)  
  9.   hist_orient[i]=-PI+i*hist_step;  
  10.     float sigma1=( ((GaussianPyr[0].Octave)[SCALESPEROCTAVE].absolute_sigma) ) / (GaussianPyr[0].subsample);//SCALESPEROCTAVE+2  
  11.  int zero_pad = (int) (max(3.0f, 2 * GAUSSKERN *sigma1 + 1.0f)*0.5+0.5);  
  12.     //Assign orientations to the keypoints.  
  13. #define ImLevels(OCTAVES,LEVELS,ROW,COL) ((float *)((GaussianPyr[(OCTAVES)].Octave[(LEVELS)].Level)->data.fl + (GaussianPyr[(OCTAVES)].Octave[(LEVELS)].Level)->step/sizeof(float) *(ROW)))[(COL)]  
  14.    
  15.     int keypoint_count = 0;  
  16.  Keypoint p = keypoints; // p指向第一個結點  
  17.    
  18.  while(p) // 沒到表尾  
  19.  {  
  20.   int i=p->octave;  
  21.   int j=p->level;  
  22.   int m=p->sy;   //行  
  23.   int n=p->sx;   //列  
  24.   if ((m>=zero_pad)&&(m<GaussianPyr[i].row-zero_pad)&&  
  25.    (n>=zero_pad)&&(n<GaussianPyr[i].col-zero_pad) )  
  26.   {  
  27.    float sigma=( ((GaussianPyr[i].Octave)[j].absolute_sigma) ) / (GaussianPyr[i].subsample);  
  28.    //產生二維高斯模板  
  29.             CvMat* mat = GaussianKernel2D( sigma );           
  30.    int dim=(int)(0.5 * (mat->rows));  
  31.    //分配用於儲存Patch幅值和方向的空間  
  32. #define MAT(ROW,COL) ((float *)(mat->data.fl + mat->step/sizeof(float) *(ROW)))[(COL)]  
  33.      
  34.    //宣告方向直方圖變數  
  35.    double* orienthist = (double *) malloc(36 * sizeof(double));  
  36.    for ( int sw = 0 ; sw < 36 ; ++sw)   
  37.    {  
  38.     orienthist[sw]=0.0;    
  39.    }  
  40.    //在特徵點的周圍統計梯度方向  
  41.             for (int x=m-dim,mm=0;x<=(m+dim);x++,mm++)   
  42.     for(int y=n-dim,nn=0;y<=(n+dim);y++,nn++)  
  43.     {       
  44.      //計算特徵點處的幅值  
  45.      double dx = 0.5*(ImLevels(i,j,x,y+1)-ImLevels(i,j,x,y-1));  //dx  
  46.      double dy = 0.5*(ImLevels(i,j,x+1,y)-ImLevels(i,j,x-1,y));  //dy  
  47.      double mag = sqrt(dx*dx+dy*dy);  //mag  
  48.      //計算方向  
  49.      double Ori =atan( 1.0*dy/dx );  
  50.      int binIdx = FindClosestRotationBin(36, Ori);                   //得到離現有方向最近的直方塊  
  51.      orienthist[binIdx] = orienthist[binIdx] + 1.0* mag * MAT(mm,nn);//利用高斯加權累加進直方圖相應的塊  
  52.     }  
  53.     // Find peaks in the orientation histogram using nonmax suppression.  
  54.     AverageWeakBins (orienthist, 36);  
  55.     // find the maximum peak in gradient orientation  
  56.     double maxGrad = 0.0;  
  57.     int maxBin = 0;  
  58.     for (int b = 0 ; b < 36 ; ++b)   
  59.     {  
  60.      if (orienthist[b] > maxGrad)   
  61.      {  
  62.       maxGrad = orienthist[b];  
  63.       maxBin = b;  
  64.      }  
  65.     }  
  66.     // First determine the real interpolated peak high at the maximum bin  
  67.     // position, which is guaranteed to be an absolute peak.  
  68.     double maxPeakValue=0.0;  
  69.     double maxDegreeCorrection=0.0;  
  70.     if ( (InterpolateOrientation ( orienthist[maxBin == 0 ? (36 - 1) : (maxBin - 1)],  
  71.                     orienthist[maxBin], orienthist[(maxBin + 1) % 36],  
  72.                     &maxDegreeCorrection, &maxPeakValue)) == false)  
  73.      printf("BUG: Parabola fitting broken");  
  74.       
  75.     // Now that we know the maximum peak value, we can find other keypoint  
  76.     // orientations, which have to fulfill two criterias:  
  77.     //  
  78.     //  1. They must be a local peak themselves. Else we might add a very  
  79.     //     similar keypoint orientation twice (imagine for example the  
  80.     //     values: 0.4 1.0 0.8, if 1.0 is maximum peak, 0.8 is still added  
  81.     //     with the default threshhold, but the maximum peak orientation  
  82.     //     was already added).  
  83.     //  2. They must have at least peakRelThresh times the maximum peak  
  84.     //     value.  
  85.     bool binIsKeypoint[36];  
  86.     for ( b = 0 ; b < 36 ; ++b)   
  87.     {  
  88.      binIsKeypoint[b] = false;  
  89.      // The maximum peak of course is  
  90.      if (b == maxBin)   
  91.      {  
  92.       binIsKeypoint[b] = true;  
  93.       continue;  
  94.      }  
  95.      // Local peaks are, too, in case they fulfill the threshhold  
  96.      if (orienthist[b] < (peakRelThresh * maxPeakValue))  
  97.       continue;  
  98.      int leftI = (b == 0) ? (36 - 1) : (b - 1);  
  99.      int rightI = (b + 1) % 36;  
  100.      if (orienthist[b] <= orienthist[leftI] || orienthist[b] <= orienthist[rightI])  
  101.       continue; // no local peak  
  102.      binIsKeypoint[b] = true;  
  103.     }  
  104.     // find other possible locations  
  105.     double oneBinRad = (2.0 * PI) / 36;  
  106.     for ( b = 0 ; b < 36 ; ++b)   
  107.     {  
  108.      if (binIsKeypoint[b] == false)  
  109.       continue;  
  110.      int bLeft = (b == 0) ? (36 - 1) : (b - 1);  
  111.      int bRight = (b + 1) % 36;  
  112.      // Get an interpolated peak direction and value guess.  
  113.      double peakValue;  
  114.      double degreeCorrection;  
  115.        
  116.      double maxPeakValue, maxDegreeCorrection;                
  117.      if (InterpolateOrientation ( orienthist[maxBin == 0 ? (36 - 1) : (maxBin - 1)],  
  118.       orienthist[maxBin], orienthist[(maxBin + 1) % 36],  
  119.       °reeCorrection, &peakValue) == false)  
  120.      {  
  121.       printf("BUG: Parabola fitting broken");  
  122.      }  
  123.        
  124.      double degree = (b + degreeCorrection) * oneBinRad - PI;  
  125.      if (degree < -PI)  
  126.       degree += 2.0 * PI;  
  127.      else if (degree > PI)  
  128.       degree -= 2.0 * PI;  
  129.      //儲存方向，可以直接利用檢測到的連結串列進行該步主方向的指定;  
  130.      //分配記憶體重新儲存特徵點  
  131.      Keypoint k;  
  132.      /* Allocate memory for the keypoint Descriptor. */  
  133.      k = (Keypoint) malloc(sizeof(struct KeypointSt));  
  134.      k->next = keyDescriptors;  
  135.      keyDescriptors = k;  
  136.      k->descrip = (float*)malloc(LEN * sizeof(float));  
  137.      k->row = p->row;  
  138.      k->col = p->col;  
  139.      k->sy = p->sy;    //行  
  140.      k->sx = p->sx;    //列  
  141.      k->octave = p->octave;  
  142.      k->level = p->level;  
  143.      k->scale = p->scale;        
  144.      k->ori = degree;  
  145.      k->mag = peakValue;    
  146.     }//for  
  147.     free(orienthist);  
  148.   }  
  149.   p=p->next;  
  150.  }   
  151. }  
  152.   
  153. //尋找與方向直方圖最近的柱，確定其index   
  154. int FindClosestRotationBin (int binCount, float angle)  
  155. {  
  156.  angle += CV_PI;  
  157.  angle /= 2.0 * CV_PI;  
  158.  // calculate the aligned bin  
  159.  angle *= binCount;  
  160.  int idx = (int) angle;  
  161.  if (idx == binCount)  
  162.   idx = 0;  
  163.  return (idx);  
  164. }  
  165.   
  166. // Average the content of the direction bins.  
  167. void AverageWeakBins (double* hist, int binCount)  
  168. {  
  169.  // TODO: make some tests what number of passes is the best. (its clear  
  170.  // one is not enough, as we may have something like  
  171.  // ( 0.4, 0.4, 0.3, 0.4, 0.4 ))  
  172.  for (int sn = 0 ; sn < 2 ; ++sn)   
  173.  {  
  174.   double firstE = hist[0];  
  175.   double last = hist[binCount-1];  
  176.   for (int sw = 0 ; sw < binCount ; ++sw)   
  177.   {  
  178.    double cur = hist[sw];  
  179.    double next = (sw == (binCount - 1)) ? firstE : hist[(sw + 1) % binCount];  
  180.    hist[sw] = (last + cur + next) / 3.0;  
  181.    last = cur;  
  182.   }  
  183.  }  
  184. }  
  185.   
  186. // Fit a parabol to the three points (-1.0 ; left), (0.0 ; middle) and  
  187. // (1.0 ; right).  
  188. // Formulas:  
  189. // f(x) = a (x - c)^2 + b  
  190. // c is the peak offset (where f'(x) is zero), b is the peak value.  
  191. // In case there is an error false is returned, otherwise a correction  
  192. // value between [-1 ; 1] is returned in 'degreeCorrection', where -1  
  193. // means the peak is located completely at the left vector, and -0.5 just  
  194. // in the middle between left and middle and > 0 to the right side. In  
  195. // 'peakValue' the maximum estimated peak value is stored.  
  196. bool InterpolateOrientation (double left, double middle,double right, double *degreeCorrection, double *peakValue)  
  197. {  
  198.  double a = ((left + right) - 2.0 * middle) / 2.0;   //拋物線捏合係數a  
  199.  // degreeCorrection = peakValue = Double.NaN;  
  200.    
  201.  // Not a parabol  
  202.  if (a == 0.0)  
  203.   return false;  
  204.  double c = (((left - middle) / a) - 1.0) / 2.0;  
  205.  double b = middle - c * c * a;  
  206.  if (c < -0.5 || c > 0.5)  
  207.   return false;  
  208.  *degreeCorrection = c;  
  209.  *peakValue = b;  
  210.  return true;  
  211. }  
  212.   
  213. //顯示特徵點處的主方向  
  214. void DisplayOrientation (IplImage* image, ImageOctaves *GaussianPyr)  
  215. {  
  216.  Keypoint p = keyDescriptors; // p指向第一個結點  
  217.  while(p) // 沒到表尾  
  218.  {  
  219.   float scale=(GaussianPyr[p->octave].Octave)[p->level].absolute_sigma;  
  220.   float autoscale = 3.0;   
  221.   float uu=autoscale*scale*cos(p->ori);  
  222.   float vv=autoscale*scale*sin(p->ori);  
  223.   float x=(p->col)+uu;  
  224.   float y=(p->row)+vv;  
  225.   cvLine( image, cvPoint((int)(p->col),(int)(p->row)),   
  226.    cvPoint((int)x,(int)y), CV_RGB(255,255,0),  
  227.    1, 8, 0 );  
  228.   // Arrow head parameters  
  229.         float alpha = 0.33; // Size of arrow head relative to the length of the vector  
  230.         float beta = 0.33;  // Width of the base of the arrow head relative to the length  
  231.     
  232.   float xx0= (p->col)+uu-alpha*(uu+beta*vv);  
  233.         float yy0= (p->row)+vv-alpha*(vv-beta*uu);  
  234.         float xx1= (p->col)+uu-alpha*(uu-beta*vv);  
  235.         float yy1= (p->row)+vv-alpha*(vv+beta*uu);  
  236.         cvLine( image, cvPoint((int)xx0,(int)yy0),   
  237.    cvPoint((int)x,(int)y), CV_RGB(255,255,0),  
  238.    1, 8, 0 );  
  239.   cvLine( image, cvPoint((int)xx1,(int)yy1),   
  240.    cvPoint((int)x,(int)y), CV_RGB(255,255,0),  
  241.    1, 8, 0 );  
  242.   p=p->next;  
  243.  }   
  244. }  

SIFT演算法第五步
    SIFT演算法第五步：抽取各個特徵點處的特徵描述字，確定特徵點的描述字。描述字是Patch網格內梯度方向的描述，旋轉網格到主方向，插值得到網格處梯度值。
一個特徵點可以用2*2*8=32維的向量，也可以用4*4*8=128維的向量更精確的進行描述。

  1. void ExtractFeatureDescriptors(int numoctaves, ImageOctaves *GaussianPyr)  
  2. {  
  3.  // The orientation histograms have 8 bins  
  4.  float orient_bin_spacing = PI/4;  
  5.     float orient_angles[8]={-PI,-PI+orient_bin_spacing,-PI*0.5, -orient_bin_spacing,  
  6.   0.0, orient_bin_spacing, PI*0.5,  PI+orient_bin_spacing};  
  7.     //產生描述字中心各點座標  
  8.     float *feat_grid=(float *) malloc( 2*16 * sizeof(float));  
  9.  for (int i=0;i<GridSpacing;i++)  
  10.  {  
  11.         for (int j=0;j<2*GridSpacing;++j,++j)  
  12.   {  
  13.    feat_grid[i*2*GridSpacing+j]=-6.0+i*GridSpacing;  
  14.    feat_grid[i*2*GridSpacing+j+1]=-6.0+0.5*j*GridSpacing;  
  15.   }  
  16.  }  
  17.     //產生網格  
  18.  float *feat_samples=(float *) malloc( 2*256 * sizeof(float));  
  19.     for ( i=0;i<4*GridSpacing;i++)  
  20.  {  
  21.         for (int j=0;j<8*GridSpacing;j+=2)  
  22.   {  
  23.    feat_samples[i*8*GridSpacing+j]=-(2*GridSpacing-0.5)+i;  
  24.    feat_samples[i*8*GridSpacing+j+1]=-(2*GridSpacing-0.5)+0.5*j;  
  25.   }  
  26.  }  
  27.  float feat_window = 2*GridSpacing;  
  28.  Keypoint p = keyDescriptors; // p指向第一個結點  
  29.  while(p) // 沒到表尾  
  30.  {  
  31.   float scale=(GaussianPyr[p->octave].Octave)[p->level].absolute_sigma;  
  32.           
  33.   float sine = sin(p->ori);  
  34.   float cosine = cos(p->ori);    
  35.         //計算中心點座標旋轉之後的位置  
  36.   float *featcenter=(float *) malloc( 2*16 * sizeof(float));  
  37.   for (int i=0;i<GridSpacing;i++)  
  38.   {  
  39.             for (int j=0;j<2*GridSpacing;j+=2)  
  40.    {  
  41.     float x=feat_grid[i*2*GridSpacing+j];  
  42.     float y=feat_grid[i*2*GridSpacing+j+1];  
  43.     featcenter[i*2*GridSpacing+j]=((cosine * x + sine * y) + p->sx);  
  44.     featcenter[i*2*GridSpacing+j+1]=((-sine * x + cosine * y) + p->sy);  
  45.    }  
  46.   }  
  47.   // calculate sample window coordinates (rotated along keypoint)  
  48.   float *feat=(float *) malloc( 2*256 * sizeof(float));  
  49.   for ( i=0;i<64*GridSpacing;i++,i++)  
  50.   {  
  51.    float x=feat_samples[i];  
  52.    float y=feat_samples[i+1];  
  53.    feat[i]=((cosine * x + sine * y) + p->sx);  
  54.    feat[i+1]=((-sine * x + cosine * y) + p->sy);  
  55.   }  
  56.   //Initialize the feature descriptor.  
  57.         float *feat_desc = (float *) malloc( 128 * sizeof(float));  
  58.         for (i=0;i<128;i++)  
  59.   {  
  60.    feat_desc[i]=0.0;  
  61.    // printf("%f  ",feat_desc[i]);    
  62.   }  
  63.         //printf("/n");  
  64.   for ( i=0;i<512;++i,++i)  
  65.   {  
  66.             float x_sample = feat[i];  
  67.             float y_sample = feat[i+1];  
  68.             // Interpolate the gradient at the sample position  
  69.       /* 
  70.       0   1   0 
  71.       1   *   1 
  72.       0   1   0   具體插值策略如圖示 
  73.    */  
  74.    float sample12=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample-1);  
  75.    float sample21=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample-1, y_sample);   
  76.    float sample22=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample);   
  77.    float sample23=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample+1, y_sample);   
  78.    float sample32=getPixelBI(((GaussianPyr[p->octave].Octave)[p->level]).Level, x_sample, y_sample+1);   
  79.    //float diff_x = 0.5*(sample23 - sample21);  
  80.             //float diff_y = 0.5*(sample32 - sample12);  
  81.    float diff_x = sample23 - sample21;  
  82.             float diff_y = sample32 - sample12;  
  83.             float mag_sample = sqrt( diff_x*diff_x + diff_y*diff_y );  
  84.             float grad_sample = atan( diff_y / diff_x );  
  85.             if(grad_sample == CV_PI)  
  86.     grad_sample = -CV_PI;  
  87.             // Compute the weighting for the x and y dimensions.  
  88.             float *x_wght=(float *) malloc( GridSpacing * GridSpacing * sizeof(float));  
  89.             float *y_wght=(float *) malloc( GridSpacing * GridSpacing * sizeof(float));  
  90.    float *pos_wght=(float *) malloc( 8*GridSpacing * GridSpacing * sizeof(float));;  
  91.    for (int m=0;m<32;++m,++m)  
  92.    {  
  93.                 float x=featcenter[m];  
  94.     float y=featcenter[m+1];  
  95.     x_wght[m/2] = max(1 - (fabs(x - x_sample)*1.0/GridSpacing), 0);  
  96.                 y_wght[m/2] = max(1 - (fabs(y - y_sample)*1.0/GridSpacing), 0);   
  97.       
  98.    }  
  99.    for ( m=0;m<16;++m)  
  100.                 for (int n=0;n<8;++n)  
  101.      pos_wght[m*8+n]=x_wght[m]*y_wght[m];  
  102.     free(x_wght);  
  103.     free(y_wght);  
  104.     //計算方向的加權，首先旋轉梯度場到主方向，然後計算差異   
  105.     float diff[8],orient_wght[128];  
  106.     for ( m=0;m<8;++m)  
  107.     {   
  108.      float angle = grad_sample-(p->ori)-orient_angles[m]+CV_PI;  
  109.      float temp = angle / (2.0 * CV_PI);  
  110.      angle -= (int)(temp) * (2.0 * CV_PI);  
  111.      diff[m]= angle - CV_PI;  
  112.     }  
  113.     // Compute the gaussian weighting.  
  114.     float x=p->sx;  
  115.     float y=p->sy;  
  116.     float g = exp(-((x_sample-x)*(x_sample-x)+(y_sample-y)*(y_sample-y))/(2*feat_window*feat_window))/(2*CV_PI*feat_window*feat_window);  
  117.       
  118.     for ( m=0;m<128;++m)  
  119.     {  
  120.      orient_wght[m] = max((1.0 - 1.0*fabs(diff[m%8])/orient_bin_spacing),0);  
  121.      feat_desc[m] = feat_desc[m] + orient_wght[m]*pos_wght[m]*g*mag_sample;  
  122.     }  
  123.     free(pos_wght);     
  124.   }  
  125.   free(feat);  
  126.   free(featcenter);  
  127.   float norm=GetVecNorm( feat_desc, 128);  
  128.   for (int m=0;m<128;m++)  
  129.   {  
  130.    feat_desc[m]/=norm;  
  131.    if (feat_desc[m]>0.2)  
  132.     feat_desc[m]=0.2;  
  133.   }  
  134.         norm=GetVecNorm( feat_desc, 128);  
  135.         for ( m=0;m<128;m++)  
  136.   {  
  137.    feat_desc[m]/=norm;  
  138.    printf("%f  ",feat_desc[m]);    
  139.   }  
  140.   printf("/n");  
  141.   p->descrip = feat_desc;  
  142.   p=p->next;  
  143.  }  
  144.  free(feat_grid);  
  145.     free(feat_samples);  
  146. }  
  147.   
  148. //為了顯示圖象金字塔，而作的影象水平拼接  
  149. CvMat* MosaicHorizen( CvMat* im1, CvMat* im2 )  
  150. {  
  151.  int row,col;  
  152.  CvMat *mosaic = cvCreateMat( max(im1->rows,im2->rows),(im1->cols+im2->cols),CV_32FC1);  
  153. #define Mosaic(ROW,COL) ((float*)(mosaic->data.fl + mosaic->step/sizeof(float)*(ROW)))[(COL)]  
  154. #define Im11Mat(ROW,COL) ((float *)(im1->data.fl + im1->step/sizeof(float) *(ROW)))[(COL)]  
  155. #define Im22Mat(ROW,COL) ((float *)(im2->data.fl + im2->step/sizeof(float) *(ROW)))[(COL)]  
  156.  cvZero(mosaic);  
  157.  /* Copy images into mosaic1. */  
  158.  for ( row = 0; row < im1->rows; row++)  
  159.   for ( col = 0; col < im1->cols; col++)  
  160.    Mosaic(row,col)=Im11Mat(row,col) ;  
  161.   for (  row = 0; row < im2->rows; row++)  
  162.    for (  col = 0; col < im2->cols; col++)  
  163.     Mosaic(row, (col+im1->cols) )= Im22Mat(row,col) ;  
  164.    return mosaic;  
  165. }  
  166.   
  167. //為了顯示圖象金字塔，而作的影象垂直拼接  
  168. CvMat* MosaicVertical( CvMat* im1, CvMat* im2 )  
  169. {  
  170.  int row,col;  
  171.  CvMat *mosaic = cvCreateMat(im1->rows+im2->rows,max(im1->cols,im2->cols), CV_32FC1);  
  172. #define Mosaic(ROW,COL) ((float*)(mosaic->data.fl + mosaic->step/sizeof(float)*(ROW)))[(COL)]  
  173. #define Im11Mat(ROW,COL) ((float *)(im1->data.fl + im1->step/sizeof(float) *(ROW)))[(COL)]  
  174. #define Im22Mat(ROW,COL) ((float *)(im2->data.fl + im2->step/sizeof(float) *(ROW)))[(COL)]  
  175.  cvZero(mosaic);  
  176.    
  177.  /* Copy images into mosaic1. */  
  178.  for ( row = 0; row < im1->rows; row++)  
  179.   for ( col = 0; col < im1->cols; col++)  
  180.    Mosaic(row,col)= Im11Mat(row,col) ;  
  181.   for ( row = 0; row < im2->rows; row++)  
  182.    for ( col = 0; col < im2->cols; col++)  
  183.     Mosaic((row+im1->rows),col)=Im22Mat(row,col) ;  
  184.      
  185.    return mosaic;  
  186. }  

 

ok，為了版述清晰，再貼一下上文所述的主函式（注，上文已貼出，此是為了版述清晰，重複造輪）：

  1. int main( void )  
  2. {  
  3.  //宣告當前幀IplImage指標  
  4.  IplImage* src = NULL;   
  5.  IplImage* image1 = NULL;   
  6.  IplImage* grey_im1 = NULL;   
  7.  IplImage* DoubleSizeImage = NULL;  
  8.    
  9.  IplImage* mosaic1 = NULL;   
  10.  IplImage* mosaic2 = NULL;   
  11.    
  12.  CvMat* mosaicHorizen1 = NULL;  
  13.  CvMat* mosaicHorizen2 = NULL;  
  14.  CvMat* mosaicVertical1 = NULL;  
  15.    
  16.  CvMat* image1Mat = NULL;  
  17.  CvMat* tempMat=NULL;  
  18.    
  19.  ImageOctaves *Gaussianpyr;  
  20.  int rows,cols;  
  21.   
  22. #define Im1Mat(ROW,COL) ((float *)(image1Mat->data.fl + image1Mat->step/sizeof(float) *(ROW)))[(COL)]  
  23.    
  24.  //灰度圖象畫素的資料結構  
  25. #define Im1B(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3]  
  26. #define Im1G(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+1]  
  27. #define Im1R(ROW,COL) ((uchar*)(image1->imageData + image1->widthStep*(ROW)))[(COL)*3+2]  
  28.    
  29.  storage = cvCreateMemStorage(0);   
  30.    
  31.  //讀取圖片  
  32.  if( (src = cvLoadImage( "street1.jpg", 1)) == 0 )  // test1.jpg einstein.pgm back1.bmp  
  33.   return -1;  
  34.   
  35.  //為影象分配記憶體   
  36.  image1 = cvCreateImage(cvSize(src->width, src->height),  IPL_DEPTH_8U,3);  
  37.  grey_im1 = cvCreateImage(cvSize(src->width, src->height),  IPL_DEPTH_8U,1);  
  38.  DoubleSizeImage = cvCreateImage(cvSize(2*(src->width), 2*(src->height)),  IPL_DEPTH_8U,3);  
  39.   
  40.  //為影象陣列分配記憶體，假設兩幅影象的大小相同，tempMat跟隨image1的大小  
  41.  image1Mat = cvCreateMat(src->height, src->width, CV_32FC1);  
  42.  //轉化成單通道影象再處理  
  43.  cvCvtColor(src, grey_im1, CV_BGR2GRAY);  
  44.  //轉換進入Mat資料結構,影象操作使用的是浮點型操作  
  45.  cvConvert(grey_im1, image1Mat);  
  46.    
  47.  double t = (double)cvGetTickCount();  
  48.  //影象歸一化  
  49.  cvConvertScale( image1Mat, image1Mat, 1.0/255, 0 );  
  50.    
  51.  int dim = min(image1Mat->rows, image1Mat->cols);  
  52.  numoctaves = (int) (log((double) dim) / log(2.0)) - 2;    //金字塔階數  
  53.  numoctaves = min(numoctaves, MAXOCTAVES);  
  54.   
  55.  //SIFT演算法第一步，預濾波除噪聲，建立金字塔底層  
  56.  tempMat = ScaleInitImage(image1Mat) ;  
  57.  //SIFT演算法第二步，建立Guassian金字塔和DOG金字塔  
  58.  Gaussianpyr = BuildGaussianOctaves(tempMat) ;  
  59.    
  60.  t = (double)cvGetTickCount() - t;  
  61.  printf( "the time of build Gaussian pyramid and DOG pyramid is %.1f/n", t/(cvGetTickFrequency()*1000.) );  
  62.    
  63. #define ImLevels(OCTAVE,LEVEL,ROW,COL) ((float *)(Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->data.fl + Gaussianpyr[(OCTAVE)].Octave[(LEVEL)].Level->step/sizeof(float) *(ROW)))[(COL)]  
  64.  //顯示高斯金字塔  
  65.  for (int i=0; i<numoctaves;i++)  
  66.  {  
  67.   if (i==0)  
  68.   {  
  69.    mosaicHorizen1=MosaicHorizen( (Gaussianpyr[0].Octave)[0].Level, (Gaussianpyr[0].Octave)[1].Level );  
  70.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
  71.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[0].Octave)[j].Level );  
  72.    for ( j=0;j<NUMSIZE;j++)  
  73.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
  74.   }  
  75.   else if (i==1)  
  76.   {  
  77.    mosaicHorizen2=MosaicHorizen( (Gaussianpyr[1].Octave)[0].Level, (Gaussianpyr[1].Octave)[1].Level );  
  78.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
  79.     mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (Gaussianpyr[1].Octave)[j].Level );  
  80.    for ( j=0;j<NUMSIZE;j++)  
  81.     mosaicHorizen2=halfSizeImage(mosaicHorizen2);  
  82.    mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );  
  83.   }  
  84.   else  
  85.   {  
  86.    mosaicHorizen1=MosaicHorizen( (Gaussianpyr[i].Octave)[0].Level, (Gaussianpyr[i].Octave)[1].Level );  
  87.    for (int j=2;j<SCALESPEROCTAVE+3;j++)  
  88.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (Gaussianpyr[i].Octave)[j].Level );  
  89.    for ( j=0;j<NUMSIZE;j++)  
  90.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
  91.    mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );  
  92.   }  
  93.  }  
  94.  mosaic1 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height),  IPL_DEPTH_8U,1);  
  95.  cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0, 0 );  
  96.  cvConvertScaleAbs( mosaicVertical1, mosaic1, 1, 0 );  
  97.    
  98.  //  cvSaveImage("GaussianPyramid of me.jpg",mosaic1);  
  99.  cvNamedWindow("mosaic1",1);  
  100.  cvShowImage("mosaic1", mosaic1);  
  101.  cvWaitKey(0);  
  102.  cvDestroyWindow("mosaic1");  
  103.  //顯示DOG金字塔  
  104.  for ( i=0; i<numoctaves;i++)  
  105.  {  
  106.   if (i==0)  
  107.   {  
  108.    mosaicHorizen1=MosaicHorizen( (DOGoctaves[0].Octave)[0].Level, (DOGoctaves[0].Octave)[1].Level );  
  109.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
  110.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[0].Octave)[j].Level );  
  111.    for ( j=0;j<NUMSIZE;j++)  
  112.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
  113.   }  
  114.   else if (i==1)  
  115.   {  
  116.    mosaicHorizen2=MosaicHorizen( (DOGoctaves[1].Octave)[0].Level, (DOGoctaves[1].Octave)[1].Level );  
  117.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
  118.     mosaicHorizen2=MosaicHorizen( mosaicHorizen2, (DOGoctaves[1].Octave)[j].Level );  
  119.    for ( j=0;j<NUMSIZE;j++)  
  120.     mosaicHorizen2=halfSizeImage(mosaicHorizen2);  
  121.    mosaicVertical1=MosaicVertical( mosaicHorizen1, mosaicHorizen2 );  
  122.   }  
  123.   else  
  124.   {  
  125.    mosaicHorizen1=MosaicHorizen( (DOGoctaves[i].Octave)[0].Level, (DOGoctaves[i].Octave)[1].Level );  
  126.    for (int j=2;j<SCALESPEROCTAVE+2;j++)  
  127.     mosaicHorizen1=MosaicHorizen( mosaicHorizen1, (DOGoctaves[i].Octave)[j].Level );  
  128.    for ( j=0;j<NUMSIZE;j++)  
  129.     mosaicHorizen1=halfSizeImage(mosaicHorizen1);  
  130.    mosaicVertical1=MosaicVertical( mosaicVertical1, mosaicHorizen1 );  
  131.   }  
  132.  }  
  133.  //考慮到DOG金字塔各層影象都會有正負，所以，必須尋找最負的，以將所有影象抬高一個臺階去顯示  
  134.  double min_val=0;  
  135.  double max_val=0;  
  136.  cvMinMaxLoc( mosaicVertical1, &min_val, &max_val,NULL, NULL, NULL );  
  137.  if ( min_val<0.0 )  
  138.   cvAddS( mosaicVertical1, cvScalarAll( (-1.0)*min_val ), mosaicVertical1, NULL );  
  139.  mosaic2 = cvCreateImage(cvSize(mosaicVertical1->width, mosaicVertical1->height),  IPL_DEPTH_8U,1);  
  140.  cvConvertScale( mosaicVertical1, mosaicVertical1, 255.0/(max_val-min_val), 0 );  
  141.  cvConvertScaleAbs( mosaicVertical1, mosaic2, 1, 0 );  
  142.    
  143.  //  cvSaveImage("DOGPyramid of me.jpg",mosaic2);  
  144.  cvNamedWindow("mosaic1",1);  
  145.  cvShowImage("mosaic1", mosaic2);  
  146.  cvWaitKey(0);  
  147.    
  148.  //SIFT演算法第三步：特徵點位置檢測，最後確定特徵點的位置  
  149.  int keycount=DetectKeypoint(numoctaves, Gaussianpyr);  
  150.  printf("the keypoints number are %d ;/n", keycount);  
  151.  cvCopy(src,image1,NULL);  
  152.  DisplayKeypointLocation( image1 ,Gaussianpyr);  
  153.    
  154.  cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
  155.  cvNamedWindow("image1",1);  
  156.  cvShowImage("image1", DoubleSizeImage);  
  157.  cvWaitKey(0);    
  158.  cvDestroyWindow("image1");  
  159.    
  160.  //SIFT演算法第四步：計算高斯影象的梯度方向和幅值，計算各個特徵點的主方向  
  161.  ComputeGrad_DirecandMag(numoctaves, Gaussianpyr);  
  162.  AssignTheMainOrientation( numoctaves, Gaussianpyr,mag_pyr,grad_pyr);  
  163.  cvCopy(src,image1,NULL);  
  164.  DisplayOrientation ( image1, Gaussianpyr);  
  165.    
  166.  //  cvPyrUp( image1, DoubleSizeImage, CV_GAUSSIAN_5x5 );  
  167.  cvNamedWindow("image1",1);  
  168.  //  cvResizeWindow("image1", 2*(image1->width), 2*(image1->height) );  
  169.  cvShowImage("image1", image1);  
  170.  cvWaitKey(0);  
  171.    
  172.  //SIFT演算法第五步：抽取各個特徵點處的特徵描述字  
  173.  ExtractFeatureDescriptors( numoctaves, Gaussianpyr);  
  174.  cvWaitKey(0);  
  175.    
  176.  //銷燬視窗  
  177.  cvDestroyWindow("image1");  
  178.  cvDestroyWindow("mosaic1");  
  179.  //釋放影象  
  180.  cvReleaseImage(&image1);  
  181.  cvReleaseImage(&grey_im1);  
  182.  cvReleaseImage(&mosaic1);  
  183.  cvReleaseImage(&mosaic2);  
  184.  return 0;  
  185. }  

 

最後，再看一下，執行效果（圖中美女為老鄉+朋友，何姐08年照）：

完。 

updated

    有很多朋友都在本文評論下要求要本程式的完整原始碼包（注：本文程式碼未貼全，複製貼上編譯肯定諸多錯誤），但由於時隔太久，這份程式碼我自己也找不到了，不過，我可以提供一份sift + KD + BBF，且可以編譯正確的程式碼供大家參考學習，有pudn帳號的朋友可以前去下載：http://www.pudn.com/downloads340/sourcecode/graph/texture_mapping/detail1486667.html（沒有pudn賬號的同學請加群：169056165，驗證資訊：sift，至群共享下載），然後用兩幅不同的圖片做了下匹配（當然，執行結果顯示是不匹配的），效果還不錯：http://weibo.com/1580904460/yDmzAEwcV#1348475194313! July、二零一二年十月十一日。


相關文章