Skip to content

Commit

Permalink
Move code for calculating residuals to QgsGCPList out of model
Browse files Browse the repository at this point in the history
  • Loading branch information
@nyalldawson
nyalldawson committed Feb 9, 2022
1 parent b174b01 commit e016f02
Show file tree
Hide file tree
Showing 3 changed files with 132 additions and 0 deletions.
57 changes: 57 additions & 0 deletions src/app/georeferencer/qgsgcplist.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
#include "qgscoordinatereferencesystem.h"
#include "qgscoordinatetransform.h"
#include "qgsproject.h"
#include "qgsgeoreftransform.h"

#include "qgsgcplist.h"
#include <QDir>
Expand Down Expand Up @@ -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<QgsPointXY> sourceCoordinates;
QVector<QgsPointXY> 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<QgsGcpPoint> QgsGCPList::asPoints() const
{
QList<QgsGcpPoint> res;
Expand Down
14 changes: 14 additions & 0 deletions src/app/georeferencer/qgsgcplist.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,14 @@
#include <QList>
#include <QVector>
#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.
Expand All @@ -50,6 +52,18 @@ class APP_EXPORT QgsGCPList : public QList<QgsGeorefDataPoint * >
*/
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.
*/
Expand Down
61 changes: 61 additions & 0 deletions tests/src/app/testqgsgeoreferencer.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ class TestQgsGeoreferencer : public QObject
void testTransformImageNoGeoference();
void testTransformImageWithExistingGeoreference();
void testRasterChangeCoords();
void testUpdateResiduals();

private:
QgisApp *mQgisApp = nullptr;
Expand Down Expand Up @@ -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"

0 comments on commit e016f02

Please sign in to comment.