Skip to content

Commit 2e7455c

Browse files
nyalldawsonnyalldawson
committed
Add some geometry utils for interpolating points on lines
1 parent 78a6118 commit 2e7455c

File tree

5 files changed

+240
-14
lines changed

5 files changed

+240
-14
lines changed

‎python/core/geometry/qgsgeometryutils.sip.in

+50
Original file line numberDiff line numberDiff line change
@@ -402,6 +402,56 @@ M value is computed if one of this point have M.
402402
# pr is a 3D point: 'PointZM (3 4 1 1)'
403403

404404
.. versionadded:: 3.0
405+
%End
406+
407+
static QgsPointXY interpolatePointOnLine( double x1, double y1, double x2, double y2, double fraction );
408+
%Docstring
409+
Interpolates the position of a point a ``fraction`` of the way along
410+
the line from (``x1``, ``y1``) to (``x2``, ``y2``).
411+
412+
Usually the ``fraction`` should be between 0 and 1, where 0 represents the
413+
point at the start of the line (``x1``, ``y1``) and 1 represents
414+
the end of the line (``x2``, ``y2``). However, it is possible to
415+
use a ``fraction`` < 0 or > 1, in which case the returned point
416+
is extrapolated from the supplied line.
417+
418+
.. versionadded:: 3.0.2
419+
420+
.. seealso:: :py:func:`interpolatePointOnLineByValue`
421+
%End
422+
423+
static QgsPoint interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, double fraction );
424+
%Docstring
425+
Interpolates the position of a point a ``fraction`` of the way along
426+
the line from ``p1`` to ``p2``.
427+
428+
Usually the ``fraction`` should be between 0 and 1, where 0 represents the
429+
point at the start of the line (``p1``) and 1 represents
430+
the end of the line (``p2``). However, it is possible to
431+
use a ``fraction`` < 0 or > 1, in which case the returned point
432+
is extrapolated from the supplied line.
433+
434+
Any Z or M values present in the input points will also be interpolated
435+
and present in the returned point.
436+
437+
.. versionadded:: 3.0.2
438+
439+
.. seealso:: :py:func:`interpolatePointOnLineByValue`
440+
%End
441+
442+
static QgsPointXY interpolatePointOnLineByValue( double x1, double y1, double v1, double x2, double y2, double v2, double value );
443+
%Docstring
444+
Interpolates the position of a point along the line from (``x1``, ``y1``)
445+
to (``x2``, ``y2``).
446+
447+
The position is interpolated using a supplied target ``value`` and the value
448+
at the start of the line (``v1``) and end of the line (``v2``). The returned
449+
point will be linearly interpolated to match position corresponding to
450+
the target ``value``.
451+
452+
.. versionadded:: 3.0.2
453+
454+
.. seealso:: :py:func:`interpolatePointOnLine`
405455
%End
406456

407457
static double gradient( const QgsPoint &pt1, const QgsPoint &pt2 );

‎src/core/geometry/qgsgeometry.cpp

+4-14
Original file line numberDiff line numberDiff line change
@@ -2709,16 +2709,6 @@ QgsGeometry QgsGeometry::smooth( const unsigned int iterations, const double off
27092709
}
27102710
}
27112711

2712-
inline QgsPoint interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, const double offset )
2713-
{
2714-
const double _offset = 1 - offset;
2715-
return QgsPoint( p1.wkbType(),
2716-
p1.x() * _offset + p2.x() * offset,
2717-
p1.y() * _offset + p2.y() * offset,
2718-
p1.is3D() ? p1.z() * _offset + p2.z() * offset : std::numeric_limits<double>::quiet_NaN(),
2719-
p1.isMeasure() ? p1.m() * _offset + p2.m() * offset : std::numeric_limits<double>::quiet_NaN() );
2720-
}
2721-
27222712
std::unique_ptr< QgsLineString > smoothCurve( const QgsLineString &line, const unsigned int iterations,
27232713
const double offset, double squareDistThreshold, double maxAngleRads,
27242714
bool isRing )
@@ -2774,21 +2764,21 @@ std::unique_ptr< QgsLineString > smoothCurve( const QgsLineString &line, const u
27742764
if ( !isRing )
27752765
{
27762766
if ( !skipFirst )
2777-
outputLine << ( i == 0 ? result->pointN( i ) : interpolatePointOnLine( p1, p2, offset ) );
2767+
outputLine << ( i == 0 ? result->pointN( i ) : QgsGeometryUtils::interpolatePointOnLine( p1, p2, offset ) );
27782768
if ( !skipLast )
2779-
outputLine << ( i == result->numPoints() - 2 ? result->pointN( i + 1 ) : interpolatePointOnLine( p1, p2, 1.0 - offset ) );
2769+
outputLine << ( i == result->numPoints() - 2 ? result->pointN( i + 1 ) : QgsGeometryUtils::interpolatePointOnLine( p1, p2, 1.0 - offset ) );
27802770
else
27812771
outputLine << p2;
27822772
}
27832773
else
27842774
{
27852775
// ring
27862776
if ( !skipFirst )
2787-
outputLine << interpolatePointOnLine( p1, p2, offset );
2777+
outputLine << QgsGeometryUtils::interpolatePointOnLine( p1, p2, offset );
27882778
else if ( i == 0 )
27892779
outputLine << p1;
27902780
if ( !skipLast )
2791-
outputLine << interpolatePointOnLine( p1, p2, 1.0 - offset );
2781+
outputLine << QgsGeometryUtils::interpolatePointOnLine( p1, p2, 1.0 - offset );
27922782
else
27932783
outputLine << p2;
27942784
}

‎src/core/geometry/qgsgeometryutils.cpp

+26
Original file line numberDiff line numberDiff line change
@@ -1175,6 +1175,32 @@ QgsPoint QgsGeometryUtils::midpoint( const QgsPoint &pt1, const QgsPoint &pt2 )
11751175
return QgsPoint( pType, x, y, z, m );
11761176
}
11771177

1178+
QgsPoint QgsGeometryUtils::interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, const double fraction )
1179+
{
1180+
const double _fraction = 1 - fraction;
1181+
return QgsPoint( p1.wkbType(),
1182+
p1.x() * _fraction + p2.x() * fraction,
1183+
p1.y() * _fraction + p2.y() * fraction,
1184+
p1.is3D() ? p1.z() * _fraction + p2.z() * fraction : std::numeric_limits<double>::quiet_NaN(),
1185+
p1.isMeasure() ? p1.m() * _fraction + p2.m() * fraction : std::numeric_limits<double>::quiet_NaN() );
1186+
}
1187+
1188+
QgsPointXY QgsGeometryUtils::interpolatePointOnLine( const double x1, const double y1, const double x2, const double y2, const double fraction )
1189+
{
1190+
const double deltaX = ( x2 - x1 ) * fraction;
1191+
const double deltaY = ( y2 - y1 ) * fraction;
1192+
return QgsPointXY( x1 + deltaX, y1 + deltaY );
1193+
}
1194+
1195+
QgsPointXY QgsGeometryUtils::interpolatePointOnLineByValue( const double x1, const double y1, const double v1, const double x2, const double y2, const double v2, const double value )
1196+
{
1197+
if ( qgsDoubleNear( v1, v2 ) )
1198+
return QgsPointXY( x1, y1 );
1199+
1200+
const double fraction = ( value - v1 ) / ( v2 - v1 );
1201+
return interpolatePointOnLine( x1, y1, x2, y2, fraction );
1202+
}
1203+
11781204
double QgsGeometryUtils::gradient( const QgsPoint &pt1, const QgsPoint &pt2 )
11791205
{
11801206
double delta_x = pt2.x() - pt1.x();

‎src/core/geometry/qgsgeometryutils.h

+47
Original file line numberDiff line numberDiff line change
@@ -434,6 +434,53 @@ class CORE_EXPORT QgsGeometryUtils
434434
*/
435435
static QgsPoint midpoint( const QgsPoint &pt1, const QgsPoint &pt2 );
436436

437+
/**
438+
* Interpolates the position of a point a \a fraction of the way along
439+
* the line from (\a x1, \a y1) to (\a x2, \a y2).
440+
*
441+
* Usually the \a fraction should be between 0 and 1, where 0 represents the
442+
* point at the start of the line (\a x1, \a y1) and 1 represents
443+
* the end of the line (\a x2, \a y2). However, it is possible to
444+
* use a \a fraction < 0 or > 1, in which case the returned point
445+
* is extrapolated from the supplied line.
446+
*
447+
* \since QGIS 3.0.2
448+
* \see interpolatePointOnLineByValue()
449+
*/
450+
static QgsPointXY interpolatePointOnLine( double x1, double y1, double x2, double y2, double fraction );
451+
452+
/**
453+
* Interpolates the position of a point a \a fraction of the way along
454+
* the line from \a p1 to \a p2.
455+
*
456+
* Usually the \a fraction should be between 0 and 1, where 0 represents the
457+
* point at the start of the line (\a p1) and 1 represents
458+
* the end of the line (\a p2). However, it is possible to
459+
* use a \a fraction < 0 or > 1, in which case the returned point
460+
* is extrapolated from the supplied line.
461+
*
462+
* Any Z or M values present in the input points will also be interpolated
463+
* and present in the returned point.
464+
*
465+
* \since QGIS 3.0.2
466+
* \see interpolatePointOnLineByValue()
467+
*/
468+
static QgsPoint interpolatePointOnLine( const QgsPoint &p1, const QgsPoint &p2, double fraction );
469+
470+
/**
471+
* Interpolates the position of a point along the line from (\a x1, \a y1)
472+
* to (\a x2, \a y2).
473+
*
474+
* The position is interpolated using a supplied target \a value and the value
475+
* at the start of the line (\a v1) and end of the line (\a v2). The returned
476+
* point will be linearly interpolated to match position corresponding to
477+
* the target \a value.
478+
*
479+
* \since QGIS 3.0.2
480+
* \see interpolatePointOnLine()
481+
*/
482+
static QgsPointXY interpolatePointOnLineByValue( double x1, double y1, double v1, double x2, double y2, double v2, double value );
483+
437484
/**
438485
* Return the gradient of a line defined by points \a pt1 and \a pt2.
439486
* \param pt1 first point.

‎tests/src/core/testqgsgeometryutils.cpp

+113
Original file line numberDiff line numberDiff line change
@@ -63,6 +63,9 @@ class TestQgsGeometryUtils: public QObject
6363
void testTangentPointAndCircle();
6464
void testCircleCircleOuterTangents();
6565
void testGml();
66+
void testInterpolatePointOnLineQgsPoint();
67+
void testInterpolatePointOnLine();
68+
void testInterpolatePointOnLineByValue();
6669
};
6770

6871

@@ -915,5 +918,115 @@ void TestQgsGeometryUtils::testGml()
915918
QGSCOMPAREGML( elemToString( elm ), expectedGML3_inverted );
916919
}
917920

921+
void TestQgsGeometryUtils::testInterpolatePointOnLineQgsPoint()
922+
{
923+
QgsPoint p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( 10, 0 ), 0 );
924+
QCOMPARE( p.x(), 0.0 );
925+
QCOMPARE( p.y(), 0.0 );
926+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( 10, 0 ), 1 );
927+
QCOMPARE( p.x(), 10.0 );
928+
QCOMPARE( p.y(), 0.0 );
929+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( 0, 10 ), 0 );
930+
QCOMPARE( p.x(), 0.0 );
931+
QCOMPARE( p.y(), 0.0 );
932+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( 0, 10 ), 1 );
933+
QCOMPARE( p.x(), 0.0 );
934+
QCOMPARE( p.y(), 10.0 );
935+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( -10, -6 ), 0.5 );
936+
QCOMPARE( p.x(), -5.0 );
937+
QCOMPARE( p.y(), -3.0 );
938+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( -10, -6 ), 0.2 );
939+
QCOMPARE( p.x(), -2.0 );
940+
QCOMPARE( p.y(), -1.2 );
941+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( -10, -6 ), 2 );
942+
QCOMPARE( p.x(), -20.0 );
943+
QCOMPARE( p.y(), -12.0 );
944+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0 ), QgsPoint( -10, -6 ), -1 );
945+
QCOMPARE( p.x(), 10.0 );
946+
QCOMPARE( p.y(), 6.0 );
947+
// with m
948+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0, 0, 5, QgsWkbTypes::PointM ), QgsPoint( -10, -6, 0, 10, QgsWkbTypes::PointM ), 0.4 );
949+
QCOMPARE( p.wkbType(), QgsWkbTypes::PointM );
950+
QCOMPARE( p.x(), -4.0 );
951+
QCOMPARE( p.y(), -2.4 );
952+
QCOMPARE( p.m(), 7.0 );
953+
// with z
954+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0, 5, 0, QgsWkbTypes::PointZ ), QgsPoint( -10, -6, 10, 0, QgsWkbTypes::PointZ ), 0.4 );
955+
QCOMPARE( p.wkbType(), QgsWkbTypes::PointZ );
956+
QCOMPARE( p.x(), -4.0 );
957+
QCOMPARE( p.y(), -2.4 );
958+
QCOMPARE( p.z(), 7.0 );
959+
// with zm
960+
p = QgsGeometryUtils::interpolatePointOnLine( QgsPoint( 0, 0, 5, 10, QgsWkbTypes::PointZM ), QgsPoint( -10, -6, 10, 5, QgsWkbTypes::PointZM ), 0.4 );
961+
QCOMPARE( p.wkbType(), QgsWkbTypes::PointZM );
962+
QCOMPARE( p.x(), -4.0 );
963+
QCOMPARE( p.y(), -2.4 );
964+
QCOMPARE( p.z(), 7.0 );
965+
QCOMPARE( p.m(), 8.0 );
966+
}
967+
968+
void TestQgsGeometryUtils::testInterpolatePointOnLine()
969+
{
970+
QgsPointXY p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, 10, 0, 0 );
971+
QCOMPARE( p.x(), 0.0 );
972+
QCOMPARE( p.y(), 0.0 );
973+
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, 10, 0, 1 );
974+
QCOMPARE( p.x(), 10.0 );
975+
QCOMPARE( p.y(), 0.0 );
976+
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, 0, 10, 0 );
977+
QCOMPARE( p.x(), 0.0 );
978+
QCOMPARE( p.y(), 0.0 );
979+
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, 0, 10, 1 );
980+
QCOMPARE( p.x(), 0.0 );
981+
QCOMPARE( p.y(), 10.0 );
982+
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, -10, -6, 0.5 );
983+
QCOMPARE( p.x(), -5.0 );
984+
QCOMPARE( p.y(), -3.0 );
985+
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, -10, -6, 0.2 );
986+
QCOMPARE( p.x(), -2.0 );
987+
QCOMPARE( p.y(), -1.2 );
988+
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, -10, -6, 2 );
989+
QCOMPARE( p.x(), -20.0 );
990+
QCOMPARE( p.y(), -12.0 );
991+
p = QgsGeometryUtils::interpolatePointOnLine( 0, 0, -10, -6, -1 );
992+
QCOMPARE( p.x(), 10.0 );
993+
QCOMPARE( p.y(), 6.0 );
994+
}
995+
996+
void TestQgsGeometryUtils::testInterpolatePointOnLineByValue()
997+
{
998+
QgsPointXY p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 0, 10, 0, 1, 0 );
999+
QCOMPARE( p.x(), 0.0 );
1000+
QCOMPARE( p.y(), 0.0 );
1001+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 0, 10, 0, 1, 1 );
1002+
QCOMPARE( p.x(), 10.0 );
1003+
QCOMPARE( p.y(), 0.0 );
1004+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 5, 0, 10, 15, 5 );
1005+
QCOMPARE( p.x(), 0.0 );
1006+
QCOMPARE( p.y(), 0.0 );
1007+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 15, 0, 10, 5, 5 );
1008+
QCOMPARE( p.x(), 0.0 );
1009+
QCOMPARE( p.y(), 10.0 );
1010+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 3, 2 );
1011+
QCOMPARE( p.x(), -5.0 );
1012+
QCOMPARE( p.y(), -3.0 );
1013+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 3, 1.4 );
1014+
QCOMPARE( p.x(), -2.0 );
1015+
QCOMPARE( p.y(), -1.2 );
1016+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 3, -10, -6, 1, 2 );
1017+
QCOMPARE( p.x(), -5.0 );
1018+
QCOMPARE( p.y(), -3.0 );
1019+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 3, -1 );
1020+
QCOMPARE( p.x(), 10.0 );
1021+
QCOMPARE( p.y(), 6.0 );
1022+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 3, 5 );
1023+
QCOMPARE( p.x(), -20.0 );
1024+
QCOMPARE( p.y(), -12.0 );
1025+
// v1 == v2, test for no crash
1026+
p = QgsGeometryUtils::interpolatePointOnLineByValue( 0, 0, 1, -10, -6, 1, 1 );
1027+
QCOMPARE( p.x(), 0.0 );
1028+
QCOMPARE( p.y(), 0.0 );
1029+
}
1030+
9181031
QGSTEST_MAIN( TestQgsGeometryUtils )
9191032
#include "testqgsgeometryutils.moc"

0 commit comments

Comments
 (0)