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