OpenCV中的SURF类
opencv中SURF被实现为一个类,包含的几个主要成员变量如下:
hessianThreshold 是特征点的判别门限
nOctaves 是图像金字塔的组数
nOctaveLayers 图像金字塔每个组内有多少层
extended为true时,特征向量是128维;否则特征向量是64维
upright为true时,所有特征点的主方向都被认为是 -90°;否则特征点的主方向通过计算得出
- class CV_EXPORTS_W SURF : public Feature2D
- {
- public:
- ...
- ...
- CV_PROP_RW double hessianThreshold;
- CV_PROP_RW int nOctaves;
- CV_PROP_RW int nOctaveLayers;
- CV_PROP_RW bool extended;
- CV_PROP_RW bool upright;
- explicit CV_WRAP SURF(double hessianThreshold,
- int nOctaves=4, int nOctaveLayers=2,
- bool extended=true, bool upright=false);
- ...
- protected:
- ...
- };
通过重载SURF类的 "()" 运算符,来实现具体的SURF操作,包括SURF特征点检测和特征点描述符的计算。
- void SURF::operator()(InputArray _img, InputArray _mask,
- CV_OUT vector<KeyPoint>& keypoints,
- OutputArray _descriptors,
- bool useProvidedKeypoints) const
- {
- Mat img = _img.getMat(), mask = _mask.getMat(), mask1, sum, msum;
- bool doDescriptors = _descriptors.needed();
- CV_Assert(!img.empty() && img.depth() == CV_8U);
- if( img.channels() > 1 )
- cvtColor(img, img, COLOR_BGR2GRAY);
- CV_Assert(mask.empty() || (mask.type() == CV_8U && mask.size() == img.size()));
- CV_Assert(hessianThreshold >= 0);
- CV_Assert(nOctaves > 0);
- CV_Assert(nOctaveLayers > 0);
- integral(img, sum, CV_32S);
- 先作一些输入合法性检查**注意输入的图像img只支持BGR 8位图像和单通道8位图像**,然后计算积分图像(存入 sum 中)。
- // Compute keypoints only if we are not asked for evaluating the descriptors are some given locations:
- if( !useProvidedKeypoints )
- {
- if( !mask.empty() )
- {
- cv::min(mask, 1, mask1);
- integral(mask1, msum, CV_32S);
- }
- fastHessianDetector( sum, msum, keypoints, nOctaves, nOctaveLayers, (float)hessianThreshold );
- }
- 如果没有指定特征点,则用fastHessianDetector()来检测特征点。如果指定了mask掩码图像,则首先将掩码图像二值化成(0,1)存入mask1,再计算mask1的积分图像msum,fastHessianDetector()函数中会用到msum。(需要使用msum的具体位置是在SURFFindInvoker::findMaximaInLayer()中,来屏蔽某些位置的特征点)
- int i, j, N = (int)keypoints.size();
- if( N > 0 )
- {
- Mat descriptors;
- bool _1d = false;
- int dcols = extended ? 128 : 64;
- size_t dsize = dcols*sizeof(float);
- dcols代表每个特征向量有多少维,dsize代表存储每个特征向量需要多少字节的空间。
- 还有一个bool型变量 _1d,它指明了输出的描述符数组是否是单列的:如果 _1d 为true,说明输出的描述数组是个 1 维单列数组(每行一个元素),每dcols行个元素构成一个特征向量;如果 _1d 为true,说明输出的描述符数组是个 2 维数组,每行有 dcols列,每行代表一个特征向量。
- if( doDescriptors )
- {
- _1d = _descriptors.kind() == _InputArray::STD_VECTOR && _descriptors.type() == CV_32F;
- if( _1d )
- {
- _descriptors.create(N*dcols, 1, CV_32F);
- descriptors = _descriptors.getMat().reshape(1, N);
- }
- else
- {
- _descriptors.create(N, dcols, CV_32F);
- descriptors = _descriptors.getMat();
- }
- }
- 如果需要计算特征点的描述符,则创建一个 N 行、dcols列 的矩阵descriptors,用于存储描述符数组。
- // we call SURFInvoker in any case, even if we do not need descriptors,
- // since it computes orientation of each feature.
- parallel_for_(Range(0, N), SURFInvoker(img, sum, keypoints, descriptors, extended, upright) );
- 无论是否需要计算描述符,都应调用SURFInvoker,因为他可以计算出特征点的主方向。
这里打断一下,前面出现了一个函数 parallel_for_(),它是什么作用呢?
如果有多组数据需要进行相似的操作,且这些操作互不影响,则这些数据可以被并行处理,在多核处理器上这可以大大增加算法的执行速度。而函数parallel_for_()的作用就是对某个区间中的多组数据进行某种相同的操作,且最大可能地将这些操作并行。它的原型为:
CV_EXPORTS void parallel_for_(const Range& range, const ParallelLoopBody& body, double nstripes=-1.);
第一个参数range代表数据的区间;第二个参数类型ParallelLoopBody,它是一个抽象接口类,封装了某种"操作",定义如下,(只重载了一个 () 操作符,且参数为 Range)
class CV_EXPORTS ParallelLoopBody
{
public:
virtual ~ParallelLoopBody();
virtual void operator() (const Range& range) const = 0;
};
对于不存在多线程机制的情况下,parallel_for_( range, body ) 等价于直接调用 body(range) .
而上面出现的 SURFInvoker 则是 ParallelLoopBody 的一个派生类。因此仅从功能上来看的话,上面那行代码的作用相当于:
SURFInvoker(img, sum, keypoints, descriptors, extended, upright).(Range(0, N)) ;
不过这样写就不能利用并行特性,运算速度会降低。
SURFInvoker主要用于提取特征点的主方向、描述符等信息。
下面继续回到 SURF::operator() 函数中,接之前的代码。
- // remove keypoints that were marked for deletion
- for( i = j = 0; i < N; i++ )
- {
- if( keypoints[i].size > 0 )
- {
- if( i > j )
- {
- keypoints[j] = keypoints[i];
- if( doDescriptors )
- memcpy( descriptors.ptr(j), descriptors.ptr(i), dsize);
- }
- j++;
- }
- }
- if( N > j )
- {
- N = j;
- keypoints.resize(N);
- if( doDescriptors )
- {
- Mat d = descriptors.rowRange(0, N);
- if( _1d )
- d = d.reshape(1, N*dcols);
- d.copyTo(_descriptors);
- }
- }
- 做一些收尾工作: 删掉所有 size <= 0 的特征点(特征点的 size<=0 表示该特征点被标记为待删除),并更新特征点的数量,对输出描述符数组的维数进行转换等。
- }
- }
SRUF特征点的检测: fastHessianDetector()
SURF特征点检测。先求每一层图像的Hessian行列式(迹),再在尺度空间中寻找极值点作为特征点。
- static void fastHessianDetector( const Mat& sum, const Mat& mask_sum, vector<KeyPoint>& keypoints,
- int nOctaves, int nOctaveLayers, float hessianThreshold )
- {
- /* Sampling step along image x and y axes at first octave. This is doubled
- for each additional octave. WARNING: Increasing this improves speed,
- however keypoint extraction becomes unreliable. */
- const int SAMPLE_STEP0 = 1;
- int nTotalLayers = (nOctaveLayers+2)*nOctaves;
- int nMiddleLayers = nOctaveLayers*nOctaves;
- 每个Octave中的图像层数应比指定的nOctaveLayers多2层,因为在尺度空间寻找极值点时,每个Octave需要额外的两层。
- nTotalLayers代表所有层的数量,包括额外层;
- nMiddleLayers代表所有“中间层”的数量,即不包含额外层;
- vector<Mat> dets(nTotalLayers);
- vector<Mat> traces(nTotalLayers);
- vector<int> sizes(nTotalLayers);
- vector<int> sampleSteps(nTotalLayers);
- vector<int> middleIndices(nMiddleLayers);
- 定义了一系列vector:
- dets : 记录每层图像的 Hessian 行列式;
- traces : 记录每层图像的 Hessian 迹;SIFT特征点算法中为了抑制边缘区域的响应所以需要计算Hessian 迹,但SURF特征点算法本身就不会对边缘区域有响应,为什么还需要计算迹?只是为了通过迹的符号(正负)来判断特征点的凹凸特性,具体参见后面SURFFindInvoker::findMaximaInLayer()函数中的相关注释。
- sizes: 记录每层图像的方框size
- sampleSteps: 记录每层图像的采样间隔
- middleIndices: 记录所有中间层图像的全局层索引,通过这个数组可以知道第 k 个中间层对应的全局层号是多少。注意前四个vector的size都是 nTotalLayers,只有这个vector的size是nMiddleLayers。
- keypoints.clear();
- // Allocate space and calculate properties of each layer
- int index = 0, middleIndex = 0, step = SAMPLE_STEP0;
- for( int octave = 0; octave < nOctaves; octave++ )
- {
- for( int layer = 0; layer < nOctaveLayers+2; layer++ )
- {
- /* The integral image sum is one pixel bigger than the source image*/
- dets[index].create( (sum.rows-1)/step, (sum.cols-1)/step, CV_32F );
- traces[index].create( (sum.rows-1)/step, (sum.cols-1)/step, CV_32F );
- sizes[index] = (SURF_HAAR_SIZE0 + SURF_HAAR_SIZE_INC*layer) << octave;
- sampleSteps[index] = step;
- if( 0 < layer && layer <= nOctaveLayers )
- middleIndices[middleIndex++] = index;
- index++;
- }
- step *= 2;
- }
- 为dets和traces分配空间,并初始化 sizes、sampleSteps 和 middleIndices 。
- 上面的step 被初始化为 SAMPLE_STEP0,后者被定义为1。且step在每次octave循环后会乘以2, "step *= 2",因此第octave组中的各层图像的采样间隔都是2^(octave);
- SURF_HAAR_SIZE0被定义为9,代表最底层图像的方框尺寸(对应尺度1.2);SURF_HAAR_SIZE_INC被定义为6,代表第0组图像中相邻两层图像的方框尺寸差,而第octave组图像中相邻两层图像的方框尺寸差则为 6*2^(octave),因此,第octave组第layer层的图像的方框尺寸size=(9+6*layer)*2^(octave)。
- middleIndices记录了所有中间层的全局层索引
- // Calculate hessian determinant and trace samples in each layer
- parallel_for_( Range(0, nTotalLayers),
- SURFBuildInvoker(sum, sizes, sampleSteps, dets, traces) );
- 利用SURFBuildInvoker为每一层图像计算行列式和迹。注意这个操作的Range是从0到nTotalLayers。
- // Find maxima in the determinant of the hessian
- parallel_for_( Range(0, nMiddleLayers),
- SURFFindInvoker(sum, mask_sum, dets, traces, sizes,
- sampleSteps, middleIndices, keypoints,
- nOctaveLayers, hessianThreshold) );
- 利用SURFFindInvoker寻找极值。注意这个操作的Range是从0到nMiddleLayers,只在中间层中寻找极值。
- std::sort(keypoints.begin(), keypoints.end(), KeypointGreater());
- 排序
- }
到这里,只剩下3个具体的操作需要研究了,这三个操作也都被实现成了三个类,他们分别是:
SURFBuildInvoker,用于计算每一层图像的hessian行列式和迹
SURFFindInvoker,用于在尺度空间寻找极值点(特征点)
SURFInvoker,用于计算特征点的主方向和描述符
下面分别介绍。从SURFInvoker(特征点描述符的计算)开始,这个最麻烦。另外两个操作都比较简单。
SURFInvoker: 特征点主方向和描述符的计算
SURFInvoker类的定义如下。我们主要关心它的构造函数和重载的()运算符。SURFInvoker 是代码最繁琐的类。
- struct SURFInvoker : ParallelLoopBody
- {
- enum { ORI_RADIUS = 6, ORI_WIN = 60, PATCH_SZ = 20 };
- 这里首先定义了几个在SURF算法中需要使用的一些常量参数,它们被定义为枚举变量。
- 先解释一下这些枚举变量的含义:
- ORI_RADIUS : 计算关键点的主方向时,邻域(圆形)区域的半径。被定义为6,且单位是关键点的尺度
- ORI_WIN : 计算关键点的主方向时,扇形窗的夹角。被定义为60°
- PATCH_SZ : 计算关键点的描述符时,邻域(方形)区域的直径。被定义为20,且单位是关键点的尺度
- // Parameters
- const Mat* img;
- const Mat* sum;
- vector<KeyPoint>* keypoints;
- Mat* descriptors;
- bool extended;
- bool upright;
- // Pre-calculated values
- int nOriSamples;
- vector<Point> apt;
- vector<float> aptw;
- vector<float> DW;
- 这里列出的是 SURFInvoker 的成员变量。可分为两部分,一部分是传入的参数,另一部分是一些需要预先计算出来的变量或数组。
- nOriSamples : 用于计算关键点主方向的邻域采样点数量。当关键点尺度 σ 较大时,位于邻域内的样点有很多,但只有坐标偏移是整数倍 σ 的样点才被用于计算主方向。因此不用担心计算量随 σ 增加而增加。为叙述方便,后面将“邻域内坐标偏移是整数倍 σ 的样点”称为“邻域网格样点”。
- apt : 存储了用于计算关键点主方向的所有邻域网格样点的坐标偏移,单位 σ 。
- aptw: 存储了用于计算关键点主方向的所有邻域网格样点的权重。高斯SIGMA = 2.5
- DW: 存储了用于计算关键点描述符的邻域采样点的权重。高斯SIGMA = 3.3
- SURFInvoker的构造函数。它主要初始化了自己所用到的参数。包括需要预先计算的变量和数组,即nOriSamples、apt、aptw 和 DW。
- SURFInvoker( const Mat& _img, const Mat& _sum,
- vector<KeyPoint>& _keypoints, Mat& _descriptors,
- bool _extended, bool _upright )
- {
- keypoints = &_keypoints;
- descriptors = &_descriptors;
- img = &_img;
- sum = &_sum;
- extended = _extended;
- upright = _upright;
- // Simple bound for number of grid points in circle of radius ORI_RADIUS
- const int nOriSampleBound = (2*ORI_RADIUS+1)*(2*ORI_RADIUS+1);
- // Allocate arrays
- apt.resize(nOriSampleBound);
- aptw.resize(nOriSampleBound);
- DW.resize(PATCH_SZ*PATCH_SZ);
- /* Coordinates and weights of samples used to calculate orientation */
- Mat G_ori = getGaussianKernel( 2*ORI_RADIUS+1, SURF_ORI_SIGMA, CV_32F );
- nOriSamples = 0;
- for( int i = -ORI_RADIUS; i <= ORI_RADIUS; i++ )
- {
- for( int j = -ORI_RADIUS; j <= ORI_RADIUS; j++ )
- {
- if( i*i + j*j <= ORI_RADIUS*ORI_RADIUS )
- {
- apt[nOriSamples] = cvPoint(i,j);
- aptw[nOriSamples++] = G_ori.at<float>(i+ORI_RADIUS,0) * G_ori.at<float>(j+ORI_RADIUS,0);
- }
- }
- }
- CV_Assert( nOriSamples <= nOriSampleBound );
- /* Gaussian used to weight descriptor samples */
- Mat G_desc = getGaussianKernel( PATCH_SZ, SURF_DESC_SIGMA, CV_32F );
- for( int i = 0; i < PATCH_SZ; i++ )
- {
- for( int j = 0; j < PATCH_SZ; j++ )
- DW[i*PATCH_SZ+j] = G_desc.at<float>(i,0) * G_desc.at<float>(j,0);
- }
- }
- 重载的()运算符实现特征点主方向和描述符的计算。
- void operator()(const Range& range) const
- {
- /* X and Y gradient wavelet data */
- const int NX=2, NY=2;
- const int dx_s[NX][5] = {{0, 0, 2, 4, -1}, {2, 0, 4, 4, 1}};
- const int dy_s[NY][5] = {{0, 0, 4, 2, 1}, {0, 2, 4, 4, -1}};
- dx_s和dy_s用于生成计算X和Y方向Haar特征时需要的特征矩形。参见后面介绍的SurfHF结构体
- // Optimisation is better using nOriSampleBound than nOriSamples for
- // array lengths. Maybe because it is a constant known at compile time
- const int nOriSampleBound =(2*ORI_RADIUS+1)*(2*ORI_RADIUS+1);
- float X[nOriSampleBound], Y[nOriSampleBound], angle[nOriSampleBound];
- uchar PATCH[PATCH_SZ+1][PATCH_SZ+1];
- float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ];
- CvMat matX = cvMat(1, nOriSampleBound, CV_32F, X);
- CvMat matY = cvMat(1, nOriSampleBound, CV_32F, Y);
- CvMat _angle = cvMat(1, nOriSampleBound, CV_32F, angle);
- Mat _patch(PATCH_SZ+1, PATCH_SZ+1, CV_8U, PATCH);
- matX、matY、_angle三个矩阵分别是所有邻域网格样点的X小波、Y小波和方向。
- _patch是计算关键点描述符时,用到的矩阵
- int dsize = extended ? 128 : 64;
- int k, k1 = range.start, k2 = range.end;
- float maxSize = 0;
- for( k = k1; k < k2; k++ )
- {
- maxSize = std::max(maxSize, (*keypoints)[k].size);
- }
- int imaxSize = std::max(cvCeil((PATCH_SZ+1)*maxSize*1.2f/9.0f), 1);
- Ptr<CvMat> winbuf = cvCreateMat( 1, imaxSize*imaxSize, CV_8U );
- [k1,k2)是需要被计算的关键点的索引区间,即只有第 k1到第(k2-1) 之间的关键点需要被计算。
- 这里先计算出了在所有需要被计算的关键点中最大的size,然后反推出了其尺度 σ :
- σ = size / 9 * 1.2
- 这个公式的依据是,size 为 9 的关键点对应的尺度 σ 是0.9,因此按照这个比例关系可求出其他size对应的尺度。
- 然后,构造了一个矩阵winbuf,它是一个 imaxSize*imaxSize 的矩阵,其中 imaxSize = (PATCH_SZ+1) * σ 。
- 显然这个矩阵是用来计算描述符的。
- 下面开始循环计算每个关键点的主方向和描述符。
- for( k = k1; k < k2; k++ )
- {
- int i, j, kk, nangle;
- float* vec;
- SurfHF dx_t[NX], dy_t[NY];
- SurfHF的定义如下:
- struct SurfHF
- {
- int p0, p1, p2, p3;
- float w;
- SurfHF(): p0(0), p1(0), p2(0), p3(0), w(0) {}
- };
- 利用积分图像计算Haar特征时,使用这个结构体可以方便地描述一个矩形区域,再结合积分图像就能很方便地求出这个矩形区域得像素和。X-haar特征等于左半边矩形减去右半边矩形,这两个矩形分别被两个SurfHF结构体表示,且其中一个具有正的权重w,另一个具有负的权重。
- p0、p1、p2、p3 分别对应矩形区域的 左上角、右上角、左下角、右下角,权重 w 等于矩形区域面积(包含像素点数量)的倒数,w的正负号取决于这个矩形是白的还是黑的。
- 每个Haar特征的计算都需要两个SurfHF结构,所以NX、NY都是2。其中 dx_t用于计算X方向的Haar特征,dy_t用于计算Y方向的Haar特征
- 利用积分图像求一个矩形区域的像素和时,可利用如下等式:
- 矩形内像素和 = 左上角积分 + 右下角积分 - 左下角积分 - 右上角积分
- KeyPoint& kp = (*keypoints)[k];
- float size = kp.size;
- Point2f center = kp.pt;
- /* The sampling intervals and wavelet sized for selecting an orientation
- and building the keypoint descriptor are defined relative to 's' */
- float s = size*1.2f/9.0f;
- 这里的 s 就是上面提到的 σ ,代表关键点的尺度:
- σ = size / 9 * 1.2
- /* To find the dominant orientation, the gradients in x and y are
- sampled in a circle of radius 6s using wavelets of size 4s.
- We ensure the gradient wavelet size is even to ensure the
- wavelet pattern is balanced and symmetric around its center */
- int grad_wav_size = 2*cvRound( 2*s );
- if( sum->rows < grad_wav_size || sum->cols < grad_wav_size )
- {
- /* when grad_wav_size is too big,
- * the sampling of gradient will be meaningless
- * mark keypoint for deletion. */
- kp.size = -1;
- continue;
- }
- 为了找到主方向,需要计算以关键点为中心、半径为6s(即6σ)的圆内的所有邻域网格样点的x、y方向的梯度小波,计算梯度小波所用的方框直径为grad_wav_size = 4s。
- 如果梯度小波的方框直径过大(比图像的尺寸还打),则认为该关键点为无效点,将其size设为-1,这样在SURF的后续处理中就会将其删掉。
- float descriptor_dir = 360.f - 90.f;
- if (upright == 0)
- {
- resizeHaarPattern( dx_s, dx_t, NX, 4, grad_wav_size, sum->cols );
- resizeHaarPattern( dy_s, dy_t, NY, 4, grad_wav_size, sum->cols );
- 这里用两个resizeHaarPattern函数分别计算出了 dx_t[NX] 和 dy_t[NY],这是计算Haar特征时需要用到的矩形 。其中NX、NY都是2。
- resizeHaarPattern 用于Haar模板的变形: 例如我们有了某个尺度下 Haar模板的特征矩形,就可以利用这个函数得到在一个新的尺度下(放缩后的)Haar模板的特征矩形。
- 下面循环计算每个邻域网格样点对主方向的贡献
- for( kk = 0, nangle = 0; kk < nOriSamples; kk++ )
- {
- int x = cvRound( center.x + apt[kk].x*s - (float)(grad_wav_size-1)/2 );
- int y = cvRound( center.y + apt[kk].y*s - (float)(grad_wav_size-1)/2 );
- 计算以第 kk 个样点为中心的梯度小波窗的左上角“最近邻取整”后的坐标 (x,y) ,其中center是关键点的坐标,s是关键点的尺度
- if( y < 0 || y >= sum->rows - grad_wav_size ||
- x < 0 || x >= sum->cols - grad_wav_size )
- continue;
- 如果梯度小波窗有一部分在图像的边界以外,则忽略该样点。
- const int* ptr = &sum->at<int>(y, x);
- float vx = calcHaarPattern( ptr, dx_t, 2 );
- float vy = calcHaarPattern( ptr, dy_t, 2 );
- X[nangle] = vx*aptw[kk];
- Y[nangle] = vy*aptw[kk];
- nangle++;
- 利用积分图像和 dx_t 、dy_t 计算出X、Y两个方向的Haar特征。
- nangle记录了关键点周围的有效样点的数量,有效样点表示该样点的梯度小波窗整个都在图像的边界之内。
- }
- if( nangle == 0 )
- {
- // No gradient could be sampled because the keypoint is too
- // near too one or more of the sides of the image. As we
- // therefore cannot find a dominant direction, we skip this
- // keypoint and mark it for later deletion from the sequence.
- kp.size = -1;
- continue;
- }
- 如果nangle是0,说明关键点(及其邻域的样点)太靠近图像的边界,这时我们标记这个关键点(size设为-1),以便后期删除之。
- matX.cols = matY.cols = _angle.cols = nangle;
- cvCartToPolar( &matX, &matY, 0, &_angle, 1 );
- 根据梯度小波矩阵计算所有邻域网格样点的方向
- float bestx = 0, besty = 0, descriptor_mod = 0;
- for( i = 0; i < 360; i += SURF_ORI_SEARCH_INC )
- {
- 计算60°的角度窗口内包含的所有样点的“组合梯度小波”:sumx和sumy,
- float sumx = 0, sumy = 0, temp_mod;
- for( j = 0; j < nangle; j++ )
- {
- int d = std::abs(cvRound(angle[j]) - i);
- if( d < ORI_WIN/2 || d > 360-ORI_WIN/2 )
- {
- sumx += X[j];
- sumy += Y[j];
- }
- }
- temp_mod = sumx*sumx + sumy*sumy;
- temp_mod代表关键点在当前角度方向上的响应,descriptor_mod代表关键点在其主方向附近的响应。在主方向附近,关键点的响应最大。下面的判断用于更新descriptor_mod以及bestx和besty: 响应越大的方向,越接近主方向。
- if( temp_mod > descriptor_mod )
- {
- descriptor_mod = temp_mod;
- bestx = sumx;
- besty = sumy;
- }
- }
- descriptor_dir = fastAtan2( -besty, bestx );
- 最后用bestx和besty计算主方向descriptor_dir。注意,关键点的主方向(角度)是按照笛卡尔坐标的习惯标定正负号,而笛卡尔坐标系与图像坐标系的y轴是相反的,因此上面计算atan2时用的是负的besty。
- }
- kp.angle = descriptor_dir;
- 将主方向descriptor_dir存储到kp.angle中。
- if( !descriptors || !descriptors->data )
- continue;
- 如果不需要计算描述符,则直接开始下一个关键点;否则,运行后面的代码来计算描述符
- /* Extract a window of pixels around the keypoint of size 20s */
- int win_size = (int)((PATCH_SZ+1)*s);
- CV_Assert( winbuf->cols >= win_size*win_size );
- Mat win(win_size, win_size, CV_8U, winbuf->data.ptr);
- 矩阵win是计算描述符用的窗口。这个矩形窗口是“斜”的,它的倾斜方向与关键点的主方向一致
- 矩阵win的size是 (PATCH_SZ+1)*s,为什么要加1?因为后面求Haar特征时使用的是 2s*2s 的窗,为了最后生成 PATCH_SZ*PATCH_SZ 即 20*20 个Haar特征,必须额外加1
- if( !upright )
- {
- descriptor_dir *= (float)(CV_PI/180);
- float sin_dir = -std::sin(descriptor_dir);
- float cos_dir = std::cos(descriptor_dir);
- /* Subpixel interpolation version (slower). Subpixel not required since
- the pixels will all get averaged when we scale down to 20 pixels */
- /*
- float w[] = { cos_dir, sin_dir, center.x,
- -sin_dir, cos_dir , center.y };
- CvMat W = cvMat(2, 3, CV_32F, w);
- cvGetQuadrangleSubPix( img, &win, &W );
- */
- float win_offset = -(float)(win_size-1)/2;
- float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;
- float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;
- uchar* WIN = win.data;
- #if 0
- // Nearest neighbour version (faster)
- for( i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir )
- {
- float pixel_x = start_x;
- float pixel_y = start_y;
- for( j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir )
- {
- int x = std::min(std::max(cvRound(pixel_x), 0), img->cols-1);
- int y = std::min(std::max(cvRound(pixel_y), 0), img->rows-1);
- WIN[i*win_size + j] = img->at<uchar>(y, x);
- }
- }
- #else
- int ncols1 = img->cols-1, nrows1 = img->rows-1;
- size_t imgstep = img->step;
- for( i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir )
- {
- double pixel_x = start_x;
- double pixel_y = start_y;
- for( j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir )
- {
- int ix = cvFloor(pixel_x), iy = cvFloor(pixel_y);
- if( (unsigned)ix < (unsigned)ncols1 &&
- (unsigned)iy < (unsigned)nrows1 )
- {
- float a = (float)(pixel_x - ix), b = (float)(pixel_y - iy);
- const uchar* imgptr = &img->at<uchar>(iy, ix);
- WIN[i*win_size + j] = (uchar)
- cvRound(imgptr[0]*(1.f - a)*(1.f - b) +
- imgptr[1]*a*(1.f - b) +
- imgptr[imgstep]*(1.f - a)*b +
- imgptr[imgstep+1]*a*b);
- }
- else
- {
- int x = std::min(std::max(cvRound(pixel_x), 0), ncols1);
- int y = std::min(std::max(cvRound(pixel_y), 0), nrows1);
- WIN[i*win_size + j] = img->at<uchar>(y, x);
- }
- }
- }
- 用双线性插值法计算“倾斜窗”win的所有像素。图像边界上无法进行双线性插值时,使用最近邻插值。
- #endif
- }
- else
- {
- // extract rect - slightly optimized version of the code above
- // TODO: find faster code, as this is simply an extract rect operation,
- // e.g. by using cvGetSubRect, problem is the border processing
- // descriptor_dir == 90 grad
- // sin_dir == 1
- // cos_dir == 0
- float win_offset = -(float)(win_size-1)/2;
- int start_x = cvRound(center.x + win_offset);
- int start_y = cvRound(center.y - win_offset);
- uchar* WIN = win.data;
- for( i = 0; i < win_size; i++, start_x++ )
- {
- int pixel_x = start_x;
- int pixel_y = start_y;
- for( j = 0; j < win_size; j++, pixel_y-- )
- {
- int x = MAX( pixel_x, 0 );
- int y = MAX( pixel_y, 0 );
- x = MIN( x, img->cols-1 );
- y = MIN( y, img->rows-1 );
- WIN[i*win_size + j] = img->at<uchar>(y, x);
- }
- }
- 主方向指定为-90°时,即(upright==true)时,直接使用最近邻插值获取倾斜窗的像素
- }
- // Scale the window to size PATCH_SZ so each pixel's size is s. This
- // makes calculating the gradients with wavelets of size 2s easy
- resize(win, _patch, _patch.size(), 0, 0, INTER_AREA);
- 调整(缩小)倾斜窗win的尺寸以便计算Haar特征(图像win缩小后,计算Haar特征使用的 2s*2s 的窗,就变成了 2*2 的窗)。
- 缩小图像时插值方法使用区域插值INTER_AREA,即计算被覆盖区域像素的平均值。
- 缩小后的图像大小为 (PATCH_SZ+1)*(PATCH_SZ+1),即 21*21。
- 后面计算Haar特征时使用的都是 2*2 的窗,因此 21*21 的图像在X、Y方向可以各生成 20*20 个Haar特征。
- // Calculate gradients in x and y with wavelets of size 2s
- for( i = 0; i < PATCH_SZ; i++ )
- for( j = 0; j < PATCH_SZ; j++ )
- {
- float dw = DW[i*PATCH_SZ + j];
- float vx = (PATCH[i][j+1] - PATCH[i][j] + PATCH[i+1][j+1] - PATCH[i+1][j])*dw;
- float vy = (PATCH[i+1][j] - PATCH[i][j] + PATCH[i+1][j+1] - PATCH[i][j+1])*dw;
- DX[i][j] = vx;
- DY[i][j] = vy;
- }
- 共能求出 PATCH_SZ * PATCH_SZ 即 20*20 个X-Haar特征和Y-Haar特征。分别存入DX和DY数组(已进行高斯加权)。
- // Construct the descriptor
- vec = descriptors->ptr<float>(k);
- for( kk = 0; kk < dsize; kk++ )
- vec[kk] = 0;
- double square_mag = 0;
- if( extended )
- {
- // 128-bin descriptor
- for( i = 0; i < 4; i++ )
- for( j = 0; j < 4; j++ )
- {
- for(int y = i*5; y < i*5+5; y++ )
- {
- for(int x = j*5; x < j*5+5; x++ )
- {
- float tx = DX[y][x], ty = DY[y][x];
- if( ty >= 0 )
- {
- vec[0] += tx;
- vec[1] += (float)fabs(tx);
- } else {
- vec[2] += tx;
- vec[3] += (float)fabs(tx);
- }
- if ( tx >= 0 )
- {
- vec[4] += ty;
- vec[5] += (float)fabs(ty);
- } else {
- vec[6] += ty;
- vec[7] += (float)fabs(ty);
- }
- }
- }
- for( kk = 0; kk < 8; kk++ )
- square_mag += vec[kk]*vec[kk];
- vec += 8;
- }
- }
- else
- {
- // 64-bin descriptor
- for( i = 0; i < 4; i++ )
- for( j = 0; j < 4; j++ )
- {
- for(int y = i*5; y < i*5+5; y++ )
- {
- for(int x = j*5; x < j*5+5; x++ )
- {
- float tx = DX[y][x], ty = DY[y][x];
- vec[0] += tx; vec[1] += ty;
- vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty);
- }
- }
- for( kk = 0; kk < 4; kk++ )
- square_mag += vec[kk]*vec[kk];
- vec+=4;
- }
- }
- 根据计算出的Haar特征来组成描述符。
- 首先将 20*20 个Haar特征等分成 4*4=16 个子区域,每个子区域有 5*5 个Haar特征。对每个子区域的Haar特征按某种方式求和,得到4或8个特征。这样每个关键点最后总共有 4*4*4 或 4*4*8 个特征。
- 此外还计算了 square_mag,它代表所有特征的平方和,用于归一化。下面几行就是进行特征向量的归一化操作
- // unit vector is essential for contrast invariance
- vec = descriptors->ptr<float>(k);
- float scale = (float)(1./(sqrt(square_mag) + DBL_EPSILON));
- for( kk = 0; kk < dsize; kk++ )
- vec[kk] *= scale;
- }
- }
- };
- 计算主方向与描述符的对比:
- 计算主方向时,搜索窗是半径为6s的圆,且计算Haar特征时使用的是 尺寸为 4s*4s 的正窗(不倾斜),由于是正窗,因此可借助积分图像快速求取Haar特征;
- 计算描述符时,搜索窗是边长为21s的方形,且是倾斜的,计算Haar特征时使用的也是 尺寸为 2s*2s 的倾斜窗,由于是倾斜窗,因此不能利用积分图像求Haar特征。
SURFBuildInvoker:行列式和迹的计算
再来看Hessian矩阵的行列式和迹是如何被计算的。
SURFBuildInvoker很简单,它重载的()运算符只是通过调用calcLayerDetAndTrace(),利用积分图像sum来求取图像的dets矩阵和traces矩阵。
calcLayerDetAndTrace() 函数的内部也是使用了前面提到的 SurfHF 结构体 (表示一个矩形区域)来计算二阶导数,二阶导数被近似成 “类Haar特征”。
- // Multi-threaded construction of the scale-space pyramid
- struct SURFBuildInvoker : ParallelLoopBody
- {
- SURFBuildInvoker( const Mat& _sum, const vector<int>& _sizes,
- const vector<int>& _sampleSteps,
- vector<Mat>& _dets, vector<Mat>& _traces )
- {
- sum = &_sum;
- sizes = &_sizes;
- sampleSteps = &_sampleSteps;
- dets = &_dets;
- traces = &_traces;
- }
- void operator()(const Range& range) const
- {
- for( int i=range.start; i<range.end; i++ )
- calcLayerDetAndTrace( *sum, (*sizes)[i], (*sampleSteps)[i], (*dets)[i], (*traces)[i] );
- }
- const Mat *sum;
- const vector<int> *sizes;
- const vector<int> *sampleSteps;
- vector<Mat>* dets;
- vector<Mat>* traces;
- };
- /*
- * Calculate the determinant and trace of the Hessian for a layer of the
- * scale-space pyramid
- */
- static void calcLayerDetAndTrace( const Mat& sum, int size, int sampleStep,
- Mat& det, Mat& trace )
- {
- const int NX=3, NY=3, NXY=4;
- const int dx_s[NX][5] = { {0, 2, 3, 7, 1}, {3, 2, 6, 7, -2}, {6, 2, 9, 7, 1} };
- const int dy_s[NY][5] = { {2, 0, 7, 3, 1}, {2, 3, 7, 6, -2}, {2, 6, 7, 9, 1} };
- const int dxy_s[NXY][5] = { {1, 1, 4, 4, 1}, {5, 1, 8, 4, -1}, {1, 5, 4, 8, -1}, {5, 5, 8, 8, 1} };
- 首先定义了 方框大小为 9 时的 3 个类 Haar 模板: dx_s、dy_s、dxy_s,分别代表 x二阶偏导、y二阶偏导、xy交叉二阶偏导。这三个模板作为基础模板。
- SurfHF Dx[NX], Dy[NY], Dxy[NXY];
- if( size > sum.rows-1 || size > sum.cols-1 )
- return;
- resizeHaarPattern( dx_s , Dx , NX , 9, size, sum.cols );
- resizeHaarPattern( dy_s , Dy , NY , 9, size, sum.cols );
- resizeHaarPattern( dxy_s, Dxy, NXY, 9, size, sum.cols );
- 根据 3 个基础模板生成对应尺度的类Haar模板: Dx、Dy、Dxy。
- /* The integral image 'sum' is one pixel bigger than the source image */
- int samples_i = 1+(sum.rows-1-size)/sampleStep;
- int samples_j = 1+(sum.cols-1-size)/sampleStep;
- /* Ignore pixels where some of the kernel is outside the image */
- int margin = (size/2)/sampleStep;
- for( int i = 0; i < samples_i; i++ )
- {
- const int* sum_ptr = sum.ptr<int>(i*sampleStep);
- float* det_ptr = &det.at<float>(i+margin, margin);
- float* trace_ptr = &trace.at<float>(i+margin, margin);
- for( int j = 0; j < samples_j; j++ )
- {
- float dx = calcHaarPattern( sum_ptr, Dx , 3 );
- float dy = calcHaarPattern( sum_ptr, Dy , 3 );
- float dxy = calcHaarPattern( sum_ptr, Dxy, 4 );
- sum_ptr += sampleStep;
- det_ptr[j] = dx*dy - 0.81f*dxy*dxy;
- 注意上面这个等式,右边第二项有个 0.81 的系数,或 0.9^2 。
- trace_ptr[j] = dx + dy;
- }
- }
- 用积分图像和类haar模板快速计算Hessian矩阵的行列式和迹。
- }
SURFFindInvoker: 搜索极值点(特征点)
最后只剩下在尺度空间寻找极值点了。在各层图像的dets矩阵中寻找极值点,并用二阶泰勒展开拟合极值点的精确位置。特征点的尺寸(size)、符号(凹凸性)也是在这里被计算。
- // Multi-threaded search of the scale-space pyramid for keypoints
- struct SURFFindInvoker : ParallelLoopBody
- {
- SURFFindInvoker( const Mat& _sum, const Mat& _mask_sum,
- const vector<Mat>& _dets, const vector<Mat>& _traces,
- const vector<int>& _sizes, const vector<int>& _sampleSteps,
- const vector<int>& _middleIndices, vector<KeyPoint>& _keypoints,
- int _nOctaveLayers, float _hessianThreshold )
- {
- sum = &_sum;
- mask_sum = &_mask_sum;
- dets = &_dets;
- traces = &_traces;
- sizes = &_sizes;
- sampleSteps = &_sampleSteps;
- middleIndices = &_middleIndices;
- keypoints = &_keypoints;
- nOctaveLayers = _nOctaveLayers;
- hessianThreshold = _hessianThreshold;
- }
- static void 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 );
- void operator()(const Range& range) const
- {
- for( int i=range.start; i<range.end; i++ )
- {
- int layer = (*middleIndices)[i];
- int octave = i / nOctaveLayers;
- findMaximaInLayer( *sum, *mask_sum, *dets, *traces, *sizes,
- *keypoints, octave, layer, hessianThreshold,
- (*sampleSteps)[layer] );
- 对每个中间层调用 findMaximaInLayer() 寻找极值点
- }
- }
- const Mat *sum;
- const Mat *mask_sum;
- const vector<Mat>* dets;
- const vector<Mat>* traces;
- const vector<int>* sizes;
- const vector<int>* sampleSteps;
- const vector<int>* middleIndices;
- vector<KeyPoint>* keypoints;
- int nOctaveLayers;
- float hessianThreshold;
- static Mutex findMaximaInLayer_m;
- };
- /*
- * Find the maxima in the determinant of the Hessian in a layer of the
- * scale-space pyramid
- */
- 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;
- 积分图像sum的宽和高都比原始图像大 1 ,计算 sum.rows 和 sum.cols 都需要减 1 。
- // Ignore pixels without a 3x3x3 neighbourhood in the layer above
- int margin = (sizes[layer+1]/2)/sampleStep+1;
- margin 是边缘区域的大小。在相邻的三层中(第layer-1、第layer和第layer+1 层)中,第layer+1层图像的边缘区域最多(因为它的滤波方框最大),所以参考第layer+1层来计算边缘区域(方框边长的一半)。
- 因为第layer层是一个中间层,所以这相邻的三层肯定都在同一个 octave 内,他们有相同的采样间隔 sampleStep。
- if( !mask_sum.empty() )
- resizeHaarPattern( dm, &Dm, NM, 9, size, mask_sum.cols );
- 根据dm的定义,可知Dm作用于某个积分图像后,仅相当于对源图像中一个大小为 size 的矩形区域求和。
- 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++ )
- {
- 循环检测极值点。不检测两边 margin 以外的边缘区域。
- float val0 = det_ptr[j];
- if( val0 > hessianThreshold )
- {
- 只对响应大于阈值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);
- (sum_i,sum_j) 代表以当前点(i,j)为中心的方框的左上角坐标。
- /* 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],
- 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] } };
- 将 3*3*3 邻域的数据提取到一个数组 N9 中,方便后面比较。
- /* 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;
- }
- calcHaarPattern( mask_ptr, &Dm, 1 ) 相当于对mask矩阵中以当前位置为中心的方框(大小为size)求和: Dm代表对一个大小为size的方框求和。
- 这一步检测当前点是否被mask: 如果以当前点为中心的方框中被mask的点数少于一半,即 mval < 0.5, 则说明当前点未被 mask,不再对其进行极值点检测。
- 为什么不直接检测方框中心点的 mask 值来判断是否跳过当前点?而非要用mask的方框区域和来判断?为什么后者更合理?
- /* Non-maxima suppression. val0 is at N9[1][4]*/
- if( val0 > N9[0][0] && val0 > N9[0][1] && val0 > N9[0][2] &&
- 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] )
- {
- 只要正的极大值( val0 > hessianThreshold 已保证了val0是正),而不要负的最小值。
- Hessian矩阵是二次型,正定或负定都将导致Hessian行列式为正。如果既非正定又非负定(这样的点不是斑点),则其Hessian行列式为负,不会被检测到。
- 正因如此,最后被检测到特征点的Hessian矩阵不是正定就是负定,因此再检测Hessian矩阵的迹的正负号即可知道到底是正定还是负定,这也决定了该特征点是凹斑点还是凸斑点
- /* 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;
- 其中 sum_i 和 sum_j 的定义为(参见前面代码):
- int sum_i = sampleStep*(i-(size/2)/sampleStep);
- int sum_j = sampleStep*(j-(size/2)/sampleStep);
- 可见二者分别代表以当前点(i,j)为中心的方框的左上角,center_i和center_j则是重新计算的该方框的中心点。
- 为什么不直接写成:
- center_i = i * sampleStep,
- center_j = j * sampleStep
- KeyPoint kpt( center_j, center_i, (float)sizes[layer],
- -1, val0, octave, CV_SIGN(trace_ptr[j]) );
- 最后一个参数是提取 hessian矩阵的迹的正负号,它反映了当前的特征点是凹的还是凸的。
- 以center_i和center_j作为关键点的坐标,以当前层的方框的size作为关键点的size。
- /* Interpolate maxima location within the 3x3x3 neighbourhood */
- int ds = size - sizes[layer-1];
- int interp_ok = interpolateKeypoint( N9, sampleStep, sampleStep, ds, kpt );
- 如果符合极值点条件,再对极值点进行二次泰勒拟合以求精确的极值点位置和size
- /* 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 );*/
- cv::AutoLock lock(findMaximaInLayer_m);
- keypoints.push_back(kpt);
- }
- }
- }
- }
- }
- }