sift算法特征描述子构建代码实现--梯度直方图生成原理及代码

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;
    }
}