QGIS API Documentation 3.41.0-Master (fda2aa46e9a)
Loading...
Searching...
No Matches
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#include "moc_qgsgcptransformer.cpp"
18
19#include <gdal.h>
20#include <gdal_alg.h>
21
22#include "qgsleastsquares.h"
23
24#include <cmath>
25
26#include <cassert>
27#include <limits>
28
29
30bool QgsGcpTransformerInterface::transform( double &x, double &y, bool inverseTransform ) const
31{
32 GDALTransformerFunc t = GDALTransformer();
33 // Fail if no transformer function was returned
34 if ( !t )
35 return false;
36
37 double z = 0.0;
38 int success = 0;
39
40 // Call GDAL transform function
41 ( *t )( GDALTransformerArgs(), inverseTransform ? 1 : 0, 1, &x, &y, &z, &success );
42 if ( !success )
43 return false;
44
45 return true;
46}
47
49{
50 switch ( method )
51 {
53 return QObject::tr( "Linear" );
55 return QObject::tr( "Helmert" );
57 return QObject::tr( "Polynomial 1" );
59 return QObject::tr( "Polynomial 2" );
61 return QObject::tr( "Polynomial 3" );
63 return QObject::tr( "Thin Plate Spline (TPS)" );
65 return QObject::tr( "Projective" );
66 default:
67 return QObject::tr( "Not set" );
68 }
69}
70
93
94QgsGcpTransformerInterface *QgsGcpTransformerInterface::createFromParameters( QgsGcpTransformerInterface::TransformMethod method, const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates )
95{
96 std::unique_ptr< QgsGcpTransformerInterface > transformer( create( method ) );
97 if ( !transformer )
98 return nullptr;
99
100 if ( !transformer->updateParametersFromGcps( sourceCoordinates, destinationCoordinates ) )
101 return nullptr;
102
103 return transformer.release();
104}
105
106
107//
108// QgsLinearGeorefTransform
109//
110
111bool QgsLinearGeorefTransform::getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const
112{
113 origin = mParameters.origin;
114 scaleX = mParameters.scaleX;
115 scaleY = mParameters.scaleY;
116 return true;
117}
118
120{
121 std::unique_ptr< QgsLinearGeorefTransform > res = std::make_unique< QgsLinearGeorefTransform >();
122 res->mParameters = mParameters;
123 return res.release();
124}
125
126bool QgsLinearGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
127{
128 if ( destinationCoordinates.size() < minimumGcpCount() )
129 return false;
130
131 mParameters.invertYAxis = invertYAxis;
132 QgsLeastSquares::linear( sourceCoordinates, destinationCoordinates, mParameters.origin, mParameters.scaleX, mParameters.scaleY );
133 return true;
134}
135
137{
138 return 2;
139}
140
142{
143 return QgsLinearGeorefTransform::linearTransform;
144}
145
147{
148 return ( void * )&mParameters;
149}
150
155
156int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
157 double *x, double *y, double *z, int *panSuccess )
158{
159 Q_UNUSED( z )
160 LinearParameters *t = static_cast<LinearParameters *>( pTransformerArg );
161 if ( !t )
162 return false;
163
164 if ( !bDstToSrc )
165 {
166 for ( int i = 0; i < nPointCount; ++i )
167 {
168 x[i] = x[i] * t->scaleX + t->origin.x();
169 y[i] = ( t->invertYAxis ? -1 : 1 ) * y[i] * t->scaleY + t->origin.y();
170 panSuccess[i] = true;
171 }
172 }
173 else
174 {
175 // Guard against division by zero
176 if ( std::fabs( t->scaleX ) < std::numeric_limits<double>::epsilon() ||
177 std::fabs( t->scaleY ) < std::numeric_limits<double>::epsilon() )
178 {
179 for ( int i = 0; i < nPointCount; ++i )
180 {
181 panSuccess[i] = false;
182 }
183 return false;
184 }
185 for ( int i = 0; i < nPointCount; ++i )
186 {
187 x[i] = ( x[i] - t->origin.x() ) / t->scaleX;
188 y[i] = ( y[i] - t->origin.y() ) / ( ( t->invertYAxis ? -1 : 1 ) * t->scaleY );
189 panSuccess[i] = true;
190 }
191 }
192
193 return true;
194}
195
196//
197// QgsHelmertGeorefTransform
198//
199bool QgsHelmertGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
200{
201 if ( destinationCoordinates.size() < minimumGcpCount() )
202 return false;
203
204 mHelmertParameters.invertYAxis = invertYAxis;
205 QgsLeastSquares::helmert( sourceCoordinates, destinationCoordinates, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle );
206 return true;
207}
208
210{
211 return 2;
212}
213
215{
216 return QgsHelmertGeorefTransform::helmertTransform;
217}
218
220{
221 return ( void * )&mHelmertParameters;
222}
223
228
229bool QgsHelmertGeorefTransform::getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const
230{
231 origin = mHelmertParameters.origin;
232 scale = mHelmertParameters.scale;
233 rotation = mHelmertParameters.angle;
234 return true;
235}
236
238{
239 std::unique_ptr< QgsHelmertGeorefTransform > res = std::make_unique< QgsHelmertGeorefTransform >();
240 res->mHelmertParameters = mHelmertParameters;
241 return res.release();
242}
243
244int QgsHelmertGeorefTransform::helmertTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
245 double *x, double *y, double *z, int *panSuccess )
246{
247 Q_UNUSED( z )
248 const HelmertParameters *t = static_cast< const HelmertParameters *>( pTransformerArg );
249 if ( !t )
250 return false;
251
252 double a = std::cos( t->angle );
253 double b = std::sin( t->angle );
254 const double x0 = t->origin.x();
255 const double y0 = t->origin.y();
256 const double s = t->scale;
257 if ( !bDstToSrc )
258 {
259 a *= s;
260 b *= s;
261 for ( int i = 0; i < nPointCount; ++i )
262 {
263 const double xT = x[i];
264 const double yT = y[i];
265
266 if ( t->invertYAxis )
267 {
268 // Because rotation parameters where estimated in a CS with negative y-axis ^= down.
269 // we need to apply the rotation matrix and a change of base:
270 // |cos a, -sin a | |1, 0| | cos a, sin a|
271 // |sin a, cos a | |0,-1| = | sin a, -cos a|
272 x[i] = x0 + ( a * xT + b * yT );
273 y[i] = y0 + ( b * xT - a * yT );
274 }
275 else
276 {
277 x[i] = x0 + ( a * xT - b * yT );
278 y[i] = y0 + ( b * xT + a * yT );
279 }
280 panSuccess[i] = true;
281 }
282 }
283 else
284 {
285 // Guard against division by zero
286 if ( std::fabs( s ) < std::numeric_limits<double>::epsilon() )
287 {
288 for ( int i = 0; i < nPointCount; ++i )
289 {
290 panSuccess[i] = false;
291 }
292 return false;
293 }
294 a /= s;
295 b /= s;
296 for ( int i = 0; i < nPointCount; ++i )
297 {
298 const double xT = x[i] - x0;
299 const double yT = y[i] - y0;
300 if ( t->invertYAxis )
301 {
302 // | cos a, sin a |^-1 |cos a, sin a|
303 // | sin a, -cos a | = |sin a, -cos a|
304 x[i] = a * xT + b * yT;
305 y[i] = b * xT - a * yT;
306 }
307 else
308 {
309 x[i] = a * xT + b * yT;
310 y[i] = -b * xT + a * yT;
311 }
312 panSuccess[i] = true;
313 }
314 }
315 return true;
316}
317
318//
319// QgsGDALGeorefTransform
320//
321
322QgsGDALGeorefTransform::QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder )
323 : mPolynomialOrder( std::min( 3u, polynomialOrder ) )
324 , mIsTPSTransform( useTPS )
325{
326}
327
329{
330 destroyGdalArgs();
331}
332
334{
335 std::unique_ptr< QgsGDALGeorefTransform > res = std::make_unique< QgsGDALGeorefTransform >( mIsTPSTransform, mPolynomialOrder );
336 res->updateParametersFromGcps( mSourceCoords, mDestCoordinates, mInvertYAxis );
337 return res.release();
338}
339
340bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
341{
342 mSourceCoords = sourceCoordinates;
343 mDestCoordinates = destinationCoordinates;
344 mInvertYAxis = invertYAxis;
345
346 assert( sourceCoordinates.size() == destinationCoordinates.size() );
347 if ( sourceCoordinates.size() != destinationCoordinates.size() )
348 return false;
349 const int n = sourceCoordinates.size();
350
351 GDAL_GCP *GCPList = new GDAL_GCP[n];
352 for ( int i = 0; i < n; i++ )
353 {
354 GCPList[i].pszId = new char[20];
355 snprintf( GCPList[i].pszId, 19, "gcp%i", i );
356 GCPList[i].pszInfo = nullptr;
357 GCPList[i].dfGCPPixel = sourceCoordinates[i].x();
358 GCPList[i].dfGCPLine = ( mInvertYAxis ? -1 : 1 ) * sourceCoordinates[i].y();
359 GCPList[i].dfGCPX = destinationCoordinates[i].x();
360 GCPList[i].dfGCPY = destinationCoordinates[i].y();
361 GCPList[i].dfGCPZ = 0;
362 }
363 destroyGdalArgs();
364
365 if ( mIsTPSTransform )
366 mGDALTransformerArgs = GDALCreateTPSTransformer( n, GCPList, false );
367 else
368 mGDALTransformerArgs = GDALCreateGCPTransformer( n, GCPList, mPolynomialOrder, false );
369
370 for ( int i = 0; i < n; i++ )
371 {
372 delete [] GCPList[i].pszId;
373 }
374 delete [] GCPList;
375
376 return nullptr != mGDALTransformerArgs;
377}
378
380{
381 if ( mIsTPSTransform )
382 return 1;
383 else
384 return ( ( mPolynomialOrder + 2 ) * ( mPolynomialOrder + 1 ) ) / 2;
385}
386
388{
389 // Fail if no arguments were calculated through updateParametersFromGCP
390 if ( !mGDALTransformerArgs )
391 return nullptr;
392
393 if ( mIsTPSTransform )
394 return GDALTPSTransform;
395 else
396 return GDALGCPTransform;
397}
398
400{
401 return mGDALTransformerArgs;
402}
403
405{
406 if ( mIsTPSTransform )
408
409 switch ( mPolynomialOrder )
410 {
411 case 1:
413 case 2:
415 case 3:
417 }
419}
420
421void QgsGDALGeorefTransform::destroyGdalArgs()
422{
423 if ( mGDALTransformerArgs )
424 {
425 if ( mIsTPSTransform )
426 GDALDestroyTPSTransformer( mGDALTransformerArgs );
427 else
428 GDALDestroyGCPTransformer( mGDALTransformerArgs );
429 }
430}
431
432//
433// QgsProjectiveGeorefTransform
434//
435
439
441{
442 std::unique_ptr< QgsProjectiveGeorefTransform > res = std::make_unique< QgsProjectiveGeorefTransform >();
443 res->mParameters = mParameters;
444 return res.release();
445}
446
447bool QgsProjectiveGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
448{
449 if ( destinationCoordinates.size() < minimumGcpCount() )
450 return false;
451
452 if ( invertYAxis )
453 {
454 // HACK: flip y coordinates, because georeferencer and gdal use different conventions
455 QVector<QgsPointXY> flippedPixelCoords;
456 flippedPixelCoords.reserve( sourceCoordinates.size() );
457 for ( const QgsPointXY &coord : sourceCoordinates )
458 {
459 flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() );
460 }
461
462 QgsLeastSquares::projective( flippedPixelCoords, destinationCoordinates, mParameters.H );
463 }
464 else
465 {
466 QgsLeastSquares::projective( sourceCoordinates, destinationCoordinates, mParameters.H );
467 }
468
469 // Invert the homography matrix using adjoint matrix
470 double *H = mParameters.H;
471
472 double adjoint[9];
473 adjoint[0] = H[4] * H[8] - H[5] * H[7];
474 adjoint[1] = -H[1] * H[8] + H[2] * H[7];
475 adjoint[2] = H[1] * H[5] - H[2] * H[4];
476
477 adjoint[3] = -H[3] * H[8] + H[5] * H[6];
478 adjoint[4] = H[0] * H[8] - H[2] * H[6];
479 adjoint[5] = -H[0] * H[5] + H[2] * H[3];
480
481 adjoint[6] = H[3] * H[7] - H[4] * H[6];
482 adjoint[7] = -H[0] * H[7] + H[1] * H[6];
483 adjoint[8] = H[0] * H[4] - H[1] * H[3];
484
485 const double det = H[0] * adjoint[0] + H[3] * adjoint[1] + H[6] * adjoint[2];
486
487 if ( std::fabs( det ) < 1024.0 * std::numeric_limits<double>::epsilon() )
488 {
489 mParameters.hasInverse = false;
490 }
491 else
492 {
493 mParameters.hasInverse = true;
494 const double oo_det = 1.0 / det;
495 for ( int i = 0; i < 9; i++ )
496 {
497 mParameters.Hinv[i] = adjoint[i] * oo_det;
498 }
499 }
500 return true;
501}
502
504{
505 return 4;
506}
507
509{
510 return QgsProjectiveGeorefTransform::projectiveTransform;
511}
512
514{
515 return ( void * )&mParameters;
516}
517
522
523int QgsProjectiveGeorefTransform::projectiveTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
524 double *x, double *y, double *z, int *panSuccess )
525{
526 Q_UNUSED( z )
527 ProjectiveParameters *t = static_cast<ProjectiveParameters *>( pTransformerArg );
528 if ( !t )
529 return false;
530
531 double *H = nullptr;
532 if ( !bDstToSrc )
533 {
534 H = t->H;
535 }
536 else
537 {
538 if ( !t->hasInverse )
539 {
540 for ( int i = 0; i < nPointCount; ++i )
541 {
542 panSuccess[i] = false;
543 }
544 return false;
545 }
546 H = t->Hinv;
547 }
548
549
550 for ( int i = 0; i < nPointCount; ++i )
551 {
552 const double Z = x[i] * H[6] + y[i] * H[7] + H[8];
553 // Projects to infinity?
554 if ( std::fabs( Z ) < 1024.0 * std::numeric_limits<double>::epsilon() )
555 {
556 panSuccess[i] = false;
557 continue;
558 }
559 const double X = ( x[i] * H[0] + y[i] * H[1] + H[2] ) / Z;
560 const double Y = ( x[i] * H[3] + y[i] * H[4] + H[5] ) / Z;
561
562 x[i] = X;
563 y[i] = Y;
564
565 panSuccess[i] = true;
566 }
567
568 return true;
569}
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.
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.
static QgsGcpTransformerInterface * createFromParameters(TransformMethod method, const QVector< QgsPointXY > &sourceCoordinates, const QVector< QgsPointXY > &destinationCoordinates)
Creates a new QgsGcpTransformerInterface subclass representing the specified transform 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:60
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.
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.