Skip to content

Commit

Permalink
Rename pixel/layer/map coordinates to "source" and "destination", to
Browse files Browse the repository at this point in the history 
make classes more generic and avoid confusion with non-raster based
GCPs

Also flip function arguments to source, destination order instead
of destination, source
  • Loading branch information
@nyalldawson
nyalldawson committed Feb 21, 2021
1 parent 30fac03 commit ce54321
Show file tree
Hide file tree
Showing 11 changed files with 1,652 additions and 156 deletions.
14 changes: 10 additions & 4 deletions python/analysis/auto_generated/georeferencing/qgsgcptransformer.sip.in
View file
Expand Up @@ -37,6 +37,9 @@ based on a transformation method and a list of GCPs.
};

QgsGcpTransformerInterface();
%Docstring
Constructor for QgsGcpTransformerInterface
%End

virtual ~QgsGcpTransformerInterface();

Expand All @@ -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
Expand Down Expand Up @@ -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.

Expand Down
129 changes: 80 additions & 49 deletions src/analysis/georeferencing/qgsgcptransformer.cpp
View file
Expand Up @@ -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();
Expand All @@ -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;
}

Expand Down Expand Up @@ -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;
}
}
Expand All @@ -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;
}
}
Expand All @@ -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;
}

Expand All @@ -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
Expand Down Expand Up @@ -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;
}
}
Expand All @@ -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;
}
}
Expand All @@ -309,30 +332,31 @@ 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++ )
{
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();
Expand Down Expand Up @@ -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;
Expand Down

0 comments on commit ce54321

Please sign in to comment.