Move classes representing GCP based transformations from app to analysis,

and make ready for exposure to public API

doc/CMakeLists.txt

@@ -110,6 +110,7 @@ if(WITH_APIDOC)
       ${CMAKE_SOURCE_DIR}/src/gui/vectortile
       ${CMAKE_SOURCE_DIR}/src/analysis
       ${CMAKE_SOURCE_DIR}/src/analysis/mesh
+      ${CMAKE_SOURCE_DIR}/src/analysis/georeferencing
       ${CMAKE_SOURCE_DIR}/src/analysis/interpolation
       ${CMAKE_SOURCE_DIR}/src/analysis/network
       ${CMAKE_SOURCE_DIR}/src/analysis/processing

python/CMakeLists.txt

@@ -245,6 +245,7 @@ if(WITH_ANALYSIS)
   include_directories(BEFORE
 
     ${CMAKE_BINARY_DIR}/src/analysis/processing
+    ${CMAKE_BINARY_DIR}/src/analysis/georeferencing
     ${CMAKE_BINARY_DIR}/src/analysis/vector
     ${CMAKE_BINARY_DIR}/src/analysis/mesh
     ${CMAKE_BINARY_DIR}/src/analysis/raster

python/analysis/analysis_auto.sip

@@ -1,5 +1,6 @@
 // Include auto-generated SIP files
 %Include auto_generated/qgsanalysis.sip
+%Include auto_generated/georeferencing/qgsgcptransformer.sip
 %Include auto_generated/interpolation/qgsgridfilewriter.sip
 %Include auto_generated/interpolation/qgsidwinterpolator.sip
 %Include auto_generated/interpolation/qgsinterpolator.sip

python/analysis/auto_generated/georeferencing/qgsgcptransformer.sip.in

@@ -0,0 +1,60 @@
+/************************************************************************
+ * This file has been generated automatically from                      *
+ *                                                                      *
+ * src/analysis/georeferencing/qgsgcptransformer.h                      *
+ *                                                                      *
+ * Do not edit manually ! Edit header and run scripts/sipify.pl again   *
+ ************************************************************************/
+
+
+
+class QgsGcpTransformerInterface /Abstract/
+{
+%Docstring
+An interface for Ground Control Points (GCP) based transformations.
+
+:py:class:`QgsGcpTransformerInterface` implementations are able to transform point locations
+based on a transformation method and a list of GCPs.
+
+.. versionadded:: 3.20
+%End
+
+%TypeHeaderCode
+#include "qgsgcptransformer.h"
+%End
+  public:
+
+    QgsGcpTransformerInterface();
+
+    virtual ~QgsGcpTransformerInterface();
+
+
+
+    virtual bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoordinates, const QVector<QgsPointXY> &pixelCoordinates ) = 0;
+%Docstring
+Fits transformation parameters using the specified Ground Control Points (GCPs) lists of map coordinates and pixel coordinates.
+
+:return: ``True`` on success, ``False`` on failure
+%End
+
+    virtual int minimumGcpCount() const = 0;
+%Docstring
+Returns the minimum number of Ground Control Points (GCPs) required for parameter fitting.
+%End
+
+
+  private:
+    QgsGcpTransformerInterface( const QgsGcpTransformerInterface &other );
+};
+
+
+
+
+
+/************************************************************************
+ * This file has been generated automatically from                      *
+ *                                                                      *
+ * src/analysis/georeferencing/qgsgcptransformer.h                      *
+ *                                                                      *
+ * Do not edit manually ! Edit header and run scripts/sipify.pl again   *
+ ************************************************************************/

src/analysis/CMakeLists.txt

@@ -3,6 +3,8 @@
 
 set(QGIS_ANALYSIS_SRCS
   qgsanalysis.cpp
+  georeferencing/qgsgcptransformer.cpp
+  georeferencing/qgsleastsquares.cpp
 
   interpolation/qgsgridfilewriter.cpp
   interpolation/qgsidwinterpolator.cpp
@@ -283,6 +285,8 @@ set(QGIS_ANALYSIS_SRCS
 set(QGIS_ANALYSIS_HDRS
   qgsanalysis.h
 
+  georeferencing/qgsgcptransformer.h
+
   interpolation/Bezier3D.h
   interpolation/CloughTocherInterpolator.h
   interpolation/qgsdualedgetriangulation.h
@@ -391,6 +395,7 @@ find_package(EXIV2 REQUIRED)
 include_directories(SYSTEM ${SPATIALITE_INCLUDE_DIR})
 include_directories(SYSTEM ${SPATIALINDEX_INCLUDE_DIR})
 include_directories(SYSTEM ${SQLITE3_INCLUDE_DIR})
+include_directories(SYSTEM ${GSL_INCLUDE_DIR})
 include_directories(BEFORE raster)
 include_directories(BEFORE mesh)
 
@@ -437,6 +442,7 @@ add_library(qgis_analysis ${LIBRARY_TYPE} ${QGIS_ANALYSIS_SRCS} ${QGIS_ANALYSIS_
 
 target_include_directories(qgis_analysis PUBLIC
   ${CMAKE_SOURCE_DIR}/src/analysis
+  ${CMAKE_SOURCE_DIR}/src/analysis/georeferencing
   ${CMAKE_SOURCE_DIR}/src/analysis/interpolation
   ${CMAKE_SOURCE_DIR}/src/analysis/mesh
   ${CMAKE_SOURCE_DIR}/src/analysis/network
@@ -497,6 +503,7 @@ target_link_libraries(
   qgis_analysis
   qgis_core
   ${EXIV2_LIBRARY}
+  ${GSL_LIBRARIES}
   )
 
 if(HAVE_OPENCL)

src/analysis/georeferencing/qgsgcptransformer.cpp

@@ -0,0 +1,398 @@
+/***************************************************************************
+    qgsgcptransformer.cpp
+     --------------------------------------
+    Date                 : 18-Feb-2009
+    Copyright            : (c) 2009 by Manuel Massing
+    Email                : m.massing at warped-space.de
+ ***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#include "qgsgcptransformer.h"
+
+#include <gdal.h>
+#include <gdal_alg.h>
+
+#include "qgsleastsquares.h"
+
+#include <cmath>
+
+#include <cassert>
+#include <limits>
+
+//
+// QgsLinearGeorefTransform
+//
+
+bool QgsLinearGeorefTransform::getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const
+{
+  origin = mParameters.origin;
+  scaleX = mParameters.scaleX;
+  scaleY = mParameters.scaleY;
+  return true;
+}
+
+bool QgsLinearGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
+{
+  if ( mapCoords.size() < minimumGcpCount() )
+    return false;
+  QgsLeastSquares::linear( mapCoords, pixelCoords, mParameters.origin, mParameters.scaleX, mParameters.scaleY );
+  return true;
+}
+
+int QgsLinearGeorefTransform::minimumGcpCount() const
+{
+  return 2;
+}
+
+GDALTransformerFunc QgsLinearGeorefTransform::GDALTransformer() const
+{
+  return QgsLinearGeorefTransform::linearTransform;
+}
+
+void *QgsLinearGeorefTransform::GDALTransformerArgs() const
+{
+  return ( void * )&mParameters;
+}
+
+int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
+    double *x, double *y, double *z, int *panSuccess )
+{
+  Q_UNUSED( z )
+  LinearParameters *t = static_cast<LinearParameters *>( pTransformerArg );
+  if ( !t )
+    return false;
+
+  if ( !bDstToSrc )
+  {
+    for ( int i = 0; i < nPointCount; ++i )
+    {
+      x[i] = x[i] * t->scaleX + t->origin.x();
+      y[i] = -y[i] * t->scaleY + t->origin.y();
+      panSuccess[i] = true;
+    }
+  }
+  else
+  {
+    // Guard against division by zero
+    if ( std::fabs( t->scaleX ) < std::numeric_limits<double>::epsilon() ||
+         std::fabs( t->scaleY ) < std::numeric_limits<double>::epsilon() )
+    {
+      for ( int i = 0; i < nPointCount; ++i )
+      {
+        panSuccess[i] = false;
+      }
+      return false;
+    }
+    for ( int i = 0; i < nPointCount; ++i )
+    {
+      x[i] = ( x[i] - t->origin.x() ) / t->scaleX;
+      y[i] = ( y[i] - t->origin.y() ) / ( -t->scaleY );
+      panSuccess[i] = true;
+    }
+  }
+
+  return true;
+}
+
+//
+// QgsHelmertGeorefTransform
+//
+bool QgsHelmertGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
+{
+  if ( mapCoords.size() < minimumGcpCount() )
+    return false;
+
+  QgsLeastSquares::helmert( mapCoords, pixelCoords, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle );
+  return true;
+}
+
+int QgsHelmertGeorefTransform::minimumGcpCount() const
+{
+  return 2;
+}
+
+
+GDALTransformerFunc QgsHelmertGeorefTransform::GDALTransformer() const
+{
+  return QgsHelmertGeorefTransform::helmert_transform;
+}
+
+void *QgsHelmertGeorefTransform::GDALTransformerArgs() const
+{
+  return ( void * )&mHelmertParameters;
+}
+
+bool QgsHelmertGeorefTransform::getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const
+{
+  origin = mHelmertParameters.origin;
+  scale = mHelmertParameters.scale;
+  rotation = mHelmertParameters.angle;
+  return true;
+}
+
+int QgsHelmertGeorefTransform::helmert_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
+    double *x, double *y, double *z, int *panSuccess )
+{
+  Q_UNUSED( z )
+  HelmertParameters *t = static_cast<HelmertParameters *>( pTransformerArg );
+  if ( !t )
+    return false;
+
+  double a = std::cos( t->angle ), b = std::sin( t->angle ), x0 = t->origin.x(), y0 = t->origin.y(), s = t->scale;
+  if ( !bDstToSrc )
+  {
+    a *= s;
+    b *= s;
+    for ( int i = 0; i < nPointCount; ++i )
+    {
+      double xT = x[i], yT = y[i];
+      // Because rotation parameters where estimated in a CS with negative y-axis ^= down.
+      // we need to apply the rotation matrix and a change of base:
+      // |cos a,-sin a| |1, 0|   | std::cos a,  std::sin a|
+      // |sin a, std::cos a| |0,-1| = | std::sin a, -cos a|
+      x[i] = x0 + ( a * xT + b * yT );
+      y[i] = y0 + ( b * xT - a * yT );
+      panSuccess[i] = true;
+    }
+  }
+  else
+  {
+    // Guard against division by zero
+    if ( std::fabs( s ) < std::numeric_limits<double>::epsilon() )
+    {
+      for ( int i = 0; i < nPointCount; ++i )
+      {
+        panSuccess[i] = false;
+      }
+      return false;
+    }
+    a /= s;
+    b /= s;
+    for ( int i = 0; i < nPointCount; ++i )
+    {
+      double xT = x[i], yT = y[i];
+      xT -= x0;
+      yT -= y0;
+      // | std::cos a,  std::sin a |^-1   |cos a,  std::sin a|
+      // | std::sin a, -cos a |    = |sin a, -cos a|
+      x[i] = a * xT + b * yT;
+      y[i] = b * xT - a * yT;
+      panSuccess[i] = true;
+    }
+  }
+  return true;
+}
+
+//
+// QgsGDALGeorefTransform
+//
+
+QgsGDALGeorefTransform::QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder )
+  : mPolynomialOrder( std::min( 3u, polynomialOrder ) )
+  , mIsTPSTransform( useTPS )
+{
+  mGDALTransformer     = nullptr;
+  mGDALTransformerArgs = nullptr;
+}
+
+QgsGDALGeorefTransform::~QgsGDALGeorefTransform()
+{
+  destroyGdalArgs();
+}
+
+bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
+{
+  assert( mapCoords.size() == pixelCoords.size() );
+  if ( mapCoords.size() != pixelCoords.size() )
+    return false;
+  int n = mapCoords.size();
+
+  GDAL_GCP *GCPList = new GDAL_GCP[n];
+  for ( int i = 0; i < n; i++ )
+  {
+    GCPList[i].pszId = new char[20];
+    snprintf( GCPList[i].pszId, 19, "gcp%i", i );
+    GCPList[i].pszInfo = nullptr;
+    GCPList[i].dfGCPPixel = pixelCoords[i].x();
+    GCPList[i].dfGCPLine  = -pixelCoords[i].y();
+    GCPList[i].dfGCPX = mapCoords[i].x();
+    GCPList[i].dfGCPY = mapCoords[i].y();
+    GCPList[i].dfGCPZ = 0;
+  }
+  destroyGdalArgs();
+
+  if ( mIsTPSTransform )
+    mGDALTransformerArgs = GDALCreateTPSTransformer( n, GCPList, false );
+  else
+    mGDALTransformerArgs = GDALCreateGCPTransformer( n, GCPList, mPolynomialOrder, false );
+
+  for ( int i = 0; i < n; i++ )
+  {
+    delete [] GCPList[i].pszId;
+  }
+  delete [] GCPList;
+
+  return nullptr != mGDALTransformerArgs;
+}
+
+int QgsGDALGeorefTransform::minimumGcpCount() const
+{
+  if ( mIsTPSTransform )
+    return 1;
+  else
+    return ( ( mPolynomialOrder + 2 ) * ( mPolynomialOrder + 1 ) ) / 2;
+}
+
+GDALTransformerFunc QgsGDALGeorefTransform::GDALTransformer() const
+{
+  // Fail if no arguments were calculated through updateParametersFromGCP
+  if ( !mGDALTransformerArgs )
+    return nullptr;
+
+  if ( mIsTPSTransform )
+    return GDALTPSTransform;
+  else
+    return GDALGCPTransform;
+}
+
+void *QgsGDALGeorefTransform::GDALTransformerArgs() const
+{
+  return mGDALTransformerArgs;
+}
+
+void QgsGDALGeorefTransform::destroyGdalArgs()
+{
+  if ( mGDALTransformerArgs )
+  {
+    if ( mIsTPSTransform )
+      GDALDestroyTPSTransformer( mGDALTransformerArgs );
+    else
+      GDALDestroyGCPTransformer( mGDALTransformerArgs );
+  }
+}
+
+//
+// QgsProjectiveGeorefTransform
+//
+
+QgsProjectiveGeorefTransform::QgsProjectiveGeorefTransform()
+  : mParameters()
+{}
+
+bool QgsProjectiveGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
+{
+  if ( mapCoords.size() < minimumGcpCount() )
+    return false;
+
+  // HACK: flip y coordinates, because georeferencer and gdal use different conventions
+  QVector<QgsPointXY> flippedPixelCoords;
+  flippedPixelCoords.reserve( pixelCoords.size() );
+  for ( const QgsPointXY &coord : pixelCoords )
+  {
+    flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() );
+  }
+
+  QgsLeastSquares::projective( mapCoords, flippedPixelCoords, mParameters.H );
+
+  // Invert the homography matrix using adjoint matrix
+  double *H = mParameters.H;
+
+  double adjoint[9];
+  adjoint[0] = H[4] * H[8] - H[5] * H[7];
+  adjoint[1] = -H[1] * H[8] + H[2] * H[7];
+  adjoint[2] = H[1] * H[5] - H[2] * H[4];
+
+  adjoint[3] = -H[3] * H[8] + H[5] * H[6];
+  adjoint[4] = H[0] * H[8] - H[2] * H[6];
+  adjoint[5] = -H[0] * H[5] + H[2] * H[3];
+
+  adjoint[6] = H[3] * H[7] - H[4] * H[6];
+  adjoint[7] = -H[0] * H[7] + H[1] * H[6];
+  adjoint[8] = H[0] * H[4] - H[1] * H[3];
+
+  double det = H[0] * adjoint[0] + H[3] * adjoint[1] + H[6] * adjoint[2];
+
+  if ( std::fabs( det ) < 1024.0 * std::numeric_limits<double>::epsilon() )
+  {
+    mParameters.hasInverse = false;
+  }
+  else
+  {
+    mParameters.hasInverse = true;
+    double oo_det = 1.0 / det;
+    for ( int i = 0; i < 9; i++ )
+    {
+      mParameters.Hinv[i] = adjoint[i] * oo_det;
+    }
+  }
+  return true;
+}
+
+int QgsProjectiveGeorefTransform::minimumGcpCount() const
+{
+  return 4;
+}
+
+GDALTransformerFunc QgsProjectiveGeorefTransform::GDALTransformer() const
+{
+  return QgsProjectiveGeorefTransform::projectiveTransform;
+}
+
+void *QgsProjectiveGeorefTransform::GDALTransformerArgs() const
+{
+  return ( void * )&mParameters;
+}
+
+int QgsProjectiveGeorefTransform::projectiveTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
+    double *x, double *y, double *z, int *panSuccess )
+{
+  Q_UNUSED( z )
+  ProjectiveParameters *t = static_cast<ProjectiveParameters *>( pTransformerArg );
+  if ( !t )
+    return false;
+
+  double *H = nullptr;
+  if ( !bDstToSrc )
+  {
+    H = t->H;
+  }
+  else
+  {
+    if ( !t->hasInverse )
+    {
+      for ( int i = 0; i < nPointCount; ++i )
+      {
+        panSuccess[i] = false;
+      }
+      return false;
+    }
+    H = t->Hinv;
+  }
+
+
+  for ( int i = 0; i < nPointCount; ++i )
+  {
+    double Z = x[i] * H[6] + y[i] * H[7] + H[8];
+    // Projects to infinity?
+    if ( std::fabs( Z ) < 1024.0 * std::numeric_limits<double>::epsilon() )
+    {
+      panSuccess[i] = false;
+      continue;
+    }
+    double X = ( x[i] * H[0] + y[i] * H[1] + H[2] ) / Z;
+    double Y = ( x[i] * H[3] + y[i] * H[4] + H[5] ) / Z;
+
+    x[i] = X;
+    y[i] = Y;
+
+    panSuccess[i] = true;
+  }
+
+  return true;
+}

src/analysis/georeferencing/qgsgcptransformer.h

@@ -0,0 +1,206 @@
+/***************************************************************************
+    qgsgcptransformer.h
+     --------------------------------------
+    Date                 : 18-Feb-2009
+    Copyright            : (c) 2009 by Manuel Massing
+    Email                : m.massing at warped-space.de
+ ***************************************************************************
+ *                                                                         *
+ *   This program is free software; you can redistribute it and/or modify  *
+ *   it under the terms of the GNU General Public License as published by  *
+ *   the Free Software Foundation; either version 2 of the License, or     *
+ *   (at your option) any later version.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef QGSGCPTRANSFORMER_H
+#define QGSGCPTRANSFORMER_H
+
+#include <gdal_alg.h>
+#include "qgspoint.h"
+#include "qgis_analysis.h"
+#include "qgis_sip.h"
+
+/**
+ * \ingroup analysis
+ * An interface for Ground Control Points (GCP) based transformations.
+ *
+ * QgsGcpTransformerInterface implementations are able to transform point locations
+ * based on a transformation method and a list of GCPs.
+ *
+ * \since QGIS 3.20
+*/
+class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT
+{
+  public:
+
+    QgsGcpTransformerInterface() = default;
+
+    virtual ~QgsGcpTransformerInterface() = default;
+
+    //! QgsGcpTransformerInterface cannot be copied
+    QgsGcpTransformerInterface( const QgsGcpTransformerInterface &other ) = delete;
+
+    //! QgsGcpTransformerInterface cannot be copied
+    QgsGcpTransformerInterface &operator=( const QgsGcpTransformerInterface &other ) = delete;
+
+    /**
+     * Fits transformation parameters using the specified Ground Control Points (GCPs) lists of map coordinates and pixel coordinates.
+     *
+     * \returns TRUE on success, FALSE on failure
+     */
+    virtual bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoordinates, const QVector<QgsPointXY> &pixelCoordinates ) = 0;
+
+    /**
+     * Returns the minimum number of Ground Control Points (GCPs) required for parameter fitting.
+     */
+    virtual int minimumGcpCount() const = 0;
+
+#ifndef SIP_RUN
+
+    /**
+     * Returns function pointer to the GDALTransformer function.
+     */
+    virtual GDALTransformerFunc GDALTransformer() const = 0;
+
+    /**
+     * Returns pointer to the GDALTransformer arguments.
+     */
+    virtual void *GDALTransformerArgs() const = 0;
+#endif
+
+  private:
+
+#ifdef SIP_RUN
+    QgsGcpTransformerInterface( const QgsGcpTransformerInterface &other )
+#endif
+};
+
+/**
+ * A simple transform which is parametrized by a translation and anistotropic scale.
+ * \ingroup analysis
+ * \note Not available in Python bindings
+ * \since QGIS 3.20
+ */
+class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP
+{
+  public:
+    QgsLinearGeorefTransform() = default;
+
+    /**
+     * Returns the origin and scale for the transform.
+     */
+    bool getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const;
+
+    bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
+    int minimumGcpCount() const override;
+    GDALTransformerFunc GDALTransformer() const override;
+    void *GDALTransformerArgs() const override;
+
+  private:
+    struct LinearParameters
+    {
+      QgsPointXY origin;
+      double scaleX, scaleY;
+    } mParameters;
+
+    static int linearTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
+                                double *x, double *y, double *z, int *panSuccess );
+
+};
+
+/**
+ * 2-dimensional helmert transform, parametrised by isotropic scale, rotation angle and translation.
+ * \ingroup analysis
+ * \note Not available in Python bindings
+ * \since QGIS 3.20
+ */
+class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP
+{
+  public:
+    QgsHelmertGeorefTransform() = default;
+
+    /**
+     * Returns the origin, scale and rotation for the transform.
+     */
+    bool getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const;
+
+    bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
+    int minimumGcpCount() const override;
+    GDALTransformerFunc GDALTransformer() const override;
+    void *GDALTransformerArgs() const override;
+
+  private:
+
+    struct HelmertParameters
+    {
+      QgsPointXY origin;
+      double scale;
+      double angle;
+    };
+    HelmertParameters mHelmertParameters;
+
+    static int helmert_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
+                                  double *x, double *y, double *z, int *panSuccess );
+
+};
+
+/**
+ * Interface to gdal thin plate splines and 1st/2nd/3rd order polynomials.
+ * \ingroup analysis
+ * \note Not available in Python bindings
+ * \since QGIS 3.20
+ */
+class ANALYSIS_EXPORT QgsGDALGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP
+{
+  public:
+    QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder );
+    ~QgsGDALGeorefTransform() override;
+
+    bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
+    int minimumGcpCount() const override;
+    GDALTransformerFunc GDALTransformer() const override;
+    void *GDALTransformerArgs() const override;
+  private:
+    void destroyGdalArgs();
+
+    const int  mPolynomialOrder;
+    const bool mIsTPSTransform;
+
+    GDALTransformerFunc  mGDALTransformer;
+    void                *mGDALTransformerArgs = nullptr;
+
+};
+
+/**
+ * A planar projective transform, expressed by a homography.
+ *
+ * Implements model fitting which minimizes algebraic error using total least squares.
+ *
+ * \ingroup analysis
+ * \note Not available in Python bindings
+ * \since QGIS 3.20
+ */
+class ANALYSIS_EXPORT QgsProjectiveGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP
+{
+  public:
+    QgsProjectiveGeorefTransform();
+
+    bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
+    int minimumGcpCount() const override;
+    GDALTransformerFunc GDALTransformer() const override;
+    void *GDALTransformerArgs() const override;
+  private:
+    struct ProjectiveParameters
+    {
+      double H[9];        // Homography
+      double Hinv[9];     // Inverted homography
+      bool hasInverse;  // Could the inverted homography be calculated?
+    } mParameters;
+
+    static int projectiveTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
+                                    double *x, double *y, double *z, int *panSuccess );
+
+};
+
+#endif //QGSGCPTRANSFORMER_H

src/analysis/georeferencing/qgsleastsquares.cpp

@@ -118,7 +118,7 @@ void QgsLeastSquares::helmert( const QVector<QgsPointXY> &mapCoords,
   rotation = std::atan2( gsl_vector_get( x, 1 ), gsl_vector_get( x, 0 ) );
 }
 
-
+#if 0
 void QgsLeastSquares::affine( QVector<QgsPointXY> mapCoords,
                               QVector<QgsPointXY> pixelCoords )
 {
@@ -171,7 +171,7 @@ void QgsLeastSquares::affine( QVector<QgsPointXY> mapCoords,
   gsl_permutation_free( p );
 
 }
-
+#endif
 
 /**
  * Scales the given coordinates so that the center of gravity is at the origin and the mean distance to the origin is sqrt(2).

src/analysis/georeferencing/qgsleastsquares.h

@@ -16,30 +16,50 @@
 #define QGSLEASTSQUARES_H
 
 #include <QVector>
-#include <cstdarg>
-#include <stdexcept>
 
+#include "qgis_analysis.h"
 #include "qgspointxy.h"
 
+#define SIP_NO_FILE
 
-class QgsLeastSquares
+/**
+ * \ingroup analysis
+ * Utilities for calculation of least squares based transformations.
+ *
+ * \note Not available in Python bindings.
+ * \since QGIS 3.20
+*/
+class ANALYSIS_EXPORT QgsLeastSquares
 {
   public:
+
+    /**
+     * Transforms the point at \a origin in-place, using a linear transformation calculated from the list of map and pixel Ground Control Points (GCPs).
+     */
     static void linear( const QVector<QgsPointXY> &mapCoords,
                         const QVector<QgsPointXY> &pixelCoords,
                         QgsPointXY &origin, double &pixelXSize, double &pixelYSize );
 
+    /**
+     * Transforms the point at \a origin in-place, using a helmert transformation calculated from the list of map and pixel Ground Control Points (GCPs).
+     */
     static void helmert( const QVector<QgsPointXY> &mapCoords,
                          const QVector<QgsPointXY> &pixelCoords,
                          QgsPointXY &origin, double &pixelSize, double &rotation );
 
+#if 0
     static void affine( QVector<QgsPointXY> mapCoords,
                         QVector<QgsPointXY> pixelCoords );
 
+#endif
+
+    /**
+     * Calculates projective parameters from the list of map and pixel Ground Control Points (GCPs).
+     */
     static void projective( QVector<QgsPointXY> mapCoords,
                             QVector<QgsPointXY> pixelCoords,
                             double H[9] );
+
 };
 
-
-#endif
+#endif // QGSLEASTSQUARES_H

src/app/CMakeLists.txt

@@ -243,7 +243,6 @@ if (WITH_GEOREFERENCER)
     georeferencer/qgsgeoreftooldeletepoint.cpp
     georeferencer/qgsgeoreftoolmovepoint.cpp
     georeferencer/qgsgeorefvalidators.cpp
-    georeferencer/qgsleastsquares.cpp
     georeferencer/qgsmapcoordsdialog.cpp
     georeferencer/qgsresidualplotitem.cpp
     georeferencer/qgstransformsettingsdialog.cpp

src/app/georeferencer/qgsgcplistmodel.cpp

@@ -89,7 +89,7 @@ void QgsGCPListModel::updateModel()
 
   if ( mGeorefTransform )
   {
-    bTransformUpdated = mGeorefTransform->updateParametersFromGCPs( mapCoords, pixelCoords );
+    bTransformUpdated = mGeorefTransform->updateParametersFromGcps( mapCoords, pixelCoords );
     mapUnitsPossible = mGeorefTransform->providesAccurateInverseTransformation();
   }
 

src/app/georeferencer/qgsgeorefmainwindow.cpp

@@ -59,7 +59,6 @@
 #include "qgsgeoreftoolmovepoint.h"
 #include "qgsgcpcanvasitem.h"
 
-#include "qgsleastsquares.h"
 #include "qgsgcplistwidget.h"
 
 #include "qgsgeorefconfigdialog.h"
@@ -1496,12 +1495,12 @@ bool QgsGeoreferencerMainWindow::calculateMeanError( double &error ) const
     }
   }
 
-  if ( nPointsEnabled == mGeorefTransform.getMinimumGCPCount() )
+  if ( nPointsEnabled == mGeorefTransform.minimumGcpCount() )
   {
     error = 0;
     return true;
   }
-  else if ( nPointsEnabled < mGeorefTransform.getMinimumGCPCount() )
+  else if ( nPointsEnabled < mGeorefTransform.minimumGcpCount() )
   {
     return false;
   }
@@ -1521,7 +1520,7 @@ bool QgsGeoreferencerMainWindow::calculateMeanError( double &error ) const
 
   // Calculate the root mean square error, adjusted for degrees of freedom of the transform
   // Caveat: The number of DoFs is assumed to be even (as each control point fixes two degrees of freedom).
-  error = std::sqrt( ( sumVxSquare + sumVySquare ) / ( nPointsEnabled - mGeorefTransform.getMinimumGCPCount() ) );
+  error = std::sqrt( ( sumVxSquare + sumVySquare ) / ( nPointsEnabled - mGeorefTransform.minimumGcpCount() ) );
   return true;
 }
 
@@ -1977,10 +1976,10 @@ bool QgsGeoreferencerMainWindow::checkReadyGeoref()
     return false;
   }
 
-  if ( mPoints.count() < static_cast<int>( mGeorefTransform.getMinimumGCPCount() ) )
+  if ( mPoints.count() < static_cast<int>( mGeorefTransform.minimumGcpCount() ) )
   {
     mMessageBar->pushMessage( tr( "Not Enough GCPs" ), tr( "%1 transformation requires at least %2 GCPs. Please define more." )
-                              .arg( convertTransformEnumToString( mTransformParam ) ).arg( mGeorefTransform.getMinimumGCPCount() )
+                              .arg( convertTransformEnumToString( mTransformParam ) ).arg( mGeorefTransform.minimumGcpCount() )
                               , Qgis::Critical );
     return false;
   }
@@ -2005,7 +2004,7 @@ bool QgsGeoreferencerMainWindow::updateGeorefTransform()
     return false;
 
   // Parametrize the transform with GCPs
-  if ( !mGeorefTransform.updateParametersFromGCPs( mapCoords, pixelCoords ) )
+  if ( !mGeorefTransform.updateParametersFromGcps( mapCoords, pixelCoords ) )
   {
     return false;
   }
@@ -2180,14 +2179,14 @@ bool QgsGeoreferencerMainWindow::equalGCPlists( const QgsGCPList &list1, const Q
 //
 //void QgsGeorefPluginGui::logRequaredGCPs()
 //{
-//  if (mGeorefTransform.getMinimumGCPCount() != 0)
+//  if (mGeorefTransform.minimumGcpCount() != 0)
 //  {
-//    if ((uint)mPoints.size() >= mGeorefTransform.getMinimumGCPCount())
+//    if ((uint)mPoints.size() >= mGeorefTransform.minimumGcpCount())
 //      showMessageInLog(tr("Info"), tr("For georeferencing requared at least %1 GCP points")
-//                       .arg(mGeorefTransform.getMinimumGCPCount()));
+//                       .arg(mGeorefTransform.minimumGcpCount()));
 //    else
 //      showMessageInLog(tr("Critical"), tr("For georeferencing requared at least %1 GCP points")
-//                       .arg(mGeorefTransform.getMinimumGCPCount()));
+//                       .arg(mGeorefTransform.minimumGcpCount()));
 //  }
 //}
 

src/app/georeferencer/qgsgeoreftransform.cpp

@@ -19,118 +19,11 @@
 #include <gdal.h>
 #include <gdal_alg.h>
 
-#include "qgsleastsquares.h"
-
 #include <cmath>
 
 #include <cassert>
 #include <limits>
 
-/**
- * A simple transform which is paremetrized by a translation and anistotropic scale.
- */
-class QgsLinearGeorefTransform : public QgsGeorefTransformInterface
-{
-  public:
-    QgsLinearGeorefTransform() = default;
-
-    bool getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const;
-
-    bool updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
-    int getMinimumGCPCount() const override;
-
-    GDALTransformerFunc  GDALTransformer()     const override { return QgsLinearGeorefTransform::linear_transform; }
-    void                *GDALTransformerArgs() const override { return ( void * )&mParameters; }
-  private:
-    struct LinearParameters
-    {
-      QgsPointXY origin;
-      double scaleX, scaleY;
-    } mParameters;
-
-    static int linear_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
-                                 double *x, double *y, double *z, int *panSuccess );
-};
-
-/**
- * 2-dimensional helmert transform, parametrised by isotropic scale, rotation angle and translation.
- */
-class QgsHelmertGeorefTransform : public QgsGeorefTransformInterface
-{
-  public:
-    QgsHelmertGeorefTransform() = default;
-    struct HelmertParameters
-    {
-      QgsPointXY origin;
-      double   scale;
-      double   angle;
-    };
-
-    bool getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const;
-    bool updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
-    int getMinimumGCPCount() const override;
-
-    GDALTransformerFunc  GDALTransformer()     const override;
-    void                *GDALTransformerArgs() const override;
-  private:
-    HelmertParameters mHelmertParameters;
-
-    static int helmert_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
-                                  double *x, double *y, double *z, int *panSuccess );
-};
-
-/**
- * Interface to gdal thin plate splines and 1st/2nd/3rd order polynomials.
- */
-class QgsGDALGeorefTransform : public QgsGeorefTransformInterface
-{
-  public:
-    QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder );
-    ~QgsGDALGeorefTransform() override;
-
-    bool updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
-    int getMinimumGCPCount() const override;
-
-    GDALTransformerFunc  GDALTransformer()     const override;
-    void                *GDALTransformerArgs() const override;
-  private:
-    void destroy_gdal_args();
-
-    const int  mPolynomialOrder;
-    const bool mIsTPSTransform;
-
-    GDALTransformerFunc  mGDALTransformer;
-    void                *mGDALTransformerArgs = nullptr;
-};
-
-/**
- * A planar projective transform, expressed by a homography.
- *
- * Implements model fitting which minimizes algebraic error using total least squares.
- */
-class QgsProjectiveGeorefTransform : public QgsGeorefTransformInterface
-{
-  public:
-    QgsProjectiveGeorefTransform() : mParameters() {}
-
-    bool updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
-    int getMinimumGCPCount() const override;
-
-    GDALTransformerFunc  GDALTransformer()     const override { return QgsProjectiveGeorefTransform::projective_transform; }
-    void                *GDALTransformerArgs() const override { return ( void * )&mParameters; }
-  private:
-    struct ProjectiveParameters
-    {
-      double H[9];        // Homography
-      double Hinv[9];     // Inverted homography
-      bool   hasInverse;  // Could the inverted homography be calculated?
-    } mParameters;
-
-    static int projective_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
-                                     double *x, double *y, double *z, int *panSuccess );
-};
-
-
 QgsGeorefTransform::QgsGeorefTransform( const QgsGeorefTransform &other )
 {
   mTransformParametrisation = InvalidTransform;
@@ -190,7 +83,7 @@ bool QgsGeorefTransform::parametersInitialized() const
   return mParametersInitialized;
 }
 
-bool QgsGeorefTransform::updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
+bool QgsGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
 {
   if ( !mGeorefTransformImplementation )
   {
@@ -200,26 +93,26 @@ bool QgsGeorefTransform::updateParametersFromGCPs( const QVector<QgsPointXY> &ma
   {
     throw ( std::domain_error( "Internal error: GCP mapping is not one-to-one" ) );
   }
-  if ( mapCoords.size() < getMinimumGCPCount() )
+  if ( mapCoords.size() < minimumGcpCount() )
   {
     return false;
   }
   if ( mRasterChangeCoords.hasCrs() )
   {
     QVector<QgsPointXY> pixelCoordsCorrect = mRasterChangeCoords.getPixelCoords( pixelCoords );
-    mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGCPs( mapCoords, pixelCoordsCorrect );
+    mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( mapCoords, pixelCoordsCorrect );
     pixelCoordsCorrect.clear();
   }
   else
   {
-    mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGCPs( mapCoords, pixelCoords );
+    mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( mapCoords, pixelCoords );
   }
   return mParametersInitialized;
 }
 
-int QgsGeorefTransform::getMinimumGCPCount() const
+int QgsGeorefTransform::minimumGcpCount() const
 {
-  return mGeorefTransformImplementation ? mGeorefTransformImplementation->getMinimumGCPCount() : 0;
+  return mGeorefTransformImplementation ? mGeorefTransformImplementation->minimumGcpCount() : 0;
 }
 
 GDALTransformerFunc QgsGeorefTransform::GDALTransformer() const
@@ -232,7 +125,7 @@ void *QgsGeorefTransform::GDALTransformerArgs() const
   return mGeorefTransformImplementation ? mGeorefTransformImplementation->GDALTransformerArgs() : nullptr;
 }
 
-QgsGeorefTransformInterface *QgsGeorefTransform::createImplementation( TransformParametrisation parametrisation )
+QgsGcpTransformerInterface *QgsGeorefTransform::createImplementation( TransformParametrisation parametrisation )
 {
   switch ( parametrisation )
   {
@@ -338,335 +231,3 @@ bool QgsGeorefTransform::gdal_transform( const QgsPointXY &src, QgsPointXY &dst,
 }
 
 
-bool QgsLinearGeorefTransform::getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const
-{
-  origin = mParameters.origin;
-  scaleX = mParameters.scaleX;
-  scaleY = mParameters.scaleY;
-  return true;
-}
-
-bool QgsLinearGeorefTransform::updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
-{
-  if ( mapCoords.size() < getMinimumGCPCount() )
-    return false;
-  QgsLeastSquares::linear( mapCoords, pixelCoords, mParameters.origin, mParameters.scaleX, mParameters.scaleY );
-  return true;
-}
-
-int QgsLinearGeorefTransform::getMinimumGCPCount() const
-{
-  return 2;
-}
-
-int QgsLinearGeorefTransform::linear_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
-    double *x, double *y, double *z, int *panSuccess )
-{
-  Q_UNUSED( z )
-  LinearParameters *t = static_cast<LinearParameters *>( pTransformerArg );
-  if ( !t )
-    return false;
-
-  if ( !bDstToSrc )
-  {
-    for ( int i = 0; i < nPointCount; ++i )
-    {
-      x[i] = x[i] * t->scaleX + t->origin.x();
-      y[i] = -y[i] * t->scaleY + t->origin.y();
-      panSuccess[i] = true;
-    }
-  }
-  else
-  {
-    // Guard against division by zero
-    if ( std::fabs( t->scaleX ) < std::numeric_limits<double>::epsilon() ||
-         std::fabs( t->scaleY ) < std::numeric_limits<double>::epsilon() )
-    {
-      for ( int i = 0; i < nPointCount; ++i )
-      {
-        panSuccess[i] = false;
-      }
-      return false;
-    }
-    for ( int i = 0; i < nPointCount; ++i )
-    {
-      x[i] = ( x[i] - t->origin.x() ) / t->scaleX;
-      y[i] = ( y[i] - t->origin.y() ) / ( -t->scaleY );
-      panSuccess[i] = true;
-    }
-  }
-
-  return true;
-}
-
-bool QgsHelmertGeorefTransform::updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
-{
-  if ( mapCoords.size() < getMinimumGCPCount() )
-    return false;
-
-  QgsLeastSquares::helmert( mapCoords, pixelCoords, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle );
-  return true;
-}
-
-int QgsHelmertGeorefTransform::getMinimumGCPCount() const
-{
-  return 2;
-}
-
-
-GDALTransformerFunc QgsHelmertGeorefTransform::GDALTransformer() const
-{
-  return QgsHelmertGeorefTransform::helmert_transform;
-}
-
-void *QgsHelmertGeorefTransform::GDALTransformerArgs() const
-{
-  return ( void * )&mHelmertParameters;
-}
-
-bool QgsHelmertGeorefTransform::getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const
-{
-  origin = mHelmertParameters.origin;
-  scale = mHelmertParameters.scale;
-  rotation = mHelmertParameters.angle;
-  return true;
-}
-
-int QgsHelmertGeorefTransform::helmert_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
-    double *x, double *y, double *z, int *panSuccess )
-{
-  Q_UNUSED( z )
-  HelmertParameters *t = static_cast<HelmertParameters *>( pTransformerArg );
-  if ( !t )
-    return false;
-
-  double a = std::cos( t->angle ), b = std::sin( t->angle ), x0 = t->origin.x(), y0 = t->origin.y(), s = t->scale;
-  if ( !bDstToSrc )
-  {
-    a *= s;
-    b *= s;
-    for ( int i = 0; i < nPointCount; ++i )
-    {
-      double xT = x[i], yT = y[i];
-      // Because rotation parameters where estimated in a CS with negative y-axis ^= down.
-      // we need to apply the rotation matrix and a change of base:
-      // |cos a,-sin a| |1, 0|   | std::cos a,  std::sin a|
-      // |sin a, std::cos a| |0,-1| = | std::sin a, -cos a|
-      x[i] = x0 + ( a * xT + b * yT );
-      y[i] = y0 + ( b * xT - a * yT );
-      panSuccess[i] = true;
-    }
-  }
-  else
-  {
-    // Guard against division by zero
-    if ( std::fabs( s ) < std::numeric_limits<double>::epsilon() )
-    {
-      for ( int i = 0; i < nPointCount; ++i )
-      {
-        panSuccess[i] = false;
-      }
-      return false;
-    }
-    a /= s;
-    b /= s;
-    for ( int i = 0; i < nPointCount; ++i )
-    {
-      double xT = x[i], yT = y[i];
-      xT -= x0;
-      yT -= y0;
-      // | std::cos a,  std::sin a |^-1   |cos a,  std::sin a|
-      // | std::sin a, -cos a |    = |sin a, -cos a|
-      x[i] = a * xT + b * yT;
-      y[i] = b * xT - a * yT;
-      panSuccess[i] = true;
-    }
-  }
-  return true;
-}
-
-QgsGDALGeorefTransform::QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder )
-  : mPolynomialOrder( std::min( 3u, polynomialOrder ) )
-  , mIsTPSTransform( useTPS )
-{
-  mGDALTransformer     = nullptr;
-  mGDALTransformerArgs = nullptr;
-}
-
-QgsGDALGeorefTransform::~QgsGDALGeorefTransform()
-{
-  destroy_gdal_args();
-}
-
-bool QgsGDALGeorefTransform::updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
-{
-  assert( mapCoords.size() == pixelCoords.size() );
-  if ( mapCoords.size() != pixelCoords.size() )
-    return false;
-  int n = mapCoords.size();
-
-  GDAL_GCP *GCPList = new GDAL_GCP[n];
-  for ( int i = 0; i < n; i++ )
-  {
-    GCPList[i].pszId = new char[20];
-    snprintf( GCPList[i].pszId, 19, "gcp%i", i );
-    GCPList[i].pszInfo = nullptr;
-    GCPList[i].dfGCPPixel = pixelCoords[i].x();
-    GCPList[i].dfGCPLine  = -pixelCoords[i].y();
-    GCPList[i].dfGCPX = mapCoords[i].x();
-    GCPList[i].dfGCPY = mapCoords[i].y();
-    GCPList[i].dfGCPZ = 0;
-  }
-  destroy_gdal_args();
-
-  if ( mIsTPSTransform )
-    mGDALTransformerArgs = GDALCreateTPSTransformer( n, GCPList, false );
-  else
-    mGDALTransformerArgs = GDALCreateGCPTransformer( n, GCPList, mPolynomialOrder, false );
-
-  for ( int i = 0; i < n; i++ )
-  {
-    delete [] GCPList[i].pszId;
-  }
-  delete [] GCPList;
-
-  return nullptr != mGDALTransformerArgs;
-}
-
-int QgsGDALGeorefTransform::getMinimumGCPCount() const
-{
-  if ( mIsTPSTransform )
-    return 1;
-  else
-    return ( ( mPolynomialOrder + 2 ) * ( mPolynomialOrder + 1 ) ) / 2;
-}
-
-GDALTransformerFunc QgsGDALGeorefTransform::GDALTransformer() const
-{
-  // Fail if no arguments were calculated through updateParametersFromGCP
-  if ( !mGDALTransformerArgs )
-    return nullptr;
-
-  if ( mIsTPSTransform )
-    return GDALTPSTransform;
-  else
-    return GDALGCPTransform;
-}
-
-void *QgsGDALGeorefTransform::GDALTransformerArgs() const
-{
-  return mGDALTransformerArgs;
-}
-
-void QgsGDALGeorefTransform::destroy_gdal_args()
-{
-  if ( mGDALTransformerArgs )
-  {
-    if ( mIsTPSTransform )
-      GDALDestroyTPSTransformer( mGDALTransformerArgs );
-    else
-      GDALDestroyGCPTransformer( mGDALTransformerArgs );
-  }
-}
-
-bool QgsProjectiveGeorefTransform::updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
-{
-  if ( mapCoords.size() < getMinimumGCPCount() )
-    return false;
-
-  // HACK: flip y coordinates, because georeferencer and gdal use different conventions
-  QVector<QgsPointXY> flippedPixelCoords;
-  flippedPixelCoords.reserve( pixelCoords.size() );
-  for ( const QgsPointXY &coord : pixelCoords )
-  {
-    flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() );
-  }
-
-  QgsLeastSquares::projective( mapCoords, flippedPixelCoords, mParameters.H );
-
-  // Invert the homography matrix using adjoint matrix
-  double *H = mParameters.H;
-
-  double adjoint[9];
-  adjoint[0] = H[4] * H[8] - H[5] * H[7];
-  adjoint[1] = -H[1] * H[8] + H[2] * H[7];
-  adjoint[2] = H[1] * H[5] - H[2] * H[4];
-
-  adjoint[3] = -H[3] * H[8] + H[5] * H[6];
-  adjoint[4] = H[0] * H[8] - H[2] * H[6];
-  adjoint[5] = -H[0] * H[5] + H[2] * H[3];
-
-  adjoint[6] = H[3] * H[7] - H[4] * H[6];
-  adjoint[7] = -H[0] * H[7] + H[1] * H[6];
-  adjoint[8] = H[0] * H[4] - H[1] * H[3];
-
-  double det = H[0] * adjoint[0] + H[3] * adjoint[1] + H[6] * adjoint[2];
-
-  if ( std::fabs( det ) < 1024.0 * std::numeric_limits<double>::epsilon() )
-  {
-    mParameters.hasInverse = false;
-  }
-  else
-  {
-    mParameters.hasInverse = true;
-    double oo_det = 1.0 / det;
-    for ( int i = 0; i < 9; i++ )
-    {
-      mParameters.Hinv[i] = adjoint[i] * oo_det;
-    }
-  }
-  return true;
-}
-
-int QgsProjectiveGeorefTransform::getMinimumGCPCount() const
-{
-  return 4;
-}
-
-int QgsProjectiveGeorefTransform::projective_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
-    double *x, double *y, double *z, int *panSuccess )
-{
-  Q_UNUSED( z )
-  ProjectiveParameters *t = static_cast<ProjectiveParameters *>( pTransformerArg );
-  if ( !t )
-    return false;
-
-  double *H = nullptr;
-  if ( !bDstToSrc )
-  {
-    H = t->H;
-  }
-  else
-  {
-    if ( !t->hasInverse )
-    {
-      for ( int i = 0; i < nPointCount; ++i )
-      {
-        panSuccess[i] = false;
-      }
-      return false;
-    }
-    H = t->Hinv;
-  }
-
-
-  for ( int i = 0; i < nPointCount; ++i )
-  {
-    double Z = x[i] * H[6] + y[i] * H[7] + H[8];
-    // Projects to infinity?
-    if ( std::fabs( Z ) < 1024.0 * std::numeric_limits<double>::epsilon() )
-    {
-      panSuccess[i] = false;
-      continue;
-    }
-    double X = ( x[i] * H[0] + y[i] * H[1] + H[2] ) / Z;
-    double Y = ( x[i] * H[3] + y[i] * H[4] + H[5] ) / Z;
-
-    x[i] = X;
-    y[i] = Y;
-
-    panSuccess[i] = true;
-  }
-
-  return true;
-}

src/app/georeferencer/qgsgeoreftransform.h

@@ -14,37 +14,14 @@
  *                                                                         *
  ***************************************************************************/
 
-#ifndef QGS_GEOREF_TRANSFORM_H
-#define QGS_GEOREF_TRANSFORM_H
-
-#include <gdal_alg.h> // just needed for GDALTransformerFunc, forward?
+#ifndef QGSGEOREFTRANSFORM_H
+#define QGSGEOREFTRANSFORM_H
 
+#include <gdal_alg.h>
 #include "qgspoint.h"
-#include <QVector>
-#include <stdexcept>
-
+#include "qgsgcptransformer.h"
 #include "qgsrasterchangecoords.h"
 
-class QgsGeorefTransformInterface
-{
-  public:
-    virtual ~QgsGeorefTransformInterface() = default;
-
-    virtual bool updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) = 0;
-
-    /**
-     * Returns the minimum number of GCPs required for parameter fitting.
-     */
-    virtual int getMinimumGCPCount() const = 0;
-
-    /**
-     * Returns function pointer to the GDALTransformer function.
-     * Used by GDALwarp.
-     */
-    virtual GDALTransformerFunc  GDALTransformer()     const = 0;
-    virtual void                *GDALTransformerArgs() const = 0;
-};
-
 /**
  * \brief Transform class for different gcp-based transform methods.
  *
@@ -56,7 +33,7 @@ class QgsGeorefTransformInterface
  * Delegates to concrete implementations of \ref QgsGeorefInterface. For exception safety,
  * this is preferred over using the subclasses directly.
  */
-class QgsGeorefTransform : public QgsGeorefTransformInterface
+class QgsGeorefTransform : public QgsGcpTransformerInterface
 {
   public:
     // GCP based transform methods implemented by subclasses
@@ -104,15 +81,8 @@ class QgsGeorefTransform : public QgsGeorefTransformInterface
     //! \returns whether the parameters of this transform have been initialized by \ref updateParametersFromGCPs
     bool parametersInitialized() const;
 
-    /**
-     * \brief Fits transformation parameters to the supplied ground control points.
-     *
-     * \returns TRUE on success, FALSE on failure
-     */
-    bool updateParametersFromGCPs( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
-
-    //! \brief Returns the minimum number of GCPs required for parameter fitting.
-    int getMinimumGCPCount() const override;
+    bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
+    int minimumGcpCount() const override;
 
     /**
      * Returns function pointer to the GDALTransformer function.
@@ -158,15 +128,15 @@ class QgsGeorefTransform : public QgsGeorefTransformInterface
     QgsGeorefTransform &operator= ( const QgsGeorefTransform & ) = delete;
 
     //! Factory function which creates an implementation for the given parametrisation.
-    static QgsGeorefTransformInterface *createImplementation( TransformParametrisation parametrisation );
+    static QgsGcpTransformerInterface *createImplementation( TransformParametrisation parametrisation );
 
     // convenience wrapper around GDALTransformerFunc
     bool gdal_transform( const QgsPointXY &src, QgsPointXY &dst, int dstToSrc ) const;
 
-    QgsGeorefTransformInterface *mGeorefTransformImplementation = nullptr;
+    QgsGcpTransformerInterface *mGeorefTransformImplementation = nullptr;
     TransformParametrisation     mTransformParametrisation;
     bool                         mParametersInitialized;
     QgsRasterChangeCoords        mRasterChangeCoords;
 };
 
-#endif
+#endif //QGSGEOREFTRANSFORM_H

src/app/georeferencer/qgstransformsettingsdialog.cpp

@@ -266,7 +266,7 @@ bool QgsTransformSettingsDialog::checkGCPpoints( int count, int &minGCPpoints )
 {
   QgsGeorefTransform georefTransform;
   georefTransform.selectTransformParametrisation( ( QgsGeorefTransform::TransformParametrisation )count );
-  minGCPpoints = georefTransform.getMinimumGCPCount();
+  minGCPpoints = georefTransform.minimumGcpCount();
   return ( mCountGCPpoints >= minGCPpoints );
 }