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 matA = 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( &matA, &_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] );
391 enum { ORI_RADIUS = 6, ORI_WIN = 60, PATCH_SZ = 20 };
393 static const int ORI_SEARCH_INC;
394 static const float ORI_SIGMA;
395 static const float DESC_SIGMA;
397 SURFInvoker( const CvSURFParams* _params,
398 CvSeq* _keypoints, CvSeq* _descriptors,
399 const CvMat* _img, const CvMat* _sum,
400 const CvPoint* _apt, const float* _aptw,
401 int _nangle0, const float* _DW )
404 keypoints = _keypoints;
405 descriptors = _descriptors;
414 void operator()(const BlockedRange& range) const
416 /* X and Y gradient wavelet data */
417 const int NX=2, NY=2;
418 int dx_s[NX][5] = {{0, 0, 2, 4, -1}, {2, 0, 4, 4, 1}};
419 int dy_s[NY][5] = {{0, 0, 4, 2, 1}, {0, 2, 4, 4, -1}};
421 const int descriptor_size = params->extended ? 128 : 64;
423 const int max_ori_samples = (2*ORI_RADIUS+1)*(2*ORI_RADIUS+1);
424 float X[max_ori_samples], Y[max_ori_samples], angle[max_ori_samples];
425 uchar PATCH[PATCH_SZ+1][PATCH_SZ+1];
426 float DX[PATCH_SZ][PATCH_SZ], DY[PATCH_SZ][PATCH_SZ];
428 CvMat matX = cvMat(1, max_ori_samples, CV_32F, X);
429 CvMat matY = cvMat(1, max_ori_samples, CV_32F, Y);
430 CvMat _angle = cvMat(1, max_ori_samples, CV_32F, angle);
431 CvMat _patch = cvMat(PATCH_SZ+1, PATCH_SZ+1, CV_8U, PATCH);
433 int k, k1 = range.begin(), k2 = range.end();
436 for( k = k1; k < k2; k++ )
437 maxSize = std::max(maxSize, ((CvSURFPoint*)cvGetSeqElem( keypoints, k ))->size);
439 maxSize = cvCeil((PATCH_SZ+1)*maxSize*1.2f/9.0f);
440 Ptr<CvMat> winbuf = cvCreateMat( 1, maxSize*maxSize, CV_8U );
442 for( k = k1; k < k2; k++ )
444 const int* sum_ptr = sum->data.i;
445 int sum_cols = sum->cols;
446 int i, j, kk, x, y, nangle;
449 CvSurfHF dx_t[NX], dy_t[NY];
451 CvSURFPoint* kp = (CvSURFPoint*)cvGetSeqElem( keypoints, k );
453 CvPoint2D32f center = kp->pt;
455 /* The sampling intervals and wavelet sized for selecting an orientation
456 and building the keypoint descriptor are defined relative to 's' */
457 float s = (float)size*1.2f/9.0f;
459 /* To find the dominant orientation, the gradients in x and y are
460 sampled in a circle of radius 6s using wavelets of size 4s.
461 We ensure the gradient wavelet size is even to ensure the
462 wavelet pattern is balanced and symmetric around its center */
463 int grad_wav_size = 2*cvRound( 2*s );
464 if ( sum->rows < grad_wav_size || sum->cols < grad_wav_size )
466 /* when grad_wav_size is too big,
467 * the sampling of gradient will be meaningless
468 * mark keypoint for deletion. */
472 icvResizeHaarPattern( dx_s, dx_t, NX, 4, grad_wav_size, sum->cols );
473 icvResizeHaarPattern( dy_s, dy_t, NY, 4, grad_wav_size, sum->cols );
474 for( kk = 0, nangle = 0; kk < nangle0; kk++ )
478 x = cvRound( center.x + apt[kk].x*s - (float)(grad_wav_size-1)/2 );
479 y = cvRound( center.y + apt[kk].y*s - (float)(grad_wav_size-1)/2 );
480 if( (unsigned)y >= (unsigned)(sum->rows - grad_wav_size) ||
481 (unsigned)x >= (unsigned)(sum->cols - grad_wav_size) )
483 ptr = sum_ptr + x + y*sum_cols;
484 vx = icvCalcHaarPattern( ptr, dx_t, 2 );
485 vy = icvCalcHaarPattern( ptr, dy_t, 2 );
486 X[nangle] = vx*aptw[kk]; Y[nangle] = vy*aptw[kk];
491 /* No gradient could be sampled because the keypoint is too
492 * near too one or more of the sides of the image. As we
493 * therefore cannot find a dominant direction, we skip this
494 * keypoint and mark it for later deletion from the sequence. */
498 matX.cols = matY.cols = _angle.cols = nangle;
499 cvCartToPolar( &matX, &matY, 0, &_angle, 1 );
501 float bestx = 0, besty = 0, descriptor_mod = 0;
502 for( i = 0; i < 360; i += ORI_SEARCH_INC )
504 float sumx = 0, sumy = 0, temp_mod;
505 for( j = 0; j < nangle; j++ )
507 int d = std::abs(cvRound(angle[j]) - i);
508 if( d < ORI_WIN/2 || d > 360-ORI_WIN/2 )
514 temp_mod = sumx*sumx + sumy*sumy;
515 if( temp_mod > descriptor_mod )
517 descriptor_mod = temp_mod;
523 float descriptor_dir = cvFastArctan( besty, bestx );
524 kp->dir = descriptor_dir;
529 descriptor_dir *= (float)(CV_PI/180);
531 /* Extract a window of pixels around the keypoint of size 20s */
532 int win_size = (int)((PATCH_SZ+1)*s);
533 CV_Assert( winbuf->cols >= win_size*win_size );
535 CvMat win = cvMat(win_size, win_size, CV_8U, winbuf->data.ptr);
536 float sin_dir = sin(descriptor_dir);
537 float cos_dir = cos(descriptor_dir) ;
539 /* Subpixel interpolation version (slower). Subpixel not required since
540 the pixels will all get averaged when we scale down to 20 pixels */
542 float w[] = { cos_dir, sin_dir, center.x,
543 -sin_dir, cos_dir , center.y };
544 CvMat W = cvMat(2, 3, CV_32F, w);
545 cvGetQuadrangleSubPix( img, &win, &W );
548 /* Nearest neighbour version (faster) */
549 float win_offset = -(float)(win_size-1)/2;
550 float start_x = center.x + win_offset*cos_dir + win_offset*sin_dir;
551 float start_y = center.y - win_offset*sin_dir + win_offset*cos_dir;
552 uchar* WIN = win.data.ptr;
553 for( i = 0; i < win_size; i++, start_x += sin_dir, start_y += cos_dir )
555 float pixel_x = start_x;
556 float pixel_y = start_y;
557 for( j = 0; j < win_size; j++, pixel_x += cos_dir, pixel_y -= sin_dir )
559 int x = std::min(std::max(cvRound(pixel_x), 0), img->cols-1);
560 int y = std::min(std::max(cvRound(pixel_y), 0), img->rows-1);
561 WIN[i*win_size + j] = img->data.ptr[y*img->step + x];
565 /* Scale the window to size PATCH_SZ so each pixel's size is s. This
566 makes calculating the gradients with wavelets of size 2s easy */
567 cvResize( &win, &_patch, CV_INTER_AREA );
569 /* Calculate gradients in x and y with wavelets of size 2s */
570 for( i = 0; i < PATCH_SZ; i++ )
571 for( j = 0; j < PATCH_SZ; j++ )
573 float dw = DW[i*PATCH_SZ + j];
574 float vx = (PATCH[i][j+1] - PATCH[i][j] + PATCH[i+1][j+1] - PATCH[i+1][j])*dw;
575 float vy = (PATCH[i+1][j] - PATCH[i][j] + PATCH[i+1][j+1] - PATCH[i][j+1])*dw;
580 /* Construct the descriptor */
581 vec = (float*)cvGetSeqElem( descriptors, k );
582 for( kk = 0; kk < (int)(descriptors->elem_size/sizeof(vec[0])); kk++ )
584 double square_mag = 0;
585 if( params->extended )
587 /* 128-bin descriptor */
588 for( i = 0; i < 4; i++ )
589 for( j = 0; j < 4; j++ )
591 for( y = i*5; y < i*5+5; y++ )
593 for( x = j*5; x < j*5+5; x++ )
595 float tx = DX[y][x], ty = DY[y][x];
599 vec[1] += (float)fabs(tx);
602 vec[3] += (float)fabs(tx);
607 vec[5] += (float)fabs(ty);
610 vec[7] += (float)fabs(ty);
614 for( kk = 0; kk < 8; kk++ )
615 square_mag += vec[kk]*vec[kk];
621 /* 64-bin descriptor */
622 for( i = 0; i < 4; i++ )
623 for( j = 0; j < 4; j++ )
625 for( y = i*5; y < i*5+5; y++ )
627 for( x = j*5; x < j*5+5; x++ )
629 float tx = DX[y][x], ty = DY[y][x];
630 vec[0] += tx; vec[1] += ty;
631 vec[2] += (float)fabs(tx); vec[3] += (float)fabs(ty);
634 for( kk = 0; kk < 4; kk++ )
635 square_mag += vec[kk]*vec[kk];
640 /* unit vector is essential for contrast invariance */
641 vec = (float*)cvGetSeqElem( descriptors, k );
642 double scale = 1./(sqrt(square_mag) + DBL_EPSILON);
643 for( kk = 0; kk < descriptor_size; kk++ )
644 vec[kk] = (float)(vec[kk]*scale);
648 const CvSURFParams* params;
659 const int SURFInvoker::ORI_SEARCH_INC = 5;
660 const float SURFInvoker::ORI_SIGMA = 2.5f;
661 const float SURFInvoker::DESC_SIGMA = 3.3f;
667 cvExtractSURF( const CvArr* _img, const CvArr* _mask,
668 CvSeq** _keypoints, CvSeq** _descriptors,
669 CvMemStorage* storage, CvSURFParams params,
670 int useProvidedKeyPts)
672 const int ORI_RADIUS = cv::SURFInvoker::ORI_RADIUS;
673 const int ORI_SIGMA = cv::SURFInvoker::ORI_SIGMA;
674 const float DESC_SIGMA = cv::SURFInvoker::DESC_SIGMA;
676 CvMat *sum = 0, *mask1 = 0, *mask_sum = 0;
678 if( _keypoints && !useProvidedKeyPts ) // If useProvidedKeyPts!=0 we'll use current contents of "*_keypoints"
683 CvSeq *keypoints, *descriptors = 0;
684 CvMat imghdr, *img = cvGetMat(_img, &imghdr);
685 CvMat maskhdr, *mask = _mask ? cvGetMat(_mask, &maskhdr) : 0;
687 const int max_ori_samples = (2*ORI_RADIUS+1)*(2*ORI_RADIUS+1);
688 int descriptor_size = params.extended ? 128 : 64;
689 const int descriptor_data_type = CV_32F;
690 const int PATCH_SZ = 20;
691 float DW[PATCH_SZ][PATCH_SZ];
692 CvMat _DW = cvMat(PATCH_SZ, PATCH_SZ, CV_32F, DW);
693 CvPoint apt[max_ori_samples];
694 float aptw[max_ori_samples];
695 int i, j, nangle0 = 0, N;
698 CV_Assert(CV_MAT_TYPE(img->type) == CV_8UC1);
699 CV_Assert(mask == 0 || (CV_ARE_SIZES_EQ(img,mask) && CV_MAT_TYPE(mask->type) == CV_8UC1));
700 CV_Assert(storage != 0);
701 CV_Assert(params.hessianThreshold >= 0);
702 CV_Assert(params.nOctaves > 0);
703 CV_Assert(params.nOctaveLayers > 0);
705 sum = cvCreateMat( img->rows+1, img->cols+1, CV_32SC1 );
706 cvIntegral( img, sum );
708 // Compute keypoints only if we are not asked for evaluating the descriptors are some given locations:
709 if (!useProvidedKeyPts)
713 mask1 = cvCreateMat( img->height, img->width, CV_8UC1 );
714 mask_sum = cvCreateMat( img->height+1, img->width+1, CV_32SC1 );
715 cvMinS( mask, 1, mask1 );
716 cvIntegral( mask1, mask_sum );
718 keypoints = icvFastHessianDetector( sum, mask_sum, storage, ¶ms );
722 CV_Assert(useProvidedKeyPts && (_keypoints != 0) && (*_keypoints != 0));
723 keypoints = *_keypoints;
726 N = keypoints->total;
729 descriptors = cvCreateSeq( 0, sizeof(CvSeq),
730 descriptor_size*CV_ELEM_SIZE(descriptor_data_type), storage );
731 cvSeqPushMulti( descriptors, 0, N );
734 /* Coordinates and weights of samples used to calculate orientation */
735 cv::Mat matG = cv::getGaussianKernel( 2*ORI_RADIUS+1, ORI_SIGMA, CV_32F );
736 const float* G = (const float*)matG.data;
738 for( i = -ORI_RADIUS; i <= ORI_RADIUS; i++ )
740 for( j = -ORI_RADIUS; j <= ORI_RADIUS; j++ )
742 if( i*i + j*j <= ORI_RADIUS*ORI_RADIUS )
744 apt[nangle0] = cvPoint(j,i);
745 aptw[nangle0++] = G[i+ORI_RADIUS]*G[j+ORI_RADIUS];
750 /* Gaussian used to weight descriptor samples */
751 double c2 = 1./(DESC_SIGMA*DESC_SIGMA*2);
753 for( i = 0; i < PATCH_SZ; i++ )
755 for( j = 0; j < PATCH_SZ; j++ )
757 double x = j - (float)(PATCH_SZ-1)/2, y = i - (float)(PATCH_SZ-1)/2;
758 double val = exp(-(x*x+y*y)*c2);
759 DW[i][j] = (float)val;
763 cvScale( &_DW, &_DW, 1./gs );
765 cv::parallel_for(cv::BlockedRange(0, N),
766 cv::SURFInvoker(¶ms, keypoints, descriptors, img, sum,
767 apt, aptw, nangle0, &DW[0][0]));
768 //cv::SURFInvoker(¶ms, keypoints, descriptors, img, sum,
769 // apt, aptw, nangle0, &DW[0][0])(cv::BlockedRange(0, N));
771 /* remove keypoints that were marked for deletion */
772 for ( i = 0; i < N; i++ )
774 CvSURFPoint* kp = (CvSURFPoint*)cvGetSeqElem( keypoints, i );
775 if ( kp->size == -1 )
777 cvSeqRemove( keypoints, i );
779 cvSeqRemove( descriptors, i );
785 if( _keypoints && !useProvidedKeyPts )
786 *_keypoints = keypoints;
788 *_descriptors = descriptors;
790 cvReleaseMat( &sum );
791 if (mask1) cvReleaseMat( &mask1 );
792 if (mask_sum) cvReleaseMat( &mask_sum );
801 hessianThreshold = 100;
807 SURF::SURF(double _threshold, int _nOctaves, int _nOctaveLayers, bool _extended)
809 hessianThreshold = _threshold;
810 extended = _extended;
811 nOctaves = _nOctaves;
812 nOctaveLayers = _nOctaveLayers;
815 int SURF::descriptorSize() const { return extended ? 128 : 64; }
818 static int getPointOctave(const CvSURFPoint& kpt, const CvSURFParams& params)
820 int octave = 0, layer = 0, best_octave = 0;
821 float min_diff = FLT_MAX;
822 for( octave = 1; octave < params.nOctaves; octave++ )
823 for( layer = 0; layer < params.nOctaveLayers; layer++ )
825 float diff = std::abs(kpt.size - (float)((HAAR_SIZE0 + HAAR_SIZE_INC*layer) << octave));
826 if( min_diff > diff )
829 best_octave = octave;
838 void SURF::operator()(const Mat& img, const Mat& mask,
839 vector<KeyPoint>& keypoints) const
841 CvMat _img = img, _mask, *pmask = 0;
843 pmask = &(_mask = mask);
844 MemStorage storage(cvCreateMemStorage(0));
846 cvExtractSURF(&_img, pmask, &kp.seq, 0, storage, *(const CvSURFParams*)this, 0);
847 Seq<CvSURFPoint>::iterator it = kp.begin();
848 size_t i, n = kp.size();
850 for( i = 0; i < n; i++, ++it )
852 const CvSURFPoint& kpt = *it;
853 keypoints[i] = KeyPoint(kpt.pt, (float)kpt.size, kpt.dir,
854 kpt.hessian, getPointOctave(kpt, *this));
858 void SURF::operator()(const Mat& img, const Mat& mask,
859 vector<KeyPoint>& keypoints,
860 vector<float>& descriptors,
861 bool useProvidedKeypoints) const
863 CvMat _img = img, _mask, *pmask = 0;
865 pmask = &(_mask = mask);
866 MemStorage storage(cvCreateMemStorage(0));
870 if( useProvidedKeypoints )
872 kp = cvCreateSeq( 0, sizeof(CvSeq), sizeof(CvSURFPoint), storage);
873 n = keypoints.size();
874 for( i = 0; i < n; i++ )
876 const KeyPoint& kpt = keypoints[i];
877 kp.push_back(cvSURFPoint(kpt.pt, 1, cvRound(kpt.size), kpt.angle, kpt.response));
881 cvExtractSURF(&_img, pmask, &kp.seq, &d, storage,
882 *(const CvSURFParams*)this, useProvidedKeypoints);
883 if( !useProvidedKeypoints )
885 Seq<CvSURFPoint>::iterator it = kp.begin();
886 size_t i, n = kp.size();
888 for( i = 0; i < n; i++, ++it )
890 const CvSURFPoint& kpt = *it;
891 keypoints[i] = KeyPoint(kpt.pt, (float)kpt.size, kpt.dir,
892 kpt.hessian, getPointOctave(kpt, *this));
895 descriptors.resize(d ? d->total*d->elem_size/sizeof(float) : 0);
896 if(descriptors.size() != 0)
897 cvCvtSeqToArray(d, &descriptors[0]);