[Backport release-3_22] Fix lots of georeferencer issues

src/app/georeferencer/qgsgcpcanvasitem.cpp

@@ -54,7 +54,7 @@ void QgsGCPCanvasItem::paint( QPainter *p )
   if ( mDataPoint )
   {
     enabled = mDataPoint->isEnabled();
-    worldCoords = mDataPoint->canvasCoords();
+    worldCoords = mDataPoint->destinationPoint();
     id = mDataPoint->id();
   }
   p->setOpacity( enabled ? 1.0 : 0.3 );
@@ -163,17 +163,21 @@ void QgsGCPCanvasItem::updatePosition()
 
   if ( mIsGCPSource )
   {
-    setPos( toCanvasCoordinates( mDataPoint->pixelCoords() ) );
-    return;
+    setPos( toCanvasCoordinates( mDataPoint->sourcePoint() ) );
   }
-  if ( mDataPoint->canvasCoords().isEmpty() )
+  else
   {
-    const QgsCoordinateReferenceSystem mapCrs = mMapCanvas->mapSettings().destinationCrs();
-    const QgsCoordinateTransform transf( mDataPoint->crs(), mapCrs, QgsProject::instance() );
-    const QgsPointXY mapCoords  = transf.transform( mDataPoint->mapCoords() );
-    mDataPoint->setCanvasCoords( mapCoords );
+    const QgsCoordinateTransform pointToCanvasTransform( mDataPoint->destinationPointCrs(), mMapCanvas->mapSettings().destinationCrs(), QgsProject::instance() );
+    try
+    {
+      const QgsPointXY canvasMapCoords = pointToCanvasTransform.transform( mDataPoint->destinationPoint() );
+      const QPointF canvasCoordinatesInPixels = toCanvasCoordinates( canvasMapCoords );
+
+      setPos( canvasCoordinatesInPixels );
+    }
+    catch ( QgsCsException & )
+    {}
   }
-  setPos( toCanvasCoordinates( mDataPoint->canvasCoords() ) );
 }
 
 void QgsGCPCanvasItem::drawResidualArrow( QPainter *p, const QgsRenderContext &context )
@@ -216,7 +220,7 @@ double QgsGCPCanvasItem::residualToScreenFactor() const
     }
   }
 
-  return 1.0 / ( mapUnitsPerScreenPixel * mapUnitsPerRasterPixel );
+  return mapUnitsPerRasterPixel / mapUnitsPerScreenPixel;
 }
 
 void QgsGCPCanvasItem::checkBoundingRectChange()

src/app/georeferencer/qgsgcplist.cpp

@@ -18,83 +18,232 @@
 #include "qgscoordinatereferencesystem.h"
 #include "qgscoordinatetransform.h"
 #include "qgsproject.h"
+#include "qgsgeoreftransform.h"
 
 #include "qgsgcplist.h"
+#include <QDir>
+#include <QTextStream>
 
-QgsGCPList::QgsGCPList( const QgsGCPList &list )
-  :  QList<QgsGeorefDataPoint *>()
+void QgsGCPList::createGCPVectors( QVector<QgsPointXY> &sourcePoints, QVector<QgsPointXY> &destinationPoints, const QgsCoordinateReferenceSystem &targetCrs, const QgsCoordinateTransformContext &context ) const
 {
-  clear();
-  QgsGCPList::const_iterator it = list.constBegin();
-  for ( ; it != list.constEnd(); ++it )
+  const int targetSize = countEnabledPoints();
+  sourcePoints.clear();
+  sourcePoints.reserve( targetSize );
+  destinationPoints.clear();
+  destinationPoints.reserve( targetSize );
+
+  for ( const QgsGeorefDataPoint *pt : std::as_const( *this ) )
   {
-    QgsGeorefDataPoint *pt = new QgsGeorefDataPoint( **it );
-    append( pt );
+    if ( !pt->isEnabled() )
+      continue;
+
+    sourcePoints.push_back( pt->sourcePoint() );
+    if ( targetCrs.isValid() )
+    {
+      destinationPoints.push_back( pt->transformedDestinationPoint( targetCrs, context ) );
+    }
+    else
+    {
+      destinationPoints.push_back( pt->destinationPoint() );
+    }
+  }
+}
+
+int QgsGCPList::countEnabledPoints() const
+{
+  if ( isEmpty() )
+    return 0;
+
+  int s = 0;
+  const_iterator it = begin();
+  while ( it != end() )
+  {
+    if ( ( *it )->isEnabled() )
+      s++;
+    ++it;
   }
+  return s;
 }
 
-void QgsGCPList::createGCPVectors( QVector<QgsPointXY> &mapCoords, QVector<QgsPointXY> &pixelCoords, const QgsCoordinateReferenceSystem targetCrs )
+void QgsGCPList::updateResiduals( QgsGeorefTransform *georefTransform, const QgsCoordinateReferenceSystem &targetCrs, const QgsCoordinateTransformContext &context, QgsUnitTypes::RenderUnit residualUnit )
 {
-  mapCoords   = QVector<QgsPointXY>( size() );
-  pixelCoords = QVector<QgsPointXY>( size() );
-  QgsPointXY transCoords;
-  for ( int i = 0, j = 0; i < sizeAll(); i++ )
+  bool bTransformUpdated = false;
+  QVector<QgsPointXY> sourceCoordinates;
+  QVector<QgsPointXY> destinationCoordinates;
+  createGCPVectors( sourceCoordinates, destinationCoordinates, targetCrs, context );
+
+  if ( georefTransform )
   {
-    QgsGeorefDataPoint *pt = at( i );
-    if ( pt->isEnabled() )
+    bTransformUpdated = georefTransform->updateParametersFromGcps( sourceCoordinates, destinationCoordinates, true );
+  }
+
+  // update residuals
+  for ( int i = 0; i < size(); ++i )
+  {
+    QgsGeorefDataPoint *p = at( i );
+
+    if ( !p )
+      continue;
+
+    p->setId( i );
+
+    const QgsPointXY transformedDestinationPoint = p->transformedDestinationPoint( targetCrs, QgsProject::instance()->transformContext() );
+
+    double dX = 0;
+    double dY = 0;
+    // Calculate residual if transform is available and up-to-date
+    if ( georefTransform && bTransformUpdated && georefTransform->parametersInitialized() )
     {
-      if ( targetCrs.isValid() )
+      QgsPointXY dst;
+      const QgsPointXY pixel = georefTransform->toSourcePixel( p->sourcePoint() );
+      if ( residualUnit == QgsUnitTypes::RenderPixels )
       {
-        try
+        // Transform from world to raster coordinate:
+        // This is the transform direction used by the warp operation.
+        // As transforms of order >=2 are not invertible, we are only
+        // interested in the residual in this direction
+        if ( georefTransform->transformWorldToRaster( transformedDestinationPoint, dst ) )
         {
-          transCoords =  QgsCoordinateTransform( pt->crs(), targetCrs,
-                                                 QgsProject::instance() ).transform( pt->mapCoords() );
-          mapCoords[j] = transCoords;
-          pt->setTransCoords( transCoords );
+          dX = ( dst.x() - pixel.x() );
+          dY = -( dst.y() - pixel.y() );
         }
-        catch ( const QgsException &e )
+      }
+      else if ( residualUnit == QgsUnitTypes::RenderMapUnits )
+      {
+        if ( georefTransform->transformRasterToWorld( pixel, dst ) )
         {
-          Q_UNUSED( e );
-          mapCoords[j] = pt->mapCoords();
+          dX = ( dst.x() - transformedDestinationPoint.x() );
+          dY = ( dst.y() - transformedDestinationPoint.y() );
         }
       }
-      else
-        mapCoords[j] = pt->mapCoords();
-      pixelCoords[j] = pt->pixelCoords();
-      j++;
     }
+    p->setResidual( QPointF( dX, dY ) );
   }
 }
 
-int QgsGCPList::size() const
+QList<QgsGcpPoint> QgsGCPList::asPoints() const
 {
-  if ( QList<QgsGeorefDataPoint *>::isEmpty() )
-    return 0;
-
-  int s = 0;
-  const_iterator it = begin();
-  while ( it != end() )
+  QList<QgsGcpPoint> res;
+  res.reserve( size() );
+  for ( QgsGeorefDataPoint *pt : *this )
   {
-    if ( ( *it )->isEnabled() )
-      s++;
-    ++it;
+    res.append( QgsGcpPoint( pt->point() ) );
   }
-  return s;
+  return res;
 }
 
-int QgsGCPList::sizeAll() const
+bool QgsGCPList::saveGcps( const QString &filePath, const QgsCoordinateReferenceSystem &targetCrs, const QgsCoordinateTransformContext &context, QString &error ) const
 {
-  return QList<QgsGeorefDataPoint *>::size();
+  error.clear();
+
+  QFile pointFile( filePath );
+  if ( pointFile.open( QIODevice::WriteOnly | QIODevice::Truncate ) )
+  {
+    QTextStream points( &pointFile );
+    if ( targetCrs.isValid() )
+    {
+      points << QStringLiteral( "#CRS: %1" ).arg( targetCrs.toWkt( QgsCoordinateReferenceSystem::WKT_PREFERRED ) );
+#if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
+      points << endl;
+#else
+      points << Qt::endl;
+#endif
+    }
+
+    points << "mapX,mapY,sourceX,sourceY,enable,dX,dY,residual";
+#if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
+    points << endl;
+#else
+    points << Qt::endl;
+#endif
+
+    for ( QgsGeorefDataPoint *pt : *this )
+    {
+      const QgsPointXY transformedDestinationPoint = pt->transformedDestinationPoint( targetCrs, context );
+      points << QStringLiteral( "%1,%2,%3,%4,%5,%6,%7,%8" )
+             .arg( qgsDoubleToString( transformedDestinationPoint.x() ),
+                   qgsDoubleToString( transformedDestinationPoint.y() ),
+                   qgsDoubleToString( pt->sourcePoint().x() ),
+                   qgsDoubleToString( pt->sourcePoint().y() ) )
+             .arg( pt->isEnabled() )
+             .arg( qgsDoubleToString( pt->residual().x() ),
+                   qgsDoubleToString( pt->residual().y() ),
+                   qgsDoubleToString( std::sqrt( pt->residual().x() * pt->residual().x() + pt->residual().y() * pt->residual().y() ) ) );
+#if QT_VERSION < QT_VERSION_CHECK(5, 14, 0)
+      points << endl;
+#else
+      points << Qt::endl;
+#endif
+    }
+    return true;
+  }
+  else
+  {
+    error = QObject::tr( "Could not write to GCP points file %1." ).arg( QDir::toNativeSeparators( filePath ) );
+    return false;
+  }
 }
 
-QgsGCPList &QgsGCPList::operator =( const QgsGCPList &list )
+QList<QgsGcpPoint> QgsGCPList::loadGcps( const QString &filePath, const QgsCoordinateReferenceSystem &defaultDestinationCrs, QgsCoordinateReferenceSystem &actualDestinationCrs, QString &error )
 {
-  clear();
-  QgsGCPList::const_iterator it = list.constBegin();
-  for ( ; it != list.constEnd(); ++it )
+  error.clear();
+  QFile pointFile( filePath );
+  if ( !pointFile.open( QIODevice::ReadOnly ) )
+  {
+    error = QObject::tr( "Could not open GCP points file %1." ).arg( QDir::toNativeSeparators( filePath ) );
+    return {};
+  }
+
+  QTextStream points( &pointFile );
+  int lineNumber = 0;
+  QString line = points.readLine();
+  lineNumber++;
+
+  int i = 0;
+  if ( line.contains( QLatin1String( "#CRS: " ) ) )
   {
-    QgsGeorefDataPoint *pt = new QgsGeorefDataPoint( **it );
-    append( pt );
+    const QString crsDef = line.remove( QStringLiteral( "#CRS: " ) );
+    if ( !crsDef.trimmed().isEmpty() )
+    {
+      actualDestinationCrs = QgsCoordinateReferenceSystem( crsDef );
+    }
+    else
+    {
+      actualDestinationCrs = defaultDestinationCrs;
+    }
+    line = points.readLine();
+    lineNumber++;
+  }
+  else
+    actualDestinationCrs = defaultDestinationCrs;
+
+  QList<QgsGcpPoint> res;
+  while ( !points.atEnd() )
+  {
+    line = points.readLine();
+    lineNumber++;
+    QStringList ls;
+    if ( line.contains( ',' ) ) // in previous format "\t" is delimiter of points in new - ","
+      ls = line.split( ',' ); // points from new georeferencer
+    else
+      ls = line.split( '\t' ); // points from prev georeferencer
+
+    if ( ls.count() < 4 )
+    {
+      error = QObject::tr( "Malformed content at line %1" ).arg( lineNumber );
+      return {};
+    }
+
+    const QgsPointXY destinationPoint( ls.at( 0 ).toDouble(), ls.at( 1 ).toDouble() ); // map x,y
+    const QgsPointXY sourcePoint( ls.at( 2 ).toDouble(), ls.at( 3 ).toDouble() ); // source x,y
+    bool enable = true;
+    if ( ls.count() >= 5 )
+    {
+      enable = ls.at( 4 ).toInt();
+    }
+    res.append( QgsGcpPoint( sourcePoint, destinationPoint, actualDestinationCrs, enable ) );
+
+    ++i;
   }
-  return *this;
+  return res;
 }