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