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