0.引言

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

總共分四個步驟：

step3 生成梯度直方圖

生成特徵點的梯度資訊，並且確定主方向和輔助主方向的關鍵點。

3.1 梯度計算

經過第二步驟，關鍵點已經有了尺度和位置資訊，缺少的梯度方向資訊。首先計算梯度。

/*****************************************************************************************************************************
 
***
*模組說明:
*        模組六---步驟1：計算關鍵點的梯度和梯度方向
*功能說明：
*        1）計算關鍵點(x,y)處的梯度幅值和梯度方向
*        2）將所計算出來的梯度幅值和梯度方向儲存在變數mag和ori中
*********************************************************************************************************************************/
bool CalcGradMagOri(const Mat& gauss, int x, int y, double& mag, double& ori)
{
    if (x > 0 && x < gauss.cols - 1 && y > 0 && y < gauss.rows - 1)
 

    {
        pixel_t *data = (pixel_t*)gauss.data;
        double dx = *(data + gauss.cols*y + (x + 1)) - (*(data + gauss.cols*y + (x - 1)));           //[1]利用X方向上的差分代替微分dx  
        double dy = *(data + gauss.cols*(y + 1) + x) - (*(data + gauss.cols*(y - 1) + x));           //[2]利用Y方向上的差分代替微分dy  

        mag = sqrt(dx*dx + dy*dy);                                          //[3]計算該關鍵點的梯度幅值  
 

        ori = atan2(dy, dx);                                                //[4]計算該關鍵點的梯度方向  弧度(0，-1)-pi~pi（-0.0001，-1）  pi/2(1,0)
        return true;
    }
    else
        return false;
}

3.2 計算梯度的方向直方圖

梯度直方圖將0~360度的方向範圍分為36個柱(bins)，其中每柱10度。

/********************************************************************************************************************************
*模組說明:
*        模組六：2.計算梯度的方向直方圖
*功能說明：
*        1）直方圖以每10度為一個柱，共36個柱，柱代表的方向為為畫素點的梯度方向，柱的長短代表了梯度幅值。
*        2）根據Lowe的建議，直方圖統計採用3*1.5*sigma
*結    論：
*        影象的關鍵點檢測完畢後，每個關鍵點就擁有三個資訊：位置、尺度、方向；同時也就使關鍵點具備平移、縮放和旋轉不變性
*********************************************************************************************************************************/
double* CalculateOrientationHistogram(const Mat&gauss,int x,int y,int bins,int radius,double sigma)
{
    double * hist = new double[bins];
    for (int i = 0; i < bins; i++)
    {
        *(hist + i) = 0.0;
    }
    double mag, ori, weight;            //[3]關鍵點的梯度幅值、方向、梯度權重
    int bin;                                //第bin個柱
    double econs = -1.0 / (2.0*sigma*sigma);//高斯平滑指數
    for (int i = 0; i < radius; i++)
    {
        for (int j = 0; j < radius; j++)
        {
            if (CalcGradMagOri(gauss,x+i,y+j,mag,ori))
            {
                weight = exp((i*i + j*j)*econs);// 使用圓形高斯函式函式進行了加權處理，也就是進行高斯平滑
                //-pi->0->pi    角度範圍
                //36->18->0     柱序號   從x軸負方向順時針轉一圈。
                bin = cvRound(bins*(CV_PI - ori) / (2 * CV_PI));    //計算第幾個bin
                bin = bin < bins ? bin : 0;     //計算在哪個塊內
                hist[bin] += mag*weight;        //統計梯度的方向直方圖  
            }
        }
    }
    return hist;
}

3.3 對方向直方圖高斯平滑

在直方圖統計時，每相鄰三個畫素點採用高斯加權，根據Lowe的建議，模板採用[0.25,0.5,0.25],並且連續加權兩次.

void GaussSmoothOriHist(double* hist, int n)
{
    double prev = hist[n - 1];
    double temp;
    double h0 = hist[0];
    for (int i = 0; i < n; i++)
    {
        temp = hist[i];
        //如果是最後一個，則跟第零個平滑
        hist[i] = 0.25*prev + 0.5*hist[i] + 0.25*(i + 1>=n ? h0 : hist[i + 1]);
        prev = temp;
    }
}

3.4 計算主方向和輔方向

方向直方圖的峰值則代表了該特徵點的方向，以直方圖中的最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，保留峰值大於主方向峰值80%的方向作為該關鍵點的輔方向。

ouble DominantDirection(double* hist,int n)
{
    double maxd = hist[0];
    for (int i = 0; i < n; i++)
    {
        if (hist[i]>maxd)
            maxd = hist[i];             //求36個柱中最大峰值
    }
    return maxd;
}

/********************************************************************************************************************************
*模組說明:
*        模組六---步驟5：計算更加精確的關鍵點主方向----拋物插值
*功能說明：
*        1）方向直方圖的峰值則代表了該特徵點的方向，以直方圖中的最大值作為該關鍵點的主方向。為了增強匹配的魯棒性，只保留峰值大於主
*           方向峰值80%的方向作為該關鍵點的輔方向。因此，對於同一梯度值得多個峰值的關鍵點位置，在相同位置和尺度將會有多個關鍵點被
*           建立但方向不同。僅有15%的關鍵點被賦予多個方向，但是可以明顯的提高關鍵點的穩定性。
*        2）在實際程式設計中，就是把該關鍵點複製成多份關鍵點，並將方向值分別賦給這些複製後的關鍵點
*        3）並且，離散的梯度直方圖要進行【插值擬合處理】，來求得更加精確的方向角度值
*********************************************************************************************************************************/
//複製關鍵點
void CopyKeypoint(const Keypoint& src, Keypoint& dst)
{
    dst.dx = src.dx;
    dst.dy = src.dy;

    dst.interval = src.interval;
    dst.octave = src.octave;
    dst.octave_scale = src.octave_scale;
    dst.offset_interval = src.offset_interval;

    dst.offset_x = src.offset_x;
    dst.offset_y = src.offset_y;

    dst.ori = src.ori;
    dst.scale = src.scale;
    dst.val = src.val;
    dst.x = src.x;
    dst.y = src.y;
}
#define Parabola_Interpolate(l,c,r) (0.5*(l-r)/(l-2*c+r))
void CalcOriFeatures(const Keypoint& keypoint,vector<Keypoint>& features,const double *hist,int n,double mag_thr)
{
    double bin;
    int l, r;
    //36個
    for (int  i = 0; i <= n; i++)
    {
        l = (i == 0) ? (n - 1) : i - 1;
        r = (i + 1) % n;
        if (hist[i]>hist[l]&&hist[i]>hist[r]&&hist[i]>mag_thr)
        {
            //插值
            bin = i + Parabola_Interpolate(hist[l],hist[i],hist[r]);
            bin = (bin < 0) ? (bin + n) : (bin >= n ? bin - n : bin);

            Keypoint new_key;
            CopyKeypoint(keypoint, new_key);
            new_key.ori = (2 * CV_PI*bin / n) - CV_PI;
            features.push_back(new_key);
        }
    }
}

3.5 生成帶有梯度資訊的關鍵點

void OrientationAssignment(vector<Keypoint>& extrema,vector<Keypoint>& features,const vector<Mat> &guass_pyr)
{
    int n = extrema.size();
    double* hist;
    for (int i = 0; i < n; i++)
    {
        //計算該關鍵點處的直方圖
        hist = CalculateOrientationHistogram(guass_pyr[extrema[i].octave*(INTERVALS + 3) + extrema[i].interval], extrema[i].x, extrema[i].y, ORI_HIST_BINS,
            cvRound(ORI_PEAK_RATIO*extrema[i].octave_scale), ORI_SIGMA_TIMES*extrema[i].octave_scale);
        //高斯平滑
        for (int j = 0; j <ORI_SMOOTH_TIMES; j++)
        {
            GaussSmoothOriHist(hist, ORI_HIST_BINS);
        }
        //取極值點
        double highest_peak = DominantDirection(hist, ORI_HIST_BINS);
        CalcOriFeatures(extrema[i], features, hist, ORI_HIST_BINS, highest_peak*ORI_PEAK_RATIO);
        delete[] hist;
    }
}