Skip to content
Permalink
Browse files
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.

git-svn-id: http://svn.osgeo.org/qgis/trunk@13987 c8812cc2-4d05-0410-92ff-de0c093fc19c
  • Loading branch information
mmassing committed Jul 30, 2010
1 parent 8cec741 commit ca918e321f108e8bc300a94a4af21ad663485481
@@ -1796,6 +1796,8 @@ QString QgsGeorefPluginGui::convertTransformEnumToString( QgsGeorefTransform::Tr
return tr( "Polynomial 3" );
case QgsGeorefTransform::ThinPlateSpline:
return tr( "Thin plate spline (TPS)" );
case QgsGeorefTransform::Projective:
return tr( "Projective" );
default:
return tr( "Not set" );
}
@@ -108,6 +108,35 @@ class QgsGDALGeorefTransform : public QgsGeorefTransformInterface
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 @@ QgsGeorefTransformInterface *QgsGeorefTransform::createImplementation( Transform
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 @@ void QgsGDALGeorefTransform::destroy_gdal_args()
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;
}
@@ -67,7 +67,8 @@ class QgsGeorefTransform : public QgsGeorefTransformInterface
PolynomialOrder2,
PolynomialOrder3,
ThinPlateSpline,
InvalidTransform
Projective,
InvalidTransform = 65535
};

QgsGeorefTransform( TransformParametrisation parametrisation );
@@ -173,3 +173,87 @@ void QgsLeastSquares::affine( std::vector<QgsPoint> mapCoords,

}


// 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 );
}


@@ -36,6 +36,10 @@ class QgsLeastSquares

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] );
};


@@ -36,6 +36,7 @@ QgsTransformSettingsDialog::QgsTransformSettingsDialog( const QString &raster, c

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 @@ void QgsTransformSettingsDialog::getTransformSettings( QgsGeorefTransform::Trans
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() );

0 comments on commit ca918e3

Please sign in to comment.