Skip to content

Commit 0532c6d

Browse files
author
mmassing
committedAug 1, 2010
Georeferencer: Normalize coordinates in order to improve numerical conditioning of homography estimation.
git-svn-id: http://svn.osgeo.org/qgis/trunk/qgis@13990 c8812cc2-4d05-0410-92ff-de0c093fc19c
·
final-1_7_1final-1_7_0
1 parent 57e207e commit 0532c6d

File tree

2 files changed

+86
-10
lines changed

2 files changed

+86
-10
lines changed
 

‎src/plugins/georeferencer/qgsgeoreftransform.cpp

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -562,8 +562,16 @@ bool QgsProjectiveGeorefTransform::updateParametersFromGCPs( const std::vector<Q
562562
{
563563
if ( mapCoords.size() < getMinimumGCPCount() )
564564
return false;
565+
566+
567+
// HACK: flip y coordinates, because georeferencer and gdal use different conventions
568+
std::vector<QgsPoint> flippedPixelCoords( pixelCoords.size() );
569+
for ( uint i = 0; i < pixelCoords.size(); i++ )
570+
{
571+
flippedPixelCoords[i] = QgsPoint( pixelCoords[i].x(), -pixelCoords[i].y() );
572+
}
565573

566-
QgsLeastSquares::projective( mapCoords, pixelCoords, mParameters.H );
574+
QgsLeastSquares::projective( mapCoords, flippedPixelCoords, mParameters.H );
567575

568576
// Invert the homography matrix using adjoint matrix
569577
double *H = mParameters.H;

‎src/plugins/georeferencer/qgsleastsquares.cpp

Lines changed: 77 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#include <stdexcept>
1818

1919
#include <gsl/gsl_linalg.h>
20+
#include <gsl/gsl_blas.h>
2021

2122
#include <QObject>
2223

@@ -174,17 +175,74 @@ void QgsLeastSquares::affine( std::vector<QgsPoint> mapCoords,
174175
}
175176

176177

178+
/**
179+
* Scales the given coordinates so that the center of gravity is at the origin and the mean distance to the origin is sqrt(2).
180+
*
181+
* Also returns 3x3 homogenous matrices which can be used to normalize and de-normalize coordinates.
182+
*/
183+
void normalizeCoordinates( const std::vector<QgsPoint> coords, std::vector<QgsPoint> &normalizedCoords,
184+
double normalizeMatrix[9], double denormalizeMatrix[9])
185+
{
186+
// Calculate center of gravity
187+
double cogX = 0.0, cogY = 0.0;
188+
for ( uint i = 0; i < coords.size(); i++ )
189+
{
190+
cogX+= coords[i].x();
191+
cogY+= coords[i].y();
192+
}
193+
cogX*= 1.0/coords.size();
194+
cogY*= 1.0/coords.size();
195+
196+
// Calculate mean distance to origin
197+
double meanDist = 0.0;
198+
for ( uint i = 0; i < coords.size(); i++ )
199+
{
200+
double X = (coords[i].x() - cogX);
201+
double Y = (coords[i].y() - cogY);
202+
meanDist+= sqrt( X*X + Y*Y );
203+
}
204+
meanDist*= 1.0/coords.size();
205+
206+
double OOD = meanDist/sqrt(2);
207+
double D = 1.0/OOD;
208+
normalizedCoords.resize(coords.size());
209+
for ( uint i = 0; i < coords.size(); i++ )
210+
{
211+
normalizedCoords[i] = QgsPoint( (coords[i].x() - cogX)*D, (coords[i].y() - cogY)*D );
212+
}
213+
214+
normalizeMatrix[0] = D; normalizeMatrix[1] = 0.0; normalizeMatrix[2] = -cogX*D;
215+
normalizeMatrix[3] = 0.0; normalizeMatrix[4] = D; normalizeMatrix[5] = -cogY*D;
216+
normalizeMatrix[6] = 0.0; normalizeMatrix[7] = 0.0; normalizeMatrix[8] = 1.0;
217+
218+
denormalizeMatrix[0] = OOD; denormalizeMatrix[1] = 0.0; denormalizeMatrix[2] = cogX;
219+
denormalizeMatrix[3] = 0.0; denormalizeMatrix[4] = OOD; denormalizeMatrix[5] = cogY;
220+
denormalizeMatrix[6] = 0.0; denormalizeMatrix[7] = 0.0; denormalizeMatrix[8] = 1.0;
221+
}
222+
177223
// Fits a homography to the given corresponding points, and
178224
// return it in H (row-major format).
179225
void QgsLeastSquares::projective( std::vector<QgsPoint> mapCoords,
180226
std::vector<QgsPoint> pixelCoords,
181227
double H[9] )
182228
{
229+
Q_ASSERT( mapCoords.size() == pixelCoords.size() );
230+
183231
if ( mapCoords.size() < 4 )
184232
{
185233
throw std::domain_error( QObject::tr( "Fitting a projective transform requires at least 4 corresponding points." ).toLocal8Bit().constData() );
186234
}
187235

236+
std::vector<QgsPoint> mapCoordsNormalized;
237+
std::vector<QgsPoint> pixelCoordsNormalized;
238+
239+
double normMap[9], denormMap[9];
240+
double normPixel[9], denormPixel[9];
241+
normalizeCoordinates( mapCoords, mapCoordsNormalized, normMap, denormMap );
242+
normalizeCoordinates( pixelCoords, pixelCoordsNormalized, normPixel, denormPixel );
243+
mapCoords = mapCoordsNormalized;
244+
pixelCoords = pixelCoordsNormalized;
245+
188246
// GSL does not support a full SVD, so we artificially add a linear dependent row
189247
// to the matrix in case the system is underconstrained.
190248
uint m = std::max( 9u, ( uint )mapCoords.size() * 2u );
@@ -193,28 +251,28 @@ void QgsLeastSquares::projective( std::vector<QgsPoint> mapCoords,
193251

194252
for ( uint i = 0; i < mapCoords.size(); i++ )
195253
{
196-
gsl_matrix_set( S, i*2, 0, pixelCoords[i].x() );
197-
gsl_matrix_set( S, i*2, 1, -pixelCoords[i].y() );
254+
gsl_matrix_set( S, i*2, 0, pixelCoords[i].x() );
255+
gsl_matrix_set( S, i*2, 1, pixelCoords[i].y() );
198256
gsl_matrix_set( S, i*2, 2, 1.0 );
199257

200258
gsl_matrix_set( S, i*2, 3, 0.0 );
201259
gsl_matrix_set( S, i*2, 4, 0.0 );
202260
gsl_matrix_set( S, i*2, 5, 0.0 );
203261

204-
gsl_matrix_set( S, i*2, 6, -mapCoords[i].x()* pixelCoords[i].x() );
205-
gsl_matrix_set( S, i*2, 7, -mapCoords[i].x()* -pixelCoords[i].y() );
262+
gsl_matrix_set( S, i*2, 6, -mapCoords[i].x()*pixelCoords[i].x() );
263+
gsl_matrix_set( S, i*2, 7, -mapCoords[i].x()*pixelCoords[i].y() );
206264
gsl_matrix_set( S, i*2, 8, -mapCoords[i].x()*1.0 );
207265

208266
gsl_matrix_set( S, i*2 + 1, 0, 0.0 );
209267
gsl_matrix_set( S, i*2 + 1, 1, 0.0 );
210268
gsl_matrix_set( S, i*2 + 1, 2, 0.0 );
211269

212270
gsl_matrix_set( S, i*2 + 1, 3, pixelCoords[i].x() );
213-
gsl_matrix_set( S, i*2 + 1, 4, -pixelCoords[i].y() );
271+
gsl_matrix_set( S, i*2 + 1, 4, pixelCoords[i].y() );
214272
gsl_matrix_set( S, i*2 + 1, 5, 1.0 );
215273

216-
gsl_matrix_set( S, i*2 + 1, 6, -mapCoords[i].y()* pixelCoords[i].x() );
217-
gsl_matrix_set( S, i*2 + 1, 7, -mapCoords[i].y()* -pixelCoords[i].y() );
274+
gsl_matrix_set( S, i*2 + 1, 6, -mapCoords[i].y()*pixelCoords[i].x() );
275+
gsl_matrix_set( S, i*2 + 1, 7, -mapCoords[i].y()*pixelCoords[i].y() );
218276
gsl_matrix_set( S, i*2 + 1, 8, -mapCoords[i].y()*1.0 );
219277
}
220278

@@ -250,10 +308,20 @@ void QgsLeastSquares::projective( std::vector<QgsPoint> mapCoords,
250308
H[i] = gsl_matrix_get( V, i, n - 1 );
251309
}
252310

311+
gsl_matrix *prodMatrix = gsl_matrix_alloc( 3, 3 );
312+
313+
gsl_matrix_view Hmatrix = gsl_matrix_view_array( H, 3, 3 );
314+
gsl_matrix_view normPixelMatrix = gsl_matrix_view_array( normPixel, 3, 3 );
315+
gsl_matrix_view denormMapMatrix = gsl_matrix_view_array( denormMap, 3, 3 );
316+
317+
// Change coordinate frame of image and pre-image from normalized to map and pixel coordinates.
318+
// H' = denormalizeMapCoords*H*normalizePixelCoords
319+
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &Hmatrix.matrix, &normPixelMatrix.matrix, 0.0, prodMatrix );
320+
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, &denormMapMatrix.matrix, prodMatrix, 0.0, &Hmatrix.matrix );
321+
322+
gsl_matrix_free( prodMatrix );
253323
gsl_matrix_free( S );
254324
gsl_matrix_free( V );
255325
gsl_vector_free( singular_values );
256326
gsl_vector_free( work );
257327
}
258-
259-

0 commit comments

Comments
 (0)
Please sign in to comment.