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