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