1 : /******************************************************************************
2 : * $Id: gdal_rpc.cpp 19870 2010-06-14 22:36:27Z rouault $
3 : *
4 : * Project: Image Warper
5 : * Purpose: Implements a rational polynomail (RPC) based transformer.
6 : * Author: Frank Warmerdam, warmerdam@pobox.com
7 : *
8 : ******************************************************************************
9 : * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
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_priv.h"
31 : #include "gdal_alg.h"
32 : #include "ogr_spatialref.h"
33 : #include "cpl_minixml.h"
34 :
35 : CPL_CVSID("$Id: gdal_rpc.cpp 19870 2010-06-14 22:36:27Z rouault $");
36 :
37 : CPL_C_START
38 : CPLXMLNode *GDALSerializeRPCTransformer( void *pTransformArg );
39 : void *GDALDeserializeRPCTransformer( CPLXMLNode *psTree );
40 : CPL_C_END
41 :
42 : /************************************************************************/
43 : /* RPCInfoToMD() */
44 : /* */
45 : /* Turn an RPCInfo structure back into it's metadata format. */
46 : /************************************************************************/
47 :
48 0 : static char ** RPCInfoToMD( GDALRPCInfo *psRPCInfo )
49 :
50 : {
51 0 : char **papszMD = NULL;
52 0 : CPLString osField, osMultiField;
53 : int i;
54 :
55 0 : osField.Printf( "%.15g", psRPCInfo->dfLINE_OFF );
56 0 : papszMD = CSLSetNameValue( papszMD, "LINE_OFF", osField );
57 :
58 0 : osField.Printf( "%.15g", psRPCInfo->dfSAMP_OFF );
59 0 : papszMD = CSLSetNameValue( papszMD, "SAMP_OFF", osField );
60 :
61 0 : osField.Printf( "%.15g", psRPCInfo->dfLAT_OFF );
62 0 : papszMD = CSLSetNameValue( papszMD, "LAT_OFF", osField );
63 :
64 0 : osField.Printf( "%.15g", psRPCInfo->dfLONG_OFF );
65 0 : papszMD = CSLSetNameValue( papszMD, "LONG_OFF", osField );
66 :
67 0 : osField.Printf( "%.15g", psRPCInfo->dfHEIGHT_OFF );
68 0 : papszMD = CSLSetNameValue( papszMD, "HEIGHT_OFF", osField );
69 :
70 0 : osField.Printf( "%.15g", psRPCInfo->dfLINE_SCALE );
71 0 : papszMD = CSLSetNameValue( papszMD, "LINE_SCALE", osField );
72 :
73 0 : osField.Printf( "%.15g", psRPCInfo->dfSAMP_SCALE );
74 0 : papszMD = CSLSetNameValue( papszMD, "SAMP_SCALE", osField );
75 :
76 0 : osField.Printf( "%.15g", psRPCInfo->dfLAT_SCALE );
77 0 : papszMD = CSLSetNameValue( papszMD, "LAT_SCALE", osField );
78 :
79 0 : osField.Printf( "%.15g", psRPCInfo->dfLONG_SCALE );
80 0 : papszMD = CSLSetNameValue( papszMD, "LONG_SCALE", osField );
81 :
82 0 : osField.Printf( "%.15g", psRPCInfo->dfHEIGHT_SCALE );
83 0 : papszMD = CSLSetNameValue( papszMD, "HEIGHT_SCALE", osField );
84 :
85 0 : osField.Printf( "%.15g", psRPCInfo->dfMIN_LONG );
86 0 : papszMD = CSLSetNameValue( papszMD, "MIN_LONG", osField );
87 :
88 0 : osField.Printf( "%.15g", psRPCInfo->dfMIN_LAT );
89 0 : papszMD = CSLSetNameValue( papszMD, "MIN_LAT", osField );
90 :
91 0 : osField.Printf( "%.15g", psRPCInfo->dfMAX_LONG );
92 0 : papszMD = CSLSetNameValue( papszMD, "MAX_LONG", osField );
93 :
94 0 : osField.Printf( "%.15g", psRPCInfo->dfMAX_LAT );
95 0 : papszMD = CSLSetNameValue( papszMD, "MAX_LAT", osField );
96 :
97 0 : for( i = 0; i < 20; i++ )
98 : {
99 0 : osField.Printf( "%.15g", psRPCInfo->adfLINE_NUM_COEFF[i] );
100 0 : if( i > 0 )
101 0 : osMultiField += " ";
102 : else
103 0 : osMultiField = "";
104 0 : osMultiField += osField;
105 : }
106 0 : papszMD = CSLSetNameValue( papszMD, "LINE_NUM_COEFF", osMultiField );
107 :
108 0 : for( i = 0; i < 20; i++ )
109 : {
110 0 : osField.Printf( "%.15g", psRPCInfo->adfLINE_DEN_COEFF[i] );
111 0 : if( i > 0 )
112 0 : osMultiField += " ";
113 : else
114 0 : osMultiField = "";
115 0 : osMultiField += osField;
116 : }
117 0 : papszMD = CSLSetNameValue( papszMD, "LINE_DEN_COEFF", osMultiField );
118 :
119 0 : for( i = 0; i < 20; i++ )
120 : {
121 0 : osField.Printf( "%.15g", psRPCInfo->adfSAMP_NUM_COEFF[i] );
122 0 : if( i > 0 )
123 0 : osMultiField += " ";
124 : else
125 0 : osMultiField = "";
126 0 : osMultiField += osField;
127 : }
128 0 : papszMD = CSLSetNameValue( papszMD, "SAMP_NUM_COEFF", osMultiField );
129 :
130 0 : for( i = 0; i < 20; i++ )
131 : {
132 0 : osField.Printf( "%.15g", psRPCInfo->adfSAMP_DEN_COEFF[i] );
133 0 : if( i > 0 )
134 0 : osMultiField += " ";
135 : else
136 0 : osMultiField = "";
137 0 : osMultiField += osField;
138 : }
139 0 : papszMD = CSLSetNameValue( papszMD, "SAMP_DEN_COEFF", osMultiField );
140 :
141 0 : return papszMD;
142 : }
143 :
144 : /************************************************************************/
145 : /* RPCComputeTerms() */
146 : /************************************************************************/
147 :
148 : static void RPCComputeTerms( double dfLong, double dfLat, double dfHeight,
149 28 : double *padfTerms )
150 :
151 : {
152 28 : padfTerms[0] = 1.0;
153 28 : padfTerms[1] = dfLong;
154 28 : padfTerms[2] = dfLat;
155 28 : padfTerms[3] = dfHeight;
156 28 : padfTerms[4] = dfLong * dfLat;
157 28 : padfTerms[5] = dfLong * dfHeight;
158 28 : padfTerms[6] = dfLat * dfHeight;
159 28 : padfTerms[7] = dfLong * dfLong;
160 28 : padfTerms[8] = dfLat * dfLat;
161 28 : padfTerms[9] = dfHeight * dfHeight;
162 :
163 28 : padfTerms[10] = dfLong * dfLat * dfHeight;
164 28 : padfTerms[11] = dfLong * dfLong * dfLong;
165 28 : padfTerms[12] = dfLong * dfLat * dfLat;
166 28 : padfTerms[13] = dfLong * dfHeight * dfHeight;
167 28 : padfTerms[14] = dfLong * dfLong * dfLat;
168 28 : padfTerms[15] = dfLat * dfLat * dfLat;
169 28 : padfTerms[16] = dfLat * dfHeight * dfHeight;
170 28 : padfTerms[17] = dfLong * dfLong * dfHeight;
171 28 : padfTerms[18] = dfLat * dfLat * dfHeight;
172 28 : padfTerms[19] = dfHeight * dfHeight * dfHeight;
173 28 : }
174 :
175 : /************************************************************************/
176 : /* RPCEvaluate() */
177 : /************************************************************************/
178 :
179 112 : static double RPCEvaluate( double *padfTerms, double *padfCoefs )
180 :
181 : {
182 112 : double dfSum = 0.0;
183 : int i;
184 :
185 2352 : for( i = 0; i < 20; i++ )
186 2240 : dfSum += padfTerms[i] * padfCoefs[i];
187 :
188 112 : return dfSum;
189 : }
190 :
191 : /************************************************************************/
192 : /* RPCTransformPoint() */
193 : /************************************************************************/
194 :
195 : static void RPCTransformPoint( GDALRPCInfo *psRPC,
196 : double dfLong, double dfLat, double dfHeight,
197 28 : double *pdfPixel, double *pdfLine )
198 :
199 : {
200 : double dfResultX, dfResultY;
201 : double adfTerms[20];
202 :
203 : RPCComputeTerms(
204 : (dfLong - psRPC->dfLONG_OFF) / psRPC->dfLONG_SCALE,
205 : (dfLat - psRPC->dfLAT_OFF) / psRPC->dfLAT_SCALE,
206 : (dfHeight - psRPC->dfHEIGHT_OFF) / psRPC->dfHEIGHT_SCALE,
207 28 : adfTerms );
208 :
209 : dfResultX = RPCEvaluate( adfTerms, psRPC->adfSAMP_NUM_COEFF )
210 28 : / RPCEvaluate( adfTerms, psRPC->adfSAMP_DEN_COEFF );
211 :
212 : dfResultY = RPCEvaluate( adfTerms, psRPC->adfLINE_NUM_COEFF )
213 28 : / RPCEvaluate( adfTerms, psRPC->adfLINE_DEN_COEFF );
214 :
215 28 : *pdfPixel = dfResultX * psRPC->dfSAMP_SCALE + psRPC->dfSAMP_OFF;
216 28 : *pdfLine = dfResultY * psRPC->dfLINE_SCALE + psRPC->dfLINE_OFF;
217 28 : }
218 :
219 : /************************************************************************/
220 : /* ==================================================================== */
221 : /* GDALRPCTransformer */
222 : /* ==================================================================== */
223 : /************************************************************************/
224 :
225 : typedef struct {
226 :
227 : GDALTransformerInfo sTI;
228 :
229 : GDALRPCInfo sRPC;
230 :
231 : double adfPLToLatLongGeoTransform[6];
232 :
233 : int bReversed;
234 :
235 : double dfPixErrThreshold;
236 :
237 : double dfHeightOffset;
238 :
239 : double dfHeightScale;
240 :
241 : char *pszDEMPath;
242 :
243 : int bHasTriedOpeningDS;
244 : GDALDataset *poDS;
245 :
246 : OGRCoordinateTransformation *poCT;
247 :
248 : double adfGeoTransform[6];
249 : double adfReverseGeoTransform[6];
250 : } GDALRPCTransformInfo;
251 :
252 : /************************************************************************/
253 : /* GDALCreateRPCTransformer() */
254 : /************************************************************************/
255 :
256 : /**
257 : * Create an RPC based transformer.
258 : *
259 : * The geometric sensor model describing the physical relationship between
260 : * image coordinates and ground coordinate is known as a Rigorous Projection
261 : * Model. A Rigorous Projection Model expresses the mapping of the image space
262 : * coordinates of rows and columns (r,c) onto the object space reference
263 : * surface geodetic coordinates (long, lat, height).
264 : *
265 : * RPC supports a generic description of the Rigorous Projection Models. The
266 : * approximation used by GDAL (RPC00) is a set of rational polynomials exp
267 : * ressing the normalized row and column values, (rn , cn), as a function of
268 : * normalized geodetic latitude, longitude, and height, (P, L, H), given a
269 : * set of normalized polynomial coefficients (LINE_NUM_COEF_n, LINE_DEN_COEF_n,
270 : * SAMP_NUM_COEF_n, SAMP_DEN_COEF_n). Normalized values, rather than actual
271 : * values are used in order to minimize introduction of errors during the
272 : * calculations. The transformation between row and column values (r,c), and
273 : * normalized row and column values (rn, cn), and between the geodetic
274 : * latitude, longitude, and height and normalized geodetic latitude,
275 : * longitude, and height (P, L, H), is defined by a set of normalizing
276 : * translations (offsets) and scales that ensure all values are contained i
277 : * the range -1 to +1.
278 : *
279 : * This function creates a GDALTransformFunc compatible transformer
280 : * for going between image pixel/line and long/lat/height coordinates
281 : * using RPCs. The RPCs are provided in a GDALRPCInfo structure which is
282 : * normally read from metadata using GDALExtractRPCInfo().
283 : *
284 : * GDAL RPC Metadata has the following entries (also described in GDAL RFC 22
285 : * and the GeoTIFF RPC document http://geotiff.maptools.org/rpc_prop.html.
286 : *
287 : * <ul>
288 : * <li>ERR_BIAS: Error - Bias. The RMS bias error in meters per horizontal axis of all points in the image (-1.0 if unknown)
289 : * <li>ERR_RAND: Error - Random. RMS random error in meters per horizontal axis of each point in the image (-1.0 if unknown)
290 : * <li>LINE_OFF: Line Offset
291 : * <li>SAMP_OFF: Sample Offset
292 : * <li>LAT_OFF: Geodetic Latitude Offset
293 : * <li>LONG_OFF: Geodetic Longitude Offset
294 : * <li>HEIGHT_OFF: Geodetic Height Offset
295 : * <li>LINE_SCALE: Line Scale
296 : * <li>SAMP_SCALE: Sample Scale
297 : * <li>LAT_SCALE: Geodetic Latitude Scale
298 : * <li>LONG_SCALE: Geodetic Longitude Scale
299 : * <li>HEIGHT_SCALE: Geodetic Height Scale
300 : * <li>LINE_NUM_COEFF (1-20): Line Numerator Coefficients. Twenty coefficients for the polynomial in the Numerator of the rn equation. (space separated)
301 : * <li>LINE_DEN_COEFF (1-20): Line Denominator Coefficients. Twenty coefficients for the polynomial in the Denominator of the rn equation. (space separated)
302 : * <li>SAMP_NUM_COEFF (1-20): Sample Numerator Coefficients. Twenty coefficients for the polynomial in the Numerator of the cn equation. (space separated)
303 : * <li>SAMP_DEN_COEFF (1-20): Sample Denominator Coefficients. Twenty coefficients for the polynomial in the Denominator of the cn equation. (space separated)
304 : * </ul>
305 : *
306 : * The transformer normally maps from pixel/line/height to long/lat/height space
307 : * as a forward transformation though in RPC terms that would be considered
308 : * an inverse transformation (and is solved by iterative approximation using
309 : * long/lat/height to pixel/line transformations). The default direction can
310 : * be reversed by passing bReversed=TRUE.
311 : *
312 : * The iterative solution of pixel/line
313 : * to lat/long/height is currently run for up to 10 iterations or until
314 : * the apparent error is less than dfPixErrThreshold pixels. Passing zero
315 : * will not avoid all error, but will cause the operation to run for the maximum
316 : * number of iterations.
317 : *
318 : * Additional options to the transformer can be supplied in papszOptions.
319 : *
320 : * Options:
321 : *
322 : * <ul>
323 : * <li> RPC_HEIGHT: a fixed height offset to be applied to all points passed
324 : * in. In this situation the Z passed into the transformation function is
325 : * assumed to be height above ground, and the RPC_HEIGHT is assumed to be
326 : * an average height above sea level for ground in the target scene.
327 : *
328 : * <li> RPC_HEIGHT_SCALE: a factor used to multiply heights above ground.
329 : * Usefull when elevation offsets of the DEM are not expressed in meters. (GDAL >= 1.8.0)
330 : *
331 : * <li> RPC_DEM: the name of a GDAL dataset (a DEM file typically) used to
332 : * extract elevation offsets from. In this situation the Z passed into the
333 : * transformation function is assumed to be height above ground. This option
334 : * should be used in replacement of RPC_HEIGHT to provide a way of defining
335 : * a non uniform ground for the target scene (GDAL >= 1.8.0)
336 : * </ul>
337 : *
338 : * @param psRPCInfo Definition of the RPC parameters.
339 : *
340 : * @param bReversed If true "forward" transformation will be lat/long to pixel/line instead of the normal pixel/line to lat/long.
341 : *
342 : * @param dfPixErrThreshold the error (measured in pixels) allowed in the
343 : * iterative solution of pixel/line to lat/long computations (the other way
344 : * is always exact given the equations).
345 : *
346 : * @param papszOptions Other transformer options (ie. RPC_HEIGHT=<z>).
347 : *
348 : * @return transformer callback data (deallocate with GDALDestroyTransformer()).
349 : */
350 :
351 : void *GDALCreateRPCTransformer( GDALRPCInfo *psRPCInfo, int bReversed,
352 : double dfPixErrThreshold,
353 3 : char **papszOptions )
354 :
355 : {
356 : GDALRPCTransformInfo *psTransform;
357 :
358 : /* -------------------------------------------------------------------- */
359 : /* Initialize core info. */
360 : /* -------------------------------------------------------------------- */
361 : psTransform = (GDALRPCTransformInfo *)
362 3 : CPLCalloc(sizeof(GDALRPCTransformInfo),1);
363 :
364 3 : memcpy( &(psTransform->sRPC), psRPCInfo, sizeof(GDALRPCInfo) );
365 3 : psTransform->bReversed = bReversed;
366 3 : psTransform->dfPixErrThreshold = dfPixErrThreshold;
367 3 : psTransform->dfHeightOffset = 0.0;
368 3 : psTransform->dfHeightScale = 1.0;
369 :
370 3 : strcpy( psTransform->sTI.szSignature, "GTI" );
371 3 : psTransform->sTI.pszClassName = "GDALRPCTransformer";
372 3 : psTransform->sTI.pfnTransform = GDALRPCTransform;
373 3 : psTransform->sTI.pfnCleanup = GDALDestroyRPCTransformer;
374 3 : psTransform->sTI.pfnSerialize = GDALSerializeRPCTransformer;
375 :
376 : /* -------------------------------------------------------------------- */
377 : /* Do we have a "average height" that we want to consider all */
378 : /* elevations to be relative to? */
379 : /* -------------------------------------------------------------------- */
380 3 : const char *pszHeight = CSLFetchNameValue( papszOptions, "RPC_HEIGHT" );
381 3 : if( pszHeight != NULL )
382 1 : psTransform->dfHeightOffset = CPLAtof(pszHeight);
383 :
384 : /* -------------------------------------------------------------------- */
385 : /* The "height scale" */
386 : /* -------------------------------------------------------------------- */
387 3 : const char *pszHeightScale = CSLFetchNameValue( papszOptions, "RPC_HEIGHT_SCALE" );
388 3 : if( pszHeightScale != NULL )
389 1 : psTransform->dfHeightScale = CPLAtof(pszHeightScale);
390 :
391 : /* -------------------------------------------------------------------- */
392 : /* The DEM file name */
393 : /* -------------------------------------------------------------------- */
394 3 : const char *pszDEMPath = CSLFetchNameValue( papszOptions, "RPC_DEM" );
395 3 : if( pszDEMPath != NULL )
396 1 : psTransform->pszDEMPath = CPLStrdup(pszDEMPath);
397 :
398 : /* -------------------------------------------------------------------- */
399 : /* Establish a reference point for calcualating an affine */
400 : /* geotransform approximate transformation. */
401 : /* -------------------------------------------------------------------- */
402 3 : double adfGTFromLL[6], dfRefPixel = -1.0, dfRefLine = -1.0;
403 3 : double dfRefLong = 0.0, dfRefLat = 0.0;
404 :
405 3 : if( psRPCInfo->dfMIN_LONG != -180 || psRPCInfo->dfMAX_LONG != 180 )
406 : {
407 0 : dfRefLong = (psRPCInfo->dfMIN_LONG + psRPCInfo->dfMAX_LONG) * 0.5;
408 0 : dfRefLat = (psRPCInfo->dfMIN_LAT + psRPCInfo->dfMAX_LAT ) * 0.5;
409 :
410 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat, 0.0,
411 0 : &dfRefPixel, &dfRefLine );
412 : }
413 :
414 : // Try with scale and offset if we don't can't use bounds or
415 : // the results seem daft.
416 3 : if( dfRefPixel < 0.0 || dfRefLine < 0.0
417 : || dfRefPixel > 100000 || dfRefLine > 100000 )
418 : {
419 3 : dfRefLong = psRPCInfo->dfLONG_OFF;
420 3 : dfRefLat = psRPCInfo->dfLAT_OFF;
421 :
422 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat, 0.0,
423 3 : &dfRefPixel, &dfRefLine );
424 : }
425 :
426 : /* -------------------------------------------------------------------- */
427 : /* Transform nearby locations to establish affine direction */
428 : /* vectors. */
429 : /* -------------------------------------------------------------------- */
430 3 : double dfRefPixelDelta, dfRefLineDelta, dfLLDelta = 0.0001;
431 :
432 : RPCTransformPoint( psRPCInfo, dfRefLong+dfLLDelta, dfRefLat, 0.0,
433 3 : &dfRefPixelDelta, &dfRefLineDelta );
434 3 : adfGTFromLL[1] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
435 3 : adfGTFromLL[2] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
436 :
437 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat+dfLLDelta, 0.0,
438 3 : &dfRefPixelDelta, &dfRefLineDelta );
439 3 : adfGTFromLL[4] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
440 3 : adfGTFromLL[5] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
441 :
442 : adfGTFromLL[0] = dfRefPixel
443 3 : - adfGTFromLL[1] * dfRefLong - adfGTFromLL[2] * dfRefLat;
444 : adfGTFromLL[3] = dfRefLine
445 3 : - adfGTFromLL[4] * dfRefLong - adfGTFromLL[5] * dfRefLat;
446 :
447 3 : GDALInvGeoTransform( adfGTFromLL, psTransform->adfPLToLatLongGeoTransform);
448 :
449 3 : return psTransform;
450 : }
451 :
452 : /************************************************************************/
453 : /* GDALDestroyReprojectionTransformer() */
454 : /************************************************************************/
455 :
456 3 : void GDALDestroyRPCTransformer( void *pTransformAlg )
457 :
458 : {
459 3 : GDALRPCTransformInfo *psTransform = (GDALRPCTransformInfo *) pTransformAlg;
460 :
461 3 : CPLFree( psTransform->pszDEMPath );
462 :
463 3 : if(psTransform->poDS)
464 1 : GDALClose(psTransform->poDS);
465 3 : if(psTransform->poCT)
466 1 : OCTDestroyCoordinateTransformation((OGRCoordinateTransformationH)psTransform->poCT);
467 :
468 3 : CPLFree( pTransformAlg );
469 3 : }
470 :
471 : /************************************************************************/
472 : /* RPCInverseTransformPoint() */
473 : /************************************************************************/
474 :
475 : static void
476 : RPCInverseTransformPoint( GDALRPCTransformInfo *psTransform,
477 : double dfPixel, double dfLine, double dfHeight,
478 5 : double *pdfLong, double *pdfLat )
479 :
480 : {
481 : double dfResultX, dfResultY;
482 : int iIter;
483 5 : GDALRPCInfo *psRPC = &(psTransform->sRPC);
484 :
485 : /* -------------------------------------------------------------------- */
486 : /* Compute an initial approximation based on linear */
487 : /* interpolation from our reference point. */
488 : /* -------------------------------------------------------------------- */
489 : dfResultX = psTransform->adfPLToLatLongGeoTransform[0]
490 : + psTransform->adfPLToLatLongGeoTransform[1] * dfPixel
491 5 : + psTransform->adfPLToLatLongGeoTransform[2] * dfLine;
492 :
493 : dfResultY = psTransform->adfPLToLatLongGeoTransform[3]
494 : + psTransform->adfPLToLatLongGeoTransform[4] * dfPixel
495 5 : + psTransform->adfPLToLatLongGeoTransform[5] * dfLine;
496 :
497 : /* -------------------------------------------------------------------- */
498 : /* Now iterate, trying to find a closer LL location that will */
499 : /* back transform to the indicated pixel and line. */
500 : /* -------------------------------------------------------------------- */
501 : double dfPixelDeltaX, dfPixelDeltaY;
502 :
503 15 : for( iIter = 0; iIter < 10; iIter++ )
504 : {
505 : double dfBackPixel, dfBackLine;
506 :
507 : RPCTransformPoint( psRPC, dfResultX, dfResultY, dfHeight,
508 15 : &dfBackPixel, &dfBackLine );
509 :
510 15 : dfPixelDeltaX = dfBackPixel - dfPixel;
511 15 : dfPixelDeltaY = dfBackLine - dfLine;
512 :
513 : dfResultX = dfResultX
514 : - dfPixelDeltaX * psTransform->adfPLToLatLongGeoTransform[1]
515 15 : - dfPixelDeltaY * psTransform->adfPLToLatLongGeoTransform[2];
516 : dfResultY = dfResultY
517 : - dfPixelDeltaX * psTransform->adfPLToLatLongGeoTransform[4]
518 15 : - dfPixelDeltaY * psTransform->adfPLToLatLongGeoTransform[5];
519 :
520 15 : if( ABS(dfPixelDeltaX) < psTransform->dfPixErrThreshold
521 : && ABS(dfPixelDeltaY) < psTransform->dfPixErrThreshold )
522 : {
523 5 : iIter = -1;
524 : //CPLDebug( "RPC", "Converged!" );
525 5 : break;
526 : }
527 :
528 : }
529 :
530 5 : if( iIter != -1 )
531 : CPLDebug( "RPC", "Iterations %d: Got: %g,%g Offset=%g,%g",
532 : iIter,
533 : dfResultX, dfResultY,
534 0 : dfPixelDeltaX, dfPixelDeltaY );
535 :
536 5 : *pdfLong = dfResultX;
537 5 : *pdfLat = dfResultY;
538 5 : }
539 :
540 : /************************************************************************/
541 : /* GDALRPCTransform() */
542 : /************************************************************************/
543 :
544 : int GDALRPCTransform( void *pTransformArg, int bDstToSrc,
545 : int nPointCount,
546 : double *padfX, double *padfY, double *padfZ,
547 8 : int *panSuccess )
548 :
549 : {
550 8 : VALIDATE_POINTER1( pTransformArg, "GDALRPCTransform", 0 );
551 :
552 8 : GDALRPCTransformInfo *psTransform = (GDALRPCTransformInfo *) pTransformArg;
553 8 : GDALRPCInfo *psRPC = &(psTransform->sRPC);
554 : int i;
555 :
556 8 : if( psTransform->bReversed )
557 0 : bDstToSrc = !bDstToSrc;
558 :
559 8 : int bands[1] = {1};
560 8 : int nRasterXSize = 0, nRasterYSize = 0;
561 :
562 : /* -------------------------------------------------------------------- */
563 : /* Lazy opening of the optionnal DEM file. */
564 : /* -------------------------------------------------------------------- */
565 8 : if(psTransform->pszDEMPath != NULL &&
566 : psTransform->bHasTriedOpeningDS == FALSE)
567 : {
568 1 : int bIsValid = FALSE;
569 1 : psTransform->bHasTriedOpeningDS = TRUE;
570 : psTransform->poDS = (GDALDataset *)
571 1 : GDALOpen( psTransform->pszDEMPath, GA_ReadOnly );
572 1 : if(psTransform->poDS != NULL && psTransform->poDS->GetRasterCount() >= 1)
573 : {
574 1 : const char* pszSpatialRef = psTransform->poDS->GetProjectionRef();
575 1 : if (pszSpatialRef != NULL && pszSpatialRef[0] != '\0')
576 : {
577 : OGRSpatialReference* poWGSSpaRef =
578 1 : new OGRSpatialReference(SRS_WKT_WGS84);
579 : OGRSpatialReference* poDSSpaRef =
580 2 : new OGRSpatialReference(pszSpatialRef);
581 1 : if(!poWGSSpaRef->IsSame(poDSSpaRef))
582 : psTransform->poCT =OGRCreateCoordinateTransformation(
583 1 : poWGSSpaRef, poDSSpaRef );
584 1 : delete poWGSSpaRef;
585 1 : delete poDSSpaRef;
586 : }
587 :
588 1 : if (psTransform->poDS->GetGeoTransform(
589 : psTransform->adfGeoTransform) == CE_None &&
590 : GDALInvGeoTransform( psTransform->adfGeoTransform,
591 : psTransform->adfReverseGeoTransform ))
592 : {
593 1 : bIsValid = TRUE;
594 : }
595 : }
596 :
597 1 : if (!bIsValid && psTransform->poDS != NULL)
598 : {
599 0 : GDALClose(psTransform->poDS);
600 0 : psTransform->poDS = NULL;
601 : }
602 : }
603 8 : if (psTransform->poDS)
604 : {
605 2 : nRasterXSize = psTransform->poDS->GetRasterXSize();
606 2 : nRasterYSize = psTransform->poDS->GetRasterYSize();
607 : }
608 :
609 : /* -------------------------------------------------------------------- */
610 : /* The simple case is transforming from lat/long to pixel/line. */
611 : /* Just apply the equations directly. */
612 : /* -------------------------------------------------------------------- */
613 8 : if( bDstToSrc )
614 : {
615 8 : for( i = 0; i < nPointCount; i++ )
616 : {
617 4 : if(psTransform->poDS)
618 : {
619 : double dfX, dfY;
620 : //check if dem is not in WGS84 and transform points padfX[i], padfY[i]
621 1 : if(psTransform->poCT)
622 : {
623 1 : double dfXOrig = padfX[i];
624 1 : double dfYOrig = padfY[i];
625 1 : double dfZOrig = padfZ[i];
626 1 : if (!psTransform->poCT->Transform(
627 : 1, &dfXOrig, &dfYOrig, &dfZOrig))
628 : {
629 0 : panSuccess[i] = FALSE;
630 0 : continue;
631 : }
632 : GDALApplyGeoTransform( psTransform->adfReverseGeoTransform,
633 1 : dfXOrig, dfYOrig, &dfX, &dfY );
634 : }
635 : else
636 : GDALApplyGeoTransform( psTransform->adfReverseGeoTransform,
637 0 : padfX[i], padfY[i], &dfX, &dfY );
638 1 : int dX = int(dfX);
639 1 : int dY = int(dfY);
640 :
641 1 : if (!(dX >= 0 && dY >= 0 &&
642 : dX+2 <= nRasterXSize && dY+2 <= nRasterYSize))
643 : {
644 0 : panSuccess[i] = FALSE;
645 0 : continue;
646 : }
647 1 : int adElevData[4] = {0,0,0,0};
648 : CPLErr eErr = psTransform->poDS->RasterIO(GF_Read, dX, dY, 2, 2,
649 : &adElevData, 2, 2,
650 1 : GDT_Int32, 1, bands, 0, 0, 0);
651 1 : if(eErr != CE_None)
652 : {
653 0 : panSuccess[i] = FALSE;
654 0 : continue;
655 : }
656 : //bilinear interpolation
657 1 : double dfDeltaX = dfX - dX;
658 1 : double dfDeltaX1 = 1.0 - dfDeltaX;
659 1 : double dfDeltaY = dfY - dY;
660 1 : double dfDeltaY1 = 1.0 - dfDeltaY;
661 :
662 1 : double dfXZ1 = adElevData[0] * dfDeltaX1 + adElevData[1] * dfDeltaX;
663 1 : double dfXZ2 = adElevData[2] * dfDeltaX1 + adElevData[3] * dfDeltaX;
664 1 : double dfYZ = dfXZ1 * dfDeltaY1 + dfXZ2 * dfDeltaY;
665 :
666 : RPCTransformPoint( psRPC, padfX[i], padfY[i],
667 : padfZ[i] + (psTransform->dfHeightOffset + dfYZ) *
668 : psTransform->dfHeightScale,
669 1 : padfX + i, padfY + i );
670 : }
671 : else
672 : RPCTransformPoint( psRPC, padfX[i], padfY[i],
673 : padfZ[i] + psTransform->dfHeightOffset *
674 : psTransform->dfHeightScale,
675 3 : padfX + i, padfY + i );
676 4 : panSuccess[i] = TRUE;
677 : }
678 :
679 4 : return TRUE;
680 : }
681 :
682 : /* -------------------------------------------------------------------- */
683 : /* Compute the inverse (pixel/line/height to lat/long). This */
684 : /* function uses an iterative method from an initial linear */
685 : /* approximation. */
686 : /* -------------------------------------------------------------------- */
687 8 : for( i = 0; i < nPointCount; i++ )
688 : {
689 : double dfResultX, dfResultY;
690 :
691 4 : if(psTransform->poDS)
692 : {
693 : RPCInverseTransformPoint( psTransform, padfX[i], padfY[i],
694 : padfZ[i] + psTransform->dfHeightOffset *
695 : psTransform->dfHeightScale,
696 1 : &dfResultX, &dfResultY );
697 :
698 : double dfX, dfY;
699 : //check if dem is not in WGS84 and transform points padfX[i], padfY[i]
700 1 : if(psTransform->poCT)
701 : {
702 1 : double dfZ = 0;
703 1 : if (!psTransform->poCT->Transform(1, &dfResultX, &dfResultY, &dfZ))
704 : {
705 0 : panSuccess[i] = FALSE;
706 0 : continue;
707 : }
708 : }
709 :
710 : GDALApplyGeoTransform( psTransform->adfReverseGeoTransform,
711 1 : dfResultX, dfResultY, &dfX, &dfY );
712 1 : int dX = int(dfX);
713 1 : int dY = int(dfY);
714 :
715 1 : if (!(dX >= 0 && dY >= 0 && dX+2 <= nRasterXSize && dY+2 <= nRasterYSize))
716 : {
717 0 : panSuccess[i] = FALSE;
718 0 : continue;
719 : }
720 1 : int adElevData[4] = {0,0,0,0};
721 : CPLErr eErr = psTransform->poDS->RasterIO(GF_Read, dX, dY, 2, 2,
722 : &adElevData, 2, 2,
723 1 : GDT_Int32, 1, bands, 0, 0, 0);
724 1 : if(eErr != CE_None)
725 : {
726 0 : panSuccess[i] = FALSE;
727 0 : continue;
728 : }
729 : //bilinear interpolation
730 1 : double dfDeltaX = dfX - dX;
731 1 : double dfDeltaX1 = 1.0 - dfDeltaX;
732 1 : double dfDeltaY = dfY - dY;
733 1 : double dfDeltaY1 = 1.0 - dfDeltaY;
734 :
735 1 : double dfXZ1 = adElevData[0] * dfDeltaX1 + adElevData[1] * dfDeltaX;
736 1 : double dfXZ2 = adElevData[2] * dfDeltaX1 + adElevData[3] * dfDeltaX;
737 1 : double dfYZ = dfXZ1 * dfDeltaY1 + dfXZ2 * dfDeltaY;
738 :
739 : RPCInverseTransformPoint( psTransform, padfX[i], padfY[i],
740 : padfZ[i] + (psTransform->dfHeightOffset + dfYZ) *
741 : psTransform->dfHeightScale,
742 1 : &dfResultX, &dfResultY );
743 : }
744 : else
745 : {
746 : RPCInverseTransformPoint( psTransform, padfX[i], padfY[i],
747 : padfZ[i] + psTransform->dfHeightOffset *
748 : psTransform->dfHeightScale,
749 3 : &dfResultX, &dfResultY );
750 :
751 : }
752 4 : padfX[i] = dfResultX;
753 4 : padfY[i] = dfResultY;
754 :
755 4 : panSuccess[i] = TRUE;
756 : }
757 :
758 4 : return TRUE;
759 : }
760 :
761 : /************************************************************************/
762 : /* GDALSerializeRPCTransformer() */
763 : /************************************************************************/
764 :
765 0 : CPLXMLNode *GDALSerializeRPCTransformer( void *pTransformArg )
766 :
767 : {
768 0 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeRPCTransformer", NULL );
769 :
770 : CPLXMLNode *psTree;
771 : GDALRPCTransformInfo *psInfo =
772 0 : (GDALRPCTransformInfo *)(pTransformArg);
773 :
774 0 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "RPCTransformer" );
775 :
776 : /* -------------------------------------------------------------------- */
777 : /* Serialize bReversed. */
778 : /* -------------------------------------------------------------------- */
779 : CPLCreateXMLElementAndValue(
780 : psTree, "Reversed",
781 0 : CPLString().Printf( "%d", psInfo->bReversed ) );
782 :
783 : /* -------------------------------------------------------------------- */
784 : /* Serialize Height Offset. */
785 : /* -------------------------------------------------------------------- */
786 : CPLCreateXMLElementAndValue(
787 : psTree, "HeightOffset",
788 0 : CPLString().Printf( "%.15g", psInfo->dfHeightOffset ) );
789 :
790 : /* -------------------------------------------------------------------- */
791 : /* Serialize Height Scale. */
792 : /* -------------------------------------------------------------------- */
793 0 : if (psInfo->dfHeightScale != 1.0)
794 : CPLCreateXMLElementAndValue(
795 : psTree, "HeightScale",
796 0 : CPLString().Printf( "%.15g", psInfo->dfHeightScale ) );
797 :
798 : /* -------------------------------------------------------------------- */
799 : /* Serialize DEM path. */
800 : /* -------------------------------------------------------------------- */
801 0 : if (psInfo->pszDEMPath != NULL)
802 : CPLCreateXMLElementAndValue(
803 : psTree, "DEMPath",
804 0 : CPLString().Printf( "%s", psInfo->pszDEMPath ) );
805 :
806 : /* -------------------------------------------------------------------- */
807 : /* Serialize pixel error threshold. */
808 : /* -------------------------------------------------------------------- */
809 : CPLCreateXMLElementAndValue(
810 : psTree, "PixErrThreshold",
811 0 : CPLString().Printf( "%.15g", psInfo->dfPixErrThreshold ) );
812 :
813 : /* -------------------------------------------------------------------- */
814 : /* RPC metadata. */
815 : /* -------------------------------------------------------------------- */
816 0 : char **papszMD = RPCInfoToMD( &(psInfo->sRPC) );
817 : CPLXMLNode *psMD= CPLCreateXMLNode( psTree, CXT_Element,
818 0 : "Metadata" );
819 :
820 0 : for( int i = 0; papszMD != NULL && papszMD[i] != NULL; i++ )
821 : {
822 : const char *pszRawValue;
823 : char *pszKey;
824 : CPLXMLNode *psMDI;
825 :
826 0 : pszRawValue = CPLParseNameValue( papszMD[i], &pszKey );
827 :
828 0 : psMDI = CPLCreateXMLNode( psMD, CXT_Element, "MDI" );
829 0 : CPLSetXMLValue( psMDI, "#key", pszKey );
830 0 : CPLCreateXMLNode( psMDI, CXT_Text, pszRawValue );
831 :
832 0 : CPLFree( pszKey );
833 : }
834 :
835 0 : CSLDestroy( papszMD );
836 :
837 0 : return psTree;
838 : }
839 :
840 : /************************************************************************/
841 : /* GDALDeserializeRPCTransformer() */
842 : /************************************************************************/
843 :
844 0 : void *GDALDeserializeRPCTransformer( CPLXMLNode *psTree )
845 :
846 : {
847 : void *pResult;
848 0 : char **papszOptions = NULL;
849 :
850 : /* -------------------------------------------------------------------- */
851 : /* Collect metadata. */
852 : /* -------------------------------------------------------------------- */
853 0 : char **papszMD = NULL;
854 : CPLXMLNode *psMDI, *psMetadata;
855 : GDALRPCInfo sRPC;
856 :
857 0 : psMetadata = CPLGetXMLNode( psTree, "Metadata" );
858 :
859 0 : if( psMetadata->eType != CXT_Element
860 : || !EQUAL(psMetadata->pszValue,"Metadata") )
861 0 : return NULL;
862 :
863 0 : for( psMDI = psMetadata->psChild; psMDI != NULL;
864 : psMDI = psMDI->psNext )
865 : {
866 0 : if( !EQUAL(psMDI->pszValue,"MDI")
867 : || psMDI->eType != CXT_Element
868 : || psMDI->psChild == NULL
869 : || psMDI->psChild->psNext == NULL
870 : || psMDI->psChild->eType != CXT_Attribute
871 : || psMDI->psChild->psChild == NULL )
872 0 : continue;
873 :
874 : papszMD =
875 : CSLSetNameValue( papszMD,
876 : psMDI->psChild->psChild->pszValue,
877 0 : psMDI->psChild->psNext->pszValue );
878 : }
879 :
880 0 : if( !GDALExtractRPCInfo( papszMD, &sRPC ) )
881 : {
882 : CPLError( CE_Failure, CPLE_AppDefined,
883 0 : "Failed to reconstitute RPC transformer." );
884 0 : CSLDestroy( papszMD );
885 0 : return NULL;
886 : }
887 :
888 0 : CSLDestroy( papszMD );
889 :
890 : /* -------------------------------------------------------------------- */
891 : /* Get other flags. */
892 : /* -------------------------------------------------------------------- */
893 : double dfPixErrThreshold;
894 : int bReversed;
895 :
896 0 : bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
897 :
898 : dfPixErrThreshold =
899 0 : CPLAtof(CPLGetXMLValue(psTree,"PixErrThreshold","0.25"));
900 :
901 : papszOptions = CSLSetNameValue( papszOptions, "RPC_HEIGHT",
902 0 : CPLGetXMLValue(psTree,"HeightOffset","0"));
903 : papszOptions = CSLSetNameValue( papszOptions, "RPC_HEIGHT_SCALE",
904 0 : CPLGetXMLValue(psTree,"HeightScale","1"));
905 0 : const char* pszDEMPath = CPLGetXMLValue(psTree,"DEMPath",NULL);
906 0 : if (pszDEMPath != NULL)
907 : papszOptions = CSLSetNameValue( papszOptions, "RPC_DEM",
908 0 : pszDEMPath);
909 :
910 : /* -------------------------------------------------------------------- */
911 : /* Generate transformation. */
912 : /* -------------------------------------------------------------------- */
913 : pResult = GDALCreateRPCTransformer( &sRPC, bReversed, dfPixErrThreshold,
914 0 : papszOptions );
915 :
916 0 : CSLDestroy( papszOptions );
917 :
918 0 : return pResult;
919 : }
|
