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
transform(double & x,double & y,bool inverseTransform) const29 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
methodToString(QgsGcpTransformerInterface::TransformMethod method)47 QString QgsGcpTransformerInterface::methodToString( QgsGcpTransformerInterface::TransformMethod method )
48 {
49 switch ( method )
50 {
51 case QgsGcpTransformerInterface::TransformMethod::Linear:
52 return QObject::tr( "Linear" );
53 case QgsGcpTransformerInterface::TransformMethod::Helmert:
54 return QObject::tr( "Helmert" );
55 case QgsGcpTransformerInterface::TransformMethod::PolynomialOrder1:
56 return QObject::tr( "Polynomial 1" );
57 case QgsGcpTransformerInterface::TransformMethod::PolynomialOrder2:
58 return QObject::tr( "Polynomial 2" );
59 case QgsGcpTransformerInterface::TransformMethod::PolynomialOrder3:
60 return QObject::tr( "Polynomial 3" );
61 case QgsGcpTransformerInterface::TransformMethod::ThinPlateSpline:
62 return QObject::tr( "Thin Plate Spline (TPS)" );
63 case QgsGcpTransformerInterface::TransformMethod::Projective:
64 return QObject::tr( "Projective" );
65 default:
66 return QObject::tr( "Not set" );
67 }
68 }
69
create(QgsGcpTransformerInterface::TransformMethod method)70 QgsGcpTransformerInterface *QgsGcpTransformerInterface::create( QgsGcpTransformerInterface::TransformMethod method )
71 {
72 switch ( method )
73 {
74 case TransformMethod::Linear:
75 return new QgsLinearGeorefTransform;
76 case TransformMethod::Helmert:
77 return new QgsHelmertGeorefTransform;
78 case TransformMethod::PolynomialOrder1:
79 return new QgsGDALGeorefTransform( false, 1 );
80 case TransformMethod::PolynomialOrder2:
81 return new QgsGDALGeorefTransform( false, 2 );
82 case TransformMethod::PolynomialOrder3:
83 return new QgsGDALGeorefTransform( false, 3 );
84 case TransformMethod::ThinPlateSpline:
85 return new QgsGDALGeorefTransform( true, 0 );
86 case TransformMethod::Projective:
87 return new QgsProjectiveGeorefTransform;
88 default:
89 return nullptr;
90 }
91 }
92
createFromParameters(QgsGcpTransformerInterface::TransformMethod method,const QVector<QgsPointXY> & sourceCoordinates,const QVector<QgsPointXY> & destinationCoordinates)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
getOriginScale(QgsPointXY & origin,double & scaleX,double & scaleY) const110 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
clone() const118 QgsGcpTransformerInterface *QgsLinearGeorefTransform::clone() const
119 {
120 std::unique_ptr< QgsLinearGeorefTransform > res = std::make_unique< QgsLinearGeorefTransform >();
121 res->mParameters = mParameters;
122 return res.release();
123 }
124
updateParametersFromGcps(const QVector<QgsPointXY> & sourceCoordinates,const QVector<QgsPointXY> & destinationCoordinates,bool invertYAxis)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
minimumGcpCount() const135 int QgsLinearGeorefTransform::minimumGcpCount() const
136 {
137 return 2;
138 }
139
GDALTransformer() const140 GDALTransformerFunc QgsLinearGeorefTransform::GDALTransformer() const
141 {
142 return QgsLinearGeorefTransform::linearTransform;
143 }
144
GDALTransformerArgs() const145 void *QgsLinearGeorefTransform::GDALTransformerArgs() const
146 {
147 return ( void * )&mParameters;
148 }
149
method() const150 QgsGcpTransformerInterface::TransformMethod QgsLinearGeorefTransform::method() const
151 {
152 return TransformMethod::Linear;
153 }
154
linearTransform(void * pTransformerArg,int bDstToSrc,int nPointCount,double * x,double * y,double * z,int * panSuccess)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 //
updateParametersFromGcps(const QVector<QgsPointXY> & sourceCoordinates,const QVector<QgsPointXY> & destinationCoordinates,bool invertYAxis)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
minimumGcpCount() const208 int QgsHelmertGeorefTransform::minimumGcpCount() const
209 {
210 return 2;
211 }
212
GDALTransformer() const213 GDALTransformerFunc QgsHelmertGeorefTransform::GDALTransformer() const
214 {
215 return QgsHelmertGeorefTransform::helmertTransform;
216 }
217
GDALTransformerArgs() const218 void *QgsHelmertGeorefTransform::GDALTransformerArgs() const
219 {
220 return ( void * )&mHelmertParameters;
221 }
222
method() const223 QgsGcpTransformerInterface::TransformMethod QgsHelmertGeorefTransform::method() const
224 {
225 return TransformMethod::Helmert;
226 }
227
getOriginScaleRotation(QgsPointXY & origin,double & scale,double & rotation) const228 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
clone() const236 QgsGcpTransformerInterface *QgsHelmertGeorefTransform::clone() const
237 {
238 std::unique_ptr< QgsHelmertGeorefTransform > res = std::make_unique< QgsHelmertGeorefTransform >();
239 res->mHelmertParameters = mHelmertParameters;
240 return res.release();
241 }
242
helmertTransform(void * pTransformerArg,int bDstToSrc,int nPointCount,double * x,double * y,double * z,int * panSuccess)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
QgsGDALGeorefTransform(bool useTPS,unsigned int polynomialOrder)321 QgsGDALGeorefTransform::QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder )
322 : mPolynomialOrder( std::min( 3u, polynomialOrder ) )
323 , mIsTPSTransform( useTPS )
324 {
325 }
326
~QgsGDALGeorefTransform()327 QgsGDALGeorefTransform::~QgsGDALGeorefTransform()
328 {
329 destroyGdalArgs();
330 }
331
clone() const332 QgsGcpTransformerInterface *QgsGDALGeorefTransform::clone() const
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
updateParametersFromGcps(const QVector<QgsPointXY> & sourceCoordinates,const QVector<QgsPointXY> & destinationCoordinates,bool invertYAxis)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 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
minimumGcpCount() const378 int QgsGDALGeorefTransform::minimumGcpCount() const
379 {
380 if ( mIsTPSTransform )
381 return 1;
382 else
383 return ( ( mPolynomialOrder + 2 ) * ( mPolynomialOrder + 1 ) ) / 2;
384 }
385
GDALTransformer() const386 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
GDALTransformerArgs() const398 void *QgsGDALGeorefTransform::GDALTransformerArgs() const
399 {
400 return mGDALTransformerArgs;
401 }
402
method() const403 QgsGcpTransformerInterface::TransformMethod QgsGDALGeorefTransform::method() const
404 {
405 if ( mIsTPSTransform )
406 return TransformMethod::ThinPlateSpline;
407
408 switch ( mPolynomialOrder )
409 {
410 case 1:
411 return TransformMethod::PolynomialOrder1;
412 case 2:
413 return TransformMethod::PolynomialOrder2;
414 case 3:
415 return TransformMethod::PolynomialOrder3;
416 }
417 return TransformMethod::InvalidTransform;
418 }
419
destroyGdalArgs()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
QgsProjectiveGeorefTransform()435 QgsProjectiveGeorefTransform::QgsProjectiveGeorefTransform()
436 : mParameters()
437 {}
438
clone() const439 QgsGcpTransformerInterface *QgsProjectiveGeorefTransform::clone() const
440 {
441 std::unique_ptr< QgsProjectiveGeorefTransform > res = std::make_unique< QgsProjectiveGeorefTransform >();
442 res->mParameters = mParameters;
443 return res.release();
444 }
445
updateParametersFromGcps(const QVector<QgsPointXY> & sourceCoordinates,const QVector<QgsPointXY> & destinationCoordinates,bool invertYAxis)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 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
minimumGcpCount() const502 int QgsProjectiveGeorefTransform::minimumGcpCount() const
503 {
504 return 4;
505 }
506
GDALTransformer() const507 GDALTransformerFunc QgsProjectiveGeorefTransform::GDALTransformer() const
508 {
509 return QgsProjectiveGeorefTransform::projectiveTransform;
510 }
511
GDALTransformerArgs() const512 void *QgsProjectiveGeorefTransform::GDALTransformerArgs() const
513 {
514 return ( void * )&mParameters;
515 }
516
method() const517 QgsGcpTransformerInterface::TransformMethod QgsProjectiveGeorefTransform::method() const
518 {
519 return TransformMethod::Projective;
520 }
521
projectiveTransform(void * pTransformerArg,int bDstToSrc,int nPointCount,double * x,double * y,double * z,int * panSuccess)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 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 }
569