30void QgsLeastSquares::linear(
const QVector<QgsPointXY> &sourceCoordinates,
const QVector<QgsPointXY> &destinationCoordinates,
QgsPointXY &origin,
double &pixelXSize,
double &pixelYSize )
32 const int n = destinationCoordinates.size();
35 throw std::domain_error( QObject::tr(
"Fit to a linear transform requires at least 2 points." ).toLocal8Bit().constData() );
38 double sumPx( 0 ), sumPy( 0 ), sumPx2( 0 ), sumPy2( 0 ), sumPxMx( 0 ), sumPyMy( 0 ), sumMx( 0 ), sumMy( 0 );
39 for (
int i = 0; i < n; ++i )
41 sumPx += sourceCoordinates.at( i ).x();
42 sumPy += sourceCoordinates.at( i ).y();
43 sumPx2 += std::pow( sourceCoordinates.at( i ).x(), 2 );
44 sumPy2 += std::pow( sourceCoordinates.at( i ).y(), 2 );
45 sumPxMx += sourceCoordinates.at( i ).x() * destinationCoordinates.at( i ).x();
46 sumPyMy += sourceCoordinates.at( i ).y() * destinationCoordinates.at( i ).y();
47 sumMx += destinationCoordinates.at( i ).x();
48 sumMy += destinationCoordinates.at( i ).y();
51 const double deltaX = n * sumPx2 - std::pow( sumPx, 2 );
52 const double deltaY = n * sumPy2 - std::pow( sumPy, 2 );
54 const double aX = ( sumPx2 * sumMx - sumPx * sumPxMx ) / deltaX;
55 const double aY = ( sumPy2 * sumMy - sumPy * sumPyMy ) / deltaY;
56 const double bX = ( n * sumPxMx - sumPx * sumMx ) / deltaX;
57 const double bY = ( n * sumPyMy - sumPy * sumMy ) / deltaY;
62 pixelXSize = std::fabs( bX );
63 pixelYSize = std::fabs( bY );
67void QgsLeastSquares::helmert(
const QVector<QgsPointXY> &sourceCoordinates,
const QVector<QgsPointXY> &destinationCoordinates,
QgsPointXY &origin,
double &pixelSize,
double &rotation )
70 ( void ) sourceCoordinates;
71 ( void ) destinationCoordinates;
75 throw QgsNotSupportedException( QObject::tr(
"Calculating a helmert transformation requires a QGIS build based GSL" ) );
77 const int n = destinationCoordinates.size();
80 throw std::domain_error( QObject::tr(
"Fit to a Helmert transform requires at least 2 points." ).toLocal8Bit().constData() );
93 for (
int i = 0; i < n; ++i )
95 A += sourceCoordinates.at( i ).x();
96 B += sourceCoordinates.at( i ).y();
97 C += destinationCoordinates.at( i ).x();
98 D += destinationCoordinates.at( i ).y();
99 E += destinationCoordinates.at( i ).x() * sourceCoordinates.at( i ).x();
100 F += destinationCoordinates.at( i ).y() * sourceCoordinates.at( i ).y();
101 G += std::pow( sourceCoordinates.at( i ).x(), 2 );
102 H += std::pow( sourceCoordinates.at( i ).y(), 2 );
103 I += destinationCoordinates.at( i ).x() * sourceCoordinates.at( i ).y();
104 J += sourceCoordinates.at( i ).x() * destinationCoordinates.at( i ).y();
112 double MData[] = { A, -B, ( double ) n, 0., B, A, 0., (
double ) n, G + H, 0., A, B, 0., G + H, -B, A };
114 double bData[] = { C, D, E + F, J - I };
117 gsl_matrix_view M = gsl_matrix_view_array( MData, 4, 4 );
118 const gsl_vector_view b = gsl_vector_view_array( bData, 4 );
119 gsl_vector *x = gsl_vector_alloc( 4 );
120 gsl_permutation *p = gsl_permutation_alloc( 4 );
122 gsl_linalg_LU_decomp( &M.matrix, p, &s );
123 gsl_linalg_LU_solve( &M.matrix, p, &b.vector, x );
124 gsl_permutation_free( p );
126 origin.
setX( gsl_vector_get( x, 2 ) );
127 origin.
setY( gsl_vector_get( x, 3 ) );
128 pixelSize = std::sqrt( std::pow( gsl_vector_get( x, 0 ), 2 ) + std::pow( gsl_vector_get( x, 1 ), 2 ) );
129 rotation = std::atan2( gsl_vector_get( x, 1 ), gsl_vector_get( x, 0 ) );
131 gsl_vector_free( x );
195void normalizeCoordinates(
const QVector<QgsPointXY> &coords, QVector<QgsPointXY> &normalizedCoords,
double normalizeMatrix[9],
double denormalizeMatrix[9] )
198 double cogX = 0.0, cogY = 0.0;
199 for (
int i = 0; i < coords.size(); i++ )
201 cogX += coords[i].x();
202 cogY += coords[i].y();
204 cogX *= 1.0 / coords.size();
205 cogY *= 1.0 / coords.size();
208 double meanDist = 0.0;
209 for (
int i = 0; i < coords.size(); i++ )
211 const double X = ( coords[i].x() - cogX );
212 const double Y = ( coords[i].y() - cogY );
213 meanDist += std::sqrt( X * X + Y * Y );
215 meanDist *= 1.0 / coords.size();
217 const double OOD = meanDist * M_SQRT1_2;
218 const double D = 1.0 / OOD;
219 normalizedCoords.resize( coords.size() );
220 for (
int i = 0; i < coords.size(); i++ )
222 normalizedCoords[i] =
QgsPointXY( ( coords[i].x() - cogX ) * D, ( coords[i].y() - cogY ) * D );
225 normalizeMatrix[0] = D;
226 normalizeMatrix[1] = 0.0;
227 normalizeMatrix[2] = -cogX * D;
228 normalizeMatrix[3] = 0.0;
229 normalizeMatrix[4] = D;
230 normalizeMatrix[5] = -cogY * D;
231 normalizeMatrix[6] = 0.0;
232 normalizeMatrix[7] = 0.0;
233 normalizeMatrix[8] = 1.0;
235 denormalizeMatrix[0] = OOD;
236 denormalizeMatrix[1] = 0.0;
237 denormalizeMatrix[2] = cogX;
238 denormalizeMatrix[3] = 0.0;
239 denormalizeMatrix[4] = OOD;
240 denormalizeMatrix[5] = cogY;
241 denormalizeMatrix[6] = 0.0;
242 denormalizeMatrix[7] = 0.0;
243 denormalizeMatrix[8] = 1.0;
251 ( void ) sourceCoordinates;
252 ( void ) destinationCoordinates;
254 throw QgsNotSupportedException( QObject::tr(
"Calculating a projective transformation requires a QGIS build based GSL" ) );
256 Q_ASSERT( sourceCoordinates.size() == destinationCoordinates.size() );
258 if ( destinationCoordinates.size() < 4 )
260 throw std::domain_error( QObject::tr(
"Fitting a projective transform requires at least 4 corresponding points." ).toLocal8Bit().constData() );
263 QVector<QgsPointXY> sourceCoordinatesNormalized;
264 QVector<QgsPointXY> destinationCoordinatesNormalized;
266 double normSource[9], denormSource[9];
267 double normDest[9], denormDest[9];
268 normalizeCoordinates( sourceCoordinates, sourceCoordinatesNormalized, normSource, denormSource );
269 normalizeCoordinates( destinationCoordinates, destinationCoordinatesNormalized, normDest, denormDest );
273 const uint m = std::max( 9u, ( uint ) destinationCoordinatesNormalized.size() * 2u );
275 gsl_matrix *S = gsl_matrix_alloc( m, n );
277 for (
int i = 0; i < destinationCoordinatesNormalized.size(); i++ )
279 gsl_matrix_set( S, i * 2, 0, sourceCoordinatesNormalized[i].x() );
280 gsl_matrix_set( S, i * 2, 1, sourceCoordinatesNormalized[i].y() );
281 gsl_matrix_set( S, i * 2, 2, 1.0 );
283 gsl_matrix_set( S, i * 2, 3, 0.0 );
284 gsl_matrix_set( S, i * 2, 4, 0.0 );
285 gsl_matrix_set( S, i * 2, 5, 0.0 );
287 gsl_matrix_set( S, i * 2, 6, -destinationCoordinatesNormalized[i].x() * sourceCoordinatesNormalized[i].x() );
288 gsl_matrix_set( S, i * 2, 7, -destinationCoordinatesNormalized[i].x() * sourceCoordinatesNormalized[i].y() );
289 gsl_matrix_set( S, i * 2, 8, -destinationCoordinatesNormalized[i].x() * 1.0 );
291 gsl_matrix_set( S, i * 2 + 1, 0, 0.0 );
292 gsl_matrix_set( S, i * 2 + 1, 1, 0.0 );
293 gsl_matrix_set( S, i * 2 + 1, 2, 0.0 );
295 gsl_matrix_set( S, i * 2 + 1, 3, sourceCoordinatesNormalized[i].x() );
296 gsl_matrix_set( S, i * 2 + 1, 4, sourceCoordinatesNormalized[i].y() );
297 gsl_matrix_set( S, i * 2 + 1, 5, 1.0 );
299 gsl_matrix_set( S, i * 2 + 1, 6, -destinationCoordinatesNormalized[i].y() * sourceCoordinatesNormalized[i].x() );
300 gsl_matrix_set( S, i * 2 + 1, 7, -destinationCoordinatesNormalized[i].y() * sourceCoordinatesNormalized[i].y() );
301 gsl_matrix_set( S, i * 2 + 1, 8, -destinationCoordinatesNormalized[i].y() * 1.0 );
304 if ( destinationCoordinatesNormalized.size() == 4 )
312 for (
int j = 0; j < 9; j++ )
314 gsl_matrix_set( S, 8, j, gsl_matrix_get( S, 7, j ) );
322 gsl_matrix *V = gsl_matrix_alloc( n, n );
323 gsl_vector *singular_values = gsl_vector_alloc( n );
324 gsl_vector *work = gsl_vector_alloc( n );
328 gsl_linalg_SV_decomp( S, V, singular_values, work );
331 for (
unsigned int i = 0; i < n; i++ )
333 H[i] = gsl_matrix_get( V, i, n - 1 );
336 gsl_matrix *prodMatrix = gsl_matrix_alloc( 3, 3 );
338 gsl_matrix_view Hmatrix = gsl_matrix_view_array( H, 3, 3 );
339 const gsl_matrix_view normSourceMatrix = gsl_matrix_view_array( normSource, 3, 3 );
340 const gsl_matrix_view denormDestMatrix = gsl_matrix_view_array( denormDest, 3, 3 );
344 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &Hmatrix.matrix, &normSourceMatrix.matrix, 0.0, prodMatrix );
345 gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &denormDestMatrix.matrix, prodMatrix, 0.0, &Hmatrix.matrix );
347 gsl_matrix_free( prodMatrix );
348 gsl_matrix_free( S );
349 gsl_matrix_free( V );
350 gsl_vector_free( singular_values );
351 gsl_vector_free( work );