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