Add some circle intersection and tangent utilities to QgsGeometryUtils

Author: nyalldawson

python/core/geometry/qgsgeometryutils.sip.in

@@ -153,6 +153,53 @@ the closest point to the initial ``intersection`` point is returned.
 @param linePoint2 a second point on the line
 @param intersection the initial point and the returned intersection point
 @return true if an intersection has been found
+%End
+
+    static int circleCircleIntersections( QgsPointXY center1, double radius1,
+                                          QgsPointXY center2, double radius2,
+                                          QgsPointXY &intersection1, QgsPointXY &intersection2 );
+%Docstring
+Calculates the intersections points between the circle with center ``center1`` and
+radius ``radius1`` and the circle with center ``center2`` and radius ``radius2``.
+
+If found, the intersection points will be stored in ``intersection1`` and ``intersection2``.
+
+:return: number of intersection points found.
+
+.. versionadded:: 3.2
+%End
+
+    static bool tangentPointAndCircle( const QgsPointXY &center, double radius,
+                                       const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2 );
+%Docstring
+Calculates the tangent points between the circle with the specified ``center`` and ``radius``
+and the point ``p``.
+
+If found, the tangent points will be stored in ``pt1`` and ``pt2``.
+
+.. versionadded:: 3.2
+%End
+
+    static int circleCircleOuterTangents(
+      const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2,
+      QgsPointXY &line1P1, QgsPointXY &line1P2,
+      QgsPointXY &line2P1, QgsPointXY &line2P2 );
+%Docstring
+Calculates the outer tangent points for two circles, centered at ``center1`` and
+``center2`` and with radii of ``radius1`` and ``radius2`` respectively.
+
+The outer tangent points correspond to the points at which the two lines
+which are drawn so that they are tangential to both circles touch
+the circles.
+
+The first tangent line is described by the points
+stored in ``line1P1`` and ``line1P2``,
+and the second line is described by the points stored in ``line2P1``
+and ``line2P2``.
+
+Returns the number of tangents (either 0 or 2).
+
+.. versionadded:: 3.2
 %End
 
     static QgsPoint projectPointOnSegment( const QgsPoint &p, const QgsPoint &s1, const QgsPoint &s2 );

src/core/geometry/qgsgeometryutils.cpp

@@ -363,6 +363,131 @@ bool QgsGeometryUtils::lineCircleIntersection( const QgsPointXY &center, const d
   }
 }
 
+// based on public domain work by 3/26/2005 Tim Voght
+// see http://paulbourke.net/geometry/circlesphere/tvoght.c
+int QgsGeometryUtils::circleCircleIntersections( QgsPointXY center0, const double r0, QgsPointXY center1, const double r1, QgsPointXY &intersection1, QgsPointXY &intersection2 )
+{
+  // determine the straight-line distance between the centers
+  const double d = center0.distance( center1 );
+
+  // check for solvability
+  if ( d > ( r0 + r1 ) )
+  {
+    // no solution. circles do not intersect.
+    return 0;
+  }
+  else if ( d < std::fabs( r0 - r1 ) )
+  {
+    // no solution. one circle is contained in the other
+    return 0;
+  }
+  else if ( qgsDoubleNear( d, 0 ) && ( qgsDoubleNear( r0, r1 ) ) )
+  {
+    // no solutions, the circles coincide
+    return 0;
+  }
+
+  /* 'point 2' is the point where the line through the circle
+   * intersection points crosses the line between the circle
+   * centers.
+  */
+
+  // Determine the distance from point 0 to point 2.
+  const double a = ( ( r0 * r0 ) - ( r1 * r1 ) + ( d * d ) ) / ( 2.0 * d ) ;
+
+  /* dx and dy are the vertical and horizontal distances between
+   * the circle centers.
+   */
+  const double dx = center1.x() - center0.x();
+  const double dy = center1.y() - center0.y();
+
+  // Determine the coordinates of point 2.
+  const double x2 = center0.x() + ( dx * a / d );
+  const double y2 = center0.y() + ( dy * a / d );
+
+  /* Determine the distance from point 2 to either of the
+   * intersection points.
+   */
+  const double h = std::sqrt( ( r0 * r0 ) - ( a * a ) );
+
+  /* Now determine the offsets of the intersection points from
+   * point 2.
+   */
+  const double rx = dy * ( h / d );
+  const double ry = dx * ( h / d );
+
+  // determine the absolute intersection points
+  intersection1 = QgsPointXY( x2 + rx, y2 - ry );
+  intersection2 = QgsPointXY( x2 - rx, y2 +  ry );
+
+  // see if we have 1 or 2 solutions
+  if ( qgsDoubleNear( d, r0 + r1 ) )
+    return 1;
+
+  return 2;
+}
+
+// Using https://stackoverflow.com/a/1351794/1861260
+// and inspired by http://csharphelper.com/blog/2014/11/find-the-tangent-lines-between-a-point-and-a-circle-in-c/
+bool QgsGeometryUtils::tangentPointAndCircle( const QgsPointXY &center, double radius, const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2 )
+{
+  // distance from point to center of circle
+  const double dx = center.x() - p.x();
+  const double dy = center.y() - p.y();
+  const double distanceSquared = dx * dx + dy * dy;
+  const double radiusSquared = radius * radius;
+  if ( distanceSquared < radiusSquared )
+  {
+    // point is inside circle!
+    return false;
+  }
+
+  // distance from point to tangent point, using pythagoras
+  const double distanceToTangent = std::sqrt( distanceSquared - radiusSquared );
+
+  // tangent points are those where the original circle intersects a circle centered
+  // on p with radius distanceToTangent
+  circleCircleIntersections( center, radius, p, distanceToTangent, pt1, pt2 );
+
+  return true;
+}
+
+// inspired by http://csharphelper.com/blog/2014/12/find-the-tangent-lines-between-two-circles-in-c/
+int QgsGeometryUtils::circleCircleOuterTangents( const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2, QgsPointXY &line1P1, QgsPointXY &line1P2, QgsPointXY &line2P1, QgsPointXY &line2P2 )
+{
+  if ( radius1 > radius2 )
+    return circleCircleOuterTangents( center2, radius2, center1, radius1, line1P1, line1P2, line2P1, line2P2 );
+
+  const double radius2a = radius2 - radius1;
+  if ( !tangentPointAndCircle( center2, radius2a, center1, line1P2, line2P2 ) )
+  {
+    // there are no tangents
+    return 0;
+  }
+
+  // get the vector perpendicular to the
+  // first tangent with length radius1
+  QgsVector v1( -( line1P2.y() - center1.y() ), line1P2.x() - center1.x() );
+  const double v1Length = v1.length();
+  v1 = v1 * ( radius1 / v1Length );
+
+  // offset the tangent vector's points
+  line1P1 = center1 + v1;
+  line1P2 = line1P2 + v1;
+
+  // get the vector perpendicular to the
+  // second tangent with length radius1
+  QgsVector v2( line2P2.y() - center1.y(), -( line2P2.x() - center1.x() ) );
+  const double v2Length = v2.length();
+  v2 = v2 * ( radius1 / v2Length );
+
+  // offset the tangent vector's points
+  line2P1 = center1 + v2;
+  line2P2 = line2P2 + v2;
+
+  return 2;
+}
+
 QVector<QgsGeometryUtils::SelfIntersection> QgsGeometryUtils::selfIntersections( const QgsAbstractGeometry *geom, int part, int ring, double tolerance )
 {
   QVector<SelfIntersection> intersections;

src/core/geometry/qgsgeometryutils.h

@@ -154,6 +154,53 @@ class CORE_EXPORT QgsGeometryUtils
                                         const QgsPointXY &linePoint1, const QgsPointXY &linePoint2,
                                         QgsPointXY &intersection SIP_INOUT );
 
+    /**
+     * Calculates the intersections points between the circle with center \a center1 and
+     * radius \a radius1 and the circle with center \a center2 and radius \a radius2.
+     *
+     * If found, the intersection points will be stored in \a intersection1 and \a intersection2.
+     *
+     * \returns number of intersection points found.
+     *
+     * \since QGIS 3.2
+     */
+    static int circleCircleIntersections( QgsPointXY center1, double radius1,
+                                          QgsPointXY center2, double radius2,
+                                          QgsPointXY &intersection1, QgsPointXY &intersection2 );
+
+    /**
+     * Calculates the tangent points between the circle with the specified \a center and \a radius
+     * and the point \a p.
+     *
+     * If found, the tangent points will be stored in \a pt1 and \a pt2.
+     *
+     * \since QGIS 3.2
+     */
+    static bool tangentPointAndCircle( const QgsPointXY &center, double radius,
+                                       const QgsPointXY &p, QgsPointXY &pt1, QgsPointXY &pt2 );
+
+    /**
+     * Calculates the outer tangent points for two circles, centered at \a center1 and
+     * \a center2 and with radii of \a radius1 and \a radius2 respectively.
+     *
+     * The outer tangent points correspond to the points at which the two lines
+     * which are drawn so that they are tangential to both circles touch
+     * the circles.
+     *
+     * The first tangent line is described by the points
+     * stored in \a line1P1 and \a line1P2,
+     * and the second line is described by the points stored in \a line2P1
+     * and \a line2P2.
+     *
+     * Returns the number of tangents (either 0 or 2).
+     *
+     * \since QGIS 3.2
+     */
+    static int circleCircleOuterTangents(
+      const QgsPointXY &center1, double radius1, const QgsPointXY &center2, double radius2,
+      QgsPointXY &line1P1, QgsPointXY &line1P2,
+      QgsPointXY &line2P1, QgsPointXY &line2P2 );
+
     /**
      * \brief Project the point on a segment
      * \param p The point

tests/src/core/testqgsgeometryutils.cpp

@@ -59,6 +59,9 @@ class TestQgsGeometryUtils: public QObject
     void testClosestPoint();
     void testSegmentIntersection();
     void testLineCircleIntersection();
+    void testCircleCircleIntersection();
+    void testTangentPointAndCircle();
+    void testCircleCircleOuterTangents();
     void testGml();
 };
 
@@ -828,6 +831,71 @@ void TestQgsGeometryUtils::testLineCircleIntersection()
   QVERIFY( !isIntersection );
 }
 
+void TestQgsGeometryUtils::testCircleCircleIntersection()
+{
+  QgsPointXY int1;
+  QgsPointXY int2;
+
+  // no intersections
+  QCOMPARE( QgsGeometryUtils::circleCircleIntersections( QgsPointXY( 0, 0 ), 1, QgsPointXY( 2, 0 ), 0.5, int1, int2 ), 0 );
+  QCOMPARE( QgsGeometryUtils::circleCircleIntersections( QgsPointXY( 0, 0 ), 1, QgsPointXY( 0.5, 0.1 ), 0.2, int1, int2 ), 0 );
+  // one intersection
+  QCOMPARE( QgsGeometryUtils::circleCircleIntersections( QgsPointXY( 0, 0 ), 1, QgsPointXY( 3, 0 ), 2, int1, int2 ), 1 );
+  QCOMPARE( int1, QgsPointXY( 1, 0 ) );
+  QCOMPARE( int2, QgsPointXY( 1, 0 ) );
+  // two intersections
+  QCOMPARE( QgsGeometryUtils::circleCircleIntersections( QgsPointXY( 5, 3 ), 2, QgsPointXY( 7, -1 ), 4, int1, int2 ), 2 );
+  QGSCOMPARENEAR( int1.x(), 3.8, 0.001 );
+  QGSCOMPARENEAR( int1.y(), 1.4, 0.001 );
+  QGSCOMPARENEAR( int2.x(), 7.0, 0.001 );
+  QGSCOMPARENEAR( int2.y(), 3.0, 0.001 );
+}
+
+void TestQgsGeometryUtils::testTangentPointAndCircle()
+{
+  QgsPointXY t1;
+  QgsPointXY t2;
+  QVERIFY( !QgsGeometryUtils::tangentPointAndCircle( QgsPointXY( 1, 2 ), 4, QgsPointXY( 1, 2 ), t1, t2 ) );
+  QVERIFY( QgsGeometryUtils::tangentPointAndCircle( QgsPointXY( 1, 2 ), 4, QgsPointXY( 8, 4 ), t1, t2 ) );
+  QGSCOMPARENEAR( t1.x(), 4.03, 0.01 );
+  QGSCOMPARENEAR( t1.y(), -0.61, 0.01 );
+  QGSCOMPARENEAR( t2.x(), 2.2, 0.01 );
+  QGSCOMPARENEAR( t2.y(), 5.82, 0.01 );
+}
+
+void TestQgsGeometryUtils::testCircleCircleOuterTangents()
+{
+  QgsPointXY l1p1;
+  QgsPointXY l1p2;
+  QgsPointXY l2p1;
+  QgsPointXY l2p2;
+
+  // no tangents
+  QCOMPARE( QgsGeometryUtils::circleCircleOuterTangents( QgsPointXY( 1, 2 ), 4, QgsPointXY( 2, 3 ), 1, l1p1, l1p2, l2p1, l2p2 ), 0 );
+
+  // tangents
+  QCOMPARE( QgsGeometryUtils::circleCircleOuterTangents( QgsPointXY( 1, 2 ), 1, QgsPointXY( 10, 3 ), 4, l1p1, l1p2, l2p1, l2p2 ), 2 );
+  QGSCOMPARENEAR( l1p1.x(), 0.566, 0.01 );
+  QGSCOMPARENEAR( l1p1.y(), 2.901, 0.01 );
+  QGSCOMPARENEAR( l1p2.x(), 8.266, 0.01 );
+  QGSCOMPARENEAR( l1p2.y(), 6.604, 0.01 );
+  QGSCOMPARENEAR( l2p1.x(), 0.7749, 0.01 );
+  QGSCOMPARENEAR( l2p1.y(), 1.025, 0.01 );
+  QGSCOMPARENEAR( l2p2.x(), 9.099, 0.01 );
+  QGSCOMPARENEAR( l2p2.y(), -0.897, 0.01 );
+
+  // larger circle first
+  QCOMPARE( QgsGeometryUtils::circleCircleOuterTangents( QgsPointXY( 10, 3 ), 4, QgsPointXY( 1, 2 ), 1, l1p1, l1p2, l2p1, l2p2 ), 2 );
+  QGSCOMPARENEAR( l1p1.x(), 0.566, 0.01 );
+  QGSCOMPARENEAR( l1p1.y(), 2.901, 0.01 );
+  QGSCOMPARENEAR( l1p2.x(), 8.266, 0.01 );
+  QGSCOMPARENEAR( l1p2.y(), 6.604, 0.01 );
+  QGSCOMPARENEAR( l2p1.x(), 0.7749, 0.01 );
+  QGSCOMPARENEAR( l2p1.y(), 1.025, 0.01 );
+  QGSCOMPARENEAR( l2p2.x(), 9.099, 0.01 );
+  QGSCOMPARENEAR( l2p2.y(), -0.897, 0.01 );
+}
+
 void TestQgsGeometryUtils::testGml()
 {
   QgsPoint point = QgsPoint( 1, 2 );