QGIS API Documentation 3.28.0-Firenze (ed3ad0430f)
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
29bool 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;
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
93QgsGcpTransformerInterface *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
110bool 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
125bool 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
141{
142 return QgsLinearGeorefTransform::linearTransform;
143}
144
146{
147 return ( void * )&mParameters;
148}
149
151{
153}
154
155int 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//
198bool 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
228bool 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
243int 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
321QgsGDALGeorefTransform::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
339bool 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 const 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
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
420void 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
446bool 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 const 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 const 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
522int 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 const 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 const double X = ( x[i] * H[0] + y[i] * H[1] + H[2] ) / Z;
559 const 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.