1.5 SIFT特征提取與描述

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

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

1.5.1 SIFT算法

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

1．尺度空間極值檢測

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

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

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

其中，（x,y）為圖像點(diǎn)的像素坐標(biāo)，I（x,y）為圖像數(shù)據(jù)。σ稱為尺度空間因子，它是高斯正態(tài)分布的方差，反映了圖像被平滑的程度，其值越小表征圖像被平滑程度越小，相應(yīng)尺度也越小。L（x,y,σ）代表了圖像的尺度空間。

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

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

金字塔中每個(gè)高斯圖像的σ為：

其中，σ0為基礎(chǔ)尺度因子；o、s分別為圖像所在的圖像組坐標(biāo)、組內(nèi)層坐標(biāo)，o∈omin+［0, …,O－1］,s∈［0, …,S－1］;omin是第一個(gè)金字塔組的坐標(biāo)，通常omin取0或者－1，當(dāng)設(shè)為－1的時(shí)候，則圖像在計(jì)算高斯尺度空間前先擴(kuò)大一倍。在Lowe的算法實(shí)現(xiàn)中，以上參數(shù)的取值為：σ0＝1.6×, omin＝－1,S＝3。

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

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

2．確定關(guān)鍵點(diǎn)位置及尺度

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

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

1）剔除對比度低的點(diǎn)

將候選特征點(diǎn)x的偏移量定義為Δx，其對比度為D（x）的絕對值D（x），對x的DoG函數(shù)表達(dá)式（1.5.3）進(jìn)行泰勒級數(shù)展開：

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

2）剔除邊緣點(diǎn)

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

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

式（1.5.8）中，Dxx、Dxy、Dyy為候選點(diǎn)鄰域?qū)?yīng)位置的像素差分。令α為H 的最大特征值，β為H的最小特征值，令γ＝，則D（x）的主曲率比值與γ成正比。由H的跡和行列式的值可得：

只與兩特征值之比有關(guān)，與特征值自身大小無關(guān)，當(dāng)兩特征值相等時(shí)最小，且隨著γ的增大而增大。設(shè)主曲率比值閾值為Tγ，則邊緣點(diǎn)的剔除公式為：

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

3）關(guān)鍵點(diǎn)方向確定

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

其中，L（x,y）所用尺度為關(guān)鍵點(diǎn)所在尺度。

1.5.2 SIFT特征描述

對每個(gè)關(guān)鍵點(diǎn)，在以其為中心的鄰域窗口內(nèi)利用直方圖的方式統(tǒng)計(jì)鄰域像素的梯度分布。此直方圖有36個(gè)柱，每柱10°，共360°。每個(gè)加入直方圖的鄰域像素樣本的權(quán)重由該像素的梯度模與高斯權(quán)重確定，此高斯窗的σ為關(guān)鍵點(diǎn)的尺度的1.5倍。加入高斯窗的目的是增強(qiáng)離關(guān)鍵點(diǎn)近的鄰域點(diǎn)對關(guān)鍵點(diǎn)的影響。

直方圖的峰值反映了關(guān)鍵點(diǎn)所處鄰域梯度的主方向。完成直方圖統(tǒng)計(jì)后，找到直方圖的最高峰值以確定關(guān)鍵點(diǎn)的方向。關(guān)鍵點(diǎn)的方向可以由離最高峰值最近的三個(gè)柱值通過拋物線插值精確得到。

在梯度方向直方圖中，當(dāng)存在一個(gè)大于等于主峰值80%能量的峰值時(shí)，則添加一個(gè)新關(guān)鍵點(diǎn)，此關(guān)鍵點(diǎn)的坐標(biāo)、尺度與當(dāng)前關(guān)鍵點(diǎn)相同，但方向?yàn)橛纱朔逯荡_定的方向。因此一個(gè)關(guān)鍵點(diǎn)可能產(chǎn)生多個(gè)坐標(biāo)、尺度相同，方向不同的關(guān)鍵點(diǎn)。這樣做的目的是增強(qiáng)匹配的魯棒性。

至此，特征點(diǎn)檢測完畢，特征描述前的準(zhǔn)備工作已經(jīng)完成。每個(gè)關(guān)鍵點(diǎn)含有三個(gè)信息：坐標(biāo)、尺度、方向。

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

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

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

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

1.5.3 實(shí)例精講

【例1.5.1】 SIFT算法實(shí)現(xiàn)的MATLAB源程序，讀者可以根據(jù)程序及相關(guān)的注釋對SIFT作進(jìn)一步的理解。

        ********************************************************************
        function[pos, scale, orient, desc]＝ SIFT(im, octaves, intervals, object＿mask, contrast
        ＿threshold, curvature＿threshold, interactive)
        % 功能：提取灰度圖像的尺度不變特征(SIFT特征)
        % 輸入：
        % im——灰度圖像，該圖像的灰度值在0～1范圍內(nèi)
        % octaves——金字塔的組數(shù)：octaves(默認(rèn)值為4)
        % intervals——該輸入?yún)?shù)決定每組金字塔的層數(shù)(默認(rèn)值為 2)
        % object＿mask——確定圖像中尺度不變特征點(diǎn)的搜索區(qū)域
        % contrast＿threshold——對比度閾值(默認(rèn)值為0.03)
        % curvature＿threshold——曲率閾值(默認(rèn)值為10.0)
        % interactive——函數(shù)運(yùn)行顯示標(biāo)志，將其設(shè)定為1，則顯示算法運(yùn)行時(shí)間和過程的相關(guān)信息
        %               如果將其設(shè)定為2，則僅顯示最終運(yùn)行結(jié)果(default＝ 1)
        % 輸出：
        % pos——Nx2 矩陣，每一行包括尺度不變特征點(diǎn)的坐標(biāo)(x, y)
        % scale——Nx3 矩陣，每一行包括尺度不變特征點(diǎn)的尺度信息
        % orient——Nx1 向量，每個(gè)元素是特征點(diǎn)的主方向，其范圍為[－pi, pi)
        % desc——Nx128 矩陣，每一行包含特征點(diǎn)的特征向量
        % 參考文獻(xiàn)：
        %[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.
        % 設(shè)定輸入量的默認(rèn)值
        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
        % 檢驗(yàn)輸入灰度圖像的像素灰度值是否已歸一化到[0,1]
        if((min(im(:))  0)｜(max(im(:))  1))
          fprintf(2, ＇Warning:image not normalized to[0,1].＼n＇);
        end
        % 將輸入圖像經(jīng)過高斯平滑處理，采用雙線性差值將其擴(kuò)大一倍
        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);
    % 記錄每一層和每一個(gè)尺度的sigma
    absolute＿sigma＝ zeros(octaves, intervals+3);
    absolute＿sigma(1,1)＝ initial＿sigma* subsample(1);
    % 記錄產(chǎn)生金字塔的濾波器的尺寸和標(biāo)準(zhǔn)差
    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)
          % 計(jì)算生成下一層幾何采樣金字塔的標(biāo)準(zhǔn)差
          % 其中，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
          %  對于擴(kuò)展的組(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);
          % 記錄濾波器的尺寸和標(biāo)準(zhǔn)差
          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
    % 下一步是查找差分高斯金字塔中的局部極值，并通過曲率和照度進(jìn)行檢驗(yàn)
    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(:))))
                        % 存儲(chǔ)對灰度大于對比度閾值的點(diǎn)的坐標(biāo)
                        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, :)];
            % 計(jì)算局部極值的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));
                          % 計(jì)算Hessian矩陣的直跡和行列式
                          Tr＿H＝ Dxx+ Dyy;
                          Det＿H＝ Dxx*Dyy－ Dxy ＾2;
                          % 計(jì)算主曲率
                          curvature＿ratio＝(Tr＿H ＾2)/Det＿H;
                    if((Det＿H ＝ 0)&(curvature＿ratio  curvature＿threshold))
                  % 存儲(chǔ)主曲率小于閾值的極值點(diǎn)的坐標(biāo)(非邊緣點(diǎn))
                  curve＿keypoints＝[curve＿keypoints; raw＿keypoints(end, :)];
                              % 將該點(diǎn)的位置的坐標(biāo)設(shè)為 1，并計(jì)算點(diǎn)的數(shù)量
                              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
    % 在交互模式下顯示特征點(diǎn)檢測的結(jié)果
    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
    %%%%%% 下一步是計(jì)算特征點(diǎn)的主方向%%%%%%
    % 在特征點(diǎn)的一個(gè)區(qū)域內(nèi)計(jì)算其梯度直方圖
    g＝ gaussian＿filter(1.5 * absolute＿sigma(1, intervals+3)/subsample(1));
    zero＿pad＝ ceil(length(g)/2);
    % 計(jì)算高斯金字塔圖像的梯度方向和幅值
    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)
          % 計(jì)算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)));
          % 計(jì)算梯度幅值
          mag＝ zeros(size(gauss＿pyr{octave, interval}));
          mag(2:(end－1),2:(end－1))＝ sqrt(diff＿x.＾ 2 + diff＿y.＾ 2);
          % 存儲(chǔ)高斯金字塔梯度幅值
          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;
          % 計(jì)算梯度主方向
          grad＝ zeros(size(gauss＿pyr{octave, interval}));
          grad(2:(end－1),2:(end－1))＝ atan2(diff＿y, diff＿x);
          grad(find(grad ＝＝ pi))＝ －pi;
          % 存儲(chǔ)高斯金字塔梯度主方向
          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
    %%%%%%下一步是確定特征點(diǎn)的主方向%%%%%%
    %%%%%% 方法：通過尋找每個(gè)關(guān)鍵點(diǎn)的子區(qū)域內(nèi)梯度直方圖的峰值來解決%%%%%%
    % 將灰度直方圖分為36等分，每隔 10度一份
    num＿bins ＝ 36;
    hist＿step＝ 2*pi/num＿bins;
    hist＿orient＝[－pi:hist＿step:(pi－h(huán)ist＿step)];
    % 初始化關(guān)鍵點(diǎn)的位置、方向和尺度信息
    pos ＝[];
    orient ＝[];
    scale ＝[];
    % 給關(guān)鍵點(diǎn)確定主方向
    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òu)造高斯加權(quán)掩模
          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－h(huán)f＿sz):(y+hf＿sz), (x－h(huán)f＿sz):(x+hf＿sz));
            grad＿window＝grad＿pyr{octave, interval}((y－h(huán)f＿sz):(y+hf＿sz), (x－h(huán)f＿sz):(x+
    hf＿sz));
            orient＿h(yuǎn)ist＝zeros(length(hist＿orient),1);
            for bin＝1:length(hist＿orient)
              diff＝ mod(grad＿window－ hist＿orient(bin)+ pi,2*pi)－ pi;
              orient＿h(yuǎn)ist(bin)＝orient＿h(yuǎn)ist(bin)+sum(sum(wght.*max(1 －abs(diff)/hist＿
    step,0)));
            end
            % 運(yùn)用非極大抑制法查找主方向直方圖的峰值
            peaks＝ orient＿h(yuǎn)ist;
            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%的直方圖也確定為特征點(diǎn)的主方向
            peak＿val ＝ max＿peak＿val;
            while(peak＿val  0.8*max＿peak＿val)
                % 最高峰值最近的三個(gè)柱值通過拋物線插值精確得到
              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＿h(yuǎn)ist(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
              % 存儲(chǔ)關(guān)鍵點(diǎn)的位置、主方向和尺度信息
              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
  % 在交互模式下顯示關(guān)鍵點(diǎn)的尺度和主方向信息
  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));
    % 將坐標(biāo)軸旋轉(zhuǎn)為關(guān)鍵點(diǎn)的方向，以確保旋轉(zhuǎn)不變性
    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);
        % 在采樣位置進(jìn)行梯度插值
        [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
        % 計(jì)算x、y方向上的權(quán)重
        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);
        % 計(jì)算高斯權(quán)重
        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
    % 將特征向量的長度歸一化，則可以進(jìn)一步去除光照變化的影響
    feat＿desc＝ feat＿desc/norm(feat＿desc);
    feat＿desc(find(feat＿desc  0.2))＝ 0.2;
    feat＿desc＝ feat＿desc/norm(feat＿desc);
    % 存儲(chǔ)特征向量
    desc＝[desc; feat＿desc];
    if(interactive  ＝ 1)&(mod(k,25)＝＝ 0)
        fprintf(2, ＇.＇);
    end
  end
  desc＿time＝ toc;
  % 調(diào)整采樣偏差
  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
  % 在交互模式下顯示運(yùn)行過程耗時(shí)
  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
  ********************************************************************

在上述實(shí)現(xiàn)SIFT算法的MATLAB程序中，調(diào)用了其他一些函數(shù)如下：

    ********************************************************************
    function g＝ gaussian＿filter(sigma)
    % 功能：生成一維高斯濾波器
    % 輸入：
    % sigma——高斯濾波器的標(biāo)準(zhǔn)差
    % 輸出：
    % 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)
    % 功能：讀入并解析一個(gè)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)
    % 功能：根據(jù)句柄重獲圖像的大小
    %  sz ＝ 圖像尺寸
    %  frac(默認(rèn)值＝ 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)
    % 功能：在原始圖像上顯示特征點(diǎn)
    % 輸入：
    % pos——特征點(diǎn)的位置矩陣
    % scale——特征點(diǎn)的尺度矩陣
    % orient——特征點(diǎn)的主方向向量
    % 輸出：
    % 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
    ********************************************************************