Branch data     Line data    Source code

       1                 :            : /***************************************************************************
       2                 :            :     qgsgcptransformer.cpp
       3                 :            :      --------------------------------------
       4                 :            :     Date                 : 18-Feb-2009
       5                 :            :     Copyright            : (c) 2009 by Manuel Massing
       6                 :            :     Email                : m.massing at warped-space.de
       7                 :            :  ***************************************************************************
       8                 :            :  *                                                                         *
       9                 :            :  *   This program is free software; you can redistribute it and/or modify  *
      10                 :            :  *   it under the terms of the GNU General Public License as published by  *
      11                 :            :  *   the Free Software Foundation; either version 2 of the License, or     *
      12                 :            :  *   (at your option) any later version.                                   *
      13                 :            :  *                                                                         *
      14                 :            :  ***************************************************************************/
      15                 :            : 
      16                 :            : #include "qgsgcptransformer.h"
      17                 :            : 
      18                 :            : #include <gdal.h>
      19                 :            : #include <gdal_alg.h>
      20                 :            : 
      21                 :            : #include "qgsleastsquares.h"
      22                 :            : 
      23                 :            : #include <cmath>
      24                 :            : 
      25                 :            : #include <cassert>
      26                 :            : #include <limits>
      27                 :            : 
      28                 :            : 
      29                 :          0 : bool QgsGcpTransformerInterface::transform( double &x, double &y, bool inverseTransform ) const
      30                 :            : {
      31                 :          0 :   GDALTransformerFunc t = GDALTransformer();
      32                 :            :   // Fail if no transformer function was returned
      33                 :          0 :   if ( !t )
      34                 :          0 :     return false;
      35                 :            : 
      36                 :          0 :   double z = 0.0;
      37                 :          0 :   int success = 0;
      38                 :            : 
      39                 :            :   // Call GDAL transform function
      40                 :          0 :   ( *t )( GDALTransformerArgs(), inverseTransform ? 1 : 0, 1,  &x, &y, &z, &success );
      41                 :          0 :   if ( !success )
      42                 :          0 :     return false;
      43                 :            : 
      44                 :          0 :   return true;
      45                 :          0 : }
      46                 :            : 
      47                 :          0 : QString QgsGcpTransformerInterface::methodToString( QgsGcpTransformerInterface::TransformMethod method )
      48                 :            : {
      49                 :          0 :   switch ( method )
      50                 :            :   {
      51                 :            :     case QgsGcpTransformerInterface::TransformMethod::Linear:
      52                 :          0 :       return QObject::tr( "Linear" );
      53                 :            :     case QgsGcpTransformerInterface::TransformMethod::Helmert:
      54                 :          0 :       return QObject::tr( "Helmert" );
      55                 :            :     case QgsGcpTransformerInterface::TransformMethod::PolynomialOrder1:
      56                 :          0 :       return QObject::tr( "Polynomial 1" );
      57                 :            :     case QgsGcpTransformerInterface::TransformMethod::PolynomialOrder2:
      58                 :          0 :       return QObject::tr( "Polynomial 2" );
      59                 :            :     case QgsGcpTransformerInterface::TransformMethod::PolynomialOrder3:
      60                 :          0 :       return QObject::tr( "Polynomial 3" );
      61                 :            :     case QgsGcpTransformerInterface::TransformMethod::ThinPlateSpline:
      62                 :          0 :       return QObject::tr( "Thin Plate Spline (TPS)" );
      63                 :            :     case QgsGcpTransformerInterface::TransformMethod::Projective:
      64                 :          0 :       return QObject::tr( "Projective" );
      65                 :            :     default:
      66                 :          0 :       return QObject::tr( "Not set" );
      67                 :            :   }
      68                 :          0 : }
      69                 :            : 
      70                 :          0 : QgsGcpTransformerInterface *QgsGcpTransformerInterface::create( QgsGcpTransformerInterface::TransformMethod method )
      71                 :            : {
      72                 :          0 :   switch ( method )
      73                 :            :   {
      74                 :            :     case TransformMethod::Linear:
      75                 :          0 :       return new QgsLinearGeorefTransform;
      76                 :            :     case TransformMethod::Helmert:
      77                 :          0 :       return new QgsHelmertGeorefTransform;
      78                 :            :     case TransformMethod::PolynomialOrder1:
      79                 :          0 :       return new QgsGDALGeorefTransform( false, 1 );
      80                 :            :     case TransformMethod::PolynomialOrder2:
      81                 :          0 :       return new QgsGDALGeorefTransform( false, 2 );
      82                 :            :     case TransformMethod::PolynomialOrder3:
      83                 :          0 :       return new QgsGDALGeorefTransform( false, 3 );
      84                 :            :     case TransformMethod::ThinPlateSpline:
      85                 :          0 :       return new QgsGDALGeorefTransform( true, 0 );
      86                 :            :     case TransformMethod::Projective:
      87                 :          0 :       return new QgsProjectiveGeorefTransform;
      88                 :            :     default:
      89                 :          0 :       return nullptr;
      90                 :            :   }
      91                 :          0 : }
      92                 :            : 
      93                 :          0 : QgsGcpTransformerInterface *QgsGcpTransformerInterface::createFromParameters( QgsGcpTransformerInterface::TransformMethod method, const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates )
      94                 :            : {
      95                 :          0 :   std::unique_ptr< QgsGcpTransformerInterface > transformer( create( method ) );
      96                 :          0 :   if ( !transformer )
      97                 :          0 :     return nullptr;
      98                 :            : 
      99                 :          0 :   if ( !transformer->updateParametersFromGcps( sourceCoordinates, destinationCoordinates ) )
     100                 :          0 :     return nullptr;
     101                 :            : 
     102                 :          0 :   return transformer.release();
     103                 :          0 : }
     104                 :            : 
     105                 :            : 
     106                 :            : //
     107                 :            : // QgsLinearGeorefTransform
     108                 :            : //
     109                 :            : 
     110                 :          0 : bool QgsLinearGeorefTransform::getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const
     111                 :            : {
     112                 :          0 :   origin = mParameters.origin;
     113                 :          0 :   scaleX = mParameters.scaleX;
     114                 :          0 :   scaleY = mParameters.scaleY;
     115                 :          0 :   return true;
     116                 :            : }
     117                 :            : 
     118                 :          0 : QgsGcpTransformerInterface *QgsLinearGeorefTransform::clone() const
     119                 :            : {
     120                 :          0 :   std::unique_ptr< QgsLinearGeorefTransform > res = std::make_unique< QgsLinearGeorefTransform >();
     121                 :          0 :   res->mParameters = mParameters;
     122                 :          0 :   return res.release();
     123                 :          0 : }
     124                 :            : 
     125                 :          0 : bool QgsLinearGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
     126                 :            : {
     127                 :          0 :   if ( destinationCoordinates.size() < minimumGcpCount() )
     128                 :          0 :     return false;
     129                 :            : 
     130                 :          0 :   mParameters.invertYAxis = invertYAxis;
     131                 :          0 :   QgsLeastSquares::linear( sourceCoordinates, destinationCoordinates, mParameters.origin, mParameters.scaleX, mParameters.scaleY );
     132                 :          0 :   return true;
     133                 :          0 : }
     134                 :            : 
     135                 :          0 : int QgsLinearGeorefTransform::minimumGcpCount() const
     136                 :            : {
     137                 :          0 :   return 2;
     138                 :            : }
     139                 :            : 
     140                 :          0 : GDALTransformerFunc QgsLinearGeorefTransform::GDALTransformer() const
     141                 :            : {
     142                 :          0 :   return QgsLinearGeorefTransform::linearTransform;
     143                 :            : }
     144                 :            : 
     145                 :          0 : void *QgsLinearGeorefTransform::GDALTransformerArgs() const
     146                 :            : {
     147                 :          0 :   return ( void * )&mParameters;
     148                 :            : }
     149                 :            : 
     150                 :          0 : QgsGcpTransformerInterface::TransformMethod QgsLinearGeorefTransform::method() const
     151                 :            : {
     152                 :          0 :   return TransformMethod::Linear;
     153                 :            : }
     154                 :            : 
     155                 :          0 : int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
     156                 :            :     double *x, double *y, double *z, int *panSuccess )
     157                 :            : {
     158                 :            :   Q_UNUSED( z )
     159                 :          0 :   LinearParameters *t = static_cast<LinearParameters *>( pTransformerArg );
     160                 :          0 :   if ( !t )
     161                 :          0 :     return false;
     162                 :            : 
     163                 :          0 :   if ( !bDstToSrc )
     164                 :            :   {
     165                 :          0 :     for ( int i = 0; i < nPointCount; ++i )
     166                 :            :     {
     167                 :          0 :       x[i] = x[i] * t->scaleX + t->origin.x();
     168                 :          0 :       y[i] = ( t->invertYAxis ? -1 : 1 ) * y[i] * t->scaleY + t->origin.y();
     169                 :          0 :       panSuccess[i] = true;
     170                 :          0 :     }
     171                 :          0 :   }
     172                 :            :   else
     173                 :            :   {
     174                 :            :     // Guard against division by zero
     175                 :          0 :     if ( std::fabs( t->scaleX ) < std::numeric_limits<double>::epsilon() ||
     176                 :          0 :          std::fabs( t->scaleY ) < std::numeric_limits<double>::epsilon() )
     177                 :            :     {
     178                 :          0 :       for ( int i = 0; i < nPointCount; ++i )
     179                 :            :       {
     180                 :          0 :         panSuccess[i] = false;
     181                 :          0 :       }
     182                 :          0 :       return false;
     183                 :            :     }
     184                 :          0 :     for ( int i = 0; i < nPointCount; ++i )
     185                 :            :     {
     186                 :          0 :       x[i] = ( x[i] - t->origin.x() ) / t->scaleX;
     187                 :          0 :       y[i] = ( y[i] - t->origin.y() ) / ( ( t->invertYAxis ? -1 : 1 ) * t->scaleY );
     188                 :          0 :       panSuccess[i] = true;
     189                 :          0 :     }
     190                 :            :   }
     191                 :            : 
     192                 :          0 :   return true;
     193                 :          0 : }
     194                 :            : 
     195                 :            : //
     196                 :            : // QgsHelmertGeorefTransform
     197                 :            : //
     198                 :          0 : bool QgsHelmertGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
     199                 :            : {
     200                 :          0 :   if ( destinationCoordinates.size() < minimumGcpCount() )
     201                 :          0 :     return false;
     202                 :            : 
     203                 :          0 :   mHelmertParameters.invertYAxis = invertYAxis;
     204                 :          0 :   QgsLeastSquares::helmert( sourceCoordinates, destinationCoordinates, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle );
     205                 :          0 :   return true;
     206                 :          0 : }
     207                 :            : 
     208                 :          0 : int QgsHelmertGeorefTransform::minimumGcpCount() const
     209                 :            : {
     210                 :          0 :   return 2;
     211                 :            : }
     212                 :            : 
     213                 :          0 : GDALTransformerFunc QgsHelmertGeorefTransform::GDALTransformer() const
     214                 :            : {
     215                 :          0 :   return QgsHelmertGeorefTransform::helmertTransform;
     216                 :            : }
     217                 :            : 
     218                 :          0 : void *QgsHelmertGeorefTransform::GDALTransformerArgs() const
     219                 :            : {
     220                 :          0 :   return ( void * )&mHelmertParameters;
     221                 :            : }
     222                 :            : 
     223                 :          0 : QgsGcpTransformerInterface::TransformMethod QgsHelmertGeorefTransform::method() const
     224                 :            : {
     225                 :          0 :   return TransformMethod::Helmert;
     226                 :            : }
     227                 :            : 
     228                 :          0 : bool QgsHelmertGeorefTransform::getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const
     229                 :            : {
     230                 :          0 :   origin = mHelmertParameters.origin;
     231                 :          0 :   scale = mHelmertParameters.scale;
     232                 :          0 :   rotation = mHelmertParameters.angle;
     233                 :          0 :   return true;
     234                 :            : }
     235                 :            : 
     236                 :          0 : QgsGcpTransformerInterface *QgsHelmertGeorefTransform::clone() const
     237                 :            : {
     238                 :          0 :   std::unique_ptr< QgsHelmertGeorefTransform > res = std::make_unique< QgsHelmertGeorefTransform >();
     239                 :          0 :   res->mHelmertParameters = mHelmertParameters;
     240                 :          0 :   return res.release();
     241                 :          0 : }
     242                 :            : 
     243                 :          0 : int QgsHelmertGeorefTransform::helmertTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
     244                 :            :     double *x, double *y, double *z, int *panSuccess )
     245                 :            : {
     246                 :            :   Q_UNUSED( z )
     247                 :          0 :   const HelmertParameters *t = static_cast< const HelmertParameters *>( pTransformerArg );
     248                 :          0 :   if ( !t )
     249                 :          0 :     return false;
     250                 :            : 
     251                 :          0 :   double a = std::cos( t->angle );
     252                 :          0 :   double b = std::sin( t->angle );
     253                 :          0 :   const double x0 = t->origin.x();
     254                 :          0 :   const double y0 = t->origin.y();
     255                 :          0 :   const double s = t->scale;
     256                 :          0 :   if ( !bDstToSrc )
     257                 :            :   {
     258                 :          0 :     a *= s;
     259                 :          0 :     b *= s;
     260                 :          0 :     for ( int i = 0; i < nPointCount; ++i )
     261                 :            :     {
     262                 :          0 :       const double xT = x[i];
     263                 :          0 :       const double yT = y[i];
     264                 :            : 
     265                 :          0 :       if ( t->invertYAxis )
     266                 :            :       {
     267                 :            :         // Because rotation parameters where estimated in a CS with negative y-axis ^= down.
     268                 :            :         // we need to apply the rotation matrix and a change of base:
     269                 :            :         // |cos a, -sin a | |1, 0|   | cos a,  sin a|
     270                 :            :         // |sin a,  cos a | |0,-1| = | sin a, -cos a|
     271                 :          0 :         x[i] = x0 + ( a * xT + b * yT );
     272                 :          0 :         y[i] = y0 + ( b * xT - a * yT );
     273                 :          0 :       }
     274                 :            :       else
     275                 :            :       {
     276                 :          0 :         x[i] = x0 + ( a * xT - b * yT );
     277                 :          0 :         y[i] = y0 + ( b * xT + a * yT );
     278                 :            :       }
     279                 :          0 :       panSuccess[i] = true;
     280                 :          0 :     }
     281                 :          0 :   }
     282                 :            :   else
     283                 :            :   {
     284                 :            :     // Guard against division by zero
     285                 :          0 :     if ( std::fabs( s ) < std::numeric_limits<double>::epsilon() )
     286                 :            :     {
     287                 :          0 :       for ( int i = 0; i < nPointCount; ++i )
     288                 :            :       {
     289                 :          0 :         panSuccess[i] = false;
     290                 :          0 :       }
     291                 :          0 :       return false;
     292                 :            :     }
     293                 :          0 :     a /= s;
     294                 :          0 :     b /= s;
     295                 :          0 :     for ( int i = 0; i < nPointCount; ++i )
     296                 :            :     {
     297                 :          0 :       const double xT = x[i] - x0;
     298                 :          0 :       const double yT = y[i] - y0;
     299                 :          0 :       if ( t->invertYAxis )
     300                 :            :       {
     301                 :            :         // | cos a,  sin a |^-1   |cos a,  sin a|
     302                 :            :         // | sin a, -cos a |    = |sin a, -cos a|
     303                 :          0 :         x[i] = a * xT + b * yT;
     304                 :          0 :         y[i] = b * xT - a * yT;
     305                 :          0 :       }
     306                 :            :       else
     307                 :            :       {
     308                 :          0 :         x[i] = a * xT + b * yT;
     309                 :          0 :         y[i] = -b * xT + a * yT;
     310                 :            :       }
     311                 :          0 :       panSuccess[i] = true;
     312                 :          0 :     }
     313                 :            :   }
     314                 :          0 :   return true;
     315                 :          0 : }
     316                 :            : 
     317                 :            : //
     318                 :            : // QgsGDALGeorefTransform
     319                 :            : //
     320                 :            : 
     321                 :          0 : QgsGDALGeorefTransform::QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder )
     322                 :          0 :   : mPolynomialOrder( std::min( 3u, polynomialOrder ) )
     323                 :          0 :   , mIsTPSTransform( useTPS )
     324                 :          0 : {
     325                 :          0 : }
     326                 :            : 
     327                 :          0 : QgsGDALGeorefTransform::~QgsGDALGeorefTransform()
     328                 :          0 : {
     329                 :          0 :   destroyGdalArgs();
     330                 :          0 : }
     331                 :            : 
     332                 :          0 : QgsGcpTransformerInterface *QgsGDALGeorefTransform::clone() const
     333                 :            : {
     334                 :          0 :   std::unique_ptr< QgsGDALGeorefTransform > res = std::make_unique< QgsGDALGeorefTransform >( mIsTPSTransform, mPolynomialOrder );
     335                 :          0 :   res->updateParametersFromGcps( mSourceCoords, mDestCoordinates, mInvertYAxis );
     336                 :          0 :   return res.release();
     337                 :          0 : }
     338                 :            : 
     339                 :          0 : bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
     340                 :            : {
     341                 :          0 :   mSourceCoords = sourceCoordinates;
     342                 :          0 :   mDestCoordinates = destinationCoordinates;
     343                 :          0 :   mInvertYAxis = invertYAxis;
     344                 :            : 
     345                 :          0 :   assert( sourceCoordinates.size() == destinationCoordinates.size() );
     346                 :          0 :   if ( sourceCoordinates.size() != destinationCoordinates.size() )
     347                 :          0 :     return false;
     348                 :          0 :   int n = sourceCoordinates.size();
     349                 :            : 
     350                 :          0 :   GDAL_GCP *GCPList = new GDAL_GCP[n];
     351                 :          0 :   for ( int i = 0; i < n; i++ )
     352                 :            :   {
     353                 :          0 :     GCPList[i].pszId = new char[20];
     354                 :          0 :     snprintf( GCPList[i].pszId, 19, "gcp%i", i );
     355                 :          0 :     GCPList[i].pszInfo = nullptr;
     356                 :          0 :     GCPList[i].dfGCPPixel = sourceCoordinates[i].x();
     357                 :          0 :     GCPList[i].dfGCPLine  = ( mInvertYAxis ? -1 : 1 ) * sourceCoordinates[i].y();
     358                 :          0 :     GCPList[i].dfGCPX = destinationCoordinates[i].x();
     359                 :          0 :     GCPList[i].dfGCPY = destinationCoordinates[i].y();
     360                 :          0 :     GCPList[i].dfGCPZ = 0;
     361                 :          0 :   }
     362                 :          0 :   destroyGdalArgs();
     363                 :            : 
     364                 :          0 :   if ( mIsTPSTransform )
     365                 :          0 :     mGDALTransformerArgs = GDALCreateTPSTransformer( n, GCPList, false );
     366                 :            :   else
     367                 :          0 :     mGDALTransformerArgs = GDALCreateGCPTransformer( n, GCPList, mPolynomialOrder, false );
     368                 :            : 
     369                 :          0 :   for ( int i = 0; i < n; i++ )
     370                 :            :   {
     371                 :          0 :     delete [] GCPList[i].pszId;
     372                 :          0 :   }
     373                 :          0 :   delete [] GCPList;
     374                 :            : 
     375                 :          0 :   return nullptr != mGDALTransformerArgs;
     376                 :          0 : }
     377                 :            : 
     378                 :          0 : int QgsGDALGeorefTransform::minimumGcpCount() const
     379                 :            : {
     380                 :          0 :   if ( mIsTPSTransform )
     381                 :          0 :     return 1;
     382                 :            :   else
     383                 :          0 :     return ( ( mPolynomialOrder + 2 ) * ( mPolynomialOrder + 1 ) ) / 2;
     384                 :          0 : }
     385                 :            : 
     386                 :          0 : GDALTransformerFunc QgsGDALGeorefTransform::GDALTransformer() const
     387                 :            : {
     388                 :            :   // Fail if no arguments were calculated through updateParametersFromGCP
     389                 :          0 :   if ( !mGDALTransformerArgs )
     390                 :          0 :     return nullptr;
     391                 :            : 
     392                 :          0 :   if ( mIsTPSTransform )
     393                 :          0 :     return GDALTPSTransform;
     394                 :            :   else
     395                 :          0 :     return GDALGCPTransform;
     396                 :          0 : }
     397                 :            : 
     398                 :          0 : void *QgsGDALGeorefTransform::GDALTransformerArgs() const
     399                 :            : {
     400                 :          0 :   return mGDALTransformerArgs;
     401                 :            : }
     402                 :            : 
     403                 :          0 : QgsGcpTransformerInterface::TransformMethod QgsGDALGeorefTransform::method() const
     404                 :            : {
     405                 :          0 :   if ( mIsTPSTransform )
     406                 :          0 :     return TransformMethod::ThinPlateSpline;
     407                 :            : 
     408                 :          0 :   switch ( mPolynomialOrder )
     409                 :            :   {
     410                 :            :     case 1:
     411                 :          0 :       return TransformMethod::PolynomialOrder1;
     412                 :            :     case 2:
     413                 :          0 :       return TransformMethod::PolynomialOrder2;
     414                 :            :     case 3:
     415                 :          0 :       return TransformMethod::PolynomialOrder3;
     416                 :            :   }
     417                 :          0 :   return TransformMethod::InvalidTransform;
     418                 :          0 : }
     419                 :            : 
     420                 :          0 : void QgsGDALGeorefTransform::destroyGdalArgs()
     421                 :            : {
     422                 :          0 :   if ( mGDALTransformerArgs )
     423                 :            :   {
     424                 :          0 :     if ( mIsTPSTransform )
     425                 :          0 :       GDALDestroyTPSTransformer( mGDALTransformerArgs );
     426                 :            :     else
     427                 :          0 :       GDALDestroyGCPTransformer( mGDALTransformerArgs );
     428                 :          0 :   }
     429                 :          0 : }
     430                 :            : 
     431                 :            : //
     432                 :            : // QgsProjectiveGeorefTransform
     433                 :            : //
     434                 :            : 
     435                 :          0 : QgsProjectiveGeorefTransform::QgsProjectiveGeorefTransform()
     436                 :          0 :   : mParameters()
     437                 :          0 : {}
     438                 :            : 
     439                 :          0 : QgsGcpTransformerInterface *QgsProjectiveGeorefTransform::clone() const
     440                 :            : {
     441                 :          0 :   std::unique_ptr< QgsProjectiveGeorefTransform > res = std::make_unique< QgsProjectiveGeorefTransform >();
     442                 :          0 :   res->mParameters = mParameters;
     443                 :          0 :   return res.release();
     444                 :          0 : }
     445                 :            : 
     446                 :          0 : bool QgsProjectiveGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
     447                 :            : {
     448                 :          0 :   if ( destinationCoordinates.size() < minimumGcpCount() )
     449                 :          0 :     return false;
     450                 :            : 
     451                 :          0 :   if ( invertYAxis )
     452                 :            :   {
     453                 :            :     // HACK: flip y coordinates, because georeferencer and gdal use different conventions
     454                 :          0 :     QVector<QgsPointXY> flippedPixelCoords;
     455                 :          0 :     flippedPixelCoords.reserve( sourceCoordinates.size() );
     456                 :          0 :     for ( const QgsPointXY &coord : sourceCoordinates )
     457                 :            :     {
     458                 :          0 :       flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() );
     459                 :            :     }
     460                 :            : 
     461                 :          0 :     QgsLeastSquares::projective( flippedPixelCoords, destinationCoordinates, mParameters.H );
     462                 :          0 :   }
     463                 :            :   else
     464                 :            :   {
     465                 :          0 :     QgsLeastSquares::projective( sourceCoordinates, destinationCoordinates, mParameters.H );
     466                 :            :   }
     467                 :            : 
     468                 :            :   // Invert the homography matrix using adjoint matrix
     469                 :          0 :   double *H = mParameters.H;
     470                 :            : 
     471                 :            :   double adjoint[9];
     472                 :          0 :   adjoint[0] = H[4] * H[8] - H[5] * H[7];
     473                 :          0 :   adjoint[1] = -H[1] * H[8] + H[2] * H[7];
     474                 :          0 :   adjoint[2] = H[1] * H[5] - H[2] * H[4];
     475                 :            : 
     476                 :          0 :   adjoint[3] = -H[3] * H[8] + H[5] * H[6];
     477                 :          0 :   adjoint[4] = H[0] * H[8] - H[2] * H[6];
     478                 :          0 :   adjoint[5] = -H[0] * H[5] + H[2] * H[3];
     479                 :            : 
     480                 :          0 :   adjoint[6] = H[3] * H[7] - H[4] * H[6];
     481                 :          0 :   adjoint[7] = -H[0] * H[7] + H[1] * H[6];
     482                 :          0 :   adjoint[8] = H[0] * H[4] - H[1] * H[3];
     483                 :            : 
     484                 :          0 :   double det = H[0] * adjoint[0] + H[3] * adjoint[1] + H[6] * adjoint[2];
     485                 :            : 
     486                 :          0 :   if ( std::fabs( det ) < 1024.0 * std::numeric_limits<double>::epsilon() )
     487                 :            :   {
     488                 :          0 :     mParameters.hasInverse = false;
     489                 :          0 :   }
     490                 :            :   else
     491                 :            :   {
     492                 :          0 :     mParameters.hasInverse = true;
     493                 :          0 :     double oo_det = 1.0 / det;
     494                 :          0 :     for ( int i = 0; i < 9; i++ )
     495                 :            :     {
     496                 :          0 :       mParameters.Hinv[i] = adjoint[i] * oo_det;
     497                 :          0 :     }
     498                 :            :   }
     499                 :          0 :   return true;
     500                 :          0 : }
     501                 :            : 
     502                 :          0 : int QgsProjectiveGeorefTransform::minimumGcpCount() const
     503                 :            : {
     504                 :          0 :   return 4;
     505                 :            : }
     506                 :            : 
     507                 :          0 : GDALTransformerFunc QgsProjectiveGeorefTransform::GDALTransformer() const
     508                 :            : {
     509                 :          0 :   return QgsProjectiveGeorefTransform::projectiveTransform;
     510                 :            : }
     511                 :            : 
     512                 :          0 : void *QgsProjectiveGeorefTransform::GDALTransformerArgs() const
     513                 :            : {
     514                 :          0 :   return ( void * )&mParameters;
     515                 :            : }
     516                 :            : 
     517                 :          0 : QgsGcpTransformerInterface::TransformMethod QgsProjectiveGeorefTransform::method() const
     518                 :            : {
     519                 :          0 :   return TransformMethod::Projective;
     520                 :            : }
     521                 :            : 
     522                 :          0 : int QgsProjectiveGeorefTransform::projectiveTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
     523                 :            :     double *x, double *y, double *z, int *panSuccess )
     524                 :            : {
     525                 :            :   Q_UNUSED( z )
     526                 :          0 :   ProjectiveParameters *t = static_cast<ProjectiveParameters *>( pTransformerArg );
     527                 :          0 :   if ( !t )
     528                 :          0 :     return false;
     529                 :            : 
     530                 :          0 :   double *H = nullptr;
     531                 :          0 :   if ( !bDstToSrc )
     532                 :            :   {
     533                 :          0 :     H = t->H;
     534                 :          0 :   }
     535                 :            :   else
     536                 :            :   {
     537                 :          0 :     if ( !t->hasInverse )
     538                 :            :     {
     539                 :          0 :       for ( int i = 0; i < nPointCount; ++i )
     540                 :            :       {
     541                 :          0 :         panSuccess[i] = false;
     542                 :          0 :       }
     543                 :          0 :       return false;
     544                 :            :     }
     545                 :          0 :     H = t->Hinv;
     546                 :            :   }
     547                 :            : 
     548                 :            : 
     549                 :          0 :   for ( int i = 0; i < nPointCount; ++i )
     550                 :            :   {
     551                 :          0 :     double Z = x[i] * H[6] + y[i] * H[7] + H[8];
     552                 :            :     // Projects to infinity?
     553                 :          0 :     if ( std::fabs( Z ) < 1024.0 * std::numeric_limits<double>::epsilon() )
     554                 :            :     {
     555                 :          0 :       panSuccess[i] = false;
     556                 :          0 :       continue;
     557                 :            :     }
     558                 :          0 :     double X = ( x[i] * H[0] + y[i] * H[1] + H[2] ) / Z;
     559                 :          0 :     double Y = ( x[i] * H[3] + y[i] * H[4] + H[5] ) / Z;
     560                 :            : 
     561                 :          0 :     x[i] = X;
     562                 :          0 :     y[i] = Y;
     563                 :            : 
     564                 :          0 :     panSuccess[i] = true;
     565                 :          0 :   }
     566                 :            : 
     567                 :          0 :   return true;
     568                 :          0 : }