Skip to content

Commit

Permalink
Add QgsGeometryUtils method to interpolate a point on an arc given a …
Browse files Browse the repository at this point in the history 
…distance
  • Loading branch information
@nyalldawson
nyalldawson committed Aug 15, 2018
1 parent 56fd4e3 commit 5b0bdbd
Show file tree
Hide file tree
Showing 4 changed files with 107 additions and 2 deletions.
15 changes: 14 additions & 1 deletion python/core/auto_generated/geometry/qgsgeometryutils.sip.in
View file
Expand Up @@ -232,6 +232,16 @@ Returns a point a specified ``distance`` toward a second point.
%End


static QgsPoint interpolatePointOnArc( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double distance );
%Docstring
Interpolates a point on an arc defined by three points, ``pt1``, ``pt2`` and ``pt3``. The arc will be
interpolated by the specified ``distance`` from ``pt1``.

Any z or m values present in the points will also be linearly interpolated in the output.

.. versionadded:: 3.4
%End

static double ccwAngle( double dy, double dx );
%Docstring
Returns the counter clockwise angle between a line with components dx, dy and the line with dx > 0 and dy = 0
Expand All @@ -245,7 +255,10 @@ Returns radius and center of the circle through pt1, pt2, pt3

static bool circleClockwise( double angle1, double angle2, double angle3 );
%Docstring
Returns true if circle is ordered clockwise
Returns true if the circle defined by three angles is ordered clockwise.

The angles are defined counter-clockwise from the origin, i.e. using
Euclidean angles as opposed to geographic "North up" angles.
%End

static bool circleAngleBetween( double angle, double angle1, double angle2, bool clockwise );
Expand Down
25 changes: 25 additions & 0 deletions src/core/geometry/qgsgeometryutils.cpp
View file
Expand Up @@ -574,6 +574,31 @@ void QgsGeometryUtils::pointOnLineWithDistance( double x1, double y1, double x2,
}
}

QgsPoint QgsGeometryUtils::interpolatePointOnArc( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double distance )
{
double centerX, centerY, radius;
circleCenterRadius( pt1, pt2, pt3, radius, centerX, centerY );

const double theta = distance / radius; // angle subtended
const double anglePt1 = std::atan2( pt1.y() - centerY, pt1.x() - centerX );
const double anglePt2 = std::atan2( pt2.y() - centerY, pt2.x() - centerX );
const double anglePt3 = std::atan2( pt3.y() - centerY, pt3.x() - centerX );
const bool isClockwise = circleClockwise( anglePt1, anglePt2, anglePt3 );
const double angleDest = anglePt1 + ( isClockwise ? -theta : theta );

const double x = centerX + radius * ( std::cos( angleDest ) );
const double y = centerY + radius * ( std::sin( angleDest ) );

const double z = pt1.is3D() ?
interpolateArcValue( angleDest, anglePt1, anglePt2, anglePt3, pt1.z(), pt2.z(), pt3.z() )
: 0;
const double m = pt1.isMeasure() ?
interpolateArcValue( angleDest, anglePt1, anglePt2, anglePt3, pt1.m(), pt2.m(), pt3.m() )
: 0;

return QgsPoint( pt1.wkbType(), x, y, z, m );
}

double QgsGeometryUtils::ccwAngle( double dy, double dx )
{
double angle = std::atan2( dy, dx ) * 180 / M_PI;
Expand Down
17 changes: 16 additions & 1 deletion src/core/geometry/qgsgeometryutils.h
View file
Expand Up @@ -264,14 +264,29 @@ class CORE_EXPORT QgsGeometryUtils
double *z1 = nullptr, double *z2 = nullptr, double *z = nullptr,
double *m1 = nullptr, double *m2 = nullptr, double *m = nullptr ) SIP_SKIP;

/**
* Interpolates a point on an arc defined by three points, \a pt1, \a pt2 and \a pt3. The arc will be
* interpolated by the specified \a distance from \a pt1.
*
* Any z or m values present in the points will also be linearly interpolated in the output.
*
* \since QGIS 3.4
*/
static QgsPoint interpolatePointOnArc( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double distance );

//! Returns the counter clockwise angle between a line with components dx, dy and the line with dx > 0 and dy = 0
static double ccwAngle( double dy, double dx );

//! Returns radius and center of the circle through pt1, pt2, pt3
static void circleCenterRadius( const QgsPoint &pt1, const QgsPoint &pt2, const QgsPoint &pt3, double &radius SIP_OUT,
double &centerX SIP_OUT, double &centerY SIP_OUT );

//! Returns true if circle is ordered clockwise
/**
* Returns true if the circle defined by three angles is ordered clockwise.
*
* The angles are defined counter-clockwise from the origin, i.e. using
* Euclidean angles as opposed to geographic "North up" angles.
*/
static bool circleClockwise( double angle1, double angle2, double angle3 );

//! Returns true if, in a circle, angle is between angle1 and angle2
Expand Down
52 changes: 52 additions & 0 deletions tests/src/core/testqgsgeometryutils.cpp
View file
Expand Up @@ -69,6 +69,7 @@ class TestQgsGeometryUtils: public QObject
void testInterpolatePointOnLine();
void testInterpolatePointOnLineByValue();
void testPointOnLineWithDistance();
void interpolatePointOnArc();
};


Expand Down Expand Up @@ -143,6 +144,8 @@ void TestQgsGeometryUtils::testCircleClockwise_data()
QTest::newRow( "circleClockwise1" ) << 10.0 << 300.0 << 270.0 << true;
QTest::newRow( "circleClockwise2" ) << 270.0 << 300.0 << 10.0 << false;
QTest::newRow( "circleClockwise3" ) << 260.0 << 245.0 << 243.0 << true;
QTest::newRow( "circleClockwise4" ) << -90.0 << 0.0 << 90.0 << false;
QTest::newRow( "circleClockwise5" ) << 0.0 << 90.0 << 180.0 << false;
}

void TestQgsGeometryUtils::testCircleClockwise()
Expand Down Expand Up @@ -1179,5 +1182,54 @@ void TestQgsGeometryUtils::testPointOnLineWithDistance()

}

void TestQgsGeometryUtils::interpolatePointOnArc()
{
QgsPoint p;
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0, 1, 2 ), QgsPoint( 11, 1, 3, 4 ), QgsPoint( 12, 0, 13, 14 ), 0 );
QGSCOMPARENEAR( p.x(), 10.0, 0.00001 );
QGSCOMPARENEAR( p.y(), 0.0, 0.00001 );
QGSCOMPARENEAR( p.z(), 1.0, 0.00001 );
QGSCOMPARENEAR( p.m(), 2.0, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0, 1, 2 ), QgsPoint( 11, 1, 3, 4 ), QgsPoint( 12, 0, 13, 14 ), 1 );
QGSCOMPARENEAR( p.x(), 10.459698, 0.00001 );
QGSCOMPARENEAR( p.y(), 0.841471, 0.00001 );
QGSCOMPARENEAR( p.z(), 2.273240, 0.00001 );
QGSCOMPARENEAR( p.m(), 3.273240, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0, 1, 2 ), QgsPoint( 11, 1, 3, 4 ), QgsPoint( 12, 0, 13, 14 ), 2 );
QGSCOMPARENEAR( p.x(), 11.416147, 0.00001 );
QGSCOMPARENEAR( p.y(), 0.909297, 0.00001 );
QGSCOMPARENEAR( p.z(), 5.732395, 0.00001 );
QGSCOMPARENEAR( p.m(), 6.732395, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0 ), QgsPoint( 11, 1 ), QgsPoint( 12, 0 ), 3.141592 );
QGSCOMPARENEAR( p.x(), 12.0, 0.00001 );
QGSCOMPARENEAR( p.y(), 0.0, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0 ), QgsPoint( 11, 1 ), QgsPoint( 12, 0 ), 3.2 );
QGSCOMPARENEAR( p.x(), 11.998295, 0.00001 );
QGSCOMPARENEAR( p.y(), -0.058374, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0 ), QgsPoint( 11, 1 ), QgsPoint( 12, 0 ), 5 );
QGSCOMPARENEAR( p.x(), 10.716338, 0.00001 );
QGSCOMPARENEAR( p.y(), -0.958924, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0, 1, 2 ), QgsPoint( 11, 1, 3, 4 ), QgsPoint( 12, 0, 13, 14 ), 3.141592 * 2 );
QGSCOMPARENEAR( p.x(), 10, 0.00001 );
QGSCOMPARENEAR( p.y(), 0, 0.00001 );
QGSCOMPARENEAR( p.z(), 32.99999, 0.00001 );
QGSCOMPARENEAR( p.m(), 33.999992, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0 ), QgsPoint( 8, 2 ), QgsPoint( 6, 0 ), 0 );
QGSCOMPARENEAR( p.x(), 10.0, 0.00001 );
QGSCOMPARENEAR( p.y(), 0.0, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0 ), QgsPoint( 8, 2 ), QgsPoint( 6, 0 ), 1 );
QGSCOMPARENEAR( p.x(), 9.755165, 0.00001 );
QGSCOMPARENEAR( p.y(), 0.958851, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0 ), QgsPoint( 8, 2 ), QgsPoint( 6, 0 ), 3.141592 );
QGSCOMPARENEAR( p.x(), 8.0, 0.00001 );
QGSCOMPARENEAR( p.y(), 2.0, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0 ), QgsPoint( 8, 2 ), QgsPoint( 6, 0 ), 3.141592 * 2 );
QGSCOMPARENEAR( p.x(), 6.0, 0.00001 );
QGSCOMPARENEAR( p.y(), 0.0, 0.00001 );
p = QgsGeometryUtils::interpolatePointOnArc( QgsPoint( 10, 0 ), QgsPoint( 8, 2 ), QgsPoint( 6, 0 ), 3.141592 * 3 );
QGSCOMPARENEAR( p.x(), 8.0, 0.00001 );
QGSCOMPARENEAR( p.y(), -2.0, 0.00001 );
}

QGSTEST_MAIN( TestQgsGeometryUtils )
#include "testqgsgeometryutils.moc"

0 comments on commit 5b0bdbd

Please sign in to comment.