Coverage Report

Created: 2025-06-09 08:44

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