LCOV - code coverage report
Current view: top level - alg - gdal_crs.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 401 453 88.5 %
Date: 2025-02-18 14:19:29 Functions: 17 17 100.0 %

          Line data    Source code
       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  Mapinfo Image Warper
       4             :  * Purpose:  Implementation of the GDALTransformer wrapper around CRS.C
       5             :  functions
       6             :  *           to build a polynomial transformation based on ground control
       7             :  *           points.
       8             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       9             :  *
      10             :  ***************************************************************************
      11             : 
      12             :     CRS.C - Center for Remote Sensing rectification routines
      13             : 
      14             :     Written By: Brian J. Buckley
      15             : 
      16             :             At: The Center for Remote Sensing
      17             :                 Michigan State University
      18             :                 302 Berkey Hall
      19             :                 East Lansing, MI  48824
      20             :                 (517)353-7195
      21             : 
      22             :     Written: 12/19/91
      23             : 
      24             :     Last Update: 12/26/91 Brian J. Buckley
      25             :     Last Update:  1/24/92 Brian J. Buckley
      26             :       Added printout of trnfile. Triggered by BDEBUG.
      27             :     Last Update:  1/27/92 Brian J. Buckley
      28             :       Fixed bug so that only the active control points were used.
      29             :     Last Update:  6/29/2011 C. F. Stallmann & R. van den Dool (South African
      30             :  National Space Agency) GCP refinement added
      31             : 
      32             :     Copyright (c) 1992, Michigan State University
      33             :  * Copyright (c) 2008-2013, Even Rouault <even dot rouault at spatialys.com>
      34             : 
      35             :     Permission is hereby granted, free of charge, to any person obtaining a
      36             :     copy of this software and associated documentation files (the "Software"),
      37             :     to deal in the Software without restriction, including without limitation
      38             :     the rights to use, copy, modify, merge, publish, distribute, sublicense,
      39             :     and/or sell copies of the Software, and to permit persons to whom the
      40             :     Software is furnished to do so, subject to the following conditions:
      41             : 
      42             :     The above copyright notice and this permission notice shall be included
      43             :     in all copies or substantial portions of the Software.
      44             : 
      45             :     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
      46             :     IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
      47             :     FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.  IN NO EVENT SHALL
      48             :     THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
      49             :     LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
      50             :     FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
      51             :     DEALINGS IN THE SOFTWARE.
      52             : 
      53             :  ****************************************************************************/
      54             : 
      55             : #include "gdal_alg.h"
      56             : #include "gdal_priv.h"
      57             : #include "cpl_conv.h"
      58             : #include "cpl_minixml.h"
      59             : #include "cpl_string.h"
      60             : #include "cpl_atomic_ops.h"
      61             : 
      62             : #include <math.h>
      63             : #include <stdlib.h>
      64             : #include <string.h>
      65             : 
      66             : #define MAXORDER 3
      67             : 
      68             : namespace
      69             : {
      70             : struct Control_Points
      71             : {
      72             :     int count;
      73             :     double *e1;
      74             :     double *n1;
      75             :     double *e2;
      76             :     double *n2;
      77             :     int *status;
      78             : };
      79             : }  // namespace
      80             : 
      81             : struct GCPTransformInfo
      82             : {
      83             :     GDALTransformerInfo sTI{};
      84             : 
      85             :     double adfToGeoX[20]{};
      86             :     double adfToGeoY[20]{};
      87             : 
      88             :     double adfFromGeoX[20]{};
      89             :     double adfFromGeoY[20]{};
      90             :     double x1_mean{};
      91             :     double y1_mean{};
      92             :     double x2_mean{};
      93             :     double y2_mean{};
      94             :     int nOrder{};
      95             :     int bReversed{};
      96             : 
      97             :     std::vector<gdal::GCP> asGCPs{};
      98             :     int bRefine{};
      99             :     int nMinimumGcps{};
     100             :     double dfTolerance{};
     101             : 
     102             :     volatile int nRefCount{};
     103             : };
     104             : 
     105             : CPL_C_START
     106             : CPLXMLNode *GDALSerializeGCPTransformer(void *pTransformArg);
     107             : void *GDALDeserializeGCPTransformer(CPLXMLNode *psTree);
     108             : CPL_C_END
     109             : 
     110             : /* crs.c */
     111             : static int CRS_georef(double, double, double *, double *, double[], double[],
     112             :                       int);
     113             : static int CRS_compute_georef_equations(GCPTransformInfo *psInfo,
     114             :                                         struct Control_Points *, double[],
     115             :                                         double[], double[], double[], int);
     116             : static int remove_outliers(GCPTransformInfo *);
     117             : 
     118             : #define MSUCCESS 1     /* SUCCESS */
     119             : #define MNPTERR 0      /* NOT ENOUGH POINTS */
     120             : #define MUNSOLVABLE -1 /* NOT SOLVABLE */
     121             : #define MMEMERR -2     /* NOT ENOUGH MEMORY */
     122             : #define MPARMERR -3    /* PARAMETER ERROR */
     123             : #define MINTERR -4     /* INTERNAL ERROR */
     124             : 
     125             : static const char *const CRS_error_message[] = {
     126             :     "Failed to compute GCP transform: Not enough points available",
     127             :     "Failed to compute GCP transform: Transform is not solvable",
     128             :     "Failed to compute GCP transform: Not enough memory",
     129             :     "Failed to compute GCP transform: Parameter error",
     130             :     "Failed to compute GCP transform: Internal error"};
     131             : 
     132             : /************************************************************************/
     133             : /*                   GDALCreateSimilarGCPTransformer()                  */
     134             : /************************************************************************/
     135             : 
     136           2 : static void *GDALCreateSimilarGCPTransformer(void *hTransformArg,
     137             :                                              double dfRatioX, double dfRatioY)
     138             : {
     139           2 :     GCPTransformInfo *psInfo = static_cast<GCPTransformInfo *>(hTransformArg);
     140             : 
     141           2 :     VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarGCPTransformer",
     142             :                       nullptr);
     143             : 
     144           2 :     if (dfRatioX == 1.0 && dfRatioY == 1.0)
     145             :     {
     146             :         /* We can just use a ref count, since using the source transformation */
     147             :         /* is thread-safe */
     148           0 :         CPLAtomicInc(&(psInfo->nRefCount));
     149             :     }
     150             :     else
     151             :     {
     152           2 :         auto newGCPs = psInfo->asGCPs;
     153           8 :         for (auto &gcp : newGCPs)
     154             :         {
     155           6 :             gcp.Pixel() /= dfRatioX;
     156           6 :             gcp.Line() /= dfRatioY;
     157             :         }
     158             :         /* As remove_outliers modifies the provided GCPs we don't need to
     159             :          * reapply it */
     160           4 :         psInfo = static_cast<GCPTransformInfo *>(GDALCreateGCPTransformer(
     161           2 :             static_cast<int>(newGCPs.size()), gdal::GCP::c_ptr(newGCPs),
     162             :             psInfo->nOrder, psInfo->bReversed));
     163             :     }
     164             : 
     165           2 :     return psInfo;
     166             : }
     167             : 
     168             : /************************************************************************/
     169             : /*                      GDALCreateGCPTransformer()                      */
     170             : /************************************************************************/
     171             : 
     172          51 : static void *GDALCreateGCPTransformerEx(int nGCPCount,
     173             :                                         const GDAL_GCP *pasGCPList,
     174             :                                         int nReqOrder, bool bReversed,
     175             :                                         bool bRefine, double dfTolerance,
     176             :                                         int nMinimumGcps)
     177             : 
     178             : {
     179             :     // If no minimumGcp parameter was passed, we  use the default value
     180             :     // according to the model
     181          51 :     if (bRefine && nMinimumGcps == -1)
     182             :     {
     183           0 :         nMinimumGcps = ((nReqOrder + 1) * (nReqOrder + 2)) / 2 + 1;
     184             :     }
     185             : 
     186          51 :     GCPTransformInfo *psInfo = nullptr;
     187          51 :     double *padfGeoX = nullptr;
     188          51 :     double *padfGeoY = nullptr;
     189          51 :     double *padfRasterX = nullptr;
     190          51 :     double *padfRasterY = nullptr;
     191          51 :     int *panStatus = nullptr;
     192          51 :     int iGCP = 0;
     193          51 :     int nCRSresult = 0;
     194             :     struct Control_Points sPoints;
     195             : 
     196          51 :     double x1_sum = 0;
     197          51 :     double y1_sum = 0;
     198          51 :     double x2_sum = 0;
     199          51 :     double y2_sum = 0;
     200             : 
     201          51 :     memset(&sPoints, 0, sizeof(sPoints));
     202             : 
     203          51 :     if (nReqOrder == 0)
     204             :     {
     205          35 :         if (nGCPCount >= 10)
     206           3 :             nReqOrder = 2; /*for now we avoid 3rd order since it is unstable*/
     207          32 :         else if (nGCPCount >= 6)
     208           1 :             nReqOrder = 2;
     209             :         else
     210          31 :             nReqOrder = 1;
     211             :     }
     212             : 
     213          51 :     psInfo = new GCPTransformInfo();
     214          51 :     psInfo->bReversed = bReversed;
     215          51 :     psInfo->nOrder = nReqOrder;
     216          51 :     psInfo->bRefine = bRefine;
     217          51 :     psInfo->dfTolerance = dfTolerance;
     218          51 :     psInfo->nMinimumGcps = nMinimumGcps;
     219             : 
     220          51 :     psInfo->nRefCount = 1;
     221             : 
     222          51 :     psInfo->asGCPs = gdal::GCP::fromC(pasGCPList, nGCPCount);
     223           1 :     if (nGCPCount == 2 && nReqOrder == 1 &&
     224          53 :         psInfo->asGCPs[0].X() != psInfo->asGCPs[1].X() &&
     225           1 :         psInfo->asGCPs[0].Y() != psInfo->asGCPs[1].Y())
     226             :     {
     227             :         // Assumes that the 2 GCPs form opposite corners of a rectangle,
     228             :         // and synthetize a 3rd corner
     229           1 :         gdal::GCP newGCP;
     230           1 :         newGCP.X() = psInfo->asGCPs[1].X();
     231           1 :         newGCP.Y() = psInfo->asGCPs[0].Y();
     232           1 :         newGCP.Pixel() = psInfo->asGCPs[1].Pixel();
     233           1 :         newGCP.Line() = psInfo->asGCPs[0].Line();
     234           1 :         psInfo->asGCPs.emplace_back(std::move(newGCP));
     235             : 
     236           1 :         nGCPCount = 3;
     237           1 :         pasGCPList = gdal::GCP::c_ptr(psInfo->asGCPs);
     238             :     }
     239             : 
     240          51 :     memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
     241             :            strlen(GDAL_GTI2_SIGNATURE));
     242          51 :     psInfo->sTI.pszClassName = "GDALGCPTransformer";
     243          51 :     psInfo->sTI.pfnTransform = GDALGCPTransform;
     244          51 :     psInfo->sTI.pfnCleanup = GDALDestroyGCPTransformer;
     245          51 :     psInfo->sTI.pfnSerialize = GDALSerializeGCPTransformer;
     246          51 :     psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarGCPTransformer;
     247             : 
     248             :     /* -------------------------------------------------------------------- */
     249             :     /*      Compute the forward and reverse polynomials.                    */
     250             :     /* -------------------------------------------------------------------- */
     251             : 
     252          51 :     if (nGCPCount == 0)
     253             :     {
     254           0 :         nCRSresult = MNPTERR;
     255             :     }
     256          51 :     else if (bRefine)
     257             :     {
     258           1 :         nCRSresult = remove_outliers(psInfo);
     259             :     }
     260             :     else
     261             :     {
     262             :         /* --------------------------------------------------------------------
     263             :          */
     264             :         /*      Allocate and initialize the working points list. */
     265             :         /* --------------------------------------------------------------------
     266             :          */
     267             :         try
     268             :         {
     269          50 :             padfGeoX = new double[nGCPCount];
     270          50 :             padfGeoY = new double[nGCPCount];
     271          50 :             padfRasterX = new double[nGCPCount];
     272          50 :             padfRasterY = new double[nGCPCount];
     273          50 :             panStatus = new int[nGCPCount];
     274         857 :             for (iGCP = 0; iGCP < nGCPCount; iGCP++)
     275             :             {
     276         807 :                 panStatus[iGCP] = 1;
     277         807 :                 padfGeoX[iGCP] = pasGCPList[iGCP].dfGCPX;
     278         807 :                 padfGeoY[iGCP] = pasGCPList[iGCP].dfGCPY;
     279         807 :                 padfRasterX[iGCP] = pasGCPList[iGCP].dfGCPPixel;
     280         807 :                 padfRasterY[iGCP] = pasGCPList[iGCP].dfGCPLine;
     281         807 :                 x1_sum += pasGCPList[iGCP].dfGCPPixel;
     282         807 :                 y1_sum += pasGCPList[iGCP].dfGCPLine;
     283         807 :                 x2_sum += pasGCPList[iGCP].dfGCPX;
     284         807 :                 y2_sum += pasGCPList[iGCP].dfGCPY;
     285             :             }
     286          50 :             psInfo->x1_mean = x1_sum / nGCPCount;
     287          50 :             psInfo->y1_mean = y1_sum / nGCPCount;
     288          50 :             psInfo->x2_mean = x2_sum / nGCPCount;
     289          50 :             psInfo->y2_mean = y2_sum / nGCPCount;
     290             : 
     291          50 :             sPoints.count = nGCPCount;
     292          50 :             sPoints.e1 = padfRasterX;
     293          50 :             sPoints.n1 = padfRasterY;
     294          50 :             sPoints.e2 = padfGeoX;
     295          50 :             sPoints.n2 = padfGeoY;
     296          50 :             sPoints.status = panStatus;
     297          50 :             nCRSresult = CRS_compute_georef_equations(
     298          50 :                 psInfo, &sPoints, psInfo->adfToGeoX, psInfo->adfToGeoY,
     299          50 :                 psInfo->adfFromGeoX, psInfo->adfFromGeoY, nReqOrder);
     300             :         }
     301           0 :         catch (const std::exception &e)
     302             :         {
     303           0 :             CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
     304           0 :             nCRSresult = MINTERR;
     305             :         }
     306          50 :         delete[] padfGeoX;
     307          50 :         delete[] padfGeoY;
     308          50 :         delete[] padfRasterX;
     309          50 :         delete[] padfRasterY;
     310          50 :         delete[] panStatus;
     311             :     }
     312             : 
     313          51 :     if (nCRSresult != 1)
     314             :     {
     315           1 :         CPLError(CE_Failure, CPLE_AppDefined, "%s",
     316           1 :                  CRS_error_message[-nCRSresult]);
     317           1 :         GDALDestroyGCPTransformer(psInfo);
     318           1 :         return nullptr;
     319             :     }
     320             :     else
     321             :     {
     322          50 :         return psInfo;
     323             :     }
     324             : }
     325             : 
     326             : /**
     327             :  * Create GCP based polynomial transformer.
     328             :  *
     329             :  * Computes least squares fit polynomials from a provided set of GCPs,
     330             :  * and stores the coefficients for later transformation of points between
     331             :  * pixel/line and georeferenced coordinates.
     332             :  *
     333             :  * The return value should be used as a TransformArg in combination with
     334             :  * the transformation function GDALGCPTransform which fits the
     335             :  * GDALTransformerFunc signature.  The returned transform argument should
     336             :  * be deallocated with GDALDestroyGCPTransformer when no longer needed.
     337             :  *
     338             :  * This function may fail (returning nullptr) if the provided set of GCPs
     339             :  * are inadequate for the requested order, the determinate is zero or they
     340             :  * are otherwise "ill conditioned".
     341             :  *
     342             :  * Note that 2nd order requires at least 6 GCPs, and 3rd order requires at
     343             :  * least 10 gcps.  If nReqOrder is 0 the highest order possible (limited to 2)
     344             :  * with the provided gcp count will be used.
     345             :  *
     346             :  * @param nGCPCount the number of GCPs in pasGCPList.
     347             :  * @param pasGCPList an array of GCPs to be used as input.
     348             :  * @param nReqOrder the requested polynomial order.  It should be 1, 2 or 3.
     349             :  * Using 3 is not recommended due to potential numeric instabilities issues.
     350             :  * @param bReversed set it to TRUE to compute the reversed transformation.
     351             :  *
     352             :  * @return the transform argument or nullptr if creation fails.
     353             :  */
     354          50 : void *GDALCreateGCPTransformer(int nGCPCount, const GDAL_GCP *pasGCPList,
     355             :                                int nReqOrder, int bReversed)
     356             : 
     357             : {
     358          50 :     return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder,
     359         100 :                                       CPL_TO_BOOL(bReversed), false, -1, -1);
     360             : }
     361             : 
     362             : /** Create GCP based polynomial transformer, with a tolerance threshold to
     363             :  * discard GCPs that transform badly.
     364             :  */
     365           1 : void *GDALCreateGCPRefineTransformer(int nGCPCount, const GDAL_GCP *pasGCPList,
     366             :                                      int nReqOrder, int bReversed,
     367             :                                      double dfTolerance, int nMinimumGcps)
     368             : 
     369             : {
     370           1 :     return GDALCreateGCPTransformerEx(nGCPCount, pasGCPList, nReqOrder,
     371           1 :                                       CPL_TO_BOOL(bReversed), true, dfTolerance,
     372           1 :                                       nMinimumGcps);
     373             : }
     374             : 
     375             : /************************************************************************/
     376             : /*                     GDALDestroyGCPTransformer()                      */
     377             : /************************************************************************/
     378             : 
     379             : /**
     380             :  * Destroy GCP transformer.
     381             :  *
     382             :  * This function is used to destroy information about a GCP based
     383             :  * polynomial transformation created with GDALCreateGCPTransformer().
     384             :  *
     385             :  * @param pTransformArg the transform arg previously returned by
     386             :  * GDALCreateGCPTransformer().
     387             :  */
     388             : 
     389          51 : void GDALDestroyGCPTransformer(void *pTransformArg)
     390             : 
     391             : {
     392          51 :     if (pTransformArg == nullptr)
     393           0 :         return;
     394             : 
     395          51 :     GCPTransformInfo *psInfo = static_cast<GCPTransformInfo *>(pTransformArg);
     396             : 
     397          51 :     if (CPLAtomicDec(&(psInfo->nRefCount)) == 0)
     398             :     {
     399          51 :         delete psInfo;
     400             :     }
     401             : }
     402             : 
     403             : /************************************************************************/
     404             : /*                          GDALGCPTransform()                          */
     405             : /************************************************************************/
     406             : 
     407             : /**
     408             :  * Transforms point based on GCP derived polynomial model.
     409             :  *
     410             :  * This function matches the GDALTransformerFunc signature, and can be
     411             :  * used to transform one or more points from pixel/line coordinates to
     412             :  * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
     413             :  *
     414             :  * @param pTransformArg return value from GDALCreateGCPTransformer().
     415             :  * @param bDstToSrc TRUE if transformation is from the destination
     416             :  * (georeferenced) coordinates to pixel/line or FALSE when transforming
     417             :  * from pixel/line to georeferenced coordinates.
     418             :  * @param nPointCount the number of values in the x, y and z arrays.
     419             :  * @param x array containing the X values to be transformed.
     420             :  * @param y array containing the Y values to be transformed.
     421             :  * @param z array containing the Z values to be transformed.
     422             :  * @param panSuccess array in which a flag indicating success (TRUE) or
     423             :  * failure (FALSE) of the transformation are placed.
     424             :  *
     425             :  * @return TRUE if all points have been successfully transformed.
     426             :  */
     427             : 
     428       17812 : int GDALGCPTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
     429             :                      double *x, double *y, CPL_UNUSED double *z,
     430             :                      int *panSuccess)
     431             : 
     432             : {
     433       17812 :     int i = 0;
     434       17812 :     GCPTransformInfo *psInfo = static_cast<GCPTransformInfo *>(pTransformArg);
     435             : 
     436       17812 :     if (psInfo->bReversed)
     437           0 :         bDstToSrc = !bDstToSrc;
     438             : 
     439       17812 :     int bRet = TRUE;
     440       58807 :     for (i = 0; i < nPointCount; i++)
     441             :     {
     442       40995 :         if (x[i] == HUGE_VAL || y[i] == HUGE_VAL)
     443             :         {
     444           0 :             bRet = FALSE;
     445           0 :             panSuccess[i] = FALSE;
     446           0 :             continue;
     447             :         }
     448             : 
     449       40995 :         if (bDstToSrc)
     450             :         {
     451       28341 :             CRS_georef(x[i] - psInfo->x2_mean, y[i] - psInfo->y2_mean, x + i,
     452       28341 :                        y + i, psInfo->adfFromGeoX, psInfo->adfFromGeoY,
     453             :                        psInfo->nOrder);
     454             :         }
     455             :         else
     456             :         {
     457       12654 :             CRS_georef(x[i] - psInfo->x1_mean, y[i] - psInfo->y1_mean, x + i,
     458       12654 :                        y + i, psInfo->adfToGeoX, psInfo->adfToGeoY,
     459             :                        psInfo->nOrder);
     460             :         }
     461       40995 :         panSuccess[i] = TRUE;
     462             :     }
     463             : 
     464       17812 :     return bRet;
     465             : }
     466             : 
     467             : /************************************************************************/
     468             : /*                    GDALSerializeGCPTransformer()                     */
     469             : /************************************************************************/
     470             : 
     471           9 : CPLXMLNode *GDALSerializeGCPTransformer(void *pTransformArg)
     472             : 
     473             : {
     474           9 :     CPLXMLNode *psTree = nullptr;
     475           9 :     GCPTransformInfo *psInfo = static_cast<GCPTransformInfo *>(pTransformArg);
     476             : 
     477           9 :     VALIDATE_POINTER1(pTransformArg, "GDALSerializeGCPTransformer", nullptr);
     478             : 
     479           9 :     psTree = CPLCreateXMLNode(nullptr, CXT_Element, "GCPTransformer");
     480             : 
     481             :     /* -------------------------------------------------------------------- */
     482             :     /*      Serialize Order and bReversed.                                  */
     483             :     /* -------------------------------------------------------------------- */
     484           9 :     CPLCreateXMLElementAndValue(psTree, "Order",
     485             :                                 CPLSPrintf("%d", psInfo->nOrder));
     486             : 
     487           9 :     CPLCreateXMLElementAndValue(psTree, "Reversed",
     488             :                                 CPLSPrintf("%d", psInfo->bReversed));
     489             : 
     490           9 :     if (psInfo->bRefine)
     491             :     {
     492           0 :         CPLCreateXMLElementAndValue(psTree, "Refine",
     493             :                                     CPLSPrintf("%d", psInfo->bRefine));
     494             : 
     495           0 :         CPLCreateXMLElementAndValue(psTree, "MinimumGcps",
     496             :                                     CPLSPrintf("%d", psInfo->nMinimumGcps));
     497             : 
     498           0 :         CPLCreateXMLElementAndValue(psTree, "Tolerance",
     499             :                                     CPLSPrintf("%f", psInfo->dfTolerance));
     500             :     }
     501             : 
     502             :     /* -------------------------------------------------------------------- */
     503             :     /*     Attach GCP List.                                                 */
     504             :     /* -------------------------------------------------------------------- */
     505           9 :     if (!psInfo->asGCPs.empty())
     506             :     {
     507           9 :         if (psInfo->bRefine)
     508             :         {
     509           0 :             remove_outliers(psInfo);
     510             :         }
     511             : 
     512           9 :         GDALSerializeGCPListToXML(psTree, psInfo->asGCPs, nullptr);
     513             :     }
     514             : 
     515           9 :     return psTree;
     516             : }
     517             : 
     518             : /************************************************************************/
     519             : /*               GDALDeserializeReprojectionTransformer()               */
     520             : /************************************************************************/
     521             : 
     522          12 : void *GDALDeserializeGCPTransformer(CPLXMLNode *psTree)
     523             : 
     524             : {
     525          12 :     std::vector<gdal::GCP> asGCPs;
     526          12 :     void *pResult = nullptr;
     527          12 :     int nReqOrder = 0;
     528          12 :     int bReversed = 0;
     529          12 :     int bRefine = 0;
     530          12 :     int nMinimumGcps = 0;
     531          12 :     double dfTolerance = 0.0;
     532             : 
     533             :     /* -------------------------------------------------------------------- */
     534             :     /*      Check for GCPs.                                                 */
     535             :     /* -------------------------------------------------------------------- */
     536          12 :     CPLXMLNode *psGCPList = CPLGetXMLNode(psTree, "GCPList");
     537             : 
     538          12 :     if (psGCPList != nullptr)
     539             :     {
     540          12 :         GDALDeserializeGCPListFromXML(psGCPList, asGCPs, nullptr);
     541             :     }
     542             : 
     543             :     /* -------------------------------------------------------------------- */
     544             :     /*      Get other flags.                                                */
     545             :     /* -------------------------------------------------------------------- */
     546          12 :     nReqOrder = atoi(CPLGetXMLValue(psTree, "Order", "3"));
     547          12 :     bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
     548          12 :     bRefine = atoi(CPLGetXMLValue(psTree, "Refine", "0"));
     549          12 :     nMinimumGcps = atoi(CPLGetXMLValue(psTree, "MinimumGcps", "6"));
     550          12 :     dfTolerance = CPLAtof(CPLGetXMLValue(psTree, "Tolerance", "1.0"));
     551             : 
     552             :     /* -------------------------------------------------------------------- */
     553             :     /*      Generate transformation.                                        */
     554             :     /* -------------------------------------------------------------------- */
     555          12 :     if (bRefine)
     556             :     {
     557           2 :         pResult = GDALCreateGCPRefineTransformer(
     558           1 :             static_cast<int>(asGCPs.size()), gdal::GCP::c_ptr(asGCPs),
     559             :             nReqOrder, bReversed, dfTolerance, nMinimumGcps);
     560             :     }
     561             :     else
     562             :     {
     563          11 :         pResult = GDALCreateGCPTransformer(static_cast<int>(asGCPs.size()),
     564             :                                            gdal::GCP::c_ptr(asGCPs), nReqOrder,
     565             :                                            bReversed);
     566             :     }
     567             : 
     568          24 :     return pResult;
     569             : }
     570             : 
     571             : /************************************************************************/
     572             : /* ==================================================================== */
     573             : /*      Everything below this point derived from the CRS.C from GRASS.  */
     574             : /* ==================================================================== */
     575             : /************************************************************************/
     576             : 
     577             : /* STRUCTURE FOR USE INTERNALLY WITH THESE FUNCTIONS.  THESE FUNCTIONS EXPECT
     578             :    SQUARE MATRICES SO ONLY ONE VARIABLE IS GIVEN (N) FOR THE MATRIX SIZE */
     579             : 
     580             : struct MATRIX
     581             : {
     582             :     int n; /* SIZE OF THIS MATRIX (N x N) */
     583             :     double *v;
     584             : };
     585             : 
     586             : /* CALCULATE OFFSET INTO ARRAY BASED ON R/C */
     587             : 
     588             : #define M(row, col) m->v[(((row)-1) * (m->n)) + (col)-1]
     589             : 
     590             : /***************************************************************************/
     591             : /*
     592             :     FUNCTION PROTOTYPES FOR STATIC (INTERNAL) FUNCTIONS
     593             : */
     594             : /***************************************************************************/
     595             : 
     596             : static int calccoef(struct Control_Points *, double, double, double *, double *,
     597             :                     int);
     598             : static int calcls(struct Control_Points *, struct MATRIX *, double, double,
     599             :                   double *, double *, double *, double *);
     600             : static int exactdet(struct Control_Points *, struct MATRIX *, double, double,
     601             :                     double *, double *, double *, double *);
     602             : static int solvemat(struct MATRIX *, double *, double *, double *, double *);
     603             : static double term(int, double, double);
     604             : 
     605             : /***************************************************************************/
     606             : /*
     607             :     TRANSFORM A SINGLE COORDINATE PAIR.
     608             : */
     609             : /***************************************************************************/
     610             : 
     611             : static int
     612       41022 : CRS_georef(double e1,  /* EASTINGS TO BE TRANSFORMED */
     613             :            double n1,  /* NORTHINGS TO BE TRANSFORMED */
     614             :            double *e,  /* EASTINGS TO BE TRANSFORMED */
     615             :            double *n,  /* NORTHINGS TO BE TRANSFORMED */
     616             :            double E[], /* EASTING COEFFICIENTS */
     617             :            double N[], /* NORTHING COEFFICIENTS */
     618             :            int order /* ORDER OF TRANSFORMATION TO BE PERFORMED, MUST MATCH THE
     619             :                      ORDER USED TO CALCULATE THE COEFFICIENTS */
     620             : )
     621             : {
     622       41022 :     double e3 = 0.0;
     623       41022 :     double e2n = 0.0;
     624       41022 :     double en2 = 0.0;
     625       41022 :     double n3 = 0.0;
     626       41022 :     double e2 = 0.0;
     627       41022 :     double en = 0.0;
     628       41022 :     double n2 = 0.0;
     629             : 
     630       41022 :     switch (order)
     631             :     {
     632       38869 :         case 1:
     633             : 
     634       38869 :             *e = E[0] + E[1] * e1 + E[2] * n1;
     635       38869 :             *n = N[0] + N[1] * e1 + N[2] * n1;
     636       38869 :             break;
     637             : 
     638        2153 :         case 2:
     639             : 
     640        2153 :             e2 = e1 * e1;
     641        2153 :             n2 = n1 * n1;
     642        2153 :             en = e1 * n1;
     643             : 
     644        2153 :             *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en +
     645        2153 :                  E[5] * n2;
     646        2153 :             *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en +
     647        2153 :                  N[5] * n2;
     648        2153 :             break;
     649             : 
     650           0 :         case 3:
     651             : 
     652           0 :             e2 = e1 * e1;
     653           0 :             en = e1 * n1;
     654           0 :             n2 = n1 * n1;
     655           0 :             e3 = e1 * e2;
     656           0 :             e2n = e2 * n1;
     657           0 :             en2 = e1 * n2;
     658           0 :             n3 = n1 * n2;
     659             : 
     660           0 :             *e = E[0] + E[1] * e1 + E[2] * n1 + E[3] * e2 + E[4] * en +
     661           0 :                  E[5] * n2 + E[6] * e3 + E[7] * e2n + E[8] * en2 + E[9] * n3;
     662           0 :             *n = N[0] + N[1] * e1 + N[2] * n1 + N[3] * e2 + N[4] * en +
     663           0 :                  N[5] * n2 + N[6] * e3 + N[7] * e2n + N[8] * en2 + N[9] * n3;
     664           0 :             break;
     665             : 
     666           0 :         default:
     667             : 
     668           0 :             return (MPARMERR);
     669             :     }
     670             : 
     671       41022 :     return (MSUCCESS);
     672             : }
     673             : 
     674             : /***************************************************************************/
     675             : /*
     676             :     COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
     677             : */
     678             : /***************************************************************************/
     679             : 
     680          53 : static int CRS_compute_georef_equations(GCPTransformInfo *psInfo,
     681             :                                         struct Control_Points *cp, double E12[],
     682             :                                         double N12[], double E21[],
     683             :                                         double N21[], int order)
     684             : {
     685          53 :     double *tempptr = nullptr;
     686          53 :     int status = 0;
     687             : 
     688          53 :     if (order < 1 || order > MAXORDER)
     689           0 :         return (MPARMERR);
     690             : 
     691             :     /* CALCULATE THE FORWARD TRANSFORMATION COEFFICIENTS */
     692             : 
     693          53 :     status = calccoef(cp, psInfo->x1_mean, psInfo->y1_mean, E12, N12, order);
     694          53 :     if (status != MSUCCESS)
     695           1 :         return (status);
     696             : 
     697             :     /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS */
     698             : 
     699          52 :     tempptr = cp->e1;
     700          52 :     cp->e1 = cp->e2;
     701          52 :     cp->e2 = tempptr;
     702          52 :     tempptr = cp->n1;
     703          52 :     cp->n1 = cp->n2;
     704          52 :     cp->n2 = tempptr;
     705             : 
     706             :     /* CALCULATE THE BACKWARD TRANSFORMATION COEFFICIENTS */
     707             : 
     708          52 :     status = calccoef(cp, psInfo->x2_mean, psInfo->y2_mean, E21, N21, order);
     709             : 
     710             :     /* SWITCH THE 1 AND 2 EASTING AND NORTHING ARRAYS BACK */
     711             : 
     712          52 :     tempptr = cp->e1;
     713          52 :     cp->e1 = cp->e2;
     714          52 :     cp->e2 = tempptr;
     715          52 :     tempptr = cp->n1;
     716          52 :     cp->n1 = cp->n2;
     717          52 :     cp->n2 = tempptr;
     718             : 
     719          52 :     return (status);
     720             : }
     721             : 
     722             : /***************************************************************************/
     723             : /*
     724             :     COMPUTE THE GEOREFFERENCING COEFFICIENTS BASED ON A SET OF CONTROL POINTS
     725             : */
     726             : /***************************************************************************/
     727             : 
     728         105 : static int calccoef(struct Control_Points *cp, double x_mean, double y_mean,
     729             :                     double E[], double N[], int order)
     730             : {
     731             :     struct MATRIX m;
     732         105 :     double *a = nullptr;
     733         105 :     double *b = nullptr;
     734         105 :     int numactive = 0; /* NUMBER OF ACTIVE CONTROL POINTS */
     735         105 :     int status = 0;
     736         105 :     int i = 0;
     737             : 
     738         105 :     memset(&m, 0, sizeof(m));
     739             : 
     740             :     /* CALCULATE THE NUMBER OF VALID CONTROL POINTS */
     741             : 
     742        1772 :     for (i = numactive = 0; i < cp->count; i++)
     743             :     {
     744        1667 :         if (cp->status[i] > 0)
     745        1667 :             numactive++;
     746             :     }
     747             : 
     748             :     /* CALCULATE THE MINIMUM NUMBER OF CONTROL POINTS NEEDED TO DETERMINE
     749             :        A TRANSFORMATION OF THIS ORDER */
     750             : 
     751         105 :     m.n = ((order + 1) * (order + 2)) / 2;
     752             : 
     753         105 :     if (numactive < m.n)
     754           1 :         return (MNPTERR);
     755             : 
     756             :     /* INITIALIZE MATRIX */
     757             : 
     758         104 :     m.v = static_cast<double *>(
     759         104 :         VSICalloc(cpl::fits_on<int>(m.n * m.n), sizeof(double)));
     760         104 :     if (m.v == nullptr)
     761             :     {
     762           0 :         return (MMEMERR);
     763             :     }
     764         104 :     a = static_cast<double *>(VSICalloc(m.n, sizeof(double)));
     765         104 :     if (a == nullptr)
     766             :     {
     767           0 :         CPLFree(m.v);
     768           0 :         return (MMEMERR);
     769             :     }
     770         104 :     b = static_cast<double *>(VSICalloc(m.n, sizeof(double)));
     771         104 :     if (b == nullptr)
     772             :     {
     773           0 :         CPLFree(m.v);
     774           0 :         CPLFree(a);
     775           0 :         return (MMEMERR);
     776             :     }
     777             : 
     778         104 :     if (numactive == m.n)
     779          24 :         status = exactdet(cp, &m, x_mean, y_mean, a, b, E, N);
     780             :     else
     781          80 :         status = calcls(cp, &m, x_mean, y_mean, a, b, E, N);
     782             : 
     783         104 :     CPLFree(m.v);
     784         104 :     CPLFree(a);
     785         104 :     CPLFree(b);
     786             : 
     787         104 :     return (status);
     788             : }
     789             : 
     790             : /***************************************************************************/
     791             : /*
     792             :     CALCULATE THE TRANSFORMATION COEFFICIENTS WITH EXACTLY THE MINIMUM
     793             :     NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.
     794             : */
     795             : /***************************************************************************/
     796             : 
     797          24 : static int exactdet(struct Control_Points *cp, struct MATRIX *m, double x_mean,
     798             :                     double y_mean, double a[], double b[],
     799             :                     double E[], /* EASTING COEFFICIENTS */
     800             :                     double N[]  /* NORTHING COEFFICIENTS */
     801             : )
     802             : {
     803          24 :     int currow = 1;
     804             : 
     805          96 :     for (int pntnow = 0; pntnow < cp->count; pntnow++)
     806             :     {
     807          72 :         if (cp->status[pntnow] > 0)
     808             :         {
     809             :             /* POPULATE MATRIX M */
     810             : 
     811         288 :             for (int j = 1; j <= m->n; j++)
     812             :             {
     813         216 :                 M(currow, j) =
     814         216 :                     term(j, cp->e1[pntnow] - x_mean, cp->n1[pntnow] - y_mean);
     815             :             }
     816             : 
     817             :             /* POPULATE MATRIX A AND B */
     818             : 
     819          72 :             a[currow - 1] = cp->e2[pntnow];
     820          72 :             b[currow - 1] = cp->n2[pntnow];
     821             : 
     822          72 :             currow++;
     823             :         }
     824             :     }
     825             : 
     826          24 :     if (currow - 1 != m->n)
     827           0 :         return (MINTERR);
     828             : 
     829          24 :     return (solvemat(m, a, b, E, N));
     830             : }
     831             : 
     832             : /***************************************************************************/
     833             : /*
     834             :     CALCULATE THE TRANSFORMATION COEFFICIENTS WITH MORE THAN THE MINIMUM
     835             :     NUMBER OF CONTROL POINTS REQUIRED FOR THIS TRANSFORMATION.  THIS
     836             :     ROUTINE USES THE LEAST SQUARES METHOD TO COMPUTE THE COEFFICIENTS.
     837             : */
     838             : /***************************************************************************/
     839             : 
     840          80 : static int calcls(struct Control_Points *cp, struct MATRIX *m, double x_mean,
     841             :                   double y_mean, double a[], double b[],
     842             :                   double E[], /* EASTING COEFFICIENTS */
     843             :                   double N[]  /* NORTHING COEFFICIENTS */
     844             : )
     845             : {
     846          80 :     int numactive = 0;
     847             : 
     848             :     /* INITIALIZE THE UPPER HALF OF THE MATRIX AND THE TWO COLUMN VECTORS */
     849             : 
     850         344 :     for (int i = 1; i <= m->n; i++)
     851             :     {
     852         864 :         for (int j = i; j <= m->n; j++)
     853         600 :             M(i, j) = 0.0;
     854         264 :         a[i - 1] = b[i - 1] = 0.0;
     855             :     }
     856             : 
     857             :     /* SUM THE UPPER HALF OF THE MATRIX AND THE COLUMN VECTORS ACCORDING TO
     858             :        THE LEAST SQUARES METHOD OF SOLVING OVER DETERMINED SYSTEMS */
     859             : 
     860        1674 :     for (int n = 0; n < cp->count; n++)
     861             :     {
     862        1594 :         if (cp->status[n] > 0)
     863             :         {
     864        1594 :             numactive++;
     865       10204 :             for (int i = 1; i <= m->n; i++)
     866             :             {
     867       37314 :                 for (int j = i; j <= m->n; j++)
     868       57408 :                     M(i, j) += term(i, cp->e1[n] - x_mean, cp->n1[n] - y_mean) *
     869       28704 :                                term(j, cp->e1[n] - x_mean, cp->n1[n] - y_mean);
     870             : 
     871        8610 :                 a[i - 1] +=
     872        8610 :                     cp->e2[n] * term(i, cp->e1[n] - x_mean, cp->n1[n] - y_mean);
     873        8610 :                 b[i - 1] +=
     874        8610 :                     cp->n2[n] * term(i, cp->e1[n] - x_mean, cp->n1[n] - y_mean);
     875             :             }
     876             :         }
     877             :     }
     878             : 
     879          80 :     if (numactive <= m->n)
     880           0 :         return (MINTERR);
     881             : 
     882             :     /* TRANSPOSE VALUES IN UPPER HALF OF M TO OTHER HALF */
     883             : 
     884         264 :     for (int i = 2; i <= m->n; i++)
     885             :     {
     886         520 :         for (int j = 1; j < i; j++)
     887         336 :             M(i, j) = M(j, i);
     888             :     }
     889             : 
     890          80 :     return (solvemat(m, a, b, E, N));
     891             : }
     892             : 
     893             : /***************************************************************************/
     894             : /*
     895             :     CALCULATE THE X/Y TERM BASED ON THE TERM NUMBER
     896             : 
     897             : ORDER\TERM   1    2    3    4    5    6    7    8    9   10
     898             :   1        e0n0 e1n0 e0n1
     899             :   2        e0n0 e1n0 e0n1 e2n0 e1n1 e0n2
     900             :   3        e0n0 e1n0 e0n1 e2n0 e1n1 e0n2 e3n0 e2n1 e1n2 e0n3
     901             : */
     902             : /***************************************************************************/
     903             : 
     904       74844 : static double term(int nTerm, double e, double n)
     905             : {
     906       74844 :     switch (nTerm)
     907             :     {
     908       13464 :         case 1:
     909       13464 :             return (1.0);
     910       13464 :         case 2:
     911       13464 :             return (e);
     912       13464 :         case 3:
     913       13464 :             return (n);
     914       11484 :         case 4:
     915       11484 :             return ((e * e));
     916       11484 :         case 5:
     917       11484 :             return ((e * n));
     918       11484 :         case 6:
     919       11484 :             return ((n * n));
     920           0 :         case 7:
     921           0 :             return ((e * e * e));
     922           0 :         case 8:
     923           0 :             return ((e * e * n));
     924           0 :         case 9:
     925           0 :             return ((e * n * n));
     926           0 :         case 10:
     927           0 :             return ((n * n * n));
     928             :     }
     929           0 :     return 0.0;
     930             : }
     931             : 
     932             : /***************************************************************************/
     933             : /*
     934             :     SOLVE FOR THE 'E' AND 'N' COEFFICIENTS BY USING A SOMEWHAT MODIFIED
     935             :     GAUSSIAN ELIMINATION METHOD.
     936             : 
     937             :     | M11 M12 ... M1n | | E0   |   | a0   |
     938             :     | M21 M22 ... M2n | | E1   | = | a1   |
     939             :     |  .   .   .   .  | | .    |   | .    |
     940             :     | Mn1 Mn2 ... Mnn | | En-1 |   | an-1 |
     941             : 
     942             :     and
     943             : 
     944             :     | M11 M12 ... M1n | | N0   |   | b0   |
     945             :     | M21 M22 ... M2n | | N1   | = | b1   |
     946             :     |  .   .   .   .  | | .    |   | .    |
     947             :     | Mn1 Mn2 ... Mnn | | Nn-1 |   | bn-1 |
     948             : */
     949             : /***************************************************************************/
     950             : 
     951         104 : static int solvemat(struct MATRIX *m, double a[], double b[], double E[],
     952             :                     double N[])
     953             : {
     954         440 :     for (int i = 1; i <= m->n; i++)
     955             :     {
     956         336 :         int j = i;
     957             : 
     958             :         /* find row with largest magnitude value for pivot value */
     959             : 
     960         336 :         double pivot =
     961         336 :             M(i, j); /* ACTUAL VALUE OF THE LARGEST PIVOT CANDIDATE */
     962         336 :         int imark = i;
     963         744 :         for (int i2 = i + 1; i2 <= m->n; i2++)
     964             :         {
     965         408 :             if (fabs(M(i2, j)) > fabs(pivot))
     966             :             {
     967          53 :                 pivot = M(i2, j);
     968          53 :                 imark = i2;
     969             :             }
     970             :         }
     971             : 
     972             :         /* if the pivot is very small then the points are nearly co-linear */
     973             :         /* co-linear points result in an undefined matrix, and nearly */
     974             :         /* co-linear points results in a solution with rounding error */
     975             : 
     976         336 :         if (pivot == 0.0)
     977           0 :             return (MUNSOLVABLE);
     978             : 
     979             :         /* if row with highest pivot is not the current row, switch them */
     980             : 
     981         336 :         if (imark != i)
     982             :         {
     983         245 :             for (int j2 = 1; j2 <= m->n; j2++)
     984             :             {
     985         198 :                 std::swap(M(imark, j2), M(i, j2));
     986             :             }
     987             : 
     988          47 :             std::swap(a[imark - 1], a[i - 1]);
     989          47 :             std::swap(b[imark - 1], b[i - 1]);
     990             :         }
     991             : 
     992             :         /* compute zeros above and below the pivot, and compute
     993             :            values for the rest of the row as well */
     994             : 
     995        1488 :         for (int i2 = 1; i2 <= m->n; i2++)
     996             :         {
     997        1152 :             if (i2 != i)
     998             :             {
     999         816 :                 const double factor = M(i2, j) / pivot;
    1000        2808 :                 for (int j2 = j; j2 <= m->n; j2++)
    1001        1992 :                     M(i2, j2) -= factor * M(i, j2);
    1002         816 :                 a[i2 - 1] -= factor * a[i - 1];
    1003         816 :                 b[i2 - 1] -= factor * b[i - 1];
    1004             :             }
    1005             :         }
    1006             :     }
    1007             : 
    1008             :     /* SINCE ALL OTHER VALUES IN THE MATRIX ARE ZERO NOW, CALCULATE THE
    1009             :        COEFFICIENTS BY DIVIDING THE COLUMN VECTORS BY THE DIAGONAL VALUES. */
    1010             : 
    1011         440 :     for (int i = 1; i <= m->n; i++)
    1012             :     {
    1013         336 :         E[i - 1] = a[i - 1] / M(i, i);
    1014         336 :         N[i - 1] = b[i - 1] / M(i, i);
    1015             :     }
    1016             : 
    1017         104 :     return (MSUCCESS);
    1018             : }
    1019             : 
    1020             : /***************************************************************************/
    1021             : /*
    1022             :   DETECTS THE WORST OUTLIER IN THE GCP LIST AND RETURNS THE INDEX OF THE
    1023             :   OUTLIER.
    1024             : 
    1025             :   THE WORST OUTLIER IS CALCULATED BASED ON THE CONTROL POINTS, COEFFICIENTS
    1026             :   AND THE ALLOWED TOLERANCE:
    1027             : 
    1028             :   sampleAdj = a0 + a1*sample + a2*line + a3*line*sample
    1029             :   lineAdj = b0 + b1*sample + b2*line + b3*line*sample
    1030             : 
    1031             :   WHERE sampleAdj AND lineAdj ARE CORRELATED GCPS
    1032             : 
    1033             :   [residualSample] = [A1][sampleCoefficients] - [b1]
    1034             :   [residualLine] = [A2][lineCoefficients] - [b2]
    1035             : 
    1036             :   sampleResidual^2 = sum( [residualSample]^2 )
    1037             :   lineResidual^2 = sum( [lineSample]^2 )
    1038             : 
    1039             :   residuals(i) = squareRoot( residualSample(i)^2 + residualLine(i)^2 )
    1040             : 
    1041             :   THE GCP WITH THE GREATEST DISTANCE residual(i) GREATER THAN THE TOLERANCE WILL
    1042             :   CONSIDERED THE WORST OUTLIER.
    1043             : 
    1044             :   IF NO OUTLIER CAN BE FOUND, -1 WILL BE RETURNED.
    1045             : */
    1046             : /***************************************************************************/
    1047           3 : static int worst_outlier(struct Control_Points *cp, double x_mean,
    1048             :                          double y_mean, int nOrder, double E[], double N[],
    1049             :                          double dfTolerance)
    1050             : {
    1051             :     // double dfSampleResidual = 0.0;
    1052             :     // double dfLineResidual = 0.0;
    1053             :     double *padfResiduals =
    1054           3 :         static_cast<double *>(CPLCalloc(sizeof(double), cp->count));
    1055             : 
    1056          30 :     for (int nI = 0; nI < cp->count; nI++)
    1057             :     {
    1058          27 :         double dfSampleRes = 0.0;
    1059          27 :         double dfLineRes = 0.0;
    1060          27 :         CRS_georef(cp->e1[nI] - x_mean, cp->n1[nI] - y_mean, &dfSampleRes,
    1061             :                    &dfLineRes, E, N, nOrder);
    1062          27 :         dfSampleRes -= cp->e2[nI];
    1063          27 :         dfLineRes -= cp->n2[nI];
    1064             :         // dfSampleResidual += dfSampleRes*dfSampleRes;
    1065             :         // dfLineResidual += dfLineRes*dfLineRes;
    1066             : 
    1067          27 :         padfResiduals[nI] =
    1068          27 :             sqrt(dfSampleRes * dfSampleRes + dfLineRes * dfLineRes);
    1069             :     }
    1070             : 
    1071           3 :     int nIndex = -1;
    1072           3 :     double dfDifference = -1.0;
    1073          30 :     for (int nI = 0; nI < cp->count; nI++)
    1074             :     {
    1075          27 :         double dfCurrentDifference = padfResiduals[nI];
    1076          27 :         if (fabs(dfCurrentDifference) < 1.19209290E-07F /*FLT_EPSILON*/)
    1077             :         {
    1078           8 :             dfCurrentDifference = 0.0;
    1079             :         }
    1080          27 :         if (dfCurrentDifference > dfDifference &&
    1081             :             dfCurrentDifference >= dfTolerance)
    1082             :         {
    1083           5 :             dfDifference = dfCurrentDifference;
    1084           5 :             nIndex = nI;
    1085             :         }
    1086             :     }
    1087           3 :     CPLFree(padfResiduals);
    1088           3 :     return nIndex;
    1089             : }
    1090             : 
    1091             : /***************************************************************************/
    1092             : /*
    1093             :   REMOVES THE WORST OUTLIERS ITERATIVELY UNTIL THE MINIMUM NUMBER OF GCPS
    1094             :   ARE REACHED OR NO OUTLIERS CAN BE DETECTED.
    1095             : 
    1096             :   1. WE CALCULATE THE COEFFICIENTS FOR ALL THE GCPS.
    1097             :   2. THE GCP LIST WILL BE SCANNED TO DETERMINE THE WORST OUTLIER USING
    1098             :      THE CALCULATED COEFFICIENTS.
    1099             :   3. THE WORST OUTLIER WILL BE REMOVED FROM THE GCP LIST.
    1100             :   4. THE COEFFICIENTS WILL BE RECALCULATED WITHOUT THE WORST OUTLIER.
    1101             :   5. STEP 1 TO 4 ARE EXECUTED UNTIL THE MINIMUM NUMBER OF GCPS WERE REACHED
    1102             :      OR IF NO GCP IS CONSIDERED AN OUTLIER WITH THE PASSED TOLERANCE.
    1103             : */
    1104             : /***************************************************************************/
    1105           1 : static int remove_outliers(GCPTransformInfo *psInfo)
    1106             : {
    1107           1 :     double *padfGeoX = nullptr;
    1108           1 :     double *padfGeoY = nullptr;
    1109           1 :     double *padfRasterX = nullptr;
    1110           1 :     double *padfRasterY = nullptr;
    1111           1 :     int *panStatus = nullptr;
    1112           1 :     int nCRSresult = 0;
    1113           1 :     int nGCPCount = 0;
    1114           1 :     int nMinimumGcps = 0;
    1115           1 :     int nReqOrder = 0;
    1116           1 :     double dfTolerance = 0;
    1117             :     struct Control_Points sPoints;
    1118             : 
    1119           1 :     double x1_sum = 0;
    1120           1 :     double y1_sum = 0;
    1121           1 :     double x2_sum = 0;
    1122           1 :     double y2_sum = 0;
    1123           1 :     memset(&sPoints, 0, sizeof(sPoints));
    1124             : 
    1125           1 :     nGCPCount = static_cast<int>(psInfo->asGCPs.size());
    1126           1 :     nMinimumGcps = psInfo->nMinimumGcps;
    1127           1 :     nReqOrder = psInfo->nOrder;
    1128           1 :     dfTolerance = psInfo->dfTolerance;
    1129             : 
    1130             :     try
    1131             :     {
    1132           1 :         padfGeoX = new double[nGCPCount];
    1133           1 :         padfGeoY = new double[nGCPCount];
    1134           1 :         padfRasterX = new double[nGCPCount];
    1135           1 :         padfRasterY = new double[nGCPCount];
    1136           1 :         panStatus = new int[nGCPCount];
    1137             : 
    1138          11 :         for (int nI = 0; nI < nGCPCount; nI++)
    1139             :         {
    1140          10 :             panStatus[nI] = 1;
    1141          10 :             padfGeoX[nI] = psInfo->asGCPs[nI].X();
    1142          10 :             padfGeoY[nI] = psInfo->asGCPs[nI].Y();
    1143          10 :             padfRasterX[nI] = psInfo->asGCPs[nI].Pixel();
    1144          10 :             padfRasterY[nI] = psInfo->asGCPs[nI].Line();
    1145          10 :             x1_sum += psInfo->asGCPs[nI].Pixel();
    1146          10 :             y1_sum += psInfo->asGCPs[nI].Line();
    1147          10 :             x2_sum += psInfo->asGCPs[nI].X();
    1148          10 :             y2_sum += psInfo->asGCPs[nI].Y();
    1149             :         }
    1150           1 :         psInfo->x1_mean = x1_sum / nGCPCount;
    1151           1 :         psInfo->y1_mean = y1_sum / nGCPCount;
    1152           1 :         psInfo->x2_mean = x2_sum / nGCPCount;
    1153           1 :         psInfo->y2_mean = y2_sum / nGCPCount;
    1154             : 
    1155           1 :         sPoints.count = nGCPCount;
    1156           1 :         sPoints.e1 = padfRasterX;
    1157           1 :         sPoints.n1 = padfRasterY;
    1158           1 :         sPoints.e2 = padfGeoX;
    1159           1 :         sPoints.n2 = padfGeoY;
    1160           1 :         sPoints.status = panStatus;
    1161             : 
    1162           1 :         nCRSresult = CRS_compute_georef_equations(
    1163           1 :             psInfo, &sPoints, psInfo->adfToGeoX, psInfo->adfToGeoY,
    1164           1 :             psInfo->adfFromGeoX, psInfo->adfFromGeoY, nReqOrder);
    1165             : 
    1166           3 :         while (sPoints.count > nMinimumGcps)
    1167             :         {
    1168           6 :             int nIndex = worst_outlier(
    1169             :                 &sPoints, psInfo->x1_mean, psInfo->y1_mean, psInfo->nOrder,
    1170           3 :                 psInfo->adfToGeoX, psInfo->adfToGeoY, dfTolerance);
    1171             : 
    1172             :             // If no outliers were detected, stop the GCP elimination
    1173           3 :             if (nIndex == -1)
    1174             :             {
    1175           1 :                 break;
    1176             :             }
    1177             : 
    1178           8 :             for (int nI = nIndex; nI < sPoints.count - 1; nI++)
    1179             :             {
    1180           6 :                 sPoints.e1[nI] = sPoints.e1[nI + 1];
    1181           6 :                 sPoints.n1[nI] = sPoints.n1[nI + 1];
    1182           6 :                 sPoints.e2[nI] = sPoints.e2[nI + 1];
    1183           6 :                 sPoints.n2[nI] = sPoints.n2[nI + 1];
    1184           6 :                 psInfo->asGCPs[nI].SetId(psInfo->asGCPs[nI + 1].Id());
    1185           6 :                 psInfo->asGCPs[nI].SetInfo(psInfo->asGCPs[nI + 1].Info());
    1186             :             }
    1187             : 
    1188           2 :             sPoints.count = sPoints.count - 1;
    1189             : 
    1190           2 :             nCRSresult = CRS_compute_georef_equations(
    1191           2 :                 psInfo, &sPoints, psInfo->adfToGeoX, psInfo->adfToGeoY,
    1192           2 :                 psInfo->adfFromGeoX, psInfo->adfFromGeoY, nReqOrder);
    1193             :         }
    1194             : 
    1195           9 :         for (int nI = 0; nI < sPoints.count; nI++)
    1196             :         {
    1197           8 :             psInfo->asGCPs[nI].X() = sPoints.e2[nI];
    1198           8 :             psInfo->asGCPs[nI].Y() = sPoints.n2[nI];
    1199           8 :             psInfo->asGCPs[nI].Pixel() = sPoints.e1[nI];
    1200           8 :             psInfo->asGCPs[nI].Line() = sPoints.n1[nI];
    1201             :         }
    1202           1 :         psInfo->asGCPs.resize(sPoints.count);
    1203             :     }
    1204           0 :     catch (const std::exception &e)
    1205             :     {
    1206           0 :         CPLError(CE_Failure, CPLE_OutOfMemory, "%s", e.what());
    1207           0 :         nCRSresult = MINTERR;
    1208             :     }
    1209           1 :     delete[] padfGeoX;
    1210           1 :     delete[] padfGeoY;
    1211           1 :     delete[] padfRasterX;
    1212           1 :     delete[] padfRasterY;
    1213           1 :     delete[] panStatus;
    1214             : 
    1215           1 :     return nCRSresult;
    1216             : }

Generated by: LCOV version 1.14