1 : /******************************************************************************
2 : * Project: GDAL
3 : * Purpose: Correlator - GDALSimpleSURF and GDALFeaturePoint classes.
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 : #include "gdal_simplesurf.h"
29 :
30 : CPL_CVSID("$Id");
31 :
32 : /************************************************************************/
33 : /* ==================================================================== */
34 : /* GDALFeaturePoint */
35 : /* ==================================================================== */
36 : /************************************************************************/
37 :
38 0 : GDALFeaturePoint::GDALFeaturePoint()
39 : {
40 0 : nX = -1;
41 0 : nY = -1;
42 0 : nScale = -1;
43 0 : nRadius = -1;
44 0 : nSign = -1;
45 :
46 0 : padfDescriptor = new double[DESC_SIZE];
47 0 : }
48 :
49 0 : GDALFeaturePoint::GDALFeaturePoint(const GDALFeaturePoint& fp)
50 : {
51 0 : nX = fp.nX;
52 0 : nY = fp.nY;
53 0 : nScale = fp.nScale;
54 0 : nRadius = fp.nRadius;
55 0 : nSign = fp.nSign;
56 :
57 0 : padfDescriptor = new double[DESC_SIZE];
58 0 : for (int i = 0; i < DESC_SIZE; i++)
59 0 : padfDescriptor[i] = fp.padfDescriptor[i];
60 0 : }
61 :
62 0 : GDALFeaturePoint::GDALFeaturePoint(int nX, int nY,
63 0 : int nScale, int nRadius, int nSign)
64 : {
65 0 : this->nX = nX;
66 0 : this->nY = nY;
67 0 : this->nScale = nScale;
68 0 : this->nRadius = nRadius;
69 0 : this->nSign = nSign;
70 :
71 0 : this->padfDescriptor = new double[DESC_SIZE];
72 0 : }
73 :
74 0 : GDALFeaturePoint& GDALFeaturePoint::operator=(const GDALFeaturePoint& point)
75 : {
76 0 : if (this != &point)
77 : {
78 0 : nX = point.nX;
79 0 : nY = point.nY;
80 0 : nScale = point.nScale;
81 0 : nRadius = point.nRadius;
82 0 : nSign = point.nSign;
83 :
84 : //Free memory
85 0 : delete[] padfDescriptor;
86 :
87 : //Copy descriptor values
88 0 : padfDescriptor = new double[DESC_SIZE];
89 0 : for (int i = 0; i < DESC_SIZE; i++)
90 0 : padfDescriptor[i] = point.padfDescriptor[i];
91 : }
92 :
93 0 : return *this;
94 : }
95 :
96 0 : int GDALFeaturePoint::GetX() { return nX; }
97 0 : void GDALFeaturePoint::SetX(int nX) { this->nX = nX; }
98 :
99 0 : int GDALFeaturePoint::GetY() { return nY; }
100 0 : void GDALFeaturePoint::SetY(int nY) { this->nY = nY; }
101 :
102 0 : int GDALFeaturePoint::GetScale() { return nScale; }
103 0 : void GDALFeaturePoint::SetScale(int nScale) { this->nScale = nScale; }
104 :
105 0 : int GDALFeaturePoint::GetRadius() { return nRadius; }
106 0 : void GDALFeaturePoint::SetRadius(int nRadius) { this->nRadius = nRadius; }
107 :
108 0 : int GDALFeaturePoint::GetSign() { return nSign; }
109 0 : void GDALFeaturePoint::SetSign(int nSign) { this->nSign = nSign; }
110 :
111 0 : double& GDALFeaturePoint::operator [] (int nIndex)
112 : {
113 0 : if (nIndex < 0 || nIndex >= DESC_SIZE)
114 : {
115 : CPLError(CE_Failure, CPLE_AppDefined,
116 0 : "Descriptor index is out of range");
117 : }
118 :
119 0 : return padfDescriptor[nIndex];
120 : }
121 :
122 0 : GDALFeaturePoint::~GDALFeaturePoint() {
123 0 : delete[] padfDescriptor;
124 0 : }
125 :
126 : /************************************************************************/
127 : /* ==================================================================== */
128 : /* GDALSimpleSurf */
129 : /* ==================================================================== */
130 : /************************************************************************/
131 :
132 0 : GDALSimpleSURF::GDALSimpleSURF(int nOctaveStart, int nOctaveEnd)
133 : {
134 0 : this->octaveStart = nOctaveStart;
135 0 : this->octaveEnd = nOctaveEnd;
136 :
137 : // Initialize Octave map with custom range
138 0 : poOctMap = new GDALOctaveMap(octaveStart, octaveEnd);
139 0 : }
140 :
141 0 : CPLErr GDALSimpleSURF::ConvertRGBToLuminosity(
142 : GDALRasterBand *red, GDALRasterBand *green, GDALRasterBand *blue,
143 : int nXSize, int nYSize, double **padfImg, int nHeight, int nWidth)
144 : {
145 0 : const double forRed = 0.21;
146 0 : const double forGreen = 0.72;
147 0 : const double forBlue = 0.07;
148 :
149 0 : if (red == NULL || green == NULL || blue == NULL)
150 : {
151 : CPLError(CE_Failure, CPLE_AppDefined,
152 0 : "Raster bands are not specified");
153 0 : return CE_Failure;
154 : }
155 :
156 0 : if (nXSize > red->GetXSize() || nYSize > red->GetYSize())
157 : {
158 : CPLError(CE_Failure, CPLE_AppDefined,
159 0 : "Red band has less size than has been requested");
160 0 : return CE_Failure;
161 : }
162 :
163 0 : if (padfImg == NULL)
164 : {
165 0 : CPLError(CE_Failure, CPLE_AppDefined, "Buffer isn't specified");
166 0 : return CE_Failure;
167 : }
168 :
169 0 : GDALDataType eRedType = red->GetRasterDataType();
170 0 : GDALDataType eGreenType = green->GetRasterDataType();
171 0 : GDALDataType eBlueType = blue->GetRasterDataType();
172 :
173 0 : int dataRedSize = GDALGetDataTypeSize(eRedType) / 8;
174 0 : int dataGreenSize = GDALGetDataTypeSize(eGreenType) / 8;
175 0 : int dataBlueSize = GDALGetDataTypeSize(eBlueType) / 8;
176 :
177 0 : void *paRedLayer = CPLMalloc(dataRedSize * nWidth * nHeight);
178 0 : void *paGreenLayer = CPLMalloc(dataGreenSize * nWidth * nHeight);
179 0 : void *paBlueLayer = CPLMalloc(dataBlueSize * nWidth * nHeight);
180 :
181 0 : red->RasterIO(GF_Read, 0, 0, nXSize, nYSize, paRedLayer, nWidth, nHeight, eRedType, 0, 0);
182 0 : green->RasterIO(GF_Read, 0, 0, nXSize, nYSize, paGreenLayer, nWidth, nHeight, eGreenType, 0, 0);
183 0 : blue->RasterIO(GF_Read, 0, 0, nXSize, nYSize, paBlueLayer, nWidth, nHeight, eBlueType, 0, 0);
184 :
185 0 : double maxValue = 255.0;
186 0 : for (int row = 0; row < nHeight; row++)
187 0 : for (int col = 0; col < nWidth; col++)
188 : {
189 : // Get RGB values
190 0 : double dfRedVal = SRCVAL(paRedLayer, eRedType,
191 : nWidth * row + col * dataRedSize);
192 0 : double dfGreenVal = SRCVAL(paGreenLayer, eGreenType,
193 : nWidth * row + col * dataGreenSize);
194 0 : double dfBlueVal = SRCVAL(paBlueLayer, eBlueType,
195 : nWidth * row + col * dataBlueSize);
196 : // Compute luminosity value
197 0 : padfImg[row][col] = (
198 : dfRedVal * forRed +
199 : dfGreenVal * forGreen +
200 0 : dfBlueVal * forBlue) / maxValue;
201 : }
202 :
203 0 : CPLFree(paRedLayer);
204 0 : CPLFree(paGreenLayer);
205 0 : CPLFree(paBlueLayer);
206 :
207 0 : return CE_None;
208 : }
209 :
210 : std::vector<GDALFeaturePoint>*
211 0 : GDALSimpleSURF::ExtractFeaturePoints(GDALIntegralImage *poImg,
212 : double dfThreshold)
213 : {
214 : std::vector<GDALFeaturePoint>* poCollection =
215 0 : new std::vector<GDALFeaturePoint>();
216 :
217 : //Calc Hessian values for layers
218 0 : poOctMap->ComputeMap(poImg);
219 :
220 : //Search for exremum points
221 0 : for (int oct = octaveStart; oct <= octaveEnd; oct++)
222 : {
223 0 : for (int k = 0; k < GDALOctaveMap::INTERVALS - 2; k++)
224 : {
225 0 : GDALOctaveLayer *bot = poOctMap->pMap[oct - 1][k];
226 0 : GDALOctaveLayer *mid = poOctMap->pMap[oct - 1][k + 1];
227 0 : GDALOctaveLayer *top = poOctMap->pMap[oct - 1][k + 2];
228 :
229 0 : for (int i = 0; i < mid->height; i++)
230 : {
231 0 : for (int j = 0; j < mid->width; j++)
232 : {
233 0 : if (poOctMap->PointIsExtremum(i, j, bot, mid, top, dfThreshold))
234 : {
235 : GDALFeaturePoint oFP(j, i, mid->scale,
236 0 : mid->radius, mid->signs[i][j]);
237 0 : SetDescriptor(&oFP, poImg);
238 0 : poCollection->push_back(oFP);
239 : }
240 : }
241 : }
242 : }
243 : }
244 :
245 0 : return poCollection;
246 : }
247 :
248 :
249 0 : double GDALSimpleSURF::GetEuclideanDistance(
250 : GDALFeaturePoint &firstPoint, GDALFeaturePoint &secondPoint)
251 : {
252 0 : double sum = 0;
253 :
254 0 : for (int i = 0; i < GDALFeaturePoint::DESC_SIZE; i++)
255 0 : sum += (firstPoint[i] - secondPoint[i]) * (firstPoint[i] - secondPoint[i]);
256 :
257 0 : return sqrt(sum);
258 : }
259 :
260 0 : void GDALSimpleSURF::NormalizeDistances(std::list<MatchedPointPairInfo> *poList)
261 : {
262 0 : double max = 0;
263 :
264 0 : std::list<MatchedPointPairInfo>::iterator i;
265 0 : for (i = poList->begin(); i != poList->end(); i++)
266 0 : if ((*i).euclideanDist > max)
267 0 : max = (*i).euclideanDist;
268 :
269 0 : if (max != 0)
270 : {
271 0 : for (i = poList->begin(); i != poList->end(); i++)
272 0 : (*i).euclideanDist /= max;
273 : }
274 0 : }
275 :
276 0 : void GDALSimpleSURF::SetDescriptor(
277 : GDALFeaturePoint *poPoint, GDALIntegralImage *poImg)
278 : {
279 : // Affects to the descriptor area
280 0 : const int haarScale = 20;
281 :
282 : // Side of the Haar wavelet
283 0 : int haarFilterSize = 2 * poPoint->GetScale();
284 :
285 : // Length of the side of the descriptor area
286 0 : int descSide = haarScale * poPoint->GetScale();
287 :
288 : // Side of the quadrant in 4x4 grid
289 0 : int quadStep = descSide / 4;
290 :
291 : // Side of the sub-quadrant in 5x5 regular grid of quadrant
292 0 : int subQuadStep = quadStep / 5;
293 :
294 0 : int leftTop_row = poPoint->GetY() - (descSide / 2);
295 0 : int leftTop_col = poPoint->GetX() - (descSide / 2);
296 :
297 0 : int count = 0;
298 :
299 0 : for (int r = leftTop_row; r < leftTop_row + descSide; r += quadStep)
300 0 : for (int c = leftTop_col; c < leftTop_col + descSide; c += quadStep)
301 : {
302 0 : double dx = 0;
303 0 : double dy = 0;
304 0 : double abs_dx = 0;
305 0 : double abs_dy = 0;
306 :
307 0 : for (int sub_r = r; sub_r < r + quadStep; sub_r += subQuadStep)
308 0 : for (int sub_c = c; sub_c < c + quadStep; sub_c += subQuadStep)
309 : {
310 : // Approximate center of sub quadrant
311 0 : int cntr_r = sub_r + subQuadStep / 2;
312 0 : int cntr_c = sub_c + subQuadStep / 2;
313 :
314 : // Left top point for Haar wavelet computation
315 0 : int cur_r = cntr_r - haarFilterSize / 2;
316 0 : int cur_c = cntr_c - haarFilterSize / 2;
317 :
318 : // Gradients
319 0 : double cur_dx = poImg->HaarWavelet_X(cur_r, cur_c, haarFilterSize);
320 0 : double cur_dy = poImg->HaarWavelet_Y(cur_r, cur_c, haarFilterSize);
321 :
322 0 : dx += cur_dx;
323 0 : dy += cur_dy;
324 0 : abs_dx += fabs(cur_dx);
325 0 : abs_dy += fabs(cur_dy);
326 : }
327 :
328 : // Fills point's descriptor
329 0 : (*poPoint)[count++] = dx;
330 0 : (*poPoint)[count++] = dy;
331 0 : (*poPoint)[count++] = abs_dx;
332 0 : (*poPoint)[count++] = abs_dy;
333 : }
334 0 : }
335 :
336 : /**
337 : * Find corresponding points (equal points in two collections).
338 : *
339 : * @param poMatched Resulting collection for matched points
340 : * @param poFirstCollection Points on the first image
341 : * @param poSecondCollection Points on the second image
342 : * @param dfThreshold Value from 0 to 1. Threshold affects to number of
343 : * matched points. If threshold is higher, amount of corresponding
344 : * points is larger, and vice versa
345 : *
346 : * @note Typical threshold's value is 0,1. BUT it's a very approximate guide.
347 : * It can be 0,001 or even 1. This threshold provides direct adjustment
348 : * of point matching.
349 : * NOTICE that if threshold is lower, matches are more robust and correct, but number of
350 : * matched points is smaller. Therefore if algorithm performs many false
351 : * detections and produces bad results, reduce threshold.
352 : * Otherwise, if algorithm finds nothing, increase threshold.
353 : *
354 : * @return CE_None or CE_Failure if error occurs.
355 : */
356 :
357 0 : CPLErr GDALSimpleSURF::MatchFeaturePoints(
358 : std::vector<GDALFeaturePoint*> *poMatchPairs,
359 : std::vector<GDALFeaturePoint> *poFirstCollect,
360 : std::vector<GDALFeaturePoint> *poSecondCollect,
361 : double dfThreshold)
362 : {
363 : /* -------------------------------------------------------------------- */
364 : /* Validate parameters. */
365 : /* -------------------------------------------------------------------- */
366 0 : if (poMatchPairs == NULL)
367 : {
368 : CPLError(CE_Failure, CPLE_AppDefined,
369 0 : "Matched points colection isn't specified");
370 0 : return CE_Failure;
371 : }
372 :
373 0 : if (poFirstCollect == NULL || poSecondCollect == NULL)
374 : {
375 : CPLError(CE_Failure, CPLE_AppDefined,
376 0 : "Feature point collections are not specified");
377 0 : return CE_Failure;
378 : }
379 :
380 : /* ==================================================================== */
381 : /* Matching algorithm. */
382 : /* ==================================================================== */
383 : // Affects to false matching pruning
384 0 : const double ratioThreshold = 0.8;
385 :
386 0 : int len_1 = poFirstCollect->size();
387 0 : int len_2 = poSecondCollect->size();
388 :
389 0 : int minLength = (len_1 < len_2) ? len_1 : len_2;
390 :
391 : // Temporary pointers. Used to swap collections
392 : std::vector<GDALFeaturePoint> *p_1;
393 : std::vector<GDALFeaturePoint> *p_2;
394 :
395 0 : bool isSwap = false;
396 :
397 : // Assign p_1 - collection with minimal number of points
398 0 : if (minLength == len_2)
399 : {
400 0 : p_1 = poSecondCollect;
401 0 : p_2 = poFirstCollect;
402 :
403 0 : int tmp = 0;
404 0 : tmp = len_1;
405 0 : len_1 = len_2;
406 0 : len_2 = tmp;
407 0 : isSwap = true;
408 : }
409 : else
410 : {
411 : // Assignment 'as is'
412 0 : p_1 = poFirstCollect;
413 0 : p_2 = poSecondCollect;
414 0 : isSwap = false;
415 : }
416 :
417 : // Stores matched point indexes and
418 : // their euclidean distances
419 : std::list<MatchedPointPairInfo> *poPairInfoList =
420 0 : new std::list<MatchedPointPairInfo>();
421 :
422 : // Flags that points in the 2nd collection are matched or not
423 0 : bool *alreadyMatched = new bool[len_2];
424 0 : for (int i = 0; i < len_2; i++)
425 0 : alreadyMatched[i] = false;
426 :
427 0 : for (int i = 0; i < len_1; i++)
428 : {
429 : // Distance to the nearest point
430 0 : double bestDist = -1;
431 : // Index of the nearest point in p_2 collection
432 0 : int bestIndex = -1;
433 :
434 : // Distance to the 2nd nearest point
435 0 : double bestDist_2 = -1;
436 :
437 : // Find the nearest and 2nd nearest points
438 0 : for (int j = 0; j < len_2; j++)
439 0 : if (!alreadyMatched[j])
440 0 : if (p_1->at(i).GetSign() ==
441 : p_2->at(j).GetSign())
442 : {
443 : // Get distance between two feature points
444 : double curDist = GetEuclideanDistance(
445 0 : p_1->at(i), p_2->at(j));
446 :
447 0 : if (bestDist == -1)
448 : {
449 0 : bestDist = curDist;
450 0 : bestIndex = j;
451 : }
452 : else
453 : {
454 0 : if (curDist < bestDist)
455 : {
456 0 : bestDist = curDist;
457 0 : bestIndex = j;
458 : }
459 : }
460 :
461 : // Findes the 2nd nearest point
462 0 : if (bestDist_2 < 0)
463 0 : bestDist_2 = curDist;
464 : else
465 0 : if (curDist > bestDist && curDist < bestDist_2)
466 0 : bestDist_2 = curDist;
467 : }
468 : /* -------------------------------------------------------------------- */
469 : /* False matching pruning. */
470 : /* If ratio bestDist to bestDist_2 greater than 0.8 => */
471 : /* consider as false detection. */
472 : /* Otherwise, add points as matched pair. */
473 : /*----------------------------------------------------------------------*/
474 0 : if (bestDist_2 > 0 && bestDist >= 0)
475 0 : if (bestDist / bestDist_2 < ratioThreshold)
476 : {
477 0 : MatchedPointPairInfo info(i, bestIndex, bestDist);
478 0 : poPairInfoList->push_back(info);
479 0 : alreadyMatched[bestIndex] = true;
480 : }
481 : }
482 :
483 :
484 : /* -------------------------------------------------------------------- */
485 : /* Pruning based on the provided threshold */
486 : /* -------------------------------------------------------------------- */
487 :
488 0 : NormalizeDistances(poPairInfoList);
489 :
490 0 : std::list<MatchedPointPairInfo>::const_iterator iter;
491 0 : for (iter = poPairInfoList->begin(); iter != poPairInfoList->end(); iter++)
492 : {
493 0 : if ((*iter).euclideanDist <= dfThreshold)
494 : {
495 0 : int i_1 = (*iter).ind_1;
496 0 : int i_2 = (*iter).ind_2;
497 :
498 : // Add copies into MatchedCollection
499 0 : if(!isSwap)
500 : {
501 0 : poMatchPairs->push_back( &(p_1->at(i_1)) );
502 0 : poMatchPairs->push_back( &(p_2->at(i_2)) );
503 : }
504 : else
505 : {
506 0 : poMatchPairs->push_back( &(p_2->at(i_2)) );
507 0 : poMatchPairs->push_back( &(p_1->at(i_1)) );
508 : }
509 : }
510 : }
511 :
512 : // Clean up
513 0 : delete[] alreadyMatched;
514 :
515 0 : return CE_None;
516 : }
517 :
518 0 : GDALSimpleSURF::~GDALSimpleSURF()
519 : {
520 0 : delete poOctMap;
521 0 : }
522 :
|
