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/qgis@13987 c8812cc2-4d05-0410-92ff-de0c093fc19c

src/plugins/georeferencer/qgsgeorefplugingui.cpp

@@ -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" );
   }

src/plugins/georeferencer/qgsgeoreftransform.cpp

@@ -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;
+}

src/plugins/georeferencer/qgsgeoreftransform.h

@@ -67,7 +67,8 @@ class QgsGeorefTransform : public QgsGeorefTransformInterface
       PolynomialOrder2,
       PolynomialOrder3,
       ThinPlateSpline,
-      InvalidTransform
+      Projective,
+      InvalidTransform = 65535
     };
 
     QgsGeorefTransform( TransformParametrisation parametrisation );

src/plugins/georeferencer/qgsleastsquares.cpp

@@ -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 );
+}
+
+

src/plugins/georeferencer/qgsleastsquares.h

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

src/plugins/georeferencer/qgstransformsettingsdialog.cpp

@@ -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() );