SIFT解析（二）特征点位置确定
最近微博上有人发起投票那篇论文是自己最受益匪浅的论文，不少人说是lowe的这篇介绍SIFT的论文。确实，在图像特征识别领域，SIFT的出现是具有重大意义的，SIFT特征以其稳定的存在，较高的区分度推进了诸多领域的发展，比如识别和配准。上一篇文章，解析了SIFT特征提取的第一步高斯金字塔的构建，并详细分析了高斯金字塔以及差分高斯金字塔如何完成一个连续的尺度空间的构建。构建高斯金字塔不是目的，目的是如何利用高斯金字塔找到极值点。
lowe在论文中阐述了为什么使用差分高斯金字塔：
1）差分高斯图像可以直接由高斯图像相减获得，简单高效
2） 差分高斯函数是尺度规范化的高斯拉普拉斯函数的近似，而高斯拉普拉斯函数的极大值和极小值点是一种非常稳定的特征点（与梯度特征、Hessian特征和Harris角点相比）
有了这些基础，我们就可以放开手脚从差分高斯金字塔中找点了。
特征点的确定主要包括两个过程：确定潜在特征点，精确确定特征点的位置和去除不稳定特征点。
确定潜在特征点
上文已经阐述，高斯拉普拉斯函数的极大值和极小值点是一种非常稳定的特征点，因此我们从差分高斯金字塔中寻找这些潜在特征点。差分高斯金字塔是一个三维空间（平面图像二维，尺度一维），因此我们在三维空间中在寻找极大值点和极小值点。具体方法是比较当前特征点的灰度值和其他26个点的灰度值的大小，这26个点包括：当前尺度下该点的8邻域以及前一尺度和后一尺度下与该点最近的9个点（9*2+8=26），如下图所示：
OpenCV该部分源码：
- void SIFT::findScaleSpaceExtrema( const vector<Mat>& gauss_pyr, const vector<Mat>& dog_pyr,
- vector<KeyPoint>& keypoints ) const
- {
- ......
- for( int o = 0; o < nOctaves; o++ )//每一个八度
- for( int i = 1; i <= nOctaveLayers; i++ )//对八度中的存在具有第1至第nOctaveLayers层高斯差分图像提取特征点
- {
- ......
- for( int r = SIFT_IMG_BORDER; r < rows-SIFT_IMG_BORDER; r++)//图像二维空间.行
- {
- ......
- for( int c = SIFT_IMG_BORDER; c < cols-SIFT_IMG_BORDER; c++)//图像二维空间.列
- {
- .......
- // 当前点与26个点比较，比较两次，分别确定是否是极大值，是否是极小值
- if( std::abs(val) > threshold &&
- ((val > 0 && val >= currptr[c-1] && val >= currptr[c+1] &&
- val >= currptr[c-step-1] && val >= currptr[c-step] && val >= currptr[c-step+1] &&
- val >= currptr[c+step-1] && val >= currptr[c+step] && val >= currptr[c+step+1] &&
- val >= nextptr[c] && val >= nextptr[c-1] && val >= nextptr[c+1] &&
- val >= nextptr[c-step-1] && val >= nextptr[c-step] && val >= nextptr[c-step+1] &&
- val >= nextptr[c+step-1] && val >= nextptr[c+step] && val >= nextptr[c+step+1] &&
- val >= prevptr[c] && val >= prevptr[c-1] && val >= prevptr[c+1] &&
- val >= prevptr[c-step-1] && val >= prevptr[c-step] && val >= prevptr[c-step+1] &&
- val >= prevptr[c+step-1] && val >= prevptr[c+step] && val >= prevptr[c+step+1]) ||
- (val < 0 && val <= currptr[c-1] && val <= currptr[c+1] &&
- val <= currptr[c-step-1] && val <= currptr[c-step] && val <= currptr[c-step+1] &&
- val <= currptr[c+step-1] && val <= currptr[c+step] && val <= currptr[c+step+1] &&
- val <= nextptr[c] && val <= nextptr[c-1] && val <= nextptr[c+1] &&
- val <= nextptr[c-step-1] && val <= nextptr[c-step] && val <= nextptr[c-step+1] &&
- val <= nextptr[c+step-1] && val <= nextptr[c+step] && val <= nextptr[c+step+1] &&
- val <= prevptr[c] && val <= prevptr[c-1] && val <= prevptr[c+1] &&
- val <= prevptr[c-step-1] && val <= prevptr[c-step] && val <= prevptr[c-step+1] &&
- val <= prevptr[c+step-1] && val <= prevptr[c+step] && val <= prevptr[c+step+1])))
- {
- ......
- }
- }
- }
- }
- }
尺度空间中的极值点已经确定出来了，下面有两个问题需要解决：
（1）这些点是最终我们确定的SIFT特征点集的超集，该超集里包含许多“间谍”-----不稳定的特征点，因此必须去掉这些不稳定的特征点。这些不稳定的特征点主要包含两类：低对比度的点（对噪声敏感）和边缘点。
（2）这一步骤中极值点的坐标还是离散的整数值，如何精确确定特征点的位置。
由于在计算上（2）问题的解决可以捎带解决（1）中低对比度点的问题，因此我们先讨论问题（2）。本部分的OpenCV源码位于sift.cpp文件的adjustLocalExtrema函数中，本文最后会贴出此部分源码，下面首先分析如何解决以上两个问题。
精确确定特征点的位置：
由于图像是一个离散的空间，特征点的位置的坐标都是整数，但是极值点的坐标并不一定就是整数，如下图所示。
因此，如何从离散空间中估计出极值点的精确位置是重要的。为了精确确定极值点坐标，Brown和Lowe使用了三元二次函数，通过迭代确定极值点的位置，具有良好的效果。
主要是根据泰勒公式，泰勒公式作用：用值已知的点A估计点A附近的某点B的值。
求上式极值，对其求导，导数等于0，得到
去除不稳定特征点
去除对比度低的点
以上求出了极值点的精确的位置，将求出的 x 带入原式，得：
我们就利用这个函数去除对比度低的点，lowe文中，当D（x）<=0.03时，去除这个特征点。
去除边缘点
差分高斯金字塔中的极值点会有许多边缘点，边缘点对一些噪声不稳定，因此需要去除这些边缘相应点。
差分高斯金字塔中会有一些不是很好的极值点，这些点的特征是：在跨越边缘的方向有较大的主曲率，在与边缘相切的方向主曲率较小。在本步骤中，需要去除这些不好的边缘相应。主曲率可以通过2阶Hessian方阵获得：
D函数中某点的主曲率和该点的H矩阵的特征值是成比例的，因此我们可以通过H矩阵的特征值来确定某点在差分高斯金字塔中的主曲率。
设矩阵H的特征值分别为α（较大）和β（较小），有如下公式：
通过以上两式，α和β就可以计算出来了，但是，不急！
如上文所述，那些不好的边缘点：跨越边缘的方向有较大的主曲率，与边缘相切的方向主曲率较小。因此，我们通过α/β的比率函数并确定阈值来体现表征那些不好的边缘点，α/β越大，说明这个点就越糟糕，就越应该被删掉，但是这样就要真真切切计算α和β的值，前面让大家不急了，是的，先不用着急计算，设定r=α/β(即 α=rβ)，使用如下公式：
以上函数是关于r的增函数（已经假设α是特征值中较大的一个），r 越大，以上函数值就越大，反之，以上函数值越大，r 就是越大的，因此我们可以通过已知的Tr（H）和Det（H）“曲线地”去判断 r的大小！所以在本步骤中，去除不好的边缘点的阈值是：
lowe论文中设定r=10。
到这里，在差分高斯金字塔中提取的特征点就完成了提纯的步骤。
下面是OpenCV源码中特征点精确位置的确定过程以及特征点提纯过程，主要实现函数为sift.cpp中adjustLocalExtrema函数：
- // Interpolates a scale-space extremum's location and scale to subpixel
- // accuracy to form an image feature. Rejects features with low contrast.
- // Based on Section 4 of Lowe's paper.
- static bool adjustLocalExtrema( const vector<Mat>& dog_pyr, KeyPoint& kpt, int octv,
- int& layer, int& r, int& c, int nOctaveLayers,
- float contrastThreshold, float edgeThreshold, float sigma )
- {
- const float img_scale = 1.f/(255*SIFT_FIXPT_SCALE);
- const float deriv_scale = img_scale*0.5f;
- const float second_deriv_scale = img_scale;
- const float cross_deriv_scale = img_scale*0.25f;
- float xi=0, xr=0, xc=0, contr=0;
- int i = 0;
- // 如上文所述，迭代计算特征点的精确位置
- for( ; i < SIFT_MAX_INTERP_STEPS; i++ )
- {
- int idx = octv*(nOctaveLayers+2) + layer;
- const Mat& img = dog_pyr[idx];
- const Mat& prev = dog_pyr[idx-1];
- const Mat& next = dog_pyr[idx+1];
- Vec3f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale,
- (img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale,
- (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale);
- float v2 = (float)img.at<sift_wt>(r, c)*2;
- float dxx = (img.at<sift_wt>(r, c+1) + img.at<sift_wt>(r, c-1) - v2)*second_deriv_scale;
- float dyy = (img.at<sift_wt>(r+1, c) + img.at<sift_wt>(r-1, c) - v2)*second_deriv_scale;
- float dss = (next.at<sift_wt>(r, c) + prev.at<sift_wt>(r, c) - v2)*second_deriv_scale;
- float dxy = (img.at<sift_wt>(r+1, c+1) - img.at<sift_wt>(r+1, c-1) -
- img.at<sift_wt>(r-1, c+1) + img.at<sift_wt>(r-1, c-1))*cross_deriv_scale;
- float dxs = (next.at<sift_wt>(r, c+1) - next.at<sift_wt>(r, c-1) -
- prev.at<sift_wt>(r, c+1) + prev.at<sift_wt>(r, c-1))*cross_deriv_scale;
- float dys = (next.at<sift_wt>(r+1, c) - next.at<sift_wt>(r-1, c) -
- prev.at<sift_wt>(r+1, c) + prev.at<sift_wt>(r-1, c))*cross_deriv_scale;
- Matx33f H(dxx, dxy, dxs,
- dxy, dyy, dys,
- dxs, dys, dss);//通过当前像素点以及周围像素点差值出H矩阵
- Vec3f X = H.solve(dD, DECOMP_LU);
- xi = -X[2];
- xr = -X[1];
- xc = -X[0];
- //有任何一个维度的偏移超过0.5，会更新当前像素点
- //如果每一个维度的偏移都没有超过0.5,当前像素的位置加上偏移就是最终的精确点
- if( std::abs(xi) < 0.5f && std::abs(xr) < 0.5f && std::abs(xc) < 0.5f )
- break;
- if( std::abs(xi) > (float)(INT_MAX/3) ||
- std::abs(xr) > (float)(INT_MAX/3) ||
- std::abs(xc) > (float)(INT_MAX/3) )
- return false;
- c += cvRound(xc);
- r += cvRound(xr);
- layer += cvRound(xi);
- if( layer < 1 || layer > nOctaveLayers ||
- c < SIFT_IMG_BORDER || c >= img.cols - SIFT_IMG_BORDER ||
- r < SIFT_IMG_BORDER || r >= img.rows - SIFT_IMG_BORDER )
- return false;
- }
- //迭代结束
- // ensure convergence of interpolation
- if( i >= SIFT_MAX_INTERP_STEPS )
- return false;
- {
- int idx = octv*(nOctaveLayers+2) + layer;
- const Mat& img = dog_pyr[idx];
- const Mat& prev = dog_pyr[idx-1];
- const Mat& next = dog_pyr[idx+1];
- Matx31f dD((img.at<sift_wt>(r, c+1) - img.at<sift_wt>(r, c-1))*deriv_scale,
- (img.at<sift_wt>(r+1, c) - img.at<sift_wt>(r-1, c))*deriv_scale,
- (next.at<sift_wt>(r, c) - prev.at<sift_wt>(r, c))*deriv_scale);
- float t = dD.dot(Matx31f(xc, xr, xi));
- contr = img.at<sift_wt>(r, c)*img_scale + t * 0.5f;
- if( std::abs( contr ) * nOctaveLayers < contrastThreshold )//去除低对比度的点
- return false;
- // principal curvatures are computed using the trace and det of Hessian
- float v2 = img.at<sift_wt>(r, c)*2.f;
- float dxx = (img.at<sift_wt>(r, c+1) + img.at<sift_wt>(r, c-1) - v2)*second_deriv_scale;
- float dyy = (img.at<sift_wt>(r+1, c) + img.at<sift_wt>(r-1, c) - v2)*second_deriv_scale;
- float dxy = (img.at<sift_wt>(r+1, c+1) - img.at<sift_wt>(r+1, c-1) -
- img.at<sift_wt>(r-1, c+1) + img.at<sift_wt>(r-1, c-1)) * cross_deriv_scale;
- float tr = dxx + dyy;
- float det = dxx * dyy - dxy * dxy;
- if( det <= 0 || tr*tr*edgeThreshold >= (edgeThreshold + 1)*(edgeThreshold + 1)*det )//去除边缘噪声点
- return false;
- }
- kpt.pt.x = (c + xc) * (1 << octv);
- kpt.pt.y = (r + xr) * (1 << octv);
- kpt.octave = octv + (layer << 8) + (cvRound((xi + 0.5)*255) << 16);
- kpt.size = sigma*powf(2.f, (layer + xi) / nOctaveLayers)*(1 << octv)*2;
- kpt.response = std::abs(contr);
- return true;
- }
以上SIFT源码均摘自OpenCV nonfree模块，lowe对SIFT拥有版权。
符合要求的特征点构建完毕，需要对该特征点进行描述了，请关注本博客SIFT系列的下一篇文章：SIFT解析（三）生成特征描述子
