diff --git a/doc/CMakeLists.txt b/doc/CMakeLists.txt index 32f432ce852..2a4289ac431 100644 --- a/doc/CMakeLists.txt +++ b/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 diff --git a/python/CMakeLists.txt b/python/CMakeLists.txt index 1e33960b081..db7953a56c4 100644 --- a/python/CMakeLists.txt +++ b/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 diff --git a/python/analysis/analysis_auto.sip b/python/analysis/analysis_auto.sip index 8a72e867a46..602360630b1 100644 --- a/python/analysis/analysis_auto.sip +++ b/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 diff --git a/python/analysis/auto_generated/georeferencing/qgsgcptransformer.sip.in b/python/analysis/auto_generated/georeferencing/qgsgcptransformer.sip.in new file mode 100644 index 00000000000..18f249be004 --- /dev/null +++ b/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 &mapCoordinates, const QVector &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 * + ************************************************************************/ diff --git a/src/analysis/CMakeLists.txt b/src/analysis/CMakeLists.txt index cb9e14153bf..3358763d818 100644 --- a/src/analysis/CMakeLists.txt +++ b/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) diff --git a/src/analysis/georeferencing/qgsgcptransformer.cpp b/src/analysis/georeferencing/qgsgcptransformer.cpp new file mode 100644 index 00000000000..55d132bcf7e --- /dev/null +++ b/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 +#include + +#include "qgsleastsquares.h" + +#include + +#include +#include + +// +// 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 &mapCoords, const QVector &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( 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::epsilon() || + std::fabs( t->scaleY ) < std::numeric_limits::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 &mapCoords, const QVector &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( 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::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 &mapCoords, const QVector &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 &mapCoords, const QVector &pixelCoords ) +{ + if ( mapCoords.size() < minimumGcpCount() ) + return false; + + // HACK: flip y coordinates, because georeferencer and gdal use different conventions + QVector 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::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( 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::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; +} diff --git a/src/analysis/georeferencing/qgsgcptransformer.h b/src/analysis/georeferencing/qgsgcptransformer.h new file mode 100644 index 00000000000..3a7fc94fdcd --- /dev/null +++ b/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 +#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 &mapCoordinates, const QVector &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 &mapCoords, const QVector &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 &mapCoords, const QVector &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 &mapCoords, const QVector &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 &mapCoords, const QVector &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 diff --git a/src/app/georeferencer/qgsleastsquares.cpp b/src/analysis/georeferencing/qgsleastsquares.cpp similarity index 99% rename from src/app/georeferencer/qgsleastsquares.cpp rename to src/analysis/georeferencing/qgsleastsquares.cpp index cbf20c7bdcb..77f6a4252d7 100644 --- a/src/app/georeferencer/qgsleastsquares.cpp +++ b/src/analysis/georeferencing/qgsleastsquares.cpp @@ -118,7 +118,7 @@ void QgsLeastSquares::helmert( const QVector &mapCoords, rotation = std::atan2( gsl_vector_get( x, 1 ), gsl_vector_get( x, 0 ) ); } - +#if 0 void QgsLeastSquares::affine( QVector mapCoords, QVector pixelCoords ) { @@ -171,7 +171,7 @@ void QgsLeastSquares::affine( QVector 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). diff --git a/src/app/georeferencer/qgsleastsquares.h b/src/analysis/georeferencing/qgsleastsquares.h similarity index 70% rename from src/app/georeferencer/qgsleastsquares.h rename to src/analysis/georeferencing/qgsleastsquares.h index f78d4166e06..3d87cd949a5 100644 --- a/src/app/georeferencer/qgsleastsquares.h +++ b/src/analysis/georeferencing/qgsleastsquares.h @@ -16,30 +16,50 @@ #define QGSLEASTSQUARES_H #include -#include -#include +#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 &mapCoords, const QVector &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 &mapCoords, const QVector &pixelCoords, QgsPointXY &origin, double &pixelSize, double &rotation ); +#if 0 static void affine( QVector mapCoords, QVector pixelCoords ); +#endif + + /** + * Calculates projective parameters from the list of map and pixel Ground Control Points (GCPs). + */ static void projective( QVector mapCoords, QVector pixelCoords, double H[9] ); + }; - -#endif +#endif // QGSLEASTSQUARES_H diff --git a/src/app/CMakeLists.txt b/src/app/CMakeLists.txt index 072e40a9bbb..60cedacd1e4 100644 --- a/src/app/CMakeLists.txt +++ b/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 diff --git a/src/app/georeferencer/qgsgcplistmodel.cpp b/src/app/georeferencer/qgsgcplistmodel.cpp index 59524587adc..016de2b3f25 100644 --- a/src/app/georeferencer/qgsgcplistmodel.cpp +++ b/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(); } diff --git a/src/app/georeferencer/qgsgeorefmainwindow.cpp b/src/app/georeferencer/qgsgeorefmainwindow.cpp index d3a27cca9cf..2236018ff54 100644 --- a/src/app/georeferencer/qgsgeorefmainwindow.cpp +++ b/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( mGeorefTransform.getMinimumGCPCount() ) ) + if ( mPoints.count() < static_cast( 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())); // } //} diff --git a/src/app/georeferencer/qgsgeoreftransform.cpp b/src/app/georeferencer/qgsgeoreftransform.cpp index 56ad4fcca1a..53b200fb4b1 100644 --- a/src/app/georeferencer/qgsgeoreftransform.cpp +++ b/src/app/georeferencer/qgsgeoreftransform.cpp @@ -19,118 +19,11 @@ #include #include -#include "qgsleastsquares.h" - #include #include #include -/** - * 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 &mapCoords, const QVector &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 &mapCoords, const QVector &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 &mapCoords, const QVector &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 &mapCoords, const QVector &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 &mapCoords, const QVector &pixelCoords ) +bool QgsGeorefTransform::updateParametersFromGcps( const QVector &mapCoords, const QVector &pixelCoords ) { if ( !mGeorefTransformImplementation ) { @@ -200,26 +93,26 @@ bool QgsGeorefTransform::updateParametersFromGCPs( const QVector &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 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 &mapCoords, const QVector &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( 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::epsilon() || - std::fabs( t->scaleY ) < std::numeric_limits::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 &mapCoords, const QVector &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( 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::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 &mapCoords, const QVector &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 &mapCoords, const QVector &pixelCoords ) -{ - if ( mapCoords.size() < getMinimumGCPCount() ) - return false; - - // HACK: flip y coordinates, because georeferencer and gdal use different conventions - QVector 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::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( 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::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; -} diff --git a/src/app/georeferencer/qgsgeoreftransform.h b/src/app/georeferencer/qgsgeoreftransform.h index c9f4df9a53d..80956a76aaa 100644 --- a/src/app/georeferencer/qgsgeoreftransform.h +++ b/src/app/georeferencer/qgsgeoreftransform.h @@ -14,37 +14,14 @@ * * ***************************************************************************/ -#ifndef QGS_GEOREF_TRANSFORM_H -#define QGS_GEOREF_TRANSFORM_H - -#include // just needed for GDALTransformerFunc, forward? +#ifndef QGSGEOREFTRANSFORM_H +#define QGSGEOREFTRANSFORM_H +#include #include "qgspoint.h" -#include -#include - +#include "qgsgcptransformer.h" #include "qgsrasterchangecoords.h" -class QgsGeorefTransformInterface -{ - public: - virtual ~QgsGeorefTransformInterface() = default; - - virtual bool updateParametersFromGCPs( const QVector &mapCoords, const QVector &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 &mapCoords, const QVector &pixelCoords ) override; - - //! \brief Returns the minimum number of GCPs required for parameter fitting. - int getMinimumGCPCount() const override; + bool updateParametersFromGcps( const QVector &mapCoords, const QVector &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 diff --git a/src/app/georeferencer/qgstransformsettingsdialog.cpp b/src/app/georeferencer/qgstransformsettingsdialog.cpp index 613f81b8419..8c456be458d 100644 --- a/src/app/georeferencer/qgstransformsettingsdialog.cpp +++ b/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 ); }