1.5 SIFT特征提取與描述

人類在識別一個物體時，不管這個物體是遠還是近，都能對它進行正確的辨識，這就是所謂的“尺度不變性”。同樣，當這個物體發生旋轉時，我們照樣可以正確地識別它，這就是所謂的“旋轉不變性”……那么，如何讓機器與人類一樣具有這種能力呢？這就是圖像局部不變性特征要解決的問題。對于圖像不變性特征提取的方法，核心是“不變性”三個字。

提到圖像局部不變性特征，有兩個人不得不提及，一個是Lindeberg，另一個是Lowe。如果將局部不變性特征方法比作一個孩子，那么，Lindeberg就是這個孩子的父親，Lowe是母親。Lindeberg奠定了局部不變性特征方法的理論基礎，播下了局部不變性特征方法的種子，而Lowe則將這顆種子孕育成一種能具體實現的方法。由于在Lindeberg的尺度空間理論中，各種高斯微分算子與哺乳動物的視網膜和視覺皮層的感受域剖面有著高度的相似性，所以，尺度空間理論經常與生物視覺相關聯，有人也稱圖像局部不變性特征方法為基于生物視覺的不變性方法。

1.5.1 SIFT算法

Lowe總結了現有的特征提取方法，并提出了一種新型高效的特征檢測描述方法——尺度不變特征變換（Scale Invariant Feature Transform, SIFT）。SIFT的主要思路是：首先建立圖像的尺度空間表示，然后在尺度空間中搜索圖像的極值點，由極值點再建立特征描述向量，最后用特征描述向量進行相似度匹配。采用SIFT方法提取的圖像特征具有放縮不變性、旋轉不變性、仿射不變性，還有一定的抗光照變化和抗視點變換性能。該特征還具有高度的可區分性，能夠在一個具有大量特征數據的數據庫中精確的匹配。SIFT還使用金字塔算法大大縮小了提取特征時的運算量。計算方法分為四個步驟：尺度空間極值提取、特征點定位、特征方向賦值和提取特征點描述。

1．尺度空間極值檢測

尺度空間理論是檢測不變性特征的基礎。Witkin（1983）提出了尺度空間理論，他主要討論了一維信號平滑處理的問題。Koenderink（1984）把這種理論擴展到二維圖像，并證明高斯卷積核是實現尺度變換的唯一變換核。

一幅二維圖像在不同尺度下的尺度空間表示可由圖像與高斯核卷積得到：

其中，G（x,y,σ）為高斯核函數：

其中，（x,y）為圖像點的像素坐標，I（x,y）為圖像數據。σ稱為尺度空間因子，它是高斯正態分布的方差，反映了圖像被平滑的程度，其值越小表征圖像被平滑程度越小，相應尺度也越小。L（x,y,σ）代表了圖像的尺度空間。

為高效地在尺度空間內檢測出穩定的特征點，Lowe使用尺度空間中差分高斯（Difference of Gaussian, DoG）極值作為判斷依據。DoG算子定義為兩個不同尺度的高斯核的差分，是歸一化（Laplacian of Gaussian, LoG）算子的近似。設k為兩個相鄰尺度間的比例因子，則DoG算子定義如下：

D（x,y,σ）構造方式如圖1.5.1所示，其建立高斯圖像（圖1.5.2中左列）與DoG（圖1.5.2中右列）兩個金字塔。高斯圖像金字塔分為多組，每組間又分為多層。一組中的多層間不同的是尺度，相鄰層間尺度相差一個比例因子k。在S個尺度間隔內變化尺度因子，如使加倍，則k應為。為了在整個金字塔內獲取DoG極值，應在高斯金字塔中生成S+3層高斯平滑圖像。下一組圖像的最底層由上一組中尺度為2σ的圖像進行因子為2的降采樣得到，其中σ為上一組中最底層圖像的尺度因子。DoG金字塔由相鄰的高斯圖像金字塔相減得到。所生成的高斯圖像金字塔與DoG金字塔如圖1.5.2所示。

金字塔中每個高斯圖像的σ為：

其中，σ0為基礎尺度因子；o、s分別為圖像所在的圖像組坐標、組內層坐標，o∈omin+［0, …,O－1］,s∈［0, …,S－1］;omin是第一個金字塔組的坐標，通常omin取0或者－1，當設為－1的時候，則圖像在計算高斯尺度空間前先擴大一倍。在Lowe的算法實現中，以上參數的取值為：σ0＝1.6×, omin＝－1,S＝3。

此外，空間坐標x是組坐標o 的函數。設x0為第0組內的空間坐標，則有：x＝2ox0,x0∈［0, …,N0－1］×［0, …,M0－1］，其中（N0,M0）是第0組中圖像的分辨率。設（N0,M0）為第0組中的圖像分辨率，則其他組的分辨率為：

金字塔構造完后開始檢測DoG局部極值。其中每個像素需要跟同一尺度的周圍鄰域8個像素和相鄰尺度對應位置的周圍鄰域9×2個像素總共26個像素進行比較，如圖1.5.3所示。僅當被檢測點的DoG值大于此26個像素點或小于此26個像素點時才將該點判定為極值點并保存以進行后續計算。

2．確定關鍵點位置及尺度

通過擬合三維二次函數以精確確定關鍵點的位置和尺度，得到原圖像SIFT候選特征點集合，如圖1.5.4所示。

在得到原圖像的SIFT候選特征點集合X0后，要從中篩選穩定的點作為該圖像的SIFT特征點，其組成的集合用X表示。由于X0中對比度較低的點對噪聲較敏感，而位于邊緣上的點難以準確定位，為確保SIFT特征點的穩定性，必須將這兩種點剔除。

1）剔除對比度低的點

將候選特征點x的偏移量定義為Δx，其對比度為D（x）的絕對值D（x），對x的DoG函數表達式（1.5.3）進行泰勒級數展開：

式中，由于x為DoG函數的極值點，所以＝0，解方程得Δx＝，通過多次迭代得到最終候選點的精確位置及尺度，將其代入式（1.5.6）求得D（），求其絕對值可得。設對比度閾值為Tc，則低對比度點的剔除公式為：

2）剔除邊緣點

由于邊緣梯度方向上主曲率值較大，而邊緣方向上曲率較小，在邊緣上得到的DoG函數的極值點與非邊緣區域的點相比，其主曲率比值較大，因此可以將主曲率比值大于一定閾值的點視為位于邊緣上的點將其剔除。

由于候選點的DoG函數D（x）的主曲率與2×2的Hessian矩陣H的特征值成正比：

式（1.5.8）中，Dxx、Dxy、Dyy為候選點鄰域對應位置的像素差分。令α為H 的最大特征值，β為H的最小特征值，令γ＝，則D（x）的主曲率比值與γ成正比。由H的跡和行列式的值可得：

只與兩特征值之比有關，與特征值自身大小無關，當兩特征值相等時最小，且隨著γ的增大而增大。設主曲率比值閾值為Tγ，則邊緣點的剔除公式為：

剔除低對比度和邊緣點后所剩的特征點如圖1.5.5所示。

3）關鍵點方向確定

通過確定關鍵點的方向，可以使特征描述符以與方向相關的方式構造，從而使算子具有旋轉不變性。關鍵點的方向利用其鄰域像素的梯度分布特性確定。對每個高斯圖像，每個點L（x,y）的梯度的模m（x,y）與方向θ（x,y）可以通過下式計算得到：

其中，L（x,y）所用尺度為關鍵點所在尺度。

1.5.2 SIFT特征描述

對每個關鍵點，在以其為中心的鄰域窗口內利用直方圖的方式統計鄰域像素的梯度分布。此直方圖有36個柱，每柱10°，共360°。每個加入直方圖的鄰域像素樣本的權重由該像素的梯度模與高斯權重確定，此高斯窗的σ為關鍵點的尺度的1.5倍。加入高斯窗的目的是增強離關鍵點近的鄰域點對關鍵點的影響。

直方圖的峰值反映了關鍵點所處鄰域梯度的主方向。完成直方圖統計后，找到直方圖的最高峰值以確定關鍵點的方向。關鍵點的方向可以由離最高峰值最近的三個柱值通過拋物線插值精確得到。

在梯度方向直方圖中，當存在一個大于等于主峰值80%能量的峰值時，則添加一個新關鍵點，此關鍵點的坐標、尺度與當前關鍵點相同，但方向為由此峰值確定的方向。因此一個關鍵點可能產生多個坐標、尺度相同，方向不同的關鍵點。這樣做的目的是增強匹配的魯棒性。

至此，特征點檢測完畢，特征描述前的準備工作已經完成。每個關鍵點含有三個信息：坐標、尺度、方向。

首先將坐標軸旋轉為關鍵點的方向，以確保旋轉不變性。

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

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

綜上所述，提取圖像SIFT特征點及其特征向量的流程如圖1.5.7所示。

1.5.3 實例精講

【例1.5.1】 SIFT算法實現的MATLAB源程序，讀者可以根據程序及相關的注釋對SIFT作進一步的理解。

        ********************************************************************
        function[pos, scale, orient, desc]＝ SIFT(im, octaves, intervals, object＿mask, contrast
        ＿threshold, curvature＿threshold, interactive)
        % 功能：提取灰度圖像的尺度不變特征(SIFT特征)
        % 輸入：
        % im——灰度圖像，該圖像的灰度值在0～1范圍內
        % octaves——金字塔的組數：octaves(默認值為4)
        % intervals——該輸入參數決定每組金字塔的層數(默認值為 2)
        % object＿mask——確定圖像中尺度不變特征點的搜索區域
        % contrast＿threshold——對比度閾值(默認值為0.03)
        % curvature＿threshold——曲率閾值(默認值為10.0)
        % interactive——函數運行顯示標志，將其設定為1，則顯示算法運行時間和過程的相關信息
        %               如果將其設定為2，則僅顯示最終運行結果(default＝ 1)
        % 輸出：
        % pos——Nx2 矩陣，每一行包括尺度不變特征點的坐標(x, y)
        % scale——Nx3 矩陣，每一行包括尺度不變特征點的尺度信息
        % orient——Nx1 向量，每個元素是特征點的主方向，其范圍為[－pi, pi)
        % desc——Nx128 矩陣，每一行包含特征點的特征向量
        % 參考文獻：
        %[1]David G.Lowe, ＂Distinctive Image Features from Sacle－Invariant Keypoints＂,
        %     accepted for publication in the International Journal of Computer
        %     Vision,2004.
        %[2]David G.Lowe, ＂Object Recognition from Local Scale－Invariant Features＂,
        %     Proc.of the International Conference on Computer Vision, Corfu,
        %     September 1999.
        % 設定輸入量的默認值
        if ～exist(＇octaves＇)
          octaves ＝ 4;
        end
        if ～exist(＇intervals＇)
          intervals ＝ 2;
        end
        if～exist(＇object＿mask＇)
          object＿mask＝ ones(size(im));
        end
        if size(object＿mask)～＝ size(im)
          object＿mask＝ ones(size(im));
        end
        if～exist(＇contrast＿threshold＇)
          contrast＿threshold＝ 0.02;
        end
        if～exist(＇curvature＿threshold＇)
          curvature＿threshold＝ 10.0;
        end
        if ～exist(＇interactive＇)
          interactive ＝ 1;
        end
        % 檢驗輸入灰度圖像的像素灰度值是否已歸一化到[0,1]
        if((min(im(:))  0)｜(max(im(:))  1))
          fprintf(2, ＇Warning:image not normalized to[0,1].＼n＇);
        end
        % 將輸入圖像經過高斯平滑處理，采用雙線性差值將其擴大一倍
        if interactive  ＝ 1
          fprintf(2, ＇Doubling image size for first octave...＼n＇);
        end
        tic;
        antialias＿sigma＝ 0.5;
        if antialias＿sigma＝＝ 0
           signal ＝ im;
        else
          g＝ gaussian＿filter(antialias＿sigma);
          if exist(＇corrsep＇)＝＝ 3
              signal ＝ corrsep(g, g, im);
          else
              signal ＝ conv2(g, g, im, ＇same＇);
          end
        end
        signal ＝ im;
        [X Y]＝ meshgrid(1:0.5:size(signal,2),1:0.5:size(signal,1));
        signal ＝ interp2(signal, X, Y, ＇*linear＇);
        %  降采樣率；
        subsample ＝[0.5];
        %%%%%%下一步是生成高斯和差分高斯(DoG)金字塔%%%%%%
        %%%%%%高斯金字塔含有s+3層，差分高斯金字塔含有s+2層%%%%%%
        if interactive  ＝ 1
          fprintf(2, ＇Prebluring image...＼n＇);
        end
        preblur＿sigma＝ sqrt(sqrt(2)＾2 －(2*antialias＿sigma)＾2);
        if preblur＿sigma＝＝ 0
          gauss＿pyr{1,1}＝ signal;
        else
          g＝ gaussian＿filter(preblur＿sigma);
          if exist(＇corrsep＇)＝＝ 3
          gauss＿pyr{1,1}＝ corrsep(g, g, signal);
      else
          gauss＿pyr{1,1}＝ conv2(g, g, signal, ＇same＇);
      end
    end
    clear signal
    pre＿time＝ toc;
    if interactive  ＝ 1
      fprintf(2, ＇Preprocessing time %.2f seconds.＼n＇, pre＿time);
    end
    % 第一組第一層的sigma
    initial＿sigma＝ sqrt((2*antialias＿sigma)＾2 + preblur＿sigma ＾2);
    % 記錄每一層和每一個尺度的sigma
    absolute＿sigma＝ zeros(octaves, intervals+3);
    absolute＿sigma(1,1)＝ initial＿sigma* subsample(1);
    % 記錄產生金字塔的濾波器的尺寸和標準差
    filter＿size＝ zeros(octaves, intervals+3);
    filter＿sigma＝ zeros(octaves, intervals+3);
    % 生成高斯和差分高斯金字塔
    if interactive  ＝ 1
      fprintf(2, ＇Expanding the Gaussian and DOG pyramids...＼n＇);
    end
    tic;
    for octave ＝ 1:octaves
      if interactive  ＝ 1
          fprintf(2, ＼＇ tProcessing octave %d:image size %d x %d subsample %.1f＼n＇, octave,
    size(gauss＿pyr{octave,1},2), size(gauss＿pyr{octave,1},1), subsample(octave));
          fprintf(2, ＇＼t＼tInterval 1 sigma %f＼n＇, absolute＿sigma(octave,1));
      end
      sigma＝ initial＿sigma;
      g＝ gaussian＿filter(sigma);
      filter＿size(octave,1)＝ length(g);
      filter＿sigma(octave,1)＝ sigma;
      DOG＿pyr{octave}＝ zeros(size(gauss＿pyr{octave,1},1), size(gauss＿pyr{octave,1},2),
    intervals+2);
      for interval ＝ 2:(intervals+3)
          % 計算生成下一層幾何采樣金字塔的標準差
          % 其中，sigma＿i+1 ＝ k*sigma.
          %  sigma＿i+1 ＾2 ＝ sigma＿f, i ＾2 + sigma＿i ＾2
          %  (k*sigma＿i)＾2 ＝ sigma＿f, i ＾2 + sigma＿i ＾2
          % 因此
          %      sigma＿f, i＝ sqrt(k ∧2 － 1)sigma＿i
          %  對于擴展的組(span the octave), k＝ 2 ＾(1/intervals)
          % 所以
          %  sigma＿f, i＝ sqrt(2 ＾(2/intervals)－ 1)sigma＿i
          sigma＿f＝ sqrt(2 ＾(2/intervals)－ 1)*sigma;
          g＝ gaussian＿filter(sigma＿f);
          sigma＝(2 ＾(1/intervals))*sigma;
          % 記錄sigma的值
          absolute＿sigma(octave, interval)＝ sigma* subsample(octave);
          % 記錄濾波器的尺寸和標準差
          filter＿size(octave, interval)＝ length(g);
          filter＿sigma(octave, interval)＝ sigma;
          if exist(＇corrsep＇)＝＝ 3
            gauss＿pyr{octave, interval}＝ corrsep(g, g, gauss＿pyr{octave, interval－1});
          else
            gauss＿pyr{octave, interval}＝ conv2(g, g, gauss＿pyr{octave, interval－1}, ＇same＇);
          end
          DOG＿pyr{octave}(:, :, interval－1) ＝ gauss＿pyr{octave, interval} － gauss＿pyr
    {octave, interval－1};
          if interactive  ＝ 1
            fprintf(2, ＇＼t＼tInterval % d sigma % f＼n ＇, interval, absolute＿sigma(octave,
    interval));
          end
        end
        if octave  octaves
          sz ＝ size(gauss＿pyr{octave, intervals+1});
          [X Y]＝ meshgrid(1:2:sz(2),1:2:sz(1));
          gauss＿pyr{octave+1,1}＝ interp2(gauss＿pyr{octave, intervals+1}, X, Y, ＇*nearest＇);
          absolute＿sigma(octave+1,1)＝ absolute＿sigma(octave, intervals+1);
          subsample ＝[subsample subsample(end)*2];
        end
    end
    pyr＿time＝ toc;
    if interactive  ＝ 1
        fprintf(2, ＇Pryamid processing time %.2f seconds.＼n＇, pyr＿time);
    end
    % 在交互模式下顯示高斯金字塔
    if interactive  ＝ 2
      sz ＝ zeros(1,2);
      sz(2)＝(intervals+3)*size(gauss＿pyr{1,1},2);
      for octave ＝ 1:octaves
          sz(1)＝ sz(1)+ size(gauss＿pyr{octave,1},1);
      end
      pic ＝ zeros(sz);
      y＝ 1;
      for octave ＝ 1:octaves
          x ＝ 1;
          sz ＝ size(gauss＿pyr{octave,1});
          for interval ＝ 1:(intervals + 3)
                pic(y:(y+sz(1)－1), x:(x+sz(2)－1))＝ gauss＿pyr{octave, interval};
            x ＝ x + sz(2);
          end
          y＝ y+ sz(1);
      end
      fig ＝ figure;
      clf;
      imshow(pic);
      resizeImageFig(fig, size(pic),0.25);
      fprintf(2, ＇The gaussian pyramid(0.25 scale).＼nPress any key to continue.＼n＇);
      pause;
      close(fig)
    end
    % 在交互模式下顯示差分高斯金字塔
    if interactive  ＝ 2
      sz ＝ zeros(1,2);
      sz(2)＝(intervals+2)*size(DOG＿pyr{1}(:, :,1),2);
      for octave ＝ 1:octaves
          sz(1)＝ sz(1)+ size(DOG＿pyr{octave}(:, :,1),1);
      end
      pic ＝ zeros(sz);
      y＝ 1;
      for octave ＝ 1:octaves
          x ＝ 1;
          sz ＝ size(DOG＿pyr{octave}(:, :,1));
          for interval ＝ 1:(intervals + 2)
                pic(y:(y+sz(1)－1), x:(x+sz(2)－1))＝ DOG＿pyr{octave}(:, :, interval);
            x ＝ x + sz(2);
          end
          y＝ y+ sz(1);
      end
      fig ＝ figure;
      clf;
      imagesc(pic);
      resizeImageFig(fig, size(pic),0.25);
      fprintf(2, ＇The DOG pyramid(0.25 scale).＼nPress any key to continue.＼n＇);
      pause;
      close(fig)
    end
    % 下一步是查找差分高斯金字塔中的局部極值，并通過曲率和照度進行檢驗
    curvature＿threshold＝((curvature＿threshold+ 1)＾2)/curvature＿threshold;
    % 二階微分核
    xx ＝[1 －2  1];
    yy＝ xx＇;
    xy＝[1 0 －1;0 0 0; －1 0 1]/4;
    raw＿keypoints＝[];
    contrast＿keypoints＝[];
    curve＿keypoints＝[];
    % 在高斯金字塔中查找局部極值
    if interactive  ＝ 1
      fprintf(2, ＇Locating keypoints...＼n＇);
    end
    tic;
    loc＝ cell(size(DOG＿pyr));
    for octave ＝ 1:octaves
      if interactive  ＝ 1
          fprintf(2, ＇＼tProcessing octave %d＼n＇, octave);
      end
      for interval ＝ 2:(intervals+1)
          keypoint＿count＝ 0;
          contrast＿mask＝ abs(DOG＿pyr{octave}(:, :, interval)) ＝ contrast＿threshold;
          loc{octave, interval}＝ zeros(size(DOG＿pyr{octave}(:, :, interval)));
          if exist(＇corrsep＇)＝＝ 3
            edge ＝ 1;
          else
            edge＝ ceil(filter＿size(octave, interval)/2);
          end
          for y＝(1+edge):(size(DOG＿pyr{octave}(:, :, interval),1)－edge)
            for x＝(1+edge):(size(DOG＿pyr{octave}(:, :, interval),2)－edge)
              if object＿mask(round(y*subsample(octave)), round(x*subsample(octave)))＝＝ 1
                  if((interactive ＝ 2)｜(contrast＿mask(y, x)＝＝ 1))
                      % 通過空間核尺度檢測最大值和最小值
                    DOG＿pyr{octave}((y－1):(y+1), (x－1):(x+1), (interval－1):(interval+1));
                    pt＿val ＝ tmp(2,2,2);
                    if((pt＿val ＝＝ min(tmp(:)))｜(pt＿val ＝＝ max(tmp(:))))
                        % 存儲對灰度大于對比度閾值的點的坐標
                        raw＿keypoints＝[raw＿keypoints; x*subsample(octave)y*
    subsample(octave)];
                        if abs(DOG＿pyr{octave}(y, x, interval)) ＝ contrast＿threshold
                          contrast＿keypoints＝[contrast＿keypoints; raw＿keypoints(end, :)];
            % 計算局部極值的Hessian矩陣
            Dxx＝ sum(DOG＿pyr{octave}(y, x－1:x+1, interval).* xx);
            Dyy＝ sum(DOG＿pyr{octave}(y－1:y+1, x, interval).* yy);
            Dxy＝ sum(sum(DOG＿pyr{octave}(y－1:y+1, x－1:x+1, interval).* xy));
                          % 計算Hessian矩陣的直跡和行列式
                          Tr＿H＝ Dxx+ Dyy;
                          Det＿H＝ Dxx*Dyy－ Dxy ＾2;
                          % 計算主曲率
                          curvature＿ratio＝(Tr＿H ＾2)/Det＿H;
                    if((Det＿H ＝ 0)&(curvature＿ratio  curvature＿threshold))
                  % 存儲主曲率小于閾值的極值點的坐標(非邊緣點)
                  curve＿keypoints＝[curve＿keypoints; raw＿keypoints(end, :)];
                              % 將該點的位置的坐標設為 1，并計算點的數量
                              loc{octave, interval}(y, x)＝ 1;
                              keypoint＿count＝ keypoint＿count+ 1;
                          end
                        end
                    end
                  end
              end
            end
          end
          if interactive  ＝ 1
            fprintf(2, ＇＼t＼t%d keypoints found on interval %d＼n＇, keypoint＿count, interval);
          end
        end
    end
    keypoint＿time＝ toc;
    if interactive  ＝ 1
        fprintf(2, ＇Keypoint location time %.2f seconds.＼n＇, keypoint＿time);
    end
    % 在交互模式下顯示特征點檢測的結果
    if interactive  ＝ 2
        fig ＝ figure;
        clf;
        imshow(im);
        hold on;
        plot(raw＿keypoints(:,1), raw＿keypoints(:,2), ＇y+ ＇);
        resizeImageFig(fig, size(im),2);
        fprintf(2, ＇DOG extrema(2x scale).＼nPress any key to continue.＼n＇);
        pause;
        close(fig);
        fig ＝ figure;
        clf;
        imshow(im);
        hold on;
        plot(contrast＿keypoints(:,1), contrast＿keypoints(:,2), ＇y+ ＇);
        resizeImageFig(fig, size(im),2);
        fprintf(2, ＇Keypoints after removing low contrast extrema(2x scale).＼nPress any key to
    continue.＼n＇);
        pause;
        close(fig);
        fig ＝ figure;
        clf;
        imshow(im);
        hold on;
        plot(curve＿keypoints(:,1), curve＿keypoints(:,2), ＇y+ ＇);
        resizeImageFig(fig, size(im),2);
        fprintf(2, ＇Keypoints after removing edge points using principal curvature filtering(2x
    scale).＼nPress any key to continue.＼n＇);
        pause;
        close(fig);
    end
    clear raw＿keypoints contrast＿keypoints curve＿keypoints
    %%%%%% 下一步是計算特征點的主方向%%%%%%
    % 在特征點的一個區域內計算其梯度直方圖
    g＝ gaussian＿filter(1.5 * absolute＿sigma(1, intervals+3)/subsample(1));
    zero＿pad＝ ceil(length(g)/2);
    % 計算高斯金字塔圖像的梯度方向和幅值
    if interactive  ＝ 1
      fprintf(2, ＇Computing gradient magnitude and orientation...＼n＇);
    end
    tic;
    mag＿thresh＝ zeros(size(gauss＿pyr));
    mag＿pyr＝ cell(size(gauss＿pyr));
    grad＿pyr＝ cell(size(gauss＿pyr));
    for octave ＝ 1:octaves
      for interval ＝ 2:(intervals+1)
          % 計算x, y的差分
          diff＿x ＝ 0.5*(gauss＿pyr{octave, interval}(2:(end－1),3:(end))－gauss＿pyr
    {octave, interval}(2:(end－1),1:(end－2)));
          diff＿y ＝ 0.5*(gauss＿pyr{octave, interval}(3:(end),2:(end－1))－gauss＿pyr
    {octave, interval}(1:(end－2),2:(end－1)));
          % 計算梯度幅值
          mag＝ zeros(size(gauss＿pyr{octave, interval}));
          mag(2:(end－1),2:(end－1))＝ sqrt(diff＿x.＾ 2 + diff＿y.＾ 2);
          % 存儲高斯金字塔梯度幅值
          mag＿pyr{octave, interval}＝ zeros(size(mag)+2*zero＿pad);
    mag＿pyr{octave, interval}((zero＿pad+1):(end－zero＿pad), (zero＿pad+1):(end－zero＿pad))
    ＝ mag;
          % 計算梯度主方向
          grad＝ zeros(size(gauss＿pyr{octave, interval}));
          grad(2:(end－1),2:(end－1))＝ atan2(diff＿y, diff＿x);
          grad(find(grad ＝＝ pi))＝ －pi;
          % 存儲高斯金字塔梯度主方向
          grad＿pyr{octave, interval}＝ zeros(size(grad)+2*zero＿pad);
    grad＿pyr{octave, interval}((zero＿pad+1):(end－zero＿pad), (zero＿pad+1):(end－zero＿pad))
    ＝ grad;
      end
    end
    clear mag grad
    grad＿time＝ toc;
    if interactive  ＝ 1
      fprintf(2, ＇Gradient calculation time %.2f seconds.＼n＇, grad＿time);
    end
    %%%%%%下一步是確定特征點的主方向%%%%%%
    %%%%%% 方法：通過尋找每個關鍵點的子區域內梯度直方圖的峰值來解決%%%%%%
    % 將灰度直方圖分為36等分，每隔 10度一份
    num＿bins ＝ 36;
    hist＿step＝ 2*pi/num＿bins;
    hist＿orient＝[－pi:hist＿step:(pi－hist＿step)];
    % 初始化關鍵點的位置、方向和尺度信息
    pos ＝[];
    orient ＝[];
    scale ＝[];
    % 給關鍵點確定主方向
    if interactive  ＝ 1
      fprintf(2, ＇Assigining keypoint orientations...＼n＇);
    end
    tic;
    for octave ＝ 1:octaves
      if interactive  ＝ 1
          fprintf(2, ＇＼tProcessing octave %d＼n＇, octave);
      end
      for interval ＝ 2:(intervals + 1)
          if interactive  ＝ 1
            fprintf(2, ＇＼t＼tProcessing interval %d ＇, interval);
          end
          keypoint＿count＝ 0;
          % 構造高斯加權掩模
          g＝ gaussian＿filter(1.5 * absolute＿sigma(octave, interval)/subsample(octave));
          hf＿sz ＝ floor(length(g)/2);
          g ＝ g＇*g;
          loc＿pad＝ zeros(size(loc{octave, interval})+2*zero＿pad);
          loc＿pad((zero＿pad+1):(end－zero＿pad), (zero＿pad+1):(end－zero＿pad))＝
    loc{octave, interval};
          [iy ix]＝find(loc＿pad＝＝1);
          for k ＝ 1:length(iy)
            x ＝ ix(k);
            y＝ iy(k);
            wght＝g.*mag＿pyr{octave, interval}((y－hf＿sz):(y+hf＿sz), (x－hf＿sz):(x+hf＿sz));
            grad＿window＝grad＿pyr{octave, interval}((y－hf＿sz):(y+hf＿sz), (x－hf＿sz):(x+
    hf＿sz));
            orient＿hist＝zeros(length(hist＿orient),1);
            for bin＝1:length(hist＿orient)
              diff＝ mod(grad＿window－ hist＿orient(bin)+ pi,2*pi)－ pi;
              orient＿hist(bin)＝orient＿hist(bin)+sum(sum(wght.*max(1 －abs(diff)/hist＿
    step,0)));
            end
            % 運用非極大抑制法查找主方向直方圖的峰值
            peaks＝ orient＿hist;
            rot＿right＝[peaks(end); peaks(1:end－1)];
            rot＿left＝[peaks(2:end); peaks(1)];
            peaks(find(peaks  rot＿right))＝ 0;
            peaks(find(peaks  rot＿left))＝ 0;
            % 提取最大峰值的值和其索引位置
            [max＿peak＿val ipeak]＝ max(peaks);
            % 將大于等于最大峰值80%的直方圖也確定為特征點的主方向
            peak＿val ＝ max＿peak＿val;
            while(peak＿val  0.8*max＿peak＿val)
                % 最高峰值最近的三個柱值通過拋物線插值精確得到
              A ＝[];
              b ＝[];
              for j ＝ －1:1
                  A＝[A; (hist＿orient(ipeak)+hist＿step*j).＾2(hist＿orient(ipeak)+
    hist＿step*j)1];
                  bin＝ mod(ipeak+ j + num＿bins－ 1, num＿bins)+ 1;
                  b＝[b; orient＿hist(bin)];
              end
              c ＝ pinv(A)*b;
              max＿orient＝ －c(2)/(2*c(1));
              while(max＿orient  －pi)
                  max＿orient＝ max＿orient+ 2*pi;
              end
              while(max＿orient ＝ pi)
                  max＿orient＝ max＿orient－ 2*pi;
              end
              % 存儲關鍵點的位置、主方向和尺度信息
              pos＝[pos; [(x－zero＿pad)(y－zero＿pad)]*subsample(octave)];
              orient＝[orient; max＿orient];
              scale＝[scale; octave interval absolute＿sigma(octave, interval)];
              keypoint＿count＝ keypoint＿count+ 1;
              peaks(ipeak)＝ 0;
              [peak＿val ipeak]＝ max(peaks);
            end
        end
        if interactive  ＝ 1
            fprintf(2, ＇(%d keypoints)＼n＇, keypoint＿count);
        end
      end
  end
  clear loc loc＿pad
  orient＿time＝ toc;
  if interactive  ＝ 1
      fprintf(2, ＇Orientation assignment time %.2f seconds.＼n＇, orient＿time);
  end
  % 在交互模式下顯示關鍵點的尺度和主方向信息
  if interactive  ＝ 2
      fig ＝ figure;
      clf;
      imshow(im);
      hold on;
      display＿keypoints(pos, scale(:,3), orient, ＇y＇);
      resizeImageFig(fig, size(im),2);
      fprintf(2, ＇Final keypoints with scale and orientation(2x scale).＼nPress any key to
  continue.＼n＇);
      pause;
      close(fig);
  end
  %%%%%% SIFT算法的最后一步是特征向量生成%%%%%%
  orient＿bin＿spacing＝ pi/4;
  orient＿angles＝[－pi:orient＿bin＿spacing:(pi－orient＿bin＿spacing)];
  grid＿spacing＝ 4;
  [x＿coords y＿coords]＝ meshgrid([－6:grid＿spacing:6]);
  feat＿grid＝[x＿coords(:)y＿coords(:)]＇;
  [x＿coords y＿coords]＝ meshgrid([－(2*grid＿spacing－0.5):(2*grid＿spacing－0.5)]);
  feat＿samples＝[x＿coords(:)y＿coords(:)]＇;
  feat＿window＝ 2*grid＿spacing;
  desc ＝[];
  if interactive  ＝ 1
    fprintf(2, ＇Computing keypoint feature descriptors for %d keypoints＇, size(pos,1));
  end
  for k ＝ 1:size(pos,1)
    x ＝ pos(k,1)/subsample(scale(k,1));
    y＝ pos(k,2)/subsample(scale(k,1));
    % 將坐標軸旋轉為關鍵點的方向，以確保旋轉不變性
    M＝[cos(orient(k))－sin(orient(k)); sin(orient(k))cos(orient(k))];
    feat＿rot＿grid＝M*feat＿grid+ repmat([x; y],1, size(feat＿grid,2));
    feat＿rot＿samples＝M*feat＿samples+ repmat([x; y],1, size(feat＿samples,2));
    % 初始化特征向量
    feat＿desc＝ zeros(1,128);
    for s＝ 1:size(feat＿rot＿samples,2)
        x＿sample＝ feat＿rot＿samples(1, s);
        y＿sample＝ feat＿rot＿samples(2, s);
        % 在采樣位置進行梯度插值
        [X Y]＝ meshgrid((x＿sample－1):(x＿sample+1), (y＿sample－1):(y＿sample+1));
        G＝ interp2(gauss＿pyr{scale(k,1), scale(k,2)}, X, Y, ＇*linear＇);
        G(find(isnan(G)))＝ 0;
        diff＿x＝ 0.5*(G(2,3)－ G(2,1));
        diff＿y＝ 0.5*(G(3,2)－ G(1,2));
        mag＿sample＝ sqrt(diff＿x ＾2 + diff＿y ＾2);
        grad＿sample＝ atan2(diff＿y, diff＿x);
        if grad＿sample＝＝ pi
          grad＿sample＝ －pi;
        end
        % 計算x、y方向上的權重
        x＿wght＝ max(1 －(abs(feat＿rot＿grid(1, :)－ x＿sample)/grid＿spacing),0);
        y＿wght＝ max(1 －(abs(feat＿rot＿grid(2, :)－ y＿sample)/grid＿spacing),0);
        pos＿wght＝ reshape(repmat(x＿wght.*y＿wght,8,1),1,128);
        diff＝ mod(grad＿sample－ orient(k)－ orient＿angles+ pi,2*pi)－ pi;
        orient＿wght＝ max(1 － abs(diff)/orient＿bin＿spacing,0);
        orient＿wght＝ repmat(orient＿wght,1,16);
        % 計算高斯權重
        g＝ exp(－((x＿sample－x)＾2+(y＿sample－y)＾2)/(2*feat＿window ＾2))/(2*pi*feat＿
  window ＾2);
        feat＿desc＝ feat＿desc+ pos＿wght.*orient＿wght*g*mag＿sample;
    end
    % 將特征向量的長度歸一化，則可以進一步去除光照變化的影響
    feat＿desc＝ feat＿desc/norm(feat＿desc);
    feat＿desc(find(feat＿desc  0.2))＝ 0.2;
    feat＿desc＝ feat＿desc/norm(feat＿desc);
    % 存儲特征向量
    desc＝[desc; feat＿desc];
    if(interactive  ＝ 1)&(mod(k,25)＝＝ 0)
        fprintf(2, ＇.＇);
    end
  end
  desc＿time＝ toc;
  % 調整采樣偏差
  sample＿offset＝ －(subsample－ 1);
  for k ＝ 1:size(pos,1)
    pos(k, :)＝ pos(k, :)+ sample＿offset(scale(k,1));
  end
  if size(pos,1)  0
      scale ＝ scale(:,3);
  end
  % 在交互模式下顯示運行過程耗時
  if interactive  ＝ 1
     fprintf(2, ＇＼nDescriptor processing time %.2f seconds.＼n＇, desc＿time);
     fprintf(2, ＇Processing time summary:＼n＇);
     fprintf(2, ＇＼tPreprocessing:＼t%.2f s＼n＇, pre＿time);
     fprintf(2, ＇＼tPyramid:＼t%.2f s＼n＇, pyr＿time);
     fprintf(2, ＇＼tKeypoints:＼t%.2f s＼n＇, keypoint＿time);
     fprintf(2, ＇＼tGradient:＼t%.2f s＼n＇, grad＿time);
     fprintf(2, ＇＼tOrientation:＼t%.2f s＼n＇, orient＿time);
     fprintf(2, ＇＼tDescriptor:＼t%.2f s＼n＇, desc＿time);
     fprintf(2, ＇Total processing time %.2f seconds.＼n＇, pre＿time+ pyr＿time+ keypoint＿
  time+ grad＿time+ orient＿time+ desc＿time);
  end
  ********************************************************************

在上述實現SIFT算法的MATLAB程序中，調用了其他一些函數如下：

    ********************************************************************
    function g＝ gaussian＿filter(sigma)
    % 功能：生成一維高斯濾波器
    % 輸入：
    % sigma——高斯濾波器的標準差
    % 輸出：
    % g——高斯濾波器
    sample ＝ 7.0/2.0;
    n ＝ 2*round(sample * sigma)+1;
    x＝1:n;
    x＝x－ceil(n/2);
    g＝ exp(－(x.＾2)/(2*sigma ＾2))/(sigma*sqrt(2*pi));
    ********************************************************************
    ********************************************************************
    function im ＝ pgmread(fname);
    % 功能：讀入．pgm圖像
    global SIZE＿LIMIT;
    if SIZE＿LIMIT
    fprintf(1, ＇SIZE＿LIMIT recognized by pgmread＼n＇);
    end
    [fid, msg]＝ fopen(fname, ＇r＇);
    if(fid ＝＝ －1)
      error(msg);
    end
    [pars type]＝ pnmReadHeader(fid);
    if(pars＝＝－1)
      fclose(fid);
      error([fname ＇:cannot parse pgm header＇]);
    end
    if ～(type ＝＝ ＇P2 ＇｜type ＝＝ ＇P5 ＇)
      fclose(fid);
      error([fname ＇:Not of type P2 or P5.＇]);
    end
    xdim ＝ pars(1);
    ydim ＝ pars(2);
    maxval ＝ pars(3);
    sz ＝ xdim * ydim;
    if SIZE＿LIMIT
      if sz  ＝ 16384
        ydim ＝ floor(16384/xdim);
        sz ＝ xdim * ydim;
        fprintf(1, ＇truncated image size:cols %d rows %d＼n＇, xdim, ydim)
      end
    end
    if(type ＝＝ ＇P2 ＇)
      [im, count]  ＝ fscanf(fid, ＇%d＇, sz);
    elseif(type ＝＝ ＇P5 ＇)
      count ＝ 0;
      im ＝[];
      stat ＝ fseek(fid, －sz, ＇eof＇);
      if ～stat
        [im, count]  ＝ fread(fid, sz, ＇uchar＇);
      end
      if(count ～＝ sz)
        fprintf(1, ＇Warning:File ended early! %s＼n＇, fname);
        fprintf(1, ＇...Padding with %d zeros.＼n＇, sz－count);
        im ＝[im ; zeros(sz－count,1)];
      end
    else
      fclose(fid);
      error([fname ＇:Not of type P2 or P5.＇]);
    end
    im ＝ reshape(im, xdim, ydim)＇;
    fclose(fid);
    ********************************************************************
    ********************************************************************
    function[pars, type]＝ pnmReadHeader(fid)
    % 功能：讀入并解析一個pnm格式圖像的頭文件
    pars ＝ －1;
    type ＝ ＇Unknown＇;
    TheLine ＝ fgetl(fid);
    szLine ＝ size(TheLine);
    endLine ＝ szLine(2);
    if(endLine  2)
      fprintf(1, [＇Unrecognized PNMfile type＼n＇]);
      pars ＝ －1;
      return;
    end
    type ＝ TheLine(1:2);
    ok ＝ 0;
    if(type(1)＝＝ ＇P＇)
      if(type(2)＝＝ ＇1 ＇｜type(2)＝＝ ＇2 ＇｜...
          type(2)＝＝ ＇3 ＇｜type(2)＝＝ ＇4 ＇｜...
          type(2)＝＝ ＇5 ＇｜type(2)＝＝ ＇6 ＇｜...
          type(2)＝＝ ＇B＇｜type(2)＝＝ ＇L＇)
        ok ＝ 1;
      else
        fprintf(1, [＇Unrecognized PNMfile type:＇type ＇＼n＇]);
      end
    end
    if(type(1)＝＝ ＇F＇)
      if(type(2)＝＝ ＇P＇｜type(2)＝＝ ＇U＇)
        ok ＝ 1;
      else
        fprintf(1, [＇Unrecognized PNMfile type:＇type ＇＼n＇]);
      end
    end
    if ～ok
      pars ＝ －1;
      return;
    end
    current ＝ 3;
    parIndex＝1;
    while(parIndex  4)
      while(current  endLine)
        TheLine ＝ fgetl(fid);
        if(TheLine ＝＝ －1)
          fprintf(1, ＇Unexpected EOF＼n＇);
          pars ＝ －1;
          return;
        end
        szLine ＝ size(TheLine);
        endLine ＝ szLine(2);
        current ＝ 1;
      end
      [token, count, errmsg, nextindex]＝ ...
          sscanf(TheLine(current:endLine), ＇%s＇,1);
      nextindex ＝ nextindex+current－1;
      if(count＝＝0)
        if(nextindex  endLine)
         current ＝ nextindex;
        else
         pars ＝ －1;
         fprintf(1, ＇Unexpected EOF＼n＇);
         return;
        end
      else
       if token(1)＝＝ ＇# ＇
        current ＝ endLine+1;
       else
        [pars(parIndex), count, errmsg, nextindex]＝ ...
          sscanf(TheLine(current:endLine), ＇%d＇,1);
        if ～(count＝＝1)
          fprintf(1, ＇Confused reading pgm header＼n＇);
          pars＝－1;
          return;
        end
        parIndex ＝ parIndex+1;
        current ＝ current+nextindex－1;
       end
      end
    end
    ********************************************************************
    ********************************************************************
    function resizeImageFig(h, sz, frac)
    % function resizeImageFig(h, sz, frac)
    % 功能：根據句柄重獲圖像的大小
    %  sz ＝ 圖像尺寸
    %  frac(默認值＝ 1)
    if(nargin  3)
      frac ＝ 1;
    end
    pos ＝ get(h, ＇Position＇);
    set(h, ＇Units＇, ＇pixels＇, ＇Position＇, ...
          [pos(1), pos(2)+pos(4)－frac*sz(1), ...
            frac*sz(2), frac*sz(1)]);
    set(gca, ＇Position＇, [0 0 1 1], ＇Visible＇, ＇off＇);
    ********************************************************************
    ********************************************************************
    function hh＝ display＿keypoints(pos, scale, orient, varargin)
    % 功能：在原始圖像上顯示特征點
    % 輸入：
    % pos——特征點的位置矩陣
    % scale——特征點的尺度矩陣
    % orient——特征點的主方向向量
    % 輸出：
    % hh——返回向量的線句柄
    hold on;
    alpha ＝ 0.33;
    beta ＝ 0.33;
    autoscale ＝ 1.5;
    plotarrows ＝ 1;
    sym ＝ ＇ ＇;
    filled ＝ 0;
    ls ＝ ＇－ ＇;
    ms ＝ ＇ ＇;
    col ＝ ＇ ＇;
    varin ＝ nargin － 3;
    while(varin  0)&isstr(varargin{varin}),
      vv ＝ varargin{varin};
      if ～isempty(vv)&strcmp(lower(vv(1)), ＇f＇)
          filled ＝ 1;
          nin ＝ nin－1;
      else
          [l, c, m, msg]＝ colstyle(vv);
          if ～isempty(msg),
            error(sprintf(＇Unknown option＂%s＂.＇, vv));
          end
          if ～isempty(l), ls ＝ l; end
          if ～isempty(c), col ＝ c; end
          if ～isempty(m), ms ＝ m; plotarrows ＝ 0; end
          if isequal(m, ＇.＇), ms ＝ ＇ ＇; end % Don＇t plot ＇.＇
          varin ＝ varin－1;
        end
    end
    if varin  0
        autoscale ＝ varargin{varin};
    end
    x ＝ pos(:,1);
    y＝ pos(:,2);
    u ＝ scale.*cos(orient);
    v ＝ scale.*sin(orient);
    if prod(size(u))＝＝1, u ＝ u(ones(size(x))); end
    if prod(size(v))＝＝1, v ＝ v(ones(size(u))); end
    if autoscale,
        u ＝ u*autoscale; v ＝ v*autoscale;
    end
    ax ＝ newplot;
    next ＝ lower(get(ax, ＇NextPlot＇));
    hold＿state＝ ishold;
    x ＝ x(:).＇; y＝ y(:).＇;
    u ＝ u(:).＇; v ＝ v(:).＇;
    uu ＝[x; x+u; repmat(NaN, size(u))];
    vv ＝[y; y+v; repmat(NaN, size(u))];
    h1 ＝ plot(uu(:), vv(:), [col ls]);
    if plotarrows,
        hu ＝[x+u－alpha*(u+beta*(v+eps)); x+u; ...
            x+u－alpha*(u－beta*(v+eps)); repmat(NaN, size(u))];
        hv ＝[y+v－alpha*(v－beta*(u+eps)); y+v; ...
            y+v－alpha*(v+beta*(u+eps)); repmat(NaN, size(v))];
        hold on
        h2 ＝ plot(hu(:), hv(:), [col ls]);
    else
        h2 ＝[];
    end
    if ～isempty(ms),
      hu ＝ x; hv ＝ y;
      hold on
      h3 ＝ plot(hu(:), hv(:), [col ms]);
      if filled, set(h3, ＇markerfacecolor＇, get(h1, ＇color＇)); end
    else
      h3 ＝[];
    end
    if～hold＿state, hold off, view(2); set(ax, ＇NextPlot＇, next); end
    if nargout  0, hh＝[h1; h2; h3]; end
    ********************************************************************