2. 程式人生 > >SIFT演算法特徵描述子構建---關鍵點定位原理及程式碼

SIFT演算法特徵描述子構建---關鍵點定位原理及程式碼

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

0.引言

sift針對區域性特徵進行特徵提取，在尺度空間尋找極值點，提取位置，尺度，旋轉不變數，生成特徵描述子。

總共分四個步驟：

step2 關鍵點/極值點提取

2.1 關鍵點位置初步探查

生成DOG金字塔後，要找到DOG空間中的區域性極值點。第一步首先根據DOG空間內，各點響應大小檢測和篩選極值點。
選擇策略如下：
這裡寫圖片描述
中間的檢測點和它同尺度的8個相鄰點和上下相鄰尺度對應的9×2個點共26個點比較，以確保在尺度空間和二維影象空間都檢測到極值點。
由於要在相鄰尺度進行比較，高斯差分金子塔S+2層，頭尾兩層無法計算，因此，定義的INTERVALS層數就是指檢測極值點的層數。

bool isExtremum(int x,int y,const vector<Mat>& dog_pyr,int index)//x為列，y為行
{
    pixel_t* data = (pixel_t*)dog_pyr[index].data;          //當前層資料的指標
    pixel_t val = *(data + dog_pyr[index].cols*y + x);      //當前點畫素
    //檢查是否為極大值
    if (val>0)
    {
        for (int i = -1; i < 2; i++)
        {
            int
 
 col_cur= dog_pyr[index+i].cols;             //當前層列數
            for (int j = -1; j < 2; j++)
            {
                for (int k = -1; k < 2; k++)
                {
                    if (val<*((pixel_t*)dog_pyr[index+i].data+col_cur*(y+j)+x+k))
                    {
                        return
 
 false;
                    }
                }
            }
        }
    }//極大值
     //檢查是否為極小值
    else
    {
        for (int i = -1; i < 2; i++)
        {
            int col_cur = dog_pyr[index + i].cols;              //當前層列數
            for (int j = -1; j < 2; j++)
            {
                for (int k = -1; k < 2; k++)
                {
                    if (val > *((pixel_t*)dog_pyr[index + i].data + col_cur*(y + j) + x + k))
                    {
                        return false;
                    }
                }
            }
        }
    }//極小值
    return true;
}

當然這樣產生離散空間的極值點並不全都是穩定的特徵點，因為某些極值點響應較弱（對比度低），同時DOG運算元會產生較強的邊緣響應。為了解決這兩個問題，需要對極值點進行插值和篩選。

2.2 關鍵點插值定位修正和邊緣響應去除

2.2.1位置修正及去除低響應

插值需要將離散的影象表示式，表達為連續的函式。以下通過擬合三維二次函式（二階的泰勒展開式）來精確確定關鍵點的位置和尺度。
這裡寫圖片描述

其中這裡寫圖片描述是位置和尺度三個變數的向量。

二次函式對稱軸就是極值點偏移量：這裡寫圖片描述

此處，為了求此公式，需要提前定義好計算逆矩陣，二階偏導，一階偏導 的函式。

/*************************************************************************************************************************
*模組說明：
*       求偏差量預備工作 【1】返回第index層的金字塔在y行x列的值
*功能說明：
*       PryAt（）---At：在該層y行x列的值
**************************************************************************************************************************/
pixel_t PyrAt(const vector<Mat>& pyr, int index, int x, int y)
{
    pixel_t *data = (pixel_t*)pyr[index].data;
    pixel_t val = *(data + y*pyr[index].cols + x);
    return val;
}
#define At(index, x, y) (PyrAt(dog_pyr, (index), (x), (y)))  
/*************************************************************************************************************************
*模組說明：
*       求偏差量預備工作，【2】對x,y,sigma的一階偏導存入dx[3]
*功能說明：
*       使用節點的一階中心差分近似導數，  取最小h=1時，最精確
*       f'(xi)=(f(xi+h)-f(xi-h))/2h
*       f'(yi)=(f(yi+h)-f(yi-h))/2h
*       f'(si)=(f(si+h)-f(si-h))/2h
**************************************************************************************************************************/
void DerivativeOf3D(int x,int y,const vector<Mat>& dog_pyr,int index,double *dx)
{
    dx[0] = (double)(At(index, x + 1, y) - At(index, x - 1, y)) / 2.0;
    dx[1] = (double)(At(index, x, y + 1) - At(index, x, y - 1)) / 2.0;
    dx[2] = (double)(At(index + 1, x, y) - At(index - 1, x, y)) / 2.0;
}
/*************************************************************************************************************************
*模組說明：
*       求偏差量預備工作，【3】對x,y,sigma的二階偏導存入Hat(3*3),對稱陣
*   |   Dxx Dxy Dxs     |
*   |   Dyx Dxy Dys     |
*   |   Dsx Dsy Dss     |
*功能說明：
*       使用節點的一階中心差分近似導數，  取最小h=1時，最精確
*       Dxx=f''(xi)=(f(xi+h)+f(xi-h)-2*f(xi))/h^2
*       Dxy=f''(yi)=(f(yi+h)+f(yi-h)-2*f(yi))/h^2
*       Dss=f''(si)=(f(si+h)+f(si-h)-2*f(si)))/h^2
*
*       Dyx=Dxy=f''(xi)(yi)=(f(xi+h,yi+h)+f(xi-h,yi-h)-f(xi+h,yi-h)-f(xi-h,yi+h))/4*h
**************************************************************************************************************************/
#define Hat(i,j) (*(H+3*(i)+(j)))
void Hessian3D(int x, int y, const vector<Mat>& dog_pyr, int index, double*H)
{
    Hat(0, 0) = At(index, x + 1, y) + At(index, x - 1, y) - 2 * At(index, x, y);//Dxx
    Hat(1, 1) = At(index, x, y + 1) + At(index, x, y - 1) - 2 * At(index, x, y);//Dyy
    Hat(1, 1) = At(index + 1, x, y) + At(index - 1, x, y) - 2 * At(index, x, y);//Dss
    Hat(1, 0) = Hat(0, 1) = (At(index, x + 1, y + 1) + At(index, x - 1, y - 1) - At(index, x + 1, y - 1) - At(index, x - 1, y + 1)) / 4.0;//Dxy
    Hat(2, 0) = Hat(0, 2) = (At(index + 1, x + 1, y) + At(index - 1, x - 1, y) - At(index - 1, x + 1, y) - At(index + 1, x - 1, y)) / 4.0;//Dxs
    Hat(2, 1) = Hat(1, 2) = (At(index + 1, x, y + 1) + At(index - 1, x, y - 1) - At(index + 1, x, y - 1) - At(index - 1, x, y + 1)) / 4.0;//Dsy
}
/*************************************************************************************************************************
*模組說明：
*       三階矩陣求逆
**************************************************************************************************************************/
#define HIat(i, j) (*(H_inve+3*(i)+j))
bool Inverse3D(const double *H,double *H_inve)
{
    //行列式
    double A =Hat(0, 0)*Hat(1, 1)*Hat(2, 2)
            + Hat(0, 1)*Hat(1, 2)*Hat(2, 0)
            + Hat(0, 2)*Hat(1, 0)*Hat(2, 1)
            - Hat(0, 0)*Hat(1, 2)*Hat(2, 1)
            - Hat(0, 1)*Hat(1, 0)*Hat(2, 2)
            - Hat(0, 2)*Hat(1, 1)*Hat(2, 0);
    if (fabs(A)<1e-10)  return false;
    HIat(0, 0) = Hat(1, 1) * Hat(2, 2) - Hat(2, 1)*Hat(1, 2);
    HIat(0, 1) = -(Hat(0, 1) * Hat(2, 2) - Hat(2, 1) * Hat(0, 2));
    HIat(0, 2) = Hat(0, 1) * Hat(1, 2) - Hat(0, 2)*Hat(1, 1);
    //伴隨矩陣第二行
    HIat(1, 0) = Hat(1, 2) * Hat(2, 0) - Hat(2, 2)*Hat(1, 0);
    HIat(1, 1) = -(Hat(0, 2) * Hat(2, 0) - Hat(0, 0) * Hat(2, 2));
    HIat(1, 2) = Hat(0, 2) * Hat(1, 0) - Hat(0, 0)*Hat(1, 2);
    //伴隨矩陣第三行
    HIat(2, 0) = Hat(1, 0) * Hat(2, 1) - Hat(1, 1)*Hat(2, 0);
    HIat(2, 1) = -(Hat(0, 0) * Hat(2, 1) - Hat(0, 1) * Hat(2, 0));
    HIat(2, 2) = Hat(0, 0) * Hat(1, 1) - Hat(0, 1)*Hat(1, 0);
    for (int i = 0; i < 9; i++)
    {
        *(H_inve + i) /= A;
    }
    return true;
}

這樣就可以計算偏差量了。

/*************************************************************************************************************************
*模組說明：計算x,y,sigma 的偏差量x^、y^、sigma^
*       
*模組內容       1.x^=-(Dxx*dx+Dxy*dy+Dxs*ds)  DOG函式擬合函式的逆矩陣的二階和擬合函式一階（為什麼這個公式）
**************************************************************************************************************************/
void GetoffsetX(int x, int y, const vector<Mat>& dog_pyr, int index, double* offset_x)
{
    //求一階偏導
    double dx[3];
    DerivativeOf3D(x, y, dog_pyr, index, dx);
    //求二階偏導
    double H[9], H_inve[9] = {0};
    Hessian3D(x, y, dog_pyr, index, H);
    if (Inverse3D(H, H_inve))
    {
        for (int i = 0; i < 3; i++)
        {
            offset_x[i] = 0.0;
            for (int j = 0; j < 3; j++)
            {
                offset_x[i] += H_inve[i * 3 + j] * dx[j];
            }
            offset_x[i]=-offset_x[i];
        }
    }
    else
        cout << "無海森逆矩陣" << endl;

得到偏移量後，首先判斷x,y,sigma任意維度的偏移量是否超過0.5，其次論文中擬合的迭代次數為5次。另外，對比度小於0.03的 極值點也要去除。
極值點處的值計算公式為：
這裡寫圖片描述

//2.D（x^）=D+0.5*D'^T(x)*x^=D()+0.5*(Dx*offset[0] +Dy* offset[1]+Ds* offset[2])
//計算D（x^）
double GetFabsDx(int x,int y,const vector<Mat>& dog_pyr,int index,const double *offset_x)
{
    double dx[3];
    DerivativeOf3D(x, y, dog_pyr, index, dx);
    double term = 0.0;
    for (int i = 0; i < 3; i++)
    {
        term += dx[i] * offset_x[i];
    }
    pixel_t val = At(index, x, y);
    return fabs(val+0.5*term);
}

整個篩選並儲存關鍵點如下：

/*************************************************************************************************************************
*模組說明：
*       模組四的第二步:修正極值點，刪除不穩定的點
*功能說明:
*       1--根據高斯差分函式產生的極值點並不全都是穩定的特徵點，因為某些極值點的響應較弱，而且DOG運算元會產生較強的邊緣響應
*       2--以上方法檢測到的極值點是離散空間的極值點，下面通過擬合三維二次函式來精確定位關鍵點的位置和尺度，同時去除對比度
*          低和不穩定的邊緣響應點(因為DOG運算元會產生較強的邊緣響應)，以增強匹配的穩定性、提高抗噪聲的能力。
*       3--修正極值點，刪除不穩定點，|D(x)| < 0.03 Lowe 2004
**************************************************************************************************************************/
Keypoint* InterpolationExtremum(int x ,int y, const vector<Mat>& dog_pyr, int index, int octave, int interval, double dxthreshold)
{
    double offset_x[3] = {0};
    const Mat &mat = dog_pyr[index];
    int idx = index;
    //int intvl=interval;
    int i = 0;
    for (; i < MAX_INTERPOLATION_STEPS; i++)
    {
        //求取偏移量
        GetoffsetX(x, y, dog_pyr, idx, offset_x); 
        //如果offset_x 的任一維度大於0.5，插值中心已經偏移到它的鄰近點上
        if (fabs(offset_x[0]) < 0.5 && fabs(offset_x[1]) < 0.5 && fabs(offset_x[2]) < 0.5)  break;  
        ////將找到的極值點對應成畫素（整數），用周圍點代替
        x += cvRound(offset_x[0]);
        y += cvRound(offset_x[1]);
        interval += cvRound(offset_x[2]);
        idx = octave* (INTERVALS+2)+ interval;
        //超出影象層或者影象邊界的範圍
        if (interval<1||interval>INTERVALS||x>mat.cols-2||x<2|| y>mat.rows - 2 || y<2)  return NULL;

    }
    //超出所設定的迭代次數
    if (i>=MAX_INTERPOLATION_STEPS) return NULL;
    //|D(x^)| < 0.03取經驗值 
    if (GetFabsDx(x, y, dog_pyr, idx, offset_x) < dxthreshold / INTERVALS)  return NULL;
    //則屬於關鍵點
    Keypoint* keypoint = new Keypoint;
    keypoint->x = x;
    keypoint->y = y;
    //keypoint->val=At()

    keypoint->offset_x = offset_x[0];
    keypoint->offset_y = offset_x[1];
    keypoint->offset_interval = offset_x[2];

    keypoint->octave = octave;
    keypoint->interval = interval;

    keypoint->dx = (x + offset_x[0])*pow(2.0, octave); //縮放成原影象大小
    keypoint->dy = (y + offset_x[1])*pow(2.0, octave);
    return keypoint;
}

2.2.2邊緣響應去除

為了剔除邊緣響應點，需要讓下比值小於一定的閾值
這裡寫圖片描述)成立時將關鍵點保留，反之剔除。
其中，在Lowe的文章中，取r＝10。在Lowe的文章中，取r＝10。
這裡寫圖片描述這裡寫圖片描述

bool passEdgeResponse(int x,int y,vector<Mat>& dog_pyr,int index,double r)
{
    //求二階偏導
    double H[9];
    Hessian3D(x, y, dog_pyr, index, H);
    double Tr_h = H[0] + H[4];//Dxx+Dyy
    double Det_h = H[0] * H[4] - H[1] * H[1];//Dxx * Dyy - Dxy * Dxy
    if (Det_h<=0)   return false;
    if (Tr_h*Tr_h / Det_h < (r + 1)*(r + 1) / r)    return true;
    return false;
}

2.2.3 關鍵點生成儲存

通過呼叫上述過程函式，將關鍵點位置尺度資訊儲存。
關鍵點資料結構：

struct Keypoint
{
    int octave;                 //【1】關鍵點所在組
    int interval;               //【2】關鍵點所在層

    double offset_interval;     //【3】調整後層的增量（是3/2嗎,怎麼是double?）

    int x, y;                   //【4】x,y座標,根據octave和interval可取的層內影象（是指所在組層？）
    double scale;               //【5】空間尺度座標scale = sigma0*pow(2.0, o+s/S)  

    double dx, dy;              //【6】特徵點座標，該座標被縮放成原影象大小
    double offset_x, offset_y;

    //============================================================  
    //特徵描述子
    //============================================================  
    double octave_scale;        //【1】offset_i
    double ori;                 //【2】方向
    double descriptor[FEATURE_ELEMENT_LENTH];//【3】特徵描述子
    int descr_length;
    double val;                 //【4】極值
};

儲存函式：

void DetectionLocalExtrema(vector<Mat>& dog_pyr,vector<Keypoint>& extrema,int intervals)
{
    int octaves = (int)dog_pyr.size() / (intervals+2);      //DOG組數
    double thresh=0.5*DXTHRESHHOLD/intervals;               
    for (int  o = 0; o < octaves; o++)
    {
        for (int i = 1; i < intervals+1; i++)
        {
            int index = o*(intervals + 2) + i;              //當前層的索引序號
            pixel_t *data =(pixel_t*) dog_pyr[index].data;  //當前層資料首地址
            //對每層進行提取特徵點
            for (int y = IMG_BORDER; y < dog_pyr[index].rows-IMG_BORDER; y++)
            {
                for (int x = IMG_BORDER; x < dog_pyr[index].cols-IMG_BORDER; x++)
                {
                    pixel_t val = *(data + y*dog_pyr[index].cols + x);
                    //當大於閾值，並且時是26箇中極值點時
                    if (fabs(val)>thresh&&isExtremum(x,y,dog_pyr,index))
                    {
                        Keypoint *extrum = InterpolationExtremum(x, y, dog_pyr, index, o, i);
                        if (extrum)
                        {
                            if (passEdgeResponse(extrum->x,extrum->y,dog_pyr,index))
                            {

                                extrum->val = At(index,extrum->x,extrum->y);
                                extrema.push_back(*extrum);
                            }
                            delete extrum;
                        }
                    }
                }
            }
        }
    }
}

2.2.4 尺度

void CalculateScale(vector<Keypoint>& features, double sigma, int intervals)
{
    double intvl = 0.0;
    for (int i = 0; i < features.size(); i++)
    {
        intvl = features[i].interval + features[i].offset_interval;
        features[i].scale = sigma * pow(2.0, features[i].octave + intvl / intervals);
        features[i].octave_scale = sigma * pow(2.0, intvl / intervals);
    }

}
//對擴大的影象特徵縮放，對應到原圖大小（非-1層）  
void HalfFeatures(vector<Keypoint>& features)
{
    for (int i = 0; i < features.size(); i++)
    {
        features[i].dx /= 2;
        features[i].dy /= 2;
        features[i].scale /= 2;
    }
}

相關推薦

SIFT演算法特徵描述子構建---關鍵點定位原理及程式碼

0.引言 sift針對區域性特徵進行特徵提取，在尺度空間尋找極值點，提取位置，尺度，旋轉不變數，生成特徵描述子。 總共分四個步驟： step2 關鍵點/極值點提取 2.1 關鍵點位置初步探查 生成DOG金字塔後，要找到DOG空間中的區域

sift演算法特徵描述子構建程式碼實現--梯度直方圖生成原理及程式碼

0.引言 sift針對區域性特徵進行特徵提取，在尺度空間尋找極值點，提取位置，尺度，旋轉不變數，生成特徵描述子。 總共分四個步驟： step3 生成梯度直方圖 生成特徵點的梯度資訊，並且確定主方向和輔助主方向的關鍵點。 3.1 梯度計算

SIFT演算法特徵描述子構建---特徵描述子構建原理及程式碼

0.引言 sift針對區域性特徵進行特徵提取，在尺度空間尋找極值點，提取位置，尺度，旋轉不變數，生成特徵描述子。 總共分四個步驟： 4 特徵描述子構建 每個關鍵點的方向、位置、尺度資訊都具備後，可以對區域性特徵進行描述，即特徵描述子。

自編碼器深度分析+定製特徵描述子構建初探

1. 關於tailored 特徵描述子         自從深度學習的繁榮，利用自編碼器、孿生網路、對抗網路採用資料驅動的形式設計tailored 特徵描述子 成為了計算機視覺領域發展的重要推動力， 這不僅大大削弱了特徵工程的壓力，而且降低了相關領域學

【機器學習--opencv3.4.1版本基於Hog特徵描述子Svm對經典手寫數字識別】

 方向梯度直方圖（Histogram of Oriented Gradient, HOG）特徵是一種在計算機視覺和影象處理中用來進行物體檢測的特徵描述子。HOG特徵通過計算和統計影象區域性區域的梯度方向直方圖來構成特徵。 #include <iostream> #inc

3D特徵描述子SHOT詳解

在3D點雲目標識別中，除了傳統的利用二維影象來輔助識別，然後對映到3D的方式；學術界在十幾年前就開始研究3D點雲的特徵描述子。 本文主要詳解義大利博洛尼亞大學教授提出的SHOT（Signature of Histogram of  Orientation）http://www.vision.dei

sift 計算機視覺——描述子

描述子實現程式碼 這裡使用開源工具包VLFeat提供的二進位制檔案來計算影象的SIFT特徵。用完整的Python實現SIFT特徵的所有步驟可能效率不是很高。VLFeat工具包可以從http://www.vlfeat.org/下載，二進位制檔案可以在所有主要的平臺上執行。VLFeat庫是用C語言來寫

OpenCv-C++-特徵描述子（BruteForce（暴力匹配））

本片文主要是通過BFMatcher()來進行特徵匹配，暴力匹配就是在一個數組或連結串列裡面進行一一查詢，迴圈匹配。如下圖： 通過SURF運算元分別找到兩張影象的特徵點，再使用BFMatcher函式對兩張影象上的特徵點進行一一檢測和計算，最後通過drawMatches函式繪製兩張影象之間

BRIEF 特徵描述子

Binary Robust Independent Elementary Features 1. BRIEF的基本原理 我們已經知道SIFT特徵採用了128維的特徵描述子，由於描述子用的浮點數，所以它將會佔用512 bytes的空間。類似地，對於SURF特徵，常見的是64維的描述子，它也將佔用256b

從特徵描述子到深度學習：計算機視覺的20年曆程 From feature descriptors to deep learning: 20 years of computer vision

We all know that deep convolutional neural networks have produced some stellar results on object detection and recognition benchmarks in the past two year

特徵描述子評測

在高斯模糊的測試中，我模擬的運動模糊（當攝像機突然移動時會發生的），在這項測試中，所有的描述子都表現出很好的結果，說結果好，我的意思是更大的高斯核大小，會導致更小的正確匹配百分比 Lighting test 在光照變化測試中，影象的不同僅僅在於整幅影象的亮度，所有的描述子都表現不錯。主要的原因是所有

影象處理-常用特徵點檢測方法以及特徵描述子

HOG(histogram of gradient) sift(Scale-invariant feature transform) sift是用於描述影象中的區域性特徵，在空間尺度中尋找極值點，並且提取出其位置、尺度、旋轉不變數，因此具有尺度和旋

角點檢測演算法及其描述子評估介紹

我以前的部落格裡面介紹了大量的角點檢測演算法和其區域性描述子設計方法，這裡將摘錄其各自的論文和相關綜述文章內的實驗估評估，然後給出自己的一些總結。（特別多的圖>O<） 一、角點檢測子的效果評估 1.1 Harris同SIFT、FAST、SURF等角點檢測方法匹

【特徵匹配】BRIEF特徵描述子原理及原始碼解析

轉載請註明出處： http://blog.csdn.net/luoshixian099/article/details/48338273   傳統的特徵點描述子如SIFT,SURF描述子，每個特徵點採用128維(SIFT)或者64維(SURF)向量去描述，每個維度上

影象特徵描述子之ORB

  ORB(Oriented FAST and Rotated BRIEF)演算法是對FAST特徵點檢測和BRIEF特徵描述子的一種結合，在原有的基礎上做了改進與優化，使得ORB特徵具備多種區域性不變性，併為實時計算提供了可能。 特徵點檢測   OR

特徵描述子（-）—HOG具體實現過程

       方向梯度直方圖（Histogram of Oriented Gradient, HOG）特徵是一種在計算機視覺和影象處理中用來進行物體檢測的特徵描述子。它通過計算和統計影象區域性區域的梯度方向直方圖來構成特徵。其提取演算法具體實現過程如下： 1.影象灰度化處

【opencv基於特徵描述子的物體檢驗尋找】

關於特徵子中經常用的的函式：findHomography，estimateRigidTransform 關於opencv的特徵描述子，我們需要了解DMatch與KPoints的資料結構，如下： DMatch結構體 /* * Struct for matching:

三維計算機視覺（五）--特徵描述子

PPF 機器人視覺中有一項重要人物就是從場景中提取物體的位置，姿態。影象處理演算法藉助Deep Learning 的東風已經在影象的物體標記領域耍的飛起了。而從三維場景中提取物體還有待研究。目前已有的思路是先提取關鍵點，再使用各種區域性特徵描述子對關鍵點進行描述，最後與待檢

特徵描述子（feature descriptor） —— HOG（方向梯度直方圖）

HOG（Histogram of Oriented Gradients），描述的是影象的區域性特徵，其命名也暗示了其計算方法，先計算影象中某一區域不同方向上梯度的值，然後累積計算頻次，得到直方圖，該直方圖便可代表該區域了，也即從影象中抽取得到的特徵向量，可以作為

影象特徵描述子之BRIEF

  BRIEF(Binary Robust Independent Elementary Features)是一種對已檢測到的特徵點進行表示和描述的特徵描述方法，和傳統的利用影象區域性鄰域的灰度直方圖或梯度直方圖提取特徵的方式不同，BRIEF是一種二進位制編