Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Add some geometry utils for interpolating points on lines
  • Loading branch information
nyalldawson committed Apr 6, 2018
1 parent 78a6118 commit 2e7455c
Show file tree
Hide file tree
Showing 5 changed files with 240 additions and 14 deletions.
50 changes: 50 additions & 0 deletions python/core/geometry/qgsgeometryutils.sip.in
Expand Up @@ -402,6 +402,56 @@ M value is computed if one of this point have M.
# pr is a 3D point: 'PointZM (3 4 1 1)'

.. versionadded:: 3.0
%End

static QgsPointXY interpolatePointOnLine( double x1, double y1, double x2, double y2, double fraction );
%Docstring
Interpolates the position of a point a ``fraction`` of the way along
the line from (``x1``, ``y1``) to (``x2``, ``y2``).

Usually the ``fraction`` should be between 0 and 1, where 0 represents the
point at the start of the line (``x1``, ``y1``) and 1 represents
the end of the line (``x2``, ``y2``). However, it is possible to
use a ``fraction`` < 0 or > 1, in which case the returned point
is extrapolated from the supplied line.

.. versionadded:: 3.0.2

.. seealso:: :py:func:`interpolatePointOnLineByValue`
%End

static QgsPoint interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, double fraction );
%Docstring
Interpolates the position of a point a ``fraction`` of the way along
the line from ``p1`` to ``p2``.

Usually the ``fraction`` should be between 0 and 1, where 0 represents the
point at the start of the line (``p1``) and 1 represents
the end of the line (``p2``). However, it is possible to
use a ``fraction`` < 0 or > 1, in which case the returned point
is extrapolated from the supplied line.

Any Z or M values present in the input points will also be interpolated
and present in the returned point.

.. versionadded:: 3.0.2

.. seealso:: :py:func:`interpolatePointOnLineByValue`
%End

static QgsPointXY interpolatePointOnLineByValue( double x1, double y1, double v1, double x2, double y2, double v2, double value );
%Docstring
Interpolates the position of a point along the line from (``x1``, ``y1``)
to (``x2``, ``y2``).

The position is interpolated using a supplied target ``value`` and the value
at the start of the line (``v1``) and end of the line (``v2``). The returned
point will be linearly interpolated to match position corresponding to
the target ``value``.

.. versionadded:: 3.0.2

.. seealso:: :py:func:`interpolatePointOnLine`
%End

static double gradient( const QgsPoint &pt1, const QgsPoint &pt2 );
Expand Down
18 changes: 4 additions & 14 deletions src/core/geometry/qgsgeometry.cpp
Expand Up @@ -2709,16 +2709,6 @@ QgsGeometry QgsGeometry::smooth( const unsigned int iterations, const double off
}
}

inline QgsPoint interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, const double offset )
{
const double _offset = 1 - offset;
return QgsPoint( p1.wkbType(),
p1.x() * _offset + p2.x() * offset,
p1.y() * _offset + p2.y() * offset,
p1.is3D() ? p1.z() * _offset + p2.z() * offset : std::numeric_limits<double>::quiet_NaN(),
p1.isMeasure() ? p1.m() * _offset + p2.m() * offset : std::numeric_limits<double>::quiet_NaN() );
}

std::unique_ptr< QgsLineString > smoothCurve( const QgsLineString &line, const unsigned int iterations,
const double offset, double squareDistThreshold, double maxAngleRads,
bool isRing )
Expand Down Expand Up @@ -2774,21 +2764,21 @@ std::unique_ptr< QgsLineString > smoothCurve( const QgsLineString &line, const u
if ( !isRing )
{
if ( !skipFirst )
outputLine << ( i == 0 ? result->pointN( i ) : interpolatePointOnLine( p1, p2, offset ) );
outputLine << ( i == 0 ? result->pointN( i ) : QgsGeometryUtils::interpolatePointOnLine( p1, p2, offset ) );
if ( !skipLast )
outputLine << ( i == result->numPoints() - 2 ? result->pointN( i + 1 ) : interpolatePointOnLine( p1, p2, 1.0 - offset ) );
outputLine << ( i == result->numPoints() - 2 ? result->pointN( i + 1 ) : QgsGeometryUtils::interpolatePointOnLine( p1, p2, 1.0 - offset ) );
else
outputLine << p2;
}
else
{
// ring
if ( !skipFirst )
outputLine << interpolatePointOnLine( p1, p2, offset );
outputLine << QgsGeometryUtils::interpolatePointOnLine( p1, p2, offset );
else if ( i == 0 )
outputLine << p1;
if ( !skipLast )
outputLine << interpolatePointOnLine( p1, p2, 1.0 - offset );
outputLine << QgsGeometryUtils::interpolatePointOnLine( p1, p2, 1.0 - offset );
else
outputLine << p2;
}
Expand Down
26 changes: 26 additions & 0 deletions src/core/geometry/qgsgeometryutils.cpp
Expand Up @@ -1175,6 +1175,32 @@ QgsPoint QgsGeometryUtils::midpoint( const QgsPoint &pt1, const QgsPoint &pt2 )
return QgsPoint( pType, x, y, z, m );
}

QgsPoint QgsGeometryUtils::interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, const double fraction )
{
const double _fraction = 1 - fraction;
return QgsPoint( p1.wkbType(),
p1.x() * _fraction + p2.x() * fraction,
p1.y() * _fraction + p2.y() * fraction,
p1.is3D() ? p1.z() * _fraction + p2.z() * fraction : std::numeric_limits<double>::quiet_NaN(),
p1.isMeasure() ? p1.m() * _fraction + p2.m() * fraction : std::numeric_limits<double>::quiet_NaN() );
}

QgsPointXY QgsGeometryUtils::interpolatePointOnLine( const double x1, const double y1, const double x2, const double y2, const double fraction )
{
const double deltaX = ( x2 - x1 ) * fraction;
const double deltaY = ( y2 - y1 ) * fraction;
return QgsPointXY( x1 + deltaX, y1 + deltaY );
}

QgsPointXY QgsGeometryUtils::interpolatePointOnLineByValue( const double x1, const double y1, const double v1, const double x2, const double y2, const double v2, const double value )
{
if ( qgsDoubleNear( v1, v2 ) )
return QgsPointXY( x1, y1 );

const double fraction = ( value - v1 ) / ( v2 - v1 );
return interpolatePointOnLine( x1, y1, x2, y2, fraction );
}

double QgsGeometryUtils::gradient( const QgsPoint &pt1, const QgsPoint &pt2 )
{
double delta_x = pt2.x() - pt1.x();
Expand Down
47 changes: 47 additions & 0 deletions src/core/geometry/qgsgeometryutils.h
Expand Up @@ -434,6 +434,53 @@ class CORE_EXPORT QgsGeometryUtils
*/
static QgsPoint midpoint( const QgsPoint &pt1, const QgsPoint &pt2 );

/**
* Interpolates the position of a point a \a fraction of the way along
* the line from (\a x1, \a y1) to (\a x2, \a y2).
*
* Usually the \a fraction should be between 0 and 1, where 0 represents the
* point at the start of the line (\a x1, \a y1) and 1 represents
* the end of the line (\a x2, \a y2). However, it is possible to
* use a \a fraction < 0 or > 1, in which case the returned point
* is extrapolated from the supplied line.
*
* \since QGIS 3.0.2
* \see interpolatePointOnLineByValue()
*/
static QgsPointXY interpolatePointOnLine( double x1, double y1, double x2, double y2, double fraction );

/**
* Interpolates the position of a point a \a fraction of the way along
* the line from \a p1 to \a p2.
*
* Usually the \a fraction should be between 0 and 1, where 0 represents the
* point at the start of the line (\a p1) and 1 represents
* the end of the line (\a p2). However, it is possible to
* use a \a fraction < 0 or > 1, in which case the returned point
* is extrapolated from the supplied line.
*
* Any Z or M values present in the input points will also be interpolated
* and present in the returned point.
*
* \since QGIS 3.0.2
* \see interpolatePointOnLineByValue()
*/
static QgsPoint interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, double fraction );

/**
* Interpolates the position of a point along the line from (\a x1, \a y1)
* to (\a x2, \a y2).
*
* The position is interpolated using a supplied target \a value and the value
* at the start of the line (\a v1) and end of the line (\a v2). The returned
* point will be linearly interpolated to match position corresponding to
* the target \a value.
*
* \since QGIS 3.0.2
* \see interpolatePointOnLine()
*/
static QgsPointXY interpolatePointOnLineByValue( double x1, double y1, double v1, double x2, double y2, double v2, double value );

/**
* Return the gradient of a line defined by points \a pt1 and \a pt2.
* \param pt1 first point.
Expand Down
113 changes: 113 additions & 0 deletions tests/src/core/testqgsgeometryutils.cpp
Expand Up @@ -63,6 +63,9 @@ class TestQgsGeometryUtils: public QObject
void testTangentPointAndCircle();
void testCircleCircleOuterTangents();
void testGml();
void testInterpolatePointOnLineQgsPoint();
void testInterpolatePointOnLine();
void testInterpolatePointOnLineByValue();
};


Expand Down Expand Up @@ -915,5 +918,115 @@ void TestQgsGeometryUtils::testGml()
QGSCOMPAREGML( elemToString( elm ), expectedGML3_inverted );
}

void TestQgsGeometryUtils::testInterpolatePointOnLineQgsPoint()
{
QgsPoint p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( 10, 0 ), 0 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( 10, 0 ), 1 );
QCOMPARE( p.x(), 10.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( 0, 10 ), 0 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( 0, 10 ), 1 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 10.0 );
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( -10, -6 ), 0.5 );
QCOMPARE( p.x(), -5.0 );
QCOMPARE( p.y(), -3.0 );
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( -10, -6 ), 0.2 );
QCOMPARE( p.x(), -2.0 );
QCOMPARE( p.y(), -1.2 );
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( -10, -6 ), 2 );
QCOMPARE( p.x(), -20.0 );
QCOMPARE( p.y(), -12.0 );
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( -10, -6 ), -1 );
QCOMPARE( p.x(), 10.0 );
QCOMPARE( p.y(), 6.0 );
// with m
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0, 0, 5, QgsWkbTypes::PointM ), QgsPoint( -10, -6, 0, 10, QgsWkbTypes::PointM ), 0.4 );
QCOMPARE( p.wkbType(), QgsWkbTypes::PointM );
QCOMPARE( p.x(), -4.0 );
QCOMPARE( p.y(), -2.4 );
QCOMPARE( p.m(), 7.0 );
// with z
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0, 5, 0, QgsWkbTypes::PointZ ), QgsPoint( -10, -6, 10, 0, QgsWkbTypes::PointZ ), 0.4 );
QCOMPARE( p.wkbType(), QgsWkbTypes::PointZ );
QCOMPARE( p.x(), -4.0 );
QCOMPARE( p.y(), -2.4 );
QCOMPARE( p.z(), 7.0 );
// with zm
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0, 5, 10, QgsWkbTypes::PointZM ), QgsPoint( -10, -6, 10, 5, QgsWkbTypes::PointZM ), 0.4 );
QCOMPARE( p.wkbType(), QgsWkbTypes::PointZM );
QCOMPARE( p.x(), -4.0 );
QCOMPARE( p.y(), -2.4 );
QCOMPARE( p.z(), 7.0 );
QCOMPARE( p.m(), 8.0 );
}

void TestQgsGeometryUtils::testInterpolatePointOnLine()
{
QgsPointXY p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, 10, 0, 0 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, 10, 0, 1 );
QCOMPARE( p.x(), 10.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, 0, 10, 0 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, 0, 10, 1 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 10.0 );
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, -10, -6, 0.5 );
QCOMPARE( p.x(), -5.0 );
QCOMPARE( p.y(), -3.0 );
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, -10, -6, 0.2 );
QCOMPARE( p.x(), -2.0 );
QCOMPARE( p.y(), -1.2 );
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, -10, -6, 2 );
QCOMPARE( p.x(), -20.0 );
QCOMPARE( p.y(), -12.0 );
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, -10, -6, -1 );
QCOMPARE( p.x(), 10.0 );
QCOMPARE( p.y(), 6.0 );
}

void TestQgsGeometryUtils::testInterpolatePointOnLineByValue()
{
QgsPointXY p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 0, 10, 0, 1, 0 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 0, 10, 0, 1, 1 );
QCOMPARE( p.x(), 10.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 5, 0, 10, 15, 5 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 0.0 );
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 15, 0, 10, 5, 5 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 10.0 );
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 3, 2 );
QCOMPARE( p.x(), -5.0 );
QCOMPARE( p.y(), -3.0 );
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 3, 1.4 );
QCOMPARE( p.x(), -2.0 );
QCOMPARE( p.y(), -1.2 );
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 3, -10, -6, 1, 2 );
QCOMPARE( p.x(), -5.0 );
QCOMPARE( p.y(), -3.0 );
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 3, -1 );
QCOMPARE( p.x(), 10.0 );
QCOMPARE( p.y(), 6.0 );
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 3, 5 );
QCOMPARE( p.x(), -20.0 );
QCOMPARE( p.y(), -12.0 );
// v1 == v2, test for no crash
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 1, 1 );
QCOMPARE( p.x(), 0.0 );
QCOMPARE( p.y(), 0.0 );
}

QGSTEST_MAIN( TestQgsGeometryUtils )
#include "testqgsgeometryutils.moc"

0 comments on commit 2e7455c

Please sign in to comment.