Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: GDAL
4 : * Purpose: GDAL Wrapper for image matching via correlation algorithm.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : * Author: Andrew Migal, migal.drew@gmail.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2012, Frank Warmerdam
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "gdal_alg.h"
15 : #include "gdal_simplesurf.h"
16 :
17 : //! @cond Doxygen_Suppress
18 :
19 : //! @endcond
20 :
21 : // TODO(schwehr): What? This below: "0,001"
22 :
23 : /**
24 : * @file
25 : * @author Andrew Migal migal.drew@gmail.com
26 : * @brief Algorithms for searching corresponding points on images.
27 : * @details This implementation is based on an simplified version
28 : * of SURF algorithm (Speeded Up Robust Features).
29 : * Provides capability for detection feature points
30 : * and finding equal points on different images.
31 : * As original, this realization is scale invariant, but sensitive to rotation.
32 : * Images should have similar rotation angles (maximum difference is up to 10-15
33 : * degrees), otherwise algorithm produces incorrect and very unstable results.
34 : */
35 :
36 : /**
37 : * Detect feature points on provided image. Please carefully read documentation
38 : * below.
39 : *
40 : * @param poDataset Image on which feature points will be detected
41 : * @param panBands Array of 3 raster bands numbers, for Red, Green, Blue bands
42 : * (in that order)
43 : * @param nOctaveStart Number of bottom octave. Octave numbers starts from one.
44 : * This value directly and strongly affects to amount of recognized points
45 : * @param nOctaveEnd Number of top octave. Should be equal or greater than
46 : * octaveStart
47 : * @param dfThreshold Value from 0 to 1. Threshold for feature point
48 : * recognition. Number of detected points is larger if threshold is lower
49 : *
50 : * @see GDALFeaturePoint, GDALSimpleSURF class for details.
51 : *
52 : * @note Every octave finds points in specific size. For small images
53 : * use small octave numbers, for high resolution - large.
54 : * For 1024x1024 images it's normal to use any octave numbers from range 1-6.
55 : * (for example, octave start - 1, octave end - 3, or octave start - 2, octave
56 : * end - 2.)
57 : * For larger images, try 1-10 range or even higher.
58 : * Pay attention that number of detected point decreases quickly per octave for
59 : * particular image. Algorithm finds more points in case of small octave number.
60 : * If method detects nothing, reduce octave start value.
61 : * In addition, if many feature points are required (the largest possible
62 : * amount), use the lowest octave start value (1) and wide octave range.
63 : *
64 : * @note Typical threshold's value is 0,001. It's pretty good for all images.
65 : * But this value depends on image's nature and may be various in each
66 : * particular case.
67 : * For example, value can be 0,002 or 0,005.
68 : * Notice that number of detected points is larger if threshold is lower.
69 : * But with high threshold feature points will be better - "stronger", more
70 : * "unique" and distinctive.
71 : *
72 : * Feel free to experiment with parameters, because character, robustness and
73 : * number of points entirely depend on provided range of octaves and threshold.
74 : *
75 : * NOTICE that every octave requires time to compute. Use a little range
76 : * or only one octave, if execution time is significant.
77 : *
78 : * @return CE_None or CE_Failure if error occurs.
79 : */
80 :
81 : static std::vector<GDALFeaturePoint> *
82 0 : GatherFeaturePoints(GDALDataset *poDataset, int *panBands, int nOctaveStart,
83 : int nOctaveEnd, double dfThreshold)
84 : {
85 0 : if (poDataset == nullptr)
86 : {
87 0 : CPLError(CE_Failure, CPLE_AppDefined, "GDALDataset isn't specified");
88 0 : return nullptr;
89 : }
90 :
91 0 : if (panBands == nullptr)
92 : {
93 0 : CPLError(CE_Failure, CPLE_AppDefined, "Raster bands are not specified");
94 0 : return nullptr;
95 : }
96 :
97 0 : if (nOctaveStart <= 0 || nOctaveEnd < 0 || nOctaveStart > nOctaveEnd)
98 : {
99 0 : CPLError(CE_Failure, CPLE_AppDefined, "Octave numbers are invalid");
100 0 : return nullptr;
101 : }
102 :
103 0 : if (dfThreshold < 0)
104 : {
105 0 : CPLError(CE_Failure, CPLE_AppDefined,
106 : "Threshold have to be greater than zero");
107 0 : return nullptr;
108 : }
109 :
110 0 : GDALRasterBand *poRstRedBand = poDataset->GetRasterBand(panBands[0]);
111 0 : GDALRasterBand *poRstGreenBand = poDataset->GetRasterBand(panBands[1]);
112 0 : GDALRasterBand *poRstBlueBand = poDataset->GetRasterBand(panBands[2]);
113 :
114 0 : const int nWidth = poRstRedBand->GetXSize();
115 0 : const int nHeight = poRstRedBand->GetYSize();
116 :
117 0 : if (nWidth == 0 || nHeight == 0)
118 : {
119 0 : CPLError(CE_Failure, CPLE_AppDefined,
120 : "Must have non-zero width and height.");
121 0 : return nullptr;
122 : }
123 :
124 : // Allocate memory for grayscale image.
125 0 : double **padfImg = new double *[nHeight];
126 0 : for (int i = 0;;)
127 : {
128 0 : padfImg[i] = new double[nWidth];
129 0 : for (int j = 0; j < nWidth; ++j)
130 0 : padfImg[i][j] = 0.0;
131 0 : ++i;
132 0 : if (i == nHeight)
133 0 : break;
134 0 : }
135 :
136 : // Create grayscale image.
137 0 : GDALSimpleSURF::ConvertRGBToLuminosity(poRstRedBand, poRstGreenBand,
138 : poRstBlueBand, nWidth, nHeight,
139 : padfImg, nHeight, nWidth);
140 :
141 : // Prepare integral image.
142 0 : GDALIntegralImage *poImg = new GDALIntegralImage();
143 0 : poImg->Initialize(const_cast<const double **>(padfImg), nHeight, nWidth);
144 :
145 : // Get feature points.
146 0 : GDALSimpleSURF *poSurf = new GDALSimpleSURF(nOctaveStart, nOctaveEnd);
147 :
148 : std::vector<GDALFeaturePoint> *poCollection =
149 0 : poSurf->ExtractFeaturePoints(poImg, dfThreshold);
150 :
151 : // Clean up.
152 0 : delete poImg;
153 0 : delete poSurf;
154 :
155 0 : for (int i = 0; i < nHeight; ++i)
156 0 : delete[] padfImg[i];
157 :
158 0 : delete[] padfImg;
159 :
160 0 : return poCollection;
161 : }
162 :
163 : /************************************************************************/
164 : /* GDALComputeMatchingPoints() */
165 : /************************************************************************/
166 :
167 : /** GDALComputeMatchingPoints. TODO document */
168 0 : GDAL_GCP CPL_DLL *GDALComputeMatchingPoints(GDALDatasetH hFirstImage,
169 : GDALDatasetH hSecondImage,
170 : char **papszOptions,
171 : int *pnGCPCount)
172 : {
173 0 : *pnGCPCount = 0;
174 :
175 : /* -------------------------------------------------------------------- */
176 : /* Override default algorithm parameters. */
177 : /* -------------------------------------------------------------------- */
178 : int nOctaveStart, nOctaveEnd;
179 : double dfSURFThreshold;
180 :
181 : nOctaveStart =
182 0 : atoi(CSLFetchNameValueDef(papszOptions, "OCTAVE_START", "2"));
183 0 : nOctaveEnd = atoi(CSLFetchNameValueDef(papszOptions, "OCTAVE_END", "2"));
184 :
185 : dfSURFThreshold =
186 0 : CPLAtof(CSLFetchNameValueDef(papszOptions, "SURF_THRESHOLD", "0.001"));
187 0 : const double dfMatchingThreshold = CPLAtof(
188 : CSLFetchNameValueDef(papszOptions, "MATCHING_THRESHOLD", "0.015"));
189 :
190 : /* -------------------------------------------------------------------- */
191 : /* Identify the bands to use. For now we are effectively */
192 : /* limited to using RGB input so if we have one band only treat */
193 : /* it as red=green=blue=band 1. Disallow non eightbit imagery. */
194 : /* -------------------------------------------------------------------- */
195 0 : int anBandMap1[3] = {1, 1, 1};
196 0 : if (GDALGetRasterCount(hFirstImage) >= 3)
197 : {
198 0 : anBandMap1[1] = 2;
199 0 : anBandMap1[2] = 3;
200 : }
201 :
202 0 : int anBandMap2[3] = {1, 1, 1};
203 0 : if (GDALGetRasterCount(hSecondImage) >= 3)
204 : {
205 0 : anBandMap2[1] = 2;
206 0 : anBandMap2[2] = 3;
207 : }
208 :
209 : /* -------------------------------------------------------------------- */
210 : /* Collect reference points on each image. */
211 : /* -------------------------------------------------------------------- */
212 : std::vector<GDALFeaturePoint> *poFPCollection1 =
213 0 : GatherFeaturePoints(GDALDataset::FromHandle(hFirstImage), anBandMap1,
214 : nOctaveStart, nOctaveEnd, dfSURFThreshold);
215 0 : if (poFPCollection1 == nullptr)
216 0 : return nullptr;
217 :
218 : std::vector<GDALFeaturePoint> *poFPCollection2 =
219 0 : GatherFeaturePoints(GDALDataset::FromHandle(hSecondImage), anBandMap2,
220 : nOctaveStart, nOctaveEnd, dfSURFThreshold);
221 :
222 0 : if (poFPCollection2 == nullptr)
223 : {
224 0 : delete poFPCollection1;
225 0 : return nullptr;
226 : }
227 :
228 : /* -------------------------------------------------------------------- */
229 : /* Try to find corresponding locations. */
230 : /* -------------------------------------------------------------------- */
231 0 : std::vector<GDALFeaturePoint *> oMatchPairs;
232 :
233 0 : if (CE_None != GDALSimpleSURF::MatchFeaturePoints(
234 : &oMatchPairs, poFPCollection1, poFPCollection2,
235 : dfMatchingThreshold))
236 : {
237 0 : delete poFPCollection1;
238 0 : delete poFPCollection2;
239 0 : return nullptr;
240 : }
241 :
242 0 : *pnGCPCount = static_cast<int>(oMatchPairs.size()) / 2;
243 :
244 : /* -------------------------------------------------------------------- */
245 : /* Translate these into GCPs - but with the output coordinate */
246 : /* system being pixel/line on the second image. */
247 : /* -------------------------------------------------------------------- */
248 : GDAL_GCP *pasGCPList =
249 0 : static_cast<GDAL_GCP *>(CPLCalloc(*pnGCPCount, sizeof(GDAL_GCP)));
250 :
251 0 : GDALInitGCPs(*pnGCPCount, pasGCPList);
252 :
253 0 : for (int i = 0; i < *pnGCPCount; i++)
254 : {
255 0 : GDALFeaturePoint *poPoint1 = oMatchPairs[i * 2];
256 0 : GDALFeaturePoint *poPoint2 = oMatchPairs[i * 2 + 1];
257 :
258 0 : pasGCPList[i].dfGCPPixel = poPoint1->GetX() + 0.5;
259 0 : pasGCPList[i].dfGCPLine = poPoint1->GetY() + 0.5;
260 :
261 0 : pasGCPList[i].dfGCPX = poPoint2->GetX() + 0.5;
262 0 : pasGCPList[i].dfGCPY = poPoint2->GetY() + 0.5;
263 0 : pasGCPList[i].dfGCPZ = 0.0;
264 : }
265 :
266 : // Cleanup the feature point lists.
267 0 : delete poFPCollection1;
268 0 : delete poFPCollection2;
269 :
270 : /* -------------------------------------------------------------------- */
271 : /* Optionally transform into the georef coordinates of the */
272 : /* output image. */
273 : /* -------------------------------------------------------------------- */
274 : const bool bGeorefOutput =
275 0 : CPLTestBool(CSLFetchNameValueDef(papszOptions, "OUTPUT_GEOREF", "NO"));
276 :
277 0 : if (bGeorefOutput)
278 : {
279 0 : double adfGeoTransform[6] = {};
280 :
281 0 : GDALGetGeoTransform(hSecondImage, adfGeoTransform);
282 :
283 0 : for (int i = 0; i < *pnGCPCount; i++)
284 : {
285 0 : GDALApplyGeoTransform(adfGeoTransform, pasGCPList[i].dfGCPX,
286 0 : pasGCPList[i].dfGCPY, &(pasGCPList[i].dfGCPX),
287 0 : &(pasGCPList[i].dfGCPY));
288 : }
289 : }
290 :
291 0 : return pasGCPList;
292 : }
|
