diff --git a/python/analysis/auto_generated/georeferencing/qgsgcptransformer.sip.in b/python/analysis/auto_generated/georeferencing/qgsgcptransformer.sip.in index 88dea23bcd9..0ec5ec3ddb1 100644 --- a/python/analysis/auto_generated/georeferencing/qgsgcptransformer.sip.in +++ b/python/analysis/auto_generated/georeferencing/qgsgcptransformer.sip.in @@ -37,6 +37,9 @@ based on a transformation method and a list of GCPs. }; QgsGcpTransformerInterface(); +%Docstring +Constructor for QgsGcpTransformerInterface +%End virtual ~QgsGcpTransformerInterface(); @@ -50,9 +53,12 @@ parameters as this one. Caller takes ownership of the returned object. %End - virtual bool updateParametersFromGcps( const QVector &mapCoordinates, const QVector &layerCoordinates ) = 0; + virtual bool updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis = false ) = 0; %Docstring -Fits transformation parameters using the specified Ground Control Points (GCPs) lists of map coordinates and layer coordinates. +Fits transformation parameters using the specified Ground Control Points (GCPs) lists of source and destination coordinates. + +If ``invertYAxis`` is set to ``True`` then the y-axis of source coordinates will be inverted, e.g. to allow for transformation of raster layers +with ascending top-to-bottom vertical axis coordinates. :return: ``True`` on success, ``False`` on failure %End @@ -88,10 +94,10 @@ Creates a new QgsGcpTransformerInterface subclass representing the specified tra Caller takes ownership of the returned object. %End - static QgsGcpTransformerInterface *createFromParameters( TransformMethod method, const QVector &mapCoordinates, const QVector &layerCoordinates ) /Factory/; + static QgsGcpTransformerInterface *createFromParameters( TransformMethod method, const QVector &sourceCoordinates, const QVector &destinationCoordinates ) /Factory/; %Docstring Creates a new QgsGcpTransformerInterface subclass representing the specified transform ``method``, initialized -using the given lists of map and layer coordinates. +using the given lists of source and destination coordinates. If the parameters cannot be fit to a transform ``None`` will be returned. diff --git a/src/analysis/georeferencing/qgsgcptransformer.cpp b/src/analysis/georeferencing/qgsgcptransformer.cpp index 5e9200701ea..8b3bb04492b 100644 --- a/src/analysis/georeferencing/qgsgcptransformer.cpp +++ b/src/analysis/georeferencing/qgsgcptransformer.cpp @@ -90,13 +90,13 @@ QgsGcpTransformerInterface *QgsGcpTransformerInterface::create( QgsGcpTransforme } } -QgsGcpTransformerInterface *QgsGcpTransformerInterface::createFromParameters( QgsGcpTransformerInterface::TransformMethod method, const QVector &mapCoordinates, const QVector &layerCoordinates ) +QgsGcpTransformerInterface *QgsGcpTransformerInterface::createFromParameters( QgsGcpTransformerInterface::TransformMethod method, const QVector &sourceCoordinates, const QVector &destinationCoordinates ) { std::unique_ptr< QgsGcpTransformerInterface > transformer( create( method ) ); if ( !transformer ) return nullptr; - if ( !transformer->updateParametersFromGcps( mapCoordinates, layerCoordinates ) ) + if ( !transformer->updateParametersFromGcps( sourceCoordinates, destinationCoordinates ) ) return nullptr; return transformer.release(); @@ -122,11 +122,13 @@ QgsGcpTransformerInterface *QgsLinearGeorefTransform::clone() const return res.release(); } -bool QgsLinearGeorefTransform::updateParametersFromGcps( const QVector &mapCoords, const QVector &layerCoords ) +bool QgsLinearGeorefTransform::updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis ) { - if ( mapCoords.size() < minimumGcpCount() ) + if ( destinationCoordinates.size() < minimumGcpCount() ) return false; - QgsLeastSquares::linear( mapCoords, layerCoords, mParameters.origin, mParameters.scaleX, mParameters.scaleY ); + + mParameters.invertYAxis = invertYAxis; + QgsLeastSquares::linear( sourceCoordinates, destinationCoordinates, mParameters.origin, mParameters.scaleX, mParameters.scaleY ); return true; } @@ -163,7 +165,7 @@ int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstTo 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(); + y[i] = ( t->invertYAxis ? -1 : 1 ) * y[i] * t->scaleY + t->origin.y(); panSuccess[i] = true; } } @@ -182,7 +184,7 @@ int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstTo 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 ); + y[i] = ( y[i] - t->origin.y() ) / ( ( t->invertYAxis ? -1 : 1 ) * t->scaleY ); panSuccess[i] = true; } } @@ -193,12 +195,13 @@ int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstTo // // QgsHelmertGeorefTransform // -bool QgsHelmertGeorefTransform::updateParametersFromGcps( const QVector &mapCoords, const QVector &layerCoords ) +bool QgsHelmertGeorefTransform::updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis ) { - if ( mapCoords.size() < minimumGcpCount() ) + if ( destinationCoordinates.size() < minimumGcpCount() ) return false; - QgsLeastSquares::helmert( mapCoords, layerCoords, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle ); + mHelmertParameters.invertYAxis = invertYAxis; + QgsLeastSquares::helmert( sourceCoordinates, destinationCoordinates, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle ); return true; } @@ -207,10 +210,9 @@ int QgsHelmertGeorefTransform::minimumGcpCount() const return 2; } - GDALTransformerFunc QgsHelmertGeorefTransform::GDALTransformer() const { - return QgsHelmertGeorefTransform::helmert_transform; + return QgsHelmertGeorefTransform::helmertTransform; } void *QgsHelmertGeorefTransform::GDALTransformerArgs() const @@ -238,28 +240,42 @@ QgsGcpTransformerInterface *QgsHelmertGeorefTransform::clone() const return res.release(); } -int QgsHelmertGeorefTransform::helmert_transform( void *pTransformerArg, int bDstToSrc, int nPointCount, +int QgsHelmertGeorefTransform::helmertTransform( void *pTransformerArg, int bDstToSrc, int nPointCount, double *x, double *y, double *z, int *panSuccess ) { Q_UNUSED( z ) - HelmertParameters *t = static_cast( pTransformerArg ); + const HelmertParameters *t = static_cast< const 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; + double a = std::cos( t->angle ); + double b = std::sin( t->angle ); + const double x0 = t->origin.x(); + const double y0 = t->origin.y(); + const double 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 ); + const double xT = x[i]; + const double yT = y[i]; + + if ( t->invertYAxis ) + { + // 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| | cos a, sin a| + // |sin a, cos a | |0,-1| = | sin a, -cos a| + x[i] = x0 + ( a * xT + b * yT ); + y[i] = y0 + ( b * xT - a * yT ); + } + else + { + x[i] = x0 + ( a * xT - b * yT ); + y[i] = y0 + ( b * xT + a * yT ); + } panSuccess[i] = true; } } @@ -278,13 +294,20 @@ int QgsHelmertGeorefTransform::helmert_transform( void *pTransformerArg, int bDs 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; + const double xT = x[i] - x0; + const double yT = y[i] - y0; + if ( t->invertYAxis ) + { + // | cos a, sin a |^-1 |cos a, sin a| + // | sin a, -cos a | = |sin a, -cos a| + x[i] = a * xT + b * yT; + y[i] = b * xT - a * yT; + } + else + { + x[i] = a * xT + b * yT; + y[i] = -b * xT + a * yT; + } panSuccess[i] = true; } } @@ -309,19 +332,20 @@ QgsGDALGeorefTransform::~QgsGDALGeorefTransform() QgsGcpTransformerInterface *QgsGDALGeorefTransform::clone() const { std::unique_ptr< QgsGDALGeorefTransform > res = qgis::make_unique< QgsGDALGeorefTransform >( mIsTPSTransform, mPolynomialOrder ); - res->updateParametersFromGcps( mMapCoords, mLayerCoords ); + res->updateParametersFromGcps( mSourceCoords, mDestCoordinates, mInvertYAxis ); return res.release(); } -bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector &mapCoords, const QVector &layerCoords ) +bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis ) { - mMapCoords = mapCoords; - mLayerCoords = layerCoords; + mSourceCoords = sourceCoordinates; + mDestCoordinates = destinationCoordinates; + mInvertYAxis = invertYAxis; - assert( mapCoords.size() == layerCoords.size() ); - if ( mapCoords.size() != layerCoords.size() ) + assert( sourceCoordinates.size() == destinationCoordinates.size() ); + if ( sourceCoordinates.size() != destinationCoordinates.size() ) return false; - int n = mapCoords.size(); + int n = sourceCoordinates.size(); GDAL_GCP *GCPList = new GDAL_GCP[n]; for ( int i = 0; i < n; i++ ) @@ -329,10 +353,10 @@ bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector GCPList[i].pszId = new char[20]; snprintf( GCPList[i].pszId, 19, "gcp%i", i ); GCPList[i].pszInfo = nullptr; - GCPList[i].dfGCPPixel = layerCoords[i].x(); - GCPList[i].dfGCPLine = -layerCoords[i].y(); - GCPList[i].dfGCPX = mapCoords[i].x(); - GCPList[i].dfGCPY = mapCoords[i].y(); + GCPList[i].dfGCPPixel = sourceCoordinates[i].x(); + GCPList[i].dfGCPLine = ( mInvertYAxis ? -1 : 1 ) * sourceCoordinates[i].y(); + GCPList[i].dfGCPX = destinationCoordinates[i].x(); + GCPList[i].dfGCPY = destinationCoordinates[i].y(); GCPList[i].dfGCPZ = 0; } destroyGdalArgs(); @@ -419,20 +443,27 @@ QgsGcpTransformerInterface *QgsProjectiveGeorefTransform::clone() const return res.release(); } -bool QgsProjectiveGeorefTransform::updateParametersFromGcps( const QVector &mapCoords, const QVector &layerCoords ) +bool QgsProjectiveGeorefTransform::updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis ) { - if ( mapCoords.size() < minimumGcpCount() ) + if ( destinationCoordinates.size() < minimumGcpCount() ) return false; - // HACK: flip y coordinates, because georeferencer and gdal use different conventions - QVector flippedPixelCoords; - flippedPixelCoords.reserve( layerCoords.size() ); - for ( const QgsPointXY &coord : layerCoords ) + if ( invertYAxis ) { - flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() ); - } + // HACK: flip y coordinates, because georeferencer and gdal use different conventions + QVector flippedPixelCoords; + flippedPixelCoords.reserve( sourceCoordinates.size() ); + for ( const QgsPointXY &coord : sourceCoordinates ) + { + flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() ); + } - QgsLeastSquares::projective( mapCoords, flippedPixelCoords, mParameters.H ); + QgsLeastSquares::projective( flippedPixelCoords, destinationCoordinates, mParameters.H ); + } + else + { + QgsLeastSquares::projective( sourceCoordinates, destinationCoordinates, mParameters.H ); + } // Invert the homography matrix using adjoint matrix double *H = mParameters.H; diff --git a/src/analysis/georeferencing/qgsgcptransformer.h b/src/analysis/georeferencing/qgsgcptransformer.h index d66f2fc66d0..f3326912d36 100644 --- a/src/analysis/georeferencing/qgsgcptransformer.h +++ b/src/analysis/georeferencing/qgsgcptransformer.h @@ -49,6 +49,7 @@ class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT InvalidTransform = 65535 //!< Invalid transform }; + //! Constructor for QgsGcpTransformerInterface QgsGcpTransformerInterface() = default; virtual ~QgsGcpTransformerInterface() = default; @@ -68,11 +69,14 @@ class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT virtual QgsGcpTransformerInterface *clone() const = 0 SIP_FACTORY; /** - * Fits transformation parameters using the specified Ground Control Points (GCPs) lists of map coordinates and layer coordinates. + * Fits transformation parameters using the specified Ground Control Points (GCPs) lists of source and destination coordinates. + * + * If \a invertYAxis is set to TRUE then the y-axis of source coordinates will be inverted, e.g. to allow for transformation of raster layers + * with ascending top-to-bottom vertical axis coordinates. * * \returns TRUE on success, FALSE on failure */ - virtual bool updateParametersFromGcps( const QVector &mapCoordinates, const QVector &layerCoordinates ) = 0; + virtual bool updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis = false ) = 0; /** * Returns the minimum number of Ground Control Points (GCPs) required for parameter fitting. @@ -107,13 +111,13 @@ class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT /** * Creates a new QgsGcpTransformerInterface subclass representing the specified transform \a method, initialized - * using the given lists of map and layer coordinates. + * using the given lists of source and destination coordinates. * * If the parameters cannot be fit to a transform NULLPTR will be returned. * * Caller takes ownership of the returned object. */ - static QgsGcpTransformerInterface *createFromParameters( TransformMethod method, const QVector &mapCoordinates, const QVector &layerCoordinates ) SIP_FACTORY; + static QgsGcpTransformerInterface *createFromParameters( TransformMethod method, const QVector &sourceCoordinates, const QVector &destinationCoordinates ) SIP_FACTORY; #ifndef SIP_RUN @@ -144,6 +148,8 @@ class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP { public: + + //! Constructor for QgsLinearGeorefTransform QgsLinearGeorefTransform() = default; /** @@ -152,7 +158,7 @@ class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterfa bool getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const; QgsGcpTransformerInterface *clone() const override; - bool updateParametersFromGcps( const QVector &mapCoords, const QVector &layerCoords ) override; + bool updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis = false ) override; int minimumGcpCount() const override; GDALTransformerFunc GDALTransformer() const override; void *GDALTransformerArgs() const override; @@ -162,7 +168,9 @@ class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterfa struct LinearParameters { QgsPointXY origin; - double scaleX, scaleY; + double scaleX = 1; + double scaleY = 1; + bool invertYAxis = false; } mParameters; static int linearTransform( void *pTransformerArg, int bDstToSrc, int nPointCount, @@ -179,6 +187,8 @@ class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterfa class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP { public: + + //! Constructor for QgsHelmertGeorefTransform QgsHelmertGeorefTransform() = default; /** @@ -187,7 +197,7 @@ class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterf bool getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const; QgsGcpTransformerInterface *clone() const override; - bool updateParametersFromGcps( const QVector &mapCoords, const QVector &layerCoords ) override; + bool updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis = false ) override; int minimumGcpCount() const override; GDALTransformerFunc GDALTransformer() const override; void *GDALTransformerArgs() const override; @@ -200,11 +210,12 @@ class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterf QgsPointXY origin; double scale; double angle; + bool invertYAxis = false; }; HelmertParameters mHelmertParameters; - static int helmert_transform( void *pTransformerArg, int bDstToSrc, int nPointCount, - double *x, double *y, double *z, int *panSuccess ); + static int helmertTransform( void *pTransformerArg, int bDstToSrc, int nPointCount, + double *x, double *y, double *z, int *panSuccess ); }; @@ -217,11 +228,13 @@ class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterf class ANALYSIS_EXPORT QgsGDALGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP { public: + + //! Constructor for QgsGDALGeorefTransform QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder ); ~QgsGDALGeorefTransform() override; QgsGcpTransformerInterface *clone() const override; - bool updateParametersFromGcps( const QVector &mapCoords, const QVector &layerCoords ) override; + bool updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis = false ) override; int minimumGcpCount() const override; GDALTransformerFunc GDALTransformer() const override; void *GDALTransformerArgs() const override; @@ -230,13 +243,13 @@ class ANALYSIS_EXPORT QgsGDALGeorefTransform : public QgsGcpTransformerInterface private: void destroyGdalArgs(); - QVector mMapCoords; - QVector mLayerCoords; + QVector mSourceCoords; + QVector mDestCoordinates; + bool mInvertYAxis = false; const int mPolynomialOrder; const bool mIsTPSTransform; - GDALTransformerFunc mGDALTransformer = nullptr; void *mGDALTransformerArgs = nullptr; }; @@ -253,10 +266,12 @@ class ANALYSIS_EXPORT QgsGDALGeorefTransform : public QgsGcpTransformerInterface class ANALYSIS_EXPORT QgsProjectiveGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP { public: + + //! Constructor for QgsProjectiveGeorefTransform QgsProjectiveGeorefTransform(); QgsGcpTransformerInterface *clone() const override; - bool updateParametersFromGcps( const QVector &mapCoords, const QVector &layerCoords ) override; + bool updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis = false ) override; int minimumGcpCount() const override; GDALTransformerFunc GDALTransformer() const override; void *GDALTransformerArgs() const override; diff --git a/src/analysis/georeferencing/qgsleastsquares.cpp b/src/analysis/georeferencing/qgsleastsquares.cpp index 77f6a4252d7..dd74bc94411 100644 --- a/src/analysis/georeferencing/qgsleastsquares.cpp +++ b/src/analysis/georeferencing/qgsleastsquares.cpp @@ -23,11 +23,11 @@ #include "qgsleastsquares.h" -void QgsLeastSquares::linear( const QVector &mapCoords, - const QVector &pixelCoords, +void QgsLeastSquares::linear( const QVector &sourceCoordinates, + const QVector &destinationCoordinates, QgsPointXY &origin, double &pixelXSize, double &pixelYSize ) { - int n = mapCoords.size(); + int n = destinationCoordinates.size(); if ( n < 2 ) { throw std::domain_error( QObject::tr( "Fit to a linear transform requires at least 2 points." ).toLocal8Bit().constData() ); @@ -36,14 +36,14 @@ void QgsLeastSquares::linear( const QVector &mapCoords, double sumPx( 0 ), sumPy( 0 ), sumPx2( 0 ), sumPy2( 0 ), sumPxMx( 0 ), sumPyMy( 0 ), sumMx( 0 ), sumMy( 0 ); for ( int i = 0; i < n; ++i ) { - sumPx += pixelCoords.at( i ).x(); - sumPy += pixelCoords.at( i ).y(); - sumPx2 += std::pow( pixelCoords.at( i ).x(), 2 ); - sumPy2 += std::pow( pixelCoords.at( i ).y(), 2 ); - sumPxMx += pixelCoords.at( i ).x() * mapCoords.at( i ).x(); - sumPyMy += pixelCoords.at( i ).y() * mapCoords.at( i ).y(); - sumMx += mapCoords.at( i ).x(); - sumMy += mapCoords.at( i ).y(); + sumPx += sourceCoordinates.at( i ).x(); + sumPy += sourceCoordinates.at( i ).y(); + sumPx2 += std::pow( sourceCoordinates.at( i ).x(), 2 ); + sumPy2 += std::pow( sourceCoordinates.at( i ).y(), 2 ); + sumPxMx += sourceCoordinates.at( i ).x() * destinationCoordinates.at( i ).x(); + sumPyMy += sourceCoordinates.at( i ).y() * destinationCoordinates.at( i ).y(); + sumMx += destinationCoordinates.at( i ).x(); + sumMy += destinationCoordinates.at( i ).y(); } double deltaX = n * sumPx2 - std::pow( sumPx, 2 ); @@ -62,30 +62,39 @@ void QgsLeastSquares::linear( const QVector &mapCoords, } -void QgsLeastSquares::helmert( const QVector &mapCoords, - const QVector &pixelCoords, +void QgsLeastSquares::helmert( const QVector &sourceCoordinates, + const QVector &destinationCoordinates, QgsPointXY &origin, double &pixelSize, double &rotation ) { - int n = mapCoords.size(); + int n = destinationCoordinates.size(); if ( n < 2 ) { throw std::domain_error( QObject::tr( "Fit to a Helmert transform requires at least 2 points." ).toLocal8Bit().constData() ); } - double A = 0, B = 0, C = 0, D = 0, E = 0, F = 0, G = 0, H = 0, I = 0, J = 0; + double A = 0; + double B = 0; + double C = 0; + double D = 0; + double E = 0; + double F = 0; + double G = 0; + double H = 0; + double I = 0; + double J = 0; for ( int i = 0; i < n; ++i ) { - A += pixelCoords.at( i ).x(); - B += pixelCoords.at( i ).y(); - C += mapCoords.at( i ).x(); - D += mapCoords.at( i ).y(); - E += mapCoords.at( i ).x() * pixelCoords.at( i ).x(); - F += mapCoords.at( i ).y() * pixelCoords.at( i ).y(); - G += std::pow( pixelCoords.at( i ).x(), 2 ); - H += std::pow( pixelCoords.at( i ).y(), 2 ); - I += mapCoords.at( i ).x() * pixelCoords.at( i ).y(); - J += pixelCoords.at( i ).x() * mapCoords.at( i ).y(); + A += sourceCoordinates.at( i ).x(); + B += sourceCoordinates.at( i ).y(); + C += destinationCoordinates.at( i ).x(); + D += destinationCoordinates.at( i ).y(); + E += destinationCoordinates.at( i ).x() * sourceCoordinates.at( i ).x(); + F += destinationCoordinates.at( i ).y() * sourceCoordinates.at( i ).y(); + G += std::pow( sourceCoordinates.at( i ).x(), 2 ); + H += std::pow( sourceCoordinates.at( i ).y(), 2 ); + I += destinationCoordinates.at( i ).x() * sourceCoordinates.at( i ).y(); + J += sourceCoordinates.at( i ).x() * destinationCoordinates.at( i ).y(); } /* The least squares fit for the parameters { a, b, x0, y0 } is the solution @@ -232,61 +241,59 @@ void normalizeCoordinates( const QVector &coords, QVector mapCoords, - QVector pixelCoords, +void QgsLeastSquares::projective( const QVector &sourceCoordinates, + const QVector &destinationCoordinates, double H[9] ) { - Q_ASSERT( mapCoords.size() == pixelCoords.size() ); + Q_ASSERT( sourceCoordinates.size() == destinationCoordinates.size() ); - if ( mapCoords.size() < 4 ) + if ( destinationCoordinates.size() < 4 ) { throw std::domain_error( QObject::tr( "Fitting a projective transform requires at least 4 corresponding points." ).toLocal8Bit().constData() ); } - QVector mapCoordsNormalized; - QVector pixelCoordsNormalized; + QVector sourceCoordinatesNormalized; + QVector destinationCoordinatesNormalized; - double normMap[9], denormMap[9]; - double normPixel[9], denormPixel[9]; - normalizeCoordinates( mapCoords, mapCoordsNormalized, normMap, denormMap ); - normalizeCoordinates( pixelCoords, pixelCoordsNormalized, normPixel, denormPixel ); - mapCoords = mapCoordsNormalized; - pixelCoords = pixelCoordsNormalized; + double normSource[9], denormSource[9]; + double normDest[9], denormDest[9]; + normalizeCoordinates( sourceCoordinates, sourceCoordinatesNormalized, normSource, denormSource ); + normalizeCoordinates( destinationCoordinates, destinationCoordinatesNormalized, normDest, denormDest ); // GSL does not support a full SVD, so we artificially add a linear dependent row // to the matrix in case the system is underconstrained. - uint m = std::max( 9u, ( uint )mapCoords.size() * 2u ); + uint m = std::max( 9u, ( uint )destinationCoordinatesNormalized.size() * 2u ); uint n = 9; gsl_matrix *S = gsl_matrix_alloc( m, n ); - for ( int i = 0; i < mapCoords.size(); i++ ) + for ( int i = 0; i < destinationCoordinatesNormalized.size(); i++ ) { - gsl_matrix_set( S, i * 2, 0, pixelCoords[i].x() ); - gsl_matrix_set( S, i * 2, 1, pixelCoords[i].y() ); + gsl_matrix_set( S, i * 2, 0, sourceCoordinatesNormalized[i].x() ); + gsl_matrix_set( S, i * 2, 1, sourceCoordinatesNormalized[i].y() ); gsl_matrix_set( S, i * 2, 2, 1.0 ); gsl_matrix_set( S, i * 2, 3, 0.0 ); gsl_matrix_set( S, i * 2, 4, 0.0 ); gsl_matrix_set( S, i * 2, 5, 0.0 ); - gsl_matrix_set( S, i * 2, 6, -mapCoords[i].x()*pixelCoords[i].x() ); - gsl_matrix_set( S, i * 2, 7, -mapCoords[i].x()*pixelCoords[i].y() ); - gsl_matrix_set( S, i * 2, 8, -mapCoords[i].x() * 1.0 ); + gsl_matrix_set( S, i * 2, 6, -destinationCoordinatesNormalized[i].x()*sourceCoordinatesNormalized[i].x() ); + gsl_matrix_set( S, i * 2, 7, -destinationCoordinatesNormalized[i].x()*sourceCoordinatesNormalized[i].y() ); + gsl_matrix_set( S, i * 2, 8, -destinationCoordinatesNormalized[i].x() * 1.0 ); gsl_matrix_set( S, i * 2 + 1, 0, 0.0 ); gsl_matrix_set( S, i * 2 + 1, 1, 0.0 ); gsl_matrix_set( S, i * 2 + 1, 2, 0.0 ); - gsl_matrix_set( S, i * 2 + 1, 3, pixelCoords[i].x() ); - gsl_matrix_set( S, i * 2 + 1, 4, pixelCoords[i].y() ); + gsl_matrix_set( S, i * 2 + 1, 3, sourceCoordinatesNormalized[i].x() ); + gsl_matrix_set( S, i * 2 + 1, 4, sourceCoordinatesNormalized[i].y() ); gsl_matrix_set( S, i * 2 + 1, 5, 1.0 ); - gsl_matrix_set( S, i * 2 + 1, 6, -mapCoords[i].y()*pixelCoords[i].x() ); - gsl_matrix_set( S, i * 2 + 1, 7, -mapCoords[i].y()*pixelCoords[i].y() ); - gsl_matrix_set( S, i * 2 + 1, 8, -mapCoords[i].y() * 1.0 ); + gsl_matrix_set( S, i * 2 + 1, 6, -destinationCoordinatesNormalized[i].y()*sourceCoordinatesNormalized[i].x() ); + gsl_matrix_set( S, i * 2 + 1, 7, -destinationCoordinatesNormalized[i].y()*sourceCoordinatesNormalized[i].y() ); + gsl_matrix_set( S, i * 2 + 1, 8, -destinationCoordinatesNormalized[i].y() * 1.0 ); } - if ( mapCoords.size() == 4 ) + if ( destinationCoordinatesNormalized.size() == 4 ) { // The GSL SVD routine only supports matrices with rows >= columns (m >= n) // Unfortunately, we can't use the SVD of the transpose (i.e. S^T = (U D V^T)^T = V D U^T) @@ -321,13 +328,13 @@ void QgsLeastSquares::projective( QVector mapCoords, gsl_matrix *prodMatrix = gsl_matrix_alloc( 3, 3 ); gsl_matrix_view Hmatrix = gsl_matrix_view_array( H, 3, 3 ); - gsl_matrix_view normPixelMatrix = gsl_matrix_view_array( normPixel, 3, 3 ); - gsl_matrix_view denormMapMatrix = gsl_matrix_view_array( denormMap, 3, 3 ); + gsl_matrix_view normSourceMatrix = gsl_matrix_view_array( normSource, 3, 3 ); + gsl_matrix_view denormDestMatrix = gsl_matrix_view_array( denormDest, 3, 3 ); - // Change coordinate frame of image and pre-image from normalized to map and pixel coordinates. + // Change coordinate frame of image and pre-image from normalized to destination and source coordinates. // H' = denormalizeMapCoords*H*normalizePixelCoords - gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &Hmatrix.matrix, &normPixelMatrix.matrix, 0.0, prodMatrix ); - gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &denormMapMatrix.matrix, prodMatrix, 0.0, &Hmatrix.matrix ); + gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &Hmatrix.matrix, &normSourceMatrix.matrix, 0.0, prodMatrix ); + gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &denormDestMatrix.matrix, prodMatrix, 0.0, &Hmatrix.matrix ); gsl_matrix_free( prodMatrix ); gsl_matrix_free( S ); diff --git a/src/analysis/georeferencing/qgsleastsquares.h b/src/analysis/georeferencing/qgsleastsquares.h index 3d87cd949a5..b09e0e37c64 100644 --- a/src/analysis/georeferencing/qgsleastsquares.h +++ b/src/analysis/georeferencing/qgsleastsquares.h @@ -34,17 +34,17 @@ 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). + * Transforms the point at \a origin in-place, using a linear transformation calculated from the list of source and destination Ground Control Points (GCPs). */ - static void linear( const QVector &mapCoords, - const QVector &pixelCoords, + static void linear( const QVector &sourceCoordinates, + const QVector &destinationCoordinates, 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). + * Transforms the point at \a origin in-place, using a helmert transformation calculated from the list of source and destination Ground Control Points (GCPs). */ - static void helmert( const QVector &mapCoords, - const QVector &pixelCoords, + static void helmert( const QVector &sourceCoordinates, + const QVector &destinationCoordinates, QgsPointXY &origin, double &pixelSize, double &rotation ); #if 0 @@ -54,10 +54,10 @@ class ANALYSIS_EXPORT QgsLeastSquares #endif /** - * Calculates projective parameters from the list of map and pixel Ground Control Points (GCPs). + * Calculates projective parameters from the list of source and destination Ground Control Points (GCPs). */ - static void projective( QVector mapCoords, - QVector pixelCoords, + static void projective( const QVector &sourceCoordinates, + const QVector &destinationCoordinates, double H[9] ); }; diff --git a/src/app/georeferencer/qgsgcplistmodel.cpp b/src/app/georeferencer/qgsgcplistmodel.cpp index 016de2b3f25..14c8920192c 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( pixelCoords, mapCoords, true ); mapUnitsPossible = mGeorefTransform->providesAccurateInverseTransformation(); } diff --git a/src/app/georeferencer/qgsgeorefmainwindow.cpp b/src/app/georeferencer/qgsgeorefmainwindow.cpp index 7038480a20a..50d9f5fdfab 100644 --- a/src/app/georeferencer/qgsgeorefmainwindow.cpp +++ b/src/app/georeferencer/qgsgeorefmainwindow.cpp @@ -2004,7 +2004,7 @@ bool QgsGeoreferencerMainWindow::updateGeorefTransform() return false; // Parametrize the transform with GCPs - if ( !mGeorefTransform.updateParametersFromGcps( mapCoords, pixelCoords ) ) + if ( !mGeorefTransform.updateParametersFromGcps( pixelCoords, mapCoords, true ) ) { return false; } diff --git a/src/app/georeferencer/qgsgeoreftransform.cpp b/src/app/georeferencer/qgsgeoreftransform.cpp index b7b9024819c..d29f1ee1816 100644 --- a/src/app/georeferencer/qgsgeoreftransform.cpp +++ b/src/app/georeferencer/qgsgeoreftransform.cpp @@ -73,36 +73,36 @@ bool QgsGeorefTransform::parametersInitialized() const QgsGcpTransformerInterface *QgsGeorefTransform::clone() const { std::unique_ptr< QgsGeorefTransform > res( new QgsGeorefTransform( *this ) ); - res->updateParametersFromGcps( mMapCoords, mLayerCoords ); + res->updateParametersFromGcps( mSourceCoordinates, mDestinationCoordinates, mInvertYAxis ); return res.release(); } -bool QgsGeorefTransform::updateParametersFromGcps( const QVector &mapCoords, const QVector &pixelCoords ) +bool QgsGeorefTransform::updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis ) { - mMapCoords = mapCoords; - mLayerCoords = pixelCoords; + mSourceCoordinates = sourceCoordinates; + mDestinationCoordinates = destinationCoordinates; + mInvertYAxis = invertYAxis; if ( !mGeorefTransformImplementation ) { return false; } - if ( mapCoords.size() != pixelCoords.size() ) // Defensive sanity check + if ( sourceCoordinates.size() != destinationCoordinates.size() ) // Defensive sanity check { throw ( std::domain_error( "Internal error: GCP mapping is not one-to-one" ) ); } - if ( mapCoords.size() < minimumGcpCount() ) + if ( sourceCoordinates.size() < minimumGcpCount() ) { return false; } if ( mRasterChangeCoords.hasCrs() ) { - QVector pixelCoordsCorrect = mRasterChangeCoords.getPixelCoords( pixelCoords ); - mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( mapCoords, pixelCoordsCorrect ); - pixelCoordsCorrect.clear(); + QVector pixelCoordsCorrect = mRasterChangeCoords.getPixelCoords( sourceCoordinates ); + mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( sourceCoordinates, pixelCoordsCorrect, invertYAxis ); } else { - mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( mapCoords, pixelCoords ); + mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( sourceCoordinates, destinationCoordinates, invertYAxis ); } return mParametersInitialized; } diff --git a/src/app/georeferencer/qgsgeoreftransform.h b/src/app/georeferencer/qgsgeoreftransform.h index 457ccf52b60..9d67f118f9e 100644 --- a/src/app/georeferencer/qgsgeoreftransform.h +++ b/src/app/georeferencer/qgsgeoreftransform.h @@ -70,7 +70,7 @@ class QgsGeorefTransform : public QgsGcpTransformerInterface bool parametersInitialized() const; QgsGcpTransformerInterface *clone() const override; - bool updateParametersFromGcps( const QVector &mapCoords, const QVector &pixelCoords ) override; + bool updateParametersFromGcps( const QVector &sourceCoordinates, const QVector &destinationCoordinates, bool invertYAxis = false ) override; int minimumGcpCount() const override; TransformMethod method() const override; GDALTransformerFunc GDALTransformer() const override; @@ -112,8 +112,9 @@ class QgsGeorefTransform : public QgsGcpTransformerInterface // convenience wrapper around GDALTransformerFunc bool gdal_transform( const QgsPointXY &src, QgsPointXY &dst, int dstToSrc ) const; - QVector mMapCoords; - QVector mLayerCoords; + QVector mSourceCoordinates; + QVector mDestinationCoordinates; + bool mInvertYAxis = false; std::unique_ptr< QgsGcpTransformerInterface > mGeorefTransformImplementation; diff --git a/tests/src/analysis/CMakeLists.txt b/tests/src/analysis/CMakeLists.txt index 487704e3399..b7410e1db2f 100644 --- a/tests/src/analysis/CMakeLists.txt +++ b/tests/src/analysis/CMakeLists.txt @@ -46,6 +46,7 @@ endmacro (ADD_QGIS_TEST) ############################################################# # Tests: set(TESTS + testqgsgcptransformer.cpp testqgsgeometrysnapper.cpp testqgsinterpolator.cpp testqgsprocessing.cpp diff --git a/tests/src/analysis/testqgsgcptransformer.cpp b/tests/src/analysis/testqgsgcptransformer.cpp index 475eccab987..a492bc88b32 100644 --- a/tests/src/analysis/testqgsgcptransformer.cpp +++ b/tests/src/analysis/testqgsgcptransformer.cpp @@ -43,29 +43,78 @@ class TestQgsGcpTransformer : public QObject void testLinear() { QgsLinearGeorefTransform transform; - transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, 1126 ) + << QgsPointXY( 2352, 716 ) + << QgsPointXY( 2067, 2398 ) + << QgsPointXY( 1102, 2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) << QgsPointXY( 322698, 129979 ) << QgsPointXY( 322501, 192061 ) - << QgsPointXY( 321803, 129198 ), - QVector< QgsPointXY >() << QgsPointXY( 288, -716 ) - << QgsPointXY( 2352, -1126 ) - << QgsPointXY( 2067, -2398 ) - << QgsPointXY( 1102, -2209 ) ); + << QgsPointXY( 321803, 192198 ) ); double x = 288; - double y = -71; + double y = 1000; QVERIFY( transform.transform( x, y ) ); QGSCOMPARENEAR( x, 321212, 1 ); - QGSCOMPARENEAR( y, 108597, 1 ); + QGSCOMPARENEAR( y, 135047, 1 ); QVERIFY( transform.transform( x, y, true ) ); QGSCOMPARENEAR( x, 288, 1 ); - QGSCOMPARENEAR( y, -71, 1 ); + QGSCOMPARENEAR( y, 1000, 1 ); + + x = 2352; + y = 1150; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 141437, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2352, 1 ); + QGSCOMPARENEAR( y, 1150, 1 ); + + x = 2067; + y = 2500; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 198947, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2067, 1 ); + QGSCOMPARENEAR( y, 2500, 1 ); + + x = 1102; + y = 2300; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 190427, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1102, 1 ); + QGSCOMPARENEAR( y, 2300, 1 ); + } + + void testLinearRasterYAxis() + { + QgsLinearGeorefTransform transform; + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, -716 ) + << QgsPointXY( 2352, -1126 ) + << QgsPointXY( 2067, -2398 ) + << QgsPointXY( 1102, -2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 288; + double y = -716; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 288, 1 ); + QGSCOMPARENEAR( y, -716, 1 ); x = 2352; y = -1126; QVERIFY( transform.transform( x, y ) ); QGSCOMPARENEAR( x, 322702, 1 ); - QGSCOMPARENEAR( y, 133775, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); QVERIFY( transform.transform( x, y, true ) ); QGSCOMPARENEAR( x, 2352, 1 ); QGSCOMPARENEAR( y, -1126, 1 ); @@ -74,7 +123,7 @@ class TestQgsGcpTransformer : public QObject y = -2398; QVERIFY( transform.transform( x, y ) ); QGSCOMPARENEAR( x, 322496, 1 ); - QGSCOMPARENEAR( y, 164131, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); QVERIFY( transform.transform( x, y, true ) ); QGSCOMPARENEAR( x, 2067, 1 ); QGSCOMPARENEAR( y, -2398, 1 ); @@ -83,12 +132,1398 @@ class TestQgsGcpTransformer : public QObject y = -2209; QVERIFY( transform.transform( x, y ) ); QGSCOMPARENEAR( x, 321800, 1 ); - QGSCOMPARENEAR( y, 159620, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); QVERIFY( transform.transform( x, y, true ) ); QGSCOMPARENEAR( x, 1102, 1 ); QGSCOMPARENEAR( y, -2209, 1 ); } + void testLinearExact() + { + QgsLinearGeorefTransform transform; + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + void testLinearExactRasterYAxis() + { + QgsLinearGeorefTransform transform; + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, -123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, -140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, -194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, -186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + + // + // Helmert + // + + void testHelmert() + { + QgsHelmertGeorefTransform transform; + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, 1126 ) + << QgsPointXY( 2352, 716 ) + << QgsPointXY( 2067, 2398 ) + << QgsPointXY( 1102, 2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 288; + double y = 1000; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 302325, 1 ); + QGSCOMPARENEAR( y, 145676, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 288, 1 ); + QGSCOMPARENEAR( y, 1000, 1 ); + + x = 2352; + y = 1150; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 340495, 1 ); + QGSCOMPARENEAR( y, 155540, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2352, 1 ); + QGSCOMPARENEAR( y, 1150, 1 ); + + x = 2067; + y = 2500; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 330540, 1 ); + QGSCOMPARENEAR( y, 179867, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2067, 1 ); + QGSCOMPARENEAR( y, 2500, 1 ); + + x = 1102; + y = 2300; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 313138, 1 ); + QGSCOMPARENEAR( y, 172822, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1102, 1 ); + QGSCOMPARENEAR( y, 2300, 1 ); + } + + void testHelmertRasterYAxis() + { + QgsHelmertGeorefTransform transform; + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, -716 ) + << QgsPointXY( 2352, -1126 ) + << QgsPointXY( 2067, -2398 ) + << QgsPointXY( 1102, -2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 288; + double y = -716; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 334590, 1 ); + QGSCOMPARENEAR( y, 115324, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 288, 1 ); + QGSCOMPARENEAR( y, -716, 1 ); + + x = 2352; + y = -1126; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 296200, 1 ); + QGSCOMPARENEAR( y, 115340, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2352, 1 ); + QGSCOMPARENEAR( y, -1126, 1 ); + + x = 2067; + y = -2398; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 296768, 1 ); + QGSCOMPARENEAR( y, 91566, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2067, 1 ); + QGSCOMPARENEAR( y, -2398, 1 ); + + x = 1102; + y = -2209; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 314707, 1 ); + QGSCOMPARENEAR( y, 91510, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1102, 1 ); + QGSCOMPARENEAR( y, -2209, 1 ); + } + + void testHelmertExact() + { + QgsHelmertGeorefTransform transform; + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + void testHelmertExactRasterYAxis() + { + QgsHelmertGeorefTransform transform; + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, -123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, -140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, -194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, -186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + + + // + // Polynomial order 1 + // + + void testPolyOrder1() + { + QgsGDALGeorefTransform transform( false, 1 ); + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, 1126 ) + << QgsPointXY( 2352, 716 ) + << QgsPointXY( 2067, 2398 ) + << QgsPointXY( 1102, 2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 288; + double y = 1000; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321209, 1 ); + QGSCOMPARENEAR( y, 128732, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 288, 1 ); + QGSCOMPARENEAR( y, 1017, 1 ); + + x = 2352; + y = 1150; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322700, 1 ); + QGSCOMPARENEAR( y, 146403, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2352, 1 ); + QGSCOMPARENEAR( y, 1164, 1 ); + + x = 2067; + y = 2500; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322501, 1 ); + QGSCOMPARENEAR( y, 202230, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2067, 1 ); + QGSCOMPARENEAR( y, 2474, 1 ); + + x = 1102; + y = 2300; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321803, 1 ); + QGSCOMPARENEAR( y, 188448, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1102, 1 ); + QGSCOMPARENEAR( y, 2279, 1 ); + } + + void testPolyOrder1RasterYAxis() + { + QgsGDALGeorefTransform transform( false, 1 ); + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, -716 ) + << QgsPointXY( 2352, -1126 ) + << QgsPointXY( 2067, -2398 ) + << QgsPointXY( 1102, -2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 288; + double y = -716; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321201, 1 ); + QGSCOMPARENEAR( y, 63463, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 288, 1 ); + QGSCOMPARENEAR( y, -716, 1 ); + + x = 2352; + y = -1126; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322686, 1 ); + QGSCOMPARENEAR( y, 24965, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2352, 1 ); + QGSCOMPARENEAR( y, -1126, 1 ); + + x = 2067; + y = -2398; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322474, 1 ); + QGSCOMPARENEAR( y, -31687, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2067, 1 ); + QGSCOMPARENEAR( y, -2398, 1 ); + + x = 1102; + y = -2209; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321780, 1 ); + QGSCOMPARENEAR( y, -13813, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1102, 1 ); + QGSCOMPARENEAR( y, -2209, 1 ); + } + + void testPolyOrder1Exact() + { + QgsGDALGeorefTransform transform( false, 1 ); + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + void testPolyOrder1ExactRasterYAxis() + { + QgsGDALGeorefTransform transform( false, 1 ); + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, -123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, -140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, -194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, -186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + + // + // Polynomial order 2 + // + + void testPolyOrder2() + { + QgsGDALGeorefTransform transform( false, 2 ); + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, 1126 ) + << QgsPointXY( 2352, 716 ) + << QgsPointXY( 2067, 2398 ) + << QgsPointXY( 2100, 1500 ) + << QgsPointXY( 1102, 2209 ) + << QgsPointXY( 300, 1550 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322550, 149979 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321310, 145979 ) ); + + double x = 288; + double y = 1000; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321175, 1 ); + QGSCOMPARENEAR( y, 126837, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 265, 1 ); + QGSCOMPARENEAR( y, 1025, 1 ); + + x = 2352; + y = 1150; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322763, 1 ); + QGSCOMPARENEAR( y, 135026, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2336, 1 ); + QGSCOMPARENEAR( y, 936, 1 ); + + x = 2067; + y = 2500; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322489, 1 ); + QGSCOMPARENEAR( y, 198409, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2097, 1 ); + QGSCOMPARENEAR( y, 2445, 1 ); + + x = 1102; + y = 2300; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321803, 1 ); + QGSCOMPARENEAR( y, 197778, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1123, 1 ); + QGSCOMPARENEAR( y, 2225, 1 ); + } + + void testPolyOrder2RasterYAxis() + { + QgsGDALGeorefTransform transform( false, 2 ); + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, -1126 ) + << QgsPointXY( 2352, -716 ) + << QgsPointXY( 2067, -2398 ) + << QgsPointXY( 2100, -1500 ) + << QgsPointXY( 1102, -2209 ) + << QgsPointXY( 300, -1550 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322550, 149979 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321310, 145979 ), true ); + + double x = 288; + double y = -716; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 320404, 1 ); + QGSCOMPARENEAR( y, 131061, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -1913, 1 ); + QGSCOMPARENEAR( y, 1625, 1 ); + + x = 2352; + y = -1126; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322022, 1 ); + QGSCOMPARENEAR( y, 176342, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1435, 1 ); + QGSCOMPARENEAR( y, 2087, 1 ); + + x = 2067; + y = -2398; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 320831, 1 ); + QGSCOMPARENEAR( y, 273376, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -310, 1 ); + QGSCOMPARENEAR( y, 134, 1 ); + + x = 1102; + y = -2209; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 319958, 1 ); + QGSCOMPARENEAR( y, 243375, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -3994, 1 ); + QGSCOMPARENEAR( y, 1158, 1 ); + } + + void testPolyOrder2Exact() + { + QgsGDALGeorefTransform transform( false, 2 ); + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322501, 152061 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321210, 152061 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322501, 152061 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321210, 152061 ) ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + void testPolyOrder2ExactRasterYAxis() + { + QgsGDALGeorefTransform transform( false, 2 ); + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322501, 152061 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321210, 152061 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322501, 152061 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321210, 152061 ), true ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, -123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, -140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, -194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, -186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + + // + // Polynomial order 3 + // + + void testPolyOrder3() + { + QgsGDALGeorefTransform transform( false, 3 ); + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, 1126 ) + << QgsPointXY( 2352, 716 ) + << QgsPointXY( 2067, 2398 ) + << QgsPointXY( 2100, 1500 ) + << QgsPointXY( 1102, 2209 ) + << QgsPointXY( 300, 1550 ) + << QgsPointXY( 300, 850 ) + << QgsPointXY( 1000, 830 ) + << QgsPointXY( 900, 1450 ) + << QgsPointXY( 700, 1550 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322550, 149979 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321310, 145979 ) + << QgsPointXY( 321310, 131979 ) + << QgsPointXY( 321703, 131579 ) + << QgsPointXY( 321603, 145979 ) + << QgsPointXY( 321003, 146179 ) + ); + + double x = 288; + double y = 1000; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321202, 1 ); + QGSCOMPARENEAR( y, 129364, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 237, 1 ); + QGSCOMPARENEAR( y, 1319, 1 ); + + x = 2352; + y = 1150; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 318595, 1 ); + QGSCOMPARENEAR( y, 101765, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 138412, 1 ); + QGSCOMPARENEAR( y, 15444, 1 ); + + x = 2067; + y = 2500; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321488, 1 ); + QGSCOMPARENEAR( y, 185766, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -1237, 1 ); + QGSCOMPARENEAR( y, 5315, 1 ); + + x = 1102; + y = 2300; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 320879, 1 ); + QGSCOMPARENEAR( y, 189391, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -4000, 1 ); + QGSCOMPARENEAR( y, 2621, 1 ); + } + + void testPolyOrder3RasterYAxis() + { + QgsGDALGeorefTransform transform( false, 3 ); + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, -1126 ) + << QgsPointXY( 2352, -716 ) + << QgsPointXY( 2067, -2398 ) + << QgsPointXY( 2100, -1500 ) + << QgsPointXY( 1102, -2209 ) + << QgsPointXY( 300, -1550 ) + << QgsPointXY( 300, -850 ) + << QgsPointXY( 1000, -830 ) + << QgsPointXY( 900, -1450 ) + << QgsPointXY( 700, -1550 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322550, 149979 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321310, 145979 ) + << QgsPointXY( 321310, 131979 ) + << QgsPointXY( 321703, 131579 ) + << QgsPointXY( 321603, 145979 ) + << QgsPointXY( 321003, 146179 ), true ); + + // these values look like nonsense... I can't verify them, except to say that at the time + // these tests were written the raster georeferencer worked correctly with polynomial order 3, so I can + // only assume they are correct...! + double x = 288; + double y = -716; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 368688, 1 ); + QGSCOMPARENEAR( y, 828425, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -476460174, 1 ); + QGSCOMPARENEAR( y, -43384555, 1 ); + + x = 2352; + y = -1126; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 456596, 1 ); + QGSCOMPARENEAR( y, 1609794, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -11099349683, 1 ); + QGSCOMPARENEAR( y, -130456708, 1 ); + + x = 2067; + y = -2398; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 731809, 1 ); + QGSCOMPARENEAR( y, 4861338, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -313632670265, 1 ); + QGSCOMPARENEAR( y, -7338651765, 1 ); + + x = 1102; + y = -2209; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 620013, 1 ); + QGSCOMPARENEAR( y, 3875762, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, -121400397631, 1 ); + QGSCOMPARENEAR( y, -5274501643, 1 ); + } + + void testPolyOrder3Exact() + { + QgsGDALGeorefTransform transform( false, 3 ); + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322501, 152061 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321210, 152061 ) + << QgsPointXY( 322010, 182061 ) + << QgsPointXY( 322010, 132061 ) + << QgsPointXY( 321050, 162061 ) + << QgsPointXY( 321050, 172061 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322501, 152061 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321210, 152061 ) + << QgsPointXY( 322010, 182061 ) + << QgsPointXY( 322010, 132061 ) + << QgsPointXY( 321050, 162061 ) + << QgsPointXY( 321050, 172061 ) ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + void testPolyOrder3ExactRasterYAxis() + { + QgsGDALGeorefTransform transform( false, 3 ); + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322501, 152061 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321210, 152061 ) + << QgsPointXY( 322010, 182061 ) + << QgsPointXY( 322010, 132061 ) + << QgsPointXY( 321050, 162061 ) + << QgsPointXY( 321050, 172061 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 322501, 152061 ) + << QgsPointXY( 321803, 192198 ) + << QgsPointXY( 321210, 152061 ) + << QgsPointXY( 322010, 182061 ) + << QgsPointXY( 322010, 132061 ) + << QgsPointXY( 321050, 162061 ) + << QgsPointXY( 321050, 172061 ), true ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, -123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, -140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, -194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, -186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + + // + // Thin plate spline + // + + void testTPSReports() + { + QgsGDALGeorefTransform transform( true, 1 ); + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, 1126 ) + << QgsPointXY( 2352, 716 ) + << QgsPointXY( 2067, 2398 ) + << QgsPointXY( 1102, 2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 288; + double y = 1000; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321209, 1 ); + QGSCOMPARENEAR( y, 124248, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 290, 1 ); + QGSCOMPARENEAR( y, 680, 1 ); + + x = 2352; + y = 1150; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322700, 1 ); + QGSCOMPARENEAR( y, 146609, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2347, 1 ); + QGSCOMPARENEAR( y, 1636, 1 ); + + x = 2067; + y = 2500; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322501, 1 ); + QGSCOMPARENEAR( y, 196160, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2068, 1 ); + QGSCOMPARENEAR( y, 2273, 1 ); + + x = 1102; + y = 2300; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321803, 1 ); + QGSCOMPARENEAR( y, 195999, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1103, 1 ); + QGSCOMPARENEAR( y, 2177, 1 ); + } + + void testTPSRasterYAxis() + { + QgsGDALGeorefTransform transform( true, 1 ); + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, -716 ) + << QgsPointXY( 2352, -1126 ) + << QgsPointXY( 2067, -2398 ) + << QgsPointXY( 1102, -2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 288; + double y = -716; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321202, 1 ); + QGSCOMPARENEAR( y, 63505, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 297, 1 ); + QGSCOMPARENEAR( y, -705, 1 ); + + x = 2352; + y = -1126; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322686, 1 ); + QGSCOMPARENEAR( y, 24959, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2361, 1 ); + QGSCOMPARENEAR( y, -1116, 1 ); + + x = 2067; + y = -2398; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322474, 1 ); + QGSCOMPARENEAR( y, -31676, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2078, 1 ); + QGSCOMPARENEAR( y, -2385, 1 ); + + x = 1102; + y = -2209; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321780, 1 ); + QGSCOMPARENEAR( y, -13787, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1114, 1 ); + QGSCOMPARENEAR( y, -2195, 1 ); + } + + void testTPSExact() + { + QgsGDALGeorefTransform transform( true, 1 ); + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + void testTPSExactRasterYAxis() + { + QgsGDALGeorefTransform transform( true, 1 ); + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, -123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, -140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, -194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, -186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + // + // Projective + // + + void testProjective() + { + QgsProjectiveGeorefTransform transform; + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, 1126 ) + << QgsPointXY( 2352, 716 ) + << QgsPointXY( 2067, 2398 ) + << QgsPointXY( 1102, 2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 288; + double y = 1000; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321256, 1 ); + QGSCOMPARENEAR( y, 123764, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 288, 1 ); + QGSCOMPARENEAR( y, 1000, 1 ); + + x = 2352; + y = 1150; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322687, 1 ); + QGSCOMPARENEAR( y, 142908, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2352, 1 ); + QGSCOMPARENEAR( y, 1150, 1 ); + + x = 2067; + y = 2500; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322494, 1 ); + QGSCOMPARENEAR( y, 197069, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2067, 1 ); + QGSCOMPARENEAR( y, 2500, 1 ); + + x = 1102; + y = 2300; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321781, 1 ); + QGSCOMPARENEAR( y, 198051, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1102, 1 ); + QGSCOMPARENEAR( y, 2300, 1 ); + } + + void testProjectiveRasterYAxis() + { + QgsProjectiveGeorefTransform transform; + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 288, -716 ) + << QgsPointXY( 2352, -1126 ) + << QgsPointXY( 2067, -2398 ) + << QgsPointXY( 1102, -2209 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 288; + double y = -716; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321203, 1 ); + QGSCOMPARENEAR( y, 63965, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 288, 1 ); + QGSCOMPARENEAR( y, -716, 1 ); + + x = 2352; + y = -1126; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322681, 1 ); + QGSCOMPARENEAR( y, 25366, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2352, 1 ); + QGSCOMPARENEAR( y, -1126, 1 ); + + x = 2067; + y = -2398; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322467, 1 ); + QGSCOMPARENEAR( y, -30635, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 2067, 1 ); + QGSCOMPARENEAR( y, -2398, 1 ); + + x = 1102; + y = -2209; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321777, 1 ); + QGSCOMPARENEAR( y, -12647, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 1102, 1 ); + QGSCOMPARENEAR( y, -2209, 1 ); + } + + void testProjectiveExact() + { + QgsProjectiveGeorefTransform transform; + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ) ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + + void testProjectiveExactRasterYAxis() + { + QgsProjectiveGeorefTransform transform; + //this is a identity transform! + transform.updateParametersFromGcps( QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), + QVector< QgsPointXY >() << QgsPointXY( 321210, 130280 ) + << QgsPointXY( 322698, 129979 ) + << QgsPointXY( 322501, 192061 ) + << QgsPointXY( 321803, 192198 ), true ); + + double x = 321212; + double y = 123003; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, -123003, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321212, 1 ); + QGSCOMPARENEAR( y, 123003, 1 ); + + x = 322702; + y = 140444; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, -140444, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322702, 1 ); + QGSCOMPARENEAR( y, 140444, 1 ); + + x = 322496; + y = 194554; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, -194554, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 322496, 1 ); + QGSCOMPARENEAR( y, 194554, 1 ); + + x = 321800; + y = 186514; + QVERIFY( transform.transform( x, y ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, -186514, 1 ); + QVERIFY( transform.transform( x, y, true ) ); + QGSCOMPARENEAR( x, 321800, 1 ); + QGSCOMPARENEAR( y, 186514, 1 ); + } + }; QGSTEST_MAIN( TestQgsGcpTransformer )