mirror of
https://github.com/qgis/QGIS.git
synced 2025-10-05 00:09:32 -04:00
Rename pixel/layer/map coordinates to "source" and "destination", to
make classes more generic and avoid confusion with non-raster based GCPs Also flip function arguments to source, destination order instead of destination, source
This commit is contained in:
Nyall Dawson
parent
30fac0375e
commit
ce54321220
@ -37,6 +37,9 @@ based on a transformation method and a list of GCPs.
|
||||
};
|
||||
|
||||
QgsGcpTransformerInterface();
|
||||
%Docstring
|
||||
Constructor for QgsGcpTransformerInterface
|
||||
%End
|
||||
|
||||
virtual ~QgsGcpTransformerInterface();
|
||||
|
||||
@ -50,9 +53,12 @@ parameters as this one.
|
||||
Caller takes ownership of the returned object.
|
||||
%End
|
||||
|
||||
virtual bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoordinates, const QVector<QgsPointXY> &layerCoordinates ) = 0;
|
||||
virtual bool updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis = false ) = 0;
|
||||
%Docstring
|
||||
Fits transformation parameters using the specified Ground Control Points (GCPs) lists of map coordinates and layer coordinates.
|
||||
Fits transformation parameters using the specified Ground Control Points (GCPs) lists of source and destination coordinates.
|
||||
|
||||
If ``invertYAxis`` is set to ``True`` then the y-axis of source coordinates will be inverted, e.g. to allow for transformation of raster layers
|
||||
with ascending top-to-bottom vertical axis coordinates.
|
||||
|
||||
:return: ``True`` on success, ``False`` on failure
|
||||
%End
|
||||
@ -88,10 +94,10 @@ Creates a new QgsGcpTransformerInterface subclass representing the specified tra
|
||||
Caller takes ownership of the returned object.
|
||||
%End
|
||||
|
||||
static QgsGcpTransformerInterface *createFromParameters( TransformMethod method, const QVector<QgsPointXY> &mapCoordinates, const QVector<QgsPointXY> &layerCoordinates ) /Factory/;
|
||||
static QgsGcpTransformerInterface *createFromParameters( TransformMethod method, const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates ) /Factory/;
|
||||
%Docstring
|
||||
Creates a new QgsGcpTransformerInterface subclass representing the specified transform ``method``, initialized
|
||||
using the given lists of map and layer coordinates.
|
||||
using the given lists of source and destination coordinates.
|
||||
|
||||
If the parameters cannot be fit to a transform ``None`` will be returned.
|
||||
|
||||
|
||||
@ -90,13 +90,13 @@ QgsGcpTransformerInterface *QgsGcpTransformerInterface::create( QgsGcpTransforme
|
||||
}
|
||||
}
|
||||
|
||||
QgsGcpTransformerInterface *QgsGcpTransformerInterface::createFromParameters( QgsGcpTransformerInterface::TransformMethod method, const QVector<QgsPointXY> &mapCoordinates, const QVector<QgsPointXY> &layerCoordinates )
|
||||
QgsGcpTransformerInterface *QgsGcpTransformerInterface::createFromParameters( QgsGcpTransformerInterface::TransformMethod method, const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates )
|
||||
{
|
||||
std::unique_ptr< QgsGcpTransformerInterface > transformer( create( method ) );
|
||||
if ( !transformer )
|
||||
return nullptr;
|
||||
|
||||
if ( !transformer->updateParametersFromGcps( mapCoordinates, layerCoordinates ) )
|
||||
if ( !transformer->updateParametersFromGcps( sourceCoordinates, destinationCoordinates ) )
|
||||
return nullptr;
|
||||
|
||||
return transformer.release();
|
||||
@ -122,11 +122,13 @@ QgsGcpTransformerInterface *QgsLinearGeorefTransform::clone() const
|
||||
return res.release();
|
||||
}
|
||||
|
||||
bool QgsLinearGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &layerCoords )
|
||||
bool QgsLinearGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
|
||||
{
|
||||
if ( mapCoords.size() < minimumGcpCount() )
|
||||
if ( destinationCoordinates.size() < minimumGcpCount() )
|
||||
return false;
|
||||
QgsLeastSquares::linear( mapCoords, layerCoords, mParameters.origin, mParameters.scaleX, mParameters.scaleY );
|
||||
|
||||
mParameters.invertYAxis = invertYAxis;
|
||||
QgsLeastSquares::linear( sourceCoordinates, destinationCoordinates, mParameters.origin, mParameters.scaleX, mParameters.scaleY );
|
||||
return true;
|
||||
}
|
||||
|
||||
@ -163,7 +165,7 @@ int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstTo
|
||||
for ( int i = 0; i < nPointCount; ++i )
|
||||
{
|
||||
x[i] = x[i] * t->scaleX + t->origin.x();
|
||||
y[i] = -y[i] * t->scaleY + t->origin.y();
|
||||
y[i] = ( t->invertYAxis ? -1 : 1 ) * y[i] * t->scaleY + t->origin.y();
|
||||
panSuccess[i] = true;
|
||||
}
|
||||
}
|
||||
@ -182,7 +184,7 @@ int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstTo
|
||||
for ( int i = 0; i < nPointCount; ++i )
|
||||
{
|
||||
x[i] = ( x[i] - t->origin.x() ) / t->scaleX;
|
||||
y[i] = ( y[i] - t->origin.y() ) / ( -t->scaleY );
|
||||
y[i] = ( y[i] - t->origin.y() ) / ( ( t->invertYAxis ? -1 : 1 ) * t->scaleY );
|
||||
panSuccess[i] = true;
|
||||
}
|
||||
}
|
||||
@ -193,12 +195,13 @@ int QgsLinearGeorefTransform::linearTransform( void *pTransformerArg, int bDstTo
|
||||
//
|
||||
// QgsHelmertGeorefTransform
|
||||
//
|
||||
bool QgsHelmertGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &layerCoords )
|
||||
bool QgsHelmertGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
|
||||
{
|
||||
if ( mapCoords.size() < minimumGcpCount() )
|
||||
if ( destinationCoordinates.size() < minimumGcpCount() )
|
||||
return false;
|
||||
|
||||
QgsLeastSquares::helmert( mapCoords, layerCoords, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle );
|
||||
mHelmertParameters.invertYAxis = invertYAxis;
|
||||
QgsLeastSquares::helmert( sourceCoordinates, destinationCoordinates, mHelmertParameters.origin, mHelmertParameters.scale, mHelmertParameters.angle );
|
||||
return true;
|
||||
}
|
||||
|
||||
@ -207,10 +210,9 @@ int QgsHelmertGeorefTransform::minimumGcpCount() const
|
||||
return 2;
|
||||
}
|
||||
|
||||
|
||||
GDALTransformerFunc QgsHelmertGeorefTransform::GDALTransformer() const
|
||||
{
|
||||
return QgsHelmertGeorefTransform::helmert_transform;
|
||||
return QgsHelmertGeorefTransform::helmertTransform;
|
||||
}
|
||||
|
||||
void *QgsHelmertGeorefTransform::GDALTransformerArgs() const
|
||||
@ -238,28 +240,42 @@ QgsGcpTransformerInterface *QgsHelmertGeorefTransform::clone() const
|
||||
return res.release();
|
||||
}
|
||||
|
||||
int QgsHelmertGeorefTransform::helmert_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
|
||||
int QgsHelmertGeorefTransform::helmertTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
|
||||
double *x, double *y, double *z, int *panSuccess )
|
||||
{
|
||||
Q_UNUSED( z )
|
||||
HelmertParameters *t = static_cast<HelmertParameters *>( pTransformerArg );
|
||||
const HelmertParameters *t = static_cast< const HelmertParameters *>( pTransformerArg );
|
||||
if ( !t )
|
||||
return false;
|
||||
|
||||
double a = std::cos( t->angle ), b = std::sin( t->angle ), x0 = t->origin.x(), y0 = t->origin.y(), s = t->scale;
|
||||
double a = std::cos( t->angle );
|
||||
double b = std::sin( t->angle );
|
||||
const double x0 = t->origin.x();
|
||||
const double y0 = t->origin.y();
|
||||
const double s = t->scale;
|
||||
if ( !bDstToSrc )
|
||||
{
|
||||
a *= s;
|
||||
b *= s;
|
||||
for ( int i = 0; i < nPointCount; ++i )
|
||||
{
|
||||
double xT = x[i], yT = y[i];
|
||||
// Because rotation parameters where estimated in a CS with negative y-axis ^= down.
|
||||
// we need to apply the rotation matrix and a change of base:
|
||||
// |cos a,-sin a| |1, 0| | std::cos a, std::sin a|
|
||||
// |sin a, std::cos a| |0,-1| = | std::sin a, -cos a|
|
||||
x[i] = x0 + ( a * xT + b * yT );
|
||||
y[i] = y0 + ( b * xT - a * yT );
|
||||
const double xT = x[i];
|
||||
const double yT = y[i];
|
||||
|
||||
if ( t->invertYAxis )
|
||||
{
|
||||
// Because rotation parameters where estimated in a CS with negative y-axis ^= down.
|
||||
// we need to apply the rotation matrix and a change of base:
|
||||
// |cos a, -sin a | |1, 0| | cos a, sin a|
|
||||
// |sin a, cos a | |0,-1| = | sin a, -cos a|
|
||||
x[i] = x0 + ( a * xT + b * yT );
|
||||
y[i] = y0 + ( b * xT - a * yT );
|
||||
}
|
||||
else
|
||||
{
|
||||
x[i] = x0 + ( a * xT - b * yT );
|
||||
y[i] = y0 + ( b * xT + a * yT );
|
||||
}
|
||||
panSuccess[i] = true;
|
||||
}
|
||||
}
|
||||
@ -278,13 +294,20 @@ int QgsHelmertGeorefTransform::helmert_transform( void *pTransformerArg, int bDs
|
||||
b /= s;
|
||||
for ( int i = 0; i < nPointCount; ++i )
|
||||
{
|
||||
double xT = x[i], yT = y[i];
|
||||
xT -= x0;
|
||||
yT -= y0;
|
||||
// | std::cos a, std::sin a |^-1 |cos a, std::sin a|
|
||||
// | std::sin a, -cos a | = |sin a, -cos a|
|
||||
x[i] = a * xT + b * yT;
|
||||
y[i] = b * xT - a * yT;
|
||||
const double xT = x[i] - x0;
|
||||
const double yT = y[i] - y0;
|
||||
if ( t->invertYAxis )
|
||||
{
|
||||
// | cos a, sin a |^-1 |cos a, sin a|
|
||||
// | sin a, -cos a | = |sin a, -cos a|
|
||||
x[i] = a * xT + b * yT;
|
||||
y[i] = b * xT - a * yT;
|
||||
}
|
||||
else
|
||||
{
|
||||
x[i] = a * xT + b * yT;
|
||||
y[i] = -b * xT + a * yT;
|
||||
}
|
||||
panSuccess[i] = true;
|
||||
}
|
||||
}
|
||||
@ -309,19 +332,20 @@ QgsGDALGeorefTransform::~QgsGDALGeorefTransform()
|
||||
QgsGcpTransformerInterface *QgsGDALGeorefTransform::clone() const
|
||||
{
|
||||
std::unique_ptr< QgsGDALGeorefTransform > res = qgis::make_unique< QgsGDALGeorefTransform >( mIsTPSTransform, mPolynomialOrder );
|
||||
res->updateParametersFromGcps( mMapCoords, mLayerCoords );
|
||||
res->updateParametersFromGcps( mSourceCoords, mDestCoordinates, mInvertYAxis );
|
||||
return res.release();
|
||||
}
|
||||
|
||||
bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &layerCoords )
|
||||
bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
|
||||
{
|
||||
mMapCoords = mapCoords;
|
||||
mLayerCoords = layerCoords;
|
||||
mSourceCoords = sourceCoordinates;
|
||||
mDestCoordinates = destinationCoordinates;
|
||||
mInvertYAxis = invertYAxis;
|
||||
|
||||
assert( mapCoords.size() == layerCoords.size() );
|
||||
if ( mapCoords.size() != layerCoords.size() )
|
||||
assert( sourceCoordinates.size() == destinationCoordinates.size() );
|
||||
if ( sourceCoordinates.size() != destinationCoordinates.size() )
|
||||
return false;
|
||||
int n = mapCoords.size();
|
||||
int n = sourceCoordinates.size();
|
||||
|
||||
GDAL_GCP *GCPList = new GDAL_GCP[n];
|
||||
for ( int i = 0; i < n; i++ )
|
||||
@ -329,10 +353,10 @@ bool QgsGDALGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY>
|
||||
GCPList[i].pszId = new char[20];
|
||||
snprintf( GCPList[i].pszId, 19, "gcp%i", i );
|
||||
GCPList[i].pszInfo = nullptr;
|
||||
GCPList[i].dfGCPPixel = layerCoords[i].x();
|
||||
GCPList[i].dfGCPLine = -layerCoords[i].y();
|
||||
GCPList[i].dfGCPX = mapCoords[i].x();
|
||||
GCPList[i].dfGCPY = mapCoords[i].y();
|
||||
GCPList[i].dfGCPPixel = sourceCoordinates[i].x();
|
||||
GCPList[i].dfGCPLine = ( mInvertYAxis ? -1 : 1 ) * sourceCoordinates[i].y();
|
||||
GCPList[i].dfGCPX = destinationCoordinates[i].x();
|
||||
GCPList[i].dfGCPY = destinationCoordinates[i].y();
|
||||
GCPList[i].dfGCPZ = 0;
|
||||
}
|
||||
destroyGdalArgs();
|
||||
@ -419,20 +443,27 @@ QgsGcpTransformerInterface *QgsProjectiveGeorefTransform::clone() const
|
||||
return res.release();
|
||||
}
|
||||
|
||||
bool QgsProjectiveGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &layerCoords )
|
||||
bool QgsProjectiveGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
|
||||
{
|
||||
if ( mapCoords.size() < minimumGcpCount() )
|
||||
if ( destinationCoordinates.size() < minimumGcpCount() )
|
||||
return false;
|
||||
|
||||
// HACK: flip y coordinates, because georeferencer and gdal use different conventions
|
||||
QVector<QgsPointXY> flippedPixelCoords;
|
||||
flippedPixelCoords.reserve( layerCoords.size() );
|
||||
for ( const QgsPointXY &coord : layerCoords )
|
||||
if ( invertYAxis )
|
||||
{
|
||||
flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() );
|
||||
}
|
||||
// HACK: flip y coordinates, because georeferencer and gdal use different conventions
|
||||
QVector<QgsPointXY> flippedPixelCoords;
|
||||
flippedPixelCoords.reserve( sourceCoordinates.size() );
|
||||
for ( const QgsPointXY &coord : sourceCoordinates )
|
||||
{
|
||||
flippedPixelCoords << QgsPointXY( coord.x(), -coord.y() );
|
||||
}
|
||||
|
||||
QgsLeastSquares::projective( mapCoords, flippedPixelCoords, mParameters.H );
|
||||
QgsLeastSquares::projective( flippedPixelCoords, destinationCoordinates, mParameters.H );
|
||||
}
|
||||
else
|
||||
{
|
||||
QgsLeastSquares::projective( sourceCoordinates, destinationCoordinates, mParameters.H );
|
||||
}
|
||||
|
||||
// Invert the homography matrix using adjoint matrix
|
||||
double *H = mParameters.H;
|
||||
|
||||
@ -49,6 +49,7 @@ class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT
|
||||
InvalidTransform = 65535 //!< Invalid transform
|
||||
};
|
||||
|
||||
//! Constructor for QgsGcpTransformerInterface
|
||||
QgsGcpTransformerInterface() = default;
|
||||
|
||||
virtual ~QgsGcpTransformerInterface() = default;
|
||||
@ -68,11 +69,14 @@ class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT
|
||||
virtual QgsGcpTransformerInterface *clone() const = 0 SIP_FACTORY;
|
||||
|
||||
/**
|
||||
* Fits transformation parameters using the specified Ground Control Points (GCPs) lists of map coordinates and layer coordinates.
|
||||
* Fits transformation parameters using the specified Ground Control Points (GCPs) lists of source and destination coordinates.
|
||||
*
|
||||
* If \a invertYAxis is set to TRUE then the y-axis of source coordinates will be inverted, e.g. to allow for transformation of raster layers
|
||||
* with ascending top-to-bottom vertical axis coordinates.
|
||||
*
|
||||
* \returns TRUE on success, FALSE on failure
|
||||
*/
|
||||
virtual bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoordinates, const QVector<QgsPointXY> &layerCoordinates ) = 0;
|
||||
virtual bool updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis = false ) = 0;
|
||||
|
||||
/**
|
||||
* Returns the minimum number of Ground Control Points (GCPs) required for parameter fitting.
|
||||
@ -107,13 +111,13 @@ class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT
|
||||
|
||||
/**
|
||||
* Creates a new QgsGcpTransformerInterface subclass representing the specified transform \a method, initialized
|
||||
* using the given lists of map and layer coordinates.
|
||||
* using the given lists of source and destination coordinates.
|
||||
*
|
||||
* If the parameters cannot be fit to a transform NULLPTR will be returned.
|
||||
*
|
||||
* Caller takes ownership of the returned object.
|
||||
*/
|
||||
static QgsGcpTransformerInterface *createFromParameters( TransformMethod method, const QVector<QgsPointXY> &mapCoordinates, const QVector<QgsPointXY> &layerCoordinates ) SIP_FACTORY;
|
||||
static QgsGcpTransformerInterface *createFromParameters( TransformMethod method, const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates ) SIP_FACTORY;
|
||||
|
||||
#ifndef SIP_RUN
|
||||
|
||||
@ -144,6 +148,8 @@ class ANALYSIS_EXPORT QgsGcpTransformerInterface SIP_ABSTRACT
|
||||
class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP
|
||||
{
|
||||
public:
|
||||
|
||||
//! Constructor for QgsLinearGeorefTransform
|
||||
QgsLinearGeorefTransform() = default;
|
||||
|
||||
/**
|
||||
@ -152,7 +158,7 @@ class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterfa
|
||||
bool getOriginScale( QgsPointXY &origin, double &scaleX, double &scaleY ) const;
|
||||
|
||||
QgsGcpTransformerInterface *clone() const override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &layerCoords ) override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis = false ) override;
|
||||
int minimumGcpCount() const override;
|
||||
GDALTransformerFunc GDALTransformer() const override;
|
||||
void *GDALTransformerArgs() const override;
|
||||
@ -162,7 +168,9 @@ class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterfa
|
||||
struct LinearParameters
|
||||
{
|
||||
QgsPointXY origin;
|
||||
double scaleX, scaleY;
|
||||
double scaleX = 1;
|
||||
double scaleY = 1;
|
||||
bool invertYAxis = false;
|
||||
} mParameters;
|
||||
|
||||
static int linearTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
|
||||
@ -179,6 +187,8 @@ class ANALYSIS_EXPORT QgsLinearGeorefTransform : public QgsGcpTransformerInterfa
|
||||
class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP
|
||||
{
|
||||
public:
|
||||
|
||||
//! Constructor for QgsHelmertGeorefTransform
|
||||
QgsHelmertGeorefTransform() = default;
|
||||
|
||||
/**
|
||||
@ -187,7 +197,7 @@ class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterf
|
||||
bool getOriginScaleRotation( QgsPointXY &origin, double &scale, double &rotation ) const;
|
||||
|
||||
QgsGcpTransformerInterface *clone() const override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &layerCoords ) override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis = false ) override;
|
||||
int minimumGcpCount() const override;
|
||||
GDALTransformerFunc GDALTransformer() const override;
|
||||
void *GDALTransformerArgs() const override;
|
||||
@ -200,11 +210,12 @@ class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterf
|
||||
QgsPointXY origin;
|
||||
double scale;
|
||||
double angle;
|
||||
bool invertYAxis = false;
|
||||
};
|
||||
HelmertParameters mHelmertParameters;
|
||||
|
||||
static int helmert_transform( void *pTransformerArg, int bDstToSrc, int nPointCount,
|
||||
double *x, double *y, double *z, int *panSuccess );
|
||||
static int helmertTransform( void *pTransformerArg, int bDstToSrc, int nPointCount,
|
||||
double *x, double *y, double *z, int *panSuccess );
|
||||
|
||||
};
|
||||
|
||||
@ -217,11 +228,13 @@ class ANALYSIS_EXPORT QgsHelmertGeorefTransform : public QgsGcpTransformerInterf
|
||||
class ANALYSIS_EXPORT QgsGDALGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP
|
||||
{
|
||||
public:
|
||||
|
||||
//! Constructor for QgsGDALGeorefTransform
|
||||
QgsGDALGeorefTransform( bool useTPS, unsigned int polynomialOrder );
|
||||
~QgsGDALGeorefTransform() override;
|
||||
|
||||
QgsGcpTransformerInterface *clone() const override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &layerCoords ) override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis = false ) override;
|
||||
int minimumGcpCount() const override;
|
||||
GDALTransformerFunc GDALTransformer() const override;
|
||||
void *GDALTransformerArgs() const override;
|
||||
@ -230,13 +243,13 @@ class ANALYSIS_EXPORT QgsGDALGeorefTransform : public QgsGcpTransformerInterface
|
||||
private:
|
||||
void destroyGdalArgs();
|
||||
|
||||
QVector<QgsPointXY> mMapCoords;
|
||||
QVector<QgsPointXY> mLayerCoords;
|
||||
QVector<QgsPointXY> mSourceCoords;
|
||||
QVector<QgsPointXY> mDestCoordinates;
|
||||
bool mInvertYAxis = false;
|
||||
|
||||
const int mPolynomialOrder;
|
||||
const bool mIsTPSTransform;
|
||||
|
||||
GDALTransformerFunc mGDALTransformer = nullptr;
|
||||
void *mGDALTransformerArgs = nullptr;
|
||||
|
||||
};
|
||||
@ -253,10 +266,12 @@ class ANALYSIS_EXPORT QgsGDALGeorefTransform : public QgsGcpTransformerInterface
|
||||
class ANALYSIS_EXPORT QgsProjectiveGeorefTransform : public QgsGcpTransformerInterface SIP_SKIP
|
||||
{
|
||||
public:
|
||||
|
||||
//! Constructor for QgsProjectiveGeorefTransform
|
||||
QgsProjectiveGeorefTransform();
|
||||
|
||||
QgsGcpTransformerInterface *clone() const override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &layerCoords ) override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis = false ) override;
|
||||
int minimumGcpCount() const override;
|
||||
GDALTransformerFunc GDALTransformer() const override;
|
||||
void *GDALTransformerArgs() const override;
|
||||
|
||||
@ -23,11 +23,11 @@
|
||||
#include "qgsleastsquares.h"
|
||||
|
||||
|
||||
void QgsLeastSquares::linear( const QVector<QgsPointXY> &mapCoords,
|
||||
const QVector<QgsPointXY> &pixelCoords,
|
||||
void QgsLeastSquares::linear( const QVector<QgsPointXY> &sourceCoordinates,
|
||||
const QVector<QgsPointXY> &destinationCoordinates,
|
||||
QgsPointXY &origin, double &pixelXSize, double &pixelYSize )
|
||||
{
|
||||
int n = mapCoords.size();
|
||||
int n = destinationCoordinates.size();
|
||||
if ( n < 2 )
|
||||
{
|
||||
throw std::domain_error( QObject::tr( "Fit to a linear transform requires at least 2 points." ).toLocal8Bit().constData() );
|
||||
@ -36,14 +36,14 @@ void QgsLeastSquares::linear( const QVector<QgsPointXY> &mapCoords,
|
||||
double sumPx( 0 ), sumPy( 0 ), sumPx2( 0 ), sumPy2( 0 ), sumPxMx( 0 ), sumPyMy( 0 ), sumMx( 0 ), sumMy( 0 );
|
||||
for ( int i = 0; i < n; ++i )
|
||||
{
|
||||
sumPx += pixelCoords.at( i ).x();
|
||||
sumPy += pixelCoords.at( i ).y();
|
||||
sumPx2 += std::pow( pixelCoords.at( i ).x(), 2 );
|
||||
sumPy2 += std::pow( pixelCoords.at( i ).y(), 2 );
|
||||
sumPxMx += pixelCoords.at( i ).x() * mapCoords.at( i ).x();
|
||||
sumPyMy += pixelCoords.at( i ).y() * mapCoords.at( i ).y();
|
||||
sumMx += mapCoords.at( i ).x();
|
||||
sumMy += mapCoords.at( i ).y();
|
||||
sumPx += sourceCoordinates.at( i ).x();
|
||||
sumPy += sourceCoordinates.at( i ).y();
|
||||
sumPx2 += std::pow( sourceCoordinates.at( i ).x(), 2 );
|
||||
sumPy2 += std::pow( sourceCoordinates.at( i ).y(), 2 );
|
||||
sumPxMx += sourceCoordinates.at( i ).x() * destinationCoordinates.at( i ).x();
|
||||
sumPyMy += sourceCoordinates.at( i ).y() * destinationCoordinates.at( i ).y();
|
||||
sumMx += destinationCoordinates.at( i ).x();
|
||||
sumMy += destinationCoordinates.at( i ).y();
|
||||
}
|
||||
|
||||
double deltaX = n * sumPx2 - std::pow( sumPx, 2 );
|
||||
@ -62,30 +62,39 @@ void QgsLeastSquares::linear( const QVector<QgsPointXY> &mapCoords,
|
||||
}
|
||||
|
||||
|
||||
void QgsLeastSquares::helmert( const QVector<QgsPointXY> &mapCoords,
|
||||
const QVector<QgsPointXY> &pixelCoords,
|
||||
void QgsLeastSquares::helmert( const QVector<QgsPointXY> &sourceCoordinates,
|
||||
const QVector<QgsPointXY> &destinationCoordinates,
|
||||
QgsPointXY &origin, double &pixelSize,
|
||||
double &rotation )
|
||||
{
|
||||
int n = mapCoords.size();
|
||||
int n = destinationCoordinates.size();
|
||||
if ( n < 2 )
|
||||
{
|
||||
throw std::domain_error( QObject::tr( "Fit to a Helmert transform requires at least 2 points." ).toLocal8Bit().constData() );
|
||||
}
|
||||
|
||||
double A = 0, B = 0, C = 0, D = 0, E = 0, F = 0, G = 0, H = 0, I = 0, J = 0;
|
||||
double A = 0;
|
||||
double B = 0;
|
||||
double C = 0;
|
||||
double D = 0;
|
||||
double E = 0;
|
||||
double F = 0;
|
||||
double G = 0;
|
||||
double H = 0;
|
||||
double I = 0;
|
||||
double J = 0;
|
||||
for ( int i = 0; i < n; ++i )
|
||||
{
|
||||
A += pixelCoords.at( i ).x();
|
||||
B += pixelCoords.at( i ).y();
|
||||
C += mapCoords.at( i ).x();
|
||||
D += mapCoords.at( i ).y();
|
||||
E += mapCoords.at( i ).x() * pixelCoords.at( i ).x();
|
||||
F += mapCoords.at( i ).y() * pixelCoords.at( i ).y();
|
||||
G += std::pow( pixelCoords.at( i ).x(), 2 );
|
||||
H += std::pow( pixelCoords.at( i ).y(), 2 );
|
||||
I += mapCoords.at( i ).x() * pixelCoords.at( i ).y();
|
||||
J += pixelCoords.at( i ).x() * mapCoords.at( i ).y();
|
||||
A += sourceCoordinates.at( i ).x();
|
||||
B += sourceCoordinates.at( i ).y();
|
||||
C += destinationCoordinates.at( i ).x();
|
||||
D += destinationCoordinates.at( i ).y();
|
||||
E += destinationCoordinates.at( i ).x() * sourceCoordinates.at( i ).x();
|
||||
F += destinationCoordinates.at( i ).y() * sourceCoordinates.at( i ).y();
|
||||
G += std::pow( sourceCoordinates.at( i ).x(), 2 );
|
||||
H += std::pow( sourceCoordinates.at( i ).y(), 2 );
|
||||
I += destinationCoordinates.at( i ).x() * sourceCoordinates.at( i ).y();
|
||||
J += sourceCoordinates.at( i ).x() * destinationCoordinates.at( i ).y();
|
||||
}
|
||||
|
||||
/* The least squares fit for the parameters { a, b, x0, y0 } is the solution
|
||||
@ -232,61 +241,59 @@ void normalizeCoordinates( const QVector<QgsPointXY> &coords, QVector<QgsPointXY
|
||||
|
||||
// Fits a homography to the given corresponding points, and
|
||||
// return it in H (row-major format).
|
||||
void QgsLeastSquares::projective( QVector<QgsPointXY> mapCoords,
|
||||
QVector<QgsPointXY> pixelCoords,
|
||||
void QgsLeastSquares::projective( const QVector<QgsPointXY> &sourceCoordinates,
|
||||
const QVector<QgsPointXY> &destinationCoordinates,
|
||||
double H[9] )
|
||||
{
|
||||
Q_ASSERT( mapCoords.size() == pixelCoords.size() );
|
||||
Q_ASSERT( sourceCoordinates.size() == destinationCoordinates.size() );
|
||||
|
||||
if ( mapCoords.size() < 4 )
|
||||
if ( destinationCoordinates.size() < 4 )
|
||||
{
|
||||
throw std::domain_error( QObject::tr( "Fitting a projective transform requires at least 4 corresponding points." ).toLocal8Bit().constData() );
|
||||
}
|
||||
|
||||
QVector<QgsPointXY> mapCoordsNormalized;
|
||||
QVector<QgsPointXY> pixelCoordsNormalized;
|
||||
QVector<QgsPointXY> sourceCoordinatesNormalized;
|
||||
QVector<QgsPointXY> destinationCoordinatesNormalized;
|
||||
|
||||
double normMap[9], denormMap[9];
|
||||
double normPixel[9], denormPixel[9];
|
||||
normalizeCoordinates( mapCoords, mapCoordsNormalized, normMap, denormMap );
|
||||
normalizeCoordinates( pixelCoords, pixelCoordsNormalized, normPixel, denormPixel );
|
||||
mapCoords = mapCoordsNormalized;
|
||||
pixelCoords = pixelCoordsNormalized;
|
||||
double normSource[9], denormSource[9];
|
||||
double normDest[9], denormDest[9];
|
||||
normalizeCoordinates( sourceCoordinates, sourceCoordinatesNormalized, normSource, denormSource );
|
||||
normalizeCoordinates( destinationCoordinates, destinationCoordinatesNormalized, normDest, denormDest );
|
||||
|
||||
// GSL does not support a full SVD, so we artificially add a linear dependent row
|
||||
// to the matrix in case the system is underconstrained.
|
||||
uint m = std::max( 9u, ( uint )mapCoords.size() * 2u );
|
||||
uint m = std::max( 9u, ( uint )destinationCoordinatesNormalized.size() * 2u );
|
||||
uint n = 9;
|
||||
gsl_matrix *S = gsl_matrix_alloc( m, n );
|
||||
|
||||
for ( int i = 0; i < mapCoords.size(); i++ )
|
||||
for ( int i = 0; i < destinationCoordinatesNormalized.size(); i++ )
|
||||
{
|
||||
gsl_matrix_set( S, i * 2, 0, pixelCoords[i].x() );
|
||||
gsl_matrix_set( S, i * 2, 1, pixelCoords[i].y() );
|
||||
gsl_matrix_set( S, i * 2, 0, sourceCoordinatesNormalized[i].x() );
|
||||
gsl_matrix_set( S, i * 2, 1, sourceCoordinatesNormalized[i].y() );
|
||||
gsl_matrix_set( S, i * 2, 2, 1.0 );
|
||||
|
||||
gsl_matrix_set( S, i * 2, 3, 0.0 );
|
||||
gsl_matrix_set( S, i * 2, 4, 0.0 );
|
||||
gsl_matrix_set( S, i * 2, 5, 0.0 );
|
||||
|
||||
gsl_matrix_set( S, i * 2, 6, -mapCoords[i].x()*pixelCoords[i].x() );
|
||||
gsl_matrix_set( S, i * 2, 7, -mapCoords[i].x()*pixelCoords[i].y() );
|
||||
gsl_matrix_set( S, i * 2, 8, -mapCoords[i].x() * 1.0 );
|
||||
gsl_matrix_set( S, i * 2, 6, -destinationCoordinatesNormalized[i].x()*sourceCoordinatesNormalized[i].x() );
|
||||
gsl_matrix_set( S, i * 2, 7, -destinationCoordinatesNormalized[i].x()*sourceCoordinatesNormalized[i].y() );
|
||||
gsl_matrix_set( S, i * 2, 8, -destinationCoordinatesNormalized[i].x() * 1.0 );
|
||||
|
||||
gsl_matrix_set( S, i * 2 + 1, 0, 0.0 );
|
||||
gsl_matrix_set( S, i * 2 + 1, 1, 0.0 );
|
||||
gsl_matrix_set( S, i * 2 + 1, 2, 0.0 );
|
||||
|
||||
gsl_matrix_set( S, i * 2 + 1, 3, pixelCoords[i].x() );
|
||||
gsl_matrix_set( S, i * 2 + 1, 4, pixelCoords[i].y() );
|
||||
gsl_matrix_set( S, i * 2 + 1, 3, sourceCoordinatesNormalized[i].x() );
|
||||
gsl_matrix_set( S, i * 2 + 1, 4, sourceCoordinatesNormalized[i].y() );
|
||||
gsl_matrix_set( S, i * 2 + 1, 5, 1.0 );
|
||||
|
||||
gsl_matrix_set( S, i * 2 + 1, 6, -mapCoords[i].y()*pixelCoords[i].x() );
|
||||
gsl_matrix_set( S, i * 2 + 1, 7, -mapCoords[i].y()*pixelCoords[i].y() );
|
||||
gsl_matrix_set( S, i * 2 + 1, 8, -mapCoords[i].y() * 1.0 );
|
||||
gsl_matrix_set( S, i * 2 + 1, 6, -destinationCoordinatesNormalized[i].y()*sourceCoordinatesNormalized[i].x() );
|
||||
gsl_matrix_set( S, i * 2 + 1, 7, -destinationCoordinatesNormalized[i].y()*sourceCoordinatesNormalized[i].y() );
|
||||
gsl_matrix_set( S, i * 2 + 1, 8, -destinationCoordinatesNormalized[i].y() * 1.0 );
|
||||
}
|
||||
|
||||
if ( mapCoords.size() == 4 )
|
||||
if ( destinationCoordinatesNormalized.size() == 4 )
|
||||
{
|
||||
// The GSL SVD routine only supports matrices with rows >= columns (m >= n)
|
||||
// Unfortunately, we can't use the SVD of the transpose (i.e. S^T = (U D V^T)^T = V D U^T)
|
||||
@ -321,13 +328,13 @@ void QgsLeastSquares::projective( QVector<QgsPointXY> mapCoords,
|
||||
gsl_matrix *prodMatrix = gsl_matrix_alloc( 3, 3 );
|
||||
|
||||
gsl_matrix_view Hmatrix = gsl_matrix_view_array( H, 3, 3 );
|
||||
gsl_matrix_view normPixelMatrix = gsl_matrix_view_array( normPixel, 3, 3 );
|
||||
gsl_matrix_view denormMapMatrix = gsl_matrix_view_array( denormMap, 3, 3 );
|
||||
gsl_matrix_view normSourceMatrix = gsl_matrix_view_array( normSource, 3, 3 );
|
||||
gsl_matrix_view denormDestMatrix = gsl_matrix_view_array( denormDest, 3, 3 );
|
||||
|
||||
// Change coordinate frame of image and pre-image from normalized to map and pixel coordinates.
|
||||
// Change coordinate frame of image and pre-image from normalized to destination and source coordinates.
|
||||
// H' = denormalizeMapCoords*H*normalizePixelCoords
|
||||
gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &Hmatrix.matrix, &normPixelMatrix.matrix, 0.0, prodMatrix );
|
||||
gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &denormMapMatrix.matrix, prodMatrix, 0.0, &Hmatrix.matrix );
|
||||
gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &Hmatrix.matrix, &normSourceMatrix.matrix, 0.0, prodMatrix );
|
||||
gsl_blas_dgemm( CblasNoTrans, CblasNoTrans, 1.0, &denormDestMatrix.matrix, prodMatrix, 0.0, &Hmatrix.matrix );
|
||||
|
||||
gsl_matrix_free( prodMatrix );
|
||||
gsl_matrix_free( S );
|
||||
|
||||
@ -34,17 +34,17 @@ class ANALYSIS_EXPORT QgsLeastSquares
|
||||
public:
|
||||
|
||||
/**
|
||||
* Transforms the point at \a origin in-place, using a linear transformation calculated from the list of map and pixel Ground Control Points (GCPs).
|
||||
* Transforms the point at \a origin in-place, using a linear transformation calculated from the list of source and destination Ground Control Points (GCPs).
|
||||
*/
|
||||
static void linear( const QVector<QgsPointXY> &mapCoords,
|
||||
const QVector<QgsPointXY> &pixelCoords,
|
||||
static void linear( const QVector<QgsPointXY> &sourceCoordinates,
|
||||
const QVector<QgsPointXY> &destinationCoordinates,
|
||||
QgsPointXY &origin, double &pixelXSize, double &pixelYSize );
|
||||
|
||||
/**
|
||||
* Transforms the point at \a origin in-place, using a helmert transformation calculated from the list of map and pixel Ground Control Points (GCPs).
|
||||
* Transforms the point at \a origin in-place, using a helmert transformation calculated from the list of source and destination Ground Control Points (GCPs).
|
||||
*/
|
||||
static void helmert( const QVector<QgsPointXY> &mapCoords,
|
||||
const QVector<QgsPointXY> &pixelCoords,
|
||||
static void helmert( const QVector<QgsPointXY> &sourceCoordinates,
|
||||
const QVector<QgsPointXY> &destinationCoordinates,
|
||||
QgsPointXY &origin, double &pixelSize, double &rotation );
|
||||
|
||||
#if 0
|
||||
@ -54,10 +54,10 @@ class ANALYSIS_EXPORT QgsLeastSquares
|
||||
#endif
|
||||
|
||||
/**
|
||||
* Calculates projective parameters from the list of map and pixel Ground Control Points (GCPs).
|
||||
* Calculates projective parameters from the list of source and destination Ground Control Points (GCPs).
|
||||
*/
|
||||
static void projective( QVector<QgsPointXY> mapCoords,
|
||||
QVector<QgsPointXY> pixelCoords,
|
||||
static void projective( const QVector<QgsPointXY> &sourceCoordinates,
|
||||
const QVector<QgsPointXY> &destinationCoordinates,
|
||||
double H[9] );
|
||||
|
||||
};
|
||||
|
||||
@ -89,7 +89,7 @@ void QgsGCPListModel::updateModel()
|
||||
|
||||
if ( mGeorefTransform )
|
||||
{
|
||||
bTransformUpdated = mGeorefTransform->updateParametersFromGcps( mapCoords, pixelCoords );
|
||||
bTransformUpdated = mGeorefTransform->updateParametersFromGcps( pixelCoords, mapCoords, true );
|
||||
mapUnitsPossible = mGeorefTransform->providesAccurateInverseTransformation();
|
||||
}
|
||||
|
||||
|
||||
@ -2004,7 +2004,7 @@ bool QgsGeoreferencerMainWindow::updateGeorefTransform()
|
||||
return false;
|
||||
|
||||
// Parametrize the transform with GCPs
|
||||
if ( !mGeorefTransform.updateParametersFromGcps( mapCoords, pixelCoords ) )
|
||||
if ( !mGeorefTransform.updateParametersFromGcps( pixelCoords, mapCoords, true ) )
|
||||
{
|
||||
return false;
|
||||
}
|
||||
|
||||
@ -73,36 +73,36 @@ bool QgsGeorefTransform::parametersInitialized() const
|
||||
QgsGcpTransformerInterface *QgsGeorefTransform::clone() const
|
||||
{
|
||||
std::unique_ptr< QgsGeorefTransform > res( new QgsGeorefTransform( *this ) );
|
||||
res->updateParametersFromGcps( mMapCoords, mLayerCoords );
|
||||
res->updateParametersFromGcps( mSourceCoordinates, mDestinationCoordinates, mInvertYAxis );
|
||||
return res.release();
|
||||
}
|
||||
|
||||
bool QgsGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords )
|
||||
bool QgsGeorefTransform::updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis )
|
||||
{
|
||||
mMapCoords = mapCoords;
|
||||
mLayerCoords = pixelCoords;
|
||||
mSourceCoordinates = sourceCoordinates;
|
||||
mDestinationCoordinates = destinationCoordinates;
|
||||
mInvertYAxis = invertYAxis;
|
||||
|
||||
if ( !mGeorefTransformImplementation )
|
||||
{
|
||||
return false;
|
||||
}
|
||||
if ( mapCoords.size() != pixelCoords.size() ) // Defensive sanity check
|
||||
if ( sourceCoordinates.size() != destinationCoordinates.size() ) // Defensive sanity check
|
||||
{
|
||||
throw ( std::domain_error( "Internal error: GCP mapping is not one-to-one" ) );
|
||||
}
|
||||
if ( mapCoords.size() < minimumGcpCount() )
|
||||
if ( sourceCoordinates.size() < minimumGcpCount() )
|
||||
{
|
||||
return false;
|
||||
}
|
||||
if ( mRasterChangeCoords.hasCrs() )
|
||||
{
|
||||
QVector<QgsPointXY> pixelCoordsCorrect = mRasterChangeCoords.getPixelCoords( pixelCoords );
|
||||
mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( mapCoords, pixelCoordsCorrect );
|
||||
pixelCoordsCorrect.clear();
|
||||
QVector<QgsPointXY> pixelCoordsCorrect = mRasterChangeCoords.getPixelCoords( sourceCoordinates );
|
||||
mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( sourceCoordinates, pixelCoordsCorrect, invertYAxis );
|
||||
}
|
||||
else
|
||||
{
|
||||
mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( mapCoords, pixelCoords );
|
||||
mParametersInitialized = mGeorefTransformImplementation->updateParametersFromGcps( sourceCoordinates, destinationCoordinates, invertYAxis );
|
||||
}
|
||||
return mParametersInitialized;
|
||||
}
|
||||
|
||||
@ -70,7 +70,7 @@ class QgsGeorefTransform : public QgsGcpTransformerInterface
|
||||
bool parametersInitialized() const;
|
||||
|
||||
QgsGcpTransformerInterface *clone() const override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &mapCoords, const QVector<QgsPointXY> &pixelCoords ) override;
|
||||
bool updateParametersFromGcps( const QVector<QgsPointXY> &sourceCoordinates, const QVector<QgsPointXY> &destinationCoordinates, bool invertYAxis = false ) override;
|
||||
int minimumGcpCount() const override;
|
||||
TransformMethod method() const override;
|
||||
GDALTransformerFunc GDALTransformer() const override;
|
||||
@ -112,8 +112,9 @@ class QgsGeorefTransform : public QgsGcpTransformerInterface
|
||||
// convenience wrapper around GDALTransformerFunc
|
||||
bool gdal_transform( const QgsPointXY &src, QgsPointXY &dst, int dstToSrc ) const;
|
||||
|
||||
QVector<QgsPointXY> mMapCoords;
|
||||
QVector<QgsPointXY> mLayerCoords;
|
||||
QVector<QgsPointXY> mSourceCoordinates;
|
||||
QVector<QgsPointXY> mDestinationCoordinates;
|
||||
bool mInvertYAxis = false;
|
||||
|
||||
std::unique_ptr< QgsGcpTransformerInterface > mGeorefTransformImplementation;
|
||||
|
||||
|
||||
@ -46,6 +46,7 @@ endmacro (ADD_QGIS_TEST)
|
||||
#############################################################
|
||||
# Tests:
|
||||
set(TESTS
|
||||
testqgsgcptransformer.cpp
|
||||
testqgsgeometrysnapper.cpp
|
||||
testqgsinterpolator.cpp
|
||||
testqgsprocessing.cpp
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
Loading…
x
Reference in New Issue
Block a user