1 /* Original code has been submitted by Liu Liu. Here is the copyright.
2 ----------------------------------------------------------------------------------
3 * An OpenCV Implementation of SURF
4 * Further Information Refer to "SURF: Speed-Up Robust Feature"
6 * liuliu.1987+opencv@gmail.com
8 * There are still serveral lacks for this experimental implementation:
9 * 1.The interpolation of sub-pixel mentioned in article was not implemented yet;
10 * 2.A comparision with original libSurf.so shows that the hessian detector is not a 100% match to their implementation;
11 * 3.Due to above reasons, I recommanded the original one for study and reuse;
13 * However, the speed of this implementation is something comparable to original one.
15 * Copyright© 2008, Liu Liu All rights reserved.
17 * Redistribution and use in source and binary forms, with or
18 * without modification, are permitted provided that the following
20 * Redistributions of source code must retain the above
21 * copyright notice, this list of conditions and the following
23 * Redistributions in binary form must reproduce the above
24 * copyright notice, this list of conditions and the following
25 * disclaimer in the documentation and/or other materials
26 * provided with the distribution.
27 * The name of Contributor may not be used to endorse or
28 * promote products derived from this software without
29 * specific prior written permission.
31 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND
32 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,
33 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
34 * MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
35 * DISCLAIMED. IN NO EVENT SHALL THE CONTRIBUTORS BE LIABLE FOR ANY
36 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
37 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
38 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA,
39 * OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
40 * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
41 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
42 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
47 The following changes have been made, comparing to the original contribution:
48 1. A lot of small optimizations, less memory allocations, got rid of global buffers
49 2. Reversed order of cvGetQuadrangleSubPix and cvResize calls; probably less accurate, but much faster
50 3. The descriptor computing part (which is most expensive) is threaded using OpenMP
51 (subpixel-accurate keypoint localization and scale estimation are still TBD)
55 KeyPoint position and scale interpolation has been implemented as described in
56 the Brown and Lowe paper cited by the SURF paper.
58 The sampling step along the x and y axes of the image for the determinant of the
59 Hessian is now the same for each layer in an octave. While this increases the
60 computation time, it ensures that a true 3x3x3 neighbourhood exists, with
61 samples calculated at the same position in the layers above and below. This
62 results in improved maxima detection and non-maxima suppression, and I think it
63 is consistent with the description in the SURF paper.
65 The wavelet size sampling interval has also been made consistent. The wavelet
66 size at the first layer of the first octave is now 9 instead of 7. Along with
67 regular position sampling steps, this makes location and scale interpolation
68 easy. I think this is consistent with the SURF paper and original
71 The scaling of the wavelet parameters has been fixed to ensure that the patterns
72 are symmetric around the centre. Previously the truncation caused by integer
73 division in the scaling ratio caused a bias towards the top left of the wavelet,
74 resulting in inconsistent keypoint positions.
76 The matrices for the determinant and trace of the Hessian are now reused in each
79 The extraction of the patch of pixels surrounding a keypoint used to build a
80 descriptor has been simplified.
82 KeyPoint descriptor normalisation has been changed from normalising each 4x4
83 cell (resulting in a descriptor of magnitude 16) to normalising the entire
84 descriptor to magnitude 1.
86 The default number of octaves has been increased from 3 to 4 to match the
87 original SURF binary default. The increase in computation time is minimal since
88 the higher octaves are sampled sparsely.
90 The default number of layers per octave has been reduced from 3 to 2, to prevent
91 redundant calculation of similar sizes in consecutive octaves. This decreases
92 computation time. The number of features extracted may be less, however the
93 additional features were mostly redundant.
95 The radius of the circle of gradient samples used to assign an orientation has
96 been increased from 4 to 6 to match the description in the SURF paper. This is
97 now defined by ORI_RADIUS, and could be made into a parameter.
99 The size of the sliding window used in orientation assignment has been reduced
100 from 120 to 60 degrees to match the description in the SURF paper. This is now
101 defined by ORI_WIN, and could be made into a parameter.
103 Other options like HAAR_SIZE0, HAAR_SIZE_INC, SAMPLE_STEP0, ORI_SEARCH_INC,
104 ORI_SIGMA and DESC_SIGMA have been separated from the code and documented.
105 These could also be made into parameters.
107 Modifications by Ian Mahon
112 CvSURFParams cvSURFParams(double threshold, int extended)
115 params.hessianThreshold = threshold;
116 params.extended = extended;
118 params.nOctaveLayers = 2;
129 icvCalcHaarPattern( const int* origin, const CvSurfHF* f, int n )
132 for( int k = 0; k < n; k++ )
133 d += (origin[f[k].p0] + origin[f[k].p3] - origin[f[k].p1] - origin[f[k].p2])*f[k].w;
138 icvResizeHaarPattern( const int src[][5], CvSurfHF* dst, int n, int oldSize, int newSize, int widthStep )
140 float ratio = (float)newSize/oldSize;
141 for( int k = 0; k < n; k++ )
143 int dx1 = cvRound( ratio*src[k][0] );
144 int dy1 = cvRound( ratio*src[k][1] );
145 int dx2 = cvRound( ratio*src[k][2] );
146 int dy2 = cvRound( ratio*src[k][3] );
147 dst[k].p0 = dy1*widthStep + dx1;
148 dst[k].p1 = dy2*widthStep + dx1;
149 dst[k].p2 = dy1*widthStep + dx2;
150 dst[k].p3 = dy2*widthStep + dx2;
151 dst[k].w = src[k][4]/((float)(dx2-dx1)*(dy2-dy1));
156 * Maxima location interpolation as described in "Invariant Features from
157 * Interest Point Groups" by Matthew Brown and David Lowe. This is performed by
158 * fitting a 3D quadratic to a set of neighbouring samples.
160 * The gradient vector and Hessian matrix at the initial keypoint location are
161 * approximated using central differences. The linear system Ax = b is then
162 * solved, where A is the Hessian, b is the negative gradient, and x is the
163 * offset of the interpolated maxima coordinates from the initial estimate.
164 * This is equivalent to an iteration of Netwon's optimisation algorithm.
166 * N9 contains the samples in the 3x3x3 neighbourhood of the maxima
167 * dx is the sampling step in x
168 * dy is the sampling step in y
169 * ds is the sampling step in size
170 * point contains the keypoint coordinates and scale to be modified
172 * Return value is 1 if interpolation was successful, 0 on failure.
175 icvInterpolateKeypoint( float N9[3][9], int dx, int dy, int ds, CvSURFPoint *point )
178 float A[9], x[3], b[3];
179 CvMat _A = cvMat(3, 3, CV_32F, A);
180 CvMat _x = cvMat(3, 1, CV_32F, x);
181 CvMat _b = cvMat(3, 1, CV_32F, b);
183 b[0] = -(N9[1][5]-N9[1][3])/2; /* Negative 1st deriv with respect to x */
184 b[1] = -(N9[1][7]-N9[1][1])/2; /* Negative 1st deriv with respect to y */
185 b[2] = -(N9[2][4]-N9[0][4])/2; /* Negative 1st deriv with respect to s */
187 A[0] = N9[1][3]-2*N9[1][4]+N9[1][5]; /* 2nd deriv x, x */
188 A[1] = (N9[1][8]-N9[1][6]-N9[1][2]+N9[1][0])/4; /* 2nd deriv x, y */
189 A[2] = (N9[2][5]-N9[2][3]-N9[0][5]+N9[0][3])/4; /* 2nd deriv x, s */
190 A[3] = A[1]; /* 2nd deriv y, x */
191 A[4] = N9[1][1]-2*N9[1][4]+N9[1][7]; /* 2nd deriv y, y */
192 A[5] = (N9[2][7]-N9[2][1]-N9[0][7]+N9[0][1])/4; /* 2nd deriv y, s */
193 A[6] = A[2]; /* 2nd deriv s, x */
194 A[7] = A[5]; /* 2nd deriv s, y */
195 A[8] = N9[0][4]-2*N9[1][4]+N9[2][4]; /* 2nd deriv s, s */
197 solve_ok = cvSolve( &_A, &_b, &_x );
200 point->pt.x += x[0]*dx;
201 point->pt.y += x[1]*dy;
202 point->size = cvRound( point->size + x[2]*ds );
208 /* Wavelet size at first layer of first octave. */
209 const int HAAR_SIZE0 = 9;
211 /* Wavelet size increment between layers. This should be an even number,
212 such that the wavelet sizes in an octave are either all even or all odd.
213 This ensures that when looking for the neighbours of a sample, the layers
214 above and below are aligned correctly. */
215 const int HAAR_SIZE_INC = 6;
218 static CvSeq* icvFastHessianDetector( const CvMat* sum, const CvMat* mask_sum,
219 CvMemStorage* storage, const CvSURFParams* params )
221 CvSeq* points = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvSURFPoint), storage );
223 /* Sampling step along image x and y axes at first octave. This is doubled
224 for each additional octave. WARNING: Increasing this improves speed,
225 however keypoint extraction becomes unreliable. */
226 const int SAMPLE_STEP0 = 1;
230 const int NX=3, NY=3, NXY=4, NM=1;
231 const int dx_s[NX][5] = { {0, 2, 3, 7, 1}, {3, 2, 6, 7, -2}, {6, 2, 9, 7, 1} };
232 const int dy_s[NY][5] = { {2, 0, 7, 3, 1}, {2, 3, 7, 6, -2}, {2, 6, 7, 9, 1} };
233 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} };
234 const int dm[NM][5] = { {0, 0, 9, 9, 1} };
235 CvSurfHF Dx[NX], Dy[NY], Dxy[NXY], Dm;
237 CvMat** dets = (CvMat**)cvStackAlloc((params->nOctaveLayers+2)*sizeof(dets[0]));
238 CvMat** traces = (CvMat**)cvStackAlloc((params->nOctaveLayers+2)*sizeof(traces[0]));
239 int *sizes = (int*)cvStackAlloc((params->nOctaveLayers+2)*sizeof(sizes[0]));
241 double dx = 0, dy = 0, dxy = 0;
242 int octave, layer, sampleStep, size, margin;
244 int i, j, sum_i, sum_j;
246 float *det_ptr, *trace_ptr;
248 /* Allocate enough space for hessian determinant and trace matrices at the
249 first octave. Clearing these initially or between octaves is not
250 required, since all values that are accessed are first calculated */
251 for( layer = 0; layer <= params->nOctaveLayers+1; layer++ )
253 dets[layer] = cvCreateMat( (sum->rows-1)/SAMPLE_STEP0, (sum->cols-1)/SAMPLE_STEP0, CV_32FC1 );
254 traces[layer] = cvCreateMat( (sum->rows-1)/SAMPLE_STEP0, (sum->cols-1)/SAMPLE_STEP0, CV_32FC1 );
257 for( octave = 0, sampleStep=SAMPLE_STEP0; octave < params->nOctaves; octave++, sampleStep*=2 )
259 /* Hessian determinant and trace sample array size in this octave */
260 rows = (sum->rows-1)/sampleStep;
261 cols = (sum->cols-1)/sampleStep;
263 /* Calculate the determinant and trace of the hessian */
264 for( layer = 0; layer <= params->nOctaveLayers+1; layer++ )
266 sizes[layer] = size = (HAAR_SIZE0+HAAR_SIZE_INC*layer)<<octave;
267 icvResizeHaarPattern( dx_s, Dx, NX, 9, size, sum->cols );
268 icvResizeHaarPattern( dy_s, Dy, NY, 9, size, sum->cols );
269 icvResizeHaarPattern( dxy_s, Dxy, NXY, 9, size, sum->cols );
270 /*printf( "octave=%d layer=%d size=%d rows=%d cols=%d\n", octave, layer, size, rows, cols );*/
272 margin = (size/2)/sampleStep;
273 for( sum_i=0, i=margin; sum_i<=(sum->rows-1)-size; sum_i+=sampleStep, i++ )
275 s_ptr = sum->data.i + sum_i*sum->cols;
276 det_ptr = dets[layer]->data.fl + i*dets[layer]->cols + margin;
277 trace_ptr = traces[layer]->data.fl + i*traces[layer]->cols + margin;
278 for( sum_j=0, j=margin; sum_j<=(sum->cols-1)-size; sum_j+=sampleStep, j++ )
280 dx = icvCalcHaarPattern( s_ptr, Dx, 3 );
281 dy = icvCalcHaarPattern( s_ptr, Dy, 3 );
282 dxy = icvCalcHaarPattern( s_ptr, Dxy, 4 );
284 *det_ptr++ = (float)(dx*dy - 0.81*dxy*dxy);
285 *trace_ptr++ = (float)(dx + dy);
290 /* Find maxima in the determinant of the hessian */
291 for( layer = 1; layer <= params->nOctaveLayers; layer++ )
294 icvResizeHaarPattern( dm, &Dm, NM, 9, size, mask_sum ? mask_sum->cols : sum->cols );
296 /* Ignore pixels without a 3x3 neighbourhood in the layer above */
297 margin = (sizes[layer+1]/2)/sampleStep+1;
298 for( i = margin; i < rows-margin; i++ )
300 det_ptr = dets[layer]->data.fl + i*dets[layer]->cols;
301 trace_ptr = traces[layer]->data.fl + i*traces[layer]->cols;
302 for( j = margin; j < cols-margin; j++ )
304 float val0 = det_ptr[j];
305 if( val0 > params->hessianThreshold )
307 /* Coordinates for the start of the wavelet in the sum image. There
308 is some integer division involved, so don't try to simplify this
309 (cancel out sampleStep) without checking the result is the same */
310 int sum_i = sampleStep*(i-(size/2)/sampleStep);
311 int sum_j = sampleStep*(j-(size/2)/sampleStep);
313 /* The 3x3x3 neighbouring samples around the maxima.
314 The maxima is included at N9[1][4] */
315 int c = dets[layer]->cols;
316 const float *det1 = dets[layer-1]->data.fl + i*c + j;
317 const float *det2 = dets[layer]->data.fl + i*c + j;
318 const float *det3 = dets[layer+1]->data.fl + i*c + j;
319 float N9[3][9] = { { det1[-c-1], det1[-c], det1[-c+1],
320 det1[-1] , det1[0] , det1[1],
321 det1[c-1] , det1[c] , det1[c+1] },
322 { det2[-c-1], det2[-c], det2[-c+1],
323 det2[-1] , det2[0] , det2[1],
324 det2[c-1] , det2[c] , det2[c+1 ] },
325 { det3[-c-1], det3[-c], det3[-c+1],
326 det3[-1 ], det3[0] , det3[1],
327 det3[c-1] , det3[c] , det3[c+1 ] } };
329 /* Check the mask - why not just check the mask at the center of the wavelet? */
332 const int* mask_ptr = mask_sum->data.i + mask_sum->cols*sum_i + sum_j;
333 float mval = icvCalcHaarPattern( mask_ptr, &Dm, 1 );
338 /* Non-maxima suppression. val0 is at N9[1][4]*/
339 if( val0 > N9[0][0] && val0 > N9[0][1] && val0 > N9[0][2] &&
340 val0 > N9[0][3] && val0 > N9[0][4] && val0 > N9[0][5] &&
341 val0 > N9[0][6] && val0 > N9[0][7] && val0 > N9[0][8] &&
342 val0 > N9[1][0] && val0 > N9[1][1] && val0 > N9[1][2] &&
343 val0 > N9[1][3] && val0 > N9[1][5] &&
344 val0 > N9[1][6] && val0 > N9[1][7] && val0 > N9[1][8] &&
345 val0 > N9[2][0] && val0 > N9[2][1] && val0 > N9[2][2] &&
346 val0 > N9[2][3] && val0 > N9[2][4] && val0 > N9[2][5] &&
347 val0 > N9[2][6] && val0 > N9[2][7] && val0 > N9[2][8] )
349 /* Calculate the wavelet center coordinates for the maxima */
350 double center_i = sum_i + (double)(size-1)/2;
351 double center_j = sum_j + (double)(size-1)/2;
353 CvSURFPoint point = cvSURFPoint( cvPoint2D32f(center_j,center_i),
354 CV_SIGN(trace_ptr[j]), sizes[layer], 0, val0 );
356 /* Interpolate maxima location within the 3x3x3 neighbourhood */
357 int ds = sizes[layer]-sizes[layer-1];
358 int interp_ok = icvInterpolateKeypoint( N9, sampleStep, sampleStep, ds, &point );
360 /* Sometimes the interpolation step gives a negative size etc. */
361 if( interp_ok && point.size >= 1 &&
362 point.pt.x >= 0 && point.pt.x <= (sum->cols-1) &&
363 point.pt.y >= 0 && point.pt.y <= (sum->rows-1) )
365 /*printf( "KeyPoint %f %f %d\n", point.pt.x, point.pt.y, point.size );*/
366 cvSeqPush( points, &point );
376 for( layer = 0; layer <= params->nOctaveLayers+1; layer++ )
378 cvReleaseMat( &dets[layer] );
379 cvReleaseMat( &traces[layer] );
387 cvExtractSURF( const CvArr* _img, const CvArr* _mask,
388 CvSeq** _keypoints, CvSeq** _descriptors,
389 CvMemStorage* storage, CvSURFParams params,
390 int useProvidedKeyPts)
392 CvMat *sum = 0, *mask1 = 0, *mask_sum = 0, **win_bufs = 0;
394 if( _keypoints && !useProvidedKeyPts ) // If useProvidedKeyPts!=0 we'll use current contents of "*_keypoints"
399 /* Radius of the circle in which to sample gradients to assign an
401 const int ORI_RADIUS = 6;
403 /* The size of the sliding window (in degrees) used to assign an
405 const int ORI_WIN = 60;
407 /* Increment used for the orientation sliding window (in degrees) */
408 const int ORI_SEARCH_INC = 5;
410 /* Standard deviation of the Gaussian used to weight the gradient samples
411 used to assign an orientation */
412 const float ORI_SIGMA = 2.5f;
414 /* Standard deviation of the Gaussian used to weight the gradient samples
415 used to build a keypoint descriptor */
416 const float DESC_SIGMA = 3.3f;
419 /* X and Y gradient wavelet data */
420 const int NX=2, NY=2;
421 int dx_s[NX][5] = {{0, 0, 2, 4, -1}, {2, 0, 4, 4, 1}};
422 int dy_s[NY][5] = {{0, 0, 4, 2, 1}, {0, 2, 4, 4, -1}};
424 CvSeq *keypoints, *descriptors = 0;
425 CvMat imghdr, *img = cvGetMat(_img, &imghdr);
426 CvMat maskhdr, *mask = _mask ? cvGetMat(_mask, &maskhdr) : 0;
428 const int max_ori_samples = (2*ORI_RADIUS+1)*(2*ORI_RADIUS+1);
429 int descriptor_size = params.extended ? 128 : 64;
430 const int descriptor_data_type = CV_32F;
431 const int PATCH_SZ = 20;
432 float DW[PATCH_SZ][PATCH_SZ];
433 CvMat _DW = cvMat(PATCH_SZ, PATCH_SZ, CV_32F, DW);
434 CvPoint apt[max_ori_samples];
435 float apt_w[max_ori_samples];
436 int i, j, k, nangle0 = 0, N;
437 int nthreads = cvGetNumThreads();
440 CV_Assert(CV_MAT_TYPE(img->type) == CV_8UC1);
441 CV_Assert(mask == 0 || (CV_ARE_SIZES_EQ(img,mask) && CV_MAT_TYPE(mask->type) == CV_8UC1));
442 CV_Assert(storage != 0);
443 CV_Assert(params.hessianThreshold >= 0);
444 CV_Assert(params.nOctaves > 0);
445 CV_Assert(params.nOctaveLayers > 0);
447 sum = cvCreateMat( img->height+1, img->width+1, CV_32SC1 );
448 cvIntegral( img, sum );
450 // Compute keypoints only if we are not asked for evaluating the descriptors are some given locations:
451 if (!useProvidedKeyPts)
455 mask1 = cvCreateMat( img->height, img->width, CV_8UC1 );
456 mask_sum = cvCreateMat( img->height+1, img->width+1, CV_32SC1 );
457 cvMinS( mask, 1, mask1 );
458 cvIntegral( mask1, mask_sum );
460 keypoints = icvFastHessianDetector( sum, mask_sum, storage, ¶ms );
464 CV_Assert(useProvidedKeyPts && (_keypoints != 0) && (*_keypoints != 0));
465 keypoints = *_keypoints;
468 N = keypoints->total;
471 descriptors = cvCreateSeq( 0, sizeof(CvSeq),
472 descriptor_size*CV_ELEM_SIZE(descriptor_data_type), storage );
473 cvSeqPushMulti( descriptors, 0, N );
476 /* Coordinates and weights of samples used to calculate orientation */
477 cv::Mat _G = cv::getGaussianKernel( 2*ORI_RADIUS+1, ORI_SIGMA, CV_32F );
478 const float* G = (const float*)_G.data;
480 for( i = -ORI_RADIUS; i <= ORI_RADIUS; i++ )
482 for( j = -ORI_RADIUS; j <= ORI_RADIUS; j++ )
484 if( i*i + j*j <= ORI_RADIUS*ORI_RADIUS )
486 apt[nangle0] = cvPoint(j,i);
487 apt_w[nangle0++] = G[i+ORI_RADIUS]*G[j+ORI_RADIUS];
492 /* Gaussian used to weight descriptor samples */
494 double c2 = 1./(DESC_SIGMA*DESC_SIGMA*2);
496 for( i = 0; i < PATCH_SZ; i++ )
498 for( j = 0; j < PATCH_SZ; j++ )
500 double x = j - (float)(PATCH_SZ-1)/2, y = i - (float)(PATCH_SZ-1)/2;
501 double val = exp(-(x*x+y*y)*c2);
502 DW[i][j] = (float)val;
506 cvScale( &_DW, &_DW, 1./gs );
509 win_bufs = (CvMat**)cvAlloc(nthreads*sizeof(win_bufs[0]));
510 for( i = 0; i < nthreads; i++ )
514 #pragma omp parallel for num_threads(nthreads) schedule(dynamic)
516 for( k = 0; k < N; k++ )
518 const int* sum_ptr = sum->data.i;
519 int sum_cols = sum->cols;
520 int i, j, kk, x, y, nangle;
521 float X[max_ori_samples], Y[max_ori_samples], angle[max_ori_samples];
522 uchar PATCH[PATCH_SZ+1][PATCH_SZ+1];
523 float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ];
524 CvMat _X = cvMat(1, max_ori_samples, CV_32F, X);
525 CvMat _Y = cvMat(1, max_ori_samples, CV_32F, Y);
526 CvMat _angle = cvMat(1, max_ori_samples, CV_32F, angle);
527 CvMat _patch = cvMat(PATCH_SZ+1, PATCH_SZ+1, CV_8U, PATCH);
529 CvSurfHF dx_t[NX], dy_t[NY];
530 int thread_idx = cvGetThreadNum();
532 CvSURFPoint* kp = (CvSURFPoint*)cvGetSeqElem( keypoints, k );
534 CvPoint2D32f center = kp->pt;
536 /* The sampling intervals and wavelet sized for selecting an orientation
537 and building the keypoint descriptor are defined relative to 's' */
538 float s = (float)size*1.2f/9.0f;
540 /* To find the dominant orientation, the gradients in x and y are
541 sampled in a circle of radius 6s using wavelets of size 4s.
542 We ensure the gradient wavelet size is even to ensure the
543 wavelet pattern is balanced and symmetric around its center */
544 int grad_wav_size = 2*cvRound( 2*s );
545 if ( sum->rows < grad_wav_size || sum->cols < grad_wav_size )
547 /* when grad_wav_size is too big,
548 * the sampling of gradient will be meaningless
549 * mark keypoint for deletion. */
553 icvResizeHaarPattern( dx_s, dx_t, NX, 4, grad_wav_size, sum->cols );
554 icvResizeHaarPattern( dy_s, dy_t, NY, 4, grad_wav_size, sum->cols );
555 for( kk = 0, nangle = 0; kk < nangle0; kk++ )
559 x = cvRound( center.x + apt[kk].x*s - (float)(grad_wav_size-1)/2 );
560 y = cvRound( center.y + apt[kk].y*s - (float)(grad_wav_size-1)/2 );
561 if( (unsigned)y >= (unsigned)(sum->rows - grad_wav_size) ||
562 (unsigned)x >= (unsigned)(sum->cols - grad_wav_size) )
564 ptr = sum_ptr + x + y*sum_cols;
565 vx = icvCalcHaarPattern( ptr, dx_t, 2 );
566 vy = icvCalcHaarPattern( ptr, dy_t, 2 );
567 X[nangle] = vx*apt_w[kk]; Y[nangle] = vy*apt_w[kk];
572 /* No gradient could be sampled because the keypoint is too
573 * near too one or more of the sides of the image. As we
574 * therefore cannot find a dominant direction, we skip this
575 * keypoint and mark it for later deletion from the sequence. */
579 _X.cols = _Y.cols = _angle.cols = nangle;
580 cvCartToPolar( &_X, &_Y, 0, &_angle, 1 );
582 float bestx = 0, besty = 0, descriptor_mod = 0;
583 for( i = 0; i < 360; i += ORI_SEARCH_INC )
585 float sumx = 0, sumy = 0, temp_mod;
586 for( j = 0; j < nangle; j++ )
588 int d = abs(cvRound(angle[j]) - i);
589 if( d < ORI_WIN/2 || d > 360-ORI_WIN/2 )
595 temp_mod = sumx*sumx + sumy*sumy;
596 if( temp_mod > descriptor_mod )
598 descriptor_mod = temp_mod;
604 float descriptor_dir = cvFastArctan( besty, bestx );
605 kp->dir = descriptor_dir;
610 descriptor_dir *= (float)(CV_PI/180);
612 /* Extract a window of pixels around the keypoint of size 20s */
613 int win_size = (int)((PATCH_SZ+1)*s);
614 if( win_bufs[thread_idx] == 0 || win_bufs[thread_idx]->cols < win_size*win_size )
616 cvReleaseMat( &win_bufs[thread_idx] );
617 win_bufs[thread_idx] = cvCreateMat( 1, win_size*win_size, CV_8U );
620 CvMat win = cvMat(win_size, win_size, CV_8U, win_bufs[thread_idx]->data.ptr);
621 float sin_dir = sin(descriptor_dir);
622 float cos_dir = cos(descriptor_dir) ;
624 /* Subpixel interpolation version (slower). Subpixel not required since
625 the pixels will all get averaged when we scale down to 20 pixels */
627 float w[] = { cos_dir, sin_dir, center.x,
628 -sin_dir, cos_dir , center.y };
629 CvMat W = cvMat(2, 3, CV_32F, w);
630 cvGetQuadrangleSubPix( img, &win, &W );
633 /* Nearest neighbour version (faster) */
634 float win_offset = -(float)(win_size-1)/2;
635 float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;
636 float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;
637 uchar* WIN = win.data.ptr;
638 for( i=0; i<win_size; i++, start_x+=sin_dir, start_y+=cos_dir )
640 float pixel_x = start_x;
641 float pixel_y = start_y;
642 for( j=0; j<win_size; j++, pixel_x+=cos_dir, pixel_y-=sin_dir )
644 int x = cvRound( pixel_x );
645 int y = cvRound( pixel_y );
648 x = MIN( x, img->cols-1 );
649 y = MIN( y, img->rows-1 );
650 WIN[i*win_size + j] = img->data.ptr[y*img->step+x];
654 /* Scale the window to size PATCH_SZ so each pixel's size is s. This
655 makes calculating the gradients with wavelets of size 2s easy */
656 cvResize( &win, &_patch, CV_INTER_AREA );
658 /* Calculate gradients in x and y with wavelets of size 2s */
659 for( i = 0; i < PATCH_SZ; i++ )
660 for( j = 0; j < PATCH_SZ; j++ )
663 float vx = (PATCH[i][j+1] - PATCH[i][j] + PATCH[i+1][j+1] - PATCH[i+1][j])*dw;
664 float vy = (PATCH[i+1][j] - PATCH[i][j] + PATCH[i+1][j+1] - PATCH[i][j+1])*dw;
669 /* Construct the descriptor */
670 vec = (float*)cvGetSeqElem( descriptors, k );
671 for( kk = 0; kk < (int)(descriptors->elem_size/sizeof(vec[0])); kk++ )
673 double square_mag = 0;
674 if( params.extended )
676 /* 128-bin descriptor */
677 for( i = 0; i < 4; i++ )
678 for( j = 0; j < 4; j++ )
680 for( y = i*5; y < i*5+5; y++ )
682 for( x = j*5; x < j*5+5; x++ )
684 float tx = DX[y][x], ty = DY[y][x];
688 vec[1] += (float)fabs(tx);
691 vec[3] += (float)fabs(tx);
696 vec[5] += (float)fabs(ty);
699 vec[7] += (float)fabs(ty);
703 for( kk = 0; kk < 8; kk++ )
704 square_mag += vec[kk]*vec[kk];
710 /* 64-bin descriptor */
711 for( i = 0; i < 4; i++ )
712 for( j = 0; j < 4; j++ )
714 for( y = i*5; y < i*5+5; y++ )
716 for( x = j*5; x < j*5+5; x++ )
718 float tx = DX[y][x], ty = DY[y][x];
719 vec[0] += tx; vec[1] += ty;
720 vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty);
723 for( kk = 0; kk < 4; kk++ )
724 square_mag += vec[kk]*vec[kk];
729 /* unit vector is essential for contrast invariance */
730 vec = (float*)cvGetSeqElem( descriptors, k );
731 double scale = 1./(sqrt(square_mag) + DBL_EPSILON);
732 for( kk = 0; kk < descriptor_size; kk++ )
733 vec[kk] = (float)(vec[kk]*scale);
736 /* remove keypoints that were marked for deletion */
737 for ( i = 0; i < N; i++ )
739 CvSURFPoint* kp = (CvSURFPoint*)cvGetSeqElem( keypoints, i );
740 if ( kp->size == -1 )
742 cvSeqRemove( keypoints, i );
744 cvSeqRemove( descriptors, i );
750 for( i = 0; i < nthreads; i++ )
751 cvReleaseMat( &win_bufs[i] );
753 if( _keypoints && !useProvidedKeyPts )
754 *_keypoints = keypoints;
756 *_descriptors = descriptors;
758 cvReleaseMat( &sum );
759 if (mask1) cvReleaseMat( &mask1 );
760 if (mask_sum) cvReleaseMat( &mask_sum );
770 hessianThreshold = 100;
776 SURF::SURF(double _threshold, int _nOctaves, int _nOctaveLayers, bool _extended)
778 hessianThreshold = _threshold;
779 extended = _extended;
780 nOctaves = _nOctaves;
781 nOctaveLayers = _nOctaveLayers;
784 int SURF::descriptorSize() const { return extended ? 128 : 64; }
787 static int getPointOctave(const CvSURFPoint& kpt, const CvSURFParams& params)
789 int octave = 0, layer = 0, best_octave = 0;
790 float min_diff = FLT_MAX;
791 for( octave = 1; octave < params.nOctaves; octave++ )
792 for( layer = 0; layer < params.nOctaveLayers; layer++ )
794 float diff = std::abs(kpt.size - (float)((HAAR_SIZE0 + HAAR_SIZE_INC*layer) << octave));
795 if( min_diff > diff )
798 best_octave = octave;
807 void SURF::operator()(const Mat& img, const Mat& mask,
808 vector<KeyPoint>& keypoints) const
810 CvMat _img = img, _mask, *pmask = 0;
812 pmask = &(_mask = mask);
813 MemStorage storage(cvCreateMemStorage(0));
815 cvExtractSURF(&_img, pmask, &kp.seq, 0, storage, *(const CvSURFParams*)this, 0);
816 Seq<CvSURFPoint>::iterator it = kp.begin();
817 size_t i, n = kp.size();
819 for( i = 0; i < n; i++, ++it )
821 const CvSURFPoint& kpt = *it;
822 keypoints[i] = KeyPoint(kpt.pt, (float)kpt.size, kpt.dir,
823 kpt.hessian, getPointOctave(kpt, *this));
827 void SURF::operator()(const Mat& img, const Mat& mask,
828 vector<KeyPoint>& keypoints,
829 vector<float>& descriptors,
830 bool useProvidedKeypoints) const
832 CvMat _img = img, _mask, *pmask = 0;
834 pmask = &(_mask = mask);
835 MemStorage storage(cvCreateMemStorage(0));
839 if( useProvidedKeypoints )
841 kp = Seq<CvSURFPoint>(storage);
842 n = keypoints.size();
843 for( i = 0; i < n; i++ )
845 const KeyPoint& kpt = keypoints[i];
846 kp.push_back(cvSURFPoint(kpt.pt, 1, cvRound(kpt.size), kpt.angle, kpt.response));
850 cvExtractSURF(&_img, pmask, &kp.seq, &d, storage,
851 *(const CvSURFParams*)this, useProvidedKeypoints);
852 if( !useProvidedKeypoints )
854 Seq<CvSURFPoint>::iterator it = kp.begin();
855 size_t i, n = kp.size();
857 for( i = 0; i < n; i++, ++it )
859 const CvSURFPoint& kpt = *it;
860 keypoints[i] = KeyPoint(kpt.pt, (float)kpt.size, kpt.dir,
861 kpt.hessian, getPointOctave(kpt, *this));
864 descriptors.resize(d ? d->total*d->elem_size/sizeof(float) : 0);
866 cvCvtSeqToArray(d, &descriptors[0]);