QGIS API Documentation  3.20.0-Odense (decaadbb31)
qgsgcptransformer.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  qgsgcptransformer.cpp
3  --------------------------------------
4  Date : 18-Feb-2009
5  Copyright : (c) 2009 by Manuel Massing
6  Email : m.massing at warped-space.de
7  ***************************************************************************
8  * *
9  * This program is free software; you can redistribute it and/or modify *
10  * it under the terms of the GNU General Public License as published by *
11  * the Free Software Foundation; either version 2 of the License, or *
12  * (at your option) any later version. *
13  * *
14  ***************************************************************************/
15 
16 #include "qgsgcptransformer.h"
17 
18 #include <gdal.h>
19 #include <gdal_alg.h>
20 
21 #include "qgsleastsquares.h"
22 
23 #include <cmath>
24 
25 #include <cassert>
26 #include <limits>
27 
28 
29 bool QgsGcpTransformerInterface::transform( double &x, double &y, bool inverseTransform ) const
30 {
31  GDALTransformerFunc t = GDALTransformer();
32  // Fail if no transformer function was returned
33  if ( !t )
34  return false;
35 
36  double z = 0.0;
37  int success = 0;
38 
39  // Call GDAL transform function
40  ( *t )( GDALTransformerArgs(), inverseTransform ? 1 : 0, 1, &x, &y, &z, &success );
41  if ( !success )
42  return false;
43 
44  return true;
45 }
46 
48 {
49  switch ( method )
50  {
52  return QObject::tr( "Linear" );
54  return QObject::tr( "Helmert" );
56  return QObject::tr( "Polynomial 1" );
58  return QObject::tr( "Polynomial 2" );
60  return QObject::tr( "Polynomial 3" );
62  return QObject::tr( "Thin Plate Spline (TPS)" );
64  return QObject::tr( "Projective" );
65  default:
66  return QObject::tr( "Not set" );
67  }
68 }
69 
71 {
72  switch ( method )
73  {
75  return new QgsLinearGeorefTransform;
77  return new QgsHelmertGeorefTransform;
79  return new QgsGDALGeorefTransform( false, 1 );
81  return new QgsGDALGeorefTransform( false, 2 );
83  return new QgsGDALGeorefTransform( false, 3 );
85  return new QgsGDALGeorefTransform( true, 0 );
88  default:
89  return nullptr;
90  }
91 }
92 
93 QgsGcpTransformerInterface *QgsGcpTransformerInterface::createFromParameters( QgsGcpTransformerInterface::TransformMethod method, const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates )
94 {
95  std::unique_ptr< QgsGcpTransformerInterface > transformer( create( method ) );
96  if ( !transformer )
97  return nullptr;
98 
99  if ( !transformer->updateParametersFromGcps( sourceCoordinates, destinationCoordinates ) )
100  return nullptr;
101 
102  return transformer.release();
103 }
104 
105 
106 //
107 // QgsLinearGeorefTransform
108 //
109 
110 bool QgsLinearGeorefTransform::getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const
111 {
112  origin = mParameters.origin;
113  scaleX = mParameters.scaleX;
114  scaleY = mParameters.scaleY;
115  return true;
116 }
117 
119 {
120  std::unique_ptr< QgsLinearGeorefTransform > res = std::make_unique< QgsLinearGeorefTransform >();
121  res->mParameters = mParameters;
122  return res.release();
123 }
124 
125 bool QgsLinearGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
126 {
127  if ( destinationCoordinates.size() < minimumGcpCount() )
128  return false;
129 
130  mParameters.invertYAxis = invertYAxis;
131  QgsLeastSquares::linear( sourceCoordinates, destinationCoordinates, mParameters.origin, mParameters.scaleX, mParameters.scaleY );
132  return true;
133 }
134 
136 {
137  return 2;
138 }
139 
140 GDALTransformerFunc QgsLinearGeorefTransform::GDALTransformer() const
141 {
142  return QgsLinearGeorefTransform::linearTransform;
143 }
144 
146 {
147  return ( void * )&mParameters;
148 }
149 
151 {
153 }
154 
155 int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
156  double *x, double *y, double *z, int *panSuccess )
157 {
158  Q_UNUSED( z )
159  LinearParameters *t = static_cast<LinearParameters *>( pTransformerArg );
160  if ( !t )
161  return false;
162 
163  if ( !bDstToSrc )
164  {
165  for ( int i = 0; i < nPointCount; ++i )
166  {
167  x[i] = x[i] * t->scaleX + t->origin.x();
168  y[i] = ( t->invertYAxis ? -1 : 1 ) * y[i] * t->scaleY + t->origin.y();
169  panSuccess[i] = true;
170  }
171  }
172  else
173  {
174  // Guard against division by zero
175  if ( std::fabs( t->scaleX ) < std::numeric_limits<double>::epsilon() ||
176  std::fabs( t->scaleY ) < std::numeric_limits<double>::epsilon() )
177  {
178  for ( int i = 0; i < nPointCount; ++i )
179  {
180  panSuccess[i] = false;
181  }
182  return false;
183  }
184  for ( int i = 0; i < nPointCount; ++i )
185  {
186  x[i] = ( x[i] - t->origin.x() ) / t->scaleX;
187  y[i] = ( y[i] - t->origin.y() ) / ( ( t->invertYAxis ? -1 : 1 ) * t->scaleY );
188  panSuccess[i] = true;
189  }
190  }
191 
192  return true;
193 }
194 
195 //
196 // QgsHelmertGeorefTransform
197 //
198 bool QgsHelmertGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
199 {
200  if ( destinationCoordinates.size() < minimumGcpCount() )
201  return false;
202 
203  mHelmertParameters.invertYAxis = invertYAxis;
204  QgsLeastSquares::helmert( sourceCoordinates, destinationCoordinates, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle );
205  return true;
206 }
207 
209 {
210  return 2;
211 }
212 
214 {
215  return QgsHelmertGeorefTransform::helmertTransform;
216 }
217 
219 {
220  return ( void * )&mHelmertParameters;
221 }
222 
224 {
226 }
227 
228 bool QgsHelmertGeorefTransform::getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const
229 {
230  origin = mHelmertParameters.origin;
231  scale = mHelmertParameters.scale;
232  rotation = mHelmertParameters.angle;
233  return true;
234 }
235 
237 {
238  std::unique_ptr< QgsHelmertGeorefTransform > res = std::make_unique< QgsHelmertGeorefTransform >();
239  res->mHelmertParameters = mHelmertParameters;
240  return res.release();
241 }
242 
243 int QgsHelmertGeorefTransform::helmertTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
244  double *x, double *y, double *z, int *panSuccess )
245 {
246  Q_UNUSED( z )
247  const HelmertParameters *t = static_cast< const HelmertParameters *>( pTransformerArg );
248  if ( !t )
249  return false;
250 
251  double a = std::cos( t->angle );
252  double b = std::sin( t->angle );
253  const double x0 = t->origin.x();
254  const double y0 = t->origin.y();
255  const double s = t->scale;
256  if ( !bDstToSrc )
257  {
258  a *= s;
259  b *= s;
260  for ( int i = 0; i < nPointCount; ++i )
261  {
262  const double xT = x[i];
263  const double yT = y[i];
264 
265  if ( t->invertYAxis )
266  {
267  // Because rotation parameters where estimated in a CS with negative y-axis ^= down.
268  // we need to apply the rotation matrix and a change of base:
269  // |cos a, -sin a | |1, 0| | cos a, sin a|
270  // |sin a, cos a | |0,-1| = | sin a, -cos a|
271  x[i] = x0 + ( a * xT + b * yT );
272  y[i] = y0 + ( b * xT - a * yT );
273  }
274  else
275  {
276  x[i] = x0 + ( a * xT - b * yT );
277  y[i] = y0 + ( b * xT + a * yT );
278  }
279  panSuccess[i] = true;
280  }
281  }
282  else
283  {
284  // Guard against division by zero
285  if ( std::fabs( s ) < std::numeric_limits<double>::epsilon() )
286  {
287  for ( int i = 0; i < nPointCount; ++i )
288  {
289  panSuccess[i] = false;
290  }
291  return false;
292  }
293  a /= s;
294  b /= s;
295  for ( int i = 0; i < nPointCount; ++i )
296  {
297  const double xT = x[i] - x0;
298  const double yT = y[i] - y0;
299  if ( t->invertYAxis )
300  {
301  // | cos a, sin a |^-1 |cos a, sin a|
302  // | sin a, -cos a | = |sin a, -cos a|
303  x[i] = a * xT + b * yT;
304  y[i] = b * xT - a * yT;
305  }
306  else
307  {
308  x[i] = a * xT + b * yT;
309  y[i] = -b * xT + a * yT;
310  }
311  panSuccess[i] = true;
312  }
313  }
314  return true;
315 }
316 
317 //
318 // QgsGDALGeorefTransform
319 //
320 
321 QgsGDALGeorefTransform::QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder )
322  : mPolynomialOrder( std::min( 3u, polynomialOrder ) )
323  , mIsTPSTransform( useTPS )
324 {
325 }
326 
328 {
329  destroyGdalArgs();
330 }
331 
333 {
334  std::unique_ptr< QgsGDALGeorefTransform > res = std::make_unique< QgsGDALGeorefTransform >( mIsTPSTransform, mPolynomialOrder );
335  res->updateParametersFromGcps( mSourceCoords, mDestCoordinates, mInvertYAxis );
336  return res.release();
337 }
338 
339 bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
340 {
341  mSourceCoords = sourceCoordinates;
342  mDestCoordinates = destinationCoordinates;
343  mInvertYAxis = invertYAxis;
344 
345  assert( sourceCoordinates.size() == destinationCoordinates.size() );
346  if ( sourceCoordinates.size() != destinationCoordinates.size() )
347  return false;
348  int n = sourceCoordinates.size();
349 
350  GDAL_GCP *GCPList = new GDAL_GCP[n];
351  for ( int i = 0; i < n; i++ )
352  {
353  GCPList[i].pszId = new char[20];
354  snprintf( GCPList[i].pszId, 19, "gcp%i", i );
355  GCPList[i].pszInfo = nullptr;
356  GCPList[i].dfGCPPixel = sourceCoordinates[i].x();
357  GCPList[i].dfGCPLine = ( mInvertYAxis ? -1 : 1 ) * sourceCoordinates[i].y();
358  GCPList[i].dfGCPX = destinationCoordinates[i].x();
359  GCPList[i].dfGCPY = destinationCoordinates[i].y();
360  GCPList[i].dfGCPZ = 0;
361  }
362  destroyGdalArgs();
363 
364  if ( mIsTPSTransform )
365  mGDALTransformerArgs = GDALCreateTPSTransformer( n, GCPList, false );
366  else
367  mGDALTransformerArgs = GDALCreateGCPTransformer( n, GCPList, mPolynomialOrder, false );
368 
369  for ( int i = 0; i < n; i++ )
370  {
371  delete [] GCPList[i].pszId;
372  }
373  delete [] GCPList;
374 
375  return nullptr != mGDALTransformerArgs;
376 }
377 
379 {
380  if ( mIsTPSTransform )
381  return 1;
382  else
383  return ( ( mPolynomialOrder + 2 ) * ( mPolynomialOrder + 1 ) ) / 2;
384 }
385 
386 GDALTransformerFunc QgsGDALGeorefTransform::GDALTransformer() const
387 {
388  // Fail if no arguments were calculated through updateParametersFromGCP
389  if ( !mGDALTransformerArgs )
390  return nullptr;
391 
392  if ( mIsTPSTransform )
393  return GDALTPSTransform;
394  else
395  return GDALGCPTransform;
396 }
397 
399 {
400  return mGDALTransformerArgs;
401 }
402 
404 {
405  if ( mIsTPSTransform )
407 
408  switch ( mPolynomialOrder )
409  {
410  case 1:
412  case 2:
414  case 3:
416  }
418 }
419 
420 void QgsGDALGeorefTransform::destroyGdalArgs()
421 {
422  if ( mGDALTransformerArgs )
423  {
424  if ( mIsTPSTransform )
425  GDALDestroyTPSTransformer( mGDALTransformerArgs );
426  else
427  GDALDestroyGCPTransformer( mGDALTransformerArgs );
428  }
429 }
430 
431 //
432 // QgsProjectiveGeorefTransform
433 //
434 
436  : mParameters()
437 {}
438 
440 {
441  std::unique_ptr< QgsProjectiveGeorefTransform > res = std::make_unique< QgsProjectiveGeorefTransform >();
442  res->mParameters = mParameters;
443  return res.release();
444 }
445 
446 bool QgsProjectiveGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
447 {
448  if ( destinationCoordinates.size() < minimumGcpCount() )
449  return false;
450 
451  if ( invertYAxis )
452  {
453  // HACK: flip y coordinates, because georeferencer and gdal use different conventions
454  QVector<QgsPointXY> flippedPixelCoords;
455  flippedPixelCoords.reserve( sourceCoordinates.size() );
456  for ( const QgsPointXY &coord : sourceCoordinates )
457  {
458  flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() );
459  }
460 
461  QgsLeastSquares::projective( flippedPixelCoords, destinationCoordinates, mParameters.H );
462  }
463  else
464  {
465  QgsLeastSquares::projective( sourceCoordinates, destinationCoordinates, mParameters.H );
466  }
467 
468  // Invert the homography matrix using adjoint matrix
469  double *H = mParameters.H;
470 
471  double adjoint[9];
472  adjoint[0] = H[4] * H[8] - H[5] * H[7];
473  adjoint[1] = -H[1] * H[8] + H[2] * H[7];
474  adjoint[2] = H[1] * H[5] - H[2] * H[4];
475 
476  adjoint[3] = -H[3] * H[8] + H[5] * H[6];
477  adjoint[4] = H[0] * H[8] - H[2] * H[6];
478  adjoint[5] = -H[0] * H[5] + H[2] * H[3];
479 
480  adjoint[6] = H[3] * H[7] - H[4] * H[6];
481  adjoint[7] = -H[0] * H[7] + H[1] * H[6];
482  adjoint[8] = H[0] * H[4] - H[1] * H[3];
483 
484  double det = H[0] * adjoint[0] + H[3] * adjoint[1] + H[6] * adjoint[2];
485 
486  if ( std::fabs( det ) < 1024.0 * std::numeric_limits<double>::epsilon() )
487  {
488  mParameters.hasInverse = false;
489  }
490  else
491  {
492  mParameters.hasInverse = true;
493  double oo_det = 1.0 / det;
494  for ( int i = 0; i < 9; i++ )
495  {
496  mParameters.Hinv[i] = adjoint[i] * oo_det;
497  }
498  }
499  return true;
500 }
501 
503 {
504  return 4;
505 }
506 
508 {
509  return QgsProjectiveGeorefTransform::projectiveTransform;
510 }
511 
513 {
514  return ( void * )&mParameters;
515 }
516 
518 {
520 }
521 
522 int QgsProjectiveGeorefTransform::projectiveTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
523  double *x, double *y, double *z, int *panSuccess )
524 {
525  Q_UNUSED( z )
526  ProjectiveParameters *t = static_cast<ProjectiveParameters *>( pTransformerArg );
527  if ( !t )
528  return false;
529 
530  double *H = nullptr;
531  if ( !bDstToSrc )
532  {
533  H = t->H;
534  }
535  else
536  {
537  if ( !t->hasInverse )
538  {
539  for ( int i = 0; i < nPointCount; ++i )
540  {
541  panSuccess[i] = false;
542  }
543  return false;
544  }
545  H = t->Hinv;
546  }
547 
548 
549  for ( int i = 0; i < nPointCount; ++i )
550  {
551  double Z = x[i] * H[6] + y[i] * H[7] + H[8];
552  // Projects to infinity?
553  if ( std::fabs( Z ) < 1024.0 * std::numeric_limits<double>::epsilon() )
554  {
555  panSuccess[i] = false;
556  continue;
557  }
558  double X = ( x[i] * H[0] + y[i] * H[1] + H[2] ) / Z;
559  double Y = ( x[i] * H[3] + y[i] * H[4] + H[5] ) / Z;
560 
561  x[i] = X;
562  y[i] = Y;
563 
564  panSuccess[i] = true;
565  }
566 
567  return true;
568 }
Interface to gdal thin plate splines and 1st/2nd/3rd order polynomials.
int minimumGcpCount() const override
Returns the minimum number of Ground Control Points (GCPs) required for parameter fitting.
TransformMethod method() const override
Returns the transformation method.
QgsGcpTransformerInterface * clone() const override
Clones the transformer, returning a new copy of the transformer with the same parameters as this one.
QgsGDALGeorefTransform(bool useTPS, unsigned int polynomialOrder)
Constructor for QgsGDALGeorefTransform.
GDALTransformerFunc GDALTransformer() const override
Returns function pointer to the GDALTransformer function.
void * GDALTransformerArgs() const override
Returns pointer to the GDALTransformer arguments.
bool updateParametersFromGcps(const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates, bool invertYAxis=false) override
Fits transformation parameters using the specified Ground Control Points (GCPs) lists of source and d...
An interface for Ground Control Points (GCP) based transformations.
static QgsGcpTransformerInterface * createFromParameters(TransformMethod method, const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates) SIP_THROW(QgsNotSupportedException)
Creates a new QgsGcpTransformerInterface subclass representing the specified transform method,...
bool transform(double &x, double &y, bool inverseTransform=false) const
Transforms the point (x, y) from source to destination coordinates.
TransformMethod
Available transformation methods.
static QgsGcpTransformerInterface * create(TransformMethod method)
Creates a new QgsGcpTransformerInterface subclass representing the specified transform method.
static QString methodToString(TransformMethod method)
Returns a translated string representing the specified transform method.
virtual void * GDALTransformerArgs() const =0
Returns pointer to the GDALTransformer arguments.
virtual GDALTransformerFunc GDALTransformer() const =0
Returns function pointer to the GDALTransformer function.
virtual TransformMethod method() const =0
Returns the transformation method.
2-dimensional helmert transform, parametrised by isotropic scale, rotation angle and translation.
TransformMethod method() const override
Returns the transformation method.
QgsGcpTransformerInterface * clone() const override
Clones the transformer, returning a new copy of the transformer with the same parameters as this one.
bool getOriginScaleRotation(QgsPointXY &origin, double &scale, double &rotation) const
Returns the origin, scale and rotation for the transform.
GDALTransformerFunc GDALTransformer() const override
Returns function pointer to the GDALTransformer function.
int minimumGcpCount() const override
Returns the minimum number of Ground Control Points (GCPs) required for parameter fitting.
bool updateParametersFromGcps(const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates, bool invertYAxis=false) override
Fits transformation parameters using the specified Ground Control Points (GCPs) lists of source and d...
void * GDALTransformerArgs() const override
Returns pointer to the GDALTransformer arguments.
static void helmert(const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates, QgsPointXY &origin, double &pixelSize, double &rotation)
Transforms the point at origin in-place, using a helmert transformation calculated from the list of s...
static void projective(const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates, double H[9])
Calculates projective parameters from the list of source and destination Ground Control Points (GCPs)...
static void linear(const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates, QgsPointXY &origin, double &pixelXSize, double &pixelYSize)
Transforms the point at origin in-place, using a linear transformation calculated from the list of so...
A simple transform which is parametrized by a translation and anistotropic scale.
bool getOriginScale(QgsPointXY &origin, double &scaleX, double &scaleY) const
Returns the origin and scale for the transform.
TransformMethod method() const override
Returns the transformation method.
void * GDALTransformerArgs() const override
Returns pointer to the GDALTransformer arguments.
QgsGcpTransformerInterface * clone() const override
Clones the transformer, returning a new copy of the transformer with the same parameters as this one.
bool updateParametersFromGcps(const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates, bool invertYAxis=false) override
Fits transformation parameters using the specified Ground Control Points (GCPs) lists of source and d...
int minimumGcpCount() const override
Returns the minimum number of Ground Control Points (GCPs) required for parameter fitting.
GDALTransformerFunc GDALTransformer() const override
Returns function pointer to the GDALTransformer function.
A class to represent a 2D point.
Definition: qgspointxy.h:59
A planar projective transform, expressed by a homography.
int minimumGcpCount() const override
Returns the minimum number of Ground Control Points (GCPs) required for parameter fitting.
QgsProjectiveGeorefTransform()
Constructor for QgsProjectiveGeorefTransform.
QgsGcpTransformerInterface * clone() const override
Clones the transformer, returning a new copy of the transformer with the same parameters as this one.
bool updateParametersFromGcps(const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates, bool invertYAxis=false) override
Fits transformation parameters using the specified Ground Control Points (GCPs) lists of source and d...
GDALTransformerFunc GDALTransformer() const override
Returns function pointer to the GDALTransformer function.
void * GDALTransformerArgs() const override
Returns pointer to the GDALTransformer arguments.
TransformMethod method() const override
Returns the transformation method.