1 : /******************************************************************************
2 : * $Id: gdal_rpc.cpp 16666 2009-03-28 12:46:49Z 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 16666 2009-03-28 12:46:49Z 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 11 : static void RPCComputeTerms( double dfLong, double dfLat, double dfHeight,
149 : double *padfTerms )
150 :
151 : {
152 11 : padfTerms[0] = 1.0;
153 11 : padfTerms[1] = dfLong;
154 11 : padfTerms[2] = dfLat;
155 11 : padfTerms[3] = dfHeight;
156 11 : padfTerms[4] = dfLong * dfLat;
157 11 : padfTerms[5] = dfLong * dfHeight;
158 11 : padfTerms[6] = dfLat * dfHeight;
159 11 : padfTerms[7] = dfLong * dfLong;
160 11 : padfTerms[8] = dfLat * dfLat;
161 11 : padfTerms[9] = dfHeight * dfHeight;
162 :
163 11 : padfTerms[10] = dfLong * dfLat * dfHeight;
164 11 : padfTerms[11] = dfLong * dfLong * dfLong;
165 11 : padfTerms[12] = dfLong * dfLat * dfLat;
166 11 : padfTerms[13] = dfLong * dfHeight * dfHeight;
167 11 : padfTerms[14] = dfLong * dfLong * dfLat;
168 11 : padfTerms[15] = dfLat * dfLat * dfLat;
169 11 : padfTerms[16] = dfLat * dfHeight * dfHeight;
170 11 : padfTerms[17] = dfLong * dfLong * dfHeight;
171 11 : padfTerms[18] = dfLat * dfLat * dfHeight;
172 11 : padfTerms[19] = dfHeight * dfHeight * dfHeight;
173 11 : }
174 :
175 : /************************************************************************/
176 : /* RPCEvaluate() */
177 : /************************************************************************/
178 :
179 44 : static double RPCEvaluate( double *padfTerms, double *padfCoefs )
180 :
181 : {
182 44 : double dfSum = 0.0;
183 : int i;
184 :
185 924 : for( i = 0; i < 20; i++ )
186 880 : dfSum += padfTerms[i] * padfCoefs[i];
187 :
188 44 : return dfSum;
189 : }
190 :
191 : /************************************************************************/
192 : /* RPCTransformPoint() */
193 : /************************************************************************/
194 :
195 11 : 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 11 : adfTerms );
208 :
209 : dfResultX = RPCEvaluate( adfTerms, psRPC->adfSAMP_NUM_COEFF )
210 11 : / RPCEvaluate( adfTerms, psRPC->adfSAMP_DEN_COEFF );
211 :
212 : dfResultY = RPCEvaluate( adfTerms, psRPC->adfLINE_NUM_COEFF )
213 11 : / RPCEvaluate( adfTerms, psRPC->adfLINE_DEN_COEFF );
214 :
215 11 : *pdfPixel = dfResultX * psRPC->dfSAMP_SCALE + psRPC->dfSAMP_OFF;
216 11 : *pdfLine = dfResultY * psRPC->dfLINE_SCALE + psRPC->dfLINE_OFF;
217 11 : }
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 : } GDALRPCTransformInfo;
240 :
241 : /************************************************************************/
242 : /* GDALCreateRPCTransformer() */
243 : /************************************************************************/
244 :
245 : /**
246 : * Create an RPC based transformer.
247 : *
248 : * The geometric sensor model describing the physical relationship between
249 : * image coordinates and ground coordinate is known as a Rigorous Projection
250 : * Model. A Rigorous Projection Model expresses the mapping of the image space
251 : * coordinates of rows and columns (r,c) onto the object space reference
252 : * surface geodetic coordinates (long, lat, height).
253 : *
254 : * RPC supports a generic description of the Rigorous Projection Models. The
255 : * approximation used by GDAL (RPC00) is a set of rational polynomials exp
256 : * ressing the normalized row and column values, (rn , cn), as a function of
257 : * normalized geodetic latitude, longitude, and height, (P, L, H), given a
258 : * set of normalized polynomial coefficients (LINE_NUM_COEF_n, LINE_DEN_COEF_n,
259 : * SAMP_NUM_COEF_n, SAMP_DEN_COEF_n). Normalized values, rather than actual
260 : * values are used in order to minimize introduction of errors during the
261 : * calculations. The transformation between row and column values (r,c), and
262 : * normalized row and column values (rn, cn), and between the geodetic
263 : * latitude, longitude, and height and normalized geodetic latitude,
264 : * longitude, and height (P, L, H), is defined by a set of normalizing
265 : * translations (offsets) and scales that ensure all values are contained i
266 : * the range -1 to +1.
267 : *
268 : * This function creates a GDALTransformFunc compatible transformer
269 : * for going between image pixel/line and long/lat/height coordinates
270 : * using RPCs. The RPCs are provided in a GDALRPCInfo structure which is
271 : * normally read from metadata using GDALExtractRPCInfo().
272 : *
273 : * GDAL RPC Metadata has the following entries (also described in GDAL RFC 22
274 : * and the GeoTIFF RPC document http://geotiff.maptools.org/rpc_prop.html.
275 : *
276 : * <ul>
277 : * <li>ERR_BIAS: Error - Bias. The RMS bias error in meters per horizontal axis of all points in the image (-1.0 if unknown)
278 : * <li>ERR_RAND: Error - Random. RMS random error in meters per horizontal axis of each point in the image (-1.0 if unknown)
279 : * <li>LINE_OFF: Line Offset
280 : * <li>SAMP_OFF: Sample Offset
281 : * <li>LAT_OFF: Geodetic Latitude Offset
282 : * <li>LONG_OFF: Geodetic Longitude Offset
283 : * <li>HEIGHT_OFF: Geodetic Height Offset
284 : * <li>LINE_SCALE: Line Scale
285 : * <li>SAMP_SCALE: Sample Scale
286 : * <li>LAT_SCALE: Geodetic Latitude Scale
287 : * <li>LONG_SCALE: Geodetic Longitude Scale
288 : * <li>HEIGHT_SCALE: Geodetic Height Scale
289 : * <li>LINE_NUM_COEFF (1-20): Line Numerator Coefficients. Twenty coefficients for the polynomial in the Numerator of the rn equation. (space separated)
290 : * <li>LINE_DEN_COEFF (1-20): Line Denominator Coefficients. Twenty coefficients for the polynomial in the Denominator of the rn equation. (space separated)
291 : * <li>SAMP_NUM_COEFF (1-20): Sample Numerator Coefficients. Twenty coefficients for the polynomial in the Numerator of the cn equation. (space separated)
292 : * <li>SAMP_DEN_COEFF (1-20): Sample Denominator Coefficients. Twenty coefficients for the polynomial in the Denominator of the cn equation. (space separated)
293 : * </ul>
294 : *
295 : * The transformer normally maps from pixel/line/height to long/lat/height space
296 : * as a forward transformation though in RPC terms that would be considered
297 : * an inverse transformation (and is solved by iterative approximation using
298 : * long/lat/height to pixel/line transformations). The default direction can
299 : * be reversed by passing bReversed=TRUE.
300 : *
301 : * The iterative solution of pixel/line
302 : * to lat/long/height is currently run for up to 10 iterations or until
303 : * the apparent error is less than dfPixErrThreshold pixels. Passing zero
304 : * will not avoid all error, but will cause the operation to run for the maximum
305 : * number of iterations.
306 : *
307 : * Additional options to the transformer can be supplied in papszOptions.
308 : * Currently only one option is supported, though in the future more may
309 : * be added, notably an option to extract elevation offsets from a DEM file.
310 : *
311 : * Options:
312 : *
313 : * <ul>
314 : * <li> RPC_HEIGHT: a fixed height offset to be applied to all points passed
315 : * in. In this situation the Z passed into the transformation function is
316 : * assumed to be height above ground, and the RPC_HEIGHT is assumed to be
317 : * an average height above sea level for ground in the target scene.
318 : * </ul>
319 : *
320 : * @param psRPCInfo Definition of the RPC parameters.
321 : *
322 : * @param bReversed If true "forward" transformation will be lat/long to pixel/line instead of the normal pixel/line to lat/long.
323 : *
324 : * @param dfPixErrThreshold the error (measured in pixels) allowed in the
325 : * iterative solution of pixel/line to lat/long computations (the other way
326 : * is always exact given the equations).
327 : *
328 : * @param papszOptions Other transformer options (ie. RPC_HEIGHT=<z>).
329 : *
330 : * @return transformer callback data (deallocate with GDALDestroyTransformer()).
331 : */
332 :
333 1 : void *GDALCreateRPCTransformer( GDALRPCInfo *psRPCInfo, int bReversed,
334 : double dfPixErrThreshold,
335 : char **papszOptions )
336 :
337 : {
338 : GDALRPCTransformInfo *psTransform;
339 :
340 : /* -------------------------------------------------------------------- */
341 : /* Initialize core info. */
342 : /* -------------------------------------------------------------------- */
343 : psTransform = (GDALRPCTransformInfo *)
344 1 : CPLCalloc(sizeof(GDALRPCTransformInfo),1);
345 :
346 1 : memcpy( &(psTransform->sRPC), psRPCInfo, sizeof(GDALRPCInfo) );
347 1 : psTransform->bReversed = bReversed;
348 1 : psTransform->dfPixErrThreshold = dfPixErrThreshold;
349 1 : psTransform->dfHeightOffset = 0.0;
350 :
351 1 : strcpy( psTransform->sTI.szSignature, "GTI" );
352 1 : psTransform->sTI.pszClassName = "GDALRPCTransformer";
353 1 : psTransform->sTI.pfnTransform = GDALRPCTransform;
354 1 : psTransform->sTI.pfnCleanup = GDALDestroyRPCTransformer;
355 1 : psTransform->sTI.pfnSerialize = GDALSerializeRPCTransformer;
356 :
357 : /* -------------------------------------------------------------------- */
358 : /* Do we have a "average height" that we want to consider all */
359 : /* elevations to be relative to? */
360 : /* -------------------------------------------------------------------- */
361 1 : const char *pszHeight = CSLFetchNameValue( papszOptions, "RPC_HEIGHT" );
362 1 : if( pszHeight != NULL )
363 0 : psTransform->dfHeightOffset = CPLAtof(pszHeight);
364 :
365 : /* -------------------------------------------------------------------- */
366 : /* Establish a reference point for calcualating an affine */
367 : /* geotransform approximate transformation. */
368 : /* -------------------------------------------------------------------- */
369 1 : double adfGTFromLL[6], dfRefPixel = -1.0, dfRefLine = -1.0;
370 1 : double dfRefLong = 0.0, dfRefLat = 0.0;
371 :
372 1 : if( psRPCInfo->dfMIN_LONG != -180 || psRPCInfo->dfMAX_LONG != 180 )
373 : {
374 0 : dfRefLong = (psRPCInfo->dfMIN_LONG + psRPCInfo->dfMAX_LONG) * 0.5;
375 0 : dfRefLat = (psRPCInfo->dfMIN_LAT + psRPCInfo->dfMAX_LAT ) * 0.5;
376 :
377 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat, 0.0,
378 0 : &dfRefPixel, &dfRefLine );
379 : }
380 :
381 : // Try with scale and offset if we don't can't use bounds or
382 : // the results seem daft.
383 1 : if( dfRefPixel < 0.0 || dfRefLine < 0.0
384 : || dfRefPixel > 100000 || dfRefLine > 100000 )
385 : {
386 1 : dfRefLong = psRPCInfo->dfLONG_OFF;
387 1 : dfRefLat = psRPCInfo->dfLAT_OFF;
388 :
389 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat, 0.0,
390 1 : &dfRefPixel, &dfRefLine );
391 : }
392 :
393 : /* -------------------------------------------------------------------- */
394 : /* Transform nearby locations to establish affine direction */
395 : /* vectors. */
396 : /* -------------------------------------------------------------------- */
397 1 : double dfRefPixelDelta, dfRefLineDelta, dfLLDelta = 0.0001;
398 :
399 : RPCTransformPoint( psRPCInfo, dfRefLong+dfLLDelta, dfRefLat, 0.0,
400 1 : &dfRefPixelDelta, &dfRefLineDelta );
401 1 : adfGTFromLL[1] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
402 1 : adfGTFromLL[2] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
403 :
404 : RPCTransformPoint( psRPCInfo, dfRefLong, dfRefLat+dfLLDelta, 0.0,
405 1 : &dfRefPixelDelta, &dfRefLineDelta );
406 1 : adfGTFromLL[4] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
407 1 : adfGTFromLL[5] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
408 :
409 : adfGTFromLL[0] = dfRefPixel
410 1 : - adfGTFromLL[1] * dfRefLong - adfGTFromLL[2] * dfRefLat;
411 : adfGTFromLL[3] = dfRefLine
412 1 : - adfGTFromLL[4] * dfRefLong - adfGTFromLL[5] * dfRefLat;
413 :
414 1 : GDALInvGeoTransform( adfGTFromLL, psTransform->adfPLToLatLongGeoTransform);
415 :
416 1 : return psTransform;
417 : }
418 :
419 : /************************************************************************/
420 : /* GDALDestroyReprojectionTransformer() */
421 : /************************************************************************/
422 :
423 1 : void GDALDestroyRPCTransformer( void *pTransformAlg )
424 :
425 : {
426 1 : CPLFree( pTransformAlg );
427 1 : }
428 :
429 : /************************************************************************/
430 : /* RPCInverseTransformPoint() */
431 : /************************************************************************/
432 :
433 : static void
434 2 : RPCInverseTransformPoint( GDALRPCTransformInfo *psTransform,
435 : double dfPixel, double dfLine, double dfHeight,
436 : double *pdfLong, double *pdfLat )
437 :
438 : {
439 : double dfResultX, dfResultY;
440 : int iIter;
441 2 : GDALRPCInfo *psRPC = &(psTransform->sRPC);
442 :
443 : /* -------------------------------------------------------------------- */
444 : /* Compute an initial approximation based on linear */
445 : /* interpolation from our reference point. */
446 : /* -------------------------------------------------------------------- */
447 2 : dfResultX = psTransform->adfPLToLatLongGeoTransform[0]
448 2 : + psTransform->adfPLToLatLongGeoTransform[1] * dfPixel
449 2 : + psTransform->adfPLToLatLongGeoTransform[2] * dfLine;
450 :
451 2 : dfResultY = psTransform->adfPLToLatLongGeoTransform[3]
452 2 : + psTransform->adfPLToLatLongGeoTransform[4] * dfPixel
453 2 : + psTransform->adfPLToLatLongGeoTransform[5] * dfLine;
454 :
455 : /* -------------------------------------------------------------------- */
456 : /* Now iterate, trying to find a closer LL location that will */
457 : /* back transform to the indicated pixel and line. */
458 : /* -------------------------------------------------------------------- */
459 : double dfPixelDeltaX, dfPixelDeltaY;
460 :
461 6 : for( iIter = 0; iIter < 10; iIter++ )
462 : {
463 : double dfBackPixel, dfBackLine;
464 :
465 : RPCTransformPoint( psRPC, dfResultX, dfResultY, dfHeight,
466 6 : &dfBackPixel, &dfBackLine );
467 :
468 6 : dfPixelDeltaX = dfBackPixel - dfPixel;
469 6 : dfPixelDeltaY = dfBackLine - dfLine;
470 :
471 : dfResultX = dfResultX
472 6 : - dfPixelDeltaX * psTransform->adfPLToLatLongGeoTransform[1]
473 6 : - dfPixelDeltaY * psTransform->adfPLToLatLongGeoTransform[2];
474 : dfResultY = dfResultY
475 6 : - dfPixelDeltaX * psTransform->adfPLToLatLongGeoTransform[4]
476 6 : - dfPixelDeltaY * psTransform->adfPLToLatLongGeoTransform[5];
477 :
478 6 : if( ABS(dfPixelDeltaX) < psTransform->dfPixErrThreshold
479 : && ABS(dfPixelDeltaY) < psTransform->dfPixErrThreshold )
480 : {
481 2 : iIter = -1;
482 : //CPLDebug( "RPC", "Converged!" );
483 2 : break;
484 : }
485 :
486 : }
487 :
488 2 : if( iIter != -1 )
489 : CPLDebug( "RPC", "Iterations %d: Got: %g,%g Offset=%g,%g",
490 : iIter,
491 : dfResultX, dfResultY,
492 0 : dfPixelDeltaX, dfPixelDeltaY );
493 :
494 2 : *pdfLong = dfResultX;
495 2 : *pdfLat = dfResultY;
496 2 : }
497 :
498 : /************************************************************************/
499 : /* GDALRPCTransform() */
500 : /************************************************************************/
501 :
502 4 : int GDALRPCTransform( void *pTransformArg, int bDstToSrc,
503 : int nPointCount,
504 : double *padfX, double *padfY, double *padfZ,
505 : int *panSuccess )
506 :
507 : {
508 4 : VALIDATE_POINTER1( pTransformArg, "GDALRPCTransform", 0 );
509 :
510 4 : GDALRPCTransformInfo *psTransform = (GDALRPCTransformInfo *) pTransformArg;
511 4 : GDALRPCInfo *psRPC = &(psTransform->sRPC);
512 : int i;
513 :
514 4 : if( psTransform->bReversed )
515 0 : bDstToSrc = !bDstToSrc;
516 :
517 : /* -------------------------------------------------------------------- */
518 : /* The simple case is transforming from lat/long to pixel/line. */
519 : /* Just apply the equations directly. */
520 : /* -------------------------------------------------------------------- */
521 4 : if( bDstToSrc )
522 : {
523 4 : for( i = 0; i < nPointCount; i++ )
524 : {
525 : RPCTransformPoint( psRPC, padfX[i], padfY[i],
526 2 : padfZ[i] + psTransform->dfHeightOffset,
527 4 : padfX + i, padfY + i );
528 2 : panSuccess[i] = TRUE;
529 : }
530 :
531 2 : return TRUE;
532 : }
533 :
534 : /* -------------------------------------------------------------------- */
535 : /* Compute the inverse (pixel/line/height to lat/long). This */
536 : /* function uses an iterative method from an initial linear */
537 : /* approximation. */
538 : /* -------------------------------------------------------------------- */
539 4 : for( i = 0; i < nPointCount; i++ )
540 : {
541 : double dfResultX, dfResultY;
542 :
543 : RPCInverseTransformPoint( psTransform, padfX[i], padfY[i],
544 2 : padfZ[i] + psTransform->dfHeightOffset,
545 2 : &dfResultX, &dfResultY );
546 :
547 2 : padfX[i] = dfResultX;
548 2 : padfY[i] = dfResultY;
549 :
550 2 : panSuccess[i] = TRUE;
551 : }
552 :
553 2 : return TRUE;
554 : }
555 :
556 : /************************************************************************/
557 : /* GDALSerializeRPCTransformer() */
558 : /************************************************************************/
559 :
560 0 : CPLXMLNode *GDALSerializeRPCTransformer( void *pTransformArg )
561 :
562 : {
563 0 : VALIDATE_POINTER1( pTransformArg, "GDALSerializeRPCTransformer", NULL );
564 :
565 : CPLXMLNode *psTree;
566 : GDALRPCTransformInfo *psInfo =
567 0 : (GDALRPCTransformInfo *)(pTransformArg);
568 :
569 0 : psTree = CPLCreateXMLNode( NULL, CXT_Element, "RPCTransformer" );
570 :
571 : /* -------------------------------------------------------------------- */
572 : /* Serialize bReversed. */
573 : /* -------------------------------------------------------------------- */
574 : CPLCreateXMLElementAndValue(
575 : psTree, "Reversed",
576 0 : CPLString().Printf( "%d", psInfo->bReversed ) );
577 :
578 : /* -------------------------------------------------------------------- */
579 : /* Serialize Height Offset. */
580 : /* -------------------------------------------------------------------- */
581 : CPLCreateXMLElementAndValue(
582 : psTree, "HeightOffset",
583 0 : CPLString().Printf( "%.15g", psInfo->dfHeightOffset ) );
584 :
585 : /* -------------------------------------------------------------------- */
586 : /* Serialize pixel error threshold. */
587 : /* -------------------------------------------------------------------- */
588 : CPLCreateXMLElementAndValue(
589 : psTree, "PixErrThreshold",
590 0 : CPLString().Printf( "%.15g", psInfo->dfPixErrThreshold ) );
591 :
592 : /* -------------------------------------------------------------------- */
593 : /* RPC metadata. */
594 : /* -------------------------------------------------------------------- */
595 0 : char **papszMD = RPCInfoToMD( &(psInfo->sRPC) );
596 : CPLXMLNode *psMD= CPLCreateXMLNode( psTree, CXT_Element,
597 0 : "Metadata" );
598 :
599 0 : for( int i = 0; papszMD != NULL && papszMD[i] != NULL; i++ )
600 : {
601 : const char *pszRawValue;
602 : char *pszKey;
603 : CPLXMLNode *psMDI;
604 :
605 0 : pszRawValue = CPLParseNameValue( papszMD[i], &pszKey );
606 :
607 0 : psMDI = CPLCreateXMLNode( psMD, CXT_Element, "MDI" );
608 0 : CPLSetXMLValue( psMDI, "#key", pszKey );
609 0 : CPLCreateXMLNode( psMDI, CXT_Text, pszRawValue );
610 :
611 0 : CPLFree( pszKey );
612 : }
613 :
614 0 : CSLDestroy( papszMD );
615 :
616 0 : return psTree;
617 : }
618 :
619 : /************************************************************************/
620 : /* GDALDeserializeRPCTransformer() */
621 : /************************************************************************/
622 :
623 0 : void *GDALDeserializeRPCTransformer( CPLXMLNode *psTree )
624 :
625 : {
626 : void *pResult;
627 0 : char **papszOptions = NULL;
628 :
629 : /* -------------------------------------------------------------------- */
630 : /* Collect metadata. */
631 : /* -------------------------------------------------------------------- */
632 0 : char **papszMD = NULL;
633 : CPLXMLNode *psMDI, *psMetadata;
634 : GDALRPCInfo sRPC;
635 :
636 0 : psMetadata = CPLGetXMLNode( psTree, "Metadata" );
637 :
638 0 : if( psMetadata->eType != CXT_Element
639 : || !EQUAL(psMetadata->pszValue,"Metadata") )
640 0 : return NULL;
641 :
642 0 : for( psMDI = psMetadata->psChild; psMDI != NULL;
643 : psMDI = psMDI->psNext )
644 : {
645 0 : if( !EQUAL(psMDI->pszValue,"MDI")
646 : || psMDI->eType != CXT_Element
647 : || psMDI->psChild == NULL
648 : || psMDI->psChild->psNext == NULL
649 : || psMDI->psChild->eType != CXT_Attribute
650 : || psMDI->psChild->psChild == NULL )
651 0 : continue;
652 :
653 : papszMD =
654 : CSLSetNameValue( papszMD,
655 : psMDI->psChild->psChild->pszValue,
656 0 : psMDI->psChild->psNext->pszValue );
657 : }
658 :
659 0 : if( !GDALExtractRPCInfo( papszMD, &sRPC ) )
660 : {
661 : CPLError( CE_Failure, CPLE_AppDefined,
662 0 : "Failed to reconstitute RPC transformer." );
663 0 : return NULL;
664 : }
665 :
666 0 : CSLDestroy( papszMD );
667 :
668 : /* -------------------------------------------------------------------- */
669 : /* Get other flags. */
670 : /* -------------------------------------------------------------------- */
671 : double dfPixErrThreshold;
672 : int bReversed;
673 :
674 0 : bReversed = atoi(CPLGetXMLValue(psTree,"Reversed","0"));
675 :
676 : dfPixErrThreshold =
677 0 : CPLAtof(CPLGetXMLValue(psTree,"PixErrThreshold","0.25"));
678 :
679 : papszOptions = CSLSetNameValue( papszOptions, "RPC_HEIGHT",
680 0 : CPLGetXMLValue(psTree,"HeightOffset","0"));
681 :
682 : /* -------------------------------------------------------------------- */
683 : /* Generate transformation. */
684 : /* -------------------------------------------------------------------- */
685 : pResult = GDALCreateRPCTransformer( &sRPC, bReversed, dfPixErrThreshold,
686 0 : papszOptions );
687 :
688 0 : CSLDestroy( papszOptions );
689 :
690 0 : return pResult;
691 : }
|
