SIFT算法特征描述子构建---特征描述子构建原理及代码

0.引言

sift针对局部特征进行特征提取，在尺度空间寻找极值点，提取位置，尺度，旋转不变量，生成特征描述子。

总共分四个步骤：

• 尺度金字塔生成
• 关键点/极值点提取
• 生成梯度直方图

• 特征描述子构建

  4 特征描述子构建

  每个关键点的方向、位置、尺度信息都具备后，可以对局部特征进行描述，即特征描述子。

  4.1 确定描述子区域

  将关键点划分为d*d(Lowe建议4)个子区域，每个子区域为一个种子点，
  每个种子点有8个方向，即128维特征。为每个子区域分配边长为3*sigma_oct的矩形采样，考虑实际计算用双线性插值，以及旋转，放大sqrt(2)*(d+1),
  最终所选图像窗口半径

4.2 坐标轴旋转为关键点的方向

旋转保持旋转不变性

4.3 将子区域内的梯度值插值加权分配到8 个方向上

//对每个种子点每个方向插值
//weight=w*dr^k*(1-dr)^（1-k）*dc^m*(1-dc)^(1-m)*do^n*(1-do)^(1-n)
//m,n,k为0，1
void InterpHistEntry(double*** hist, double xbin, double ybin, double obin, double mag, int bins, int d)
{
    //对邻近两行的贡献因子为dr和1 - dr，对邻近两列的贡献因子为dc和1 - dc，对邻近两个方向的贡献因子为do和1 - do
    double d_r, d_c, d_o;
    int r0, c0, o0;
    //r0取不大于xbin的正整数时。
    //r0 + 0 <= xbin <= r0 + 1
    //mag在区间[r0, r1]上做插值
    r0 = cvFloor(ybin);
    c0 = cvFloor(xbin);
    o0 = cvFloor(obin);
    d_r = ybin - r0; // 计算偏差
    d_c = xbin - c0;
    d_o = obin - o0;

    double v_r, v_c, v_o;//v_r为mag在行上、列加权，最后方向的加权
    int rb, cb, ob;//r,c,o分别取0，1时候的行列
    double** row, *h;
    for (int r = 0; r <= 1; r++)
    {
        rb = r + r0;
        if (rb >= 0 && rb<d)//在4*4描述子内
        {
            v_r = mag*((r == 0) ? 1 - d_r : d_r);//分别计算r取0，1时的值，最后相加
            row = hist[rb];//固定行，进入hist的列
            for (int c = 0; c <= 1; c++)
            {
                cb = c0 + c;
                if (cb >= 0 && cb<d)
                {
                    v_c = v_r*((c == 0) ? 1 - d_c : d_c);
                    h = row[cb];
                    for (int o = 0; o <= 1; o++)
                    {
                        ob = (o0 + o) % bins;
                        v_o = v_c*((o == 0) ? 1 - d_o : d_o);
                        h[ob] += v_o;//累加各层取0，1时候的加权值
                    }
                }

            }
        }

    }
}

double*** CalculateDescrHist(const Mat&gauss, int x, int y, double octave_scale, double ori, int bins, int width)
{
//二维数组指针,width 为子区域尺寸d=4*4
double *hist = new double [width];
for (int i = 0; i < width; i++)
{
//长度为width的一维数组指针
hist[i] = new double *[width];
for (int j = 0; j < width; j++)
{
//每个hist处是一个36维数组
hist[i][j] = new double[bins];
}
}
//初始化
for (int r = 0; r < width; r++)
for (int c = 0; c < width; c++)
for (int o = 0; o < bins; o++)
hist[r][c][o] = 0.0;

//高斯权值，Lowe建议子区域的像素的梯度大小按sigma=0.5*d的高斯加权计算，即2
double sigma = 0.5*width;
double conste = -1.0 / (2 * sigma*sigma);

double sub_hist_width = DESCR_SCALE_ADJUST*octave_scale;//每个子区域尺寸为mσ个像元 尺度特征点的尺度值3*sig_oct
                                                        //子区域半径
int radius = cvRound((sqrt(2)*sub_hist_width*(width + 1)) / 2.0);
double grad_ori;
double grad_mag;
/*计算采样区域点坐标旋转
|x`|  |cos -sin| |x|
|  | =|        |*| |
|y`|  |sin  cos| |y|
子区域下标
| x``|  1               |x|
|    |=--------------*  | | +1/d
| y``| sub_hist_width   |y|
*/
double cos_ori = cos(ori);
double sin_ori = sin(ori);
for (int i = -radius; i < radius; i++)
{
    for (int j = -radius; j < radius; j++)
    {
        double rot_x = (cos_ori*j - sin_ori + i);
        double rot_y = (sin_ori*j + cos_ori + i);   
        double xbin = rot_x / sub_hist_width + width / 2 - 0.5;                     //xbin, ybin为落在4 * 4窗口中的下标值
        double ybin = rot_y / sub_hist_width + width / 2 - 0.5;
        if (xbin>-1.0&&xbin<width&&ybin>-1 && ybin<width)
        {
            //计算关键点的梯度
            if (CalcGradMagOri(gauss, x + j, y + i, grad_mag, grad_ori))
            {
                //梯度方向夹角
                grad_ori = (CV_PI - grad_ori) - ori;
                while (grad_ori<0)
                {
                    grad_ori += CV_PI * 2;
                }
                while (grad_ori >= CV_PI * 2)
                {
                    grad_ori -= 2 * CV_PI;
                }
                double obin = grad_ori*(bins / (2 * CV_PI));//种子点所在子窗口的方向 
                                                            //公式子区域像素梯度进行高斯加权：exp(-（(x`2)+（y`2））/(2*(0.5d)^2))
                double weight = exp(conste*(rot_x*rot_x + rot_y*rot_y));

                //插值计算每个种子点处的梯度
                InterpHistEntry(hist, xbin, ybin, obin, grad_mag*weight, bins, width);
            }
        }

    }
}


return hist;

}

4.4 归一化及门限

如上统计的4*4*8=128个梯度信息即为该关键点的特征向量。特征向量形成后，为了去除光照变化的影响，需要对它们进行归一化处理，对于图像灰度值整体漂移，图像各点的梯度是邻域像素相减得到，所以也能去除。

void NormalizeDescr(Keypoint& feat)
{
    double len_sq = 0;
    for (int i = 0; i < feat.descr_length; i++)
    {
        len_sq += feat.descriptor[i] * feat.descriptor[i];
    }
    len_sq = sqrt(len_sq);
    for (int i = 0; i < feat.descr_length; i++)
    {
        feat.descriptor[i] = feat.descriptor[i] / len_sq;
    }
}

设置门限
非线性光照，相机饱和度变化对造成某些方向的梯度值过大，而对方向的影响微弱。因此设置门限值(向量归一化后，一般取0.2)截断较大的梯度值。然后，再进行一次归一化处理，提高特征的鉴别性。

//直方图存入特征描述子
void HistToDescriptor(double***hist, int width, int bins, Keypoint& feature)
{
    int k = 0;
    for (int r = 0; r < width; r++)
    {
        for (int c = 0; c < width; c++)
        {
            for (int o = 0; o < bins; o++)
            {
                feature.descriptor[k++] = hist[r][c][o];//放进128维特征描述子内
            }
        }
    }
    feature.descr_length = k;
    //描述子归一化
    NormalizeDescr(feature);
    //描述子门限
    for (int i = 0; i < k; i++)
    {
        if (feature.descriptor[i]>DESCR_MAG_THR)
        {
            feature.descriptor[i] = DESCR_MAG_THR;
        }
    }
    //第二次归一化
    NormalizeDescr(feature);
    int int_val;//整形值  将double型转为整形描述子
    for (int i = 0; i < k; i++)
    {
        int_val = INT_DESCR_FCTR*feature.descriptor[i];
        feature.descriptor[i] = min(255, int_val);
    }

}

4.5 生成描述子

void DescriptorRepresentation(vector<Keypoint>&features, const vector<Mat>& guass_pyr, int bins, int width)
{
    double ***hist;
    for (int i = 0; i < features.size(); i++)
    {
        hist = CalculateDescrHist(guass_pyr[features[i].octave*(INTERVALS + 3) + features[i].interval], features[i].x, features[i].y, features[i].octave_scale, features[i].ori, bins, width);
        HistToDescriptor(hist, width, bins, features[i]);
        //释放空间
        for (int j = 0; j < width; j++)
        {
            for (int k = 0; k < width; k++)
            {
                delete[] hist[j][k];
            }
            delete[] hist[j];
        }
        delete[] hist;
    }
}