Line data Source code
1 : /******************************************************************************
2 : * $Id$
3 : * Project: GDAL
4 : * Purpose: Correlator
5 : * Author: Andrew Migal, migal.drew@gmail.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2012, Andrew Migal
9 : *
10 : * Permission is hereby granted, free of charge, to any person obtaining a
11 : * copy of this software and associated documentation files (the "Software"),
12 : * to deal in the Software without restriction, including without limitation
13 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
14 : * and/or sell copies of the Software, and to permit persons to whom the
15 : * Software is furnished to do so, subject to the following conditions:
16 : *
17 : * The above copyright notice and this permission notice shall be included
18 : * in all copies or substantial portions of the Software.
19 : *
20 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26 : * DEALINGS IN THE SOFTWARE.
27 : ****************************************************************************/
28 :
29 : /**
30 : * @file
31 : * @author Andrew Migal migal.drew@gmail.com
32 : * @brief Class for searching corresponding points on images.
33 : */
34 :
35 : #ifndef GDALSIMPLESURF_H_
36 : #define GDALSIMPLESURF_H_
37 :
38 : #include "gdal_priv.h"
39 : #include "cpl_conv.h"
40 : #include <list>
41 :
42 : /**
43 : * @brief Class of "feature point" in raster. Used by SURF-based algorithm.
44 : *
45 : * @details This point presents coordinates of distinctive pixel in image.
46 : * In computer vision, feature points - the most "strong" and "unique"
47 : * pixels (or areas) in picture, which can be distinguished from others.
48 : * For more details, see FAST corner detector, SIFT, SURF and similar
49 : * algorithms.
50 : */
51 : class GDALFeaturePoint
52 : {
53 : public:
54 : /**
55 : * Standard constructor. Initializes all parameters with negative numbers
56 : * and allocates memory for descriptor.
57 : */
58 : GDALFeaturePoint();
59 :
60 : /**
61 : * Copy constructor
62 : * @param fp Copied instance of GDALFeaturePoint class
63 : */
64 : GDALFeaturePoint(const GDALFeaturePoint &fp);
65 :
66 : /**
67 : * Create instance of GDALFeaturePoint class
68 : *
69 : * @param nX X-coordinate (pixel)
70 : * @param nY Y-coordinate (line)
71 : * @param nScale Scale which contains this point (2, 4, 8, 16 and so on)
72 : * @param nRadius Half of the side of descriptor area
73 : * @param nSign Sign of Hessian determinant for this point
74 : *
75 : * @note This constructor normally is invoked by SURF-based algorithm,
76 : * which provides all necessary parameters.
77 : */
78 : GDALFeaturePoint(int nX, int nY, int nScale, int nRadius, int nSign);
79 : virtual ~GDALFeaturePoint();
80 :
81 : /** Assignment operator */
82 : GDALFeaturePoint &operator=(const GDALFeaturePoint &point);
83 :
84 : /**
85 : * Provide access to point's descriptor.
86 : *
87 : * @param nIndex Position of descriptor's value.
88 : * nIndex should be within range from 0 to DESC_SIZE (in current version -
89 : * 64)
90 : *
91 : * @return Reference to value of descriptor in 'nIndex' position.
92 : * If index is out of range then behavior is undefined.
93 : */
94 : double &operator[](int nIndex);
95 :
96 : /**
97 : * Provide access to point's descriptor.
98 : *
99 : * @param nIndex Position of descriptor's value.
100 : * nIndex should be within range from 0 to DESC_SIZE (in current version -
101 : * 64)
102 : *
103 : * @return Reference to value of descriptor in 'nIndex' position.
104 : * If index is out of range then behavior is undefined.
105 : */
106 : double operator[](int nIndex) const;
107 :
108 : /** Descriptor length */
109 : static const int DESC_SIZE = 64;
110 :
111 : /**
112 : * Fetch X-coordinate (pixel) of point
113 : *
114 : * @return X-coordinate in pixels
115 : */
116 : int GetX() const;
117 :
118 : /**
119 : * Set X coordinate of point
120 : *
121 : * @param nX X coordinate in pixels
122 : */
123 : void SetX(int nX);
124 :
125 : /**
126 : * Fetch Y-coordinate (line) of point.
127 : *
128 : * @return Y-coordinate in pixels.
129 : */
130 : int GetY() const;
131 :
132 : /**
133 : * Set Y coordinate of point.
134 : *
135 : * @param nY Y coordinate in pixels.
136 : */
137 : void SetY(int nY);
138 :
139 : /**
140 : * Fetch scale of point.
141 : *
142 : * @return Scale for this point.
143 : */
144 : int GetScale() const;
145 :
146 : /**
147 : * Set scale of point.
148 : *
149 : * @param nScale Scale for this point.
150 : */
151 : void SetScale(int nScale);
152 :
153 : /**
154 : * Fetch radius of point.
155 : *
156 : * @return Radius for this point.
157 : */
158 : int GetRadius() const;
159 :
160 : /**
161 : * Set radius of point.
162 : *
163 : * @param nRadius Radius for this point.
164 : */
165 : void SetRadius(int nRadius);
166 :
167 : /**
168 : * Fetch sign of Hessian determinant of point.
169 : *
170 : * @return Sign for this point.
171 : */
172 : int GetSign() const;
173 :
174 : /**
175 : * Set sign of point.
176 : *
177 : * @param nSign Sign of Hessian determinant for this point.
178 : */
179 : void SetSign(int nSign);
180 :
181 : private:
182 : // Coordinates of point in image
183 : int nX;
184 : int nY;
185 : // --------------------
186 : int nScale;
187 : int nRadius;
188 : int nSign;
189 : // Descriptor array
190 : double *padfDescriptor;
191 : };
192 :
193 : /**
194 : * @author Andrew Migal migal.drew@gmail.com
195 : * @brief Integral image class (summed area table).
196 : * @details Integral image is a table for fast computing the sum of
197 : * values in rectangular subarea. In more detail, for 2-dimensional array
198 : * of numbers this class provides capability to get sum of values in
199 : * rectangular arbitrary area with any size in constant time.
200 : * Integral image is constructed from grayscale picture.
201 : */
202 0 : class GDALIntegralImage
203 : {
204 : CPL_DISALLOW_COPY_ASSIGN(GDALIntegralImage)
205 :
206 : public:
207 : GDALIntegralImage();
208 : virtual ~GDALIntegralImage();
209 :
210 : /**
211 : * Compute integral image for specified array. Result is stored internally.
212 : *
213 : * @param padfImg Pointer to 2-dimensional array of values
214 : * @param nHeight Number of rows in array
215 : * @param nWidth Number of columns in array
216 : */
217 : void Initialize(const double **padfImg, int nHeight, int nWidth);
218 :
219 : /**
220 : * Fetch value of specified position in integral image.
221 : *
222 : * @param nRow Row of this position
223 : * @param nCol Column of this position
224 : *
225 : * @return Value in specified position or zero if parameters are out of
226 : * range.
227 : */
228 : double GetValue(int nRow, int nCol);
229 :
230 : /**
231 : * Get sum of values in specified rectangular grid. Rectangle is constructed
232 : * from left top point.
233 : *
234 : * @param nRow Row of left top point of rectangle
235 : * @param nCol Column of left top point of rectangle
236 : * @param nWidth Width of rectangular area (number of columns)
237 : * @param nHeight Height of rectangular area (number of rows)
238 : *
239 : * @return Sum of values in specified grid.
240 : */
241 : double GetRectangleSum(int nRow, int nCol, int nWidth, int nHeight);
242 :
243 : /**
244 : * Get value of horizontal Haar wavelet in specified square grid.
245 : *
246 : * @param nRow Row of left top point of square
247 : * @param nCol Column of left top point of square
248 : * @param nSize Side of the square
249 : *
250 : * @return Value of horizontal Haar wavelet in specified square grid.
251 : */
252 : double HaarWavelet_X(int nRow, int nCol, int nSize);
253 :
254 : /**
255 : * Get value of vertical Haar wavelet in specified square grid.
256 : *
257 : * @param nRow Row of left top point of square
258 : * @param nCol Column of left top point of square
259 : * @param nSize Side of the square
260 : *
261 : * @return Value of vertical Haar wavelet in specified square grid.
262 : */
263 : double HaarWavelet_Y(int nRow, int nCol, int nSize);
264 :
265 : /**
266 : * Fetch height of integral image.
267 : *
268 : * @return Height of integral image (number of rows).
269 : */
270 : int GetHeight();
271 :
272 : /**
273 : * Fetch width of integral image.
274 : *
275 : * @return Width of integral image (number of columns).
276 : */
277 : int GetWidth();
278 :
279 : private:
280 : double **pMatrix = nullptr;
281 : int nWidth = 0;
282 : int nHeight = 0;
283 : };
284 :
285 : /**
286 : * @author Andrew Migal migal.drew@gmail.com
287 : * @brief Class for computation and storage of Hessian values in SURF-based
288 : * algorithm.
289 : *
290 : * @details SURF-based algorithm normally uses this class for searching
291 : * feature points on raster images. Class also contains traces of Hessian
292 : * matrices to provide fast computations.
293 : */
294 : class GDALOctaveLayer
295 : {
296 : CPL_DISALLOW_COPY_ASSIGN(GDALOctaveLayer)
297 :
298 : public:
299 : GDALOctaveLayer();
300 :
301 : /**
302 : * Create instance with provided parameters.
303 : *
304 : * @param nOctave Number of octave which contains this layer
305 : * @param nInterval Number of position in octave
306 : *
307 : * @note Normally constructor is invoked only by SURF-based algorithm.
308 : */
309 : GDALOctaveLayer(int nOctave, int nInterval);
310 : virtual ~GDALOctaveLayer();
311 :
312 : /**
313 : * Perform calculation of Hessian determinants and their signs
314 : * for specified integral image. Result is stored internally.
315 : *
316 : * @param poImg Integral image object, which provides all necessary
317 : * data for computation
318 : *
319 : * @note Normally method is invoked only by SURF-based algorithm.
320 : */
321 : void ComputeLayer(GDALIntegralImage *poImg);
322 :
323 : /**
324 : * Octave which contains this layer (1,2,3...)
325 : */
326 : int octaveNum;
327 : /**
328 : * Length of the side of filter
329 : */
330 : int filterSize;
331 : /**
332 : * Length of the border
333 : */
334 : int radius;
335 : /**
336 : * Scale for this layer
337 : */
338 : int scale;
339 : /**
340 : * Image width in pixels
341 : */
342 : int width;
343 : /**
344 : * Image height in pixels
345 : */
346 : int height;
347 : /**
348 : * Hessian values for image pixels
349 : */
350 : double **detHessians;
351 : /**
352 : * Hessian signs for speeded matching
353 : */
354 : int **signs;
355 : };
356 :
357 : /**
358 : * @author Andrew Migal migal.drew@gmail.com
359 : * @brief Class for handling octave layers in SURF-based algorithm.
360 : * @details Class contains OctaveLayers and provides capability to construct
361 : * octave space and distinguish feature points. Normally this class is used only
362 : * by SURF-based algorithm.
363 : */
364 : class GDALOctaveMap
365 : {
366 : CPL_DISALLOW_COPY_ASSIGN(GDALOctaveMap)
367 :
368 : public:
369 : /**
370 : * Create octave space. Octave numbers are start with one. (1, 2, 3, 4, ...
371 : * )
372 : *
373 : * @param nOctaveStart Number of bottom octave
374 : * @param nOctaveEnd Number of top octave. Should be equal or greater than
375 : * OctaveStart
376 : */
377 : GDALOctaveMap(int nOctaveStart, int nOctaveEnd);
378 : virtual ~GDALOctaveMap();
379 :
380 : /**
381 : * Calculate Hessian values for octave space
382 : * (for all stored octave layers) using specified integral image
383 : * @param poImg Integral image instance which provides necessary data
384 : * @see GDALOctaveLayer
385 : */
386 : void ComputeMap(GDALIntegralImage *poImg);
387 :
388 : /**
389 : * Method makes decision that specified point
390 : * in middle octave layer is maximum among all points
391 : * from 3x3x3 neighbourhood (surrounding points in
392 : * bottom, middle and top layers). Provided layers should be from the same
393 : * octave's interval. Detects feature points.
394 : *
395 : * @param row Row of point, which is candidate to be feature point
396 : * @param col Column of point, which is candidate to be feature point
397 : * @param bot Bottom octave layer
398 : * @param mid Middle octave layer
399 : * @param top Top octave layer
400 : * @param threshold Threshold for feature point recognition. Detected
401 : * feature point will have Hessian value greater than this provided
402 : * threshold.
403 : *
404 : * @return TRUE if candidate was evaluated as feature point or FALSE
405 : * otherwise.
406 : */
407 : static bool PointIsExtremum(int row, int col, GDALOctaveLayer *bot,
408 : GDALOctaveLayer *mid, GDALOctaveLayer *top,
409 : double threshold);
410 :
411 : /**
412 : * 2-dimensional array of octave layers
413 : */
414 : GDALOctaveLayer ***pMap;
415 :
416 : /**
417 : * Value for constructing internal octave space
418 : */
419 : static const int INTERVALS = 4;
420 :
421 : /**
422 : * Number of bottom octave
423 : */
424 : int octaveStart;
425 :
426 : /**
427 : * Number of top octave. Should be equal or greater than OctaveStart
428 : */
429 : int octaveEnd;
430 : };
431 :
432 : /**
433 : * @author Andrew Migal migal.drew@gmail.com
434 : * @brief Class for searching corresponding points on images.
435 : * @details Provides capability for detection feature points
436 : * and finding equal points on different images.
437 : * Class implements simplified version of SURF algorithm (Speeded Up Robust
438 : * Features). As original, this realization is scale invariant, but sensitive to
439 : * rotation. Images should have similar rotation angles (maximum difference is
440 : * up to 10-15 degrees), otherwise algorithm produces incorrect and very
441 : * unstable results.
442 : */
443 :
444 : class GDALSimpleSURF
445 : {
446 : private:
447 : /**
448 : * Class stores indexes of pair of point
449 : * and distance between them.
450 : */
451 : class MatchedPointPairInfo
452 : {
453 : public:
454 0 : MatchedPointPairInfo(int nInd_1, int nInd_2, double dfDist)
455 0 : : ind_1(nInd_1), ind_2(nInd_2), euclideanDist(dfDist)
456 : {
457 0 : }
458 :
459 : int ind_1;
460 : int ind_2;
461 : double euclideanDist;
462 : };
463 :
464 : CPL_DISALLOW_COPY_ASSIGN(GDALSimpleSURF)
465 :
466 : public:
467 : /**
468 : * Prepare class according to specified parameters. Octave numbers affects
469 : * to amount of detected points and their robustness.
470 : * Range between bottom and top octaves also affects to required time of
471 : * detection points (if range is large, algorithm should perform more
472 : * operations).
473 : * @param nOctaveStart Number of bottom octave. Octave numbers starts with
474 : * one
475 : * @param nOctaveEnd Number of top octave. Should be equal or greater than
476 : * OctaveStart
477 : *
478 : * @note
479 : * Every octave finds points with specific size. For small images
480 : * use small octave numbers, for high resolution - large.
481 : * For 1024x1024 images it's normal to use any octave numbers from range
482 : * 1-6. (for example, octave start - 1, octave end - 3, or octave start - 2,
483 : * octave end - 2.) For larger images, try 1-10 range or even higher. Pay
484 : * attention that number of detected point decreases quickly per octave for
485 : * particular image. Algorithm finds more points in case of small octave
486 : * numbers. If method detects nothing, reduce bottom bound of octave range.
487 : *
488 : * NOTICE that every octave requires time to compute. Use a little range
489 : * or only one octave if execution time is significant.
490 : */
491 : GDALSimpleSURF(int nOctaveStart, int nOctaveEnd);
492 : virtual ~GDALSimpleSURF();
493 :
494 : /**
495 : * Convert image with RGB channels to grayscale using "luminosity" method.
496 : * Result is used in SURF-based algorithm, but may be used anywhere where
497 : * grayscale images with nice contrast are required.
498 : *
499 : * @param red Image's red channel
500 : * @param green Image's green channel
501 : * @param blue Image's blue channel
502 : * @param nXSize Width of initial image
503 : * @param nYSize Height of initial image
504 : * @param padfImg Array for resulting grayscale image
505 : * @param nHeight Height of resulting image
506 : * @param nWidth Width of resulting image
507 : *
508 : * @return CE_None or CE_Failure if error occurs.
509 : */
510 : static CPLErr ConvertRGBToLuminosity(GDALRasterBand *red,
511 : GDALRasterBand *green,
512 : GDALRasterBand *blue, int nXSize,
513 : int nYSize, double **padfImg,
514 : int nHeight, int nWidth);
515 :
516 : /**
517 : * Find feature points using specified integral image.
518 : *
519 : * @param poImg Integral image to be used
520 : * @param dfThreshold Threshold for feature point recognition. Detected
521 : * feature point will have Hessian value greater than this provided
522 : * threshold.
523 : *
524 : * @note Typical threshold's value is 0,001. But this value
525 : * can be various in each case and depends on image's nature.
526 : * For example, value can be 0.002 or 0.005.
527 : * Fill free to experiment with it.
528 : * If threshold is high, than number of detected feature points is small,
529 : * and vice versa.
530 : */
531 : std::vector<GDALFeaturePoint> *
532 : ExtractFeaturePoints(GDALIntegralImage *poImg, double dfThreshold);
533 :
534 : /**
535 : * Find corresponding points (equal points in two collections).
536 : *
537 : * @param poMatchPairs Resulting collection for matched points
538 : * @param poFirstCollect Points on the first image
539 : * @param poSecondCollect Points on the second image
540 : * @param dfThreshold Value from 0 to 1. Threshold affects to number of
541 : * matched points. If threshold is higher, amount of corresponding
542 : * points is larger, and vice versa
543 : *
544 : * @note Typical threshold's value is 0,1. BUT it's a very approximate
545 : * guide. It can be 0,001 or even 1. This threshold provides direct
546 : * adjustment of point matching. NOTICE that if threshold is lower, matches
547 : * are more robust and correct, but number of matched points is smaller.
548 : * Therefore if algorithm performs many false detections and produces bad
549 : * results, reduce threshold. Otherwise, if algorithm finds nothing,
550 : * increase threshold.
551 : *
552 : * @return CE_None or CE_Failure if error occurs.
553 : */
554 : static CPLErr
555 : MatchFeaturePoints(std::vector<GDALFeaturePoint *> *poMatchPairs,
556 : std::vector<GDALFeaturePoint> *poFirstCollect,
557 : std::vector<GDALFeaturePoint> *poSecondCollect,
558 : double dfThreshold);
559 :
560 : private:
561 : /**
562 : * Compute euclidean distance between descriptors of two feature points.
563 : * It's used in comparison and matching of points.
564 : *
565 : * @param firstPoint First feature point to be compared
566 : * @param secondPoint Second feature point to be compared
567 : *
568 : * @return Euclidean distance between descriptors.
569 : */
570 : static double GetEuclideanDistance(const GDALFeaturePoint &firstPoint,
571 : const GDALFeaturePoint &secondPoint);
572 :
573 : /**
574 : * Set provided distance values to range from 0 to 1.
575 : *
576 : * @param poList List of distances to be normalized
577 : */
578 : static void NormalizeDistances(std::list<MatchedPointPairInfo> *poList);
579 :
580 : /**
581 : * Compute descriptor for specified feature point.
582 : *
583 : * @param poPoint Feature point instance
584 : * @param poImg image where feature point was found
585 : */
586 : static void SetDescriptor(GDALFeaturePoint *poPoint,
587 : GDALIntegralImage *poImg);
588 :
589 : private:
590 : int octaveStart;
591 : int octaveEnd;
592 : GDALOctaveMap *poOctMap;
593 : };
594 :
595 : #endif /* GDALSIMPLESURF_H_ */
|
