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