1 : /******************************************************************************
2 : * $Id: gdal_rpc.cpp 24099 2012-03-09 09:24:57Z bishop $
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 24099 2012-03-09 09:24:57Z bishop $");
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 56 : static void RPCComputeTerms( double dfLong, double dfLat, double dfHeight,
149 : double *padfTerms )
150 :
151 : {
152 56 : padfTerms[0] = 1.0;
153 56 : padfTerms[1] = dfLong;
154 56 : padfTerms[2] = dfLat;
155 56 : padfTerms[3] = dfHeight;
156 56 : padfTerms[4] = dfLong * dfLat;
157 56 : padfTerms[5] = dfLong * dfHeight;
158 56 : padfTerms[6] = dfLat * dfHeight;
159 56 : padfTerms[7] = dfLong * dfLong;
160 56 : padfTerms[8] = dfLat * dfLat;
161 56 : padfTerms[9] = dfHeight * dfHeight;
162 :
163 56 : padfTerms[10] = dfLong * dfLat * dfHeight;
164 56 : padfTerms[11] = dfLong * dfLong * dfLong;
165 56 : padfTerms[12] = dfLong * dfLat * dfLat;
166 56 : padfTerms[13] = dfLong * dfHeight * dfHeight;
167 56 : padfTerms[14] = dfLong * dfLong * dfLat;
168 56 : padfTerms[15] = dfLat * dfLat * dfLat;
169 56 : padfTerms[16] = dfLat * dfHeight * dfHeight;
170 56 : padfTerms[17] = dfLong * dfLong * dfHeight;
171 56 : padfTerms[18] = dfLat * dfLat * dfHeight;
172 56 : padfTerms[19] = dfHeight * dfHeight * dfHeight;
173 56 : }
174 :
175 : /************************************************************************/
176 : /* RPCEvaluate() */
177 : /************************************************************************/
178 :
179 224 : static double RPCEvaluate( double *padfTerms, double *padfCoefs )
180 :
181 : {
182 224 : double dfSum = 0.0;
183 : int i;
184 :
185 4704 : for( i = 0; i < 20; i++ )
186 4480 : dfSum += padfTerms[i] * padfCoefs[i];
187 :
188 224 : return dfSum;
189 : }
190 :
191 : /************************************************************************/
192 : /* RPCTransformPoint() */
193 : /************************************************************************/
194 :
195 56 : static void RPCTransformPoint( GDALRPCInfo *psRPC,
196 : double dfLong, double dfLat, double dfHeight,
197 : 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 56 : adfTerms );
208 :
209 : dfResultX = RPCEvaluate( adfTerms, psRPC->adfSAMP_NUM_COEFF )
210 56 : / RPCEvaluate( adfTerms, psRPC->adfSAMP_DEN_COEFF );
211 :
212 : dfResultY = RPCEvaluate( adfTerms, psRPC->adfLINE_NUM_COEFF )
213 56 : / RPCEvaluate( adfTerms, psRPC->adfLINE_DEN_COEFF );
214 :
215 56 : *pdfPixel = dfResultX * psRPC->dfSAMP_SCALE + psRPC->dfSAMP_OFF;
216 56 : *pdfLine = dfResultY * psRPC->dfLINE_SCALE + psRPC->dfLINE_OFF;
217 56 : }
218 :
219 : /************************************************************************/
220 : /* ==================================================================== */
221 : /* GDALRPCTransformer */
222 : /* ==================================================================== */
223 : /************************************************************************/
224 :
225 : /*! DEM Resampling Algorithm */
226 : typedef enum {
227 : /*! Nearest neighbour (select on one input pixel) */ DRA_NearestNeighbour=0,
228 : /*! Bilinear (2x2 kernel) */ DRA_Bilinear=1,
229 : /*! Cubic Convolution Approximation (4x4 kernel) */ DRA_Cubic=2
230 : } DEMResampleAlg;
231 :
232 : typedef struct {
233 :
234 : GDALTransformerInfo sTI;
235 :
236 : GDALRPCInfo sRPC;
237 :
238 : double adfPLToLatLongGeoTransform[6];
239 :
240 : int bReversed;
241 :
242 : double dfPixErrThreshold;
243 :
244 : double dfHeightOffset;
245 :
246 : double dfHeightScale;
247 :
248 : char *pszDEMPath;
249 :
250 : DEMResampleAlg eResampleAlg;
251 :
252 : int bHasTriedOpeningDS;
253 : GDALDataset *poDS;
254 :
255 : OGRCoordinateTransformation *poCT;
256 :
257 : double adfGeoTransform[6];
258 : double adfReverseGeoTransform[6];
259 : } GDALRPCTransformInfo;
260 :
261 : /************************************************************************/
262 : /* GDALCreateRPCTransformer() */
263 : /************************************************************************/
264 :
265 : /**
266 : * Create an RPC based transformer.
267 : *
268 : * The geometric sensor model describing the physical relationship between
269 : * image coordinates and ground coordinate is known as a Rigorous Projection
270 : * Model. A Rigorous Projection Model expresses the mapping of the image space
271 : * coordinates of rows and columns (r,c) onto the object space reference
272 : * surface geodetic coordinates (long, lat, height).
273 : *
274 : * RPC supports a generic description of the Rigorous Projection Models. The
275 : * approximation used by GDAL (RPC00) is a set of rational polynomials exp
276 : * ressing the normalized row and column values, (rn , cn), as a function of
277 : * normalized geodetic latitude, longitude, and height, (P, L, H), given a
278 : * set of normalized polynomial coefficients (LINE_NUM_COEF_n, LINE_DEN_COEF_n,
279 : * SAMP_NUM_COEF_n, SAMP_DEN_COEF_n). Normalized values, rather than actual
280 : * values are used in order to minimize introduction of errors during the
281 : * calculations. The transformation between row and column values (r,c), and
282 : * normalized row and column values (rn, cn), and between the geodetic
283 : * latitude, longitude, and height and normalized geodetic latitude,
284 : * longitude, and height (P, L, H), is defined by a set of normalizing
285 : * translations (offsets) and scales that ensure all values are contained i
286 : * the range -1 to +1.
287 : *
288 : * This function creates a GDALTransformFunc compatible transformer
289 : * for going between image pixel/line and long/lat/height coordinates
290 : * using RPCs. The RPCs are provided in a GDALRPCInfo structure which is
291 : * normally read from metadata using GDALExtractRPCInfo().
292 : *
293 : * GDAL RPC Metadata has the following entries (also described in GDAL RFC 22
294 : * and the GeoTIFF RPC document http://geotiff.maptools.org/rpc_prop.html.
295 : *
296 : * <ul>
297 : * <li>ERR_BIAS: Error - Bias. The RMS bias error in meters per horizontal axis of all points in the image (-1.0 if unknown)
298 : * <li>ERR_RAND: Error - Random. RMS random error in meters per horizontal axis of each point in the image (-1.0 if unknown)
299 : * <li>LINE_OFF: Line Offset
300 : * <li>SAMP_OFF: Sample Offset
301 : * <li>LAT_OFF: Geodetic Latitude Offset
302 : * <li>LONG_OFF: Geodetic Longitude Offset
303 : * <li>HEIGHT_OFF: Geodetic Height Offset
304 : * <li>LINE_SCALE: Line Scale
305 : * <li>SAMP_SCALE: Sample Scale
306 : * <li>LAT_SCALE: Geodetic Latitude Scale
307 : * <li>LONG_SCALE: Geodetic Longitude Scale
308 : * <li>HEIGHT_SCALE: Geodetic Height Scale
309 : * <li>LINE_NUM_COEFF (1-20): Line Numerator Coefficients. Twenty coefficients for the polynomial in the Numerator of the rn equation. (space separated)
310 : * <li>LINE_DEN_COEFF (1-20): Line Denominator Coefficients. Twenty coefficients for the polynomial in the Denominator of the rn equation. (space separated)
311 : * <li>SAMP_NUM_COEFF (1-20): Sample Numerator Coefficients. Twenty coefficients for the polynomial in the Numerator of the cn equation. (space separated)
312 : * <li>SAMP_DEN_COEFF (1-20): Sample Denominator Coefficients. Twenty coefficients for the polynomial in the Denominator of the cn equation. (space separated)
313 : * </ul>
314 : *
315 : * The transformer normally maps from pixel/line/height to long/lat/height space
316 : * as a forward transformation though in RPC terms that would be considered
317 : * an inverse transformation (and is solved by iterative approximation using
318 : * long/lat/height to pixel/line transformations). The default direction can
319 : * be reversed by passing bReversed=TRUE.
320 : *
321 : * The iterative solution of pixel/line
322 : * to lat/long/height is currently run for up to 10 iterations or until
323 : * the apparent error is less than dfPixErrThreshold pixels. Passing zero
324 : * will not avoid all error, but will cause the operation to run for the maximum
325 : * number of iterations.
326 : *
327 : * Additional options to the transformer can be supplied in papszOptions.
328 : *
329 : * Options:
330 : *
331 : * <ul>
332 : * <li> RPC_HEIGHT: a fixed height offset to be applied to all points passed
333 : * in. In this situation the Z passed into the transformation function is
334 : * assumed to be height above ground, and the RPC_HEIGHT is assumed to be
335 : * an average height above sea level for ground in the target scene.
336 : *
337 : * <li> RPC_HEIGHT_SCALE: a factor used to multiply heights above ground.
338 : * Usefull when elevation offsets of the DEM are not expressed in meters. (GDAL >= 1.8.0)
339 : *
340 : * <li> RPC_DEM: the name of a GDAL dataset (a DEM file typically) used to
341 : * extract elevation offsets from. In this situation the Z passed into the
342 : * transformation function is assumed to be height above ground. This option
343 : * should be used in replacement of RPC_HEIGHT to provide a way of defining
344 : * a non uniform ground for the target scene (GDAL >= 1.8.0)
345 : *
346 : * <li> RPC_DEMINTERPOLATION: the DEM interpolation (near, bilinear or cubic)
347 : * </ul>
348 : *
349 : * @param psRPCInfo Definition of the RPC parameters.
350 : *
351 : * @param bReversed If true "forward" transformation will be lat/long to pixel/line instead of the normal pixel/line to lat/long.
352 : *
353 : * @param dfPixErrThreshold the error (measured in pixels) allowed in the
354 : * iterative solution of pixel/line to lat/long computations (the other way
355 : * is always exact given the equations).
356 : *
357 : * @param papszOptions Other transformer options (ie. RPC_HEIGHT=<z>).
358 : *
359 : * @return transformer callback data (deallocate with GDALDestroyTransformer()).
360 : */
361 :
362 6 : void *GDALCreateRPCTransformer( GDALRPCInfo *psRPCInfo, int bReversed,
363 : double dfPixErrThreshold,
364 : char **papszOptions )
365 :
366 : {
367 : GDALRPCTransformInfo *psTransform;
368 :
369 : /* -------------------------------------------------------------------- */
370 : /* Initialize core info. */
371 : /* -------------------------------------------------------------------- */
372 : psTransform = (GDALRPCTransformInfo *)
373 6 : CPLCalloc(sizeof(GDALRPCTransformInfo),1);
374 :
375 6 : memcpy( &(psTransform->sRPC), psRPCInfo, sizeof(GDALRPCInfo) );
376 6 : psTransform->bReversed = bReversed;
377 6 : psTransform->dfPixErrThreshold = dfPixErrThreshold;
378 6 : psTransform->dfHeightOffset = 0.0;
379 6 : psTransform->dfHeightScale = 1.0;
380 :
381 6 : strcpy( psTransform->sTI.szSignature, "GTI" );
382 6 : psTransform->sTI.pszClassName = "GDALRPCTransformer";
383 6 : psTransform->sTI.pfnTransform = GDALRPCTransform;
384 6 : psTransform->sTI.pfnCleanup = GDALDestroyRPCTransformer;
385 6 : psTransform->sTI.pfnSerialize = GDALSerializeRPCTransformer;
386 :
387 : /* -------------------------------------------------------------------- */
388 : /* Do we have a "average height" that we want to consider all */
389 : /* elevations to be relative to? */
390 : /* -------------------------------------------------------------------- */
391 6 : const char *pszHeight = CSLFetchNameValue( papszOptions, "RPC_HEIGHT" );
392 6 : if( pszHeight != NULL )
393 2 : psTransform->dfHeightOffset = CPLAtof(pszHeight);
394 :
395 : /* -------------------------------------------------------------------- */
396 : /* The "height scale" */
397 : /* -------------------------------------------------------------------- */
398 6 : const char *pszHeightScale = CSLFetchNameValue( papszOptions, "RPC_HEIGHT_SCALE" );
399 6 : if( pszHeightScale != NULL )
400 2 : psTransform->dfHeightScale = CPLAtof(pszHeightScale);
401 :
402 : /* -------------------------------------------------------------------- */
403 : /* The DEM file name */
404 : /* -------------------------------------------------------------------- */
405 6 : const char *pszDEMPath = CSLFetchNameValue( papszOptions, "RPC_DEM" );
406 6 : if( pszDEMPath != NULL )
407 2 : psTransform->pszDEMPath = CPLStrdup(pszDEMPath);
408 :
409 : /* -------------------------------------------------------------------- */
410 : /* The DEM interpolation */
411 : /* -------------------------------------------------------------------- */
412 6 : const char *pszDEMInterpolation = CSLFetchNameValueDef( papszOptions, "RPC_DEMINTERPOLATION", "bilinear" );
413 6 : if(EQUAL(pszDEMInterpolation, "near" ))
414 0 : psTransform->eResampleAlg = DRA_NearestNeighbour;
415 6 : else if(EQUAL(pszDEMInterpolation, "bilinear" ))
416 6 : psTransform->eResampleAlg = DRA_Bilinear;
417 0 : else if(EQUAL(pszDEMInterpolation, "cubic" ))
418 0 : psTransform->eResampleAlg = DRA_Cubic;
419 : else
420 0 : psTransform->eResampleAlg = DRA_Bilinear;
421 :
422 : /* -------------------------------------------------------------------- */
423 : /* Establish a reference point for calcualating an affine */
424 : /* geotransform approximate transformation. */
425 : /* -------------------------------------------------------------------- */
426 6 : double adfGTFromLL[6], dfRefPixel = -1.0, dfRefLine = -1.0;
427 6 : double dfRefLong = 0.0, dfRefLat = 0.0;
428 :
429 6 : if( psRPCInfo->dfMIN_LONG != -180 || psRPCInfo->dfMAX_LONG != 180 )
430 : {
431 0 : dfRefLong = (psRPCInfo->dfMIN_LONG + psRPCInfo->dfMAX_LONG) * 0.5;
432 0 : dfRefLat = (psRPCInfo->dfMIN_LAT + psRPCInfo->dfMAX_LAT ) * 0.5;
433 :
434 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat, 0.0,
435 0 : &dfRefPixel, &dfRefLine );
436 : }
437 :
438 : // Try with scale and offset if we don't can't use bounds or
439 : // the results seem daft.
440 6 : if( dfRefPixel < 0.0 || dfRefLine < 0.0
441 : || dfRefPixel > 100000 || dfRefLine > 100000 )
442 : {
443 6 : dfRefLong = psRPCInfo->dfLONG_OFF;
444 6 : dfRefLat = psRPCInfo->dfLAT_OFF;
445 :
446 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat, 0.0,
447 6 : &dfRefPixel, &dfRefLine );
448 : }
449 :
450 : /* -------------------------------------------------------------------- */
451 : /* Transform nearby locations to establish affine direction */
452 : /* vectors. */
453 : /* -------------------------------------------------------------------- */
454 6 : double dfRefPixelDelta, dfRefLineDelta, dfLLDelta = 0.0001;
455 :
456 : RPCTransformPoint( psRPCInfo, dfRefLong+dfLLDelta, dfRefLat, 0.0,
457 6 : &dfRefPixelDelta, &dfRefLineDelta );
458 6 : adfGTFromLL[1] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
459 6 : adfGTFromLL[2] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
460 :
461 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat+dfLLDelta, 0.0,
462 6 : &dfRefPixelDelta, &dfRefLineDelta );
463 6 : adfGTFromLL[4] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
464 6 : adfGTFromLL[5] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
465 :
466 : adfGTFromLL[0] = dfRefPixel
467 6 : - adfGTFromLL[1] * dfRefLong - adfGTFromLL[2] * dfRefLat;
468 : adfGTFromLL[3] = dfRefLine
469 6 : - adfGTFromLL[4] * dfRefLong - adfGTFromLL[5] * dfRefLat;
470 :
471 6 : GDALInvGeoTransform( adfGTFromLL, psTransform->adfPLToLatLongGeoTransform);
472 :
473 6 : return psTransform;
474 : }
475 :
476 : /************************************************************************/
477 : /* GDALDestroyReprojectionTransformer() */
478 : /************************************************************************/
479 :
480 6 : void GDALDestroyRPCTransformer( void *pTransformAlg )
481 :
482 : {
483 6 : GDALRPCTransformInfo *psTransform = (GDALRPCTransformInfo *) pTransformAlg;
484 :
485 6 : CPLFree( psTransform->pszDEMPath );
486 :
487 6 : if(psTransform->poDS)
488 2 : GDALClose(psTransform->poDS);
489 6 : if(psTransform->poCT)
490 2 : OCTDestroyCoordinateTransformation((OGRCoordinateTransformationH)psTransform->poCT);
491 :
492 6 : CPLFree( pTransformAlg );
493 6 : }
494 :
495 : /************************************************************************/
496 : /* RPCInverseTransformPoint() */
497 : /************************************************************************/
498 :
499 : static void
500 10 : RPCInverseTransformPoint( GDALRPCTransformInfo *psTransform,
501 : double dfPixel, double dfLine, double dfHeight,
502 : double *pdfLong, double *pdfLat )
503 :
504 : {
505 : double dfResultX, dfResultY;
506 : int iIter;
507 10 : GDALRPCInfo *psRPC = &(psTransform->sRPC);
508 :
509 : /* -------------------------------------------------------------------- */
510 : /* Compute an initial approximation based on linear */
511 : /* interpolation from our reference point. */
512 : /* -------------------------------------------------------------------- */
513 10 : dfResultX = psTransform->adfPLToLatLongGeoTransform[0]
514 10 : + psTransform->adfPLToLatLongGeoTransform[1] * dfPixel
515 10 : + psTransform->adfPLToLatLongGeoTransform[2] * dfLine;
516 :
517 10 : dfResultY = psTransform->adfPLToLatLongGeoTransform[3]
518 10 : + psTransform->adfPLToLatLongGeoTransform[4] * dfPixel
519 10 : + psTransform->adfPLToLatLongGeoTransform[5] * dfLine;
520 :
521 : /* -------------------------------------------------------------------- */
522 : /* Now iterate, trying to find a closer LL location that will */
523 : /* back transform to the indicated pixel and line. */
524 : /* -------------------------------------------------------------------- */
525 10 : double dfPixelDeltaX=0.0, dfPixelDeltaY=0.0;
526 :
527 30 : for( iIter = 0; iIter < 10; iIter++ )
528 : {
529 : double dfBackPixel, dfBackLine;
530 :
531 : RPCTransformPoint( psRPC, dfResultX, dfResultY, dfHeight,
532 30 : &dfBackPixel, &dfBackLine );
533 :
534 30 : dfPixelDeltaX = dfBackPixel - dfPixel;
535 30 : dfPixelDeltaY = dfBackLine - dfLine;
536 :
537 : dfResultX = dfResultX
538 30 : - dfPixelDeltaX * psTransform->adfPLToLatLongGeoTransform[1]
539 30 : - dfPixelDeltaY * psTransform->adfPLToLatLongGeoTransform[2];
540 : dfResultY = dfResultY
541 30 : - dfPixelDeltaX * psTransform->adfPLToLatLongGeoTransform[4]
542 30 : - dfPixelDeltaY * psTransform->adfPLToLatLongGeoTransform[5];
543 :
544 30 : if( ABS(dfPixelDeltaX) < psTransform->dfPixErrThreshold
545 : && ABS(dfPixelDeltaY) < psTransform->dfPixErrThreshold )
546 : {
547 10 : iIter = -1;
548 : //CPLDebug( "RPC", "Converged!" );
549 10 : break;
550 : }
551 :
552 : }
553 :
554 10 : if( iIter != -1 )
555 : CPLDebug( "RPC", "Iterations %d: Got: %g,%g Offset=%g,%g",
556 : iIter,
557 : dfResultX, dfResultY,
558 0 : dfPixelDeltaX, dfPixelDeltaY );
559 :
560 10 : *pdfLong = dfResultX;
561 10 : *pdfLat = dfResultY;
562 10 : }
563 :
564 :
565 0 : double BiCubicKernel(double dfVal)
566 : {
567 0 : if ( dfVal > 2.0 )
568 0 : return 0.0;
569 :
570 : double a, b, c, d;
571 0 : double xm1 = dfVal - 1.0;
572 0 : double xp1 = dfVal + 1.0;
573 0 : double xp2 = dfVal + 2.0;
574 :
575 0 : a = ( xp2 <= 0.0 ) ? 0.0 : xp2 * xp2 * xp2;
576 0 : b = ( xp1 <= 0.0 ) ? 0.0 : xp1 * xp1 * xp1;
577 0 : c = ( dfVal <= 0.0 ) ? 0.0 : dfVal * dfVal * dfVal;
578 0 : d = ( xm1 <= 0.0 ) ? 0.0 : xm1 * xm1 * xm1;
579 :
580 0 : return ( 0.16666666666666666667 * ( a - ( 4.0 * b ) + ( 6.0 * c ) - ( 4.0 * d ) ) );
581 : }
582 :
583 : /************************************************************************/
584 : /* GDALRPCTransform() */
585 : /************************************************************************/
586 :
587 16 : int GDALRPCTransform( void *pTransformArg, int bDstToSrc,
588 : int nPointCount,
589 : double *padfX, double *padfY, double *padfZ,
590 : int *panSuccess )
591 :
592 : {
593 16 : VALIDATE_POINTER1( pTransformArg, "GDALRPCTransform", 0 );
594 :
595 16 : GDALRPCTransformInfo *psTransform = (GDALRPCTransformInfo *) pTransformArg;
596 16 : GDALRPCInfo *psRPC = &(psTransform->sRPC);
597 : int i;
598 :
599 16 : if( psTransform->bReversed )
600 0 : bDstToSrc = !bDstToSrc;
601 :
602 16 : int bands[1] = {1};
603 16 : int nRasterXSize = 0, nRasterYSize = 0;
604 :
605 : /* -------------------------------------------------------------------- */
606 : /* Lazy opening of the optionnal DEM file. */
607 : /* -------------------------------------------------------------------- */
608 16 : if(psTransform->pszDEMPath != NULL &&
609 : psTransform->bHasTriedOpeningDS == FALSE)
610 : {
611 2 : int bIsValid = FALSE;
612 2 : psTransform->bHasTriedOpeningDS = TRUE;
613 : psTransform->poDS = (GDALDataset *)
614 2 : GDALOpen( psTransform->pszDEMPath, GA_ReadOnly );
615 2 : if(psTransform->poDS != NULL && psTransform->poDS->GetRasterCount() >= 1)
616 : {
617 2 : const char* pszSpatialRef = psTransform->poDS->GetProjectionRef();
618 2 : if (pszSpatialRef != NULL && pszSpatialRef[0] != '\0')
619 : {
620 : OGRSpatialReference* poWGSSpaRef =
621 2 : new OGRSpatialReference(SRS_WKT_WGS84);
622 : OGRSpatialReference* poDSSpaRef =
623 4 : new OGRSpatialReference(pszSpatialRef);
624 2 : if(!poWGSSpaRef->IsSame(poDSSpaRef))
625 : psTransform->poCT =OGRCreateCoordinateTransformation(
626 2 : poWGSSpaRef, poDSSpaRef );
627 2 : delete poWGSSpaRef;
628 2 : delete poDSSpaRef;
629 : }
630 :
631 4 : if (psTransform->poDS->GetGeoTransform(
632 2 : psTransform->adfGeoTransform) == CE_None &&
633 : GDALInvGeoTransform( psTransform->adfGeoTransform,
634 : psTransform->adfReverseGeoTransform ))
635 : {
636 2 : bIsValid = TRUE;
637 : }
638 : }
639 :
640 2 : if (!bIsValid && psTransform->poDS != NULL)
641 : {
642 0 : GDALClose(psTransform->poDS);
643 0 : psTransform->poDS = NULL;
644 : }
645 : }
646 16 : if (psTransform->poDS)
647 : {
648 4 : nRasterXSize = psTransform->poDS->GetRasterXSize();
649 4 : nRasterYSize = psTransform->poDS->GetRasterYSize();
650 : }
651 :
652 : /* -------------------------------------------------------------------- */
653 : /* The simple case is transforming from lat/long to pixel/line. */
654 : /* Just apply the equations directly. */
655 : /* -------------------------------------------------------------------- */
656 16 : if( bDstToSrc )
657 : {
658 16 : for( i = 0; i < nPointCount; i++ )
659 : {
660 8 : if(psTransform->poDS)
661 : {
662 : double dfX, dfY;
663 : //check if dem is not in WGS84 and transform points padfX[i], padfY[i]
664 2 : if(psTransform->poCT)
665 : {
666 2 : double dfXOrig = padfX[i];
667 2 : double dfYOrig = padfY[i];
668 2 : double dfZOrig = padfZ[i];
669 4 : if (!psTransform->poCT->Transform(
670 2 : 1, &dfXOrig, &dfYOrig, &dfZOrig))
671 : {
672 0 : panSuccess[i] = FALSE;
673 0 : continue;
674 : }
675 : GDALApplyGeoTransform( psTransform->adfReverseGeoTransform,
676 2 : dfXOrig, dfYOrig, &dfX, &dfY );
677 : }
678 : else
679 : GDALApplyGeoTransform( psTransform->adfReverseGeoTransform,
680 0 : padfX[i], padfY[i], &dfX, &dfY );
681 2 : int dX = int(dfX);
682 2 : int dY = int(dfY);
683 :
684 2 : if (!(dX >= 0 && dY >= 0 &&
685 : dX+2 <= nRasterXSize && dY+2 <= nRasterYSize))
686 : {
687 0 : panSuccess[i] = FALSE;
688 0 : continue;
689 : }
690 :
691 2 : double dfDEMH(0);
692 2 : double dfDeltaX = dfX - dX;
693 2 : double dfDeltaY = dfY - dY;
694 :
695 2 : if(psTransform->eResampleAlg == DRA_Cubic)
696 : {
697 0 : int dXNew = dX - 1;
698 0 : int dYNew = dY - 1;
699 0 : if (!(dXNew >= 0 && dYNew >= 0 && dXNew + 4 <= nRasterXSize && dYNew + 4 <= nRasterYSize))
700 : {
701 0 : panSuccess[i] = FALSE;
702 0 : continue;
703 : }
704 : //cubic interpolation
705 0 : int adElevData[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
706 : CPLErr eErr = psTransform->poDS->RasterIO(GF_Read, dXNew, dYNew, 4, 4,
707 : &adElevData, 4, 4,
708 0 : GDT_Int32, 1, bands, 0, 0, 0);
709 0 : if(eErr != CE_None)
710 : {
711 0 : panSuccess[i] = FALSE;
712 0 : continue;
713 : }
714 :
715 0 : double dfSumH(0);
716 0 : for ( int i = 0; i < 5; i++ )
717 : {
718 : // Loop across the X axis
719 0 : for ( int j = 0; j < 5; j++ )
720 : {
721 : // Calculate the weight for the specified pixel according
722 : // to the bicubic b-spline kernel we're using for
723 : // interpolation
724 0 : int dKernIndX = j - 1;
725 0 : int dKernIndY = i - 1;
726 0 : double dfPixelWeight = BiCubicKernel(dKernIndX - dfDeltaX) * BiCubicKernel(dKernIndY - dfDeltaY);
727 :
728 : // Create a sum of all values
729 : // adjusted for the pixel's calculated weight
730 0 : dfSumH += adElevData[j + i * 4] * dfPixelWeight;
731 : }
732 : }
733 0 : dfDEMH = dfSumH;
734 : }
735 2 : else if(psTransform->eResampleAlg == DRA_Bilinear)
736 : {
737 2 : if (!(dX >= 0 && dY >= 0 && dX + 2 <= nRasterXSize && dY + 2 <= nRasterYSize))
738 : {
739 0 : panSuccess[i] = FALSE;
740 0 : continue;
741 : }
742 : //bilinear interpolation
743 2 : int anElevData[4] = {0,0,0,0};
744 : CPLErr eErr = psTransform->poDS->RasterIO(GF_Read, dX, dY, 2, 2,
745 : &anElevData, 2, 2,
746 2 : GDT_Int32, 1, bands, 0, 0, 0);
747 2 : if(eErr != CE_None)
748 : {
749 0 : panSuccess[i] = FALSE;
750 0 : continue;
751 : }
752 2 : double dfDeltaX1 = 1.0 - dfDeltaX;
753 2 : double dfDeltaY1 = 1.0 - dfDeltaY;
754 :
755 2 : double dfXZ1 = anElevData[0] * dfDeltaX1 + anElevData[1] * dfDeltaX;
756 2 : double dfXZ2 = anElevData[2] * dfDeltaX1 + anElevData[3] * dfDeltaX;
757 2 : double dfYZ = dfXZ1 * dfDeltaY1 + dfXZ2 * dfDeltaY;
758 2 : dfDEMH = dfYZ;
759 : }
760 : else
761 : {
762 0 : if (!(dX >= 0 && dY >= 0 && dX <= nRasterXSize && dY <= nRasterYSize))
763 : {
764 0 : panSuccess[i] = FALSE;
765 0 : continue;
766 : }
767 : CPLErr eErr = psTransform->poDS->RasterIO(GF_Read, dX, dY, 1, 1,
768 : &dfDEMH, 1, 1,
769 0 : GDT_Int32, 1, bands, 0, 0, 0);
770 0 : if(eErr != CE_None)
771 : {
772 0 : panSuccess[i] = FALSE;
773 0 : continue;
774 : }
775 : }
776 :
777 : RPCTransformPoint( psRPC, padfX[i], padfY[i],
778 2 : padfZ[i] + (psTransform->dfHeightOffset + dfDEMH) *
779 : psTransform->dfHeightScale,
780 4 : padfX + i, padfY + i );
781 : }
782 : else
783 : RPCTransformPoint( psRPC, padfX[i], padfY[i],
784 6 : padfZ[i] + psTransform->dfHeightOffset *
785 : psTransform->dfHeightScale,
786 12 : padfX + i, padfY + i );
787 8 : panSuccess[i] = TRUE;
788 : }
789 :
790 8 : return TRUE;
791 : }
792 :
793 : /* -------------------------------------------------------------------- */
794 : /* Compute the inverse (pixel/line/height to lat/long). This */
795 : /* function uses an iterative method from an initial linear */
796 : /* approximation. */
797 : /* -------------------------------------------------------------------- */
798 16 : for( i = 0; i < nPointCount; i++ )
799 : {
800 : double dfResultX, dfResultY;
801 :
802 8 : if(psTransform->poDS)
803 : {
804 : RPCInverseTransformPoint( psTransform, padfX[i], padfY[i],
805 2 : padfZ[i] + psTransform->dfHeightOffset *
806 : psTransform->dfHeightScale,
807 2 : &dfResultX, &dfResultY );
808 :
809 : double dfX, dfY;
810 : //check if dem is not in WGS84 and transform points padfX[i], padfY[i]
811 2 : if(psTransform->poCT)
812 : {
813 2 : double dfZ = 0;
814 2 : if (!psTransform->poCT->Transform(1, &dfResultX, &dfResultY, &dfZ))
815 : {
816 0 : panSuccess[i] = FALSE;
817 0 : continue;
818 : }
819 : }
820 :
821 : GDALApplyGeoTransform( psTransform->adfReverseGeoTransform,
822 2 : dfResultX, dfResultY, &dfX, &dfY );
823 2 : int dX = int(dfX);
824 2 : int dY = int(dfY);
825 :
826 2 : double dfDEMH(0);
827 2 : double dfDeltaX = dfX - dX;
828 2 : double dfDeltaY = dfY - dY;
829 :
830 2 : if(psTransform->eResampleAlg == DRA_Cubic)
831 : {
832 0 : int dXNew = dX - 1;
833 0 : int dYNew = dY - 1;
834 0 : if (!(dXNew >= 0 && dYNew >= 0 && dXNew + 4 <= nRasterXSize && dYNew + 4 <= nRasterYSize))
835 : {
836 0 : panSuccess[i] = FALSE;
837 0 : continue;
838 : }
839 : //cubic interpolation
840 0 : int adElevData[16] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
841 : CPLErr eErr = psTransform->poDS->RasterIO(GF_Read, dXNew, dYNew, 4, 4,
842 : &adElevData, 4, 4,
843 0 : GDT_Int32, 1, bands, 0, 0, 0);
844 0 : if(eErr != CE_None)
845 : {
846 0 : panSuccess[i] = FALSE;
847 0 : continue;
848 : }
849 :
850 0 : double dfSumH(0);
851 0 : for ( int i = 0; i < 5; i++ )
852 : {
853 : // Loop across the X axis
854 0 : for ( int j = 0; j < 5; j++ )
855 : {
856 : // Calculate the weight for the specified pixel according
857 : // to the bicubic b-spline kernel we're using for
858 : // interpolation
859 0 : int dKernIndX = j - 1;
860 0 : int dKernIndY = i - 1;
861 0 : double dfPixelWeight = BiCubicKernel(dKernIndX - dfDeltaX) * BiCubicKernel(dKernIndY - dfDeltaY);
862 :
863 : // Create a sum of all values
864 : // adjusted for the pixel's calculated weight
865 0 : dfSumH += adElevData[j + i * 4] * dfPixelWeight;
866 : }
867 : }
868 0 : dfDEMH = dfSumH;
869 : }
870 2 : else if(psTransform->eResampleAlg == DRA_Bilinear)
871 : {
872 2 : if (!(dX >= 0 && dY >= 0 && dX + 2 <= nRasterXSize && dY + 2 <= nRasterYSize))
873 : {
874 0 : panSuccess[i] = FALSE;
875 0 : continue;
876 : }
877 : //bilinear interpolation
878 2 : int adElevData[4] = {0,0,0,0};
879 : CPLErr eErr = psTransform->poDS->RasterIO(GF_Read, dX, dY, 2, 2,
880 : &adElevData, 2, 2,
881 2 : GDT_Int32, 1, bands, 0, 0, 0);
882 2 : if(eErr != CE_None)
883 : {
884 0 : panSuccess[i] = FALSE;
885 0 : continue;
886 : }
887 2 : double dfDeltaX1 = 1.0 - dfDeltaX;
888 2 : double dfDeltaY1 = 1.0 - dfDeltaY;
889 :
890 2 : double dfXZ1 = adElevData[0] * dfDeltaX1 + adElevData[1] * dfDeltaX;
891 2 : double dfXZ2 = adElevData[2] * dfDeltaX1 + adElevData[3] * dfDeltaX;
892 2 : double dfYZ = dfXZ1 * dfDeltaY1 + dfXZ2 * dfDeltaY;
893 2 : dfDEMH = dfYZ;
894 : }
895 : else
896 : {
897 0 : if (!(dX >= 0 && dY >= 0 && dX <= nRasterXSize && dY <= nRasterYSize))
898 : {
899 0 : panSuccess[i] = FALSE;
900 0 : continue;
901 : }
902 : CPLErr eErr = psTransform->poDS->RasterIO(GF_Read, dX, dY, 1, 1,
903 : &dfDEMH, 1, 1,
904 0 : GDT_Int32, 1, bands, 0, 0, 0);
905 0 : if(eErr != CE_None)
906 : {
907 0 : panSuccess[i] = FALSE;
908 0 : continue;
909 : }
910 : }
911 :
912 : RPCInverseTransformPoint( psTransform, padfX[i], padfY[i],
913 2 : padfZ[i] + (psTransform->dfHeightOffset + dfDEMH) *
914 : psTransform->dfHeightScale,
915 2 : &dfResultX, &dfResultY );
916 : }
917 : else
918 : {
919 : RPCInverseTransformPoint( psTransform, padfX[i], padfY[i],
920 6 : padfZ[i] + psTransform->dfHeightOffset *
921 : psTransform->dfHeightScale,
922 6 : &dfResultX, &dfResultY );
923 :
924 : }
925 8 : padfX[i] = dfResultX;
926 8 : padfY[i] = dfResultY;
927 :
928 8 : panSuccess[i] = TRUE;
929 : }
930 :
931 8 : return TRUE;
932 : }
933 :
934 : /************************************************************************/
935 : /* GDALSerializeRPCTransformer() */
936 : /************************************************************************/
937 :
938 0 : CPLXMLNode *GDALSerializeRPCTransformer( void *pTransformArg )
939 :
940 : {
941 0 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeRPCTransformer", NULL );
942 :
943 : CPLXMLNode *psTree;
944 : GDALRPCTransformInfo *psInfo =
945 0 : (GDALRPCTransformInfo *)(pTransformArg);
946 :
947 0 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "RPCTransformer" );
948 :
949 : /* -------------------------------------------------------------------- */
950 : /* Serialize bReversed. */
951 : /* -------------------------------------------------------------------- */
952 : CPLCreateXMLElementAndValue(
953 : psTree, "Reversed",
954 0 : CPLString().Printf( "%d", psInfo->bReversed ) );
955 :
956 : /* -------------------------------------------------------------------- */
957 : /* Serialize Height Offset. */
958 : /* -------------------------------------------------------------------- */
959 : CPLCreateXMLElementAndValue(
960 : psTree, "HeightOffset",
961 0 : CPLString().Printf( "%.15g", psInfo->dfHeightOffset ) );
962 :
963 : /* -------------------------------------------------------------------- */
964 : /* Serialize Height Scale. */
965 : /* -------------------------------------------------------------------- */
966 0 : if (psInfo->dfHeightScale != 1.0)
967 : CPLCreateXMLElementAndValue(
968 : psTree, "HeightScale",
969 0 : CPLString().Printf( "%.15g", psInfo->dfHeightScale ) );
970 :
971 : /* -------------------------------------------------------------------- */
972 : /* Serialize DEM path. */
973 : /* -------------------------------------------------------------------- */
974 0 : if (psInfo->pszDEMPath != NULL)
975 : CPLCreateXMLElementAndValue(
976 : psTree, "DEMPath",
977 0 : CPLString().Printf( "%s", psInfo->pszDEMPath ) );
978 :
979 : /* -------------------------------------------------------------------- */
980 : /* Serialize DEM interpolation */
981 : /* -------------------------------------------------------------------- */
982 0 : CPLString soDEMInterpolation;
983 0 : switch(psInfo->eResampleAlg)
984 : {
985 : case DRA_NearestNeighbour:
986 0 : soDEMInterpolation = "near";
987 0 : break;
988 : case DRA_Cubic:
989 0 : soDEMInterpolation = "cubic";
990 0 : break;
991 : default:
992 : case DRA_Bilinear:
993 0 : soDEMInterpolation = "bilinear";
994 : }
995 : CPLCreateXMLElementAndValue(
996 0 : psTree, "DEMInterpolation", soDEMInterpolation );
997 :
998 : /* -------------------------------------------------------------------- */
999 : /* Serialize pixel error threshold. */
1000 : /* -------------------------------------------------------------------- */
1001 : CPLCreateXMLElementAndValue(
1002 : psTree, "PixErrThreshold",
1003 0 : CPLString().Printf( "%.15g", psInfo->dfPixErrThreshold ) );
1004 :
1005 : /* -------------------------------------------------------------------- */
1006 : /* RPC metadata. */
1007 : /* -------------------------------------------------------------------- */
1008 0 : char **papszMD = RPCInfoToMD( &(psInfo->sRPC) );
1009 : CPLXMLNode *psMD= CPLCreateXMLNode( psTree, CXT_Element,
1010 0 : "Metadata" );
1011 :
1012 0 : for( int i = 0; papszMD != NULL && papszMD[i] != NULL; i++ )
1013 : {
1014 : const char *pszRawValue;
1015 : char *pszKey;
1016 : CPLXMLNode *psMDI;
1017 :
1018 0 : pszRawValue = CPLParseNameValue( papszMD[i], &pszKey );
1019 :
1020 0 : psMDI = CPLCreateXMLNode( psMD, CXT_Element, "MDI" );
1021 0 : CPLSetXMLValue( psMDI, "#key", pszKey );
1022 0 : CPLCreateXMLNode( psMDI, CXT_Text, pszRawValue );
1023 :
1024 0 : CPLFree( pszKey );
1025 : }
1026 :
1027 0 : CSLDestroy( papszMD );
1028 :
1029 0 : return psTree;
1030 : }
1031 :
1032 : /************************************************************************/
1033 : /* GDALDeserializeRPCTransformer() */
1034 : /************************************************************************/
1035 :
1036 0 : void *GDALDeserializeRPCTransformer( CPLXMLNode *psTree )
1037 :
1038 : {
1039 : void *pResult;
1040 0 : char **papszOptions = NULL;
1041 :
1042 : /* -------------------------------------------------------------------- */
1043 : /* Collect metadata. */
1044 : /* -------------------------------------------------------------------- */
1045 0 : char **papszMD = NULL;
1046 : CPLXMLNode *psMDI, *psMetadata;
1047 : GDALRPCInfo sRPC;
1048 :
1049 0 : psMetadata = CPLGetXMLNode( psTree, "Metadata" );
1050 :
1051 0 : if( psMetadata == NULL
1052 : || psMetadata->eType != CXT_Element
1053 : || !EQUAL(psMetadata->pszValue,"Metadata") )
1054 0 : return NULL;
1055 :
1056 0 : for( psMDI = psMetadata->psChild; psMDI != NULL;
1057 : psMDI = psMDI->psNext )
1058 : {
1059 0 : if( !EQUAL(psMDI->pszValue,"MDI")
1060 : || psMDI->eType != CXT_Element
1061 : || psMDI->psChild == NULL
1062 : || psMDI->psChild->psNext == NULL
1063 : || psMDI->psChild->eType != CXT_Attribute
1064 : || psMDI->psChild->psChild == NULL )
1065 0 : continue;
1066 :
1067 : papszMD =
1068 : CSLSetNameValue( papszMD,
1069 : psMDI->psChild->psChild->pszValue,
1070 0 : psMDI->psChild->psNext->pszValue );
1071 : }
1072 :
1073 0 : if( !GDALExtractRPCInfo( papszMD, &sRPC ) )
1074 : {
1075 : CPLError( CE_Failure, CPLE_AppDefined,
1076 0 : "Failed to reconstitute RPC transformer." );
1077 0 : CSLDestroy( papszMD );
1078 0 : return NULL;
1079 : }
1080 :
1081 0 : CSLDestroy( papszMD );
1082 :
1083 : /* -------------------------------------------------------------------- */
1084 : /* Get other flags. */
1085 : /* -------------------------------------------------------------------- */
1086 : double dfPixErrThreshold;
1087 : int bReversed;
1088 :
1089 0 : bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
1090 :
1091 : dfPixErrThreshold =
1092 0 : CPLAtof(CPLGetXMLValue(psTree,"PixErrThreshold","0.25"));
1093 :
1094 : papszOptions = CSLSetNameValue( papszOptions, "RPC_HEIGHT",
1095 0 : CPLGetXMLValue(psTree,"HeightOffset","0"));
1096 : papszOptions = CSLSetNameValue( papszOptions, "RPC_HEIGHT_SCALE",
1097 0 : CPLGetXMLValue(psTree,"HeightScale","1"));
1098 0 : const char* pszDEMPath = CPLGetXMLValue(psTree,"DEMPath",NULL);
1099 0 : if (pszDEMPath != NULL)
1100 : papszOptions = CSLSetNameValue( papszOptions, "RPC_DEM",
1101 0 : pszDEMPath);
1102 :
1103 0 : const char* pszDEMInterpolation = CPLGetXMLValue(psTree,"DEMInterpolation", "bilinear");
1104 0 : if (pszDEMInterpolation != NULL)
1105 : papszOptions = CSLSetNameValue( papszOptions, "RPC_DEMINTERPOLATION",
1106 0 : pszDEMInterpolation);
1107 :
1108 : /* -------------------------------------------------------------------- */
1109 : /* Generate transformation. */
1110 : /* -------------------------------------------------------------------- */
1111 : pResult = GDALCreateRPCTransformer( &sRPC, bReversed, dfPixErrThreshold,
1112 0 : papszOptions );
1113 :
1114 0 : CSLDestroy( papszOptions );
1115 :
1116 0 : return pResult;
1117 : }
|
