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