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