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