Skip to content

Commit 8b0c4b4

Browse files
author
mmassing
committedMay 23, 2010
Applied patch #2731 by mhugent: show georeferencer residuals in map units when possible.
git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@13552 c8812cc2-4d05-0410-92ff-de0c093fc19c
·
final-1_7_1final-1_7_0
1 parent 738b4ca commit 8b0c4b4

File tree

3 files changed

+177
-92
lines changed

3 files changed

+177
-92
lines changed
 

‎src/plugins/georeferencer/qgsgcplistmodel.cpp

Lines changed: 36 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616

1717
#include "qgsgcplist.h"
1818
#include "qgsgcplistmodel.h"
19-
19+
#include "qgis.h"
2020
#include "qgsgeorefdatapoint.h"
2121
#include "qgsgeoreftransform.h"
2222

@@ -85,23 +85,43 @@ void QgsGCPListModel::updateModel()
8585
if ( !mGCPList )
8686
return;
8787

88-
// // Setup table header
89-
QStringList itemLabels;
90-
// // Set column headers
91-
itemLabels << "on/off" << "id" << "srcX" << "srcY" << "dstX" << "dstY" << "dX" << "dY" << "residual";
92-
// setColumnCount(itemLabels.size());
93-
setHorizontalHeaderLabels( itemLabels );
94-
setRowCount( mGCPList->size() );
95-
9688
bool bTransformUpdated = false;
89+
bool wldTransform = false;
90+
double wldScaleX, wldScaleY, rotation;
91+
QgsPoint origin;
92+
9793
if ( mGeorefTransform )
9894
{
9995
vector<QgsPoint> mapCoords, pixelCoords;
10096
mGCPList->createGCPVectors( mapCoords, pixelCoords );
10197

10298
// TODO: the parameters should probable be updated externally (by user interaction)
10399
bTransformUpdated = mGeorefTransform->updateParametersFromGCPs( mapCoords, pixelCoords );
100+
//transformation that involves only scaling and rotation (linear or helmert) ?
101+
wldTransform = mGeorefTransform->getOriginScaleRotation( origin, wldScaleX, wldScaleY, rotation );
102+
if ( wldTransform && !doubleNear( rotation, 0.0 ) )
103+
{
104+
wldScaleX *= cos( rotation );
105+
wldScaleY *= cos( rotation );
106+
}
107+
if ( wldTransform )
108+
{
109+
110+
}
111+
}
112+
113+
// // Setup table header
114+
QStringList itemLabels;
115+
if ( wldTransform )
116+
{
117+
itemLabels << "on/off" << "id" << "srcX" << "srcY" << "dstX" << "dstY" << QString( "dX[" ) + tr( "map units" ) + "]" << QString( "dY[" ) + tr( "map units" ) + "]" << "residual";
104118
}
119+
else
120+
{
121+
itemLabels << "on/off" << "id" << "srcX" << "srcY" << "dstX" << "dstY" << QString( "dX[" ) + tr( "pixels" ) + "]" << QString( "dY[" ) + tr( "pixels" ) + "]" << "residual";
122+
}
123+
setHorizontalHeaderLabels( itemLabels );
124+
setRowCount( mGCPList->size() );
105125

106126
for ( int i = 0; i < mGCPList->sizeAll(); ++i )
107127
{
@@ -135,8 +155,13 @@ void QgsGCPListModel::updateModel()
135155
// As transforms of order >=2 are not invertible, we are only
136156
// interested in the residual in this direction
137157
mGeorefTransform->transformWorldToRaster( p->mapCoords(), dst );
138-
dX = ( dst.x() - p->pixelCoords().x() );
158+
dX = ( dst.x() - p->pixelCoords().x() );
139159
dY = -( dst.y() - p->pixelCoords().y() );
160+
if ( wldTransform )
161+
{
162+
dX *= wldScaleX;
163+
dY *= wldScaleY;
164+
}
140165
residual = sqrt( dX * dX + dY * dY );
141166
}
142167
else
@@ -152,7 +177,7 @@ void QgsGCPListModel::updateModel()
152177
if ( residual >= 0.f )
153178
{
154179
setItem( i, j++, QGSSTANDARDITEM( dX ) /*create_item<double>(dX)*/ );
155-
setItem( i, j++, QGSSTANDARDITEM( dY ) /*create_item<double>(-dY)*/);
180+
setItem( i, j++, QGSSTANDARDITEM( dY ) /*create_item<double>(-dY)*/ );
156181
setItem( i, j++, QGSSTANDARDITEM( residual ) /*create_item<double>(residual)*/ );
157182
}
158183
else

‎src/plugins/georeferencer/qgsgeorefplugingui.cpp

Lines changed: 128 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -136,7 +136,6 @@ void QgsGeorefPluginGui::dockThisWindow( bool dock )
136136

137137
QgsGeorefPluginGui::~QgsGeorefPluginGui()
138138
{
139-
QgsTransformSettingsDialog::resetSettings();
140139
clearGCPData();
141140

142141
// delete layer (and don't signal it as it's our private layer)
@@ -336,6 +335,7 @@ bool QgsGeorefPluginGui::getTransformSettings()
336335
mActionLinkQGisToGeoref->setEnabled( false );
337336
}
338337

338+
updateTransformParamLabel();
339339
return true;
340340
}
341341

@@ -472,6 +472,7 @@ void QgsGeorefPluginGui::addPoint( const QgsPoint& pixelCoords, const QgsPoint&
472472
}
473473

474474
connect( mCanvas, SIGNAL( extentsChanged() ), pnt, SLOT( updateCoords() ) );
475+
updateGeorefTransform();
475476

476477
// if (verbose)
477478
// logRequaredGCPs();
@@ -498,13 +499,15 @@ void QgsGeorefPluginGui::deleteDataPoint( const QPoint &coords )
498499
break;
499500
}
500501
}
502+
updateGeorefTransform();
501503
}
502504

503505
void QgsGeorefPluginGui::deleteDataPoint( int index )
504506
{
505507
mGCPListWidget->model()->removeRow( index );
506508
delete mPoints.takeAt( index );
507509
mGCPListWidget->updateGCPList();
510+
updateGeorefTransform();
508511
}
509512

510513
void QgsGeorefPluginGui::selectPoint( const QPoint &p )
@@ -930,12 +933,24 @@ void QgsGeorefPluginGui::createDockWidgets()
930933
this, SLOT( replaceDataPoint( QgsGeorefDataPoint*, int ) ) );
931934
connect( mGCPListWidget, SIGNAL( deleteDataPoint( int ) ),
932935
this, SLOT( deleteDataPoint( int ) ) );
936+
connect( mGCPListWidget, SIGNAL( pointEnabled( QgsGeorefDataPoint*, int ) ), this, SLOT( updateGeorefTransform() ) );
933937
}
934938

935939
void QgsGeorefPluginGui::createStatusBar()
936940
{
937941
QFont myFont( "Arial", 9 );
938942

943+
mTransformParamLabel = new QLabel( statusBar() );
944+
mTransformParamLabel->setFont( myFont );
945+
mTransformParamLabel->setMinimumWidth( 10 );
946+
mTransformParamLabel->setMaximumHeight( 20 );
947+
mTransformParamLabel->setMargin( 3 );
948+
mTransformParamLabel->setAlignment( Qt::AlignCenter );
949+
mTransformParamLabel->setFrameStyle( QFrame::NoFrame );
950+
mTransformParamLabel->setText( tr( "Transform: " ) + convertTransformEnumToString( mTransformParam ) );
951+
mTransformParamLabel->setToolTip( tr( "Current transform parametrisation" ) );
952+
statusBar()->addPermanentWidget( mTransformParamLabel, 0 );
953+
939954
mCoordsLabel = new QLabel( QString(), statusBar() );
940955
mCoordsLabel->setFont( myFont );
941956
mCoordsLabel->setMinimumWidth( 10 );
@@ -947,17 +962,6 @@ void QgsGeorefPluginGui::createStatusBar()
947962
mCoordsLabel->setText( tr( "Coordinate: " ) );
948963
mCoordsLabel->setToolTip( tr( "Current map coordinate" ) );
949964
statusBar()->addPermanentWidget( mCoordsLabel, 0 );
950-
951-
mTransformParamLabel = new QLabel( statusBar() );
952-
mTransformParamLabel->setFont( myFont );
953-
mTransformParamLabel->setMinimumWidth( 10 );
954-
mTransformParamLabel->setMaximumHeight( 20 );
955-
mTransformParamLabel->setMargin( 3 );
956-
mTransformParamLabel->setAlignment( Qt::AlignCenter );
957-
mTransformParamLabel->setFrameStyle( QFrame::NoFrame );
958-
mTransformParamLabel->setText( tr( "Transform: " ) + convertTransformEnumToString( mTransformParam ) );
959-
mTransformParamLabel->setToolTip( tr( "Current transform parametrisation" ) );
960-
statusBar()->addPermanentWidget( mTransformParamLabel, 0 );
961965
}
962966

963967
void QgsGeorefPluginGui::setupConnections()
@@ -1208,6 +1212,53 @@ bool QgsGeorefPluginGui::writeWorldFile( QgsPoint origin, double pixelXSize, dou
12081212
return true;
12091213
}
12101214

1215+
bool QgsGeorefPluginGui::calculateMeanError( double& error ) const
1216+
{
1217+
if ( mGeorefTransform.transformParametrisation() == QgsGeorefTransform::InvalidTransform )
1218+
{
1219+
return false;
1220+
}
1221+
1222+
int nPointsEnabled = 0;
1223+
QgsGCPList::const_iterator gcpIt = mPoints.constBegin();
1224+
for ( ; gcpIt != mPoints.constEnd(); ++gcpIt )
1225+
{
1226+
if (( *gcpIt )->isEnabled() )
1227+
{
1228+
++nPointsEnabled;
1229+
}
1230+
}
1231+
1232+
if ( nPointsEnabled == mGeorefTransform.getMinimumGCPCount() )
1233+
{
1234+
error = 0;
1235+
return true;
1236+
}
1237+
else if ( nPointsEnabled < mGeorefTransform.getMinimumGCPCount() )
1238+
{
1239+
return false;
1240+
}
1241+
1242+
double sumVxSquare = 0;
1243+
double sumVySquare = 0;
1244+
double resXMap, resYMap;
1245+
1246+
gcpIt = mPoints.constBegin();
1247+
for ( ; gcpIt != mPoints.constEnd(); ++gcpIt )
1248+
{
1249+
if (( *gcpIt )->isEnabled() )
1250+
{
1251+
sumVxSquare += (( *gcpIt )->residual().x() * ( *gcpIt )->residual().x() );
1252+
sumVySquare += (( *gcpIt )->residual().y() * ( *gcpIt )->residual().y() );
1253+
}
1254+
}
1255+
1256+
// Calculate the root mean square error, adjusted for degrees of freedom of the transform
1257+
// Caveat: The number of DoFs is assumed to be even (as each control point fixes two degrees of freedom).
1258+
error = sqrt(( sumVxSquare + sumVySquare ) / ( nPointsEnabled - mGeorefTransform.getMinimumGCPCount() ));
1259+
return true;
1260+
}
1261+
12111262
bool QgsGeorefPluginGui::writePDFReportFile( const QString& fileName, const QgsGeorefTransform& transform )
12121263
{
12131264
if ( !mCanvas )
@@ -1250,7 +1301,24 @@ bool QgsGeorefPluginGui::writePDFReportFile( const QString& fileName, const QgsG
12501301
titleLabel->setFrame( false );
12511302

12521303
//composer map
1253-
QgsComposerMap* composerMap = new QgsComposerMap( composition, 2, titleLabel->rect().bottom() + titleLabel->transform().dy(), 206, 277 );
1304+
QgsRectangle canvasExtent = mCanvas->extent();
1305+
//calculate width and height considering extent aspect ratio and max Width 206, maxHeight 70
1306+
double widthExtentRatio = 206 / canvasExtent.width();
1307+
double heightExtentRatio = 70 / canvasExtent.height();
1308+
double mapWidthMM = 0;
1309+
double mapHeightMM = 0;
1310+
if ( widthExtentRatio < heightExtentRatio )
1311+
{
1312+
mapWidthMM = 206;
1313+
mapHeightMM = 206 / canvasExtent.width() * canvasExtent.height();
1314+
}
1315+
else
1316+
{
1317+
mapHeightMM = 70;
1318+
mapWidthMM = 70 / canvasExtent.height() * canvasExtent.width();
1319+
}
1320+
1321+
QgsComposerMap* composerMap = new QgsComposerMap( composition, 2, titleLabel->rect().bottom() + titleLabel->transform().dy(), mapWidthMM, mapHeightMM );
12541322
composerMap->setLayerSet( canvasRenderer->layerSet() );
12551323
composerMap->setNewExtent( mCanvas->extent() );
12561324
composerMap->setMapCanvas( mCanvas );
@@ -1266,26 +1334,9 @@ bool QgsGeorefPluginGui::writePDFReportFile( const QString& fileName, const QgsG
12661334
//transformation that involves only scaling and rotation (linear or helmert) ?
12671335
bool wldTransform = transform.getOriginScaleRotation( origin, scaleX, scaleY, rotation );
12681336

1269-
//consider rotation in scale parameter
1270-
double wldScaleX = scaleX;
1271-
double wldScaleY = scaleY;
1272-
if ( wldTransform && !doubleNear( rotation, 0.0 ) )
1273-
{
1274-
wldScaleX *= cos( rotation );
1275-
wldScaleY *= cos( rotation );
1276-
}
1277-
12781337
if ( wldTransform )
12791338
{
1280-
QString parameterTitle = tr( "Transformation parameters" );
1281-
if ( transform.transformParametrisation() == QgsGeorefTransform::Helmert )
1282-
{
1283-
parameterTitle += " (Helmert)";
1284-
}
1285-
else if ( transform.transformParametrisation() == QgsGeorefTransform::Linear )
1286-
{
1287-
parameterTitle += " (Linear)";
1288-
}
1339+
QString parameterTitle = tr( "Transformation parameters" ) + QString( " (" ) + convertTransformEnumToString( transform.transformParametrisation() ) + QString( ")" );
12891340
parameterLabel = new QgsComposerLabel( composition );
12901341
parameterLabel->setFont( titleFont );
12911342
parameterLabel->setText( parameterTitle );
@@ -1304,29 +1355,9 @@ bool QgsGeorefPluginGui::writePDFReportFile( const QString& fileName, const QgsG
13041355
}
13051356
}
13061357

1307-
//calculate mean error (in map units)
1358+
//calculate mean error
13081359
double meanError = 0;
1309-
if ( nPointsEnabled > 2 )
1310-
{
1311-
double sumVxSquare = 0;
1312-
double sumVySquare = 0;
1313-
double resXMap, resYMap;
1314-
1315-
QgsGCPList::const_iterator gcpIt = mPoints.constBegin();
1316-
for ( ; gcpIt != mPoints.constEnd(); ++gcpIt )
1317-
{
1318-
if (( *gcpIt )->isEnabled() )
1319-
{
1320-
resXMap = ( *gcpIt )->residual().x() * wldScaleX;
1321-
resYMap = ( *gcpIt )->residual().y() * wldScaleY;
1322-
sumVxSquare += ( resXMap * resXMap );
1323-
sumVySquare += ( resYMap * resYMap );
1324-
}
1325-
}
1326-
1327-
meanError = sqrt(( sumVxSquare + sumVySquare ) / ( 2 * nPointsEnabled - 4 ) ) * sqrt( 2.0 );
1328-
}
1329-
1360+
calculateMeanError( meanError );
13301361

13311362
parameterTable = new QgsComposerTextTable( composition );
13321363
parameterTable->setHeaderFont( tableHeaderFont );
@@ -1366,7 +1397,7 @@ bool QgsGeorefPluginGui::writePDFReportFile( const QString& fileName, const QgsG
13661397
//convert residual scale bar plot to map units if scaling is equal in x- and y-direction (e.g. helmert)
13671398
if ( wldTransform )
13681399
{
1369-
if ( doubleNear( wldScaleX, wldScaleY ) )
1400+
if ( doubleNear( scaleX, scaleX ) )
13701401
{
13711402
resPlotItem->setPixelToMapUnits( scaleX );
13721403
resPlotItem->setConvertScaleToMapUnits( true );
@@ -1394,17 +1425,7 @@ bool QgsGeorefPluginGui::writePDFReportFile( const QString& fileName, const QgsG
13941425
{
13951426
QStringList currentGCPStrings;
13961427
QPointF residual = ( *gcpIt )->residual();
1397-
double residualX = residual.x();
1398-
if ( wldTransform )
1399-
{
1400-
residualX *= wldScaleX;
1401-
}
1402-
double residualY = residual.y();
1403-
if ( wldTransform )
1404-
{
1405-
residualY *= wldScaleY;
1406-
}
1407-
double residualTot = sqrt( residualX * residualX + residualY * residualY );
1428+
double residualTot = sqrt( residual.x() * residual.x() + residual.y() * residual.y() );
14081429

14091430
currentGCPStrings << QString::number(( *gcpIt )->id() );
14101431
if (( *gcpIt )->isEnabled() )
@@ -1416,7 +1437,7 @@ bool QgsGeorefPluginGui::writePDFReportFile( const QString& fileName, const QgsG
14161437
currentGCPStrings << tr( "no" );
14171438
}
14181439
currentGCPStrings << QString::number(( *gcpIt )->pixelCoords().x(), 'f', 2 ) << QString::number(( *gcpIt )->pixelCoords().y(), 'f', 2 ) << QString::number(( *gcpIt )->mapCoords().x(), 'f', 2 )\
1419-
<< QString::number(( *gcpIt )->mapCoords().y(), 'f', 2 ) << QString::number( residualX ) << QString::number( residualY ) << QString::number( residualTot );
1440+
<< QString::number(( *gcpIt )->mapCoords().y(), 'f', 2 ) << QString::number( residual.x() ) << QString::number( residual.y() ) << QString::number( residualTot );
14201441
gcpTable->addRow( currentGCPStrings );
14211442
}
14221443

@@ -1443,6 +1464,35 @@ bool QgsGeorefPluginGui::writePDFReportFile( const QString& fileName, const QgsG
14431464
return true;
14441465
}
14451466

1467+
void QgsGeorefPluginGui::updateTransformParamLabel()
1468+
{
1469+
if ( !mTransformParamLabel )
1470+
{
1471+
return;
1472+
}
1473+
1474+
QString transformName = convertTransformEnumToString( mGeorefTransform.transformParametrisation() );
1475+
QString labelString = tr( "Transform: " ) + transformName;
1476+
1477+
QgsPoint origin;
1478+
double scaleX, scaleY, rotation;
1479+
if ( mGeorefTransform.getOriginScaleRotation( origin, scaleX, scaleY, rotation ) )
1480+
{
1481+
labelString += " ";
1482+
labelString += tr( "Translation (%1, %2)" ).arg( origin.x() ).arg( origin.y() ); labelString += " ";
1483+
labelString += tr( "Scale (%1, %2)" ).arg( scaleX ).arg( scaleY ); labelString += " ";
1484+
labelString += tr( "Rotation: %1" ).arg( rotation );
1485+
}
1486+
1487+
double meanError = 0;
1488+
if ( calculateMeanError( meanError ) )
1489+
{
1490+
labelString += " ";
1491+
labelString += tr( "Mean error: %1" ).arg( meanError );
1492+
}
1493+
mTransformParamLabel->setText( labelString );
1494+
}
1495+
14461496
// Gdal script
14471497
void QgsGeorefPluginGui::showGDALScript( int argNum... )
14481498
{
@@ -1605,22 +1655,20 @@ bool QgsGeorefPluginGui::checkReadyGeoref()
16051655

16061656
bool QgsGeorefPluginGui::updateGeorefTransform()
16071657
{
1608-
if ( mGCPsDirty || !mGeorefTransform.parametersInitialized() )
1609-
{
1610-
std::vector<QgsPoint> mapCoords, pixelCoords;
1611-
if ( mGCPListWidget->gcpList() )
1612-
mGCPListWidget->gcpList()->createGCPVectors( mapCoords, pixelCoords );
1613-
else
1614-
return false;
1615-
1616-
// Parametrize the transform with GCPs
1617-
if ( !mGeorefTransform.updateParametersFromGCPs( mapCoords, pixelCoords ) )
1618-
{
1619-
return false;
1620-
}
1658+
std::vector<QgsPoint> mapCoords, pixelCoords;
1659+
if ( mGCPListWidget->gcpList() )
1660+
mGCPListWidget->gcpList()->createGCPVectors( mapCoords, pixelCoords );
1661+
else
1662+
return false;
16211663

1622-
mGCPsDirty = false;
1664+
// Parametrize the transform with GCPs
1665+
if ( !mGeorefTransform.updateParametersFromGCPs( mapCoords, pixelCoords ) )
1666+
{
1667+
return false;
16231668
}
1669+
1670+
mGCPsDirty = false;
1671+
updateTransformParamLabel();
16241672
return true;
16251673
}
16261674

‎src/plugins/georeferencer/qgsgeorefplugingui.h

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -102,6 +102,8 @@ class QgsGeorefPluginGui : public QMainWindow, private Ui::QgsGeorefPluginGuiBas
102102
void showMouseCoords( const QgsPoint pt );
103103
void updateMouseCoordinatePrecision();
104104

105+
bool updateGeorefTransform();
106+
105107
private:
106108
enum SaveGCPs
107109
{
@@ -133,6 +135,7 @@ class QgsGeorefPluginGui : public QMainWindow, private Ui::QgsGeorefPluginGuiBas
133135
bool georeference();
134136
bool writeWorldFile( QgsPoint origin, double pixelXSize, double pixelYSize, double rotation );
135137
bool writePDFReportFile( const QString& fileName, const QgsGeorefTransform& transform );
138+
void updateTransformParamLabel();
136139

137140
// gdal script
138141
void showGDALScript( int argNum... );
@@ -145,7 +148,6 @@ class QgsGeorefPluginGui : public QMainWindow, private Ui::QgsGeorefPluginGuiBas
145148

146149
// utils
147150
bool checkReadyGeoref();
148-
bool updateGeorefTransform();
149151
QgsRectangle transformViewportBoundingBox( const QgsRectangle &canvasExtent, const QgsGeorefTransform &t,
150152
bool rasterToWorld = true, uint numSamples = 4 );
151153
QString convertTransformEnumToString( QgsGeorefTransform::TransformParametrisation transform );
@@ -159,6 +161,16 @@ class QgsGeorefPluginGui : public QMainWindow, private Ui::QgsGeorefPluginGuiBas
159161
void logRequaredGCPs();
160162
void clearGCPData();
161163

164+
/**
165+
* Calculates root mean squared error for the currently active
166+
* ground control points and transform method.
167+
* Note that he RMSE measure is adjusted for the degrees of freedom of the
168+
* used polynomial transform.
169+
* @param error out: the mean error
170+
* @return true in case of success
171+
*/
172+
bool calculateMeanError( double& error ) const;
173+
162174
/**Docks / undocks this window*/
163175
void dockThisWindow( bool dock );
164176

0 commit comments

Comments
 (0)
Please sign in to comment.