diff --git a/src/app/georeferencer/qgsgcplist.cpp b/src/app/georeferencer/qgsgcplist.cpp index 7ff6d159127a..3498ad3b8f7f 100644 --- a/src/app/georeferencer/qgsgcplist.cpp +++ b/src/app/georeferencer/qgsgcplist.cpp @@ -18,6 +18,7 @@ #include "qgscoordinatereferencesystem.h" #include "qgscoordinatetransform.h" #include "qgsproject.h" +#include "qgsgeoreftransform.h" #include "qgsgcplist.h" #include @@ -64,6 +65,62 @@ int QgsGCPList::countEnabledPoints() const return s; } +void QgsGCPList::updateResiduals( QgsGeorefTransform *georefTransform, const QgsCoordinateReferenceSystem &targetCrs, const QgsCoordinateTransformContext &context, QgsUnitTypes::RenderUnit residualUnit ) +{ + bool bTransformUpdated = false; + QVector sourceCoordinates; + QVector destinationCoordinates; + createGCPVectors( sourceCoordinates, destinationCoordinates, targetCrs, context ); + + if ( georefTransform ) + { + bTransformUpdated = georefTransform->updateParametersFromGcps( sourceCoordinates, destinationCoordinates, true ); + } + + // update residuals + for ( int i = 0; i < size(); ++i ) + { + QgsGeorefDataPoint *p = at( i ); + + if ( !p ) + continue; + + p->setId( i ); + + const QgsPointXY transformedDestinationPoint = p->transformedDestinationPoint( targetCrs, QgsProject::instance()->transformContext() ); + + double dX = 0; + double dY = 0; + // Calculate residual if transform is available and up-to-date + if ( georefTransform && bTransformUpdated && georefTransform->parametersInitialized() ) + { + QgsPointXY dst; + const QgsPointXY pixel = georefTransform->toSourcePixel( p->sourcePoint() ); + if ( residualUnit == QgsUnitTypes::RenderPixels ) + { + // Transform from world to raster coordinate: + // This is the transform direction used by the warp operation. + // As transforms of order >=2 are not invertible, we are only + // interested in the residual in this direction + if ( georefTransform->transformWorldToRaster( transformedDestinationPoint, dst ) ) + { + dX = ( dst.x() - pixel.x() ); + dY = -( dst.y() - pixel.y() ); + } + } + else if ( residualUnit == QgsUnitTypes::RenderMapUnits ) + { + if ( georefTransform->transformRasterToWorld( pixel, dst ) ) + { + dX = ( dst.x() - transformedDestinationPoint.x() ); + dY = ( dst.y() - transformedDestinationPoint.y() ); + } + } + } + p->setResidual( QPointF( dX, dY ) ); + } +} + QList QgsGCPList::asPoints() const { QList res; diff --git a/src/app/georeferencer/qgsgcplist.h b/src/app/georeferencer/qgsgcplist.h index 8513dd8282c8..6662a95e97e2 100644 --- a/src/app/georeferencer/qgsgcplist.h +++ b/src/app/georeferencer/qgsgcplist.h @@ -19,12 +19,14 @@ #include #include #include "qgis_app.h" +#include "qgsunittypes.h" class QgsGeorefDataPoint; class QgsGcpPoint; class QgsPointXY; class QgsCoordinateReferenceSystem; class QgsCoordinateTransformContext; +class QgsGeorefTransform; /** * A container for GCP data points. @@ -50,6 +52,18 @@ class APP_EXPORT QgsGCPList : public QList */ int countEnabledPoints() const; + /** + * Updates the stored residual sizes for points in the list. + * + * \param georefTransform transformation to use for residual calculation + * \param targetCrs georeference output CRS + * \param context transform context + * \param residualUnit units for residual calculation. Supported values are QgsUnitTypes::RenderPixels or QgsUnitTypes::RenderMapUnits + */ + void updateResiduals( QgsGeorefTransform *georefTransform, + const QgsCoordinateReferenceSystem &targetCrs, const QgsCoordinateTransformContext &context, + QgsUnitTypes::RenderUnit residualUnit ); + /** * Returns the container as a list of GCP points. */ diff --git a/tests/src/app/testqgsgeoreferencer.cpp b/tests/src/app/testqgsgeoreferencer.cpp index c9a6ad1090de..d42b76770696 100644 --- a/tests/src/app/testqgsgeoreferencer.cpp +++ b/tests/src/app/testqgsgeoreferencer.cpp @@ -49,6 +49,7 @@ class TestQgsGeoreferencer : public QObject void testTransformImageNoGeoference(); void testTransformImageWithExistingGeoreference(); void testRasterChangeCoords(); + void testUpdateResiduals(); private: QgisApp *mQgisApp = nullptr; @@ -483,5 +484,65 @@ void TestQgsGeoreferencer::testRasterChangeCoords() QCOMPARE( transform.toColumnLine( QgsPointXY( 100, 200 ) ).y(), 200.0 ); } +void TestQgsGeoreferencer::testUpdateResiduals() +{ + // test updating residuals + QgsGeorefTransform transform( QgsGcpTransformerInterface::TransformMethod::Linear ); + transform.loadRaster( QStringLiteral( TEST_DATA_DIR ) + QStringLiteral( "/landsat.tif" ) ); + QVERIFY( transform.hasExistingGeoreference() ); + + QgsGCPList list; + QgsMapCanvas c1; + QgsMapCanvas c2; + list.append( new QgsGeorefDataPoint( &c1, &c2, + QgsPointXY( 781662.375, 3350923.125 ), QgsPointXY( -30, 40 ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), + true ) ); + list.append( new QgsGeorefDataPoint( &c1, &c2, + QgsPointXY( 787362.375, 3350923.125 ), QgsPointXY( 16697923, -3503549 ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ), + true ) ); + list.append( new QgsGeorefDataPoint( &c1, &c2, + QgsPointXY( 787362.375, 3362323.125 ), QgsPointXY( -35, 42 ), QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), + true ) ); + + list.updateResiduals( &transform, QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsProject::instance()->transformContext(), QgsUnitTypes::RenderPixels ); + QGSCOMPARENEAR( list.at( 0 )->residual().x(), 0, 0.00001 ); + QGSCOMPARENEAR( list.at( 0 )->residual().y(), -189.189, 0.1 ); + QGSCOMPARENEAR( list.at( 1 )->residual().x(), 105.7142, 0.1 ); + QGSCOMPARENEAR( list.at( 1 )->residual().y(), 189.189, 0.1 ); + QGSCOMPARENEAR( list.at( 2 )->residual().x(), -105.7142, 0.1 ); + QGSCOMPARENEAR( list.at( 2 )->residual().y(), 0, 0.00001 ); + + // in map units + list.updateResiduals( &transform, QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:4326" ) ), QgsProject::instance()->transformContext(), QgsUnitTypes::RenderMapUnits ); + QGSCOMPARENEAR( list.at( 0 )->residual().x(), 0, 0.00001 ); + QGSCOMPARENEAR( list.at( 0 )->residual().y(), -34.999, 0.1 ); + QGSCOMPARENEAR( list.at( 1 )->residual().x(), -92.499, 0.1 ); + QGSCOMPARENEAR( list.at( 1 )->residual().y(), 34.99999, 0.1 ); + QGSCOMPARENEAR( list.at( 2 )->residual().x(), 92.4999972, 0.1 ); + QGSCOMPARENEAR( list.at( 2 )->residual().y(), 0, 0.00001 ); + + // different target CRS + list.updateResiduals( &transform, QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ), QgsProject::instance()->transformContext(), QgsUnitTypes::RenderPixels ); + QGSCOMPARENEAR( list.at( 0 )->residual().x(), 0, 0.00001 ); + QGSCOMPARENEAR( list.at( 0 )->residual().y(), -186.828, 0.1 ); + QGSCOMPARENEAR( list.at( 1 )->residual().x(), 105.7142, 0.1 ); + QGSCOMPARENEAR( list.at( 1 )->residual().y(), 186.828, 0.1 ); + QGSCOMPARENEAR( list.at( 2 )->residual().x(), -105.7142, 0.1 ); + QGSCOMPARENEAR( list.at( 2 )->residual().y(), 0, 0.00001 ); + + // projective transform -- except 0 residuals here + QgsGeorefTransform projective( QgsGcpTransformerInterface::TransformMethod::Projective ); + projective.loadRaster( QStringLiteral( TEST_DATA_DIR ) + QStringLiteral( "/landsat.tif" ) ); + QVERIFY( projective.hasExistingGeoreference() ); + + list.updateResiduals( &projective, QgsCoordinateReferenceSystem( QStringLiteral( "EPSG:3857" ) ), QgsProject::instance()->transformContext(), QgsUnitTypes::RenderPixels ); + QGSCOMPARENEAR( list.at( 0 )->residual().x(), 0, 0.00001 ); + QGSCOMPARENEAR( list.at( 0 )->residual().y(), 0, 0.00001 ); + QGSCOMPARENEAR( list.at( 1 )->residual().x(), 0, 0.00001 ); + QGSCOMPARENEAR( list.at( 1 )->residual().y(), 0, 0.00001 ); + QGSCOMPARENEAR( list.at( 2 )->residual().x(), 0, 0.00001 ); + QGSCOMPARENEAR( list.at( 2 )->residual().y(), 0, 0.00001 ); +} + QGSTEST_MAIN( TestQgsGeoreferencer ) #include "testqgsgeoreferencer.moc"