opencv安装目录/doc 下有个文件：opencv_tutorials.pdf 这个文件给我们提供很多学习代码。。有时间可以好好看下。。很有用的。。。我还没认真看。。。
现在进入正题
surf特征提取方法：
- std::vector<cv::Keypoint>keypoints;
- cv::SurfFeatureDetector surf(2500.); //阀值
- surf.detect(image,keypoints);
那么我们看下surf特征点是如何提取的？
在feature2d.hpp中有
typedef SURF SurfFeatureDetector对于cv::SurfFeatureDetector surf(2500.); -> // SURF类实例化对象surf 传入参数2500
SURF(double hessianThreshold, // = 2500
int nOctaves=4, int nOctaveLayers=2,
bool extended=true, bool upright=false); ->
SURF::SURF(double _threshold, int _nOctaves, int _nOctaveLayers, bool _extended, bool _upright)
{
hessianThreshold = _threshold;
extended = _extended;
upright = _upright;
nOctaves = _nOctaves;
nOctaveLayers = _nOctaveLayers;
}
对于surf.detect(image,keypoints);
SURF类中不存在detect方法，我们发现SURF类是继承Feature2D类，我们可以到Feature2D类中查找相关方法；
class CV_EXPORTS_W SURF : public Feature2D->
class CV_EXPORTS_W Feature2D : public FeatureDetector, public DescriptorExtractor ->
class CV_EXPORTS_W FeatureDetector : public virtual Algorithm
{
......
CV_WRAP void detect( const Mat& image, CV_OUT vector<KeyPoint>& keypoints, const Mat& mask=Mat() ) const;->
......
virtual void detectImpl( const Mat& image, vector<KeyPoint>& keypoints, const Mat& mask=Mat() ) const = 0;
}
detectImpl( image, keypoints, mask );
由于FeatureDector中的detectImpl为纯虚函数 则调用子类中的detectImpl 即调用到
void SURF::detectImpl( const Mat& image, vector<KeyPoint>& keypoints, const Mat& mask) const
(*this)(image, mask, keypoints, noArray(), false);->
void SURF::operator()(InputArray _img, InputArray _mask,
CV_OUT vector<KeyPoint>& keypoints,
OutputArray _descriptors,
bool useProvidedKeypoints) const->
fastHessianDetector( sum, msum, keypoints, nOctaves, nOctaveLayers, (float)hessianThreshold ); //这一步进行了角点提取
fastHessianDetector 开始创建尺度空间 接着调用
parallel_for( BlockedRange(0, nTotalLayers),SURFBuildInvoker(sum, sizes, sampleSteps, dets, traces) );->
body(range); ->
struct SURFBuildInvoker{
void operator()(const BlockedRange& range) const->
}
for( int i=range.begin(); i<range.end(); i++ )
calcLayerDetAndTrace( *sum, (*sizes)[i], (*sampleSteps)[i], (*dets)[i], (*traces)[i] ); //这里开始计算矩阵的行列式和迹了
(这个函数式通过对每个像素点（边界点除外）在它9*9领域内的盒子滤波 dxdx ，dydy ，dxdy ) 计算好后通过surf论文公式就能计算出每个像素点的行列式和迹。
parallel_for( BlockedRange(0, nMiddleLayers),
SURFFindInvoker(sum, mask_sum, dets, traces, sizes,
sampleSteps, middleIndices, keypoints,
nOctaveLayers, hessianThreshold) ); ->
body(range); ->
findMaximaInLayer( *sum, *mask_sum, *dets, *traces, *sizes,
*keypoints, octave, layer, hessianThreshold,
(*sampleSteps)[layer] ); ->
这个函数是读取每个像素的行列式的值，如果值大于threshold，则与该像素领域3*3*3的每个像素比较。。看是否为最大值。。如果为最大值则进行差值，精确极值点的位置
最后保存该极值点
对于计算特征点位置精确计算函数如下
- void SURFFindInvoker::<span style="color:#FF0000;">findMaximaInLayer</span>( const Mat& sum, const Mat& mask_sum,
- const vector<Mat>& dets, const vector<Mat>& traces,
- const vector<int>& sizes, vector<KeyPoint>& keypoints,
- int octave, int layer, float hessianThreshold, int sampleStep )
- {
- // Wavelet Data
- const int NM=1;
- const int dm[NM][5] = { {0, 0, 9, 9, 1} };
- SurfHF Dm;
- int size = sizes[layer];
- // The integral image 'sum' is one pixel bigger than the source image
- int layer_rows = (sum.rows-1)/sampleStep;
- int layer_cols = (sum.cols-1)/sampleStep;
- // Ignore pixels without a 3x3x3 neighbourhood in the layer above
- int margin = (sizes[layer+1]/2)/sampleStep+1;
- if( !mask_sum.empty() )
- resizeHaarPattern( dm, &Dm, NM, 9, size, mask_sum.cols );
- int step = (int)(dets[layer].step/dets[layer].elemSize());
- for( int i = margin; i < layer_rows - margin; i++ )
- {
- const float* det_ptr = dets[layer].ptr<float>(i);
- const float* trace_ptr = traces[layer].ptr<float>(i);
- for( int j = margin; j < layer_cols-margin; j++ )
- {
- float val0 = det_ptr[j];
- if( val0 > hessianThreshold ) //判断每个像素是否超过阀值
- {
- /* Coordinates for the start of the wavelet in the sum image. There
- is some integer division involved, so don't try to simplify this
- (cancel out sampleStep) without checking the result is the same */
- int sum_i = sampleStep*(i-(size/2)/sampleStep);
- int sum_j = sampleStep*(j-(size/2)/sampleStep);
- /* The 3x3x3 neighbouring samples around the maxima.
- The maxima is included at N9[1][4] */
- const float *det1 = &dets[layer-1].at<float>(i, j);
- const float *det2 = &dets[layer].at<float>(i, j);
- const float *det3 = &dets[layer+1].at<float>(i, j);
- float N9[3][9] = { { det1[-step-1], det1[-step], det1[-step+1], //以该点为中心的3层，每层9个像素
- det1[-1] , det1[0] , det1[1],
- det1[step-1] , det1[step] , det1[step+1] },
- { det2[-step-1], det2[-step], det2[-step+1],
- det2[-1] , det2[0] , det2[1],
- det2[step-1] , det2[step] , det2[step+1] },
- { det3[-step-1], det3[-step], det3[-step+1],
- det3[-1] , det3[0] , det3[1],
- det3[step-1] , det3[step] , det3[step+1] } };
- /* Check the mask - why not just check the mask at the center of the wavelet? */
- if( !mask_sum.empty() )
- {
- const int* mask_ptr = &mask_sum.at<int>(sum_i, sum_j);
- float mval = calcHaarPattern( mask_ptr, &Dm, 1 );
- if( mval < 0.5 )
- continue;
- }
- /* Non-maxima suppression. val0 is at N9[1][4]*/
- if( val0 > N9[0][0] && val0 > N9[0][1] && val0 > N9[0][2] && //计算val0是否为极大值点
- val0 > N9[0][3] && val0 > N9[0][4] && val0 > N9[0][5] &&
- val0 > N9[0][6] && val0 > N9[0][7] && val0 > N9[0][8] &&
- val0 > N9[1][0] && val0 > N9[1][1] && val0 > N9[1][2] &&
- val0 > N9[1][3] && val0 > N9[1][5] &&
- val0 > N9[1][6] && val0 > N9[1][7] && val0 > N9[1][8] &&
- val0 > N9[2][0] && val0 > N9[2][1] && val0 > N9[2][2] &&
- val0 > N9[2][3] && val0 > N9[2][4] && val0 > N9[2][5] &&
- val0 > N9[2][6] && val0 > N9[2][7] && val0 > N9[2][8] )
- {
- /* Calculate the wavelet center coordinates for the maxima */
- float center_i = sum_i + (size-1)*0.5f;
- float center_j = sum_j + (size-1)*0.5f;
- KeyPoint kpt( center_j, center_i, (float)sizes[layer],
- -1, val0, octave, CV_SIGN(trace_ptr[j]) );
- /* Interpolate maxima location within the 3x3x3 neighbourhood */
- int ds = size - sizes[layer-1];
- int interp_ok = <span style="color:#FF0000;">interpolateKeypoint</span>( N9, sampleStep, sampleStep, ds, kpt ); //若为极大值点还要进行插值计算求精
- /* Sometimes the interpolation step gives a negative size etc. */
- if( interp_ok )
- {
- /*printf( "KeyPoint %f %f %d\n", point.pt.x, point.pt.y, point.size );*/
- #ifdef HAVE_TBB
- tbb::mutex::scoped_lock lock(findMaximaInLayer_m);
- #endif
- keypoints.push_back(kpt); //保存该点
- }
- }
- }
- }
- }
- }
- <span style="color:#FF0000;">static int interpolateKeypoint</span>( float N9[3][9], int dx, int dy, int ds, KeyPoint& kpt )
- {
- Vec3f b(-(N9[1][5]-N9[1][3])/2, // Negative 1st deriv with respect to x //对该像素点求一阶导。。
- -(N9[1][7]-N9[1][1])/2, // Negative 1st deriv with respect to y
- -(N9[2][4]-N9[0][4])/2); // Negative 1st deriv with respect to s
- Matx33f A(
- N9[1][3]-2*N9[1][4]+N9[1][5], // 2nd deriv x, x //对该像素点求2阶导
- (N9[1][8]-N9[1][6]-N9[1][2]+N9[1][0])/4, // 2nd deriv x, y
- (N9[2][5]-N9[2][3]-N9[0][5]+N9[0][3])/4, // 2nd deriv x, s
- (N9[1][8]-N9[1][6]-N9[1][2]+N9[1][0])/4, // 2nd deriv x, y
- N9[1][1]-2*N9[1][4]+N9[1][7], // 2nd deriv y, y
- (N9[2][7]-N9[2][1]-N9[0][7]+N9[0][1])/4, // 2nd deriv y, s
- (N9[2][5]-N9[2][3]-N9[0][5]+N9[0][3])/4, // 2nd deriv x, s
- (N9[2][7]-N9[2][1]-N9[0][7]+N9[0][1])/4, // 2nd deriv y, s
- N9[0][4]-2*N9[1][4]+N9[2][4]); // 2nd deriv s, s
- Vec3f x = A.solve(b, DECOMP_LU); //计算极值点的精确位置
- bool ok = (x[0] != 0 || x[1] != 0 || x[2] != 0) &&
- std::abs(x[0]) <= 1 && std::abs(x[1]) <= 1 && std::abs(x[2]) <= 1;
- if( ok )
- {
- kpt.pt.x += x[0]*dx;
- kpt.pt.y += x[1]*dy;
- kpt.size = (float)cvRound( kpt.size + x[2]*ds );
- }
- return ok;
- };
理论模型如下图
扫一扫
2181
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
