sift计算描述子代码详解_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，1void 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

{

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

{

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 

{

for (int 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-1 && ybin

{

//计算关键点的梯度

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 

{

len_sq += feat.descriptor[i] * feat.descriptor[i];

}

len_sq = sqrt(len_sq);    for (int i = 0; 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 

{        for (int c = 0; c 

{            for (int o = 0; o 

{

feature.descriptor[k++] = hist[r][c][o];//放进128维特征描述子内

}

}

}

feature.descr_length = k;    //描述子归一化

NormalizeDescr(feature);    //描述子门限

for (int i = 0; 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 

{

int_val = INT_DESCR_FCTR*feature.descriptor[i];

feature.descriptor[i] = min(255, int_val);

}

}

4.5 生成描述子void DescriptorRepresentation(vector&features, const vector& guass_pyr, int bins, int width)

{    double ***hist;    for (int i = 0; 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 

{            for (int k = 0; k 

{                delete[] hist[j][k];

}            delete[] hist[j];

}        delete[] hist;

}

}