1
2 /* autopano-sift, Automatic panorama image creation
3 * Copyright (C) 2004-2005 -- Sebastian Nowozin
4 *
5 * This program is free software released under the GNU General Public
6 * License, which is included in this software package (doc/LICENSE).
7 */
8
9 /* ScaleSpace.cs
10 *
11 * Key feature identification functionality, octave and scale-space handling.
12 *
13 * (C) Copyright 2004-2005 -- Sebastian Nowozin (nowozin@cs.tu-berlin.de)
14 *
15 * "The University of British Columbia has applied for a patent on the SIFT
16 * algorithm in the United States. Commercial applications of this software
17 * may require a license from the University of British Columbia."
18 * For more information, see the LICENSE file supplied with the distribution.
19 */
20
21 #include "AutoPanoSift.h"
22
OctavePyramid_new0()23 OctavePyramid* OctavePyramid_new0 ()
24 {
25 OctavePyramid* self = (OctavePyramid*)malloc(sizeof(OctavePyramid));
26 self->octaves = NULL;
27 self->Verbose = true;
28 return self;
29 }
30
OctavePyramid_delete(OctavePyramid * self)31 void OctavePyramid_delete(OctavePyramid* self)
32 {
33 if (self) {
34 ArrayList_delete(self->octaves);
35 memset(self, 0, sizeof(OctavePyramid));
36 free(self);
37 }
38 }
39
OctavePyramid_Count(OctavePyramid * self)40 int OctavePyramid_Count(OctavePyramid* self)
41 {
42 return ArrayList_Count(self->octaves);
43 }
44
OctavePyramid_GetScaleSpace(OctavePyramid * self,int idx)45 DScaleSpace* OctavePyramid_GetScaleSpace(OctavePyramid* self, int idx)
46 {
47 return (DScaleSpace*) ArrayList_GetItem(self->octaves, idx);
48 }
49
50 // Build the largest possible number of octaves, each holding
51 // 'levelsPerOctave' scales in scale-space. Each octave is downscaled by
52 // 0.5 and the scales in each octave represent a sigma change of
53 // 'octaveSigma' to 2.0 * 'octaveSigma'. If minSize is greater zero, every
54 // scale-space with less than minSize effective pixels in x-dimension will
55 // be discarded.
56 //
57 // Return the number of octaves build
OctavePyramid_BuildOctaves(OctavePyramid * self,ImageMap * source,double scale,int levelsPerOctave,double octaveSigma,int minSize)58 int OctavePyramid_BuildOctaves (OctavePyramid* self,
59 ImageMap* source, double scale,
60 int levelsPerOctave, double octaveSigma, int minSize)
61 {
62 self->octaves = ArrayList_new0 (DScaleSpace_delete);
63
64 DScaleSpace* downSpace = NULL;
65 ImageMap* prev = source;
66
67 while (prev != NULL && prev->xDim >= minSize && prev->yDim >= minSize) {
68 DScaleSpace* dsp = DScaleSpace_new0 ();
69 dsp->Verbose = self->Verbose;
70
71 if (self->Verbose)
72 WriteLine ("Building octave, (%d, %d)", prev->xDim, prev->yDim);
73
74 // Create both the gaussian filtered images and the DoG maps
75 DScaleSpace_BuildGaussianMaps (dsp, prev, scale, levelsPerOctave, octaveSigma);
76 DScaleSpace_BuildDiffMaps (dsp);
77
78 ArrayList_AddItem (self->octaves , dsp);
79
80 prev = ImageMap_ScaleHalf(DScaleSpace_LastGaussianMap(dsp));
81
82 if (downSpace != NULL)
83 downSpace->Up = dsp;
84
85 dsp->Down = downSpace;
86 downSpace = dsp;
87
88 scale *= 2.0;
89 }
90 if (prev != NULL) {
91 ImageMap_delete(prev);
92 }
93
94 return ArrayList_Count(self->octaves);
95 }
96
DScaleSpace_new0()97 DScaleSpace* DScaleSpace_new0()
98 {
99 DScaleSpace* self = (DScaleSpace*)malloc(sizeof(DScaleSpace));
100 self->Verbose = true;
101 self->Down = NULL;
102 self->Up = NULL;
103 self->baseImg = NULL;
104 self->magnitudes = NULL;
105 self->directions = NULL;
106 self->imgScaled = NULL;
107 self->spaces = NULL;
108 return self;
109 }
110
DScaleSpace_delete(DScaleSpace * self)111 void DScaleSpace_delete(DScaleSpace* self)
112 {
113 ArrayList_delete(self->magnitudes);
114 ArrayList_delete(self->directions);
115 ArrayList_delete(self->imgScaled);
116 ArrayList_delete(self->spaces);
117 free(self);
118 }
119
DScaleSpace_GetGaussianMap(DScaleSpace * self,int idx)120 ImageMap* DScaleSpace_GetGaussianMap (DScaleSpace* self, int idx)
121 {
122 return (ImageMap*) ArrayList_GetItem(self->imgScaled, idx);
123 }
124
125 // The last created gaussian map, with 2 \sigma blur. (For use with next
126 // octave).
DScaleSpace_LastGaussianMap(DScaleSpace * self)127 ImageMap* DScaleSpace_LastGaussianMap(DScaleSpace* self)
128 {
129 // position s = Length - 2 has: D(k^s * sigma) = D(2 * sigma)
130 if (ArrayList_Count(self->imgScaled) < 2)
131 FatalError ("bu keneng: too few gaussian maps");
132
133 return (ImageMap*) ArrayList_GetItem(self->imgScaled, ArrayList_Count(self->imgScaled) - 2);
134 }
135
DScaleSpace_Count(DScaleSpace * self)136 int DScaleSpace_Count(DScaleSpace* self)
137 {
138 return ArrayList_Count(self->spaces);
139 }
140
141 // Return a single DoG map
DScaleSpace_GetMap(DScaleSpace * self,int idx)142 ImageMap* DScaleSpace_GetMap(DScaleSpace* self, int idx)
143 {
144 return (ImageMap*) ArrayList_GetItem(self->spaces, idx);
145 }
146
147 // Generate keypoints from localized peak list.
148 // TODO: description
149 // scaleCount: number of scales (3)
DScaleSpace_GenerateKeypoints(DScaleSpace * self,ArrayList * localizedPeaks,int scaleCount,double octaveSigma)150 ArrayList* DScaleSpace_GenerateKeypoints(DScaleSpace* self,
151 ArrayList* localizedPeaks,
152 int scaleCount, double octaveSigma)
153 {
154 ArrayList* keypoints = ArrayList_new0 (NULL);
155
156 int i;
157 for (i=0; i<ArrayList_Count(localizedPeaks); i++) {
158 ScalePoint* sp = (ScalePoint*) ArrayList_GetItem(localizedPeaks, i);
159
160 // Generate zero or more keypoints from the scale point locations.
161 // TODO: make the values configurable
162 ArrayList* selfPointKeys = DScaleSpace_GenerateKeypointSingle (self, self->basePixScale,
163 sp, 36, 0.8, scaleCount, octaveSigma);
164
165 // Generate the feature descriptor.
166 selfPointKeys = DScaleSpace_CreateDescriptors (self, selfPointKeys,
167 (ImageMap*) ArrayList_GetItem(self->magnitudes, sp->level),
168 (ImageMap*) ArrayList_GetItem(self->directions, sp->level), 2.0, 4, 8, 0.2);
169
170 // Only copy over those keypoints that have been successfully
171 // assigned a descriptor (feature vector).
172 int j;
173 for(j=0; j<ArrayList_Count(selfPointKeys); j++) {
174 Keypoint* kp = (Keypoint*) ArrayList_GetItem(selfPointKeys, j);
175
176 if (kp->hasFV == false)
177 FatalError ("should not happen");
178
179 // Transform the this level image relative coordinate system
180 // to the original image coordinates by multiplying with the
181 // current img scale (which starts with either 0.5 or 1.0 and
182 // then always doubles: 2.0, 4.0, ..)
183 // Note that the kp coordinates are not used for processing by
184 // the detection methods and this has to be the last step.
185 // Also transform the relative-to-image scale to an
186 // absolute-to-original-image scale.
187 kp->x *= kp->imgScale;
188 kp->y *= kp->imgScale;
189 kp->scale *= kp->imgScale;
190
191 ArrayList_AddItem (keypoints, kp);
192 }
193 ArrayList_delete(selfPointKeys);
194 }
195
196 return (keypoints);
197 }
198
199 // Assign each feature point one or more standardized orientations.
200 // (section 5 in Lowe's paper)
201 //
202 // We use an orientation histogram with 36 bins, 10 degrees each. For
203 // this, every pixel (x,y) lieing in a circle of 'squareDim' diameter
204 // within in a 'squareDim' sized field within the image L ('gaussImg') is
205 // examined and two measures calculated:
206 //
207 // m = \sqrt{ (L_{x+1,y} - L_{x-1,y})^2 + (L_{x,y+1} - L_{x,y-1})^2 }
208 // theta = tan^{-1} ( \frac{ L_{x,y+1} - L_{x,y-1} }
209 // { L_{x+1,y} - L_{x-1,y} } )
210 //
211 // Where m is the gradient magnitude around the pixel and theta is the
212 // gradient orientation. The 'imgScale' value is the octave scale,
213 // starting with 1.0 at the finest-detail octave, and doubling every
214 // octave. The gradient orientations are discreetized to 'binCount'
215 // directions (should be: 36). For every peak orientation that lies within
216 // 'peakRelThresh' of the maximum peak value, a keypoint location is
217 // added (should be: 0.8).
218 //
219 // Note that 'space' is the gaussian smoothed original image, not the
220 // difference-of-gaussian one used for peak-search.
DScaleSpace_GenerateKeypointSingle(DScaleSpace * self,double imgScale,ScalePoint * point,int binCount,double peakRelThresh,int scaleCount,double octaveSigma)221 ArrayList* DScaleSpace_GenerateKeypointSingle (DScaleSpace* self,
222 double imgScale, ScalePoint* point,
223 int binCount, double peakRelThresh, int scaleCount,
224 double octaveSigma)
225 {
226 // The relative estimated keypoint scale. The actual absolute keypoint
227 // scale to the original image is yielded by the product of imgScale.
228 // But as we operate in the current octave, the size relative to the
229 // anchoring images is missing the imgScale factor.
230 double kpScale = octaveSigma *
231 pow (2.0, (point->level + point->local->scaleAdjust) / scaleCount);
232
233 // Lowe03, "A gaussian-weighted circular window with a \sigma three
234 // times that of the scale of the keypoint".
235 //
236 // With \sigma = 3.0 * kpScale, the square dimension we have to
237 // consider is (3 * \sigma) (until the weight becomes very small).
238 double sigma = 3.0 * kpScale;
239 int radius = (int) (3.0 * sigma / 2.0 + 0.5);
240 int radiusSq = radius * radius;
241
242 ImageMap* magnitude = (ImageMap*) ArrayList_GetItem(self->magnitudes, point->level);
243 ImageMap* direction = (ImageMap*) ArrayList_GetItem(self->directions, point->level);
244
245 // As the point may lie near the border, build the rectangle
246 // coordinates we can still reach, minus the border pixels, for which
247 // we do not have gradient information available.
248 int xMin = max (point->x - radius, 1);
249 int xMax = min (point->x + radius, magnitude->xDim - 1);
250 int yMin = max (point->y - radius, 1);
251 int yMax = min (point->y + radius, magnitude->yDim - 1);
252
253 // Precompute 1D gaussian divisor (2 \sigma^2) in:
254 // G(r) = e^{-\frac{r^2}{2 \sigma^2}}
255 double gaussianSigmaFactor = 2.0 * sigma * sigma;
256
257 double* bins = (double*)calloc(binCount, sizeof(double));
258
259 // Build the direction histogram
260 int y;
261 for (y = yMin ; y < yMax ; ++y) {
262 int x;
263 for ( x = xMin ; x < xMax ; ++x) {
264 // Only consider pixels in the circle, else we might skew the
265 // orientation histogram by considering more pixels into the
266 // corner directions
267 int relX = x - point->x;
268 int relY = y - point->y;
269 if (DScaleSpace_IsInCircle (relX, relY, radiusSq) == false)
270 continue;
271
272 // The gaussian weight factor.
273 double gaussianWeight = exp
274 (- ((relX * relX + relY * relY) / gaussianSigmaFactor));
275
276 // find the closest bin and add the direction
277 int binIdx = DScaleSpace_FindClosestRotationBin (self, binCount, ImageMap_GetPixel(direction, x, y));
278 bins[binIdx] += ImageMap_GetPixel(magnitude, x, y) * gaussianWeight;
279 }
280 }
281
282 // As there may be succeeding histogram entries like this:
283 // ( ..., 0.4, 0.3, 0.4, ... ) where the real peak is located at the
284 // middle of this three entries, we can improve the distinctiveness of
285 // the bins by applying an averaging pass.
286 //
287 // TODO: is this really the best method? (we also loose a bit of
288 // information. Maybe there is a one-step method that conserves more)
289 DScaleSpace_AverageWeakBins (self, bins, binCount);
290
291 // find the maximum peak in gradient orientation
292 double maxGrad = 0.0;
293 int maxBin = 0;
294 int b;
295 for (b = 0 ; b < binCount ; ++b) {
296 if (bins[b] > maxGrad) {
297 maxGrad = bins[b];
298 maxBin = b;
299 }
300 }
301
302 // First determine the real interpolated peak high at the maximum bin
303 // position, which is guaranteed to be an absolute peak.
304 //
305 // XXX: should we use the estimated peak value as reference for the
306 // 0.8 check or the original bin-value?
307 double maxPeakValue, maxDegreeCorrection;
308 DScaleSpace_InterpolateOrientation (self,
309 bins[maxBin == 0 ? (binCount - 1) : (maxBin - 1)],
310 bins[maxBin], bins[(maxBin + 1) % binCount],
311 &maxDegreeCorrection, &maxPeakValue);
312
313 // Now that we know the maximum peak value, we can find other keypoint
314 // orientations, which have to fulfill two criterias:
315 //
316 // 1. They must be a local peak themselves. Else we might add a very
317 // similar keypoint orientation twice (imagine for example the
318 // values: 0.4 1.0 0.8, if 1.0 is maximum peak, 0.8 is still added
319 // with the default threshhold, but the maximum peak orientation
320 // was already added).
321 // 2. They must have at least peakRelThresh times the maximum peak
322 // value.
323 bool* binIsKeypoint = (bool*)malloc(binCount*sizeof(bool));
324 for (b = 0 ; b < binCount ; ++b) {
325 binIsKeypoint[b] = false;
326
327 // The maximum peak of course is
328 if (b == maxBin) {
329 binIsKeypoint[b] = true;
330 continue;
331 }
332
333 // Local peaks are, too, in case they fulfill the threshhold
334 if (bins[b] < (peakRelThresh * maxPeakValue))
335 continue;
336
337 int leftI = (b == 0) ? (binCount - 1) : (b - 1);
338 int rightI = (b + 1) % binCount;
339 if (bins[b] <= bins[leftI] || bins[b] <= bins[rightI])
340 continue; // no local peak
341
342 binIsKeypoint[b] = true;
343 }
344
345 // All the valid keypoint bins are now marked in binIsKeypoint, now
346 // build them.
347 ArrayList* keypoints = ArrayList_new0 (NULL);
348
349 // find other possible locations
350 double oneBinRad = (2.0 * M_PI) / binCount;
351
352 for (b = 0 ; b < binCount ; ++b) {
353 if (binIsKeypoint[b] == false)
354 continue;
355
356 int bLeft = (b == 0) ? (binCount - 1) : (b - 1);
357 int bRight = (b + 1) % binCount;
358
359 // Get an interpolated peak direction and value guess.
360 double peakValue;
361 double degreeCorrection;
362
363 if (DScaleSpace_InterpolateOrientation (self, bins[bLeft], bins[b], bins[bRight],
364 °reeCorrection, &peakValue) == false)
365 {
366 FatalError("BUG: Parabola fitting broken");
367 }
368
369 // [-1.0 ; 1.0] -> [0 ; binrange], and add the fixed absolute bin
370 // position.
371 // We subtract PI because bin 0 refers to 0, binCount-1 bin refers
372 // to a bin just below 2PI, so -> [-PI ; PI]. Note that at this
373 // point we determine the canonical descriptor anchor angle. It
374 // does not matter where we set it relative to the peak degree,
375 // but it has to be constant. Also, if the output of this
376 // implementation is to be matched with other implementations it
377 // must be the same constant angle (here: -PI).
378 double degree = (b + degreeCorrection) * oneBinRad - M_PI;
379
380 if (degree < -M_PI)
381 degree += 2.0 * M_PI;
382 else if (degree > M_PI)
383 degree -= 2.0 * M_PI;
384
385 Keypoint* kp = Keypoint_new ((ImageMap*)ArrayList_GetItem(self->imgScaled, point->level),
386 point->x + point->local->fineX,
387 point->y + point->local->fineY,
388 imgScale, kpScale, degree);
389 ArrayList_AddItem (keypoints, kp);
390 }
391
392 free(binIsKeypoint);
393 free(bins);
394 return (keypoints);
395 }
396
397 // Fit a parabol to the three points (-1.0 ; left), (0.0 ; middle) and
398 // (1.0 ; right).
399 //
400 // Formulas:
401 // f(x) = a (x - c)^2 + b
402 //
403 // c is the peak offset (where f'(x) is zero), b is the peak value.
404 //
405 // In case there is an error false is returned, otherwise a correction
406 // value between [-1 ; 1] is returned in 'degreeCorrection', where -1
407 // means the peak is located completely at the left vector, and -0.5 just
408 // in the middle between left and middle and > 0 to the right side. In
409 // 'peakValue' the maximum estimated peak value is stored.
DScaleSpace_InterpolateOrientation(DScaleSpace * self,double left,double middle,double right,double * degreeCorrection,double * peakValue)410 bool DScaleSpace_InterpolateOrientation (DScaleSpace* self,
411 double left, double middle,
412 double right, double* degreeCorrection, double* peakValue)
413 {
414 double a = ((left + right) - 2.0 * middle) / 2.0;
415 *degreeCorrection = *peakValue = -1;
416
417 // Not a parabol
418 if (a == 0.0) {
419 *degreeCorrection = 0;
420 return (true);
421 }
422
423 double c = (((left - middle) / a) - 1.0) / 2.0;
424 double b = middle - c * c * a;
425
426 if (c < -0.5 || c > 0.5)
427 FatalError
428 ("InterpolateOrientation: off peak ]-0.5 ; 0.5[");
429
430 *degreeCorrection = c;
431 *peakValue = b;
432
433 return (true);
434 }
435
436 // Find the bin out of 'binCount' bins that matches the 'angle' closest.
437 // 'angle' fulfills -PI <= angle <= PI. Bin 0 is assigned to -PI, the
438 // binCount-1 bin refers to just below PI.
439 //
440 // Return the index of the closest bin.
DScaleSpace_FindClosestRotationBin(DScaleSpace * self,int binCount,double angle)441 int DScaleSpace_FindClosestRotationBin (DScaleSpace* self, int binCount, double angle)
442 {
443 angle += M_PI;
444 angle /= 2.0 * M_PI;
445
446 // calculate the aligned bin
447 angle *= binCount;
448
449 int idx = (int) angle;
450 if (idx == binCount)
451 idx = 0;
452
453 return (idx);
454 }
455
456 // Average the content of the direction bins.
DScaleSpace_AverageWeakBins(DScaleSpace * self,double * bins,int binCount)457 void DScaleSpace_AverageWeakBins (DScaleSpace* self, double* bins, int binCount)
458 {
459 // TODO: make some tests what number of passes is the best. (its clear
460 // one is not enough, as we may have something like
461 // ( 0.4, 0.4, 0.3, 0.4, 0.4 ))
462 int sn;
463 for ( sn = 0 ; sn < 4 ; ++sn) {
464 double firstE = bins[0];
465 double last = bins[binCount - 1];
466
467 int sw;
468 for ( sw = 0 ; sw < binCount ; ++sw) {
469 double cur = bins[sw];
470 double next = (sw == (binCount - 1)) ?
471 firstE : bins[(sw + 1) % binCount];
472
473 bins[sw] = (last + cur + next) / 3.0;
474 last = cur;
475 }
476 }
477 }
478
479 // Create the descriptor vector for a list of keypoints.
480 //
481 // keypoints: The list of keypoints to be processed. Everything but the
482 // descriptor must be filled in already.
483 // magnitude/direction: The precomputed gradient magnitude and direction
484 // maps.
485 // considerScaleFactor: The downscale factor, which describes the amount
486 // of pixels in the circular region relative to the keypoint scale.
487 // Low values means few pixels considered, large values extend the
488 // range. (Use values between 1.0 and 6.0)
489 // descDim: The dimension size of the output descriptor. There will be
490 // descDim * descDim * directionCount elements in the feature vector.
491 // directionCount: The dimensionality of the low level gradient vectors.
492 // fvGradHicap: The feature vector gradient length hi-cap threshhold.
493 // (Should be: 0.2)
494 //
495 // Some parts modelled after Alexandre Jenny's Matlab implementation.
496 //
497 // Return a list of survivors, which a descriptor was created for
498 // successfully.
DScaleSpace_CreateDescriptors(DScaleSpace * self,ArrayList * keypoints,ImageMap * magnitude,ImageMap * direction,double considerScaleFactor,int descDim,int directionCount,double fvGradHicap)499 ArrayList* DScaleSpace_CreateDescriptors (DScaleSpace* self,
500 ArrayList* keypoints,
501 ImageMap* magnitude, ImageMap* direction,
502 double considerScaleFactor, int descDim, int directionCount,
503 double fvGradHicap)
504 {
505 if (ArrayList_Count(keypoints) <= 0)
506 return (keypoints);
507
508 considerScaleFactor *= ((Keypoint*) ArrayList_GetItem(keypoints, 0))->scale;
509 double dDim05 = ((double) descDim) / 2.0;
510
511 // Now calculate the radius: We consider pixels in a square with
512 // dimension 'descDim' plus 0.5 in each direction. As the feature
513 // vector elements at the diagonal borders are most distant from the
514 // center pixel we have scale up with sqrt(2).
515 int radius = (int) (((descDim + 1.0) / 2) *
516 sqrt (2.0) * considerScaleFactor + 0.5);
517
518 // Instead of modifying the original list, we just copy the keypoints
519 // that received a descriptor.
520 ArrayList* survivors = ArrayList_new0 (NULL);
521
522 // Precompute the sigma for the "center-most, border-less" gaussian
523 // weighting.
524 // (We are operating to dDim05, CV book tells us G(x), x > 3 \sigma
525 // negligible, but this range seems much shorter!?)
526 //
527 // In Lowe03, page 15 it says "A Gaussian weighting function with
528 // \sigma equal to one half the width of the descriptor window is
529 // used", so we just use his advice.
530 double sigma2Sq = 2.0 * dDim05 * dDim05;
531
532 int i;
533 for (i=0; i<ArrayList_Count(keypoints); i++) {
534 Keypoint* kp = (Keypoint*) ArrayList_GetItem(keypoints,i);
535 // The angle to rotate with: negate the orientation.
536 double angle = -kp->orientation;
537
538 Keypoint_CreateVector (kp, descDim, descDim, directionCount);
539 //Console.WriteLine (" FV allocated");
540
541 int y;
542 for (y = -radius ; y < radius ; ++y) {
543 int x;
544 for (x = -radius ; x < radius ; ++x) {
545 // Rotate and scale
546 double yR = sin (angle) * x +
547 cos (angle) * y;
548 double xR = cos (angle) * x -
549 sin (angle) * y;
550
551 yR /= considerScaleFactor;
552 xR /= considerScaleFactor;
553
554 // Now consider all (xR, yR) that are anchored within
555 // (- descDim/2 - 0.5 ; -descDim/2 - 0.5) to
556 // (descDim/2 + 0.5 ; descDim/2 + 0.5),
557 // as only those can influence the FV.
558 if (yR >= (dDim05 + 0.5) || xR >= (dDim05 + 0.5) ||
559 xR <= -(dDim05 + 0.5) || yR <= -(dDim05 + 0.5))
560 continue;
561
562 int currentX = (int) (x + kp->x + 0.5);
563 int currentY = (int) (y + kp->y + 0.5);
564 if (currentX < 1 || currentX >= (magnitude->xDim - 1) ||
565 currentY < 1 || currentY >= (magnitude->yDim - 1))
566 continue;
567
568 /*
569 WriteLine (" (%d,%d) by angle {%lf} -> (%lf,%lf)",
570 x, y, angle, xR, yR);
571 */
572
573 // Weight the magnitude relative to the center of the
574 // whole FV. We do not need a normalizing factor now, as
575 // we normalize the whole FV later anyway (see below).
576 // xR, yR are each in -(dDim05 + 0.5) to (dDim05 + 0.5)
577 // range
578 double magW = exp (-(xR * xR + yR * yR) / sigma2Sq) *
579 ImageMap_GetPixel(magnitude, currentX, currentY);
580
581 // Anchor to (-1.0, -1.0)-(dDim + 1.0, dDim + 1.0), where
582 // the FV points are located at (x, y)
583 yR += dDim05 - 0.5;
584 xR += dDim05 - 0.5;
585
586 // Build linear interpolation weights:
587 // A B
588 // C D
589 //
590 // The keypoint is located between A, B, C and D.
591 int xIdx[2] = {0,0};
592 int yIdx[2] = {0,0};
593 int dirIdx[2] ={0,0};
594 double xWeight[2] = {0.0,0.0};
595 double yWeight[2] = {0.0,0.0};
596 double dirWeight[2] = {0.0,0.0};
597
598 if (xR >= 0) {
599 xIdx[0] = (int) xR;
600 xWeight[0] = (1.0 - (xR - xIdx[0]));
601 }
602 if (yR >= 0) {
603 yIdx[0] = (int) yR;
604 yWeight[0] = (1.0 - (yR - yIdx[0]));
605 }
606
607 if (xR < (descDim - 1)) {
608 xIdx[1] = (int) (xR + 1.0);
609 xWeight[1] = xR - xIdx[1] + 1.0;
610 }
611 if (yR < (descDim - 1)) {
612 yIdx[1] = (int) (yR + 1.0);
613 yWeight[1] = yR - yIdx[1] + 1.0;
614 }
615
616 // Rotate the gradient direction by the keypoint
617 // orientation, then normalize to [-pi ; pi] range.
618 double dir = ImageMap_GetPixel(direction, currentX, currentY) - kp->orientation;
619 if (dir <= -M_PI)
620 dir += M_PI;
621 if (dir > M_PI)
622 dir -= M_PI;
623
624 double idxDir = (dir * directionCount) /
625 (2.0 * M_PI);
626 if (idxDir < 0.0)
627 idxDir += directionCount;
628
629 dirIdx[0] = (int) idxDir;
630 dirIdx[1] = (dirIdx[0] + 1) % directionCount;
631 dirWeight[0] = 1.0 - (idxDir - dirIdx[0]);
632 dirWeight[1] = idxDir - dirIdx[0];
633
634 /*
635 WriteLine (" (%lf,%lf) yields:", xR, yR);
636 WriteLine (" x<%d,%d>*(%lf,%lf)",
637 xIdx[0], xIdx[1], xWeight[0], xWeight[1]);
638 WriteLine (" y<%d,%d>*(%lf,%lf)",
639 yIdx[0], yIdx[1], yWeight[0], yWeight[1]);
640 WriteLine (" dir<%d,%d>*(%lf,%lf)",
641 dirIdx[0], dirIdx[1], dirWeight[0], dirWeight[1]);
642 WriteLine (" weighting m * w: %lf * %lf",
643 ImageMap_GetPixel(magnitude,currentX, currentY), exp (-(xR * xR +
644 yR * yR) / sigma2Sq));
645 */
646 int iy;
647 for (iy = 0 ; iy < 2 ; ++iy) {
648 int ix;
649 for (ix = 0 ; ix < 2 ; ++ix) {
650 int id;
651 for (id = 0 ; id < 2 ; ++id) {
652 Keypoint_FVSet (kp, xIdx[ix], yIdx[iy], dirIdx[id],
653 Keypoint_FVGet (kp, xIdx[ix], yIdx[iy], dirIdx[id]) +
654 xWeight[ix] * yWeight[iy] * dirWeight[id] * magW);
655 }
656 }
657 }
658 }
659 }
660
661 // Normalize and hicap the feature vector, as recommended on page
662 // 16 in Lowe03.
663 DScaleSpace_CapAndNormalizeFV (self, kp, fvGradHicap);
664
665 ArrayList_AddItem (survivors, kp);
666 }
667
668 ArrayList_delete(keypoints);
669 return (survivors);
670 }
671
672
673 // Threshhold and normalize feature vector.
674 // Note that the feature vector as a whole is normalized (Lowe's paper is
675 // a bit unclear at that point).
DScaleSpace_CapAndNormalizeFV(DScaleSpace * self,Keypoint * kp,double fvGradHicap)676 void DScaleSpace_CapAndNormalizeFV (DScaleSpace* self, Keypoint* kp, double fvGradHicap)
677 {
678 // Straight normalization
679 double norm = 0.0;
680 int n;
681 for (n = 0 ; n < Keypoint_FVLinearDim(kp) ; ++n)
682 norm += pow (Keypoint_FVLinearGet (kp, n), 2.0);
683
684 norm = sqrt (norm);
685 if (norm == 0.0)
686 return;
687
688 for (n = 0 ; n < Keypoint_FVLinearDim(kp) ; ++n)
689 Keypoint_FVLinearSet (kp, n, Keypoint_FVLinearGet (kp, n) / norm);
690
691 // Hicap after normalization
692 for (n = 0 ; n < Keypoint_FVLinearDim(kp) ; ++n) {
693 if (Keypoint_FVLinearGet (kp, n) > fvGradHicap) {
694 Keypoint_FVLinearSet (kp, n, fvGradHicap);
695 }
696 }
697
698 // Renormalize again
699 norm = 0.0;
700 for ( n = 0 ; n < Keypoint_FVLinearDim(kp) ; ++n)
701 norm += pow (Keypoint_FVLinearGet (kp, n), 2.0);
702
703 norm = sqrt (norm);
704
705 for ( n = 0 ; n < Keypoint_FVLinearDim(kp) ; ++n)
706 Keypoint_FVLinearSet (kp, n, Keypoint_FVLinearGet (kp, n) / norm);
707 }
708
709 // Simple helper predicate to tell if (rX, rY) is within a circle of
710 // \sqrt{radiusSq} radius, assuming the circle center is (0, 0).
DScaleSpace_IsInCircle(int rX,int rY,int radiusSq)711 bool DScaleSpace_IsInCircle (int rX, int rY, int radiusSq)
712 {
713 rX *= rX;
714 rY *= rY;
715 if ((rX + rY) <= radiusSq)
716 return (true);
717
718 return (false);
719 }
720
721 // Remove peaks by peak magnitude and peak edge response. Find the
722 // sub-pixel local offset by interpolation.
723 //
724 // Sub-pixel localization and peak magnitude:
725 // After this method returns, every peak has a relative localization
726 // offset and its peak magnitude available in 'peak.Local'. The peak
727 // magnitude value must be above 'dValueLoThresh' for the point to
728 // survive. Usual values might lie in the range 0.0 (no filtering) to
729 // 0.03 (Lowe/Brown's recommendation). We normally use a value around
730 // 0.0001 to 0.00025 (and Brown's values seem quite large to me). The
731 // scaleAdjustThresh value is explained in LoweDetector.cs.
732 //
733 // Edge filtering:
734 // 'edgeRatio' denotes the required hi-threshhold for the ratio between
735 // the principle curvatures. Small values (1.5 to 3.0) will filter most
736 // points, leaving only the most corner-like points. Larger values (3.0 to
737 // 10.0) will remove the points which lie on a straight edge, whose
738 // position might be more vulnerable to noise.
739 //
740 // Return a filtered list of ScalePoint elements, with only the remaining
741 // survivors.
DScaleSpace_FilterAndLocalizePeaks(DScaleSpace * self,ArrayList * peaks,double edgeRatio,double dValueLoThresh,double scaleAdjustThresh,int relocationMaximum)742 ArrayList* DScaleSpace_FilterAndLocalizePeaks (DScaleSpace* self, ArrayList* peaks, double edgeRatio,
743 double dValueLoThresh, double scaleAdjustThresh, int relocationMaximum)
744 {
745 ArrayList* filtered = ArrayList_new0 (NULL);
746
747 ImageMap* space0 = (ImageMap*) ArrayList_GetItem(self->spaces, 0);
748 int** processed = IntMap_new(space0->xDim, space0->yDim);
749
750 int i;
751 for(i=0; i<ArrayList_Count(peaks); i++) {
752 ScalePoint* peak = (ScalePoint*)ArrayList_GetItem(peaks, i);
753
754 if (DScaleSpace_IsTooEdgelike (self, (ImageMap*) ArrayList_GetItem(self->spaces, peak->level), peak->x, peak->y, edgeRatio))
755 continue;
756
757 // When the localization hits some problem, i.e. while moving the
758 // point a border is reached, then skip this point.
759 if (DScaleSpace_LocalizeIsWeak (self, peak, relocationMaximum, processed))
760 continue;
761
762 // Now we approximated the exact sub-pixel peak position.
763 // Comment the following line out to get a number of image files
764 // which show the located peak in the closest DoG scale.
765 /*DEBUGSaveRectangle (spaces[peak.Level], peak.X, peak.Y,
766 String.Format ("rect_{0}.png", peak.Local.DValue);
767 */
768
769 /*WriteLine ("peak.Local.ScaleAdjust = %f",
770 peak->local->scaleAdjust);*/
771 if (abs (peak->local->scaleAdjust) > scaleAdjustThresh)
772 continue;
773
774 // Additional local pixel information is now available, threshhold
775 // the D(^x)
776 /*WriteLine ("%d %d %f # == DVALUE", peak->y, peak->x, peak->local->dValue);*/
777 if (abs (peak->local->dValue) <= dValueLoThresh)
778 continue;
779
780 /*WriteLine ("%d %d %f %f # FILTERLOCALIZE",
781 peak->y, peak->x, peak->local->scaleAdjust, peak->local->dValue);*/
782
783 // its edgy enough, add it
784 ArrayList_AddItem (filtered, peak);
785 }
786
787 IntMap_delete(processed);
788 return (filtered);
789 }
790
791 // Return true if the point is not suitable, either because it lies on a
792 // border pixel or the Hessian matrix cannot be inverted.
793 // If false is returned, the pixel is suitable for localization and
794 // additional localization information has been put into 'point.Local'.
795 // No more than 'steps' corrections are made.
DScaleSpace_LocalizeIsWeak(DScaleSpace * self,ScalePoint * point,int steps,int ** processed)796 bool DScaleSpace_LocalizeIsWeak (DScaleSpace* self, ScalePoint* point, int steps, int** processed)
797 {
798 bool needToAdjust = true;
799 int adjusted = steps;
800
801 while (needToAdjust) {
802 int x = point->x;
803 int y = point->y;
804
805 // Points we cannot say anything about, as they lie on the border
806 // of the scale space
807 if (point->level <= 0 || point->level >= (ArrayList_Count(self->spaces) - 1))
808 return (true);
809
810 ImageMap* space = (ImageMap*) ArrayList_GetItem(self->spaces,point->level);
811 if (x <= 0 || x >= (space->xDim - 1))
812 return (true);
813 if (y <= 0 || y >= (space->yDim - 1))
814 return (true);
815
816 double dp;
817 SimpleMatrix* adj = DScaleSpace_GetAdjustment (self, point, point->level, x, y, &dp);
818
819 // Get adjustments and check if we require further adjustments due
820 // to pixel level moves. If so, turn the adjustments into real
821 // changes and continue the loop. Do not adjust the plane, as we
822 // are usually quite low on planes in thie space and could not do
823 // further adjustments from the top/bottom planes.
824 double adjS = SimpleMatrix_GetValue(adj, 0, 0);
825 double adjY = SimpleMatrix_GetValue(adj, 1, 0);
826 double adjX = SimpleMatrix_GetValue(adj, 2, 0);
827 SimpleMatrix_delete(adj);
828 if (abs (adjX) > 0.5 || abs (adjY) > 0.5) {
829 // Already adjusted the last time, give up
830 if (adjusted == 0) {
831 //WriteLine ("too many adjustments, returning");
832 return (true);
833 }
834
835 adjusted -= 1;
836
837 // Check that just one pixel step is needed, otherwise discard
838 // the point
839 double distSq = adjX * adjX + adjY * adjY;
840 if (distSq > 2.0)
841 return (true);
842
843 point->x = (int) (point->x + adjX + 0.5);
844 point->y = (int) (point->y + adjY + 0.5);
845 //point->level = (int) (point->level + adjS + 0.5);
846
847 /*WriteLine ("moved point by ({0},{1}: {2}) to ({3},{4}: {5})",
848 adjX, adjY, adjS, point.X, point.Y, point.Level);*/
849 continue;
850 }
851
852 /* for processing with gnuplot
853 *
854 Console.WriteLine ("{0} {1} # POINT LEVEL {2}", point.X,
855 point.Y, basePixScale);
856 Console.WriteLine ("{0} {1} {2} # ADJ POINT LEVEL {3}",
857 adjS, adjX, adjY, basePixScale);
858 */
859
860 // Check if we already have a keypoint within this octave for this
861 // pixel position in order to avoid dupes. (Maybe we can move this
862 // check earlier after any adjustment, so we catch dupes earlier).
863 // If its not in there, mark it for later searches.
864 //
865 // FIXME: check why there does not seem to be a dupe at all
866 if (processed[point->x][point->y] != 0)
867 return (true);
868
869 processed[point->x][point->y] = 1;
870
871 // Save final sub-pixel adjustments.
872 PointLocalInformation* local = PointLocalInformation_new3 (adjS, adjX, adjY);
873 //local.DValue = dp;
874 local->dValue = ImageMap_GetPixel(space, point->x, point->y) + 0.5 * dp;
875 point->local = local;
876
877 needToAdjust = false;
878 }
879
880 return (false);
881 }
882
DScaleSpace_IsTooEdgelike(DScaleSpace * self,ImageMap * space,int x,int y,double r)883 bool DScaleSpace_IsTooEdgelike (DScaleSpace* self, ImageMap* space, int x, int y, double r)
884 {
885 double D_xx, D_yy, D_xy;
886
887 // Calculate the Hessian H elements [ D_xx, D_xy ; D_xy , D_yy ]
888 D_xx = space->values[x + 1][y] + space->values[x - 1][y] - 2.0 * space->values[x][y];
889 D_yy = space->values[x][y + 1] + space->values[x][y - 1] - 2.0 * space->values[x][y];
890 D_xy = 0.25 * ((space->values[x + 1][y + 1] - space->values[x + 1][y - 1]) -
891 (space->values[x - 1][y + 1] - space->values[x - 1][y - 1]));
892
893 // page 13 in Lowe's paper
894 double TrHsq = D_xx + D_yy;
895 TrHsq *= TrHsq;
896 double DetH = D_xx * D_yy - (D_xy * D_xy);
897
898 double r1sq = (r + 1.0);
899 r1sq *= r1sq;
900
901 // BUG: this can invert < to >, uhh: if ((TrHsq * r) < (DetH * r1sq))
902 if ((TrHsq / DetH) < (r1sq / r)) {
903 /*Console.WriteLine ("{0} {1} {2} {3} {4} # EDGETEST",
904 y, x, (TrHsq * r), (DetH * r1sq),
905 (TrHsq / DetH) / (r1sq / r));*/
906 return (false);
907 }
908
909 return (true);
910 }
911
912 // Return adjustment (scale, y, x) on success,
913 // return null on failure
914 // TODO: integrate this
DScaleSpace_GetAdjustment(DScaleSpace * self,ScalePoint * point,int level,int x,int y,double * dp)915 SimpleMatrix* DScaleSpace_GetAdjustment (DScaleSpace* self, ScalePoint* point,
916 int level, int x, int y, double* dp)
917 {
918 /*WriteLine ("GetAdjustment (point, %d, %d, %d, out double dp)",
919 level, x, y);*/
920 *dp = 0.0;
921 if (point->level <= 0 || point->level >= (ArrayList_Count(self->spaces) - 1))
922 FatalError ("point.Level is not within [bottom-1;top-1] range");
923
924 ImageMap* below = (ImageMap*) ArrayList_GetItem(self->spaces, level - 1);
925 ImageMap* current = (ImageMap*) ArrayList_GetItem(self->spaces, level);
926 ImageMap* above = (ImageMap*) ArrayList_GetItem(self->spaces, level + 1);
927
928 SimpleMatrix* H = SimpleMatrix_new (3, 3);
929 H->values[0][0] = below->values[x][y] - 2 * current->values[x][y] + above->values[x][y];
930 H->values[0][1] = H->values[1][0] = 0.25 * (above->values[x][y + 1] - above->values[x][y - 1] -
931 (below->values[x][y + 1] - below->values[x][y - 1]));
932 H->values[0][2] = H->values[2][0] = 0.25 * (above->values[x + 1][y] - above->values[x - 1][y] -
933 (below->values[x + 1][y] - below->values[x - 1][y]));
934 H->values[1][1] = current->values[x][y - 1] - 2 * current->values[x][y] + current->values[x][y + 1];
935 H->values[1][2] = H->values[2][1] = 0.25 * (current->values[x + 1][y + 1] - current->values[x - 1][y + 1] -
936 (current->values[x + 1][y - 1] - current->values[x - 1][y - 1]));
937 H->values[2][2] = current->values[x - 1][y] - 2 * current->values[x][y] + current->values[x + 1][y];
938
939 SimpleMatrix* d = SimpleMatrix_new (3, 1);
940 d->values[0][0] = 0.5 * (above->values[x][y] - below->values[x][y]);
941 d->values[1][0] = 0.5 * (current->values[x][y + 1] - current->values[x][y - 1]);
942 d->values[2][0] = 0.5 * (current->values[x + 1][y] - current->values[x - 1][y]);
943
944 SimpleMatrix* b = SimpleMatrix_clone (d);
945 SimpleMatrix_Negate (b);
946
947 // Solve: A x = b
948 SimpleMatrix_SolveLinear (H, b);
949
950 SimpleMatrix_delete(H);
951 *dp = SimpleMatrix_Dot (b, d);
952
953 SimpleMatrix_delete(d);
954
955 return (b);
956 }
957
958
959 // Peak localization in scale-space.
960 //
961 // From Lowe's paper its not really clear whether we always need three
962 // neighbourhood spaces or should also search only one or two spaces. As
963 // figure 3 might suggest the later, we do it like this.
964 //
965 // Return an arraylist holding ScalePoint elements.
DScaleSpace_FindPeaks(DScaleSpace * self,double dogThresh)966 ArrayList* DScaleSpace_FindPeaks (DScaleSpace* self, double dogThresh)
967 {
968 if (self->Verbose)
969 WriteLine (" FindPeaks: scale %02.2f, testing %d levels",
970 self->basePixScale, DScaleSpace_Count(self) - 2);
971
972 ArrayList* peaks = ArrayList_new0 (ScalePoint_delete);
973
974 ImageMap *current, *above, *below;
975
976 // Search the D(k * sigma) to D(2 * sigma) spaces
977 int level;
978 for (level = 1 ; level < (DScaleSpace_Count(self) - 1) ; ++level)
979 {
980 current = DScaleSpace_GetMap(self, level);
981 below = DScaleSpace_GetMap(self, level - 1);
982 above = DScaleSpace_GetMap(self, level + 1);
983
984 //// WriteLine ("peak-search at level %d", level);
985 /*
986 WriteLine ("below/current/above: %s %s %s",
987 below == NULL ? "-" : "X",
988 current == NULL ? "-" : "X",
989 above == NULL ? "-" : "X");
990 WriteLine ("peak-search at level %d", level);
991 */
992 ArrayList* newPeaks;
993 ArrayList_AddRange (peaks, newPeaks = DScaleSpace_FindPeaksThreeLevel (self, below, current, above,
994 level, dogThresh));
995 ArrayList_delete(newPeaks);
996
997 below = current;
998 }
999
1000 return (peaks);
1001 }
1002
DScaleSpace_FindPeaksThreeLevel(DScaleSpace * self,ImageMap * below,ImageMap * current,ImageMap * above,int curLev,double dogThresh)1003 ArrayList* DScaleSpace_FindPeaksThreeLevel (DScaleSpace* self, ImageMap* below, ImageMap* current,
1004 ImageMap* above, int curLev, double dogThresh)
1005 {
1006 ArrayList* peaks = ArrayList_new0 (NULL);
1007
1008 int x, y;
1009 for ( y = 1 ; y < (current->yDim - 1) ; ++y) {
1010 for ( x = 1 ; x < (current->xDim - 1) ; ++x) {
1011 bool cIsMax = true;
1012 bool cIsMin = true;
1013
1014 double c = current->values[x][y]; // Center value
1015
1016 /* If the magnitude is below the threshhold, skip it early.
1017 */
1018 if (abs (c) <= dogThresh)
1019 continue;
1020
1021 DScaleSpace_CheckMinMax (self, current, c, x, y, &cIsMin, &cIsMax, true);
1022 DScaleSpace_CheckMinMax (self, below, c, x, y, &cIsMin, &cIsMax, false);
1023 DScaleSpace_CheckMinMax (self, above, c, x, y, &cIsMin, &cIsMax, false);
1024 if (cIsMin == false && cIsMax == false)
1025 continue;
1026
1027 //WriteLine ("%d %d %f # DOG", y, x, c);
1028
1029 /* Add the peak that survived the first checks, to the peak
1030 * list.
1031 */
1032 ArrayList_AddItem (peaks, ScalePoint_new3 (x, y, curLev));
1033 }
1034 }
1035
1036 return (peaks);
1037 }
1038
1039 // Check if a pixel ('x', 'y') with value 'c' is minimum or maximum in the
1040 // 'layer' image map. Except for the center, and its above/below planes
1041 // corresponding pixels, use a strong > and < check (because the if is
1042 // inverted it looks like >= and <=).
DScaleSpace_CheckMinMax(DScaleSpace * self,ImageMap * layer,double c,int x,int y,bool * IsMin,bool * IsMax,bool cLayer)1043 void DScaleSpace_CheckMinMax (DScaleSpace* self, ImageMap* layer, double c, int x, int y,
1044 bool* IsMin, bool* IsMax, bool cLayer)
1045 {
1046 if (layer == NULL)
1047 return;
1048
1049 if (*IsMin == true) {
1050 if (layer->values[x - 1][y - 1] <= c ||
1051 layer->values[x][y - 1] <= c ||
1052 layer->values[x + 1][y - 1] <= c ||
1053 layer->values[x - 1][y] <= c ||
1054 // note here its just < instead of <=
1055 (cLayer ? false : (layer->values[x][y] < c)) ||
1056 layer->values[x + 1][y] <= c ||
1057 layer->values[x - 1][y + 1] <= c ||
1058 layer->values[x][y + 1] <= c ||
1059 layer->values[x + 1][y + 1] <= c)
1060 *IsMin = false;
1061 }
1062 if (*IsMax == true) {
1063 if (layer->values[x - 1][y - 1] >= c ||
1064 layer->values[x][y - 1] >= c ||
1065 layer->values[x + 1][y - 1] >= c ||
1066 layer->values[x - 1][y] >= c ||
1067 // note here its just > instead of >=
1068 (cLayer ? false : (layer->values[x][y] > c)) ||
1069 layer->values[x + 1][y] >= c ||
1070 layer->values[x - 1][y + 1] >= c ||
1071 layer->values[x][y + 1] >= c ||
1072 layer->values[x + 1][y + 1] >= c)
1073 *IsMax = false;
1074 }
1075 }
1076
DScaleSpace_SToK(int s)1077 double DScaleSpace_SToK (int s)
1078 {
1079 return (pow (2.0, 1.0 / s));
1080 }
1081
1082 // Precompute all gradient magnitude and direction planes for one octave.
DScaleSpace_GenerateMagnitudeAndDirectionMaps(DScaleSpace * self)1083 void DScaleSpace_GenerateMagnitudeAndDirectionMaps (DScaleSpace* self)
1084 {
1085 // We leave the first entry to null, and ommit the last. This way, the
1086 // magnitudes and directions maps have the same index as the
1087 // imgScaled maps they below to.
1088 int Count = DScaleSpace_Count(self);
1089 self->magnitudes = ArrayList_new(Count - 1, ImageMap_delete);
1090 self->directions = ArrayList_new(Count - 1, ImageMap_delete);
1091
1092 // Build the maps, omitting the border pixels, as we cannot build
1093 // gradient information there.
1094 int s;
1095 for ( s = 1 ; s < (Count - 1) ; ++s) {
1096 ImageMap* imgScaledAtS = (ImageMap*) ArrayList_GetItem(self->imgScaled,s);
1097
1098 ImageMap* magnitudeAtS = ImageMap_new (imgScaledAtS->xDim, imgScaledAtS->yDim);
1099 ImageMap* directionsAtS = ImageMap_new (imgScaledAtS->xDim, imgScaledAtS->yDim);
1100
1101 ArrayList_SetItem(self->magnitudes, s, magnitudeAtS);
1102 ArrayList_SetItem(self->directions, s, directionsAtS);
1103
1104 int x, y;
1105 for ( y = 1 ; y < (imgScaledAtS->yDim - 1) ; ++y) {
1106 for ( x = 1 ; x < (imgScaledAtS->xDim - 1) ; ++x) {
1107 // gradient magnitude m
1108 magnitudeAtS->values[x][y] = sqrt (
1109 pow ((double)(imgScaledAtS->values[x + 1][y] -
1110 imgScaledAtS->values[x - 1][y]), 2.0) +
1111 pow ((double)(imgScaledAtS->values[x][y + 1] -
1112 imgScaledAtS->values[x][y - 1]), 2.0));
1113
1114 // gradient direction theta
1115 directionsAtS->values[x][y] = atan2
1116 (imgScaledAtS->values[x][y + 1] - imgScaledAtS->values[x][y - 1],
1117 imgScaledAtS->values[x + 1][y] - imgScaledAtS->values[x - 1][y]);
1118 }
1119 }
1120 }
1121 }
1122
DScaleSpace_ClearMagnitudeAndDirectionMaps(DScaleSpace * self)1123 void DScaleSpace_ClearMagnitudeAndDirectionMaps (DScaleSpace* self)
1124 {
1125 // TODO: this should free up temporary memory, but this results in a heap corruption ???
1126
1127 //ArrayList_delete(self->magnitudes, NULL/*ImageMap_delete*/);
1128 //ArrayList_delete(self->directions, NULL/*ImageMap_delete*/);
1129 //self->magnitudes = self->directions = NULL;
1130 }
1131
1132 // Build a set of Difference-of-Gaussian (DoG) maps from the gaussian
1133 // blurred images.
1134 // This method has to be called after BuildGaussianMaps has completed.
DScaleSpace_BuildDiffMaps(DScaleSpace * self)1135 void DScaleSpace_BuildDiffMaps (DScaleSpace* self)
1136 {
1137 // Generate DoG maps. The maps are organized like this:
1138 // 0: D(sigma)
1139 // 1: D(k * sigma)
1140 // 2: D(k^2 * sigma)
1141 // ...
1142 // s: D(k^s * sigma) = D(2 * sigma)
1143 // s+1: D(k * 2 * sigma)
1144 //
1145 // So, we can start peak searching at 1 to s, and have a DoG map into
1146 // each direction.
1147 self->spaces = ArrayList_new(ArrayList_Count(self->imgScaled) - 1, ImageMap_delete);
1148
1149 // After the loop completes, we have used (s + 1) maps, yielding s
1150 // D-gaussian maps in the range of sigma to 2*sigma, as k^s = 2, which
1151 // is defined as one octave.
1152 int sn;
1153 for ( sn = 0 ; sn < ArrayList_Count(self->spaces) ; ++sn) {
1154 // XXX: order correct? It should not really matter as we search
1155 // for both minimums and maximums, but still, what is recommended?
1156 // (otherwise maybe the gradient directions are inverted?)
1157 ImageMap* imgScaledAtSnPlus1 = (ImageMap*) ArrayList_GetItem(self->imgScaled, sn+1);
1158 ImageMap* imgScaledAtSn = (ImageMap*) ArrayList_GetItem(self->imgScaled, sn);
1159 ImageMap* imgDiff = (ImageMap*) ImageMap_Sub(imgScaledAtSnPlus1, imgScaledAtSn);
1160 ArrayList_SetItem(self->spaces, sn, imgDiff);
1161 }
1162 }
1163
1164 // Incrementally blur the input image first so it reaches the next octave.
DScaleSpace_BuildGaussianMaps(DScaleSpace * self,ImageMap * first,double firstScale,int scales,double sigma)1165 void DScaleSpace_BuildGaussianMaps (DScaleSpace* self, ImageMap* first, double firstScale,
1166 int scales, double sigma)
1167 {
1168 // We need one more gaussian blurred image than the number of DoG
1169 // maps. But for the minima/maxima pixel search, we need two more. See
1170 // BuildDiffMaps.
1171 self->imgScaled = ArrayList_new(scales + 1 + 1 + 1, ImageMap_delete);
1172 self->basePixScale = firstScale;
1173
1174 // Ln1(x, y, k^{p+1}) = G(x, y, \sqrt{k^2-1}) * Ln0(x, y, k^p).
1175 ImageMap* prev = first;
1176 ArrayList_SetItem(self->imgScaled, 0, first);
1177
1178 /* Many thanks to Liu for this explanation and fix:
1179 *
1180 * Gaussian G(sigma), with relation
1181 * G(sigma_1) * G(sigma_2) = G(sqrt(sigma_1^2 + * sigma_2^2))
1182 *
1183 * Then, we have:
1184 *
1185 * G(k^{p+1}) = G(k^p) * G(sigma),
1186 * and our goal is to compute every iterations sigma value so self
1187 * equation iteratively produces the next level. Hence:
1188 *
1189 * sigma = \sqrt{\left(k^{p+1}\right)^2 - \left(k^p\right)^2}
1190 * = \sqrt{k^{2p+2} - k^{2p}}
1191 * = \sqrt{k^2p (k^2 - 1)}
1192 * = k^p \sqrt{k^2 - 1}
1193 *
1194 * In the code below, 'w' is the running k^p, where p increases by one
1195 * each iteration. kTerm is the constant \sqrt{k^2 - 1} term.
1196 */
1197 double w = sigma;
1198 double kTerm = sqrt (pow (DScaleSpace_SToK (scales), 2.0) - 1.0);
1199 int scI;
1200 for (scI = 1 ; scI < ArrayList_Count(self->imgScaled) ; ++scI) {
1201 GaussianConvolution* gauss = GaussianConvolution_new1 (w * kTerm);
1202 prev = GaussianConvolution_Convolve (gauss, prev);
1203 GaussianConvolution_delete(gauss);
1204 ArrayList_SetItem(self->imgScaled, scI, prev);
1205 w *= DScaleSpace_SToK (scales);
1206 }
1207
1208 }
1209
1210
ScalePoint_new0()1211 ScalePoint* ScalePoint_new0 ()
1212 {
1213 ScalePoint * self = (ScalePoint*)malloc(sizeof(ScalePoint));
1214 self->local = NULL;
1215 return self;
1216 }
1217
ScalePoint_delete(ScalePoint * self)1218 void ScalePoint_delete (ScalePoint* self)
1219 {
1220 if (self) {
1221 PointLocalInformation_delete(self->local);
1222 free(self);
1223 }
1224 }
1225
ScalePoint_new3(int x,int y,int level)1226 ScalePoint* ScalePoint_new3 (int x, int y, int level)
1227 {
1228 ScalePoint* self = ScalePoint_new0();
1229 self->x = x;
1230 self->y = y;
1231 self->level = level;
1232 self->local = NULL;
1233 return self;
1234 }
1235
PointLocalInformation_new0()1236 PointLocalInformation* PointLocalInformation_new0 ()
1237 {
1238 PointLocalInformation* self = (PointLocalInformation*)malloc(sizeof(PointLocalInformation));
1239 return self;
1240 }
1241
PointLocalInformation_delete(PointLocalInformation * self)1242 void PointLocalInformation_delete (PointLocalInformation* self)
1243 {
1244 if (self) {
1245 free( self );
1246 }
1247 }
1248
PointLocalInformation_new3(double fineS,double fineX,double fineY)1249 PointLocalInformation* PointLocalInformation_new3 (double fineS, double fineX, double fineY)
1250 {
1251 PointLocalInformation* self = PointLocalInformation_new0();
1252 self->fineX = fineX;
1253 self->fineY = fineY;
1254 self->scaleAdjust = fineS;
1255 self->dValue = 0;
1256 return self;
1257 }
1258
1259
1260
1261
Keypoint_FVGet(Keypoint * self,int xI,int yI,int oI)1262 double Keypoint_FVGet (Keypoint* self, int xI, int yI, int oI)
1263 {
1264 return (self->featureVector[(xI * self->yDim * self->oDim) + (yI * self->oDim) + oI]);
1265 }
Keypoint_FVSet(Keypoint * self,int xI,int yI,int oI,double value)1266 void Keypoint_FVSet (Keypoint* self, int xI, int yI, int oI, double value)
1267 {
1268 self->featureVector[(xI * self->yDim * self->oDim) + (yI * self->oDim) + oI] = value;
1269 }
1270
Keypoint_FVLinearDim(Keypoint * self)1271 int Keypoint_FVLinearDim(Keypoint* self) {
1272 return (self->featureVectorLength);
1273 }
1274
Keypoint_FVLinearGet(Keypoint * self,int idx)1275 double Keypoint_FVLinearGet (Keypoint* self, int idx)
1276 {
1277 return self->featureVector[idx];
1278 }
1279
Keypoint_FVLinearSet(Keypoint * self,int idx,double value)1280 void Keypoint_FVLinearSet (Keypoint* self, int idx, double value)
1281 {
1282 self->featureVector[idx] = value;
1283 }
1284
Keypoint_CreateLinearVector(Keypoint * self,int dim)1285 void Keypoint_CreateLinearVector (Keypoint* self, int dim)
1286 {
1287 self->featureVector = (double*)malloc( sizeof(double)*dim);
1288 self->featureVectorLength = dim;
1289 }
1290
Keypoint_CreateVector(Keypoint * self,int xDim,int yDim,int oDim)1291 void Keypoint_CreateVector (Keypoint* self, int xDim, int yDim, int oDim)
1292 {
1293 self->hasFV = true;
1294 self->xDim = xDim;
1295 self->yDim = yDim;
1296 self->oDim = oDim;
1297 self->featureVectorLength = yDim * xDim * oDim;
1298 self->featureVector = (double*)calloc(self->featureVectorLength, sizeof(double));
1299 }
1300
Keypoint_new0()1301 Keypoint* Keypoint_new0() {
1302 Keypoint* self = (Keypoint*)malloc(sizeof(Keypoint));
1303 self->featureVector = NULL;
1304 self->hasFV = false;
1305 self->featureVectorLength = 0;
1306 return self;
1307 }
1308
Keypoint_delete(Keypoint * self)1309 void Keypoint_delete(Keypoint* self) {
1310 if (self) {
1311 if (self->featureVector) {
1312 free(self->featureVector);
1313 }
1314 free(self);
1315 }
1316 }
1317
1318 // Keypoint constructor.
1319 //
1320 // image: The smoothed gaussian image the keypoint was located in.
1321 // x, y: The subpixel level coordinates of the keypoint.
1322 // imgScale: The scale of the gaussian image, with 1.0 being the original
1323 // detail scale (source image), and doubling at each octave.
1324 // kpScale: The scale of the keypoint.
1325 // orientation: Orientation degree in the range of [-PI ; PI] of the
1326 // keypoint.
1327 //
1328 // First add a keypoint, then use 'MakeDescriptor' to generate the local
1329 // image descriptor for this keypoint.
Keypoint_new(ImageMap * image,double x,double y,double imgScale,double kpScale,double orientation)1330 Keypoint* Keypoint_new (ImageMap* image, double x, double y, double imgScale,
1331 double kpScale, double orientation)
1332 {
1333 Keypoint* self = Keypoint_new0();
1334 self->image = image;
1335 self->x = x;
1336 self->y = y;
1337 self->imgScale = imgScale;
1338 self->scale = kpScale;
1339 self->orientation = orientation;
1340 self->hasFV = false;
1341 self->featureVector = NULL;
1342 return self;
1343 }
1344