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