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