Line data    Source code

       1             : /******************************************************************************
       2             :  *
       3             :  * Project:  High Performance Image Reprojector
       4             :  * Purpose:  Thin Plate Spline transformer (GDAL wrapper portion)
       5             :  * Author:   Frank Warmerdam, warmerdam@pobox.com
       6             :  *
       7             :  ******************************************************************************
       8             :  * Copyright (c) 2004, Frank Warmerdam <warmerdam@pobox.com>
       9             :  * Copyright (c) 2011-2013, Even Rouault <even dot rouault at spatialys.com>
      10             :  *
      11             :  * SPDX-License-Identifier: MIT
      12             :  ****************************************************************************/
      13             : 
      14             : #include "cpl_port.h"
      15             : #include "thinplatespline.h"
      16             : 
      17             : #include <stdlib.h>
      18             : #include <string.h>
      19             : #include <map>
      20             : #include <utility>
      21             : 
      22             : #include "cpl_atomic_ops.h"
      23             : #include "cpl_conv.h"
      24             : #include "cpl_error.h"
      25             : #include "cpl_minixml.h"
      26             : #include "cpl_multiproc.h"
      27             : #include "cpl_string.h"
      28             : #include "gdal.h"
      29             : #include "gdal_alg.h"
      30             : #include "gdal_alg_priv.h"
      31             : #include "gdal_priv.h"
      32             : #include "gdalgenericinverse.h"
      33             : 
      34             : CPL_C_START
      35             : CPLXMLNode *GDALSerializeTPSTransformer(void *pTransformArg);
      36             : void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree);
      37             : CPL_C_END
      38             : 
      39             : struct TPSTransformInfo
      40             : {
      41             :     GDALTransformerInfo sTI{};
      42             : 
      43             :     VizGeorefSpline2D *poForward{};
      44             :     VizGeorefSpline2D *poReverse{};
      45             :     bool bForwardSolved{};
      46             :     bool bReverseSolved{};
      47             :     double dfSrcApproxErrorReverse{};
      48             : 
      49             :     bool bReversed{};
      50             : 
      51             :     std::vector<gdal::GCP> asGCPs{};
      52             : 
      53             :     volatile int nRefCount{};
      54             : };
      55             : 
      56             : /************************************************************************/
      57             : /*                   GDALCreateSimilarTPSTransformer()                  */
      58             : /************************************************************************/
      59             : 
      60           3 : static void *GDALCreateSimilarTPSTransformer(void *hTransformArg,
      61             :                                              double dfRatioX, double dfRatioY)
      62             : {
      63           3 :     VALIDATE_POINTER1(hTransformArg, "GDALCreateSimilarTPSTransformer",
      64             :                       nullptr);
      65             : 
      66           3 :     TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(hTransformArg);
      67             : 
      68           3 :     if (dfRatioX == 1.0 && dfRatioY == 1.0)
      69             :     {
      70             :         // We can just use a ref count, since using the source transformation
      71             :         // is thread-safe.
      72           0 :         CPLAtomicInc(&(psInfo->nRefCount));
      73             :     }
      74             :     else
      75             :     {
      76           3 :         auto newGCPs = psInfo->asGCPs;
      77          21 :         for (auto &gcp : newGCPs)
      78             :         {
      79          18 :             gcp.Pixel() /= dfRatioX;
      80          18 :             gcp.Line() /= dfRatioY;
      81             :         }
      82           3 :         psInfo = static_cast<TPSTransformInfo *>(GDALCreateTPSTransformer(
      83           3 :             static_cast<int>(newGCPs.size()), gdal::GCP::c_ptr(newGCPs),
      84           3 :             psInfo->bReversed));
      85             :     }
      86             : 
      87           3 :     return psInfo;
      88             : }
      89             : 
      90             : /************************************************************************/
      91             : /*                      GDALCreateTPSTransformer()                      */
      92             : /************************************************************************/
      93             : 
      94             : /**
      95             :  * Create Thin Plate Spline transformer from GCPs.
      96             :  *
      97             :  * The thin plate spline transformer produces exact transformation
      98             :  * at all control points and smoothly varying transformations between
      99             :  * control points with greatest influence from local control points.
     100             :  * It is suitable for for many applications not well modeled by polynomial
     101             :  * transformations.
     102             :  *
     103             :  * Creating the TPS transformer involves solving systems of linear equations
     104             :  * related to the number of control points involved.  This solution is
     105             :  * computed within this function call.  It can be quite an expensive operation
     106             :  * for large numbers of GCPs.  For instance, for reference, it takes on the
     107             :  * order of 10s for 400 GCPs on a 2GHz Athlon processor.
     108             :  *
     109             :  * TPS Transformers are serializable.
     110             :  *
     111             :  * The GDAL Thin Plate Spline transformer is based on code provided by
     112             :  * Gilad Ronnen on behalf of VIZRT Inc (http://www.visrt.com).  Incorporation
     113             :  * of the algorithm into GDAL was supported by the Centro di Ecologia Alpina
     114             :  * (http://www.cealp.it).
     115             :  *
     116             :  * @param nGCPCount the number of GCPs in pasGCPList.
     117             :  * @param pasGCPList an array of GCPs to be used as input.
     118             :  * @param bReversed set it to TRUE to compute the reversed transformation.
     119             :  *
     120             :  * @return the transform argument or NULL if creation fails.
     121             :  */
     122             : 
     123           5 : void *GDALCreateTPSTransformer(int nGCPCount, const GDAL_GCP *pasGCPList,
     124             :                                int bReversed)
     125             : {
     126           5 :     return GDALCreateTPSTransformerInt(nGCPCount, pasGCPList, bReversed,
     127           5 :                                        nullptr);
     128             : }
     129             : 
     130           0 : static void GDALTPSComputeForwardInThread(void *pData)
     131             : {
     132           0 :     TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pData);
     133           0 :     psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
     134           0 : }
     135             : 
     136          19 : void *GDALCreateTPSTransformerInt(int nGCPCount, const GDAL_GCP *pasGCPList,
     137             :                                   int bReversed, char **papszOptions)
     138             : 
     139             : {
     140             :     /* -------------------------------------------------------------------- */
     141             :     /*      Allocate transform info.                                        */
     142             :     /* -------------------------------------------------------------------- */
     143          19 :     TPSTransformInfo *psInfo = new TPSTransformInfo();
     144             : 
     145          19 :     psInfo->asGCPs = gdal::GCP::fromC(pasGCPList, nGCPCount);
     146             : 
     147          19 :     psInfo->bReversed = CPL_TO_BOOL(bReversed);
     148          19 :     psInfo->poForward = new VizGeorefSpline2D(2);
     149          19 :     psInfo->poReverse = new VizGeorefSpline2D(2);
     150             : 
     151          19 :     memcpy(psInfo->sTI.abySignature, GDAL_GTI2_SIGNATURE,
     152             :            strlen(GDAL_GTI2_SIGNATURE));
     153          19 :     psInfo->sTI.pszClassName = "GDALTPSTransformer";
     154          19 :     psInfo->sTI.pfnTransform = GDALTPSTransform;
     155          19 :     psInfo->sTI.pfnCleanup = GDALDestroyTPSTransformer;
     156          19 :     psInfo->sTI.pfnSerialize = GDALSerializeTPSTransformer;
     157          19 :     psInfo->sTI.pfnCreateSimilar = GDALCreateSimilarTPSTransformer;
     158             : 
     159             :     /* -------------------------------------------------------------------- */
     160             :     /*      Attach (non-redundant) points to the transformation.            */
     161             :     /* -------------------------------------------------------------------- */
     162          38 :     std::map<std::pair<double, double>, int> oMapPixelLineToIdx;
     163          38 :     std::map<std::pair<double, double>, int> oMapXYToIdx;
     164        2426 :     for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
     165             :     {
     166        2407 :         const double afPL[2] = {pasGCPList[iGCP].dfGCPPixel,
     167        2407 :                                 pasGCPList[iGCP].dfGCPLine};
     168        2407 :         const double afXY[2] = {pasGCPList[iGCP].dfGCPX,
     169        2407 :                                 pasGCPList[iGCP].dfGCPY};
     170             : 
     171             :         std::map<std::pair<double, double>, int>::iterator oIter(
     172             :             oMapPixelLineToIdx.find(
     173        2407 :                 std::pair<double, double>(afPL[0], afPL[1])));
     174        2407 :         if (oIter != oMapPixelLineToIdx.end())
     175             :         {
     176           3 :             if (afXY[0] == pasGCPList[oIter->second].dfGCPX &&
     177           1 :                 afXY[1] == pasGCPList[oIter->second].dfGCPY)
     178             :             {
     179           1 :                 continue;
     180             :             }
     181             :             else
     182             :             {
     183           1 :                 CPLError(CE_Warning, CPLE_AppDefined,
     184             :                          "GCP %d and %d have same (pixel,line)=(%f,%f), "
     185             :                          "but different (X,Y): (%f,%f) vs (%f,%f)",
     186           1 :                          iGCP + 1, oIter->second, afPL[0], afPL[1], afXY[0],
     187           1 :                          afXY[1], pasGCPList[oIter->second].dfGCPX,
     188           1 :                          pasGCPList[oIter->second].dfGCPY);
     189             :             }
     190             :         }
     191             :         else
     192             :         {
     193        2405 :             oMapPixelLineToIdx[std::pair<double, double>(afPL[0], afPL[1])] =
     194             :                 iGCP;
     195             :         }
     196             : 
     197        2406 :         oIter = oMapXYToIdx.find(std::pair<double, double>(afXY[0], afXY[1]));
     198        2406 :         if (oIter != oMapXYToIdx.end())
     199             :         {
     200           1 :             CPLError(CE_Warning, CPLE_AppDefined,
     201             :                      "GCP %d and %d have same (x,y)=(%f,%f), "
     202             :                      "but different (pixel,line): (%f,%f) vs (%f,%f)",
     203           1 :                      iGCP + 1, oIter->second, afXY[0], afXY[1], afPL[0],
     204           1 :                      afPL[1], pasGCPList[oIter->second].dfGCPPixel,
     205           1 :                      pasGCPList[oIter->second].dfGCPLine);
     206             :         }
     207             :         else
     208             :         {
     209        2405 :             oMapXYToIdx[std::pair<double, double>(afXY[0], afXY[1])] = iGCP;
     210             :         }
     211             : 
     212        2406 :         bool bOK = true;
     213        2406 :         if (bReversed)
     214             :         {
     215           0 :             bOK &= psInfo->poReverse->add_point(afPL[0], afPL[1], afXY);
     216           0 :             bOK &= psInfo->poForward->add_point(afXY[0], afXY[1], afPL);
     217             :         }
     218             :         else
     219             :         {
     220        2406 :             bOK &= psInfo->poForward->add_point(afPL[0], afPL[1], afXY);
     221        2406 :             bOK &= psInfo->poReverse->add_point(afXY[0], afXY[1], afPL);
     222             :         }
     223        2406 :         if (!bOK)
     224             :         {
     225           0 :             GDALDestroyTPSTransformer(psInfo);
     226           0 :             return nullptr;
     227             :         }
     228             :     }
     229             : 
     230          19 :     psInfo->nRefCount = 1;
     231             : 
     232          19 :     psInfo->dfSrcApproxErrorReverse = CPLAtof(
     233             :         CSLFetchNameValueDef(papszOptions, "SRC_APPROX_ERROR_IN_PIXEL", "0"));
     234             : 
     235          19 :     int nThreads = 1;
     236          19 :     if (nGCPCount > 100)
     237             :     {
     238             :         const char *pszWarpThreads =
     239           2 :             CSLFetchNameValue(papszOptions, "NUM_THREADS");
     240           2 :         if (pszWarpThreads == nullptr)
     241           2 :             pszWarpThreads = CPLGetConfigOption("GDAL_NUM_THREADS", "1");
     242           2 :         if (EQUAL(pszWarpThreads, "ALL_CPUS"))
     243           0 :             nThreads = CPLGetNumCPUs();
     244             :         else
     245           2 :             nThreads = atoi(pszWarpThreads);
     246             :     }
     247             : 
     248          19 :     if (nThreads > 1)
     249             :     {
     250             :         // Compute direct and reverse transforms in parallel.
     251             :         CPLJoinableThread *hThread =
     252           0 :             CPLCreateJoinableThread(GDALTPSComputeForwardInThread, psInfo);
     253           0 :         psInfo->bReverseSolved = psInfo->poReverse->solve() != 0;
     254           0 :         if (hThread != nullptr)
     255           0 :             CPLJoinThread(hThread);
     256             :         else
     257           0 :             psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
     258             :     }
     259             :     else
     260             :     {
     261          19 :         psInfo->bForwardSolved = psInfo->poForward->solve() != 0;
     262          19 :         psInfo->bReverseSolved = psInfo->poReverse->solve() != 0;
     263             :     }
     264             : 
     265          19 :     if (!psInfo->bForwardSolved || !psInfo->bReverseSolved)
     266             :     {
     267           2 :         GDALDestroyTPSTransformer(psInfo);
     268           2 :         return nullptr;
     269             :     }
     270             : 
     271          17 :     return psInfo;
     272             : }
     273             : 
     274             : /************************************************************************/
     275             : /*                     GDALDestroyTPSTransformer()                      */
     276             : /************************************************************************/
     277             : 
     278             : /**
     279             :  * Destroy TPS transformer.
     280             :  *
     281             :  * This function is used to destroy information about a GCP based
     282             :  * polynomial transformation created with GDALCreateTPSTransformer().
     283             :  *
     284             :  * @param pTransformArg the transform arg previously returned by
     285             :  * GDALCreateTPSTransformer().
     286             :  */
     287             : 
     288          19 : void GDALDestroyTPSTransformer(void *pTransformArg)
     289             : 
     290             : {
     291          19 :     if (pTransformArg == nullptr)
     292           0 :         return;
     293             : 
     294          19 :     TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
     295             : 
     296          19 :     if (CPLAtomicDec(&(psInfo->nRefCount)) == 0)
     297             :     {
     298          19 :         delete psInfo->poForward;
     299          19 :         delete psInfo->poReverse;
     300             : 
     301          19 :         delete psInfo;
     302             :     }
     303             : }
     304             : 
     305             : /************************************************************************/
     306             : /*                          GDALTPSTransform()                          */
     307             : /************************************************************************/
     308             : 
     309             : /**
     310             :  * Transforms point based on GCP derived polynomial model.
     311             :  *
     312             :  * This function matches the GDALTransformerFunc signature, and can be
     313             :  * used to transform one or more points from pixel/line coordinates to
     314             :  * georeferenced coordinates (SrcToDst) or vice versa (DstToSrc).
     315             :  *
     316             :  * @param pTransformArg return value from GDALCreateTPSTransformer().
     317             :  * @param bDstToSrc TRUE if transformation is from the destination
     318             :  * (georeferenced) coordinates to pixel/line or FALSE when transforming
     319             :  * from pixel/line to georeferenced coordinates.
     320             :  * @param nPointCount the number of values in the x, y and z arrays.
     321             :  * @param x array containing the X values to be transformed.
     322             :  * @param y array containing the Y values to be transformed.
     323             :  * @param z array containing the Z values to be transformed.
     324             :  * @param panSuccess array in which a flag indicating success (TRUE) or
     325             :  * failure (FALSE) of the transformation are placed.
     326             :  *
     327             :  * @return TRUE.
     328             :  */
     329             : 
     330       34790 : int GDALTPSTransform(void *pTransformArg, int bDstToSrc, int nPointCount,
     331             :                      double *x, double *y, CPL_UNUSED double *z,
     332             :                      int *panSuccess)
     333             : {
     334       34790 :     VALIDATE_POINTER1(pTransformArg, "GDALTPSTransform", 0);
     335             : 
     336       34790 :     TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
     337             : 
     338      112088 :     for (int i = 0; i < nPointCount; i++)
     339             :     {
     340       77298 :         double xy_out[2] = {0.0, 0.0};
     341             : 
     342       77298 :         if (bDstToSrc)
     343             :         {
     344             :             // Compute initial guess
     345       70066 :             psInfo->poReverse->get_point(x[i], y[i], xy_out);
     346             : 
     347      461315 :             const auto ForwardTransformer = [](double xIn, double yIn,
     348             :                                                double &xOut, double &yOut,
     349             :                                                void *pUserData)
     350             :             {
     351      461315 :                 double xyOut[2] = {0.0, 0.0};
     352      461315 :                 TPSTransformInfo *l_psInfo =
     353             :                     static_cast<TPSTransformInfo *>(pUserData);
     354      461315 :                 l_psInfo->poForward->get_point(xIn, yIn, xyOut);
     355      461315 :                 xOut = xyOut[0];
     356      461315 :                 yOut = xyOut[1];
     357      461315 :                 return true;
     358             :             };
     359             : 
     360             :             // Refine the initial guess
     361       70066 :             GDALGenericInverse2D(
     362       70066 :                 x[i], y[i], xy_out[0], xy_out[1], ForwardTransformer, psInfo,
     363             :                 xy_out[0], xy_out[1],
     364             :                 /* computeJacobianMatrixOnlyAtFirstIter = */ true,
     365             :                 /* toleranceOnOutputCoordinates = */ 0,
     366             :                 psInfo->dfSrcApproxErrorReverse);
     367       70066 :             x[i] = xy_out[0];
     368       70066 :             y[i] = xy_out[1];
     369             :         }
     370             :         else
     371             :         {
     372        7232 :             psInfo->poForward->get_point(x[i], y[i], xy_out);
     373        7232 :             x[i] = xy_out[0];
     374        7232 :             y[i] = xy_out[1];
     375             :         }
     376       77298 :         panSuccess[i] = TRUE;
     377             :     }
     378             : 
     379       34790 :     return TRUE;
     380             : }
     381             : 
     382             : /************************************************************************/
     383             : /*                    GDALSerializeTPSTransformer()                     */
     384             : /************************************************************************/
     385             : 
     386           2 : CPLXMLNode *GDALSerializeTPSTransformer(void *pTransformArg)
     387             : 
     388             : {
     389           2 :     VALIDATE_POINTER1(pTransformArg, "GDALSerializeTPSTransformer", nullptr);
     390             : 
     391           2 :     TPSTransformInfo *psInfo = static_cast<TPSTransformInfo *>(pTransformArg);
     392             : 
     393             :     CPLXMLNode *psTree =
     394           2 :         CPLCreateXMLNode(nullptr, CXT_Element, "TPSTransformer");
     395             : 
     396             :     /* -------------------------------------------------------------------- */
     397             :     /*      Serialize bReversed.                                            */
     398             :     /* -------------------------------------------------------------------- */
     399           2 :     CPLCreateXMLElementAndValue(
     400             :         psTree, "Reversed",
     401           4 :         CPLString().Printf("%d", static_cast<int>(psInfo->bReversed)));
     402             : 
     403             :     /* -------------------------------------------------------------------- */
     404             :     /*      Attach GCP List.                                                */
     405             :     /* -------------------------------------------------------------------- */
     406           2 :     if (!psInfo->asGCPs.empty())
     407             :     {
     408           2 :         GDALSerializeGCPListToXML(psTree, psInfo->asGCPs, nullptr);
     409             :     }
     410             : 
     411           2 :     if (psInfo->dfSrcApproxErrorReverse > 0)
     412             :     {
     413           2 :         CPLCreateXMLElementAndValue(
     414             :             psTree, "SrcApproxErrorInPixel",
     415           4 :             CPLString().Printf("%g", psInfo->dfSrcApproxErrorReverse));
     416             :     }
     417             : 
     418           2 :     return psTree;
     419             : }
     420             : 
     421             : /************************************************************************/
     422             : /*                   GDALDeserializeTPSTransformer()                    */
     423             : /************************************************************************/
     424             : 
     425           3 : void *GDALDeserializeTPSTransformer(CPLXMLNode *psTree)
     426             : 
     427             : {
     428             :     /* -------------------------------------------------------------------- */
     429             :     /*      Check for GCPs.                                                 */
     430             :     /* -------------------------------------------------------------------- */
     431           3 :     CPLXMLNode *psGCPList = CPLGetXMLNode(psTree, "GCPList");
     432             : 
     433           6 :     std::vector<gdal::GCP> asGCPs;
     434           3 :     if (psGCPList != nullptr)
     435             :     {
     436           3 :         GDALDeserializeGCPListFromXML(psGCPList, asGCPs, nullptr);
     437             :     }
     438             : 
     439             :     /* -------------------------------------------------------------------- */
     440             :     /*      Get other flags.                                                */
     441             :     /* -------------------------------------------------------------------- */
     442           3 :     const int bReversed = atoi(CPLGetXMLValue(psTree, "Reversed", "0"));
     443             : 
     444           3 :     CPLStringList aosOptions;
     445             :     aosOptions.SetNameValue(
     446             :         "SRC_APPROX_ERROR_IN_PIXEL",
     447           3 :         CPLGetXMLValue(psTree, "SrcApproxErrorInPixel", nullptr));
     448             : 
     449             :     /* -------------------------------------------------------------------- */
     450             :     /*      Generate transformation.                                        */
     451             :     /* -------------------------------------------------------------------- */
     452           3 :     void *pResult = GDALCreateTPSTransformerInt(static_cast<int>(asGCPs.size()),
     453             :                                                 gdal::GCP::c_ptr(asGCPs),
     454             :                                                 bReversed, aosOptions.List());
     455             : 
     456           6 :     return pResult;
     457             : }