MVE - Multi-View Environment mve-devel
Loading...
Searching...
No Matches
surf.cc
Go to the documentation of this file.
1/*
2 * Copyright (C) 2015, Simon Fuhrmann
3 * TU Darmstadt - Graphics, Capture and Massively Parallel Computing
4 * All rights reserved.
5 *
6 * This software may be modified and distributed under the terms
7 * of the BSD 3-Clause license. See the LICENSE.txt file for details.
8 */
9
10#include <iostream>
11
12#include "util/timer.h"
13#include "math/functions.h"
14#include "math/vector.h"
15#include "math/matrix.h"
16#include "math/matrix_tools.h"
17#include "mve/image.h"
18#include "mve/image_tools.h"
19#include "mve/image_drawing.h"
20#include "sfm/defines.h"
21#include "sfm/surf.h"
22
24
25namespace
26{
27 /* Kernel sizes per octave (given in 1/3 of the full size). */
28 int const kernel_sizes[4][4] =
29 {
30 { 3, 5, 7, 9 }, // 9 15 21 27
31 { 5, 9, 13, 17 }, // 15 27 39 51
32 { 9, 17, 25, 33 }, // 27 51 75 99
33 { 17, 33, 49, 65 } // 51 99 147 195
34 };
35} // namespace
36
37/* ---------------------------------------------------------------- */
38
39Surf::Surf (Options const& options)
40 : options(options)
41{
42 if (this->options.debug_output)
43 this->options.verbose_output = true;
44}
45
46/* ---------------------------------------------------------------- */
47
48void
50{
51 this->keypoints.clear();
52 this->descriptors.clear();
53 this->octaves.clear();
54
55 util::WallTimer timer, total_timer;
56
57 /* Compute Hessian response maps and find SS maxima (SURF 3.3). */
58 if (this->options.verbose_output)
59 {
60 std::cout << "SURF: Creating 4 octaves (0 to 4)..." << std::endl;
61 }
62 timer.reset();
63 this->create_octaves();
64 if (this->options.debug_output)
65 {
66 std::cout << "SURF: Creating octaves took "
67 << timer.get_elapsed() << " ms." << std::endl;
68 }
69
70 /* Detect local extrema in the SS of Hessian response maps. */
71 if (this->options.debug_output)
72 {
73 std::cout << "SURF: Detecting local extrema..." << std::endl;
74 }
75 timer.reset();
76 this->extrema_detection();
77 if (this->options.debug_output)
78 {
79 std::cout << "SURF: Extrema detection took "
80 << timer.get_elapsed() << " ms." << std::endl;
81 }
82
83 /* Sub-pixel keypoint localization and filtering of weak keypoints. */
84 if (this->options.debug_output)
85 {
86 std::cout << "SURF: Localizing and filtering keypoints..." << std::endl;
87 }
88 timer.reset();
90 this->octaves.clear();
91 if (this->options.debug_output)
92 {
93 std::cout << "SURF: Localization and filtering took "
94 << timer.get_elapsed() << " ms." << std::endl;
95 }
96
97 /* Compute the SURF descriptor for the keypoint location. */
98 if (this->options.verbose_output)
99 {
100 std::cout << "SURF: Generating keypoint descriptors..." << std::endl;
101 }
102 timer.reset();
103 this->descriptor_assignment();
104 if (this->options.debug_output)
105 {
106 std::cout << "SURF: Generated " << this->descriptors.size()
107 << " descriptors, took " << timer.get_elapsed() << "ms."
108 << std::endl;
109 }
110 if (this->options.verbose_output)
111 {
112 std::cout << "SURF: Generated " << this->descriptors.size()
113 << " descriptors from " << this->keypoints.size() << " keypoints,"
114 << " took " << total_timer.get_elapsed() << "ms." << std::endl;
115 }
116
117 /* Cleanup. */
118 this->sat.reset();
119}
120
121/* ---------------------------------------------------------------- */
122
123void
125{
126 /* Desaturate input image. */
127 if (image->channels() != 1 && image->channels() != 3)
128 throw std::invalid_argument("Gray or color image expected");
129 if (image->channels() == 3)
130 image = mve::image::desaturate<uint8_t>(image,
132
133 /* Build summed area table (SAT). */
134 //util::WallTimer timer;
135 this->sat = mve::image::integral_image<uint8_t,SatType>(image);
136 //std::cout << "Creating SAT image took "
137 // << timer.get_elapsed() << " ms." << std::endl;
138}
139
140/* ---------------------------------------------------------------- */
141
142void
144{
145 /* Prepare octaves. */
146 this->octaves.resize(4);
147 for (int i = 0; i < 4; ++i)
148 this->octaves[i].imgs.resize(4);
149
150 /* Create octaves. */
151 for (int octave = 0; octave < 4; ++octave)
152 for (int sample = 0; sample < 4; ++sample)
153 this->create_response_map(octave, sample);
154}
155
156/* ---------------------------------------------------------------- */
157
158void
160{
161 /*
162 * In order to create the Hessian response map for filter size 'fs',
163 * we need to convolute the image with second order Gaussian partial
164 * derivative filters Dxx, Dyy and Dxy and compute the response
165 * as det(H) = Dxx * Dyy - w * Dxy * Dxy, where w = 0.81.
166 * So we need to compute Dxx, Dyy and Dxy with filter size 'fs'.
167 * Note: filter size 'fs' is defined as filter width / 3.
168 * For details, see SURF 3.2.
169 */
170 typedef Octave::RespType RespType;
171
172 /* Filter size. The actual kernel side length is 3 * fs. */
173 int const fs = kernel_sizes[o][k];
174 /* The sample spacing for the octaves. */
175 int const step = math::fastpow(2, o);
176 /* Weight to balance between real gaussian kernel and approximated one. */
177 RespType const weight = 0.912; // See SURF 3.3 (4).
178 /* NormalizationKernel width is 3 * fs, height is 2 * fs - 1. */
179 RespType const inv_karea = 1.0 / (fs * (2 * fs - 1));
180
181 /* Original dimensions and octave dimensions. */
182 int const w = this->sat->width();
183 int const h = this->sat->height();
184 int ow = w;
185 int oh = h;
186 for (int i = 0; i < o; ++i)
187 {
188 ow = (ow + 1) >> 1;
189 oh = (oh + 1) >> 1;
190 }
191
192 /* Generate the response maps. */
194 int const border = fs + fs / 2 + 1;
195 for (int y = 0, i = 0; y < h; y += step)
196 for (int x = 0; x < w; x += step, ++i)
197 {
198 if (x < border || x + border >= w || y < border || y + border >= h)
199 {
200 img->at(i) = 0.0f;
201 continue;
202 }
203
204 SatType dxx = this->filter_dxx(fs, x, y);
205 SatType dyy = this->filter_dyy(fs, x, y);
206 SatType dxy = this->filter_dxy(fs, x, y);
207 RespType dxx_t = static_cast<RespType>(dxx) * inv_karea;
208 RespType dyy_t = static_cast<RespType>(dyy) * inv_karea;
209 RespType dxy_t = static_cast<RespType>(dxy) * inv_karea;
210 /* Compute the determinant of the hessian. */
211 img->at(i) = dxx_t * dyy_t - weight * dxy_t * dxy_t;
212 /* The laplacian can be computed as dxx_t + dyy_t. */
213 // float laplacian = dxx_t + dyy_t;
214 }
215
216 this->octaves[o].imgs[k] = img;
217}
218
219/* ---------------------------------------------------------------- */
220
222Surf::filter_dxx (int fs, int x, int y)
223{
224 int const w = this->sat->width();
225 int const fs2 = fs / 2;
226 int const row1 = (x - fs - fs2 - 1) + w * (y - fs);
227 int const row2 = row1 + w * (fs + fs - 1);
228
229 Surf::SatType const v0 = this->sat->at(row1);
230 Surf::SatType const v1 = this->sat->at(row1 + fs);
231 Surf::SatType const v2 = this->sat->at(row1 + 2 * fs);
232 Surf::SatType const v3 = this->sat->at(row1 + 3 * fs);
233 Surf::SatType const v4 = this->sat->at(row2);
234 Surf::SatType const v5 = this->sat->at(row2 + fs);
235 Surf::SatType const v6 = this->sat->at(row2 + 2 * fs);
236 Surf::SatType const v7 = this->sat->at(row2 + 3 * fs);
237
238 Surf::SatType ret = 0;
239 ret += v5 + v0 - v4 - v1;
240 ret -= 2 * (v6 + v1 - v5 - v2);
241 ret += v7 + v2 - v6 - v3;
242 return ret;
243}
244
246Surf::filter_dyy (int fs, int x, int y)
247{
248 int const w = this->sat->width();
249 int const fs2 = fs / 2;
250 int const row1 = (x - fs) + w * (y - fs - fs2 - 1);
251 int const row2 = row1 + w * fs;
252 int const row3 = row2 + w * fs;
253 int const row4 = row3 + w * fs;
254
255 Surf::SatType const v0 = this->sat->at(row1);
256 Surf::SatType const v1 = this->sat->at(row2);
257 Surf::SatType const v2 = this->sat->at(row3);
258 Surf::SatType const v3 = this->sat->at(row4);
259 Surf::SatType const v4 = this->sat->at(row1 + fs + fs - 1);
260 Surf::SatType const v5 = this->sat->at(row2 + fs + fs - 1);
261 Surf::SatType const v6 = this->sat->at(row3 + fs + fs - 1);
262 Surf::SatType const v7 = this->sat->at(row4 + fs + fs - 1);
263
264 Surf::SatType ret = 0;
265 ret += v5 + v0 - v1 - v4;
266 ret -= 2 * (v6 + v1 - v2 - v5);
267 ret += v7 + v2 - v3 - v6;
268 return ret;
269}
270
272Surf::filter_dxy (int fs, int x, int y)
273{
274 int const w = this->sat->width();
275 int const row1 = (x - fs - 1) + w * (y - fs - 1);
276 int const row2 = row1 + w * fs;
277 int const row3 = row2 + w;
278 int const row4 = row3 + w * fs;
279
280 Surf::SatType ret = 0;
281 Surf::SatType v0 = this->sat->at(row1);
282 Surf::SatType v1 = this->sat->at(row1 + fs);
283 Surf::SatType v2 = this->sat->at(row2);
284 Surf::SatType v3 = this->sat->at(row2 + fs);
285 ret += v3 + v0 - v2 - v1;
286
287 v0 = this->sat->at(row1 + fs + 1);
288 v1 = this->sat->at(row1 + fs + fs + 1);
289 v2 = this->sat->at(row2 + fs + 1);
290 v3 = this->sat->at(row2 + fs + fs + 1);
291 ret -= v3 + v0 - v2 - v1;
292
293 v0 = this->sat->at(row3);
294 v1 = this->sat->at(row3 + fs);
295 v2 = this->sat->at(row4);
296 v3 = this->sat->at(row4 + fs);
297 ret -= v3 + v0 - v2 - v1;
298
299 v0 = this->sat->at(row3 + fs + 1);
300 v1 = this->sat->at(row3 + fs + fs + 1);
301 v2 = this->sat->at(row4 + fs + 1);
302 v3 = this->sat->at(row4 + fs + fs + 1);
303 ret += v3 + v0 - v2 - v1;
304
305 return ret;
306}
307
308/* ---------------------------------------------------------------- */
309
310void
312{
313 /*
314 * At this stage each octave contains 4 scale space samples and local
315 * maxima/minima in the approximated DoG function need to be found.
316 * To this end, a simple non-maximum suppression technique is applied.
317 */
318 for (std::size_t o = 0; o < this->octaves.size(); ++o)
319 {
320 int const width = this->octaves[o].imgs[0]->width();
321 int const height = this->octaves[o].imgs[0]->height();
322 for (std::size_t s = 1; s < 3; ++s)
323 {
324 Octave::RespImage::ConstPtr ss = this->octaves[o].imgs[s];
325 Octave::RespType const* ptr = &ss->at(0);
326 for (int y = 0; y < height; ++y, ptr += width)
327 for (int x = 1; x + 1 < width; ++x)
328 if (ptr[x] > ptr[x-1] && ptr[x] > ptr[x + 1])
329 this->check_maximum(o, s, x, y);
330 }
331 }
332}
333
334/* ---------------------------------------------------------------- */
335
336void
337Surf::check_maximum(int o, int s, int x, int y)
338{
339 /*
340 * Assumes that given coordinates are within bounds and a 1 pixel
341 * boundary for comarisons in x, y and s direction.
342 */
343 int const w = this->octaves[o].imgs[s]->width();
344 int off1[8] = { -w - 1, -w, -w + 1, -1, +1, w - 1, w, w + 1 };
345 int off2[9] = { -w - 1, -w, -w + 1, -1, 0, +1, w - 1, w, w + 1 };
346
347 Octave::RespType const* ptr;
348 int const off = x + y * w;
349
350 /* Perform NMS check on the candidate sample first. */
351 ptr = this->octaves[o].imgs[s]->get_data_pointer() + off;
352 Octave::RespType value = ptr[0];
353 for (int i = 0; i < 8; ++i)
354 if (ptr[off1[i]] >= value)
355 return;
356 /* Perform NMS check on the above sample. */
357 ptr = this->octaves[o].imgs[s - 1]->get_data_pointer() + off;
358 for (int i = 0; i < 9; ++i)
359 if (ptr[off2[i]] >= value)
360 return;
361 /* Perform NMS check on the below sample. */
362 ptr = this->octaves[o].imgs[s + 1]->get_data_pointer() + off;
363 for (int i = 0; i < 9; ++i)
364 if (ptr[off2[i]] >= value)
365 return;
366
367 /* Seems like we found a keypoint. */
368 Keypoint kp;
369 kp.octave = o;
370 kp.sample = static_cast<float>(s);
371 kp.x = static_cast<float>(x);
372 kp.y = static_cast<float>(y);
373 this->keypoints.push_back(kp);
374}
375
376/* ---------------------------------------------------------------- */
377
378void
380{
381 Keypoints filtered_keypoints;
382 filtered_keypoints.reserve(this->keypoints.size());
383 for (std::size_t i = 0; i < this->keypoints.size(); ++i)
384 {
385 Keypoint kp = this->keypoints[i];
386 bool solve_good = this->keypoint_localization(&kp);
387 if (!solve_good)
388 continue;
389 filtered_keypoints.push_back(kp);
390 }
391 std::swap(this->keypoints, filtered_keypoints);
392}
393
394/* ---------------------------------------------------------------- */
395
396bool
398{
399 int const sample = static_cast<int>(kp->sample);
400 int const w = this->octaves[kp->octave].imgs[sample]->width();
401 //int const h = this->octaves[kp->octave].imgs[sample]->height();
402 int off = static_cast<int>(kp->x) + static_cast<int>(kp->y) * w;
403 Octave::RespType const *s0, *s1, *s2;
404 s0 = this->octaves[kp->octave].imgs[sample - 1]->get_data_pointer() + off;
405 s1 = this->octaves[kp->octave].imgs[sample]->get_data_pointer() + off;
406 s2 = this->octaves[kp->octave].imgs[sample + 1]->get_data_pointer() + off;
407
408 float N9[3][9] =
409 {
410 { s0[-1-w], s0[-w], s0[1-w], s0[-1], s0[0], s0[1], s0[w-1], s0[w], s0[1+w] },
411 { s1[-1-w], s1[-w], s1[1-w], s1[-1], s1[0], s1[1], s1[w-1], s1[w], s1[1+w] },
412 { s2[-1-w], s2[-w], s2[1-w], s2[-1], s2[0], s2[1], s2[w-1], s2[w], s2[1+w] }
413 };
414
415 /* Switch to processing in double. Determinant can be very large. */
416 math::Vec3d vec_b;
417 math::Matrix3d mat_a;
418 vec_b[0] = -(N9[1][5] - N9[1][3]) * 0.5; // 1st deriv x.
419 vec_b[1] = -(N9[1][7] - N9[1][1]) * 0.5; // 1st deriv y.
420 vec_b[2] = -(N9[2][4] - N9[0][4]) * 0.5; // 1st deriv s.
421
422 mat_a[0] = N9[1][3] - 2.0 * N9[1][4] + N9[1][5]; // xx
423 mat_a[1] = (N9[1][8] - N9[1][6] - N9[1][2] + N9[1][0]) * 0.25; // xy
424 mat_a[2] = (N9[2][5] - N9[2][3] - N9[0][5] + N9[0][3]) * 0.25; // xs
425 mat_a[3] = mat_a[1]; // yx
426 mat_a[4] = N9[1][1] - 2.0 * N9[1][4] + N9[1][7]; // yy
427 mat_a[5] = (N9[2][7] - N9[2][1] - N9[0][7] + N9[0][1]) * 0.25; // ys
428 mat_a[6] = mat_a[2]; // sx
429 mat_a[7] = mat_a[5]; // sy
430 mat_a[8] = N9[0][4] - 2.0 * N9[1][4] + N9[2][4]; // ss
431
432 /* Compute determinant to detect singular matrix. */
433 double det_a = math::matrix_determinant(mat_a);
434 //std::cout << "Determinant: " << det_a << std::endl;
435 if (MATH_EPSILON_EQ(det_a, 0.0f, 1.0e-5f))
436 {
437 //std::cout << "Rejection keypoint: det = " << det_a << std::endl;
438 return false;
439 }
440
441 /* Invert the matrix to get the accurate keypoint. */
442 mat_a = math::matrix_inverse(mat_a, det_a);
443 math::Vec3d vec_x = mat_a * vec_b;
444
445 /* Reject keypoint if location is too far off original point. */
446 if (vec_x.maximum() > 0.5 || vec_x.minimum() < -0.5f)
447 return false;
448
449 /* Compute actual DoG value at accurate keypoint x. */
450 float dog_value = N9[1][4] - 0.5 * vec_b.dot(vec_x);
451 if (dog_value < this->options.contrast_threshold)
452 {
453 //std::cout << "Rejection keypoint: Contrast Thres" << std::endl;
454 return false;
455 }
456
457 /* Update keypoint. */
458 float sampling = static_cast<float>(math::fastpow(2, kp->octave));
459 kp->x = (kp->x + vec_x[0]) * sampling;
460 kp->y = (kp->y + vec_x[1]) * sampling;
461 // OpenSURF code: ipt.scale = static_cast<float>((0.1333f) * (m->filter + xi * filterStep));
462 kp->sample += vec_x[2]; // Is this correct? No! FIXME!
463
464 if (kp->x < 0.0f || kp->x + 1.0f > this->sat->width()
465 || kp->y < 0.0f || kp->y + 1.0f > this->sat->height())
466 {
467 //std::cout << "SURF: Rejection keypoint: OOB "
468 // << kp->x << " " << kp->y << " / " << vec_x[0] << " "
469 // << vec_x[1] << std::endl;
470 return false;
471 }
472
473 return true;
474}
475
476/* ---------------------------------------------------------------- */
477
478void
480{
481 this->descriptors.clear();
482 this->descriptors.reserve(keypoints.size());
483 for (std::size_t i = 0; i < this->keypoints.size(); ++i)
484 {
485 Keypoint const& kp = this->keypoints[i];
486
487 /* Copy over the basic information to the descriptor. */
488 Descriptor descr;
489 descr.x = kp.x;
490 descr.y = kp.y;
491
492 /*
493 * The scale is obtained from the filter size. The smallest filter in
494 * SURF has size 9 and corresponds to a scale of 1.2. Thus, the
495 * scale of a filter with size X has a scale of X * 1.2 / 9.
496 */
497 int sample = static_cast<int>(kp.sample + 0.5f);
498 descr.scale = 3 * kernel_sizes[kp.octave][sample] * 1.2f / 9.0f;
499
500 /* Find the orientation of the keypoint. */
501 if (!this->descriptor_orientation(&descr))
502 continue;
503
504 /* Compute descriptor relative to orientation. */
505 if (!this->descriptor_computation(&descr,
506 this->options.use_upright_descriptor))
507 continue;
508
509 this->descriptors.push_back(descr);
510 }
511}
512
513/* ---------------------------------------------------------------- */
514
515bool
517{
518 int const descr_x = static_cast<int>(descr->x + 0.5f);
519 int const descr_y = static_cast<int>(descr->y + 0.5f);
520 int const descr_scale = static_cast<int>(descr->scale);
521 int const width = this->sat->width();
522 int const height = this->sat->height();
523
524 /*
525 * Pre-computed gaussian weights for the (circular) window. The gaussian
526 * is computed as exp((dx^2 + dy^2) / (2 * sigma^2) with sigma = 2.5.
527 */
528 float const gaussian[109] = { 0.0658748, 0.0982736, 0.12493, 0.135335,
529 0.12493, 0.0982736, 0.0658748, 0.0773047, 0.135335, 0.201897, 0.256661,
530 0.278037, 0.256661, 0.201897, 0.135335, 0.0773047, 0.0658748, 0.135335,
531 0.236928, 0.353455, 0.449329, 0.486752, 0.449329, 0.353455, 0.236928,
532 0.135335, 0.0658748, 0.0982736, 0.201897, 0.353455, 0.527292, 0.67032,
533 0.726149, 0.67032, 0.527292, 0.353455, 0.201897, 0.0982736, 0.12493,
534 0.256661, 0.449329, 0.67032, 0.852144, 0.923116, 0.852144, 0.67032,
535 0.449329, 0.256661, 0.12493, 0.135335, 0.278037, 0.486752, 0.726149,
536 0.923116, 1, 0.923116, 0.726149, 0.486752, 0.278037, 0.135335, 0.12493,
537 0.256661, 0.449329, 0.67032, 0.852144, 0.923116, 0.852144, 0.67032,
538 0.449329, 0.256661, 0.12493, 0.0982736, 0.201897, 0.353455, 0.527292,
539 0.67032, 0.726149, 0.67032, 0.527292, 0.353455, 0.201897, 0.0982736,
540 0.0658748, 0.135335, 0.236928, 0.353455, 0.449329, 0.486752, 0.449329,
541 0.353455, 0.236928, 0.135335, 0.0658748, 0.0773047, 0.135335, 0.201897,
542 0.256661, 0.278037, 0.256661, 0.201897, 0.135335, 0.0773047, 0.0658748,
543 0.0982736, 0.12493, 0.135335, 0.12493, 0.0982736, 0.0658748 };
544
545 /*
546 * At least a 12 * scale pixel kernel for the support circle is needed.
547 * Additionally, computing the Haar Wavelet response uses a kernel size
548 * of 4 * scale pixel. That makes a total of (6 + 2) * scale pixel on
549 * either side of the kernel spacing. One additional pixel is needed
550 * to the upper-left side to simplify integral image access.
551 */
552 int const spacing = 8 * descr_scale + 1;
553 if (descr_x < spacing || descr_y < spacing
554 || descr_x + spacing >= width || descr_y + spacing >= height)
555 return false;
556
557 /* The number of samples is constant, depending on radius. */
558#define NUM_SAMPLES 109
559 float dx[NUM_SAMPLES], dy[NUM_SAMPLES];
560
561 /*
562 * Iterate over the pixels of a circle with radius 6 * scale
563 * and compute Haar Wavelet responses in x- and y-direction.
564 */
565 for (int ry = -5, index = 0; ry <= 5; ++ry)
566 for (int rx = -5; rx <= 5; ++rx)
567 {
568 if (rx * rx + ry * ry >= 36)
569 continue;
570 this->filter_dx_dy(descr_x + rx * descr_scale,
571 descr_y + ry * descr_scale, 2 * descr_scale,
572 dx + index, dy + index);
573 dx[index] *= gaussian[index];
574 dy[index] *= gaussian[index];
575 index += 1;
576 }
577
578 /*
579 * Iterate in a sliding window over the 109 extracted samples
580 * in order to find a dominant orientation for the keypoint.
581 */
582 double const window_increment = MATH_PI / 8.0;
583 double const window_halfsize = MATH_PI / 6.0;
584 double best_dx = 0.0, best_dy = 0.0;
585 double best_length = 0.0f;
586 for (double deg = -MATH_PI; deg < MATH_PI; deg += window_increment)
587 {
588 /* Evaluate samples for this window. */
589 double const win_start = deg - window_halfsize;
590 double const win_end = deg + window_halfsize;
591 double sum_dx = 0.0f, sum_dy = 0.0f;
592 for (int i = 0; i < NUM_SAMPLES; ++i)
593 {
594 double const sample_deg = std::atan2(dy[i], dx[i]);
595 double const sample_deg_p = sample_deg + 2.0 * MATH_PI;
596 double const sample_deg_m = sample_deg - 2.0 * MATH_PI;
597 if ((sample_deg > win_start && sample_deg < win_end)
598 || (sample_deg_p > win_start && sample_deg_p < win_end)
599 || (sample_deg_m > win_start && sample_deg_m < win_end))
600 {
601 sum_dx += dx[i];
602 sum_dy += dy[i];
603 }
604 }
605
606 /* Total vector length of dx/dy sums defines dominance. */
607 double const length = sum_dx * sum_dx + sum_dy * sum_dy;
608 if (length > best_length)
609 {
610 best_dx = sum_dx;
611 best_dy = sum_dy;
612 best_length = length;
613 }
614 }
615
616 /* TODO: Threshold on the vector length? */
617 descr->orientation = std::atan2(best_dy, best_dx);
618 return true;
619}
620
621/* ---------------------------------------------------------------- */
622
623void
624Surf::filter_dx_dy (int x, int y, int fs, float* dx, float* dy)
625{
626 int const width = this->sat->width();
627 SatType const* ptr = this->sat->begin();
628
629 /*
630 * To have a center pixel, filter size needs to be odd.
631 * To ensure symmetry in the filters, we include the center row in
632 * both sides of the dy filter, and the center column in both sides
633 * of the dx filter, which cancels them out. However, this costs
634 * four additional lookups (12 instead of 8).
635 */
636 SatType const* iter = ptr + (x - fs - 1) + (y - fs - 1) * width;
637 SatType x1 = *iter; iter += fs;
638 SatType x2 = *iter; iter += 1;
639 SatType x3 = *iter; iter += fs;
640 SatType x4 = *iter; iter += (width - 1) * (2 * fs + 1);
641
642 SatType x5 = *iter; iter += fs;
643 SatType x6 = *iter; iter += 1;
644 SatType x7 = *iter; iter += fs;
645 SatType x8 = *iter;
646
647 iter = ptr + (x - fs - 1) + (y - 1) * width;
648 SatType y1 = *iter; iter += 2 * fs + 1;
649 SatType y2 = *iter; iter += width - 2 * fs - 1;
650 SatType y3 = *iter; iter += 2 * fs + 1;
651 SatType y4 = *iter;
652
653 /*
654 * Normalize filter by size "(2 * fs + 1) * fs" and normalize discrete
655 * derivative with distance between the Wavelet box centers "fs + 1".
656 */
657 float norm = static_cast<float>((2 * fs + 1) * fs * (fs + 1));
658 *dx = static_cast<float>((x8 + x2 - x4 - x6) - (x7 + x1 - x3 - x5)) / norm;
659 *dy = static_cast<float>((x8 + y1 - x5 - y2) - (y4 + x1 - y3 - x4)) / norm;
660}
661
662/* ---------------------------------------------------------------- */
663
664bool
666{
667 int const descr_scale = static_cast<int>(descr->scale);
668 int const width = this->sat->width();
669 int const height = this->sat->height();
670
671 /*
672 * Size of the descriptor is 20s (10s to each side). The Wavelet filter
673 * have size 2s (1s to each side), plus one additional pixel for simpler
674 * integral image lookup. Since the window can be rotated, 4s additional
675 * pixel are needed.
676 */
677 float const spacing = static_cast<float>(15 * descr_scale + 1);
678 if (descr->x < spacing || descr->y < spacing
679 || descr->x + spacing >= width || descr->y + spacing > height)
680 return false;
681
682 float sin_ori = 0.0f, cos_ori = 1.0f;
683 if (!upright)
684 {
685 sin_ori = std::sin(descr->orientation);
686 cos_ori = std::cos(descr->orientation);
687 }
688
689 /* Interest point region has size 20 * scale. */
690 descr->data.fill(0.0f);
691 float* descr_iter = &descr->data[0];
692 for (int y = -10; y < 10; ++y)
693 {
694 for (int x = -10; x < 10; ++x)
695 {
696 /* Rotate sample coordinate. */
697 int rot_x = static_cast<int>(math::round(descr->x
698 + (cos_ori * (x + 0.5f) - sin_ori * (y + 0.5f))
699 * descr_scale));
700 int rot_y = static_cast<int>(math::round(descr->y
701 + (sin_ori * (x + 0.5f) + cos_ori * (y + 0.5f))
702 * descr_scale));
703
704 /* Obtain and rotate gradient. */
705 float dx, dy;
706 this->filter_dx_dy(rot_x, rot_y, descr_scale, &dx, &dy);
707 float ori_dx = cos_ori * dx + sin_ori * dy;
708 float ori_dy = -sin_ori * dx + cos_ori * dy;
709
710 /* Keypoints are weighted with a Gaussian, sigma = 3.3*s. */
711 float weight = std::exp(-static_cast<float>(x * x + y * y)
712 / math::fastpow(2.0f * 3.3f, 2));
713 *(descr_iter++) += weight * ori_dx;
714 *(descr_iter++) += weight * ori_dy;
715 *(descr_iter++) += weight * std::abs(ori_dx);
716 *(descr_iter++) += weight * std::abs(ori_dy);
717
718 if ((x + 10) % 5 != 4)
719 descr_iter -= 4;
720 }
721 if ((y + 10) % 5 != 4)
722 descr_iter -= 4 * 4;
723 }
724
725 /* Normalize descriptor, reject too small descriptors. */
726 float norm = descr->data.square_norm();
727 if (MATH_EPSILON_EQ(norm, 0.0f, 1e-8)) // TODO: Play with this.
728 return false;
729 descr->data /= std::sqrt(norm);
730
731 return true;
732}
733
T square_norm(void) const
Computes the squared norm of the vector (much cheaper).
Definition vector.h:441
T minimum(void) const
Returns the smallest element in the vector.
Definition vector.h:397
T dot(Vector< T, N > const &other) const
Dot (or scalar) product between self and another vector.
Definition vector.h:542
Vector< T, N > & fill(T const &value)
Fills all vector elements with the given value.
Definition vector.h:381
T maximum(void) const
Returns the largest element in the vector.
Definition vector.h:404
std::shared_ptr< Image< T > > Ptr
Definition image.h:42
std::shared_ptr< Image< T > const > ConstPtr
Definition image.h:43
static Ptr create(void)
Smart pointer image constructor.
Definition image.h:191
std::vector< Keypoint > Keypoints
Definition surf.h:98
void create_octaves(void)
Definition surf.cc:143
void descriptor_assignment(void)
Definition surf.cc:479
bool keypoint_localization(Surf::Keypoint *kp)
Definition surf.cc:397
void keypoint_localization_and_filtering(void)
Definition surf.cc:379
SatType filter_dyy(int fs, int x, int y)
Definition surf.cc:246
void filter_dx_dy(int x, int y, int fs, float *dx, float *dy)
Definition surf.cc:624
void process(void)
Starts Surf keypoint detection and descriptor extraction.
Definition surf.cc:49
SatType filter_dxy(int fs, int x, int y)
Definition surf.cc:272
int64_t SatType
Signed type for the SAT image values.
Definition surf.h:128
void set_image(mve::ByteImage::ConstPtr image)
Sets the input image.
Definition surf.cc:124
bool descriptor_orientation(Descriptor *descr)
Definition surf.cc:516
bool descriptor_computation(Descriptor *descr, bool upright)
Definition surf.cc:665
void check_maximum(int o, int s, int x, int y)
Definition surf.cc:337
void extrema_detection(void)
Definition surf.cc:311
void create_response_map(int o, int k)
Definition surf.cc:159
SatType filter_dxx(int fs, int x, int y)
Definition surf.cc:222
Cross-platform high-resolution real-time timer.
Definition timer.h:30
std::size_t get_elapsed(void) const
Returns the milli seconds since last reset.
Definition timer.h:94
void reset(void)
Definition timer.h:88
#define MATH_EPSILON_EQ(x, v, eps)
Definition defines.h:96
#define MATH_PI
Definition defines.h:42
T round(T const &x)
Removes the fractional part of the value to the closest integer.
Definition functions.h:70
T fastpow(T const &base, unsigned int exp)
Takes base to the integer power of 'exp'.
Definition functions.h:235
Matrix< T, N, N > matrix_inverse(Matrix< T, N, N > const &mat)
Calculates the inverse of the given matrix.
T matrix_determinant(Matrix< T, N, N > const &mat)
Calculates the determinant of the given matrix.
@ DESATURATE_LIGHTNESS
(max(R,G,B) + min(R,G,B)) * 1/2
void swap(mve::Image< T > &a, mve::Image< T > &b)
Specialization of std::swap for efficient image swapping.
Definition image.h:478
#define SFM_NAMESPACE_END
Definition defines.h:14
#define SFM_NAMESPACE_BEGIN
Definition defines.h:13
Representation of a SURF descriptor.
Definition surf.h:84
float scale
The scale (or sigma value) of the keypoint.
Definition surf.h:90
float orientation
The orientation of the image keypoint in [-PI, PI].
Definition surf.h:92
float y
The sub-pixel y-coordinate of the image keypoint.
Definition surf.h:88
math::Vector< float, 64 > data
The descriptor data, elements are signed in [-1.0, 1.0].
Definition surf.h:94
float x
The sub-pixel x-coordinate of the image keypoint.
Definition surf.h:86
Representation of a SURF keypoint.
Definition surf.h:71
int octave
Octave index of the keypoint.
Definition surf.h:72
float x
Detected keypoint X coordinate.
Definition surf.h:74
float y
Detected keypoint Y coordinate.
Definition surf.h:75
float sample
Scale space sample index within octave in [0, 3].
Definition surf.h:73
float RespType
Type for the Hessian response value.
Definition surf.h:121
SURF options.
Definition surf.h:47
bool verbose_output
Produce status messages on the console.
Definition surf.h:59
float contrast_threshold
Sets the hessian threshold, defaults to 500.0.
Definition surf.h:51
bool debug_output
Produce even more messages on the console.
Definition surf.h:64
bool use_upright_descriptor
Trade rotation invariance for speed.
Definition surf.h:54
#define NUM_SAMPLES