SIFT特徵點演算法原始碼分析（opencv）

阿新 • • 發佈：2019-01-24

SIFT特徵點演算法原始碼分析

SIFT演算法在opencv中被實現為一個類： SIFT ，主要的操作都在這個類過載的"()"運算子中實現。下面介紹這個類，以及其中呼叫的一些關鍵的函式。

SIFT類的建構函式：初始化演算法引數

SIFT::SIFT(int nfeatures=0, int nOctaveLayers=3, double contrastThreshold=0.04, double edgeThreshold=  
    10, double sigma=1.6)

引數：

nfeatures：特徵點數目（演算法對檢測出的特徵點排名，返回最好的nfeatures個特徵點）。
nOctaveLayers：金字塔中每組的層數（演算法中會自己計算這個值，後面會介紹）。
contrastThreshold：特徵點的判別閾值。這個引數的大小怎麼把握？
edgeThreshold：過濾掉邊緣效應的閾值。
sigma：金字塔第0層影象高斯濾波系

SIFT的“()”運算子：SIFT演算法的主要流程

void SIFT::operator()(InputArray _image, InputArray _mask,
                      vector<KeyPoint>& keypoints,
                      OutputArray _descriptors,
                      bool useProvidedKeypoints) const
{

    int firstOctave = -1, actualNOctaves = 0, actualNLayers = 0;
firstOctave 代表最底層一組的組號。如果為0，則源影象所在的組就是最底層一組；如果為 -1，則需要將源影象放大 1 倍（先與sigma=0.5的高斯核卷積再雙線性插值放大），放大後的影象位於最底層一組；
actualNOctaves 實際所需金字塔共多少組；actualNLayers 每組金字塔實際需要多少層。只有當特徵點keypoints被直接輸入時，這兩個引數才用得著。如果是從從影象中找keypoint，可以忽略這兩引數。


    Mat image = _image.getMat(), mask = _mask.getMat();

    if( image.empty() || image.depth() != CV_8U )
        CV_Error( CV_StsBadArg, "image is empty or has incorrect depth (!=CV_8U)" );

    if( !mask.empty() && mask.type() != CV_8UC1 )
        CV_Error( CV_StsBadArg, "mask has incorrect type (!=CV_8UC1)" );


    if( useProvidedKeypoints )
    {
        firstOctave = 0;
        int maxOctave = INT_MIN;
        for( size_t i = 0; i < keypoints.size(); i++ )
        {
            int octave, layer;
            float scale;
            unpackOctave(keypoints[i], octave, layer, scale);
            firstOctave = std::min(firstOctave, octave);
            maxOctave = std::max(maxOctave, octave);
            actualNLayers = std::max(actualNLayers, layer-2);

        }

        firstOctave = std::min(firstOctave, 0);
        CV_Assert( firstOctave >= -1 && actualNLayers <= nOctaveLayers );
        actualNOctaves = maxOctave - firstOctave + 1;
    }
如果特徵點是直接輸入的，則依次讀取輸入的特徵點，看看這些特徵點需要多少層、多少組金字塔。
約束： 
最底層的金字塔組號至少大於 -1 ，即源影象最多被放大一倍；
實際每組金字塔的層數不得大於nOctaveLayers，nOctaveLayers是在初始化SIFT時指定的；
實際的金字塔組數根據特徵點中的組號確定（按需分配，節省空間和提高效率）。


    Mat base = createInitialImage(image, firstOctave < 0, (float)sigma);
    vector<Mat> gpyr, dogpyr;
    int nOctaves = actualNOctaves > 0 ? actualNOctaves : cvRound(log( (double)std::min( base.cols, base.rows ) ) / log(2.) - 2) - firstOctave;
對原始影象進行 方差=sigma^2 的高斯濾波，得到基礎影象base：金字塔最底層的影象。
計算金字塔的組數 nOctaves： 使得最高的一組金字塔的寬（或高）只有8個畫素（or 4個畫素？）。

    //double t, tf = getTickFrequency();
    //t = (double)getTickCount();
    buildGaussianPyramid(base, gpyr, nOctaves);
    buildDoGPyramid(gpyr, dogpyr);
建立差分高斯金字塔DOG

    //t = (double)getTickCount() - t;
    //printf("pyramid construction time: %g\n", t*1000./tf);

    if( !useProvidedKeypoints )
    {
        //t = (double)getTickCount();
        findScaleSpaceExtrema(gpyr, dogpyr, keypoints);
        KeyPointsFilter::removeDuplicated( keypoints );

        if( nfeatures > 0 )
            KeyPointsFilter::retainBest(keypoints, nfeatures);
        //t = (double)getTickCount() - t;
        //printf("keypoint detection time: %g\n", t*1000./tf);

        if( firstOctave < 0 )
            for( size_t i = 0; i < keypoints.size(); i++ )
            {
                KeyPoint& kpt = keypoints[i];
                float scale = 1.f/(float)(1 << -firstOctave);
                kpt.octave = (kpt.octave & ~255) | ((kpt.octave + firstOctave) & 255);
                kpt.pt *= scale;
                kpt.size *= scale;
            }

        if( !mask.empty() )
            KeyPointsFilter::runByPixelsMask( keypoints, mask );
    }
    else
    {
        // filter keypoints by mask
        //KeyPointsFilter::runByPixelsMask( keypoints, mask );
    }
如果特徵點不是直接輸入的（false=useProvidedKeypoints），則需要先呼叫findScaleSpaceExtrema()找到特徵點，即尺度空間內的極值點（正的極大值和負的極小值），並呼叫removeDuplicated()去掉一些重複的特徵點（座標pt、鄰域大小size、方向angle都一致的特徵點被認為是重複）。如果特徵點個數超出上限，還要用retainBest() 來挑出最好的nfeatures個特徵點(響應response越高的特徵點被認為越"好")。
另外，由於findScaleSpaceExtrema()找特徵點時預設最底層的金字塔是第0組，因此如果實際的最底層是第 -1 組的話，需要將求出的關鍵點的組號減一，並將關鍵點的座標pt、鄰域尺寸size等乘以scale(=0.5，縮小一倍)。關鍵點的座標和鄰域尺寸都是以第0組的影象為基準的。


    if( _descriptors.needed() )
    {
        //t = (double)getTickCount();
        int dsize = descriptorSize();

descriptorSize()被定義為  SIFT_DESCR_WIDTH*SIFT_DESCR_WIDTH*SIFT_DESCR_HIST_BINS , 即 4*4*8，共4*4個種子點，每個種子點需要8個bin(梯度直方圖)。

        _descriptors.create((int)keypoints.size(), dsize, CV_32F);
        Mat descriptors = _descriptors.getMat();

        calcDescriptors(gpyr, keypoints, descriptors, nOctaveLayers, firstOctave);
        //t = (double)getTickCount() - t;
        //printf("descriptor extraction time: %g\n", t*1000./tf);
    }
}

呼叫calcDescriptors()為每個特徵點計算特徵向量，特徵向量是128維的。

子流程 createInitialImage() ：建立初始影象（第0層影象）

static Mat createInitialImage( const Mat& img, bool doubleImageSize, float sigma )
{
    Mat gray, gray_fpt;
    if( img.channels() == 3 || img.channels() == 4 )
        cvtColor(img, gray, COLOR_BGR2GRAY);
    else
        img.copyTo(gray);
    gray.convertTo(gray_fpt, DataType<sift_wt>::type, SIFT_FIXPT_SCALE, 0);

    float sig_diff;

    if( doubleImageSize )
    {
        sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA * 4, 0.01f) );
        Mat dbl;
        resize(gray_fpt, dbl, Size(gray.cols*2, gray.rows*2), 0, 0, INTER_LINEAR);
        GaussianBlur(dbl, dbl, Size(), sig_diff, sig_diff);
        return dbl;
    }
    else
    {
        sig_diff = sqrtf( std::max(sigma * sigma - SIFT_INIT_SIGMA * SIFT_INIT_SIGMA, 0.01f) );
        GaussianBlur(gray_fpt, gray_fpt, Size(), sig_diff, sig_diff);
        return gray_fpt;
    }
}

如果需要將影象放大一倍(doubleImageSize為true)，則現將源影象用雙線性插值放大一倍，再進行方差 = ( sigma^2 - 4*0.5^2 ) 的高斯濾波，生成基礎影象；
如果不需要放大，則直接將源影象進行方差 = ( sigma^2 - 0.5^2 ) 的高斯濾波，生成基礎影象；（0.5^2不乘4）
0.5是哪裡來的？輸入的源影象被認為是由某個假想的原始影象經過方差 = 0.5^2 的高斯濾波後得來的。如果需要進行放大，那放大後的源影象就是由某個原始影象放大後經過方差=(2*0.5)^2=4*0.5^2 高斯濾波得來。所有的 sigma 都是以假想的原始影象為基準的，所以做高斯濾波時，方差是 sigma^2 減去 4*0.5^2 或 0.5^2，而不是直接用sigma^2 。（兩個高斯核的疊加作用仍等效於一個新的高斯核，新高斯核的方差等於前兩個核的方差之和）

子流程 findScaleSpaceExtrema() ：尋找尺度空間極值點（即特徵點）

void SIFT::findScaleSpaceExtrema( const vector<Mat>& gauss_pyr, const vector<Mat>& dog_pyr,
                                  vector<KeyPoint>& keypoints ) const
{
    int nOctaves = (int)gauss_pyr.size()/(nOctaveLayers + 3);
    int threshold = cvFloor(0.5 * contrastThreshold / nOctaveLayers * 255 * SIFT_FIXPT_SCALE);
    const int n = SIFT_ORI_HIST_BINS;
    float hist[n];
    KeyPoint kpt;

    keypoints.clear();

    for( int o = 0; o < nOctaves; o++ )
        for( int i = 1; i <= nOctaveLayers; i++ )
        {
            int idx = o*(nOctaveLayers+2)+i;
            const Mat& img = dog_pyr[idx];
            const Mat& prev = dog_pyr[idx-1];
            const Mat& next = dog_pyr[idx+1];
            int step = (int)img.step1();
            int rows = img.rows, cols = img.cols;

            for( int r = SIFT_IMG_BORDER; r < rows-SIFT_IMG_BORDER; r++)
            {
                const sift_wt* currptr = img.ptr<sift_wt>(r);
                const sift_wt* prevptr = prev.ptr<sift_wt>(r);
                const sift_wt* nextptr = next.ptr<sift_wt>(r);

                for( int c = SIFT_IMG_BORDER; c < cols-SIFT_IMG_BORDER; c++)
                {
                    sift_wt val = currptr[c];

                    // find local extrema with pixel accuracy
                    if( std::abs(val) > threshold &&
                       ((val > 0 && val >= currptr[c-1] && val >= currptr[c+1] &&
                         val >= currptr[c-step-1] && val >= currptr[c-step] && val >= currptr[c-step+1] &&
                         val >= currptr[c+step-1] && val >= currptr[c+step] && val >= currptr[c+step+1] &&
                         val >= nextptr[c] && val >= nextptr[c-1] && val >= nextptr[c+1] &&
                         val >= nextptr[c-step-1] && val >= nextptr[c-step] && val >= nextptr[c-step+1] &&
                         val >= nextptr[c+step-1] && val >= nextptr[c+step] && val >= nextptr[c+step+1] &&
                         val >= prevptr[c] && val >= prevptr[c-1] && val >= prevptr[c+1] &&
                         val >= prevptr[c-step-1] && val >= prevptr[c-step] && val >= prevptr[c-step+1] &&
                         val >= prevptr[c+step-1] && val >= prevptr[c+step] && val >= prevptr[c+step+1]) ||
                        (val < 0 && val <= currptr[c-1] && val <= currptr[c+1] &&
                         val <= currptr[c-step-1] && val <= currptr[c-step] && val <= currptr[c-step+1] &&
                         val <= currptr[c+step-1] && val <= currptr[c+step] && val <= currptr[c+step+1] &&
                         val <= nextptr[c] && val <= nextptr[c-1] && val <= nextptr[c+1] &&
                         val <= nextptr[c-step-1] && val <= nextptr[c-step] && val <= nextptr[c-step+1] &&
                         val <= nextptr[c+step-1] && val <= nextptr[c+step] && val <= nextptr[c+step+1] &&
                         val <= prevptr[c] && val <= prevptr[c-1] && val <= prevptr[c+1] &&
                         val <= prevptr[c-step-1] && val <= prevptr[c-step] && val <= prevptr[c-step+1] &&
                         val <= prevptr[c+step-1] && val <= prevptr[c+step] && val <= prevptr[c+step+1])))

這個這麼長的if語句尋找的是正的極大值和負的極小值

                    {
                        int r1 = r, c1 = c, layer = i;
                        if( !adjustLocalExtrema(dog_pyr, kpt, o, layer, r1, c1,
                                                nOctaveLayers, (float)contrastThreshold,
                                                (float)edgeThreshold, (float)sigma) )
                            continue;

adjustLocalExtrema() 用來精確求解角點位置，並過濾掉邊緣點。後面會解釋這個函式。



                        float scl_octv = kpt.size*0.5f/(1 << o);
                        float omax = calcOrientationHist(gauss_pyr[o*(nOctaveLayers+3) + layer],
                                                         Point(c1, r1),
                                                         cvRound(SIFT_ORI_RADIUS * scl_octv),
                                                         SIFT_ORI_SIG_FCTR * scl_octv,
                                                         hist, n);

求出關鍵點鄰域內的梯度方向直方圖並求出最大的bin包含的元素數量 omax。此直方圖在與關鍵點尺度相對應的金字塔層中求取。
scl_octv是關鍵點的組內尺度座標，SIFT_ORI_RADIUS被定義為3，因此鄰域的半徑為 3*關鍵點組內尺度；
求梯度直方圖時對鄰域內的點會依其距關鍵點的距離進行高斯加權，加權的Sigma=1.5*關鍵點組內尺度。（SIFT_ORI_SIG_FCTR被定義為1.5）。
到現在為止只是求出了直方圖，但還並未提取出關鍵點的主方向和輔方向。提取關鍵點方向由下面的幾行實現。

下面的SIFT_ORI_PEAK_RATIO被定義為0.8。元素數量多於 0.8 倍omax 的bin都被認為是關鍵點的主方向或輔方向。每個主/輔方向都用一個單獨的特徵點表示。

                        float mag_thr = (float)(omax * SIFT_ORI_PEAK_RATIO);
                        for( int j = 0; j < n; j++ )
                        {
                            int l = j > 0 ? j - 1 : n - 1;
                            int r2 = j < n-1 ? j + 1 : 0;

                            if( hist[j] > hist[l]  &&  hist[j] > hist[r2]  &&  hist[j] >= mag_thr )
                            {
                                float bin = j + 0.5f * (hist[l]-hist[r2]) / (hist[l] - 2*hist[j] + hist[r2]);
j是整數，但直方圖的極值點一般並不是準確地位於整數位置，因此這裡進行了極值點擬合。對直方圖一維二階泰勒展開，再取導數為0，即可擬合出精度更高的極值點，極值點的位置等於 (負的一階導數/二階導數)。

                                bin = bin < 0 ? n + bin : bin >= n ? bin - n : bin;
                                kpt.angle = 360.f - (float)((360.f/n) * bin);
                                if(std::abs(kpt.angle - 360.f) < FLT_EPSILON)
                                    kpt.angle = 0.f;
                                keypoints.push_back(kpt);
                            }
                        }

                    }
                }
            }
        }
}

上面出現了一個新函式：adjustLocalExtrema() 。它可以用來精確求解角點位置（二階泰勒級數擬合），並過濾掉邊緣點。下面闡述該函式的原理。

在尺度空間（三維）的 (r1,c1,layer) 點處，進行三維二階泰勒展開，可以求取極值的精確位置。
三維二階泰勒展開的原理如下：
△f = f ( x+△x, y+△y, z+△z ) - f ( x, y, z )
= f ( x+△x, y+△y, z+△z ) - f ( x+△x, y+△y, z ) +
f ( x+△x, y+△y, z ) - f ( x+△x, y, z ) +
f ( x+△x, y, z ) - f ( x, y, z )

// 進行泰勒展開近似
≈ Fz (x+△x, y+△y, z) * △z + 0.5 * Fzz(x, y, z) * △z^2 +
Fy (x+△x, y, z) * △y + 0.5 * Fyy(x, y, z) * △y^2 +
Fx (x, y, z) * △x + 0.5 * Fxx(x, y, z) * △x^2

≈ [ Fz (x, y, z) + Fzx * △x + Fyz * △y ] * △z + 0.5 * Fzz(x, y, z) * △z^2 +
[ Fy (x, y, z) + Fxy * △x ] * △y + 0.5 * Fyy(x, y, z) * △y^2 +
Fx (x, y, z) * △x + 0.5 * Fxx(x, y, z) * △x^2

上面推導的基本思路是， f ( x, y, z ) 到 f ( x+△x, y+△y, z+△z ) 之間的增量，可以分解為：

" f ( x, y, z ) 到 f ( x+△x, y, z ) 的增量 " +

" f ( x+△x, y, z ) 到 f ( x+△x, y+△y, z ) 的增量 " +

" f ( x+△x, y+△y, z ) 到 f ( x+△x, y+△y, z+△z ) 的增量 "。

而這三個分量分別都可以用單變數（一維）的二階泰勒展開求解，由此組合起來就可以得到三維的二階泰勒展開式。

整理一下就是：
△f ≈ G0' * △ + 0.5* △'* H0 * △
其中G0代表點(x,y,z)處的梯度向量，單引號上角標代表轉置，△代表偏移( △x, △y, △z )，H0 代表點(x,y,z)處的Hessian矩陣。對兩邊求導並取導數為0的點，[△(△f)]/[△] = 0，得到的解就是極值點。極值點的方程為： -G0 = H0 * △，或 G0 = H0 * (-△)
這個公式以及上面的推理方法同樣適用於更高維的泰勒展開。

另外，上式的第二項 △'* H0 * △，構成了一個二次型：
2*(△f - G0' * △) = △'* H0 * △
H 的特徵向量反映了二階分量的變化方向，特徵值反映了對應的方向上二階分量變化的快慢和正負。在邊緣區域，垂直邊緣的方向上，一、二階分量變化都很劇烈，而在平行於邊緣的方向上，一、二階分量變化都很小。因此通過比較Hessian矩陣的兩個特徵值還可辨別出位於邊緣的點：大特徵值與小特徵值之比遠大於1的點是邊緣點（這些邊緣點要被幹掉）。設兩個特徵值分別是α和β，則 Tr(H)^2/det(H) = (α+β)^2/αβ = (1+γ)^2/γ 的值就可以反映出大、小特徵值之比。利用這個關係可以避免特徵值的求取。
這跟 Harris 角點檢測中去除邊緣點的原理類似，只不過Harris 角點檢測中用到的2X2特徵矩陣是由一階導數的二次式構成，而這裡的 Hession 矩陣是由二階導數構成。

adjustLocalExtrema() 的實現程式碼如下：

static bool adjustLocalExtrema( const vector<Mat>& dog_pyr, KeyPoint& kpt, int octv,
                                int& layer, int& r, int& c, int nOctaveLayers,
                                float contrastThreshold, float edgeThreshold, float sigma )
{
    const float img_scale = 1.f/(255*SIFT_FIXPT_SCALE);
    const float deriv_scale = img_scale*0.5f;
    const float second_deriv_scale = img_scale;
    const float cross_deriv_scale = img_scale*0.25f;

    float xi=0, xr=0, xc=0, contr=0;
    int i = 0;

    for( ; i < SIFT_MAX_INTERP_STEPS; i++ )
    {
        int idx = octv*(nOctaveLayers+2) + layer;
        const Mat& img = dog_pyr[idx];
        const Mat& prev = dog_pyr[idx-1];
        const Mat& next = dog_pyr[idx+1];

        Vec3f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale,
                 (img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale,
                 (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale);

        float v2 = (float)img.at<sift_wt>(r, c)*2;
        float dxx = (img.at<sift_wt>(r, c+1) + img.at<sift_wt>(r, c-1) - v2)*second_deriv_scale;
        float dyy = (img.at<sift_wt>(r+1, c) + img.at<sift_wt>(r-1, c) - v2)*second_deriv_scale;
        float dss = (next.at<sift_wt>(r, c) + prev.at<sift_wt>(r, c) - v2)*second_deriv_scale;
        float dxy = (img.at<sift_wt>(r+1, c+1) - img.at<sift_wt>(r+1, c-1) -
                     img.at<sift_wt>(r-1, c+1) + img.at<sift_wt>(r-1, c-1))*cross_deriv_scale;
        float dxs = (next.at<sift_wt>(r, c+1) - next.at<sift_wt>(r, c-1) -
                     prev.at<sift_wt>(r, c+1) + prev.at<sift_wt>(r, c-1))*cross_deriv_scale;
        float dys = (next.at<sift_wt>(r+1, c) - next.at<sift_wt>(r-1, c) -
                     prev.at<sift_wt>(r+1, c) + prev.at<sift_wt>(r-1, c))*cross_deriv_scale;

        Matx33f H(dxx, dxy, dxs,
                  dxy, dyy, dys,
                  dxs, dys, dss);

        Vec3f X = H.solve(dD, DECOMP_LU);

H.solve(dD)用於求解線性方程組，該方程組是  dD = H * X 。
前面提到極值點的方程是 ： G0 = H0 * (-△) ，這裡G0就是dD,H0就是H，(-△)就是X。因此為求得△，只需對方程解X求反。如下面的xi、xr和xc。

        xi = -X[2];
        xr = -X[1];
        xc = -X[0];

        if( std::abs(xi) < 0.5f && std::abs(xr) < 0.5f && std::abs(xc) < 0.5f )
            break;

        if( std::abs(xi) > (float)(INT_MAX/3) ||
            std::abs(xr) > (float)(INT_MAX/3) ||
            std::abs(xc) > (float)(INT_MAX/3) )
            return false;

        c += cvRound(xc);
        r += cvRound(xr);
        layer += cvRound(xi);


如果擬合的極值點與當前的整數點偏差大於0.5，則說明真正的極值點偏離當前整數點很多，需要換一個更好的整數點，然後重新擬合。知道偏差小於0.5或迭代超過一定次數。



        if( layer < 1 || layer > nOctaveLayers ||
            c < SIFT_IMG_BORDER || c >= img.cols - SIFT_IMG_BORDER  ||
            r < SIFT_IMG_BORDER || r >= img.rows - SIFT_IMG_BORDER )
            return false;
    }

    // ensure convergence of interpolation
    if( i >= SIFT_MAX_INTERP_STEPS )
        return false;

    {
        int idx = octv*(nOctaveLayers+2) + layer;
        const Mat& img = dog_pyr[idx];
        const Mat& prev = dog_pyr[idx-1];
        const Mat& next = dog_pyr[idx+1];
        Matx31f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale,
                   (img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale,
                   (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale);
        float t = dD.dot(Matx31f(xc, xr, xi));

        contr = img.at<sift_wt>(r, c)*img_scale + t * 0.5f;
        if( std::abs( contr ) * nOctaveLayers < contrastThreshold )
            return false;

        // principal curvatures are computed using the trace and det of Hessian
        float v2 = img.at<sift_wt>(r, c)*2.f;
        float dxx = (img.at<sift_wt>(r, c+1) + img.at<sift_wt>(r, c-1) - v2)*second_deriv_scale;
        float dyy = (img.at<sift_wt>(r+1, c) + img.at<sift_wt>(r-1, c) - v2)*second_deriv_scale;
        float dxy = (img.at<sift_wt>(r+1, c+1) - img.at<sift_wt>(r+1, c-1) -
                     img.at<sift_wt>(r-1, c+1) + img.at<sift_wt>(r-1, c-1)) * cross_deriv_scale;
        float tr = dxx + dyy;
        float det = dxx * dyy - dxy * dxy;

求出了Hessian矩陣的軌跡和行列式，以便做邊緣判斷。

        if( det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det )
            return false;
    }


    kpt.pt.x = (c + xc) * (1 << octv);
    kpt.pt.y = (r + xr) * (1 << octv);
    kpt.octave = octv + (layer << 8) + (cvRound((xi + 0.5)*255) << 16);
    kpt.size = sigma*powf(2.f, (layer + xi) / nOctaveLayers)*(1 << octv)*2;
    kpt.response = std::abs(contr);

    return true;
}

注意，關鍵點kpt的座標pt和鄰域size都是以第0組影象為基準的。size代表鄰域的直徑，它的大小為 2*Sigma, Sigma是關鍵點所在的尺度。鄰域的半徑正好是 Sigma 。

子流程 calcDescriptors() ：特徵點描述符的計算

static void calcDescriptors(const vector<Mat>& gpyr, const vector<KeyPoint>& keypoints,
                            Mat& descriptors, int nOctaveLayers, int firstOctave )
{
    int d = SIFT_DESCR_WIDTH, n = SIFT_DESCR_HIST_BINS;

SIFT_DESCR_WIDTH 是4，SIFT_DESCR_HIST_BINS是8。


    for( size_t i = 0; i < keypoints.size(); i++ )
    {
        KeyPoint kpt = keypoints[i];
        int octave, layer;
        float scale;
        unpackOctave(kpt, octave, layer, scale);

        CV_Assert(octave >= firstOctave && layer <= nOctaveLayers+2);
        float size=kpt.size*scale;
        Point2f ptf(kpt.pt.x*scale, kpt.pt.y*scale);
        const Mat& img = gpyr[(octave - firstOctave)*(nOctaveLayers + 3) + layer];

解析出關鍵點金字塔組號、層號和比例。再根據組號和層號找到對應層的金字塔影象，根據scale計算出在該層金字塔影象上關鍵點的位置和鄰域大小。（關鍵點kpt內部儲存的座標pt和size都是以第0組第0層影象為基準的，因此要想在其他層中使用，需要先乘以scale變換一下）

        float angle = 360.f - kpt.angle;
        if(std::abs(angle - 360.f) < FLT_EPSILON)
            angle = 0.f;
        calcSIFTDescriptor(img, ptf, angle, size*0.5f, d, n, descriptors.ptr<float>((int)i));
計算關鍵點的描述符（一個128維的向量）

    }
}






static void calcSIFTDescriptor( const Mat& img, Point2f ptf, float ori, float scl,
                               int d, int n, float* dst )
{
    Point pt(cvRound(ptf.x), cvRound(ptf.y));
    float cos_t = cosf(ori*(float)(CV_PI/180));
    float sin_t = sinf(ori*(float)(CV_PI/180));
    float bins_per_rad = n / 360.f;
    float exp_scale = -1.f/(d * d * 0.5f);

exp_scale = -1.f / [2*(d/2)^2], 所以高斯加權的Sigma=d/2 ;
d=4,n=8。bins_per_rad表示每1°包含多少個bin。

    float hist_width = SIFT_DESCR_SCL_FCTR * scl;
    int radius = cvRound(hist_width * 1.4142135623730951f * (d + 1) * 0.5f);

scl等於 size*0.5,正好對應關鍵點的尺度。
SIFT_DESCR_SCL_FCTR = 3 ，hist_width是每個種子點（一個矩形區域）的寬度，共4*4個種子點。
radius 理論上應為 hist_width*d/2，但考慮到雙線性差值的需要，應將radius額外擴充套件一些，所以乘以 d+1；又考慮到旋轉時正方形的對角線，所以再乘以 根號2。


    // Clip the radius to the diagonal of the image to avoid autobuffer too large exception
    radius = std::min(radius, (int) sqrt((double) img.cols*img.cols + img.rows*img.rows));
    cos_t /= hist_width;
    sin_t /= hist_width;
這時 cos_t 和 sin_t 已經不單單是正餘弦值了，他們除了一個係數。

    int i, j, k, len = (radius*2+1)*(radius*2+1), histlen = (d+2)*(d+2)*(n+2);
    int rows = img.rows, cols = img.cols;

    AutoBuffer<float> buf(len*6 + histlen);
    float *X = buf, *Y = X + len, *Mag = Y, *Ori = Mag + len, *W = Ori + len;
    float *RBin = W + len, *CBin = RBin + len, *hist = CBin + len;

    for( i = 0; i < d+2; i++ )
    {
        for( j = 0; j < d+2; j++ )
            for( k = 0; k < n+2; k++ )
                hist[(i*(d+2) + j)*(n+2) + k] = 0.;
    }

    for( i = -radius, k = 0; i <= radius; i++ )
        for( j = -radius; j <= radius; j++ )
        {
            // Calculate sample's histogram array coords rotated relative to ori.
            // Subtract 0.5 so samples that fall e.g. in the center of row 1 (i.e.
            // r_rot = 1.5) have full weight placed in row 1 after interpolation.
            float c_rot = j * cos_t - i * sin_t;
            float r_rot = j * sin_t + i * cos_t;

由於cos_t 和 sin_t除過了hist_width，所以這時旋轉後的座標的 c_rot 和 r_rot 的單位直接就是種子點了。下面的 rbin 和 cbin 計算的就是種子點的編號。種子點的編號從(-1,-1)到(d,d)。減去0.5是為了方便差值時計算權重。

            float rbin = r_rot + d/2 - 0.5f;
            float cbin = c_rot + d/2 - 0.5f;
            int r = pt.y + i, c = pt.x + j;

            if( rbin > -1 && rbin < d && cbin > -1 && cbin < d &&
                r > 0 && r < rows - 1 && c > 0 && c < cols - 1 )
超出窗邊界或影象邊界的點不予考慮。注意這裡 rbin 和 cbin 的下線都是 -1，也就是說它們可能為負
            {
                float dx = (float)(img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1));
                float dy = (float)(img.at<sift_wt>(r-1, c) - img.at<sift_wt>(r+1, c));
                X[k] = dx; Y[k] = dy; RBin[k] = rbin; CBin[k] = cbin;
                W[k] = (c_rot * c_rot + r_rot * r_rot)*exp_scale;
                k++;
            }
        }

    len = k;
    fastAtan2(Y, X, Ori, len, true);
    magnitude(X, Y, Mag, len);
    exp(W, W, len);
Ori 存放角度陣列，Mag存放梯度的能量陣列，W存放高斯權重

    for( k = 0; k < len; k++ )
    {
        float rbin = RBin[k], cbin = CBin[k];
        float obin = (Ori[k] - ori)*bins_per_rad;
        float mag = Mag[k]*W[k];
mag等於高斯權重乘以梯度能量

        int r0 = cvFloor( rbin );
        int c0 = cvFloor( cbin );
        int o0 = cvFloor( obin );
由於 rbin 可能為負，所以 r0 可能為 -1，這樣r0的取值就是 [-1,d-1]，而後面雙線性插值時還會擴充套件到d,因此r0的取值是[-1,d]，共d+2個取值，c0也一樣，所以hist陣列的r維和c維的size都是 d+2。


        rbin -= r0;
        cbin -= c0;
        obin -= o0;

只保留rbin、cbin和obin的小數部分。每個樣點對與之相鄰的各個bin的貢獻比重與該樣點到每個bin的距離成反比。


        if( o0 < 0 )
            o0 += n;
        if( o0 >= n )
            o0 -= n;
角度的bin是迴圈的，所以處理一下 o0。這樣o0的範圍就是 [0,n-1] ，但後面雙線性差值時可能還會擴充套件到n，因此o0的範圍就是 [0,n],共n+1個值，所以 hist 陣列o維的長度應為 n+1 。( ?? 這裡的程式中o維的長度定義為了 n+2，是我分析錯了還是opencv確實保留得多了? )。

        // histogram update using tri-linear interpolation
        float v_r1 = mag*rbin, v_r0 = mag - v_r1;	
        float v_rc11 = v_r1*cbin, v_rc10 = v_r1 - v_rc11;
        float v_rc01 = v_r0*cbin, v_rc00 = v_r0 - v_rc01;
        float v_rco111 = v_rc11*obin, v_rco110 = v_rc11 - v_rco111;
        float v_rco101 = v_rc10*obin, v_rco100 = v_rc10 - v_rco101;
        float v_rco011 = v_rc01*obin, v_rco010 = v_rc01 - v_rco011;
        float v_rco001 = v_rc00*obin, v_rco000 = v_rc00 - v_rco001;

對於一個三維的直方圖， 樣點對與之相鄰的位於(r,c,o)處的bin的貢獻是：
mag * (dis(r,rbin)) * (dis(c,cbin)) * (dis(o,obin))
其中 dis(r,rbin) 表示r到rbin的距離（差）。注意r,c,o需要滿足：dis(r,rbin)<1 && dis(c,cbin)<1 && dis(o,obin)<1 ,才算是與樣點相鄰。


        int idx = ((r0+1)*(d+2) + c0+1)*(n+2) + o0;
        hist[idx] += v_rco000;
        hist[idx+1] += v_rco001;
        hist[idx+(n+2)] += v_rco010;
        hist[idx+(n+3)] += v_rco011;
        hist[idx+(d+2)*(n+2)] += v_rco100;
        hist[idx+(d+2)*(n+2)+1] += v_rco101;
        hist[idx+(d+3)*(n+2)] += v_rco110;
        hist[idx+(d+3)*(n+2)+1] += v_rco111;
    }

    // finalize histogram, since the orientation histograms are circular
    for( i = 0; i < d; i++ )
        for( j = 0; j < d; j++ )
        {
            int idx = ((i+1)*(d+2) + (j+1))*(n+2);
            hist[idx] += hist[idx+n];
            hist[idx+1] += hist[idx+n+1];
            for( k = 0; k < n; k++ )
                dst[(i*d + j)*n + k] = hist[idx+k];
        }
將 hist 拷貝到 dst
    // copy histogram to the descriptor,
    // apply hysteresis thresholding
    // and scale the result, so that it can be easily converted
    // to byte array
    float nrm2 = 0;
    len = d*d*n;
    for( k = 0; k < len; k++ )
        nrm2 += dst[k]*dst[k];
    float thr = std::sqrt(nrm2)*SIFT_DESCR_MAG_THR;
    for( i = 0, nrm2 = 0; i < k; i++ )
    {
        float val = std::min(dst[i], thr);
        dst[i] = val;
        nrm2 += val*val;
    }
    nrm2 = SIFT_INT_DESCR_FCTR/std::max(std::sqrt(nrm2), FLT_EPSILON);
抑制暴走的特徵：每個特徵不能大於特徵向量總模長的 0.2 倍。

#if 1
    for( k = 0; k < len; k++ )
    {
        dst[k] = saturate_cast<uchar>(dst[k]*nrm2);
    }
歸一化並將每個特徵都轉化為8位深度。
#else
    float nrm1 = 0;
    for( k = 0; k < len; k++ )
    {
        dst[k] *= nrm2;
        nrm1 += dst[k];
    }
    nrm1 = 1.f/std::max(nrm1, FLT_EPSILON);
    for( k = 0; k < len; k++ )
    {
        dst[k] = std::sqrt(dst[k] * nrm1);//saturate_cast<uchar>(std::sqrt(dst[k] * nrm1)*SIFT_INT_DESCR_FCTR);
    }
#endif
}

SIFT特徵點演算法原始碼分析（opencv）

SIFT類的建構函式：初始化演算法引數

SIFT的“()”運算子：SIFT演算法的主要流程

子流程 createInitialImage() ：建立初始影象（第0層影象）

子流程 findScaleSpaceExtrema() ：尋找尺度空間極值點（即特徵點）

子流程 calcDescriptors() ：特徵點描述符的計算

SIFT特徵點演算法原始碼分析（opencv）

Mahout協同過濾演算法原始碼分析（3）--parallelALS

支援向量機—SMO演算法原始碼分析（1）

SIFT特徵點檢測學習一（轉載）

OpenCV學習筆記（31）KAZE 演算法原理與原始碼分析（五）KAZE的原始碼優化及與SIFT的比較

OpenCV學習筆記（30）KAZE 演算法原理與原始碼分析（四）KAZE特徵的效能分析與比較

OpenCV角點檢測原始碼分析（Harris和ShiTomasi角點）

cornerSubPixel亞畫素角點檢測原始碼分析（最小二乘法）

Weka演算法Classifier-tree-J48原始碼分析（一）演算法和基本資料結構

第二人生原始碼分析（3）程式入口點

TensorFlow學習筆記之原始碼分析（1）----最近演算法nearest_neighbor

多執行緒高併發程式設計(12) -- 阻塞演算法實現ArrayBlockingQueue原始碼分析（1）

Android ADB 原始碼分析（三）

Mybatis 原始碼分析（2）—— 引數處理

Mybatis 原始碼分析（9）—— 事物管理

Mybatis 原始碼分析（8）—— 一二級快取

Mybatis原始碼分析（7）—— 結果集處理

Mybatis原始碼分析（6）—— 從JDBC看Mybatis的設計

Mybatis原始碼分析（5）—— 外掛的原理

Mybatis原始碼分析（4）—— Mapper的建立和獲取

SIFT特徵點演算法原始碼分析（opencv）

SIFT類的建構函式：初始化演算法引數

SIFT的“()”運算子：SIFT演算法的主要流程

子流程 createInitialImage() ： 建立初始影象（第0層影象）

子流程 findScaleSpaceExtrema() ：尋找尺度空間極值點（即特徵點）

子流程 calcDescriptors() ： 特徵點描述符的計算

相關推薦

子流程 createInitialImage() ：建立初始影象（第0層影象）

子流程 calcDescriptors() ：特徵點描述符的計算