[QGIS Commit] r13987 - trunk/qgis/src/plugins/georeferencer
svn_qgis at osgeo.org
svn_qgis at osgeo.org
Fri Jul 30 10:05:55 EDT 2010
Author: mmassing
Date: 2010-07-30 14:05:55 +0000 (Fri, 30 Jul 2010)
New Revision: 13987
Modified:
trunk/qgis/src/plugins/georeferencer/qgsgeorefplugingui.cpp
trunk/qgis/src/plugins/georeferencer/qgsgeoreftransform.cpp
trunk/qgis/src/plugins/georeferencer/qgsgeoreftransform.h
trunk/qgis/src/plugins/georeferencer/qgsleastsquares.cpp
trunk/qgis/src/plugins/georeferencer/qgsleastsquares.h
trunk/qgis/src/plugins/georeferencer/qgstransformsettingsdialog.cpp
Log:
Add support for projective transformations, which can relate
perspective images of a plane.
Fits a homography to the GCP using total least squares, which
minimizes algebraic error.
Modified: trunk/qgis/src/plugins/georeferencer/qgsgeorefplugingui.cpp
===================================================================
--- trunk/qgis/src/plugins/georeferencer/qgsgeorefplugingui.cpp 2010-07-30 13:24:13 UTC (rev 13986)
+++ trunk/qgis/src/plugins/georeferencer/qgsgeorefplugingui.cpp 2010-07-30 14:05:55 UTC (rev 13987)
@@ -1796,6 +1796,8 @@
return tr( "Polynomial 3" );
case QgsGeorefTransform::ThinPlateSpline:
return tr( "Thin plate spline (TPS)" );
+ case QgsGeorefTransform::Projective:
+ return tr( "Projective" );
default:
return tr( "Not set" );
}
Modified: trunk/qgis/src/plugins/georeferencer/qgsgeoreftransform.cpp
===================================================================
--- trunk/qgis/src/plugins/georeferencer/qgsgeoreftransform.cpp 2010-07-30 13:24:13 UTC (rev 13986)
+++ trunk/qgis/src/plugins/georeferencer/qgsgeoreftransform.cpp 2010-07-30 14:05:55 UTC (rev 13987)
@@ -108,6 +108,35 @@
void *mGDALTransformerArgs;
};
+/**
+ * A planar projective transform, expressed by a homography.
+ *
+ * Implements model fitting which minimizes algebraic error using total least squares.
+ */
+class QgsProjectiveGeorefTransform : public QgsGeorefTransformInterface
+{
+ public:
+ QgsProjectiveGeorefTransform() { }
+ ~QgsProjectiveGeorefTransform() { }
+
+ bool updateParametersFromGCPs( const std::vector<QgsPoint> &mapCoords, const std::vector<QgsPoint> &pixelCoords );
+ uint getMinimumGCPCount() const;
+
+ GDALTransformerFunc GDALTransformer() const { return QgsProjectiveGeorefTransform::projective_transform; }
+ void *GDALTransformerArgs() const { return ( void * )&mParameters; }
+ private:
+ struct ProjectiveParameters
+ {
+ double H[9]; // Homography
+ double Hinv[9]; // Inverted homography
+ bool hasInverse; // Could the inverted homography be calculated?
+ } mParameters;
+
+ static int projective_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
+ double *x, double *y, double *z, int *panSuccess );
+};
+
+
QgsGeorefTransform::QgsGeorefTransform( const QgsGeorefTransform &other )
{
mTransformParametrisation = InvalidTransform;
@@ -216,6 +245,7 @@
case PolynomialOrder2: return new QgsGDALGeorefTransform( false, 2 );
case PolynomialOrder3: return new QgsGDALGeorefTransform( false, 3 );
case ThinPlateSpline: return new QgsGDALGeorefTransform( true, 0 );
+ case Projective: return new QgsProjectiveGeorefTransform;
default: return NULL;
break;
}
@@ -527,3 +557,109 @@
GDALDestroyGCPTransformer( mGDALTransformerArgs );
}
}
+
+bool QgsProjectiveGeorefTransform::updateParametersFromGCPs( const std::vector<QgsPoint> &mapCoords, const std::vector<QgsPoint> &pixelCoords )
+{
+ if ( mapCoords.size() < getMinimumGCPCount() )
+ return false;
+
+ QgsLeastSquares::projective( mapCoords, pixelCoords, mParameters.H );
+
+ // Invert the homography matrix using adjoint matrix
+ double *H = mParameters.H;
+
+ double adjoint[9];
+ adjoint[0] = H[4]*H[8] - H[5]*H[7];
+ adjoint[1] = -H[1]*H[8] + H[2]*H[7];
+ adjoint[2] = H[1]*H[5] - H[2]*H[4];
+
+ adjoint[3] = -H[3]*H[8] + H[5]*H[6];
+ adjoint[4] = H[0]*H[8] - H[2]*H[6];
+ adjoint[5] = -H[0]*H[5] + H[2]*H[3];
+
+ adjoint[6] = H[3]*H[7] - H[4]*H[6];
+ adjoint[7] = -H[0]*H[7] + H[1]*H[6];
+ adjoint[8] = H[0]*H[4] - H[1]*H[3];
+
+ double det = H[0]*adjoint[0] + H[3]*adjoint[1] + H[6]*adjoint[2];
+
+ if (std::abs(det) < 1024.0*std::numeric_limits<double>::epsilon())
+ {
+ mParameters.hasInverse = false;
+ }
+ else
+ {
+ mParameters.hasInverse = true;
+ double oo_det = 1.0/det;
+ for ( int i = 0; i < 9; i++ )
+ {
+ mParameters.Hinv[i] = adjoint[i]*oo_det;
+ }
+
+ for (int i = 0; i < 3; i++) {
+ for (int j = 0; j < 3; j++) {
+ double sum = 0.0;
+ for (int k = 0; k < 3; k++) {
+ sum+= H[j*3 + k]*mParameters.Hinv[i + k*3];
+ }
+ std::cout<<sum<<", ";
+ }
+ std::cout<<std::endl;
+ }
+ }
+ return true;
+}
+
+uint QgsProjectiveGeorefTransform::getMinimumGCPCount() const
+{
+ return 4;
+}
+
+int QgsProjectiveGeorefTransform::projective_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
+ double *x, double *y, double *z, int *panSuccess )
+{
+ ProjectiveParameters* t = static_cast<ProjectiveParameters*>( pTransformerArg );
+ if ( t == NULL )
+ {
+ return false;
+ }
+
+ double *H;
+ if ( !bDstToSrc )
+ {
+ H = t->H;
+ }
+ else
+ {
+ if ( !t->hasInverse )
+ {
+ for (int i = 0; i < nPointCount; ++i )
+ {
+ panSuccess[i] = false;
+ }
+ return false;
+ }
+ H = t->Hinv;
+ }
+
+
+ for (int i = 0; i < nPointCount; ++i )
+ {
+ double Z = x[i]*H[6] + y[i]*H[7] + H[8];
+ // Projects to infinity?
+ if ( std::abs(Z) < 1024.0*std::numeric_limits<double>::epsilon() )
+ {
+ panSuccess[i] = false;
+ continue;
+ }
+ double X = (x[i]*H[0] + y[i]*H[1] + H[2]) / Z;
+ double Y = (x[i]*H[3] + y[i]*H[4] + H[5]) / Z;
+
+ x[i] = X;
+ y[i] = Y;
+
+ panSuccess[i] = true;
+ }
+
+ return true;
+}
Modified: trunk/qgis/src/plugins/georeferencer/qgsgeoreftransform.h
===================================================================
--- trunk/qgis/src/plugins/georeferencer/qgsgeoreftransform.h 2010-07-30 13:24:13 UTC (rev 13986)
+++ trunk/qgis/src/plugins/georeferencer/qgsgeoreftransform.h 2010-07-30 14:05:55 UTC (rev 13987)
@@ -67,7 +67,8 @@
PolynomialOrder2,
PolynomialOrder3,
ThinPlateSpline,
- InvalidTransform
+ Projective,
+ InvalidTransform = 65535
};
QgsGeorefTransform( TransformParametrisation parametrisation );
Modified: trunk/qgis/src/plugins/georeferencer/qgsleastsquares.cpp
===================================================================
--- trunk/qgis/src/plugins/georeferencer/qgsleastsquares.cpp 2010-07-30 13:24:13 UTC (rev 13986)
+++ trunk/qgis/src/plugins/georeferencer/qgsleastsquares.cpp 2010-07-30 14:05:55 UTC (rev 13987)
@@ -173,3 +173,87 @@
}
+
+// Fits a homography to the given corresponding points, and
+// return it in H (row-major format).
+void QgsLeastSquares::projective( std::vector<QgsPoint> mapCoords,
+ std::vector<QgsPoint> pixelCoords,
+ double H[9] )
+{
+ if ( mapCoords.size() < 4 )
+ {
+ throw std::domain_error( QObject::tr( "Fitting a projective transform requires at least 4 corresponding points." ).toLocal8Bit().constData() );
+ }
+
+ // 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 n = 9;
+ gsl_matrix *S = gsl_matrix_alloc ( m, n );
+
+ for ( uint i = 0; i < mapCoords.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, 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+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, 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 );
+ }
+
+ if (mapCoords.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)
+ // to work around this, because the solution lies in the right nullspace of S, and
+ // gsl only supports a thin SVD of S^T, which does not return these vectors.
+
+ // HACK: duplicate last row to get a 9x9 equation system
+ for (int j = 0; j < 9; j++)
+ {
+ gsl_matrix_set( S, 8, j, gsl_matrix_get( S, 7, j) );
+ }
+ }
+
+ // Solve Sh = 0 in the total least squares sense, i.e.
+ // with Sh = min and |h|=1. The solution "h" is given by the
+ // right singular eigenvector of S corresponding, to the smallest
+ // singular value (via SVD).
+ gsl_matrix *V = gsl_matrix_alloc (n, n);
+ gsl_vector *singular_values = gsl_vector_alloc(n);
+ gsl_vector *work = gsl_vector_alloc(n);
+
+ // V = n x n
+ // U = m x n (thin SVD) U D V^T
+ gsl_linalg_SV_decomp(S, V, singular_values, work);
+
+ double eigen[9];
+ // Columns of V store the right singular vectors of S
+ for (int i = 0; i < n; i++)
+ {
+ H[i] = gsl_matrix_get( V, i, n - 1 );
+ }
+
+ gsl_matrix_free( S );
+ gsl_matrix_free( V );
+ gsl_vector_free( singular_values );
+ gsl_vector_free( work );
+}
+
+
Modified: trunk/qgis/src/plugins/georeferencer/qgsleastsquares.h
===================================================================
--- trunk/qgis/src/plugins/georeferencer/qgsleastsquares.h 2010-07-30 13:24:13 UTC (rev 13986)
+++ trunk/qgis/src/plugins/georeferencer/qgsleastsquares.h 2010-07-30 14:05:55 UTC (rev 13987)
@@ -36,6 +36,10 @@
static void affine( std::vector<QgsPoint> mapCoords,
std::vector<QgsPoint> pixelCoords );
+
+ static void projective( std::vector<QgsPoint> mapCoords,
+ std::vector<QgsPoint> pixelCoords,
+ double H[9] );
};
Modified: trunk/qgis/src/plugins/georeferencer/qgstransformsettingsdialog.cpp
===================================================================
--- trunk/qgis/src/plugins/georeferencer/qgstransformsettingsdialog.cpp 2010-07-30 13:24:13 UTC (rev 13986)
+++ trunk/qgis/src/plugins/georeferencer/qgstransformsettingsdialog.cpp 2010-07-30 14:05:55 UTC (rev 13987)
@@ -36,6 +36,7 @@
cmbTransformType->addItem( tr( "Linear" ) , ( int )QgsGeorefTransform::Linear ) ;
cmbTransformType->addItem( tr( "Helmert" ), ( int )QgsGeorefTransform::Helmert );
+ cmbTransformType->addItem( tr( "Projective" ), ( int )QgsGeorefTransform::Projective );
cmbTransformType->addItem( tr( "Polynomial 1" ), ( int )QgsGeorefTransform::PolynomialOrder1 );
cmbTransformType->addItem( tr( "Polynomial 2" ), ( int )QgsGeorefTransform::PolynomialOrder2 );
cmbTransformType->addItem( tr( "Polynomial 3" ), ( int )QgsGeorefTransform::PolynomialOrder3 );
@@ -96,7 +97,7 @@
if ( cmbTransformType->currentIndex() == -1 )
tp = QgsGeorefTransform::InvalidTransform;
else
- tp = ( QgsGeorefTransform::TransformParametrisation )cmbTransformType->currentIndex();
+ tp = ( QgsGeorefTransform::TransformParametrisation )cmbTransformType->itemData( cmbTransformType->currentIndex() ).toInt();
rm = ( QgsImageWarper::ResamplingMethod )cmbResampling->currentIndex();
comprMethod = mListCompression.at( cmbCompressionComboBox->currentIndex() );
More information about the QGIS-commit
mailing list