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::findMaximaInLayer( 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 = interpolateKeypoint( 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); //儲存該點
}
}
}
}
}
}
對于計算特征點位置精确計算函數如下 static int interpolateKeypoint( 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;
};
理論模型如下圖
