《SIFT原理与源码分析》系列文章索引：http://www.cnblogs.com/tianyalu/p/5467813.html

由前一篇《关键点搜索与定位》，我们已经找到了关键点。为了实现图像旋转不变性，需要根据检测到的关键点局部图像结构为特征点方向赋值。也就是在findScaleSpaceExtrema()函数里看到的alcOrientationHist()语句：

// 计算梯度直方图

float omax = calcOrientationHist(gauss_pyr[o*(nOctaveLayers+) + layer],

                                                Point(c1, r1),

                                                cvRound(SIFT_ORI_RADIUS * scl_octv),

                                                SIFT_ORI_SIG_FCTR * scl_octv,

                                                hist, n);

我们使用图像的梯度直方图法求关键点局部结构的稳定方向。

梯度方向和幅值

在前文中，精确定位关键点后也找到改特征点的尺度值σ，根据这一尺度值，得到最接近这一尺度值的高斯图像：

使用有限差分，计算以关键点为中心，以3×1.5σ为半径的区域内图像梯度的幅角和幅值，公式如下：

梯度直方图

在完成关键点邻域内高斯图像梯度计算后，使用直方图统计邻域内像素对应的梯度方向和幅值。

有关直方图的基础知识可以参考《数字图像直方图》，可以看做是离散点的概率表示形式。此处方向直方图的核心是统计以关键点为原点，一定区域内的图像像素点对关键点方向生成所作的贡献。

梯度方向直方图的横轴是梯度方向角，纵轴是剃度方向角对应的梯度幅值累加值。梯度方向直方图将0°~360°的范围分为36个柱，每10°为一个柱。下图是从高斯图像上求取梯度，再由梯度得到梯度方向直方图的例图。

在计算直方图时，每个加入直方图的采样点都使用圆形高斯函数函数进行了加权处理，也就是进行高斯平滑。这主要是因为SIFT算法只考虑了尺度和旋转不变形，没有考虑仿射不变性。通过高斯平滑，可以使关键点附近的梯度幅值有较大权重，从而部分弥补没考虑仿射不变形产生的特征点不稳定。

通常离散的梯度直方图要进行插值拟合处理，以求取更精确的方向角度值。（这和《关键点搜索与定位》中插值的思路是一样的）。

关键点方向

直方图峰值代表该关键点处邻域内图像梯度的主方向，也就是该关键点的主方向。在梯度方向直方图中，当存在另一个相当于主峰值 80%能量的峰值时，则将这个方向认为是该关键点的辅方向。所以一个关键点可能检测得到多个方向，这可以增强匹配的鲁棒性。Lowe的论文指出大概有15%关键点具有多方向，但这些点对匹配的稳定性至为关键。

获得图像关键点主方向后，每个关键点有三个信息(x,y,σ,θ)：位置、尺度、方向。由此我们可以确定一个SIFT特征区域。通常使用一个带箭头的圆或直接使用箭头表示SIFT区域的三个值：中心表示特征点位置，半径表示关键点尺度（r=2.5σ）,箭头表示主方向。具有多个方向的关键点可以复制成多份，然后将方向值分别赋给复制后的关键点。如下图：

源码

// Computes a gradient orientation histogram at a specified pixel

// 计算特定点的梯度方向直方图

static float calcOrientationHist( const Mat& img, Point pt, int radius,

                                  float sigma, float* hist, int n )

{

    //len：2r+1也就是以r为半径的圆（正方形）像素个数

    int i, j, k, len = (radius*+)*(radius*+);  

    float expf_scale = -.f/(.f * sigma * sigma);

    AutoBuffer<float> buf(len* + n+);

    float *X = buf, *Y = X + len, *Mag = X, *Ori = Y + len, *W = Ori + len;

    float* temphist = W + len + ;  

    for( i = ; i < n; i++ )

        temphist[i] = .f;  

    // 图像梯度直方图统计的像素范围

    for( i = -radius, k = ; i <= radius; i++ )

    {

        int y = pt.y + i;

        if( y <=  || y >= img.rows -  )

            continue;

        for( j = -radius; j <= radius; j++ )

        {

            int x = pt.x + j;

            if( x <=  || x >= img.cols -  )

                continue;  

            float dx = (float)(img.at<short>(y, x+) - img.at<short>(y, x-));

            float dy = (float)(img.at<short>(y-, x) - img.at<short>(y+, x));  

            X[k] = dx; Y[k] = dy; W[k] = (i*i + j*j)*expf_scale;

            k++;

        }

    }  

    len = k;  

    // compute gradient values, orientations and the weights over the pixel neighborhood

    exp(W, W, len);

    fastAtan2(Y, X, Ori, len, true);

    magnitude(X, Y, Mag, len);   

    // 计算直方图的每个bin

    for( k = ; k < len; k++ )

    {

        int bin = cvRound((n/.f)*Ori[k]);

        if( bin >= n )

            bin -= n;

        if( bin <  )

            bin += n;

        temphist[bin] += W[k]*Mag[k];

    }  

    // smooth the histogram

    // 高斯平滑

    temphist[-] = temphist[n-];

    temphist[-] = temphist[n-];

    temphist[n] = temphist[];

    temphist[n+] = temphist[];

    for( i = ; i < n; i++ )

    {

        hist[i] = (temphist[i-] + temphist[i+])*(.f/.f) +

            (temphist[i-] + temphist[i+])*(.f/.f) +

            temphist[i]*(.f/.f);

    }  

    // 得到主方向

    float maxval = hist[];

    for( i = ; i < n; i++ )

        maxval = std::max(maxval, hist[i]);  

    return maxval;

}