Line data Source code
1 : /******************************************************************************
2 : *
3 : * Project: Image Warper
4 : * Purpose: Implements a rational polynomial (RPC) based transformer.
5 : * Author: Frank Warmerdam, warmerdam@pobox.com
6 : *
7 : ******************************************************************************
8 : * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>
9 : * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>
10 : *
11 : * SPDX-License-Identifier: MIT
12 : ****************************************************************************/
13 :
14 : #include "cpl_port.h"
15 : #include "gdal_alg.h"
16 :
17 : #include <cmath>
18 : #include <cstddef>
19 : #include <cstdlib>
20 : #include <cstring>
21 :
22 : #include <algorithm>
23 : #include <limits>
24 : #include <string>
25 :
26 : #include "cpl_conv.h"
27 : #include "cpl_error.h"
28 : #include "cpl_mem_cache.h"
29 : #include "cpl_minixml.h"
30 : #include "cpl_string.h"
31 : #include "cpl_vsi.h"
32 : #include "gdal.h"
33 : #include "gdal_interpolateatpoint.h"
34 : #include "gdal_mdreader.h"
35 : #include "gdal_alg_priv.h"
36 : #include "gdal_priv.h"
37 :
38 : #ifdef USE_NEON_OPTIMIZATIONS
39 : #define USE_SSE2
40 : #elif defined(__x86_64) || defined(_M_X64)
41 : #define USE_SSE2
42 : #endif
43 :
44 : #ifdef USE_SSE2
45 : #include "gdalsse_priv.h"
46 : #define USE_SSE2_OPTIM
47 : #endif
48 :
49 : #include "ogr_api.h"
50 : #include "ogr_geometry.h"
51 : #include "ogr_spatialref.h"
52 : #include "ogr_srs_api.h"
53 : #include "gdalresamplingkernels.h"
54 :
55 : // #define DEBUG_VERBOSE_EXTRACT_DEM
56 :
57 : CPL_C_START
58 : CPLXMLNode *GDALSerializeRPCTransformer(void *pTransformArg);
59 : void *GDALDeserializeRPCTransformer(CPLXMLNode *psTree);
60 : CPL_C_END
61 :
62 : constexpr int MAX_ABS_VALUE_WARNINGS = 20;
63 : constexpr double DEFAULT_PIX_ERR_THRESHOLD = 0.1;
64 :
65 : /************************************************************************/
66 : /* RPCInfoToMD() */
67 : /* */
68 : /* Turn an RPCInfo structure back into its metadata format. */
69 : /************************************************************************/
70 :
71 0 : char **RPCInfoV1ToMD(GDALRPCInfoV1 *psRPCInfo)
72 :
73 : {
74 : GDALRPCInfoV2 sRPCInfo;
75 0 : memcpy(&sRPCInfo, psRPCInfo, sizeof(GDALRPCInfoV1));
76 0 : sRPCInfo.dfERR_BIAS = std::numeric_limits<double>::quiet_NaN();
77 0 : sRPCInfo.dfERR_RAND = std::numeric_limits<double>::quiet_NaN();
78 0 : return RPCInfoV2ToMD(&sRPCInfo);
79 : }
80 :
81 1 : char **RPCInfoV2ToMD(GDALRPCInfoV2 *psRPCInfo)
82 :
83 : {
84 1 : char **papszMD = nullptr;
85 2 : CPLString osField, osMultiField;
86 :
87 1 : if (!std::isnan(psRPCInfo->dfERR_BIAS))
88 : {
89 1 : osField.Printf("%.15g", psRPCInfo->dfERR_BIAS);
90 1 : papszMD = CSLSetNameValue(papszMD, RPC_ERR_BIAS, osField);
91 : }
92 :
93 1 : if (!std::isnan(psRPCInfo->dfERR_RAND))
94 : {
95 1 : osField.Printf("%.15g", psRPCInfo->dfERR_RAND);
96 1 : papszMD = CSLSetNameValue(papszMD, RPC_ERR_RAND, osField);
97 : }
98 :
99 1 : osField.Printf("%.15g", psRPCInfo->dfLINE_OFF);
100 1 : papszMD = CSLSetNameValue(papszMD, RPC_LINE_OFF, osField);
101 :
102 1 : osField.Printf("%.15g", psRPCInfo->dfSAMP_OFF);
103 1 : papszMD = CSLSetNameValue(papszMD, RPC_SAMP_OFF, osField);
104 :
105 1 : osField.Printf("%.15g", psRPCInfo->dfLAT_OFF);
106 1 : papszMD = CSLSetNameValue(papszMD, RPC_LAT_OFF, osField);
107 :
108 1 : osField.Printf("%.15g", psRPCInfo->dfLONG_OFF);
109 1 : papszMD = CSLSetNameValue(papszMD, RPC_LONG_OFF, osField);
110 :
111 1 : osField.Printf("%.15g", psRPCInfo->dfHEIGHT_OFF);
112 1 : papszMD = CSLSetNameValue(papszMD, RPC_HEIGHT_OFF, osField);
113 :
114 1 : osField.Printf("%.15g", psRPCInfo->dfLINE_SCALE);
115 1 : papszMD = CSLSetNameValue(papszMD, RPC_LINE_SCALE, osField);
116 :
117 1 : osField.Printf("%.15g", psRPCInfo->dfSAMP_SCALE);
118 1 : papszMD = CSLSetNameValue(papszMD, RPC_SAMP_SCALE, osField);
119 :
120 1 : osField.Printf("%.15g", psRPCInfo->dfLAT_SCALE);
121 1 : papszMD = CSLSetNameValue(papszMD, RPC_LAT_SCALE, osField);
122 :
123 1 : osField.Printf("%.15g", psRPCInfo->dfLONG_SCALE);
124 1 : papszMD = CSLSetNameValue(papszMD, RPC_LONG_SCALE, osField);
125 :
126 1 : osField.Printf("%.15g", psRPCInfo->dfHEIGHT_SCALE);
127 1 : papszMD = CSLSetNameValue(papszMD, RPC_HEIGHT_SCALE, osField);
128 :
129 1 : osField.Printf("%.15g", psRPCInfo->dfMIN_LONG);
130 1 : papszMD = CSLSetNameValue(papszMD, RPC_MIN_LONG, osField);
131 :
132 1 : osField.Printf("%.15g", psRPCInfo->dfMIN_LAT);
133 1 : papszMD = CSLSetNameValue(papszMD, RPC_MIN_LAT, osField);
134 :
135 1 : osField.Printf("%.15g", psRPCInfo->dfMAX_LONG);
136 1 : papszMD = CSLSetNameValue(papszMD, RPC_MAX_LONG, osField);
137 :
138 1 : osField.Printf("%.15g", psRPCInfo->dfMAX_LAT);
139 1 : papszMD = CSLSetNameValue(papszMD, RPC_MAX_LAT, osField);
140 :
141 21 : for (int i = 0; i < 20; i++)
142 : {
143 20 : osField.Printf("%.15g", psRPCInfo->adfLINE_NUM_COEFF[i]);
144 20 : if (i > 0)
145 19 : osMultiField += " ";
146 : else
147 1 : osMultiField = "";
148 20 : osMultiField += osField;
149 : }
150 1 : papszMD = CSLSetNameValue(papszMD, "LINE_NUM_COEFF", osMultiField);
151 :
152 21 : for (int i = 0; i < 20; i++)
153 : {
154 20 : osField.Printf("%.15g", psRPCInfo->adfLINE_DEN_COEFF[i]);
155 20 : if (i > 0)
156 19 : osMultiField += " ";
157 : else
158 1 : osMultiField = "";
159 20 : osMultiField += osField;
160 : }
161 1 : papszMD = CSLSetNameValue(papszMD, "LINE_DEN_COEFF", osMultiField);
162 :
163 21 : for (int i = 0; i < 20; i++)
164 : {
165 20 : osField.Printf("%.15g", psRPCInfo->adfSAMP_NUM_COEFF[i]);
166 20 : if (i > 0)
167 19 : osMultiField += " ";
168 : else
169 1 : osMultiField = "";
170 20 : osMultiField += osField;
171 : }
172 1 : papszMD = CSLSetNameValue(papszMD, "SAMP_NUM_COEFF", osMultiField);
173 :
174 21 : for (int i = 0; i < 20; i++)
175 : {
176 20 : osField.Printf("%.15g", psRPCInfo->adfSAMP_DEN_COEFF[i]);
177 20 : if (i > 0)
178 19 : osMultiField += " ";
179 : else
180 1 : osMultiField = "";
181 20 : osMultiField += osField;
182 : }
183 1 : papszMD = CSLSetNameValue(papszMD, "SAMP_DEN_COEFF", osMultiField);
184 :
185 2 : return papszMD;
186 : }
187 :
188 : /************************************************************************/
189 : /* RPCComputeTerms() */
190 : /************************************************************************/
191 :
192 460269 : static void RPCComputeTerms(double dfLong, double dfLat, double dfHeight,
193 : double *padfTerms)
194 :
195 : {
196 460269 : padfTerms[0] = 1.0;
197 460269 : padfTerms[1] = dfLong;
198 460269 : padfTerms[2] = dfLat;
199 460269 : padfTerms[3] = dfHeight;
200 460269 : padfTerms[4] = dfLong * dfLat;
201 460269 : padfTerms[5] = dfLong * dfHeight;
202 460269 : padfTerms[6] = dfLat * dfHeight;
203 460269 : padfTerms[7] = dfLong * dfLong;
204 460269 : padfTerms[8] = dfLat * dfLat;
205 460269 : padfTerms[9] = dfHeight * dfHeight;
206 :
207 460269 : padfTerms[10] = dfLong * dfLat * dfHeight;
208 460269 : padfTerms[11] = dfLong * dfLong * dfLong;
209 460269 : padfTerms[12] = dfLong * dfLat * dfLat;
210 460269 : padfTerms[13] = dfLong * dfHeight * dfHeight;
211 460269 : padfTerms[14] = dfLong * dfLong * dfLat;
212 460269 : padfTerms[15] = dfLat * dfLat * dfLat;
213 460269 : padfTerms[16] = dfLat * dfHeight * dfHeight;
214 460269 : padfTerms[17] = dfLong * dfLong * dfHeight;
215 460269 : padfTerms[18] = dfLat * dfLat * dfHeight;
216 460269 : padfTerms[19] = dfHeight * dfHeight * dfHeight;
217 460269 : }
218 :
219 : /************************************************************************/
220 : /* ==================================================================== */
221 : /* GDALRPCTransformer */
222 : /* ==================================================================== */
223 : /************************************************************************/
224 :
225 : /*! DEM Resampling Algorithm */
226 : typedef enum
227 : {
228 : /*! Nearest neighbour (select on one input pixel) */ DRA_NearestNeighbour =
229 : 0,
230 : /*! Bilinear (2x2 kernel) */ DRA_Bilinear = 1,
231 : /*! Cubic Convolution Approximation (4x4 kernel) */ DRA_CubicSpline = 2
232 : } DEMResampleAlg;
233 :
234 : typedef struct
235 : {
236 :
237 : GDALTransformerInfo sTI;
238 :
239 : GDALRPCInfoV2 sRPC;
240 :
241 : double adfPLToLatLongGeoTransform[6];
242 : double dfRefZ;
243 :
244 : int bReversed;
245 :
246 : double dfPixErrThreshold;
247 :
248 : double dfHeightOffset;
249 :
250 : double dfHeightScale;
251 :
252 : char *pszDEMPath;
253 :
254 : DEMResampleAlg eResampleAlg;
255 :
256 : int bHasDEMMissingValue;
257 : double dfDEMMissingValue;
258 : char *pszDEMSRS;
259 : int bApplyDEMVDatumShift;
260 :
261 : GDALDataset *poDS;
262 : // the key is (nYBlock << 32) | nXBlock)
263 : lru11::Cache<uint64_t, std::shared_ptr<std::vector<double>>> *poCacheDEM;
264 :
265 : OGRCoordinateTransformation *poCT;
266 :
267 : int nMaxIterations;
268 :
269 : double adfDEMGeoTransform[6];
270 : double adfDEMReverseGeoTransform[6];
271 :
272 : #ifdef USE_SSE2_OPTIM
273 : double adfDoubles[20 * 4 + 1];
274 : // LINE_NUM_COEFF, LINE_DEN_COEFF, SAMP_NUM_COEFF and then SAMP_DEN_COEFF.
275 : double *padfCoeffs;
276 : #endif
277 :
278 : bool bRPCInverseVerbose;
279 : char *pszRPCInverseLog;
280 :
281 : char *pszRPCFootprint;
282 : OGRGeometry *poRPCFootprintGeom;
283 : OGRPreparedGeometry *poRPCFootprintPreparedGeom;
284 :
285 : } GDALRPCTransformInfo;
286 :
287 : static bool GDALRPCOpenDEM(GDALRPCTransformInfo *psTransform);
288 :
289 : /************************************************************************/
290 : /* RPCEvaluate() */
291 : /************************************************************************/
292 : #ifdef USE_SSE2_OPTIM
293 :
294 460269 : static void RPCEvaluate4(const double *padfTerms, const double *padfCoefs,
295 : double &dfSum1, double &dfSum2, double &dfSum3,
296 : double &dfSum4)
297 :
298 : {
299 460269 : XMMReg2Double sum1 = XMMReg2Double::Zero();
300 460269 : XMMReg2Double sum2 = XMMReg2Double::Zero();
301 460269 : XMMReg2Double sum3 = XMMReg2Double::Zero();
302 460269 : XMMReg2Double sum4 = XMMReg2Double::Zero();
303 5062960 : for (int i = 0; i < 20; i += 2)
304 : {
305 : const XMMReg2Double terms =
306 4602690 : XMMReg2Double::Load2ValAligned(padfTerms + i);
307 :
308 : // LINE_NUM_COEFF.
309 : const XMMReg2Double coefs1 =
310 4602690 : XMMReg2Double::Load2ValAligned(padfCoefs + i);
311 :
312 : // LINE_DEN_COEFF.
313 : const XMMReg2Double coefs2 =
314 4602690 : XMMReg2Double::Load2ValAligned(padfCoefs + i + 20);
315 :
316 : // SAMP_NUM_COEFF.
317 : const XMMReg2Double coefs3 =
318 4602690 : XMMReg2Double::Load2ValAligned(padfCoefs + i + 40);
319 :
320 : // SAMP_DEN_COEFF.
321 : const XMMReg2Double coefs4 =
322 4602690 : XMMReg2Double::Load2ValAligned(padfCoefs + i + 60);
323 :
324 4602690 : sum1 += terms * coefs1;
325 4602690 : sum2 += terms * coefs2;
326 4602690 : sum3 += terms * coefs3;
327 4602690 : sum4 += terms * coefs4;
328 : }
329 460269 : dfSum1 = sum1.GetHorizSum();
330 460269 : dfSum2 = sum2.GetHorizSum();
331 460269 : dfSum3 = sum3.GetHorizSum();
332 460269 : dfSum4 = sum4.GetHorizSum();
333 460269 : }
334 :
335 : #else
336 :
337 : static double RPCEvaluate(const double *padfTerms, const double *padfCoefs)
338 :
339 : {
340 : double dfSum1 = 0.0;
341 : double dfSum2 = 0.0;
342 :
343 : for (int i = 0; i < 20; i += 2)
344 : {
345 : dfSum1 += padfTerms[i] * padfCoefs[i];
346 : dfSum2 += padfTerms[i + 1] * padfCoefs[i + 1];
347 : }
348 :
349 : return dfSum1 + dfSum2;
350 : }
351 :
352 : #endif
353 :
354 : /************************************************************************/
355 : /* RPCTransformPoint() */
356 : /************************************************************************/
357 :
358 460269 : static void RPCTransformPoint(const GDALRPCTransformInfo *psRPCTransformInfo,
359 : double dfLong, double dfLat, double dfHeight,
360 : double *pdfPixel, double *pdfLine)
361 :
362 : {
363 460269 : double adfTermsWithMargin[20 + 1] = {};
364 : // Make padfTerms aligned on 16-byte boundary for SSE2 aligned loads.
365 460269 : double *padfTerms =
366 460269 : adfTermsWithMargin +
367 : (reinterpret_cast<GUIntptr_t>(adfTermsWithMargin) % 16) / 8;
368 :
369 : // Avoid dateline issues.
370 460269 : double diffLong = dfLong - psRPCTransformInfo->sRPC.dfLONG_OFF;
371 460269 : if (diffLong < -270)
372 : {
373 2 : diffLong += 360;
374 : }
375 460267 : else if (diffLong > 270)
376 : {
377 6 : diffLong -= 360;
378 : }
379 :
380 460269 : const double dfNormalizedLong =
381 460269 : diffLong / psRPCTransformInfo->sRPC.dfLONG_SCALE;
382 460269 : const double dfNormalizedLat =
383 460269 : (dfLat - psRPCTransformInfo->sRPC.dfLAT_OFF) /
384 460269 : psRPCTransformInfo->sRPC.dfLAT_SCALE;
385 460269 : const double dfNormalizedHeight =
386 460269 : (dfHeight - psRPCTransformInfo->sRPC.dfHEIGHT_OFF) /
387 460269 : psRPCTransformInfo->sRPC.dfHEIGHT_SCALE;
388 :
389 : // The absolute values of the 3 above normalized values are supposed to be
390 : // below 1. Warn (as debug message) if it is not the case. We allow for some
391 : // margin above 1 (1.5, somewhat arbitrary chosen) before warning.
392 : static int nCountWarningsAboutAboveOneNormalizedValues = 0;
393 460269 : if (nCountWarningsAboutAboveOneNormalizedValues < MAX_ABS_VALUE_WARNINGS)
394 : {
395 99804 : bool bWarned = false;
396 99804 : if (fabs(dfNormalizedLong) > 1.5)
397 : {
398 31 : bWarned = true;
399 31 : CPLDebug(
400 : "RPC",
401 : "Normalized %s for (lon,lat,height)=(%f,%f,%f) is %f, "
402 : "i.e. with an absolute value of > 1, which may cause numeric "
403 : "stability problems",
404 : "longitude", dfLong, dfLat, dfHeight, dfNormalizedLong);
405 : }
406 99804 : if (fabs(dfNormalizedLat) > 1.5)
407 : {
408 21 : bWarned = true;
409 21 : CPLDebug(
410 : "RPC",
411 : "Normalized %s for (lon,lat,height)=(%f,%f,%f) is %f, "
412 : "ie with an absolute value of > 1, which may cause numeric "
413 : "stability problems",
414 : "latitude", dfLong, dfLat, dfHeight, dfNormalizedLat);
415 : }
416 99804 : if (fabs(dfNormalizedHeight) > 1.5)
417 : {
418 39 : bWarned = true;
419 39 : CPLDebug(
420 : "RPC",
421 : "Normalized %s for (lon,lat,height)=(%f,%f,%f) is %f, "
422 : "i.e. with an absolute value of > 1, which may cause numeric "
423 : "stability problems",
424 : "height", dfLong, dfLat, dfHeight, dfNormalizedHeight);
425 : }
426 99804 : if (bWarned)
427 : {
428 : // Limit the number of warnings.
429 60 : nCountWarningsAboutAboveOneNormalizedValues++;
430 60 : if (nCountWarningsAboutAboveOneNormalizedValues ==
431 : MAX_ABS_VALUE_WARNINGS)
432 : {
433 3 : CPLDebug("RPC", "No more such debug warnings will be emitted");
434 : }
435 : }
436 : }
437 :
438 460269 : RPCComputeTerms(dfNormalizedLong, dfNormalizedLat, dfNormalizedHeight,
439 : padfTerms);
440 :
441 : #ifdef USE_SSE2_OPTIM
442 460269 : double dfSampNum = 0.0;
443 460269 : double dfSampDen = 0.0;
444 460269 : double dfLineNum = 0.0;
445 460269 : double dfLineDen = 0.0;
446 460269 : RPCEvaluate4(padfTerms, psRPCTransformInfo->padfCoeffs, dfLineNum,
447 : dfLineDen, dfSampNum, dfSampDen);
448 460269 : const double dfResultX = dfSampNum / dfSampDen;
449 460269 : const double dfResultY = dfLineNum / dfLineDen;
450 : #else
451 : const double dfResultX =
452 : RPCEvaluate(padfTerms, psRPCTransformInfo->sRPC.adfSAMP_NUM_COEFF) /
453 : RPCEvaluate(padfTerms, psRPCTransformInfo->sRPC.adfSAMP_DEN_COEFF);
454 :
455 : const double dfResultY =
456 : RPCEvaluate(padfTerms, psRPCTransformInfo->sRPC.adfLINE_NUM_COEFF) /
457 : RPCEvaluate(padfTerms, psRPCTransformInfo->sRPC.adfLINE_DEN_COEFF);
458 : #endif
459 :
460 : // RPCs are using the center of upper left pixel = 0,0 convention
461 : // convert to top left corner = 0,0 convention used in GDAL.
462 460269 : *pdfPixel = dfResultX * psRPCTransformInfo->sRPC.dfSAMP_SCALE +
463 460269 : psRPCTransformInfo->sRPC.dfSAMP_OFF + 0.5;
464 460269 : *pdfLine = dfResultY * psRPCTransformInfo->sRPC.dfLINE_SCALE +
465 460269 : psRPCTransformInfo->sRPC.dfLINE_OFF + 0.5;
466 460269 : }
467 :
468 : /************************************************************************/
469 : /* GDALSerializeRPCDEMResample() */
470 : /************************************************************************/
471 :
472 0 : static const char *GDALSerializeRPCDEMResample(DEMResampleAlg eResampleAlg)
473 : {
474 0 : switch (eResampleAlg)
475 : {
476 0 : case DRA_NearestNeighbour:
477 0 : return "near";
478 0 : case DRA_CubicSpline:
479 0 : return "cubic";
480 0 : default:
481 : case DRA_Bilinear:
482 0 : return "bilinear";
483 : }
484 : }
485 :
486 : /************************************************************************/
487 : /* GDALCreateSimilarRPCTransformer() */
488 : /************************************************************************/
489 :
490 1 : static void *GDALCreateSimilarRPCTransformer(void *hTransformArg,
491 : double dfRatioX, double dfRatioY)
492 : {
493 1 : VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarRPCTransformer",
494 : nullptr);
495 :
496 1 : GDALRPCTransformInfo *psInfo =
497 : static_cast<GDALRPCTransformInfo *>(hTransformArg);
498 :
499 : GDALRPCInfoV2 sRPC;
500 1 : memcpy(&sRPC, &(psInfo->sRPC), sizeof(GDALRPCInfoV2));
501 :
502 1 : if (dfRatioX != 1.0 || dfRatioY != 1.0)
503 : {
504 1 : sRPC.dfLINE_OFF /= dfRatioY;
505 1 : sRPC.dfLINE_SCALE /= dfRatioY;
506 1 : sRPC.dfSAMP_OFF /= dfRatioX;
507 1 : sRPC.dfSAMP_SCALE /= dfRatioX;
508 : }
509 :
510 1 : char **papszOptions = nullptr;
511 1 : papszOptions = CSLSetNameValue(papszOptions, "RPC_HEIGHT",
512 : CPLSPrintf("%.17g", psInfo->dfHeightOffset));
513 1 : papszOptions = CSLSetNameValue(papszOptions, "RPC_HEIGHT_SCALE",
514 : CPLSPrintf("%.17g", psInfo->dfHeightScale));
515 1 : if (psInfo->pszDEMPath != nullptr)
516 : {
517 : papszOptions =
518 0 : CSLSetNameValue(papszOptions, "RPC_DEM", psInfo->pszDEMPath);
519 : papszOptions =
520 0 : CSLSetNameValue(papszOptions, "RPC_DEMINTERPOLATION",
521 : GDALSerializeRPCDEMResample(psInfo->eResampleAlg));
522 0 : if (psInfo->bHasDEMMissingValue)
523 : papszOptions =
524 0 : CSLSetNameValue(papszOptions, "RPC_DEM_MISSING_VALUE",
525 : CPLSPrintf("%.17g", psInfo->dfDEMMissingValue));
526 : papszOptions =
527 0 : CSLSetNameValue(papszOptions, "RPC_DEM_APPLY_VDATUM_SHIFT",
528 0 : (psInfo->bApplyDEMVDatumShift) ? "TRUE" : "FALSE");
529 : }
530 1 : papszOptions = CSLSetNameValue(papszOptions, "RPC_MAX_ITERATIONS",
531 : CPLSPrintf("%d", psInfo->nMaxIterations));
532 :
533 : GDALRPCTransformInfo *psNewInfo =
534 1 : static_cast<GDALRPCTransformInfo *>(GDALCreateRPCTransformerV2(
535 : &sRPC, psInfo->bReversed, psInfo->dfPixErrThreshold, papszOptions));
536 1 : CSLDestroy(papszOptions);
537 :
538 1 : return psNewInfo;
539 : }
540 :
541 : /************************************************************************/
542 : /* GDALRPCGetHeightAtLongLat() */
543 : /************************************************************************/
544 :
545 : static int GDALRPCGetDEMHeight(GDALRPCTransformInfo *psTransform,
546 : const double dfXIn, const double dfYIn,
547 : double *pdfDEMH);
548 :
549 304118 : static bool GDALRPCGetHeightAtLongLat(GDALRPCTransformInfo *psTransform,
550 : const double dfXIn, const double dfYIn,
551 : double *pdfHeight,
552 : double *pdfDEMPixel = nullptr,
553 : double *pdfDEMLine = nullptr)
554 : {
555 304118 : double dfVDatumShift = 0.0;
556 304118 : double dfDEMH = 0.0;
557 304118 : if (psTransform->poDS)
558 : {
559 296262 : double dfX = 0.0;
560 296262 : double dfY = 0.0;
561 296262 : double dfXTemp = dfXIn;
562 296262 : double dfYTemp = dfYIn;
563 : // Check if dem is not in WGS84 and transform points padfX[i], padfY[i].
564 296262 : if (psTransform->poCT)
565 : {
566 189912 : double dfZ = 0.0;
567 189912 : if (!psTransform->poCT->Transform(1, &dfXTemp, &dfYTemp, &dfZ))
568 : {
569 0 : return false;
570 : }
571 :
572 : // We must take the opposite since poCT transforms from
573 : // WGS84 to geoid. And we are going to do the reverse:
574 : // take an elevation over the geoid and transforms it to WGS84.
575 189912 : if (psTransform->bApplyDEMVDatumShift)
576 189912 : dfVDatumShift = -dfZ;
577 : }
578 :
579 296262 : bool bRetried = false;
580 296270 : retry:
581 296270 : GDALApplyGeoTransform(psTransform->adfDEMReverseGeoTransform, dfXTemp,
582 : dfYTemp, &dfX, &dfY);
583 296270 : if (pdfDEMPixel)
584 43197 : *pdfDEMPixel = dfX;
585 296270 : if (pdfDEMLine)
586 43197 : *pdfDEMLine = dfY;
587 :
588 296270 : if (!GDALRPCGetDEMHeight(psTransform, dfX, dfY, &dfDEMH))
589 : {
590 : // Try to handle the case where the DEM is in LL WGS84 and spans
591 : // over [-180,180], (or very close to it ), presumably with much
592 : // hole in the middle if using VRT, and the longitude goes beyond
593 : // that interval.
594 44616 : if (!bRetried && psTransform->poCT == nullptr &&
595 16677 : (dfXIn >= 180.0 || dfXIn <= -180.0))
596 : {
597 8 : const int nRasterXSize = psTransform->poDS->GetRasterXSize();
598 8 : const double dfMinDEMLong = psTransform->adfDEMGeoTransform[0];
599 8 : const double dfMaxDEMLong =
600 8 : psTransform->adfDEMGeoTransform[0] +
601 8 : nRasterXSize * psTransform->adfDEMGeoTransform[1];
602 8 : if (fabs(dfMinDEMLong - -180) < 0.1 &&
603 8 : fabs(dfMaxDEMLong - 180) < 0.1)
604 : {
605 8 : if (dfXIn >= 180)
606 : {
607 4 : dfXTemp = dfXIn - 360;
608 4 : dfYTemp = dfYIn;
609 : }
610 : else
611 : {
612 4 : dfXTemp = dfXIn + 360;
613 4 : dfYTemp = dfYIn;
614 : }
615 8 : bRetried = true;
616 8 : goto retry;
617 : }
618 : }
619 :
620 44608 : if (psTransform->bHasDEMMissingValue)
621 11832 : dfDEMH = psTransform->dfDEMMissingValue;
622 : else
623 : {
624 32776 : return false;
625 : }
626 : }
627 : #ifdef DEBUG_VERBOSE_EXTRACT_DEM
628 : CPLDebug("RPC_DEM", "X=%f, Y=%f -> Z=%f", dfX, dfY, dfDEMH);
629 : #endif
630 : }
631 :
632 271342 : *pdfHeight = dfVDatumShift + (psTransform->dfHeightOffset +
633 271342 : dfDEMH * psTransform->dfHeightScale);
634 271342 : return true;
635 : }
636 :
637 : /************************************************************************/
638 : /* GDALCreateRPCTransformer() */
639 : /************************************************************************/
640 :
641 0 : void *GDALCreateRPCTransformerV1(GDALRPCInfoV1 *psRPCInfo, int bReversed,
642 : double dfPixErrThreshold, char **papszOptions)
643 :
644 : {
645 : GDALRPCInfoV2 sRPCInfo;
646 0 : memcpy(&sRPCInfo, psRPCInfo, sizeof(GDALRPCInfoV1));
647 0 : sRPCInfo.dfERR_BIAS = std::numeric_limits<double>::quiet_NaN();
648 0 : sRPCInfo.dfERR_RAND = std::numeric_limits<double>::quiet_NaN();
649 0 : return GDALCreateRPCTransformerV2(&sRPCInfo, bReversed, dfPixErrThreshold,
650 0 : papszOptions);
651 : }
652 :
653 : /**
654 : * Create an RPC based transformer.
655 : *
656 : * The geometric sensor model describing the physical relationship between
657 : * image coordinates and ground coordinates is known as a Rigorous Projection
658 : * Model. A Rigorous Projection Model expresses the mapping of the image space
659 : * coordinates of rows and columns (r,c) onto the object space reference
660 : * surface geodetic coordinates (long, lat, height).
661 : *
662 : * A RPC supports a generic description of the Rigorous Projection Models. The
663 : * approximation used by GDAL (RPC00) is a set of rational polynomials
664 : * expressing the normalized row and column values, (rn , cn), as a function of
665 : * normalized geodetic latitude, longitude, and height, (P, L, H), given a
666 : * set of normalized polynomial coefficients (LINE_NUM_COEF_n, LINE_DEN_COEF_n,
667 : * SAMP_NUM_COEF_n, SAMP_DEN_COEF_n). Normalized values, rather than actual
668 : * values are used in order to minimize introduction of errors during the
669 : * calculations. The transformation between row and column values (r,c), and
670 : * normalized row and column values (rn, cn), and between the geodetic
671 : * latitude, longitude, and height and normalized geodetic latitude,
672 : * longitude, and height (P, L, H), is defined by a set of normalizing
673 : * translations (offsets) and scales that ensure all values are contained in
674 : * the range -1 to +1.
675 : *
676 : * This function creates a GDALTransformFunc compatible transformer
677 : * for going between image pixel/line and long/lat/height coordinates
678 : * using RPCs. The RPCs are provided in a GDALRPCInfo structure which is
679 : * normally read from metadata using GDALExtractRPCInfo().
680 : *
681 : * GDAL RPC Metadata has the following entries (also described in GDAL RFC 22
682 : * and the GeoTIFF RPC document http://geotiff.maptools.org/rpc_prop.html .
683 : *
684 : * <ul>
685 : * <li>ERR_BIAS: Error - Bias. The RMS bias error in meters per horizontal axis
686 : * of all points in the image (-1.0 if unknown)
687 : * <li>ERR_RAND: Error - Random. RMS random error in meters per horizontal axis
688 : * of each point in the image (-1.0 if unknown)
689 : * <li>LINE_OFF: Line Offset
690 : * <li>SAMP_OFF: Sample Offset
691 : * <li>LAT_OFF: Geodetic Latitude Offset
692 : * <li>LONG_OFF: Geodetic Longitude Offset
693 : * <li>HEIGHT_OFF: Geodetic Height Offset
694 : * <li>LINE_SCALE: Line Scale
695 : * <li>SAMP_SCALE: Sample Scale
696 : * <li>LAT_SCALE: Geodetic Latitude Scale
697 : * <li>LONG_SCALE: Geodetic Longitude Scale
698 : * <li>HEIGHT_SCALE: Geodetic Height Scale
699 :
700 : * <li>LINE_NUM_COEFF (1-20): Line Numerator Coefficients. Twenty coefficients
701 : * for the polynomial in the Numerator of the rn equation. (space separated)
702 : * <li>LINE_DEN_COEFF (1-20): Line Denominator Coefficients. Twenty coefficients
703 : * for the polynomial in the Denominator of the rn equation. (space separated)
704 : * <li>SAMP_NUM_COEFF (1-20): Sample Numerator Coefficients. Twenty coefficients
705 : * for the polynomial in the Numerator of the cn equation. (space separated)
706 : * <li>SAMP_DEN_COEFF (1-20): Sample Denominator Coefficients. Twenty
707 : * coefficients for the polynomial in the Denominator of the cn equation. (space
708 : * separated)
709 : * </ul>
710 : *
711 : * The transformer normally maps from pixel/line/height to long/lat/height space
712 : * as a forward transformation though in RPC terms that would be considered
713 : * an inverse transformation (and is solved by iterative approximation using
714 : * long/lat/height to pixel/line transformations). The default direction can
715 : * be reversed by passing bReversed=TRUE.
716 : *
717 : * The iterative solution of pixel/line
718 : * to lat/long/height is currently run for up to 10 iterations or until
719 : * the apparent error is less than dfPixErrThreshold pixels. Passing zero
720 : * will not avoid all error, but will cause the operation to run for the maximum
721 : * number of iterations.
722 : *
723 : * Starting with GDAL 2.1, debugging of the RPC inverse transformer can be done
724 : * by setting the RPC_INVERSE_VERBOSE configuration option to YES (in which case
725 : * extra debug information will be displayed in the "RPC" debug category, so
726 : * requiring CPL_DEBUG to be also set) and/or by setting RPC_INVERSE_LOG to a
727 : * filename that will contain the content of iterations (this last option only
728 : * makes sense when debugging point by point, since each time
729 : * RPCInverseTransformPoint() is called, the file is rewritten).
730 : *
731 : * Additional options to the transformer can be supplied in papszOptions.
732 : *
733 : * Options:
734 : *
735 : * <ul>
736 : * <li> RPC_HEIGHT: a fixed height offset to be applied to all points passed
737 : * in. In this situation the Z passed into the transformation function is
738 : * assumed to be height above ground, and the RPC_HEIGHT is assumed to be
739 : * an average height above sea level for ground in the target scene.</li>
740 : *
741 : * <li> RPC_HEIGHT_SCALE: a factor used to multiply heights above ground.
742 : * Useful when elevation offsets of the DEM are not expressed in meters.</li>
743 : *
744 : * <li> RPC_DEM: the name of a GDAL dataset (a DEM file typically) used to
745 : * extract elevation offsets from. In this situation the Z passed into the
746 : * transformation function is assumed to be height above ground. This option
747 : * should be used in replacement of RPC_HEIGHT to provide a way of defining
748 : * a non uniform ground for the target scene</li>
749 : *
750 : * <li> RPC_DEMINTERPOLATION: the DEM interpolation ("near", "bilinear" or
751 : "cubic").
752 : * Default is "bilinear".</li>
753 : *
754 : * <li> RPC_DEM_MISSING_VALUE: value of DEM height that must be used in case
755 : * the DEM has nodata value at the sampling point, or if its extent does not
756 : * cover the requested coordinate. When not specified, missing values will cause
757 : * a failed transform.</li>
758 : *
759 : * <li> RPC_DEM_SRS: (GDAL >= 3.2) WKT SRS, or any string recognized by
760 : * OGRSpatialReference::SetFromUserInput(), to be used as an override for DEM
761 : SRS.
762 : * Useful if DEM SRS does not have an explicit vertical component. </li>
763 : *
764 : * <li> RPC_DEM_APPLY_VDATUM_SHIFT: whether the vertical component of a compound
765 : * SRS for the DEM should be used (when it is present). This is useful so as to
766 : * be able to transform the "raw" values from the DEM expressed with respect to
767 : * a geoid to the heights with respect to the WGS84 ellipsoid. When this is
768 : * enabled, the GTIFF_REPORT_COMPD_CS configuration option will be also set
769 : * temporarily so as to get the vertical information from GeoTIFF
770 : * files. Defaults to TRUE. (GDAL >= 2.1.0)</li>
771 : *
772 : * <li> RPC_PIXEL_ERROR_THRESHOLD: overrides the dfPixErrThreshold parameter, ie
773 : the error (measured in pixels) allowed in the
774 : * iterative solution of pixel/line to lat/long computations (the other way
775 : * is always exact given the equations). (GDAL >= 2.1.0)</li>
776 : *
777 : * <li> RPC_MAX_ITERATIONS: maximum number of iterations allowed in the
778 : * iterative solution of pixel/line to lat/long computations. Default value is
779 : * 10 in the absence of a DEM, or 20 if there is a DEM. (GDAL >= 2.1.0)</li>
780 : *
781 : * <li> RPC_FOOTPRINT: WKT or GeoJSON polygon (in long / lat coordinate space)
782 : * with a validity footprint for the RPC. Any coordinate transformation that
783 : * goes from or arrive outside this footprint will be considered invalid. This
784 : * is useful in situations where the RPC values become highly unstable outside
785 : * of the area on which they have been computed for, potentially leading to
786 : * undesirable "echoes" / false positives. This requires GDAL to be built
787 : against
788 : * GEOS.</li>
789 : *
790 : * </ul>
791 : *
792 : * @param psRPCInfo Definition of the RPC parameters.
793 : *
794 : * @param bReversed If true "forward" transformation will be lat/long to
795 : * pixel/line instead of the normal pixel/line to lat/long.
796 : *
797 : * @param dfPixErrThreshold the error (measured in pixels) allowed in the
798 : * iterative solution of pixel/line to lat/long computations (the other way
799 : * is always exact given the equations). Starting with GDAL 2.1, this may also
800 : * be set through the RPC_PIXEL_ERROR_THRESHOLD transformer option.
801 : * If a negative or null value is provided, then this defaults to 0.1 pixel.
802 : *
803 : * @param papszOptions Other transformer options (i.e. RPC_HEIGHT=z).
804 : *
805 : * @return transformer callback data (deallocate with GDALDestroyTransformer()).
806 : */
807 :
808 46 : void *GDALCreateRPCTransformerV2(const GDALRPCInfoV2 *psRPCInfo, int bReversed,
809 : double dfPixErrThreshold, char **papszOptions)
810 :
811 : {
812 : /* -------------------------------------------------------------------- */
813 : /* Initialize core info. */
814 : /* -------------------------------------------------------------------- */
815 : GDALRPCTransformInfo *psTransform = static_cast<GDALRPCTransformInfo *>(
816 46 : CPLCalloc(sizeof(GDALRPCTransformInfo), 1));
817 :
818 46 : memcpy(&(psTransform->sRPC), psRPCInfo, sizeof(GDALRPCInfoV2));
819 46 : psTransform->bReversed = bReversed;
820 : const char *pszPixErrThreshold =
821 46 : CSLFetchNameValue(papszOptions, "RPC_PIXEL_ERROR_THRESHOLD");
822 46 : if (pszPixErrThreshold != nullptr)
823 1 : psTransform->dfPixErrThreshold = CPLAtof(pszPixErrThreshold);
824 45 : else if (dfPixErrThreshold > 0)
825 1 : psTransform->dfPixErrThreshold = dfPixErrThreshold;
826 : else
827 44 : psTransform->dfPixErrThreshold = DEFAULT_PIX_ERR_THRESHOLD;
828 46 : psTransform->dfHeightOffset = 0.0;
829 46 : psTransform->dfHeightScale = 1.0;
830 :
831 46 : memcpy(psTransform->sTI.abySignature, GDAL_GTI2_SIGNATURE,
832 : strlen(GDAL_GTI2_SIGNATURE));
833 46 : psTransform->sTI.pszClassName = GDAL_RPC_TRANSFORMER_CLASS_NAME;
834 46 : psTransform->sTI.pfnTransform = GDALRPCTransform;
835 46 : psTransform->sTI.pfnCleanup = GDALDestroyRPCTransformer;
836 46 : psTransform->sTI.pfnSerialize = GDALSerializeRPCTransformer;
837 46 : psTransform->sTI.pfnCreateSimilar = GDALCreateSimilarRPCTransformer;
838 :
839 : #ifdef USE_SSE2_OPTIM
840 : // Make sure padfCoeffs is aligned on a 16-byte boundary for SSE2 aligned
841 : // loads.
842 46 : psTransform->padfCoeffs =
843 46 : psTransform->adfDoubles +
844 46 : (reinterpret_cast<GUIntptr_t>(psTransform->adfDoubles) % 16) / 8;
845 46 : memcpy(psTransform->padfCoeffs, psRPCInfo->adfLINE_NUM_COEFF,
846 : 20 * sizeof(double));
847 46 : memcpy(psTransform->padfCoeffs + 20, psRPCInfo->adfLINE_DEN_COEFF,
848 : 20 * sizeof(double));
849 46 : memcpy(psTransform->padfCoeffs + 40, psRPCInfo->adfSAMP_NUM_COEFF,
850 : 20 * sizeof(double));
851 46 : memcpy(psTransform->padfCoeffs + 60, psRPCInfo->adfSAMP_DEN_COEFF,
852 : 20 * sizeof(double));
853 : #endif
854 :
855 : /* -------------------------------------------------------------------- */
856 : /* Do we have a "average height" that we want to consider all */
857 : /* elevations to be relative to? */
858 : /* -------------------------------------------------------------------- */
859 46 : const char *pszHeight = CSLFetchNameValue(papszOptions, "RPC_HEIGHT");
860 46 : if (pszHeight != nullptr)
861 6 : psTransform->dfHeightOffset = CPLAtof(pszHeight);
862 :
863 : /* -------------------------------------------------------------------- */
864 : /* The "height scale" */
865 : /* -------------------------------------------------------------------- */
866 : const char *pszHeightScale =
867 46 : CSLFetchNameValue(papszOptions, "RPC_HEIGHT_SCALE");
868 46 : if (pszHeightScale != nullptr)
869 7 : psTransform->dfHeightScale = CPLAtof(pszHeightScale);
870 :
871 : /* -------------------------------------------------------------------- */
872 : /* The DEM file name */
873 : /* -------------------------------------------------------------------- */
874 46 : const char *pszDEMPath = CSLFetchNameValue(papszOptions, "RPC_DEM");
875 46 : if (pszDEMPath != nullptr)
876 : {
877 33 : psTransform->pszDEMPath = CPLStrdup(pszDEMPath);
878 : }
879 :
880 : /* -------------------------------------------------------------------- */
881 : /* The DEM interpolation */
882 : /* -------------------------------------------------------------------- */
883 : const char *pszDEMInterpolation =
884 46 : CSLFetchNameValueDef(papszOptions, "RPC_DEMINTERPOLATION", "bilinear");
885 46 : if (EQUAL(pszDEMInterpolation, "near"))
886 : {
887 3 : psTransform->eResampleAlg = DRA_NearestNeighbour;
888 : }
889 43 : else if (EQUAL(pszDEMInterpolation, "bilinear"))
890 : {
891 40 : psTransform->eResampleAlg = DRA_Bilinear;
892 : }
893 3 : else if (EQUAL(pszDEMInterpolation, "cubic"))
894 : {
895 3 : psTransform->eResampleAlg = DRA_CubicSpline;
896 : }
897 : else
898 : {
899 0 : CPLDebug("RPC", "Unknown interpolation %s. Defaulting to bilinear",
900 : pszDEMInterpolation);
901 0 : psTransform->eResampleAlg = DRA_Bilinear;
902 : }
903 :
904 : /* -------------------------------------------------------------------- */
905 : /* The DEM missing value */
906 : /* -------------------------------------------------------------------- */
907 : const char *pszDEMMissingValue =
908 46 : CSLFetchNameValue(papszOptions, "RPC_DEM_MISSING_VALUE");
909 46 : if (pszDEMMissingValue != nullptr)
910 : {
911 6 : psTransform->bHasDEMMissingValue = TRUE;
912 6 : psTransform->dfDEMMissingValue = CPLAtof(pszDEMMissingValue);
913 : }
914 :
915 : /* -------------------------------------------------------------------- */
916 : /* The DEM SRS override */
917 : /* -------------------------------------------------------------------- */
918 46 : const char *pszDEMSRS = CSLFetchNameValue(papszOptions, "RPC_DEM_SRS");
919 46 : if (pszDEMSRS != nullptr)
920 : {
921 1 : psTransform->pszDEMSRS = CPLStrdup(pszDEMSRS);
922 : }
923 :
924 : /* -------------------------------------------------------------------- */
925 : /* Whether to apply vdatum shift */
926 : /* -------------------------------------------------------------------- */
927 46 : psTransform->bApplyDEMVDatumShift =
928 46 : CPLFetchBool(papszOptions, "RPC_DEM_APPLY_VDATUM_SHIFT", true);
929 :
930 46 : psTransform->nMaxIterations =
931 46 : atoi(CSLFetchNameValueDef(papszOptions, "RPC_MAX_ITERATIONS", "0"));
932 :
933 : /* -------------------------------------------------------------------- */
934 : /* Debug */
935 : /* -------------------------------------------------------------------- */
936 46 : psTransform->bRPCInverseVerbose =
937 46 : CPLTestBool(CPLGetConfigOption("RPC_INVERSE_VERBOSE", "NO"));
938 : const char *pszRPCInverseLog =
939 46 : CPLGetConfigOption("RPC_INVERSE_LOG", nullptr);
940 46 : if (pszRPCInverseLog != nullptr)
941 1 : psTransform->pszRPCInverseLog = CPLStrdup(pszRPCInverseLog);
942 :
943 : /* -------------------------------------------------------------------- */
944 : /* Footprint */
945 : /* -------------------------------------------------------------------- */
946 46 : const char *pszFootprint = CSLFetchNameValue(papszOptions, "RPC_FOOTPRINT");
947 46 : if (pszFootprint != nullptr)
948 : {
949 1 : psTransform->pszRPCFootprint = CPLStrdup(pszFootprint);
950 1 : if (pszFootprint[0] == '{')
951 : {
952 0 : psTransform->poRPCFootprintGeom =
953 0 : OGRGeometryFactory::createFromGeoJson(pszFootprint);
954 : }
955 : else
956 : {
957 1 : OGRGeometryFactory::createFromWkt(
958 : pszFootprint, nullptr, &(psTransform->poRPCFootprintGeom));
959 : }
960 1 : if (psTransform->poRPCFootprintGeom)
961 : {
962 1 : if (OGRHasPreparedGeometrySupport())
963 : {
964 1 : psTransform->poRPCFootprintPreparedGeom =
965 1 : OGRCreatePreparedGeometry(
966 : OGRGeometry::ToHandle(psTransform->poRPCFootprintGeom));
967 : }
968 : else
969 : {
970 0 : CPLError(CE_Warning, CPLE_AppDefined,
971 : "GEOS not available. RPC_FOOTPRINT will be ignored");
972 : }
973 : }
974 : }
975 :
976 : /* -------------------------------------------------------------------- */
977 : /* Open DEM if needed. */
978 : /* -------------------------------------------------------------------- */
979 :
980 46 : if (psTransform->pszDEMPath != nullptr && !GDALRPCOpenDEM(psTransform))
981 : {
982 1 : GDALDestroyRPCTransformer(psTransform);
983 1 : return nullptr;
984 : }
985 :
986 : /* -------------------------------------------------------------------- */
987 : /* Establish a reference point for calcualating an affine */
988 : /* geotransform approximate transformation. */
989 : /* -------------------------------------------------------------------- */
990 45 : double adfGTFromLL[6] = {};
991 45 : double dfRefPixel = -1.0;
992 45 : double dfRefLine = -1.0;
993 45 : double dfRefLong = 0.0;
994 45 : double dfRefLat = 0.0;
995 :
996 45 : if (psRPCInfo->dfMIN_LONG != -180 || psRPCInfo->dfMAX_LONG != 180)
997 : {
998 1 : dfRefLong = (psRPCInfo->dfMIN_LONG + psRPCInfo->dfMAX_LONG) * 0.5;
999 1 : dfRefLat = (psRPCInfo->dfMIN_LAT + psRPCInfo->dfMAX_LAT) * 0.5;
1000 :
1001 1 : double dfX = dfRefLong;
1002 1 : double dfY = dfRefLat;
1003 1 : double dfZ = 0.0;
1004 1 : int nSuccess = 0;
1005 : // Try with DEM first.
1006 1 : if (GDALRPCTransform(psTransform, !(psTransform->bReversed), 1, &dfX,
1007 2 : &dfY, &dfZ, &nSuccess) &&
1008 1 : nSuccess)
1009 : {
1010 1 : dfRefPixel = dfX;
1011 1 : dfRefLine = dfY;
1012 : }
1013 : else
1014 : {
1015 0 : RPCTransformPoint(psTransform, dfRefLong, dfRefLat, 0.0,
1016 : &dfRefPixel, &dfRefLine);
1017 : }
1018 : }
1019 :
1020 : // Try with scale and offset if we don't can't use bounds or
1021 : // the results seem daft.
1022 45 : if (dfRefPixel < 0.0 || dfRefLine < 0.0 || dfRefPixel > 100000 ||
1023 1 : dfRefLine > 100000)
1024 : {
1025 44 : dfRefLong = psRPCInfo->dfLONG_OFF;
1026 44 : dfRefLat = psRPCInfo->dfLAT_OFF;
1027 :
1028 44 : double dfX = dfRefLong;
1029 44 : double dfY = dfRefLat;
1030 44 : double dfZ = 0.0;
1031 44 : int nSuccess = 0;
1032 : // Try with DEM first.
1033 44 : if (GDALRPCTransform(psTransform, !(psTransform->bReversed), 1, &dfX,
1034 88 : &dfY, &dfZ, &nSuccess) &&
1035 44 : nSuccess)
1036 : {
1037 28 : dfRefPixel = dfX;
1038 28 : dfRefLine = dfY;
1039 : }
1040 : else
1041 : {
1042 16 : RPCTransformPoint(psTransform, dfRefLong, dfRefLat, 0.0,
1043 : &dfRefPixel, &dfRefLine);
1044 : }
1045 : }
1046 :
1047 45 : psTransform->dfRefZ = 0.0;
1048 45 : GDALRPCGetHeightAtLongLat(psTransform, dfRefLong, dfRefLat,
1049 : &psTransform->dfRefZ);
1050 :
1051 : /* -------------------------------------------------------------------- */
1052 : /* Transform nearby locations to establish affine direction */
1053 : /* vectors. */
1054 : /* -------------------------------------------------------------------- */
1055 45 : double dfRefPixelDelta = 0.0;
1056 45 : double dfRefLineDelta = 0.0;
1057 45 : double dfLLDelta = 0.0001;
1058 :
1059 45 : RPCTransformPoint(psTransform, dfRefLong + dfLLDelta, dfRefLat,
1060 : psTransform->dfRefZ, &dfRefPixelDelta, &dfRefLineDelta);
1061 45 : adfGTFromLL[1] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
1062 45 : adfGTFromLL[4] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
1063 :
1064 45 : RPCTransformPoint(psTransform, dfRefLong, dfRefLat + dfLLDelta,
1065 : psTransform->dfRefZ, &dfRefPixelDelta, &dfRefLineDelta);
1066 45 : adfGTFromLL[2] = (dfRefPixelDelta - dfRefPixel) / dfLLDelta;
1067 45 : adfGTFromLL[5] = (dfRefLineDelta - dfRefLine) / dfLLDelta;
1068 :
1069 45 : adfGTFromLL[0] =
1070 45 : dfRefPixel - adfGTFromLL[1] * dfRefLong - adfGTFromLL[2] * dfRefLat;
1071 45 : adfGTFromLL[3] =
1072 45 : dfRefLine - adfGTFromLL[4] * dfRefLong - adfGTFromLL[5] * dfRefLat;
1073 :
1074 45 : if (!GDALInvGeoTransform(adfGTFromLL,
1075 45 : psTransform->adfPLToLatLongGeoTransform))
1076 : {
1077 0 : CPLError(CE_Failure, CPLE_AppDefined, "Cannot invert geotransform");
1078 0 : GDALDestroyRPCTransformer(psTransform);
1079 0 : return nullptr;
1080 : }
1081 :
1082 45 : return psTransform;
1083 : }
1084 :
1085 : /************************************************************************/
1086 : /* GDALDestroyReprojectionTransformer() */
1087 : /************************************************************************/
1088 :
1089 : /** Destroy RPC transformer */
1090 46 : void GDALDestroyRPCTransformer(void *pTransformAlg)
1091 :
1092 : {
1093 46 : if (pTransformAlg == nullptr)
1094 0 : return;
1095 :
1096 46 : GDALRPCTransformInfo *psTransform =
1097 : static_cast<GDALRPCTransformInfo *>(pTransformAlg);
1098 :
1099 46 : CPLFree(psTransform->pszDEMPath);
1100 46 : CPLFree(psTransform->pszDEMSRS);
1101 :
1102 46 : if (psTransform->poDS)
1103 32 : GDALClose(psTransform->poDS);
1104 46 : delete psTransform->poCacheDEM;
1105 46 : if (psTransform->poCT)
1106 11 : OCTDestroyCoordinateTransformation(
1107 11 : reinterpret_cast<OGRCoordinateTransformationH>(psTransform->poCT));
1108 46 : CPLFree(psTransform->pszRPCInverseLog);
1109 :
1110 46 : CPLFree(psTransform->pszRPCFootprint);
1111 46 : delete psTransform->poRPCFootprintGeom;
1112 46 : OGRDestroyPreparedGeometry(psTransform->poRPCFootprintPreparedGeom);
1113 :
1114 46 : CPLFree(pTransformAlg);
1115 : }
1116 :
1117 : /************************************************************************/
1118 : /* RPCInverseTransformPoint() */
1119 : /************************************************************************/
1120 :
1121 13614 : static bool RPCInverseTransformPoint(GDALRPCTransformInfo *psTransform,
1122 : double dfPixel, double dfLine,
1123 : double dfUserHeight, double *pdfLong,
1124 : double *pdfLat)
1125 :
1126 : {
1127 : // Memo:
1128 : // Known to work with 40 iterations with DEM on all points (int coord and
1129 : // +0.5,+0.5 shift) of flock1.20160216_041050_0905.tif, especially on (0,0).
1130 :
1131 : /* -------------------------------------------------------------------- */
1132 : /* Compute an initial approximation based on linear */
1133 : /* interpolation from our reference point. */
1134 : /* -------------------------------------------------------------------- */
1135 13614 : double dfResultX = psTransform->adfPLToLatLongGeoTransform[0] +
1136 13614 : psTransform->adfPLToLatLongGeoTransform[1] * dfPixel +
1137 13614 : psTransform->adfPLToLatLongGeoTransform[2] * dfLine;
1138 :
1139 13614 : double dfResultY = psTransform->adfPLToLatLongGeoTransform[3] +
1140 13614 : psTransform->adfPLToLatLongGeoTransform[4] * dfPixel +
1141 13614 : psTransform->adfPLToLatLongGeoTransform[5] * dfLine;
1142 :
1143 13614 : if (psTransform->bRPCInverseVerbose)
1144 : {
1145 1 : CPLDebug("RPC", "Computing inverse transform for (pixel,line)=(%f,%f)",
1146 : dfPixel, dfLine);
1147 : }
1148 13614 : VSILFILE *fpLog = nullptr;
1149 13614 : if (psTransform->pszRPCInverseLog)
1150 : {
1151 1 : fpLog = VSIFOpenL(
1152 2 : CPLResetExtensionSafe(psTransform->pszRPCInverseLog, "csvt")
1153 : .c_str(),
1154 : "wb");
1155 1 : if (fpLog != nullptr)
1156 : {
1157 1 : VSIFPrintfL(fpLog, "Integer,Real,Real,Real,String,Real,Real\n");
1158 1 : VSIFCloseL(fpLog);
1159 : }
1160 1 : fpLog = VSIFOpenL(psTransform->pszRPCInverseLog, "wb");
1161 1 : if (fpLog != nullptr)
1162 1 : VSIFPrintfL(
1163 : fpLog,
1164 : "iter,long,lat,height,WKT,error_pixel_x,error_pixel_y\n");
1165 : }
1166 :
1167 : /* -------------------------------------------------------------------- */
1168 : /* Now iterate, trying to find a closer LL location that will */
1169 : /* back transform to the indicated pixel and line. */
1170 : /* -------------------------------------------------------------------- */
1171 13614 : double dfPixelDeltaX = 0.0;
1172 13614 : double dfPixelDeltaY = 0.0;
1173 13614 : double dfLastResultX = 0.0;
1174 13614 : double dfLastResultY = 0.0;
1175 13614 : double dfLastPixelDeltaX = 0.0;
1176 13614 : double dfLastPixelDeltaY = 0.0;
1177 13614 : bool bLastPixelDeltaValid = false;
1178 27228 : const int nMaxIterations = (psTransform->nMaxIterations > 0)
1179 26237 : ? psTransform->nMaxIterations
1180 12623 : : (psTransform->poDS != nullptr) ? 20
1181 : : 10;
1182 13614 : int nCountConsecutiveErrorBelow2 = 0;
1183 :
1184 13614 : int iIter = 0; // Used after for.
1185 47918 : for (; iIter < nMaxIterations; iIter++)
1186 : {
1187 47388 : double dfBackPixel = 0.0;
1188 47388 : double dfBackLine = 0.0;
1189 :
1190 : // Update DEMH.
1191 47388 : double dfDEMH = 0.0;
1192 47388 : double dfDEMPixel = 0.0;
1193 47388 : double dfDEMLine = 0.0;
1194 47388 : if (!GDALRPCGetHeightAtLongLat(psTransform, dfResultX, dfResultY,
1195 : &dfDEMH, &dfDEMPixel, &dfDEMLine))
1196 : {
1197 20334 : if (psTransform->poDS)
1198 : {
1199 20334 : CPLDebug("RPC", "DEM (pixel, line) = (%g, %g)", dfDEMPixel,
1200 : dfDEMLine);
1201 : }
1202 :
1203 : // The first time, the guess might be completely out of the
1204 : // validity of the DEM, so pickup the "reference Z" as the
1205 : // first guess or the closest point of the DEM by snapping to it.
1206 20334 : if (iIter == 0)
1207 : {
1208 10262 : bool bUseRefZ = true;
1209 10262 : if (psTransform->poDS)
1210 : {
1211 10262 : if (dfDEMPixel >= psTransform->poDS->GetRasterXSize())
1212 4311 : dfDEMPixel = psTransform->poDS->GetRasterXSize() - 0.5;
1213 5951 : else if (dfDEMPixel < 0)
1214 1385 : dfDEMPixel = 0.5;
1215 10262 : if (dfDEMLine >= psTransform->poDS->GetRasterYSize())
1216 21 : dfDEMLine = psTransform->poDS->GetRasterYSize() - 0.5;
1217 10241 : else if (dfDEMPixel < 0)
1218 0 : dfDEMPixel = 0.5;
1219 10262 : if (GDALRPCGetDEMHeight(psTransform, dfDEMPixel, dfDEMLine,
1220 10262 : &dfDEMH))
1221 : {
1222 1900 : bUseRefZ = false;
1223 1900 : CPLDebug("RPC",
1224 : "Iteration %d for (pixel, line) = (%g, %g): "
1225 : "No elevation value at %.15g %.15g. "
1226 : "Using elevation %g at DEM (pixel, line) = "
1227 : "(%g, %g) (snapping to boundaries) instead",
1228 : iIter, dfPixel, dfLine, dfResultX, dfResultY,
1229 : dfDEMH, dfDEMPixel, dfDEMLine);
1230 : }
1231 : }
1232 10262 : if (bUseRefZ)
1233 : {
1234 8362 : dfDEMH = psTransform->dfRefZ;
1235 8362 : CPLDebug("RPC",
1236 : "Iteration %d for (pixel, line) = (%g, %g): "
1237 : "No elevation value at %.15g %.15g. "
1238 : "Using elevation %g of reference point instead",
1239 : iIter, dfPixel, dfLine, dfResultX, dfResultY,
1240 : dfDEMH);
1241 : }
1242 : }
1243 : else
1244 : {
1245 10072 : CPLDebug("RPC",
1246 : "Iteration %d for (pixel, line) = (%g, %g): "
1247 : "No elevation value at %.15g %.15g. Erroring out",
1248 : iIter, dfPixel, dfLine, dfResultX, dfResultY);
1249 10072 : if (fpLog)
1250 0 : VSIFCloseL(fpLog);
1251 10072 : return false;
1252 : }
1253 : }
1254 :
1255 37316 : RPCTransformPoint(psTransform, dfResultX, dfResultY,
1256 : dfUserHeight + dfDEMH, &dfBackPixel, &dfBackLine);
1257 :
1258 37316 : dfPixelDeltaX = dfBackPixel - dfPixel;
1259 37316 : dfPixelDeltaY = dfBackLine - dfLine;
1260 :
1261 37316 : if (psTransform->bRPCInverseVerbose)
1262 : {
1263 8 : CPLDebug("RPC",
1264 : "Iter %d: dfPixelDeltaX=%.02f, dfPixelDeltaY=%.02f, "
1265 : "long=%f, lat=%f, height=%f",
1266 : iIter, dfPixelDeltaX, dfPixelDeltaY, dfResultX, dfResultY,
1267 : dfUserHeight + dfDEMH);
1268 : }
1269 37316 : if (fpLog != nullptr)
1270 : {
1271 8 : VSIFPrintfL(fpLog,
1272 : "%d,%.12f,%.12f,%f,\"POINT(%.12f %.12f)\",%f,%f\n",
1273 : iIter, dfResultX, dfResultY, dfUserHeight + dfDEMH,
1274 : dfResultX, dfResultY, dfPixelDeltaX, dfPixelDeltaY);
1275 : }
1276 :
1277 : const double dfError =
1278 37316 : std::max(std::abs(dfPixelDeltaX), std::abs(dfPixelDeltaY));
1279 37316 : if (dfError < psTransform->dfPixErrThreshold)
1280 : {
1281 3012 : iIter = -1;
1282 3012 : if (psTransform->bRPCInverseVerbose)
1283 : {
1284 1 : CPLDebug("RPC", "Converged!");
1285 : }
1286 3012 : break;
1287 : }
1288 34304 : else if (psTransform->poDS != nullptr && bLastPixelDeltaValid &&
1289 21058 : dfPixelDeltaX * dfLastPixelDeltaX < 0 &&
1290 106 : dfPixelDeltaY * dfLastPixelDeltaY < 0)
1291 : {
1292 : // When there is a DEM, if the error changes sign, we might
1293 : // oscillate forever, so take a mean position as a new guess.
1294 40 : if (psTransform->bRPCInverseVerbose)
1295 : {
1296 1 : CPLDebug("RPC",
1297 : "Oscillation detected. "
1298 : "Taking mean of 2 previous results as new guess");
1299 : }
1300 40 : dfResultX = (fabs(dfPixelDeltaX) * dfLastResultX +
1301 40 : fabs(dfLastPixelDeltaX) * dfResultX) /
1302 40 : (fabs(dfPixelDeltaX) + fabs(dfLastPixelDeltaX));
1303 40 : dfResultY = (fabs(dfPixelDeltaY) * dfLastResultY +
1304 40 : fabs(dfLastPixelDeltaY) * dfResultY) /
1305 40 : (fabs(dfPixelDeltaY) + fabs(dfLastPixelDeltaY));
1306 40 : bLastPixelDeltaValid = false;
1307 40 : nCountConsecutiveErrorBelow2 = 0;
1308 40 : continue;
1309 : }
1310 :
1311 34264 : double dfBoostFactor = 1.0;
1312 34264 : if (psTransform->poDS != nullptr && nCountConsecutiveErrorBelow2 >= 5 &&
1313 : dfError < 2)
1314 : {
1315 : // When there is a DEM, if we remain below a given threshold
1316 : // (somewhat arbitrarily set to 2 pixels) for some time, apply a
1317 : // "boost factor" for the new guessed result, in the hope we will go
1318 : // out of the somewhat current stuck situation.
1319 13 : dfBoostFactor = 10;
1320 13 : if (psTransform->bRPCInverseVerbose)
1321 : {
1322 1 : CPLDebug("RPC", "Applying boost factor 10");
1323 : }
1324 : }
1325 :
1326 34264 : if (dfError < 2)
1327 1340 : nCountConsecutiveErrorBelow2++;
1328 : else
1329 32924 : nCountConsecutiveErrorBelow2 = 0;
1330 :
1331 34264 : const double dfNewResultX =
1332 34264 : dfResultX -
1333 34264 : (dfPixelDeltaX * psTransform->adfPLToLatLongGeoTransform[1] *
1334 : dfBoostFactor) -
1335 34264 : (dfPixelDeltaY * psTransform->adfPLToLatLongGeoTransform[2] *
1336 : dfBoostFactor);
1337 34264 : const double dfNewResultY =
1338 34264 : dfResultY -
1339 34264 : (dfPixelDeltaX * psTransform->adfPLToLatLongGeoTransform[4] *
1340 : dfBoostFactor) -
1341 34264 : (dfPixelDeltaY * psTransform->adfPLToLatLongGeoTransform[5] *
1342 : dfBoostFactor);
1343 :
1344 34264 : dfLastResultX = dfResultX;
1345 34264 : dfLastResultY = dfResultY;
1346 34264 : dfResultX = dfNewResultX;
1347 34264 : dfResultY = dfNewResultY;
1348 34264 : dfLastPixelDeltaX = dfPixelDeltaX;
1349 34264 : dfLastPixelDeltaY = dfPixelDeltaY;
1350 34264 : bLastPixelDeltaValid = true;
1351 : }
1352 3542 : if (fpLog != nullptr)
1353 1 : VSIFCloseL(fpLog);
1354 :
1355 3542 : if (iIter != -1)
1356 : {
1357 530 : CPLDebug("RPC", "Failed Iterations %d: Got: %.16g,%.16g Offset=%g,%g",
1358 : iIter, dfResultX, dfResultY, dfPixelDeltaX, dfPixelDeltaY);
1359 530 : return false;
1360 : }
1361 :
1362 3012 : *pdfLong = dfResultX;
1363 3012 : *pdfLat = dfResultY;
1364 3012 : return true;
1365 : }
1366 :
1367 : /************************************************************************/
1368 : /* GDALRPCGetDEMHeight() */
1369 : /************************************************************************/
1370 :
1371 306532 : static int GDALRPCGetDEMHeight(GDALRPCTransformInfo *psTransform,
1372 : const double dfXIn, const double dfYIn,
1373 : double *pdfDEMH)
1374 : {
1375 306532 : GDALRIOResampleAlg eResample = GDALRIOResampleAlg::GRIORA_NearestNeighbour;
1376 306532 : switch (psTransform->eResampleAlg)
1377 : {
1378 14 : case DEMResampleAlg::DRA_NearestNeighbour:
1379 14 : eResample = GDALRIOResampleAlg::GRIORA_NearestNeighbour;
1380 14 : break;
1381 306504 : case DEMResampleAlg::DRA_Bilinear:
1382 306504 : eResample = GDALRIOResampleAlg::GRIORA_Bilinear;
1383 306504 : break;
1384 14 : case DEMResampleAlg::DRA_CubicSpline:
1385 14 : eResample = GDALRIOResampleAlg::GRIORA_CubicSpline;
1386 14 : break;
1387 : }
1388 :
1389 306532 : std::unique_ptr<DoublePointsCache> cacheDEM{psTransform->poCacheDEM};
1390 : int res =
1391 306532 : GDALInterpolateAtPoint(psTransform->poDS->GetRasterBand(1), eResample,
1392 306532 : cacheDEM, dfXIn, dfYIn, pdfDEMH, nullptr);
1393 306532 : psTransform->poCacheDEM = cacheDEM.release();
1394 613064 : return res;
1395 : }
1396 :
1397 : /************************************************************************/
1398 : /* RPCIsValidLongLat() */
1399 : /************************************************************************/
1400 :
1401 547280 : static bool RPCIsValidLongLat(const GDALRPCTransformInfo *psTransform,
1402 : double dfLong, double dfLat)
1403 : {
1404 547280 : if (!psTransform->poRPCFootprintPreparedGeom)
1405 427167 : return true;
1406 :
1407 240226 : OGRPoint p(dfLong, dfLat);
1408 240226 : return CPL_TO_BOOL(OGRPreparedGeometryContains(
1409 240226 : psTransform->poRPCFootprintPreparedGeom, OGRGeometry::ToHandle(&p)));
1410 : }
1411 :
1412 : /************************************************************************/
1413 : /* GDALRPCTransformWholeLineWithDEM() */
1414 : /************************************************************************/
1415 :
1416 : static int
1417 2192 : GDALRPCTransformWholeLineWithDEM(const GDALRPCTransformInfo *psTransform,
1418 : int nPointCount, double *padfX, double *padfY,
1419 : double *padfZ, int *panSuccess, int nXLeft,
1420 : int nXWidth, int nYTop, int nYHeight)
1421 : {
1422 : double *padfDEMBuffer = static_cast<double *>(
1423 2192 : VSI_MALLOC3_VERBOSE(sizeof(double), nXWidth, nYHeight));
1424 2192 : if (padfDEMBuffer == nullptr)
1425 : {
1426 0 : for (int i = 0; i < nPointCount; i++)
1427 0 : panSuccess[i] = FALSE;
1428 0 : return FALSE;
1429 : }
1430 2192 : CPLErr eErr = psTransform->poDS->GetRasterBand(1)->RasterIO(
1431 : GF_Read, nXLeft, nYTop, nXWidth, nYHeight, padfDEMBuffer, nXWidth,
1432 : nYHeight, GDT_Float64, 0, 0, nullptr);
1433 2192 : if (eErr != CE_None)
1434 : {
1435 0 : for (int i = 0; i < nPointCount; i++)
1436 0 : panSuccess[i] = FALSE;
1437 0 : VSIFree(padfDEMBuffer);
1438 0 : return FALSE;
1439 : }
1440 :
1441 2192 : int bGotNoDataValue = FALSE;
1442 : const double dfNoDataValue =
1443 2192 : psTransform->poDS->GetRasterBand(1)->GetNoDataValue(&bGotNoDataValue);
1444 :
1445 : // dfY in pixel center convention.
1446 2192 : const double dfY = psTransform->adfDEMReverseGeoTransform[3] +
1447 2192 : padfY[0] * psTransform->adfDEMReverseGeoTransform[5] -
1448 : 0.5;
1449 2192 : const int nY = static_cast<int>(dfY);
1450 2192 : const double dfDeltaY = dfY - nY;
1451 :
1452 180781 : for (int i = 0; i < nPointCount; i++)
1453 : {
1454 178589 : if (padfX[i] == HUGE_VAL)
1455 0 : continue;
1456 :
1457 178589 : double dfDEMH = 0.0;
1458 178589 : const double dfZ_i = padfZ ? padfZ[i] : 0.0;
1459 :
1460 178589 : if (psTransform->eResampleAlg == DRA_CubicSpline)
1461 : {
1462 : // dfX in pixel center convention.
1463 10 : const double dfX =
1464 10 : psTransform->adfDEMReverseGeoTransform[0] +
1465 10 : padfX[i] * psTransform->adfDEMReverseGeoTransform[1] - 0.5;
1466 10 : const int nX = static_cast<int>(dfX);
1467 10 : const double dfDeltaX = dfX - nX;
1468 :
1469 10 : const int nXNew = nX - 1;
1470 :
1471 10 : double dfSumH = 0.0;
1472 10 : double dfSumWeight = 0.0;
1473 50 : for (int k_i = 0; k_i < 4; k_i++)
1474 : {
1475 : // Loop across the X axis.
1476 200 : for (int k_j = 0; k_j < 4; k_j++)
1477 : {
1478 : // Calculate the weight for the specified pixel according
1479 : // to the bicubic b-spline kernel we're using for
1480 : // interpolation.
1481 160 : const int dKernIndX = k_j - 1;
1482 160 : const int dKernIndY = k_i - 1;
1483 : const double dfPixelWeight =
1484 160 : CubicSplineKernel(dKernIndX - dfDeltaX) *
1485 160 : CubicSplineKernel(dKernIndY - dfDeltaY);
1486 :
1487 : // Create a sum of all values
1488 : // adjusted for the pixel's calculated weight.
1489 160 : const double dfElev =
1490 160 : padfDEMBuffer[k_i * nXWidth + nXNew - nXLeft + k_j];
1491 160 : if (bGotNoDataValue &&
1492 160 : ARE_REAL_EQUAL(dfNoDataValue, dfElev))
1493 0 : continue;
1494 :
1495 160 : dfSumH += dfElev * dfPixelWeight;
1496 160 : dfSumWeight += dfPixelWeight;
1497 : }
1498 : }
1499 10 : if (dfSumWeight == 0.0)
1500 : {
1501 0 : if (psTransform->bHasDEMMissingValue)
1502 0 : dfDEMH = psTransform->dfDEMMissingValue;
1503 : else
1504 : {
1505 0 : panSuccess[i] = FALSE;
1506 0 : continue;
1507 : }
1508 : }
1509 : else
1510 10 : dfDEMH = dfSumH / dfSumWeight;
1511 : }
1512 178579 : else if (psTransform->eResampleAlg == DRA_Bilinear)
1513 : {
1514 : // dfX in pixel center convention.
1515 178569 : const double dfX =
1516 178569 : psTransform->adfDEMReverseGeoTransform[0] +
1517 178569 : padfX[i] * psTransform->adfDEMReverseGeoTransform[1] - 0.5;
1518 178569 : const int nX = static_cast<int>(dfX);
1519 178569 : const double dfDeltaX = dfX - nX;
1520 :
1521 : // Bilinear interpolation.
1522 178569 : double adfElevData[4] = {};
1523 178569 : memcpy(adfElevData, padfDEMBuffer + nX - nXLeft,
1524 : 2 * sizeof(double));
1525 178569 : memcpy(adfElevData + 2, padfDEMBuffer + nXWidth + nX - nXLeft,
1526 : 2 * sizeof(double));
1527 :
1528 178569 : int bFoundNoDataElev = FALSE;
1529 178569 : if (bGotNoDataValue)
1530 : {
1531 46350 : int k_valid_sample = -1;
1532 231750 : for (int k_i = 0; k_i < 4; k_i++)
1533 : {
1534 185400 : if (ARE_REAL_EQUAL(dfNoDataValue, adfElevData[k_i]))
1535 : {
1536 0 : bFoundNoDataElev = TRUE;
1537 : }
1538 185400 : else if (k_valid_sample < 0)
1539 : {
1540 46350 : k_valid_sample = k_i;
1541 : }
1542 : }
1543 46350 : if (bFoundNoDataElev)
1544 : {
1545 0 : if (k_valid_sample >= 0)
1546 : {
1547 0 : if (!RPCIsValidLongLat(psTransform, padfX[i], padfY[i]))
1548 : {
1549 0 : panSuccess[i] = FALSE;
1550 0 : padfX[i] = HUGE_VAL;
1551 0 : padfY[i] = HUGE_VAL;
1552 0 : continue;
1553 : }
1554 0 : dfDEMH = adfElevData[k_valid_sample];
1555 0 : RPCTransformPoint(
1556 0 : psTransform, padfX[i], padfY[i],
1557 0 : dfZ_i + (psTransform->dfHeightOffset + dfDEMH) *
1558 0 : psTransform->dfHeightScale,
1559 0 : padfX + i, padfY + i);
1560 :
1561 0 : panSuccess[i] = TRUE;
1562 0 : continue;
1563 : }
1564 0 : else if (psTransform->bHasDEMMissingValue)
1565 : {
1566 0 : if (!RPCIsValidLongLat(psTransform, padfX[i], padfY[i]))
1567 : {
1568 0 : panSuccess[i] = FALSE;
1569 0 : padfX[i] = HUGE_VAL;
1570 0 : padfY[i] = HUGE_VAL;
1571 0 : continue;
1572 : }
1573 0 : dfDEMH = psTransform->dfDEMMissingValue;
1574 0 : RPCTransformPoint(
1575 0 : psTransform, padfX[i], padfY[i],
1576 0 : dfZ_i + (psTransform->dfHeightOffset + dfDEMH) *
1577 0 : psTransform->dfHeightScale,
1578 0 : padfX + i, padfY + i);
1579 :
1580 0 : panSuccess[i] = TRUE;
1581 0 : continue;
1582 : }
1583 : else
1584 : {
1585 0 : panSuccess[i] = FALSE;
1586 0 : padfX[i] = HUGE_VAL;
1587 0 : padfY[i] = HUGE_VAL;
1588 0 : continue;
1589 : }
1590 : }
1591 : }
1592 178569 : const double dfDeltaX1 = 1.0 - dfDeltaX;
1593 178569 : const double dfDeltaY1 = 1.0 - dfDeltaY;
1594 :
1595 178569 : const double dfXZ1 =
1596 178569 : adfElevData[0] * dfDeltaX1 + adfElevData[1] * dfDeltaX;
1597 178569 : const double dfXZ2 =
1598 178569 : adfElevData[2] * dfDeltaX1 + adfElevData[3] * dfDeltaX;
1599 178569 : const double dfYZ = dfXZ1 * dfDeltaY1 + dfXZ2 * dfDeltaY;
1600 178569 : dfDEMH = dfYZ;
1601 : }
1602 : else
1603 : {
1604 10 : const double dfX =
1605 10 : psTransform->adfDEMReverseGeoTransform[0] +
1606 10 : padfX[i] * psTransform->adfDEMReverseGeoTransform[1];
1607 10 : const int nX = int(dfX);
1608 :
1609 10 : dfDEMH = padfDEMBuffer[nX - nXLeft];
1610 10 : if (bGotNoDataValue && ARE_REAL_EQUAL(dfNoDataValue, dfDEMH))
1611 : {
1612 0 : if (psTransform->bHasDEMMissingValue)
1613 0 : dfDEMH = psTransform->dfDEMMissingValue;
1614 : else
1615 : {
1616 0 : panSuccess[i] = FALSE;
1617 0 : padfX[i] = HUGE_VAL;
1618 0 : padfY[i] = HUGE_VAL;
1619 0 : continue;
1620 : }
1621 : }
1622 : }
1623 :
1624 178589 : if (!RPCIsValidLongLat(psTransform, padfX[i], padfY[i]))
1625 : {
1626 0 : panSuccess[i] = FALSE;
1627 0 : padfX[i] = HUGE_VAL;
1628 0 : padfY[i] = HUGE_VAL;
1629 0 : continue;
1630 : }
1631 178589 : RPCTransformPoint(psTransform, padfX[i], padfY[i],
1632 178589 : dfZ_i + (psTransform->dfHeightOffset + dfDEMH) *
1633 178589 : psTransform->dfHeightScale,
1634 178589 : padfX + i, padfY + i);
1635 :
1636 178589 : panSuccess[i] = TRUE;
1637 : }
1638 :
1639 2192 : VSIFree(padfDEMBuffer);
1640 :
1641 2192 : return TRUE;
1642 : }
1643 :
1644 : /************************************************************************/
1645 : /* GDALRPCOpenDEM() */
1646 : /************************************************************************/
1647 :
1648 33 : static bool GDALRPCOpenDEM(GDALRPCTransformInfo *psTransform)
1649 : {
1650 33 : CPLAssert(psTransform->pszDEMPath != nullptr);
1651 :
1652 33 : bool bIsValid = false;
1653 :
1654 66 : CPLString osPrevValueConfigOption;
1655 33 : if (psTransform->bApplyDEMVDatumShift)
1656 : {
1657 : osPrevValueConfigOption =
1658 32 : CPLGetThreadLocalConfigOption("GTIFF_REPORT_COMPD_CS", "");
1659 32 : CPLSetThreadLocalConfigOption("GTIFF_REPORT_COMPD_CS", "YES");
1660 : }
1661 33 : CPLConfigOptionSetter oSetter("CPL_ALLOW_VSISTDIN", "NO", true);
1662 33 : psTransform->poDS =
1663 33 : GDALDataset::FromHandle(GDALOpen(psTransform->pszDEMPath, GA_ReadOnly));
1664 65 : if (psTransform->poDS != nullptr &&
1665 32 : psTransform->poDS->GetRasterCount() >= 1)
1666 : {
1667 64 : OGRSpatialReference oDEMSRS;
1668 32 : if (psTransform->pszDEMSRS != nullptr)
1669 : {
1670 1 : oDEMSRS.SetFromUserInput(psTransform->pszDEMSRS);
1671 1 : oDEMSRS.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1672 : }
1673 :
1674 32 : auto poDSSpaRefSrc = psTransform->pszDEMSRS != nullptr
1675 32 : ? &oDEMSRS
1676 31 : : psTransform->poDS->GetSpatialRef();
1677 32 : if (poDSSpaRefSrc)
1678 : {
1679 32 : auto poDSSpaRef = poDSSpaRefSrc->Clone();
1680 :
1681 32 : if (!psTransform->bApplyDEMVDatumShift)
1682 1 : poDSSpaRef->StripVertical();
1683 :
1684 32 : auto wkt_EPSG_4979 =
1685 : "GEODCRS[\"WGS 84\",\n"
1686 : " DATUM[\"World Geodetic System 1984\",\n"
1687 : " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n"
1688 : " LENGTHUNIT[\"metre\",1]]],\n"
1689 : " PRIMEM[\"Greenwich\",0,\n"
1690 : " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
1691 : " CS[ellipsoidal,3],\n"
1692 : " AXIS[\"geodetic latitude (Lat)\",north,\n"
1693 : " ORDER[1],\n"
1694 : " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
1695 : " AXIS[\"geodetic longitude (Lon)\",east,\n"
1696 : " ORDER[2],\n"
1697 : " ANGLEUNIT[\"degree\",0.0174532925199433]],\n"
1698 : " AXIS[\"ellipsoidal height (h)\",up,\n"
1699 : " ORDER[3],\n"
1700 : " LENGTHUNIT[\"metre\",1]],\n"
1701 : " AREA[\"World (by country)\"],\n"
1702 : " BBOX[-90,-180,90,180],\n"
1703 : " ID[\"EPSG\",4979]]";
1704 : OGRSpatialReference *poWGSSpaRef = new OGRSpatialReference(
1705 32 : poDSSpaRef->IsCompound() ? wkt_EPSG_4979
1706 32 : : SRS_WKT_WGS84_LAT_LONG);
1707 32 : poWGSSpaRef->SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1708 :
1709 32 : if (!poWGSSpaRef->IsSame(poDSSpaRef))
1710 12 : psTransform->poCT =
1711 12 : OGRCreateCoordinateTransformation(poWGSSpaRef, poDSSpaRef);
1712 :
1713 32 : if (psTransform->poCT != nullptr && !poDSSpaRef->IsCompound())
1714 : {
1715 : // Empiric attempt to guess if the coordinate transformation
1716 : // to WGS84 is a no-op. For example for NED13 datasets in
1717 : // NAD83.
1718 10 : double adfX[] = {-179.0, 179.0, 179.0, -179.0, 0.0, 0.0};
1719 10 : double adfY[] = {89.0, 89.0, -89.0, -89.0, 0.0, 0.0};
1720 10 : double adfZ[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1721 :
1722 : // Also test with a "reference point" from the RPC values.
1723 10 : double dfRefLong = 0.0;
1724 10 : double dfRefLat = 0.0;
1725 10 : if (psTransform->sRPC.dfMIN_LONG != -180 ||
1726 10 : psTransform->sRPC.dfMAX_LONG != 180)
1727 : {
1728 0 : dfRefLong = (psTransform->sRPC.dfMIN_LONG +
1729 0 : psTransform->sRPC.dfMAX_LONG) *
1730 : 0.5;
1731 0 : dfRefLat = (psTransform->sRPC.dfMIN_LAT +
1732 0 : psTransform->sRPC.dfMAX_LAT) *
1733 : 0.5;
1734 : }
1735 : else
1736 : {
1737 10 : dfRefLong = psTransform->sRPC.dfLONG_OFF;
1738 10 : dfRefLat = psTransform->sRPC.dfLAT_OFF;
1739 : }
1740 10 : adfX[5] = dfRefLong;
1741 10 : adfY[5] = dfRefLat;
1742 :
1743 10 : if (psTransform->poCT->Transform(6, adfX, adfY, adfZ) &&
1744 10 : fabs(adfX[0] - -179.0) < 1.0e-12 &&
1745 1 : fabs(adfY[0] - 89.0) < 1.0e-12 &&
1746 1 : fabs(adfX[1] - 179.0) < 1.0e-12 &&
1747 1 : fabs(adfY[1] - 89.0) < 1.0e-12 &&
1748 1 : fabs(adfX[2] - 179.0) < 1.0e-12 &&
1749 1 : fabs(adfY[2] - -89.0) < 1.0e-12 &&
1750 1 : fabs(adfX[3] - -179.0) < 1.0e-12 &&
1751 1 : fabs(adfY[3] - -89.0) < 1.0e-12 &&
1752 1 : fabs(adfX[4] - 0.0) < 1.0e-12 &&
1753 1 : fabs(adfY[4] - 0.0) < 1.0e-12 &&
1754 21 : fabs(adfX[5] - dfRefLong) < 1.0e-12 &&
1755 1 : fabs(adfY[5] - dfRefLat) < 1.0e-12)
1756 : {
1757 1 : CPLDebug("RPC",
1758 : "Short-circuiting coordinate transformation "
1759 : "from DEM SRS to WGS 84 due to apparent nop");
1760 1 : delete psTransform->poCT;
1761 1 : psTransform->poCT = nullptr;
1762 : }
1763 : }
1764 :
1765 32 : delete poWGSSpaRef;
1766 32 : delete poDSSpaRef;
1767 : }
1768 :
1769 96 : if (psTransform->poDS->GetGeoTransform(
1770 64 : psTransform->adfDEMGeoTransform) == CE_None &&
1771 32 : GDALInvGeoTransform(psTransform->adfDEMGeoTransform,
1772 32 : psTransform->adfDEMReverseGeoTransform))
1773 : {
1774 32 : bIsValid = true;
1775 : }
1776 : }
1777 :
1778 33 : if (psTransform->bApplyDEMVDatumShift)
1779 : {
1780 32 : CPLSetThreadLocalConfigOption("GTIFF_REPORT_COMPD_CS",
1781 32 : !osPrevValueConfigOption.empty()
1782 0 : ? osPrevValueConfigOption.c_str()
1783 : : nullptr);
1784 : }
1785 :
1786 66 : return bIsValid;
1787 : }
1788 :
1789 : /************************************************************************/
1790 : /* GDALRPCTransform() */
1791 : /************************************************************************/
1792 :
1793 : /** RPC transform */
1794 9386 : int GDALRPCTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
1795 : double *padfX, double *padfY, double *padfZ,
1796 : int *panSuccess)
1797 :
1798 : {
1799 9386 : VALIDATE_POINTER1(pTransformArg, "GDALRPCTransform", 0);
1800 :
1801 9386 : GDALRPCTransformInfo *psTransform =
1802 : static_cast<GDALRPCTransformInfo *>(pTransformArg);
1803 :
1804 9386 : if (psTransform->bReversed)
1805 0 : bDstToSrc = !bDstToSrc;
1806 :
1807 : /* -------------------------------------------------------------------- */
1808 : /* The simple case is transforming from lat/long to pixel/line. */
1809 : /* Just apply the equations directly. */
1810 : /* -------------------------------------------------------------------- */
1811 9386 : if (bDstToSrc)
1812 : {
1813 : // Optimization to avoid doing too many picking in DEM in the particular
1814 : // case where each point to transform is on a single line of the DEM.
1815 : // To make it simple and fast we check that all input latitudes are
1816 : // identical, that the DEM is in WGS84 geodetic and that it has no
1817 : // rotation. Such case is for example triggered when doing gdalwarp
1818 : // with a target SRS of EPSG:4326 or EPSG:3857.
1819 3473 : if (nPointCount >= 10 && psTransform->poDS != nullptr &&
1820 3446 : psTransform->poCT == nullptr &&
1821 2884 : padfY[0] == padfY[nPointCount - 1] &&
1822 2826 : padfY[0] == padfY[nPointCount / 2] &&
1823 2822 : psTransform->adfDEMReverseGeoTransform[1] > 0.0 &&
1824 2822 : psTransform->adfDEMReverseGeoTransform[2] == 0.0 &&
1825 13948 : psTransform->adfDEMReverseGeoTransform[4] == 0.0 &&
1826 2822 : CPLTestBool(CPLGetConfigOption("GDAL_RPC_DEM_OPTIM", "YES")))
1827 : {
1828 2822 : bool bUseOptimized = true;
1829 2822 : double dfMinX = padfX[0];
1830 2822 : double dfMaxX = padfX[0];
1831 327074 : for (int i = 1; i < nPointCount; i++)
1832 : {
1833 324252 : if (padfY[i] != padfY[0])
1834 : {
1835 0 : bUseOptimized = false;
1836 0 : break;
1837 : }
1838 324252 : if (padfX[i] < dfMinX)
1839 0 : dfMinX = padfX[i];
1840 324252 : if (padfX[i] > dfMaxX)
1841 324225 : dfMaxX = padfX[i];
1842 : }
1843 2822 : if (bUseOptimized)
1844 : {
1845 2822 : double dfX1 = 0.0;
1846 2822 : double dfY1 = 0.0;
1847 2822 : double dfX2 = 0.0;
1848 2822 : double dfY2 = 0.0;
1849 2822 : GDALApplyGeoTransform(psTransform->adfDEMReverseGeoTransform,
1850 : dfMinX, padfY[0], &dfX1, &dfY1);
1851 2822 : GDALApplyGeoTransform(psTransform->adfDEMReverseGeoTransform,
1852 : dfMaxX, padfY[0], &dfX2, &dfY2);
1853 :
1854 : // Convert to center of pixel convention for reading the image
1855 : // data.
1856 2822 : if (psTransform->eResampleAlg != DRA_NearestNeighbour)
1857 : {
1858 2821 : dfX1 -= 0.5;
1859 2821 : dfY1 -= 0.5;
1860 2821 : dfX2 -= 0.5;
1861 : // dfY2 -= 0.5;
1862 : }
1863 2822 : int nXLeft = static_cast<int>(floor(dfX1));
1864 2822 : int nXRight = static_cast<int>(floor(dfX2));
1865 2822 : int nXWidth = nXRight - nXLeft + 1;
1866 2822 : int nYTop = static_cast<int>(floor(dfY1));
1867 : int nYHeight;
1868 2822 : if (psTransform->eResampleAlg == DRA_CubicSpline)
1869 : {
1870 1 : nXLeft--;
1871 1 : nXWidth += 3;
1872 1 : nYTop--;
1873 1 : nYHeight = 4;
1874 : }
1875 2821 : else if (psTransform->eResampleAlg == DRA_Bilinear)
1876 : {
1877 2820 : nXWidth++;
1878 2820 : nYHeight = 2;
1879 : }
1880 : else
1881 : {
1882 1 : nYHeight = 1;
1883 : }
1884 5022 : if (nXLeft >= 0 &&
1885 2200 : nXLeft + nXWidth <= psTransform->poDS->GetRasterXSize() &&
1886 5022 : nYTop >= 0 &&
1887 2200 : nYTop + nYHeight <= psTransform->poDS->GetRasterYSize())
1888 : {
1889 : static bool bOnce = false;
1890 2192 : if (!bOnce)
1891 : {
1892 2 : bOnce = true;
1893 2 : CPLDebug("RPC",
1894 : "Using GDALRPCTransformWholeLineWithDEM");
1895 : }
1896 2192 : return GDALRPCTransformWholeLineWithDEM(
1897 : psTransform, nPointCount, padfX, padfY, padfZ,
1898 2192 : panSuccess, nXLeft, nXWidth, nYTop, nYHeight);
1899 : }
1900 : }
1901 : }
1902 :
1903 371140 : for (int i = 0; i < nPointCount; i++)
1904 : {
1905 365679 : if (!RPCIsValidLongLat(psTransform, padfX[i], padfY[i]))
1906 : {
1907 108994 : panSuccess[i] = FALSE;
1908 108994 : padfX[i] = HUGE_VAL;
1909 108994 : padfY[i] = HUGE_VAL;
1910 121421 : continue;
1911 : }
1912 256685 : double dfHeight = 0.0;
1913 256685 : if (!GDALRPCGetHeightAtLongLat(psTransform, padfX[i], padfY[i],
1914 : &dfHeight))
1915 : {
1916 12427 : panSuccess[i] = FALSE;
1917 12427 : padfX[i] = HUGE_VAL;
1918 12427 : padfY[i] = HUGE_VAL;
1919 12427 : continue;
1920 : }
1921 :
1922 244258 : RPCTransformPoint(psTransform, padfX[i], padfY[i],
1923 244258 : (padfZ ? padfZ[i] : 0.0) + dfHeight, padfX + i,
1924 244258 : padfY + i);
1925 244258 : panSuccess[i] = TRUE;
1926 : }
1927 :
1928 5461 : return TRUE;
1929 : }
1930 :
1931 1733 : if (padfZ == nullptr)
1932 : {
1933 0 : CPLError(CE_Failure, CPLE_NotSupported,
1934 : "Z array should be provided for reverse RPC computation");
1935 0 : return FALSE;
1936 : }
1937 :
1938 : /* -------------------------------------------------------------------- */
1939 : /* Compute the inverse (pixel/line/height to lat/long). This */
1940 : /* function uses an iterative method from an initial linear */
1941 : /* approximation. */
1942 : /* -------------------------------------------------------------------- */
1943 15347 : for (int i = 0; i < nPointCount; i++)
1944 : {
1945 13614 : double dfResultX = 0.0;
1946 13614 : double dfResultY = 0.0;
1947 :
1948 13614 : if (!RPCInverseTransformPoint(psTransform, padfX[i], padfY[i], padfZ[i],
1949 : &dfResultX, &dfResultY))
1950 : {
1951 10602 : panSuccess[i] = FALSE;
1952 10602 : padfX[i] = HUGE_VAL;
1953 10602 : padfY[i] = HUGE_VAL;
1954 10602 : continue;
1955 : }
1956 3012 : if (!RPCIsValidLongLat(psTransform, padfX[i], padfY[i]))
1957 : {
1958 0 : panSuccess[i] = FALSE;
1959 0 : padfX[i] = HUGE_VAL;
1960 0 : padfY[i] = HUGE_VAL;
1961 0 : continue;
1962 : }
1963 :
1964 3012 : padfX[i] = dfResultX;
1965 3012 : padfY[i] = dfResultY;
1966 :
1967 3012 : panSuccess[i] = TRUE;
1968 : }
1969 :
1970 1733 : return TRUE;
1971 : }
1972 :
1973 : /************************************************************************/
1974 : /* GDALSerializeRPCTransformer() */
1975 : /************************************************************************/
1976 :
1977 1 : CPLXMLNode *GDALSerializeRPCTransformer(void *pTransformArg)
1978 :
1979 : {
1980 1 : VALIDATE_POINTER1(pTransformArg, "GDALSerializeRPCTransformer", nullptr);
1981 :
1982 1 : GDALRPCTransformInfo *psInfo =
1983 : static_cast<GDALRPCTransformInfo *>(pTransformArg);
1984 :
1985 : CPLXMLNode *psTree =
1986 1 : CPLCreateXMLNode(nullptr, CXT_Element, "RPCTransformer");
1987 :
1988 : /* -------------------------------------------------------------------- */
1989 : /* Serialize bReversed. */
1990 : /* -------------------------------------------------------------------- */
1991 1 : CPLCreateXMLElementAndValue(psTree, "Reversed",
1992 2 : CPLString().Printf("%d", psInfo->bReversed));
1993 :
1994 : /* -------------------------------------------------------------------- */
1995 : /* Serialize Height Offset. */
1996 : /* -------------------------------------------------------------------- */
1997 1 : CPLCreateXMLElementAndValue(
1998 : psTree, "HeightOffset",
1999 2 : CPLString().Printf("%.15g", psInfo->dfHeightOffset));
2000 :
2001 : /* -------------------------------------------------------------------- */
2002 : /* Serialize Height Scale. */
2003 : /* -------------------------------------------------------------------- */
2004 1 : if (psInfo->dfHeightScale != 1.0)
2005 0 : CPLCreateXMLElementAndValue(
2006 : psTree, "HeightScale",
2007 0 : CPLString().Printf("%.15g", psInfo->dfHeightScale));
2008 :
2009 : /* -------------------------------------------------------------------- */
2010 : /* Serialize DEM path. */
2011 : /* -------------------------------------------------------------------- */
2012 1 : if (psInfo->pszDEMPath != nullptr)
2013 : {
2014 0 : CPLCreateXMLElementAndValue(
2015 0 : psTree, "DEMPath", CPLString().Printf("%s", psInfo->pszDEMPath));
2016 :
2017 : /* --------------------------------------------------------------------
2018 : */
2019 : /* Serialize DEM interpolation */
2020 : /* --------------------------------------------------------------------
2021 : */
2022 0 : CPLCreateXMLElementAndValue(
2023 : psTree, "DEMInterpolation",
2024 : GDALSerializeRPCDEMResample(psInfo->eResampleAlg));
2025 :
2026 0 : if (psInfo->bHasDEMMissingValue)
2027 : {
2028 0 : CPLCreateXMLElementAndValue(
2029 : psTree, "DEMMissingValue",
2030 : CPLSPrintf("%.17g", psInfo->dfDEMMissingValue));
2031 : }
2032 :
2033 0 : CPLCreateXMLElementAndValue(psTree, "DEMApplyVDatumShift",
2034 0 : psInfo->bApplyDEMVDatumShift ? "true"
2035 : : "false");
2036 :
2037 : /* --------------------------------------------------------------------
2038 : */
2039 : /* Serialize DEM SRS */
2040 : /* --------------------------------------------------------------------
2041 : */
2042 0 : if (psInfo->pszDEMSRS != nullptr)
2043 : {
2044 0 : CPLCreateXMLElementAndValue(psTree, "DEMSRS", psInfo->pszDEMSRS);
2045 : }
2046 : }
2047 :
2048 : /* -------------------------------------------------------------------- */
2049 : /* Serialize pixel error threshold. */
2050 : /* -------------------------------------------------------------------- */
2051 1 : CPLCreateXMLElementAndValue(
2052 : psTree, "PixErrThreshold",
2053 2 : CPLString().Printf("%.15g", psInfo->dfPixErrThreshold));
2054 :
2055 : /* -------------------------------------------------------------------- */
2056 : /* RPC metadata. */
2057 : /* -------------------------------------------------------------------- */
2058 1 : char **papszMD = RPCInfoV2ToMD(&(psInfo->sRPC));
2059 1 : CPLXMLNode *psMD = CPLCreateXMLNode(psTree, CXT_Element, "Metadata");
2060 :
2061 21 : for (int i = 0; papszMD != nullptr && papszMD[i] != nullptr; i++)
2062 : {
2063 20 : char *pszKey = nullptr;
2064 :
2065 20 : const char *pszRawValue = CPLParseNameValue(papszMD[i], &pszKey);
2066 :
2067 20 : CPLXMLNode *psMDI = CPLCreateXMLNode(psMD, CXT_Element, "MDI");
2068 20 : CPLSetXMLValue(psMDI, "#key", pszKey);
2069 20 : CPLCreateXMLNode(psMDI, CXT_Text, pszRawValue);
2070 :
2071 20 : CPLFree(pszKey);
2072 : }
2073 :
2074 1 : CSLDestroy(papszMD);
2075 :
2076 1 : return psTree;
2077 : }
2078 :
2079 : /************************************************************************/
2080 : /* GDALDeserializeRPCTransformer() */
2081 : /************************************************************************/
2082 :
2083 0 : void *GDALDeserializeRPCTransformer(CPLXMLNode *psTree)
2084 :
2085 : {
2086 0 : char **papszOptions = nullptr;
2087 :
2088 : /* -------------------------------------------------------------------- */
2089 : /* Collect metadata. */
2090 : /* -------------------------------------------------------------------- */
2091 0 : CPLXMLNode *psMetadata = CPLGetXMLNode(psTree, "Metadata");
2092 :
2093 0 : if (psMetadata == nullptr || psMetadata->eType != CXT_Element ||
2094 0 : !EQUAL(psMetadata->pszValue, "Metadata"))
2095 0 : return nullptr;
2096 :
2097 0 : char **papszMD = nullptr;
2098 0 : for (CPLXMLNode *psMDI = psMetadata->psChild; psMDI != nullptr;
2099 0 : psMDI = psMDI->psNext)
2100 : {
2101 0 : if (!EQUAL(psMDI->pszValue, "MDI") || psMDI->eType != CXT_Element ||
2102 0 : psMDI->psChild == nullptr || psMDI->psChild->psNext == nullptr ||
2103 0 : psMDI->psChild->eType != CXT_Attribute ||
2104 0 : psMDI->psChild->psChild == nullptr)
2105 0 : continue;
2106 :
2107 0 : papszMD = CSLSetNameValue(papszMD, psMDI->psChild->psChild->pszValue,
2108 0 : psMDI->psChild->psNext->pszValue);
2109 : }
2110 :
2111 : GDALRPCInfoV2 sRPC;
2112 0 : if (!GDALExtractRPCInfoV2(papszMD, &sRPC))
2113 : {
2114 0 : CPLError(CE_Failure, CPLE_AppDefined,
2115 : "Failed to reconstitute RPC transformer.");
2116 0 : CSLDestroy(papszMD);
2117 0 : return nullptr;
2118 : }
2119 :
2120 0 : CSLDestroy(papszMD);
2121 :
2122 : /* -------------------------------------------------------------------- */
2123 : /* Get other flags. */
2124 : /* -------------------------------------------------------------------- */
2125 0 : const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
2126 :
2127 : const double dfPixErrThreshold =
2128 0 : CPLAtof(CPLGetXMLValue(psTree, "PixErrThreshold",
2129 : CPLSPrintf("%f", DEFAULT_PIX_ERR_THRESHOLD)));
2130 :
2131 0 : papszOptions = CSLSetNameValue(papszOptions, "RPC_HEIGHT",
2132 : CPLGetXMLValue(psTree, "HeightOffset", "0"));
2133 0 : papszOptions = CSLSetNameValue(papszOptions, "RPC_HEIGHT_SCALE",
2134 : CPLGetXMLValue(psTree, "HeightScale", "1"));
2135 0 : const char *pszDEMPath = CPLGetXMLValue(psTree, "DEMPath", nullptr);
2136 0 : if (pszDEMPath != nullptr)
2137 0 : papszOptions = CSLSetNameValue(papszOptions, "RPC_DEM", pszDEMPath);
2138 :
2139 : const char *pszDEMInterpolation =
2140 0 : CPLGetXMLValue(psTree, "DEMInterpolation", "bilinear");
2141 0 : if (pszDEMInterpolation != nullptr)
2142 0 : papszOptions = CSLSetNameValue(papszOptions, "RPC_DEMINTERPOLATION",
2143 : pszDEMInterpolation);
2144 :
2145 : const char *pszDEMMissingValue =
2146 0 : CPLGetXMLValue(psTree, "DEMMissingValue", nullptr);
2147 0 : if (pszDEMMissingValue != nullptr)
2148 0 : papszOptions = CSLSetNameValue(papszOptions, "RPC_DEM_MISSING_VALUE",
2149 : pszDEMMissingValue);
2150 :
2151 : const char *pszDEMApplyVDatumShift =
2152 0 : CPLGetXMLValue(psTree, "DEMApplyVDatumShift", nullptr);
2153 0 : if (pszDEMApplyVDatumShift != nullptr)
2154 0 : papszOptions = CSLSetNameValue(
2155 : papszOptions, "RPC_DEM_APPLY_VDATUM_SHIFT", pszDEMApplyVDatumShift);
2156 0 : const char *pszDEMSRS = CPLGetXMLValue(psTree, "DEMSRS", nullptr);
2157 0 : if (pszDEMSRS != nullptr)
2158 0 : papszOptions = CSLSetNameValue(papszOptions, "RPC_DEM_SRS", pszDEMSRS);
2159 :
2160 : /* -------------------------------------------------------------------- */
2161 : /* Generate transformation. */
2162 : /* -------------------------------------------------------------------- */
2163 0 : void *pResult = GDALCreateRPCTransformerV2(&sRPC, bReversed,
2164 : dfPixErrThreshold, papszOptions);
2165 :
2166 0 : CSLDestroy(papszOptions);
2167 :
2168 0 : return pResult;
2169 : }
|
