/******************************************************************************

 *

 * Project:  Image Warper

 * Purpose:  Implements a rational polynomial (RPC) based transformer.

 * Author:   Frank Warmerdam, warmerdam@pobox.com

 *

 ******************************************************************************

 * Copyright (c) 2003, Frank Warmerdam <warmerdam@pobox.com>

 * Copyright (c) 2009-2013, Even Rouault <even dot rouault at spatialys.com>

 *

 * SPDX-License-Identifier: MIT

 ****************************************************************************/

#include "cpl_port.h"

#include "gdal_alg.h"

#include <cmath>

#include <cstddef>

#include <cstdlib>

#include <cstring>

#include <algorithm>

#include <limits>

#include <string>

#include "cpl_conv.h"

#include "cpl_error.h"

#include "cpl_mem_cache.h"

#include "cpl_minixml.h"

#include "cpl_string.h"

#include "cpl_vsi.h"

#include "gdal.h"

#include "gdal_interpolateatpoint.h"

#include "gdal_mdreader.h"

#include "gdal_alg_priv.h"

#include "gdal_priv.h"

#ifdef USE_NEON_OPTIMIZATIONS

#define USE_SSE2

#elif defined(__x86_64) || defined(_M_X64)

#define USE_SSE2

#endif

#ifdef USE_SSE2

#include "gdalsse_priv.h"

#define USE_SSE2_OPTIM

#endif

#include "ogr_api.h"

#include "ogr_geometry.h"

#include "ogr_spatialref.h"

#include "ogr_srs_api.h"

#include "gdalresamplingkernels.h"

// #define DEBUG_VERBOSE_EXTRACT_DEM

CPL_C_START

CPLXMLNode *GDALSerializeRPCTransformer(void *pTransformArg);

void *GDALDeserializeRPCTransformer(CPLXMLNode *psTree);

CPL_C_END

constexpr int MAX_ABS_VALUE_WARNINGS = 20;

constexpr double DEFAULT_PIX_ERR_THRESHOLD = 0.1;

/************************************************************************/

/*                            RPCInfoToMD()                             */

/*                                                                      */

/*      Turn an RPCInfo structure back into its metadata format.        */

/************************************************************************/

char **RPCInfoV1ToMD(GDALRPCInfoV1 *psRPCInfo)

{

    GDALRPCInfoV2 sRPCInfo;

    memcpy(&sRPCInfo, psRPCInfo, sizeof(GDALRPCInfoV1));

    sRPCInfo.dfERR_BIAS = std::numeric_limits<double>::quiet_NaN();

    sRPCInfo.dfERR_RAND = std::numeric_limits<double>::quiet_NaN();

    return RPCInfoV2ToMD(&sRPCInfo);

}

char **RPCInfoV2ToMD(GDALRPCInfoV2 *psRPCInfo)

{

    char **papszMD = nullptr;

    CPLString osField, osMultiField;

    if (!std::isnan(psRPCInfo->dfERR_BIAS))

    {

        osField.Printf("%.15g", psRPCInfo->dfERR_BIAS);

        papszMD = CSLSetNameValue(papszMD, RPC_ERR_BIAS, osField);

    }

    if (!std::isnan(psRPCInfo->dfERR_RAND))

    {

        osField.Printf("%.15g", psRPCInfo->dfERR_RAND);

        papszMD = CSLSetNameValue(papszMD, RPC_ERR_RAND, osField);

    }

    osField.Printf("%.15g", psRPCInfo->dfLINE_OFF);

    papszMD = CSLSetNameValue(papszMD, RPC_LINE_OFF, osField);

    osField.Printf("%.15g", psRPCInfo->dfSAMP_OFF);

    papszMD = CSLSetNameValue(papszMD, RPC_SAMP_OFF, osField);

    osField.Printf("%.15g", psRPCInfo->dfLAT_OFF);

    papszMD = CSLSetNameValue(papszMD, RPC_LAT_OFF, osField);

    osField.Printf("%.15g", psRPCInfo->dfLONG_OFF);

    papszMD = CSLSetNameValue(papszMD, RPC_LONG_OFF, osField);

    osField.Printf("%.15g", psRPCInfo->dfHEIGHT_OFF);

    papszMD = CSLSetNameValue(papszMD, RPC_HEIGHT_OFF, osField);

    osField.Printf("%.15g", psRPCInfo->dfLINE_SCALE);

    papszMD = CSLSetNameValue(papszMD, RPC_LINE_SCALE, osField);

    osField.Printf("%.15g", psRPCInfo->dfSAMP_SCALE);

    papszMD = CSLSetNameValue(papszMD, RPC_SAMP_SCALE, osField);

    osField.Printf("%.15g", psRPCInfo->dfLAT_SCALE);

    papszMD = CSLSetNameValue(papszMD, RPC_LAT_SCALE, osField);

    osField.Printf("%.15g", psRPCInfo->dfLONG_SCALE);

    papszMD = CSLSetNameValue(papszMD, RPC_LONG_SCALE, osField);

    osField.Printf("%.15g", psRPCInfo->dfHEIGHT_SCALE);

    papszMD = CSLSetNameValue(papszMD, RPC_HEIGHT_SCALE, osField);

    osField.Printf("%.15g", psRPCInfo->dfMIN_LONG);

    papszMD = CSLSetNameValue(papszMD, RPC_MIN_LONG, osField);

    osField.Printf("%.15g", psRPCInfo->dfMIN_LAT);

    papszMD = CSLSetNameValue(papszMD, RPC_MIN_LAT, osField);

    osField.Printf("%.15g", psRPCInfo->dfMAX_LONG);

    papszMD = CSLSetNameValue(papszMD, RPC_MAX_LONG, osField);

    osField.Printf("%.15g", psRPCInfo->dfMAX_LAT);

    papszMD = CSLSetNameValue(papszMD, RPC_MAX_LAT, osField);

    for (int i = 0; i < 20; i++)

    {

        osField.Printf("%.15g", psRPCInfo->adfLINE_NUM_COEFF[i]);

        if (i > 0)

            osMultiField += " ";

        else

            osMultiField = "";

        osMultiField += osField;

    }

    papszMD = CSLSetNameValue(papszMD, "LINE_NUM_COEFF", osMultiField);

    for (int i = 0; i < 20; i++)

    {

        osField.Printf("%.15g", psRPCInfo->adfLINE_DEN_COEFF[i]);

        if (i > 0)

            osMultiField += " ";

        else

            osMultiField = "";

        osMultiField += osField;

    }

    papszMD = CSLSetNameValue(papszMD, "LINE_DEN_COEFF", osMultiField);

    for (int i = 0; i < 20; i++)

    {

        osField.Printf("%.15g", psRPCInfo->adfSAMP_NUM_COEFF[i]);

        if (i > 0)

            osMultiField += " ";

        else

            osMultiField = "";

        osMultiField += osField;

    }

    papszMD = CSLSetNameValue(papszMD, "SAMP_NUM_COEFF", osMultiField);

    for (int i = 0; i < 20; i++)

    {

        osField.Printf("%.15g", psRPCInfo->adfSAMP_DEN_COEFF[i]);

        if (i > 0)

            osMultiField += " ";

        else

            osMultiField = "";

        osMultiField += osField;

    }

    papszMD = CSLSetNameValue(papszMD, "SAMP_DEN_COEFF", osMultiField);

    return papszMD;

}

/************************************************************************/

/*                          RPCComputeTerms()                           */

/************************************************************************/

static void RPCComputeTerms(double dfLong, double dfLat, double dfHeight,

                            double *padfTerms)

{

    padfTerms[0] = 1.0;

    padfTerms[1] = dfLong;

    padfTerms[2] = dfLat;

    padfTerms[3] = dfHeight;

    padfTerms[4] = dfLong * dfLat;

    padfTerms[5] = dfLong * dfHeight;

    padfTerms[6] = dfLat * dfHeight;

    padfTerms[7] = dfLong * dfLong;

    padfTerms[8] = dfLat * dfLat;

    padfTerms[9] = dfHeight * dfHeight;

    padfTerms[10] = dfLong * dfLat * dfHeight;

    padfTerms[11] = dfLong * dfLong * dfLong;

    padfTerms[12] = dfLong * dfLat * dfLat;

    padfTerms[13] = dfLong * dfHeight * dfHeight;

    padfTerms[14] = dfLong * dfLong * dfLat;

    padfTerms[15] = dfLat * dfLat * dfLat;

    padfTerms[16] = dfLat * dfHeight * dfHeight;

    padfTerms[17] = dfLong * dfLong * dfHeight;

    padfTerms[18] = dfLat * dfLat * dfHeight;

    padfTerms[19] = dfHeight * dfHeight * dfHeight;

}

/************************************************************************/

/* ==================================================================== */

/*                           GDALRPCTransformer                         */

/* ==================================================================== */

/************************************************************************/

/*! DEM Resampling Algorithm */

typedef enum

{

    /*! Nearest neighbour (select on one input pixel) */ DRA_NearestNeighbour =

        0,

    /*! Bilinear (2x2 kernel) */ DRA_Bilinear = 1,

    /*! Cubic Convolution Approximation (4x4 kernel) */ DRA_CubicSpline = 2

} DEMResampleAlg;

typedef struct

{

    GDALTransformerInfo sTI;

    GDALRPCInfoV2 sRPC;

    double adfPLToLatLongGeoTransform[6];

    double dfRefZ;

    int bReversed;

    double dfPixErrThreshold;

    double dfHeightOffset;

    double dfHeightScale;

    char *pszDEMPath;

    DEMResampleAlg eResampleAlg;

    int bHasDEMMissingValue;

    double dfDEMMissingValue;

    char *pszDEMSRS;

    int bApplyDEMVDatumShift;

    GDALDataset *poDS;

    // the key is (nYBlock << 32) | nXBlock)

    lru11::Cache<uint64_t, std::shared_ptr<std::vector<double>>> *poCacheDEM;

    OGRCoordinateTransformation *poCT;

    int nMaxIterations;

    double adfDEMGeoTransform[6];

    double adfDEMReverseGeoTransform[6];

#ifdef USE_SSE2_OPTIM

    double adfDoubles[20 * 4 + 1];

    // LINE_NUM_COEFF, LINE_DEN_COEFF, SAMP_NUM_COEFF and then SAMP_DEN_COEFF.

    double *padfCoeffs;

#endif

    bool bRPCInverseVerbose;

    char *pszRPCInverseLog;

    char *pszRPCFootprint;

    OGRGeometry *poRPCFootprintGeom;

    OGRPreparedGeometry *poRPCFootprintPreparedGeom;

} GDALRPCTransformInfo;

static bool GDALRPCOpenDEM(GDALRPCTransformInfo *psTransform);

/************************************************************************/

/*                            RPCEvaluate()                             */

/************************************************************************/

#ifdef USE_SSE2_OPTIM

static void RPCEvaluate4(const double *padfTerms, const double *padfCoefs,

                         double &dfSum1, double &dfSum2, double &dfSum3,

                         double &dfSum4)

{

    XMMReg2Double sum1 = XMMReg2Double::Zero();

    XMMReg2Double sum2 = XMMReg2Double::Zero();

    XMMReg2Double sum3 = XMMReg2Double::Zero();

    XMMReg2Double sum4 = XMMReg2Double::Zero();

    for (int i = 0; i < 20; i += 2)

    {

        const XMMReg2Double terms =

            XMMReg2Double::Load2ValAligned(padfTerms + i);

        // LINE_NUM_COEFF.

        const XMMReg2Double coefs1 =

            XMMReg2Double::Load2ValAligned(padfCoefs + i);

        // LINE_DEN_COEFF.

        const XMMReg2Double coefs2 =

            XMMReg2Double::Load2ValAligned(padfCoefs + i + 20);

        // SAMP_NUM_COEFF.

        const XMMReg2Double coefs3 =

            XMMReg2Double::Load2ValAligned(padfCoefs + i + 40);

        // SAMP_DEN_COEFF.

        const XMMReg2Double coefs4 =

            XMMReg2Double::Load2ValAligned(padfCoefs + i + 60);

        sum1 += terms * coefs1;

        sum2 += terms * coefs2;

        sum3 += terms * coefs3;

        sum4 += terms * coefs4;

    }

    dfSum1 = sum1.GetHorizSum();

    dfSum2 = sum2.GetHorizSum();

    dfSum3 = sum3.GetHorizSum();

    dfSum4 = sum4.GetHorizSum();

}

#else

static double RPCEvaluate(const double *padfTerms, const double *padfCoefs)

{

    double dfSum1 = 0.0;

    double dfSum2 = 0.0;

    for (int i = 0; i < 20; i += 2)

    {

        dfSum1 += padfTerms[i] * padfCoefs[i];

        dfSum2 += padfTerms[i + 1] * padfCoefs[i + 1];

    }

    return dfSum1 + dfSum2;

}

#endif

/************************************************************************/

/*                         RPCTransformPoint()                          */

/************************************************************************/

static void RPCTransformPoint(const GDALRPCTransformInfo *psRPCTransformInfo,

                              double dfLong, double dfLat, double dfHeight,

                              double *pdfPixel, double *pdfLine)

{

    double adfTermsWithMargin[20 + 1] = {};

    // Make padfTerms aligned on 16-byte boundary for SSE2 aligned loads.

    double *padfTerms =

        adfTermsWithMargin +

        (reinterpret_cast<GUIntptr_t>(adfTermsWithMargin) % 16) / 8;

    // Avoid dateline issues.

    double diffLong = dfLong - psRPCTransformInfo->sRPC.dfLONG_OFF;

    if (diffLong < -270)

    {

        diffLong += 360;

    }

    else if (diffLong > 270)

    {

        diffLong -= 360;

    }

    const double dfNormalizedLong =

        diffLong / psRPCTransformInfo->sRPC.dfLONG_SCALE;

    const double dfNormalizedLat =

        (dfLat - psRPCTransformInfo->sRPC.dfLAT_OFF) /

        psRPCTransformInfo->sRPC.dfLAT_SCALE;

    const double dfNormalizedHeight =

        (dfHeight - psRPCTransformInfo->sRPC.dfHEIGHT_OFF) /

        psRPCTransformInfo->sRPC.dfHEIGHT_SCALE;

    // The absolute values of the 3 above normalized values are supposed to be

    // below 1. Warn (as debug message) if it is not the case. We allow for some

    // margin above 1 (1.5, somewhat arbitrary chosen) before warning.

    static int nCountWarningsAboutAboveOneNormalizedValues = 0;

    if (nCountWarningsAboutAboveOneNormalizedValues < MAX_ABS_VALUE_WARNINGS)

    {

        bool bWarned = false;

        if (fabs(dfNormalizedLong) > 1.5)

        {

            bWarned = true;

            CPLDebug(

                "RPC",

                "Normalized %s for (lon,lat,height)=(%f,%f,%f) is %f, "

                "i.e. with an absolute value of > 1, which may cause numeric "

                "stability problems",

                "longitude", dfLong, dfLat, dfHeight, dfNormalizedLong);

        }

        if (fabs(dfNormalizedLat) > 1.5)

        {

            bWarned = true;

            CPLDebug(

                "RPC",

                "Normalized %s for (lon,lat,height)=(%f,%f,%f) is %f, "

                "ie with an absolute value of > 1, which may cause numeric "

                "stability problems",

                "latitude", dfLong, dfLat, dfHeight, dfNormalizedLat);

        }

        if (fabs(dfNormalizedHeight) > 1.5)

        {

            bWarned = true;

            CPLDebug(

                "RPC",

                "Normalized %s for (lon,lat,height)=(%f,%f,%f) is %f, "

                "i.e. with an absolute value of > 1, which may cause numeric "

                "stability problems",

                "height", dfLong, dfLat, dfHeight, dfNormalizedHeight);

        }

        if (bWarned)

        {

            // Limit the number of warnings.

            nCountWarningsAboutAboveOneNormalizedValues++;

            if (nCountWarningsAboutAboveOneNormalizedValues ==

                MAX_ABS_VALUE_WARNINGS)

            {

                CPLDebug("RPC", "No more such debug warnings will be emitted");

            }

        }

    }

    RPCComputeTerms(dfNormalizedLong, dfNormalizedLat, dfNormalizedHeight,

                    padfTerms);

#ifdef USE_SSE2_OPTIM

    double dfSampNum = 0.0;

    double dfSampDen = 0.0;

    double dfLineNum = 0.0;

    double dfLineDen = 0.0;

    RPCEvaluate4(padfTerms, psRPCTransformInfo->padfCoeffs, dfLineNum,

                 dfLineDen, dfSampNum, dfSampDen);

    const double dfResultX = dfSampNum / dfSampDen;

    const double dfResultY = dfLineNum / dfLineDen;

#else

    const double dfResultX =

        RPCEvaluate(padfTerms, psRPCTransformInfo->sRPC.adfSAMP_NUM_COEFF) /

        RPCEvaluate(padfTerms, psRPCTransformInfo->sRPC.adfSAMP_DEN_COEFF);

    const double dfResultY =

        RPCEvaluate(padfTerms, psRPCTransformInfo->sRPC.adfLINE_NUM_COEFF) /

        RPCEvaluate(padfTerms, psRPCTransformInfo->sRPC.adfLINE_DEN_COEFF);

#endif

    // RPCs are using the center of upper left pixel = 0,0 convention

    // convert to top left corner = 0,0 convention used in GDAL.

    *pdfPixel = dfResultX * psRPCTransformInfo->sRPC.dfSAMP_SCALE +

                psRPCTransformInfo->sRPC.dfSAMP_OFF + 0.5;

    *pdfLine = dfResultY * psRPCTransformInfo->sRPC.dfLINE_SCALE +

               psRPCTransformInfo->sRPC.dfLINE_OFF + 0.5;

}

/************************************************************************/

/*                     GDALSerializeRPCDEMResample()                    */

/************************************************************************/

static const char *GDALSerializeRPCDEMResample(DEMResampleAlg eResampleAlg)

{

    switch (eResampleAlg)

    {

        case DRA_NearestNeighbour:

            return "near";

        case DRA_CubicSpline:

            return "cubic";

        default:

        case DRA_Bilinear:

            return "bilinear";

    }

}

/************************************************************************/

/*                   GDALCreateSimilarRPCTransformer()                  */

/************************************************************************/

static void *GDALCreateSimilarRPCTransformer(void *hTransformArg,

                                             double dfRatioX, double dfRatioY)

{

    VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarRPCTransformer",

                      nullptr);

    GDALRPCTransformInfo *psInfo =

        static_cast<GDALRPCTransformInfo *>(hTransformArg);

    GDALRPCInfoV2 sRPC;

    memcpy(&sRPC, &(psInfo->sRPC), sizeof(GDALRPCInfoV2));

    if (dfRatioX != 1.0 || dfRatioY != 1.0)

    {

        sRPC.dfLINE_OFF /= dfRatioY;

        sRPC.dfLINE_SCALE /= dfRatioY;

        sRPC.dfSAMP_OFF /= dfRatioX;

        sRPC.dfSAMP_SCALE /= dfRatioX;

    }

    char **papszOptions = nullptr;

    papszOptions = CSLSetNameValue(papszOptions, "RPC_HEIGHT",

                                   CPLSPrintf("%.17g", psInfo->dfHeightOffset));

    papszOptions = CSLSetNameValue(papszOptions, "RPC_HEIGHT_SCALE",

                                   CPLSPrintf("%.17g", psInfo->dfHeightScale));

    if (psInfo->pszDEMPath != nullptr)

    {

        papszOptions =

            CSLSetNameValue(papszOptions, "RPC_DEM", psInfo->pszDEMPath);

        papszOptions =

            CSLSetNameValue(papszOptions, "RPC_DEMINTERPOLATION",

                            GDALSerializeRPCDEMResample(psInfo->eResampleAlg));

        if (psInfo->bHasDEMMissingValue)

            papszOptions =

                CSLSetNameValue(papszOptions, "RPC_DEM_MISSING_VALUE",

                                CPLSPrintf("%.17g", psInfo->dfDEMMissingValue));

        papszOptions =

            CSLSetNameValue(papszOptions, "RPC_DEM_APPLY_VDATUM_SHIFT",

                            (psInfo->bApplyDEMVDatumShift) ? "TRUE" : "FALSE");

    }

    papszOptions = CSLSetNameValue(papszOptions, "RPC_MAX_ITERATIONS",

                                   CPLSPrintf("%d", psInfo->nMaxIterations));

    GDALRPCTransformInfo *psNewInfo =

        static_cast<GDALRPCTransformInfo *>(GDALCreateRPCTransformerV2(

            &sRPC, psInfo->bReversed, psInfo->dfPixErrThreshold, papszOptions));

    CSLDestroy(papszOptions);

    return psNewInfo;

}

/************************************************************************/

/*                      GDALRPCGetHeightAtLongLat()                     */

/************************************************************************/

static int GDALRPCGetDEMHeight(GDALRPCTransformInfo *psTransform,

                               const double dfXIn, const double dfYIn,

                               double *pdfDEMH);

static bool GDALRPCGetHeightAtLongLat(GDALRPCTransformInfo *psTransform,

                                      const double dfXIn, const double dfYIn,

                                      double *pdfHeight,

                                      double *pdfDEMPixel = nullptr,

                                      double *pdfDEMLine = nullptr)

{

    double dfVDatumShift = 0.0;

    double dfDEMH = 0.0;

    if (psTransform->poDS)

    {

        double dfX = 0.0;

        double dfY = 0.0;

        double dfXTemp = dfXIn;

        double dfYTemp = dfYIn;

        // Check if dem is not in WGS84 and transform points padfX[i], padfY[i].

        if (psTransform->poCT)

        {

            double dfZ = 0.0;

            if (!psTransform->poCT->Transform(1, &dfXTemp, &dfYTemp, &dfZ))

            {

                return false;

            }

            // We must take the opposite since poCT transforms from

            // WGS84 to geoid. And we are going to do the reverse:

            // take an elevation over the geoid and transforms it to WGS84.

            if (psTransform->bApplyDEMVDatumShift)

                dfVDatumShift = -dfZ;

        }

        bool bRetried = false;

    retry:

        GDALApplyGeoTransform(psTransform->adfDEMReverseGeoTransform, dfXTemp,

                              dfYTemp, &dfX, &dfY);

        if (pdfDEMPixel)

            *pdfDEMPixel = dfX;

        if (pdfDEMLine)

            *pdfDEMLine = dfY;

        if (!GDALRPCGetDEMHeight(psTransform, dfX, dfY, &dfDEMH))

        {

            // Try to handle the case where the DEM is in LL WGS84 and spans

            // over [-180,180], (or very close to it ), presumably with much

            // hole in the middle if using VRT, and the longitude goes beyond

            // that interval.

            if (!bRetried && psTransform->poCT == nullptr &&

                (dfXIn >= 180.0 || dfXIn <= -180.0))

            {

                const int nRasterXSize = psTransform->poDS->GetRasterXSize();

                const double dfMinDEMLong = psTransform->adfDEMGeoTransform[0];

                const double dfMaxDEMLong =

                    psTransform->adfDEMGeoTransform[0] +

                    nRasterXSize * psTransform->adfDEMGeoTransform[1];

                if (fabs(dfMinDEMLong - -180) < 0.1 &&

                    fabs(dfMaxDEMLong - 180) < 0.1)

                {

                    if (dfXIn >= 180)

                    {

                        dfXTemp = dfXIn - 360;

                        dfYTemp = dfYIn;

                    }

                    else

                    {

                        dfXTemp = dfXIn + 360;

                        dfYTemp = dfYIn;

                    }

                    bRetried = true;

                    goto retry;

                }

            }

            if (psTransform->bHasDEMMissingValue)

                dfDEMH = psTransform->dfDEMMissingValue;

            else

            {

                return false;

            }

        }

#ifdef DEBUG_VERBOSE_EXTRACT_DEM

        CPLDebug("RPC_DEM", "X=%f, Y=%f -> Z=%f", dfX, dfY, dfDEMH);

#endif

    }

    *pdfHeight = dfVDatumShift + (psTransform->dfHeightOffset +

                                  dfDEMH * psTransform->dfHeightScale);

    return true;

}

/************************************************************************/

/*                      GDALCreateRPCTransformer()                      */

/************************************************************************/

void *GDALCreateRPCTransformerV1(GDALRPCInfoV1 *psRPCInfo, int bReversed,

                                 double dfPixErrThreshold, char **papszOptions)

{

    GDALRPCInfoV2 sRPCInfo;

    memcpy(&sRPCInfo, psRPCInfo, sizeof(GDALRPCInfoV1));

    sRPCInfo.dfERR_BIAS = std::numeric_limits<double>::quiet_NaN();

    sRPCInfo.dfERR_RAND = std::numeric_limits<double>::quiet_NaN();

    return GDALCreateRPCTransformerV2(&sRPCInfo, bReversed, dfPixErrThreshold,

                                      papszOptions);

}

/**

 * Create an RPC based transformer.

 *

 * The geometric sensor model describing the physical relationship between

 * image coordinates and ground coordinates is known as a Rigorous Projection

 * Model. A Rigorous Projection Model expresses the mapping of the image space

 * coordinates of rows and columns (r,c) onto the object space reference

 * surface geodetic coordinates (long, lat, height).

 *

 * A RPC supports a generic description of the Rigorous Projection Models. The

 * approximation used by GDAL (RPC00) is a set of rational polynomials

 * expressing the normalized row and column values, (rn , cn), as a function of

 * normalized geodetic latitude, longitude, and height, (P, L, H), given a

 * set of normalized polynomial coefficients (LINE_NUM_COEF_n, LINE_DEN_COEF_n,

 * SAMP_NUM_COEF_n, SAMP_DEN_COEF_n). Normalized values, rather than actual

 * values are used in order to minimize introduction of errors during the

 * calculations. The transformation between row and column values (r,c), and

 * normalized row and column values (rn, cn), and between the geodetic

 * latitude, longitude, and height and normalized geodetic latitude,

 * longitude, and height (P, L, H), is defined by a set of normalizing

 * translations (offsets) and scales that ensure all values are contained in

 * the range -1 to +1.

 *

 * This function creates a GDALTransformFunc compatible transformer

 * for going between image pixel/line and long/lat/height coordinates

 * using RPCs.  The RPCs are provided in a GDALRPCInfo structure which is

 * normally read from metadata using GDALExtractRPCInfo().

 *

 * GDAL RPC Metadata has the following entries (also described in GDAL RFC 22

 * and the GeoTIFF RPC document http://geotiff.maptools.org/rpc_prop.html .

 *

 * <ul>

 * <li>ERR_BIAS: Error - Bias. The RMS bias error in meters per horizontal axis

 * of all points in the image (-1.0 if unknown)

 * <li>ERR_RAND: Error - Random. RMS random error in meters per horizontal axis

 * of each point in the image (-1.0 if unknown)

 * <li>LINE_OFF: Line Offset

 * <li>SAMP_OFF: Sample Offset

 * <li>LAT_OFF: Geodetic Latitude Offset

 * <li>LONG_OFF: Geodetic Longitude Offset

 * <li>HEIGHT_OFF: Geodetic Height Offset

 * <li>LINE_SCALE: Line Scale

 * <li>SAMP_SCALE: Sample Scale

 * <li>LAT_SCALE: Geodetic Latitude Scale

 * <li>LONG_SCALE: Geodetic Longitude Scale

 * <li>HEIGHT_SCALE: Geodetic Height Scale

 * <li>LINE_NUM_COEFF (1-20): Line Numerator Coefficients. Twenty coefficients

 * for the polynomial in the Numerator of the rn equation. (space separated)

 * <li>LINE_DEN_COEFF (1-20): Line Denominator Coefficients. Twenty coefficients

 * for the polynomial in the Denominator of the rn equation. (space separated)

 * <li>SAMP_NUM_COEFF (1-20): Sample Numerator Coefficients. Twenty coefficients

 * for the polynomial in the Numerator of the cn equation. (space separated)

 * <li>SAMP_DEN_COEFF (1-20): Sample Denominator Coefficients. Twenty

 * coefficients for the polynomial in the Denominator of the cn equation. (space

 * separated)

 * </ul>

 *

 * Some drivers (such as DIMAP) may also fill a HEIGHT_DEFAULT item that can be

 * used by GDALCreateGenImgProjTransformer2() to initialize the below RPC_HEIGHT

 * transformer option if none of RPC_HEIGHT and RPC_DEM are specified.

 * Otherwise, if none of RPC_HEIGHT and RPC_DEM are specified as transformer

 * options and if HEIGHT_DEFAULT is no available, a height of 0 will be used.

 *

 * The transformer normally maps from pixel/line/height to long/lat/height space

 * as a forward transformation though in RPC terms that would be considered

 * an inverse transformation (and is solved by iterative approximation using

 * long/lat/height to pixel/line transformations).  The default direction can

 * be reversed by passing bReversed=TRUE.

 *

 * The iterative solution of pixel/line

 * to lat/long/height is currently run for up to 10 iterations or until

 * the apparent error is less than dfPixErrThreshold pixels.  Passing zero

 * will not avoid all error, but will cause the operation to run for the maximum

 * number of iterations.

 *

 * Starting with GDAL 2.1, debugging of the RPC inverse transformer can be done

 * by setting the RPC_INVERSE_VERBOSE configuration option to YES (in which case

 * extra debug information will be displayed in the "RPC" debug category, so

 * requiring CPL_DEBUG to be also set) and/or by setting RPC_INVERSE_LOG to a

 * filename that will contain the content of iterations (this last option only

 * makes sense when debugging point by point, since each time

 * RPCInverseTransformPoint() is called, the file is rewritten).

 *

 * Additional options to the transformer can be supplied in papszOptions.

 *

 * Options:

 *

 * <ul>

 * <li> RPC_HEIGHT: a fixed height offset to be applied to all points passed

 * in.  In this situation the Z passed into the transformation function is

 * assumed to be height above ground, and the RPC_HEIGHT is assumed to be

 * an average height above sea level for ground in the target scene.</li>

 *

 * <li> RPC_HEIGHT_SCALE: a factor used to multiply heights above ground.

 * Useful when elevation offsets of the DEM are not expressed in meters.</li>

 *

 * <li> RPC_DEM: the name of a GDAL dataset (a DEM file typically) used to

 * extract elevation offsets from. In this situation the Z passed into the

 * transformation function is assumed to be height above ground. This option

 * should be used in replacement of RPC_HEIGHT to provide a way of defining

 * a non uniform ground for the target scene</li>

 *

 * <li> RPC_DEMINTERPOLATION: the DEM interpolation ("near", "bilinear" or

 "cubic").

 *      Default is "bilinear".</li>

 *

 * <li> RPC_DEM_MISSING_VALUE: value of DEM height that must be used in case

 * the DEM has nodata value at the sampling point, or if its extent does not

 * cover the requested coordinate. When not specified, missing values will cause

 * a failed transform.</li>

 *

 * <li> RPC_DEM_SRS: (GDAL >= 3.2) WKT SRS, or any string recognized by

 * OGRSpatialReference::SetFromUserInput(), to be used as an override for DEM

 SRS.

 * Useful if DEM SRS does not have an explicit vertical component. </li>

 *

 * <li> RPC_DEM_APPLY_VDATUM_SHIFT: whether the vertical component of a compound

 * SRS for the DEM should be used (when it is present). This is useful so as to

 * be able to transform the "raw" values from the DEM expressed with respect to

 * a geoid to the heights with respect to the WGS84 ellipsoid. When this is

 * enabled, the GTIFF_REPORT_COMPD_CS configuration option will be also set

 * temporarily so as to get the vertical information from GeoTIFF

 * files. Defaults to TRUE. (GDAL >= 2.1.0)</li>

 *

 * <li> RPC_PIXEL_ERROR_THRESHOLD: overrides the dfPixErrThreshold parameter, ie

  the error (measured in pixels) allowed in the

 * iterative solution of pixel/line to lat/long computations (the other way

 * is always exact given the equations).  (GDAL >= 2.1.0)</li>

 *

 * <li> RPC_MAX_ITERATIONS: maximum number of iterations allowed in the

 * iterative solution of pixel/line to lat/long computations. Default value is

 * 10 in the absence of a DEM, or 20 if there is a DEM.  (GDAL >= 2.1.0)</li>

 *

 * <li> RPC_FOOTPRINT: WKT or GeoJSON polygon (in long / lat coordinate space)

 * with a validity footprint for the RPC. Any coordinate transformation that

 * goes from or arrive outside this footprint will be considered invalid. This

 * is useful in situations where the RPC values become highly unstable outside

 * of the area on which they have been computed for, potentially leading to

 * undesirable "echoes" / false positives. This requires GDAL to be built

 against

 * GEOS.</li>

 *

 * </ul>

 *

 * @param psRPCInfo Definition of the RPC parameters.

 *

 * @param bReversed If true "forward" transformation will be lat/long to

 * pixel/line instead of the normal pixel/line to lat/long.

 *

 * @param dfPixErrThreshold the error (measured in pixels) allowed in the

 * iterative solution of pixel/line to lat/long computations (the other way

 * is always exact given the equations). Starting with GDAL 2.1, this may also

 * be set through the RPC_PIXEL_ERROR_THRESHOLD transformer option.

 * If a negative or null value is provided, then this defaults to 0.1 pixel.

 *

 * @param papszOptions Other transformer options (i.e. RPC_HEIGHT=z).

 *

 * @return transformer callback data (deallocate with GDALDestroyTransformer()).

 */

void *GDALCreateRPCTransformerV2(const GDALRPCInfoV2 *psRPCInfo, int bReversed,

                                 double dfPixErrThreshold,

                                 CSLConstList papszOptions)

{

    /* -------------------------------------------------------------------- */

    /*      Initialize core info.                                           */

    /* -------------------------------------------------------------------- */

    GDALRPCTransformInfo *psTransform = static_cast<GDALRPCTransformInfo *>(

        CPLCalloc(sizeof(GDALRPCTransformInfo), 1));

    memcpy(&(psTransform->sRPC), psRPCInfo, sizeof(GDALRPCInfoV2));

    psTransform->bReversed = bReversed;

    const char *pszPixErrThreshold =

        CSLFetchNameValue(papszOptions, "RPC_PIXEL_ERROR_THRESHOLD");

    if (pszPixErrThreshold != nullptr)

        psTransform->dfPixErrThreshold = CPLAtof(pszPixErrThreshold);

    else if (dfPixErrThreshold > 0)

        psTransform->dfPixErrThreshold = dfPixErrThreshold;

    else

        psTransform->dfPixErrThreshold = DEFAULT_PIX_ERR_THRESHOLD;

    psTransform->dfHeightOffset = 0.0;

    psTransform->dfHeightScale = 1.0;

    memcpy(psTransform->sTI.abySignature, GDAL_GTI2_SIGNATURE,

           strlen(GDAL_GTI2_SIGNATURE));

    psTransform->sTI.pszClassName = GDAL_RPC_TRANSFORMER_CLASS_NAME;

    psTransform->sTI.pfnTransform = GDALRPCTransform;

    psTransform->sTI.pfnCleanup = GDALDestroyRPCTransformer;

    psTransform->sTI.pfnSerialize = GDALSerializeRPCTransformer;

    psTransform->sTI.pfnCreateSimilar = GDALCreateSimilarRPCTransformer;

#ifdef USE_SSE2_OPTIM

    // Make sure padfCoeffs is aligned on a 16-byte boundary for SSE2 aligned

    // loads.

    psTransform->padfCoeffs =

        psTransform->adfDoubles +

        (reinterpret_cast<GUIntptr_t>(psTransform->adfDoubles) % 16) / 8;

    memcpy(psTransform->padfCoeffs, psRPCInfo->adfLINE_NUM_COEFF,

           20 * sizeof(double));

    memcpy(psTransform->padfCoeffs + 20, psRPCInfo->adfLINE_DEN_COEFF,

           20 * sizeof(double));

    memcpy(psTransform->padfCoeffs + 40, psRPCInfo->adfSAMP_NUM_COEFF,

           20 * sizeof(double));

    memcpy(psTransform->padfCoeffs + 60, psRPCInfo->adfSAMP_DEN_COEFF,

           20 * sizeof(double));

#endif

    /* -------------------------------------------------------------------- */

    /*      Do we have a "average height" that we want to consider all      */

    /*      elevations to be relative to?                                   */

    /* -------------------------------------------------------------------- */

    const char *pszHeight = CSLFetchNameValue(papszOptions, "RPC_HEIGHT");

    if (pszHeight != nullptr)

        psTransform->dfHeightOffset = CPLAtof(pszHeight);

    /* -------------------------------------------------------------------- */

    /*                       The "height scale"                             */

    /* -------------------------------------------------------------------- */

    const char *pszHeightScale =

        CSLFetchNameValue(papszOptions, "RPC_HEIGHT_SCALE");

    if (pszHeightScale != nullptr)

        psTransform->dfHeightScale = CPLAtof(pszHeightScale);

    /* -------------------------------------------------------------------- */

    /*                       The DEM file name                              */

    /* -------------------------------------------------------------------- */

    const char *pszDEMPath = CSLFetchNameValue(papszOptions, "RPC_DEM");

    if (pszDEMPath != nullptr)

    {

        psTransform->pszDEMPath = CPLStrdup(pszDEMPath);

    }

    /* -------------------------------------------------------------------- */

    /*                      The DEM interpolation                           */

    /* -------------------------------------------------------------------- */

    const char *pszDEMInterpolation =

        CSLFetchNameValueDef(papszOptions, "RPC_DEMINTERPOLATION", "bilinear");

    if (EQUAL(pszDEMInterpolation, "near"))

    {

        psTransform->eResampleAlg = DRA_NearestNeighbour;

    }

    else if (EQUAL(pszDEMInterpolation, "bilinear"))

    {

        psTransform->eResampleAlg = DRA_Bilinear;

    }

    else if (EQUAL(pszDEMInterpolation, "cubic"))

    {

        psTransform->eResampleAlg = DRA_CubicSpline;

    }

    else

    {

        CPLDebug("RPC", "Unknown interpolation %s. Defaulting to bilinear",

                 pszDEMInterpolation);

        psTransform->eResampleAlg = DRA_Bilinear;

    }

    /* -------------------------------------------------------------------- */

    /*                       The DEM missing value                          */

    /* -------------------------------------------------------------------- */

    const char *pszDEMMissingValue =

        CSLFetchNameValue(papszOptions, "RPC_DEM_MISSING_VALUE");

    if (pszDEMMissingValue != nullptr)

    {

        psTransform->bHasDEMMissingValue = TRUE;

        psTransform->dfDEMMissingValue = CPLAtof(pszDEMMissingValue);

    }

    /* -------------------------------------------------------------------- */

    /*                        The DEM SRS override                          */

    /* -------------------------------------------------------------------- */

    const char *pszDEMSRS = CSLFetchNameValue(papszOptions, "RPC_DEM_SRS");

    if (pszDEMSRS != nullptr)

    {

        psTransform->pszDEMSRS = CPLStrdup(pszDEMSRS);

    }

    /* -------------------------------------------------------------------- */

    /*      Whether to apply vdatum shift                                   */

    /* -------------------------------------------------------------------- */

    psTransform->bApplyDEMVDatumShift =

        CPLFetchBool(papszOptions, "RPC_DEM_APPLY_VDATUM_SHIFT", true);

    psTransform->nMaxIterations =

        atoi(CSLFetchNameValueDef(papszOptions, "RPC_MAX_ITERATIONS", "0"));

    /* -------------------------------------------------------------------- */

    /*      Debug                                                           */

    /* -------------------------------------------------------------------- */

    psTransform->bRPCInverseVerbose =

        CPLTestBool(CPLGetConfigOption("RPC_INVERSE_VERBOSE", "NO"));

    const char *pszRPCInverseLog =

        CPLGetConfigOption("RPC_INVERSE_LOG", nullptr);

    if (pszRPCInverseLog != nullptr)

        psTransform->pszRPCInverseLog = CPLStrdup(pszRPCInverseLog);

    /* -------------------------------------------------------------------- */

    /*      Footprint                                                       */

    /* -------------------------------------------------------------------- */

    const char *pszFootprint = CSLFetchNameValue(papszOptions, "RPC_FOOTPRINT");

    if (pszFootprint != nullptr)

    {

        psTransform->pszRPCFootprint = CPLStrdup(pszFootprint);

        if (pszFootprint[0] == '{')

        {

            psTransform->poRPCFootprintGeom =

                OGRGeometryFactory::createFromGeoJson(pszFootprint);

        }

        else

        {

            OGRGeometryFactory::createFromWkt(

                pszFootprint, nullptr, &(psTransform->poRPCFootprintGeom));

        }