[geos-commits] [SCM] GEOS branch main updated. 789196b9e19eed3299f47ea5437ac089ac841e7d
git at osgeo.org
git at osgeo.org
Wed Apr 29 07:00:56 PDT 2026
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, main has been updated
via 789196b9e19eed3299f47ea5437ac089ac841e7d (commit)
via 8d42241e9c557ef9ca6bc97ecc71dd636ea82973 (commit)
via e5beda552e1c949d234a4c1e45457ee0fd032e04 (commit)
via 2ef81880b9d0c9032902960f350e58591661f7bd (commit)
from 7445d96c1d10df18418b73f85a1b6c4d05381498 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 789196b9e19eed3299f47ea5437ac089ac841e7d
Author: Daniel Baston <dbaston at gmail.com>
Date: Sun Apr 19 21:20:54 2026 -0400
DistanceOp: Support curved geometries
diff --git a/include/geos/algorithm/CircularArcIntersector.h b/include/geos/algorithm/CircularArcIntersector.h
index 8affd3c65..b3a7a41ca 100644
--- a/include/geos/algorithm/CircularArcIntersector.h
+++ b/include/geos/algorithm/CircularArcIntersector.h
@@ -130,19 +130,6 @@ private:
void addArcSegmentIntersectionPoint(const CoordinateXY& computedIntPt, const CircularArc& lhs,
const geom::CoordinateSequence& seq, std::size_t pos0, std::size_t pos1, bool useSegEndpoints);
- /** Determines whether and where two circles intersect a line segment.
- *
- * @param center The center point of the circle
- * @param r The radius of the circle
- * @param p0 The first point of the line segment
- * @param p1 The second point of the line segment
- * @param isect0 Set to the first intersection point, if it exists
- * @param isect1 Set to the second intersection point, if it exists
- * @return The number of intersection points
- */
- static int
- circleIntersects(const CoordinateXY& center, double r, const CoordinateXY& p0, const CoordinateXY& p1, CoordinateXY& isect0, CoordinateXY& isect1);
-
void computeCocircularIntersection(const CircularArc& arc1, const CircularArc& arc2);
/** Checks whether the provided point has already been recorded as an intersection point. */
diff --git a/include/geos/algorithm/CircularArcs.h b/include/geos/algorithm/CircularArcs.h
index e1d72ef91..d8793d561 100644
--- a/include/geos/algorithm/CircularArcs.h
+++ b/include/geos/algorithm/CircularArcs.h
@@ -14,12 +14,13 @@
#pragma once
-#include <array>
-
#include <geos/export.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/Envelope.h>
+#include <array>
+#include <optional>
+
namespace geos {
namespace algorithm {
@@ -45,7 +46,194 @@ public:
/// Return the point defined by a circle center, radius, and angle
static geom::CoordinateXY createPoint(const geom::CoordinateXY& center, double radius, double theta);
-};
+ /** Determines whether and where a circle intersects a line segment.
+ *
+ * @param center The center point of the circle
+ * @param r The radius of the circle
+ * @param p0 The first point of the line segment
+ * @param p1 The second point of the line segment
+ * @param isect0 Set to the first intersection point, if it exists
+ * @param isect1 Set to the second intersection point, if it exists
+ * @return The number of intersection points
+ */
+ static int
+ circleIntersectsSegment(const geom::CoordinateXY& center, double r,
+ const geom::CoordinateXY& p0, const geom::CoordinateXY& p1,
+ geom::CoordinateXY& isect0, geom::CoordinateXY& isect1);
+ /** Returns whether an angle is included between two other angles.
+ *
+ * @param t0 starting angle (radians)
+ * @param t2 ending engle (radians)
+ * @param isCCW whether the interval from t0 to t2 is counterclockwise
+ * @param theta angle to test (radians)
+ * @return true if theta is included between t0 and t2
+ */
+ static bool
+ containsAngle(double t0, double t2, bool isCCW, double theta);
+
+ /** Returns whether a point on a circle is included between two arc endpoints.
+ * It is assumed that the point lies on the circle; this is not tested.
+ *
+ * @param center The center coordinate of the circle
+ * @param p0 The first arc endpoint
+ * @param p2 The second arc endpoint
+ * @param isCCW Whether the arc between p0 and p2 is counterclockwise
+ * @param q The point to check for inclusion
+ * @return true if q is on the arc, false otherwise
+ */
+ static bool
+ containsPointOnCircle(const geom::CoordinateXY& center,
+ const geom::CoordinateXY& p0, const geom::CoordinateXY& p2, bool isCCW,
+ const geom::CoordinateXY& q);
+
+ /** Calculates the distance between an arc and a point.
+ *
+ * @param c The center coordinate of the circle
+ * @param r The radius of the circle
+ * @param p0 The first arc endpoint
+ * @param p2 The second arc endpoint
+ * @param isCCW Whether the arc between p0 and p2 is counterclockwise
+ * @param q The point to calculate distance to
+ * @return the distance between the arc and the point
+ */
+ static double
+ distanceArcPoint(const geom::CoordinateXY& c, double r,
+ const geom::CoordinateXY& p0, const geom::CoordinateXY& p2, bool isCCW,
+ const geom::CoordinateXY& q);
+
+ /** Calculates the distance between an arc and a line segment.
+ * @param c The center coordinate of the circle
+ * @param r The radius of the circle
+ * @param p0 The first arc endpoint
+ * @param p2 The second arc endpoint
+ * @param isCCW Whether the arc between p0 and p2 is counterclockwise
+ * @param q0 The first point of the line segment
+ * @param q1 The second point of the line segment
+ * @return the distance between the arc and the line segment
+ */
+ static double
+ distanceArcSegment(const geom::CoordinateXY& c, double r,
+ const geom::CoordinateXY& p0, const geom::CoordinateXY& p2, bool isCCW,
+ const geom::CoordinateXY& q0, const geom::CoordinateXY& q1);
+
+ /** Calculates the distance between two arcs.
+ *
+ * @param ca the center coordinate of the first arc's circle
+ * @param ra the radius of the first arc's circle
+ * @param a0 the first endpoint of the first arc
+ * @param a2 the second endpoint of the first arc
+ * @param aCCW whether the first arc is counterclockwise
+ * @param cb the center coordinate of the second arc's circle
+ * @param rb the radius of the second arc's circle
+ * @param b0 the first endpoint of the second arc
+ * @param b2 the second endpoint of the second arc
+ * @param bCCW whether the second arc is counterclockwise
+ * @return the distance between the two arcs
+ */
+ static double
+ distanceArcArc(const geom::CoordinateXY& ca, double ra,
+ const geom::CoordinateXY& a0, const geom::CoordinateXY& a2, bool aCCW,
+ const geom::CoordinateXY& cb, double rb,
+ const geom::CoordinateXY& b0, const geom::CoordinateXY& b2, bool bCCW);
+
+ /** Return the closest point on an arc to a specified point.
+ *
+ * @param c the center of the arc
+ * @param r the radius of the arc
+ * @param p0 the first endpoint of the arc
+ * @param p2 the second endpoint of the arc
+ * @param isCCW whether the arc is counterclockwise
+ * @param q the point to find the closest point to
+ * @return the closest point on the arc to the specified point
+ */
+ static geom::CoordinateXY
+ closestPointArcPoint(const geom::CoordinateXY& c, double r,
+ const geom::CoordinateXY& p0, const geom::CoordinateXY& p2,
+ bool isCCW, const geom::CoordinateXY &q);
+
+ /** Return the pair of points having the shortest distance between an arc and a line segment
+ *
+ * @param c the center of the arc
+ * @param r the radius of the arc
+ * @param p0 the first endpoint of the arc
+ * @param p2 the second endpoint of the arc
+ * @param isCCW whether the arc is counterclockwise
+ * @param q0 the first endpoint of the line segment
+ * @param q1 the second endpoint of the line segment
+ * @return array whose first element is a point on the arc and second element is a point on the line segment
+ */
+ static std::array<geom::CoordinateXY, 2>
+ closestPointsArcSegment(const geom::CoordinateXY& c, double r,
+ const geom::CoordinateXY& p0, const geom::CoordinateXY& p2, bool isCCW, const geom::CoordinateXY& q0, const geom::CoordinateXY& q1);
+
+
+ /** Return the pair of points having the shortest distance between two arcs
+ *
+ * @param ca the center of the first arc
+ * @param ra the radius of the first arc
+ * @param a0 the first endpoint of the first arc
+ * @param a2 the second endpoint of the first arc
+ * @param aCCW whether the first arc is counterclockwise
+ * @param cb the center of the second arc
+ * @param rb the radius of the second arc
+ * @param b0 the first endpoint of the second arc
+ * @param b2 the second endpoint of the second arc
+ * @param bCCW whether the second arc is counterclockwise
+ * @return array whose first element is a point on the first arc and second element is a point on the second arc
+ */
+ static std::array<geom::CoordinateXY, 2> closestPointsArcArc(const geom::CoordinateXY& ca,
+ double ra,
+ const geom::CoordinateXY& a0,
+ const geom::CoordinateXY& a2,
+ bool aCCW,
+ const geom::CoordinateXY& cb,
+ double rb,
+ const geom::CoordinateXY& b0,
+ const geom::CoordinateXY& b2,
+ bool bCCW);
+
+ /** Returns a single XY point of intersection between two cocircular arcs, if one exists.
+ *
+ * To identify all intersection point(s) and/or arc(s), with handling of Z and M ordinates, use CircularArcIntersector.
+ *
+ * @param center the center coordinate of the arcs
+ * @param radius the radius of the arcs
+ * @param a0 the first endpoint of the first arc
+ * @param a2 the second endpoint of the first arc
+ * @param aCCW whether the first arc is counterclockwise
+ * @param b0 the first endpoint of the second arc
+ * @param b2 the second endpoint of the second arc
+ * @param bCCW whether the second arc is counterclockwise
+ * @return a point of intersection between the two arcs, if one exists
+ */
+ static std::optional<geom::CoordinateXY>
+ cocircularArcsIntersectionPoint(const geom::CoordinateXY& center, double radius,
+ const geom::CoordinateXY& a0, const geom::CoordinateXY& a2, bool aCCW,
+ const geom::CoordinateXY& b0, const geom::CoordinateXY& b2, bool bCCW);
+
+ /** Returns a single point of intersection between two arcs, if one exists.
+ *
+ * To identify all intersection point(s) and/or arc(s), with handling of Z and M ordinates, use CircularArcIntersector.
+ *
+ * @param ca the center of the first arc
+ * @param ra the radius of the first arc
+ * @param a0 the first endpoint of the first arc
+ * @param a2 the second endpoint of the first arc
+ * @param aCCW whether the first arc is counterclockwise
+ * @param cb the center of the second arc
+ * @param rb the radius of the second arc
+ * @param b0 the first endpoint of the second arc
+ * @param b2 the second endpoint of the second arc
+ * @param bCCW whether the second arc is counterclockwise
+ * @return a point of intersection between the two arcs, if one exists
+ */
+ static std::optional<geom::CoordinateXY>
+ arcIntersectionPoint(const geom::CoordinateXY& ca, double ra,
+ const geom::CoordinateXY& a0, const geom::CoordinateXY& a2, bool aCCW,
+ const geom::CoordinateXY& cb, double rb,
+ const geom::CoordinateXY& b0, const geom::CoordinateXY& b2, bool bCCW);
+
+};
}
}
diff --git a/include/geos/algorithm/PointLocator.h b/include/geos/algorithm/PointLocator.h
index 6221b5db8..be504c2b0 100644
--- a/include/geos/algorithm/PointLocator.h
+++ b/include/geos/algorithm/PointLocator.h
@@ -26,11 +26,13 @@
namespace geos {
namespace geom {
class CoordinateXY;
+class Curve;
class Geometry;
class LinearRing;
class LineString;
class Polygon;
class Point;
+class Surface;
}
}
@@ -96,9 +98,9 @@ private:
geom::Location locate(const geom::CoordinateXY& p, const geom::LineString* l);
- geom::Location locateInPolygonRing(const geom::CoordinateXY& p, const geom::LinearRing* ring);
+ static geom::Location locateInPolygonRing(const geom::CoordinateXY& p, const geom::Curve* ring);
- geom::Location locate(const geom::CoordinateXY& p, const geom::Polygon* poly);
+ static geom::Location locate(const geom::CoordinateXY& p, const geom::Surface* poly);
};
diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h
index bb3d3d2c2..9f9bbc693 100644
--- a/include/geos/geom/CircularArc.h
+++ b/include/geos/geom/CircularArc.h
@@ -190,6 +190,14 @@ public:
return containsAngle(theta);
}
+ CoordinateXY closestPoint(const CoordinateXY& p) const;
+ std::array<CoordinateXY, 2> closestPoints(const CoordinateXY& p1, const CoordinateXY& p2) const;
+ std::array<CoordinateXY, 2> closestPoints(const CircularArc& other) const;
+
+ double distance(const CoordinateXY& p) const;
+ double distance(const CoordinateXY& p1, const CoordinateXY& p2) const;
+ double distance(const CircularArc& other) const;
+
/// Return true if the arc is pointing positive in the y direction
/// at the location of a specified point. The point is assumed to
/// be on the arc.
diff --git a/include/geos/operation/distance/DistanceOp.h b/include/geos/operation/distance/DistanceOp.h
index 03d3d2c92..f45079d72 100644
--- a/include/geos/operation/distance/DistanceOp.h
+++ b/include/geos/operation/distance/DistanceOp.h
@@ -37,11 +37,14 @@
// Forward declarations
namespace geos {
namespace geom {
+class CircularString;
class Coordinate;
class Polygon;
class LineString;
class Point;
class Geometry;
+class SimpleCurve;
+class Surface;
}
}
@@ -191,7 +194,7 @@ private:
void computeContainmentDistance();
void computeInside(std::vector<GeometryLocation> & locs,
- const std::vector<const geom::Polygon*>& polys,
+ const std::vector<const geom::Surface*>& polys,
std::array<GeometryLocation, 2> & locPtPoly);
@@ -202,8 +205,8 @@ private:
void computeFacetDistance();
void computeMinDistanceLines(
- const std::vector<const geom::LineString*>& lines0,
- const std::vector<const geom::LineString*>& lines1,
+ const std::vector<const geom::SimpleCurve*>& lines0,
+ const std::vector<const geom::SimpleCurve*>& lines1,
std::array<GeometryLocation, 2> & locGeom);
void computeMinDistancePoints(
@@ -212,17 +215,33 @@ private:
std::array<GeometryLocation, 2> & locGeom);
void computeMinDistanceLinesPoints(
- const std::vector<const geom::LineString*>& lines0,
+ const std::vector<const geom::SimpleCurve*>& lines0,
const std::vector<const geom::Point*>& points1,
std::array<GeometryLocation, 2> & locGeom);
+ void computeMinDistance(const geom::CircularString* cs0,
+ const geom::LineString* line1,
+ std::array<GeometryLocation, 2> & locGeom);
+
+ void computeMinDistance(const geom::LineString* line0,
+ const geom::CircularString* cs1,
+ std::array<GeometryLocation, 2> & locGeom);
+
void computeMinDistance(const geom::LineString* line0,
const geom::LineString* line1,
std::array<GeometryLocation, 2> & locGeom);
+ void computeMinDistance(const geom::CircularString* cs0,
+ const geom::CircularString* cs1,
+ std::array<GeometryLocation, 2> & locGeom);
+
void computeMinDistance(const geom::LineString* line,
const geom::Point* pt,
std::array<GeometryLocation, 2> & locGeom);
+
+ void computeMinDistance(const geom::CircularString* cs,
+ const geom::Point* pt,
+ std::array<GeometryLocation, 2> & locGeom);
};
diff --git a/src/algorithm/CircularArcIntersector.cpp b/src/algorithm/CircularArcIntersector.cpp
index 1b3bedcaf..31a6eea7c 100644
--- a/src/algorithm/CircularArcIntersector.cpp
+++ b/src/algorithm/CircularArcIntersector.cpp
@@ -190,72 +190,7 @@ interpolateZM(const CircularArc& arc0,
}
}
-int
-CircularArcIntersector::circleIntersects(const CoordinateXY& center,
- double r,
- const CoordinateXY& p0,
- const CoordinateXY& p1,
- CoordinateXY& ret0,
- CoordinateXY& ret1)
-{
- const double& x0 = center.x;
- const double& y0 = center.y;
- Envelope segEnv(p0, p1);
-
- CoordinateXY isect0, isect1;
- int n = 0;
-
- if (p1.x == p0.x) {
- // vertical line
- double x = p1.x;
-
- double A = 1;
- double B = 2*y0;
- double C = x*x - 2*x*x0 + x0*x0 + y0*y0 - r*r;
-
- double d = std::sqrt(B*B - 4*A*C);
- double Y1 = (-B + d)/(2*A);
- double Y2 = (-B - d)/(2*A);
-
- isect0 = {x, Y1};
- isect1 = {x, Y2};
- }
- else {
- double m = (p1.y - p0.y) / (p1.x - p0.x);
- double b = p1.y - p1.x*m;
-
- // Ax^2 + Bx + C = 0
- double A = 1 + m*m;
- double B = -2*x0 + 2*m*b - 2*m*y0;
- double C = x0*x0 + b*b - 2*b*y0 + y0*y0 - r*r;
-
- double d = std::sqrt(B*B - 4*A*C);
- double X1 = (-B + d)/(2*A);
- double X2 = (-B - d)/(2*A);
-
- // TODO use robust quadratic solver such as https://github.com/archermarx/quadratic ?
- // auto [X1, X2] = quadratic::solve(A, B, C);
-
- isect0 = {X1, m* X1 + b};
- isect1 = {X2, m* X2 + b};
- }
-
- if (segEnv.intersects(isect0)) {
- ret0 = isect0;
- if (segEnv.intersects(isect1) && !isect1.equals2D(isect0)) {
- ret1 = isect1;
- n = 2;
- } else {
- n = 1;
- }
- } else if (segEnv.intersects(isect1)) {
- ret0 = isect1;
- n = 1;
- }
-
- return n;
-}
bool
CircularArcIntersector::hasIntersection(const geom::CoordinateXY &p) const {
@@ -287,7 +222,7 @@ CircularArcIntersector::intersects(const CircularArc& arc, const CoordinateSeque
const double r = arc.getRadius();
CoordinateXYZM isect0, isect1;
- auto n = circleIntersects(c, r, seq.getAt<CoordinateXY>(segPos0), seq.getAt<CoordinateXY>(segPos1), isect0, isect1);
+ auto n = CircularArcs::circleIntersectsSegment(c, r, seq.getAt<CoordinateXY>(segPos0), seq.getAt<CoordinateXY>(segPos1), isect0, isect1);
if (n > 0 && arc.containsPointOnCircle(isect0)) {
addArcSegmentIntersectionPoint(isect0, arc, seq, segPos0, segPos1, useSegEndpoints);
diff --git a/src/algorithm/CircularArcs.cpp b/src/algorithm/CircularArcs.cpp
index aee3b7b5d..ee28c96c4 100644
--- a/src/algorithm/CircularArcs.cpp
+++ b/src/algorithm/CircularArcs.cpp
@@ -3,7 +3,7 @@
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
- * Copyright (C) 2024 ISciences, LLC
+ * Copyright (C) 2024-2026 ISciences, LLC
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU Lesser General Public Licence as published
@@ -16,10 +16,13 @@
#include <geos/algorithm/CircularArcs.h>
#include <geos/algorithm/Orientation.h>
#include <geos/geom/Envelope.h>
+#include <geos/geom/LineSegment.h>
#include <geos/geom/Quadrant.h>
#include <geos/math/DD.h>
using geos::geom::CoordinateXY;
+using geos::geom::Envelope;
+using geos::geom::LineSegment;
namespace geos {
namespace algorithm {
@@ -198,5 +201,423 @@ CircularArcs::expandEnvelope(geom::Envelope& e, const geom::CoordinateXY& p0, co
}
}
+int
+CircularArcs::circleIntersectsSegment(const CoordinateXY& center, double r,
+ const CoordinateXY& p0, const CoordinateXY& p1,
+ CoordinateXY& ret0, CoordinateXY& ret1)
+{
+ const double& x0 = center.x;
+ const double& y0 = center.y;
+
+ Envelope segEnv(p0, p1);
+
+ CoordinateXY isect0, isect1;
+ int n = 0;
+
+ if (p1.x == p0.x) {
+ // vertical line
+ double x = p1.x;
+
+ double A = 1;
+ double B = -2*y0;
+ double C = x*x - 2*x*x0 + x0*x0 + y0*y0 - r*r;
+
+ double d = std::sqrt(B*B - 4*A*C);
+ double Y1 = (-B + d)/(2*A);
+ double Y2 = (-B - d)/(2*A);
+
+ isect0 = {x, Y1};
+ isect1 = {x, Y2};
+ }
+ else {
+ double m = (p1.y - p0.y) / (p1.x - p0.x);
+ double b = p1.y - p1.x*m;
+
+ // Ax^2 + Bx + C = 0
+ double A = 1 + m*m;
+ double B = -2*x0 + 2*m*b - 2*m*y0;
+ double C = x0*x0 + b*b - 2*b*y0 + y0*y0 - r*r;
+
+ double d = std::sqrt(B*B - 4*A*C);
+ double X1 = (-B + d)/(2*A);
+ double X2 = (-B - d)/(2*A);
+
+ // TODO use robust quadratic solver such as https://github.com/archermarx/quadratic ?
+ // auto [X1, X2] = quadratic::solve(A, B, C);
+
+ isect0 = {X1, m* X1 + b};
+ isect1 = {X2, m* X2 + b};
+ }
+
+ if (segEnv.intersects(isect0)) {
+ ret0 = isect0;
+ if (segEnv.intersects(isect1) && !isect1.equals2D(isect0)) {
+ ret1 = isect1;
+ n = 2;
+ } else {
+ n = 1;
+ }
+ } else if (segEnv.intersects(isect1)) {
+ ret0 = isect1;
+ n = 1;
+ }
+
+ return n;
+}
+
+bool
+CircularArcs::containsAngle(double t0, double t2, bool isCCW, double theta)
+{
+ if (theta == t0 || theta == t2) {
+ return true;
+ }
+
+ if (isCCW) {
+ std::swap(t0, t2);
+ }
+
+ t2 -= t0;
+ theta -= t0;
+
+ if (t2 < 0) {
+ t2 += 2*MATH_PI;
+ }
+ if (theta < 0) {
+ theta += 2*MATH_PI;
+ }
+
+ return theta >= t2;
+
+}
+
+bool
+CircularArcs::containsPointOnCircle(const CoordinateXY& center,
+ const CoordinateXY& p0, const CoordinateXY& p2, bool isCCW,
+ const CoordinateXY& q)
+{
+ if (p0.equals2D(p2)) {
+ // closed circle contains all points
+ return true;
+ }
+
+ const double t0 = getAngle(p0, center);
+ const double t2 = getAngle(p2, center);
+ const double theta = getAngle(q, center);
+
+ return containsAngle(t0, t2, isCCW, theta);
+}
+
+double
+CircularArcs::distanceArcPoint(const CoordinateXY& c, double r,
+ const CoordinateXY& p0, const CoordinateXY& p2, bool isCCW,
+ const CoordinateXY &q)
+{
+ const CoordinateXY p = closestPointArcPoint(c, r, p0, p2, isCCW, q);
+ return p.distance(q);
+}
+
+CoordinateXY
+CircularArcs::closestPointArcPoint(const CoordinateXY& c, double r,
+ const CoordinateXY& p0, const CoordinateXY& p2, bool isCCW,
+ const CoordinateXY &q)
+{
+ if (!std::isfinite(r)) {
+ return LineSegment::closestPoint(p0, p2, q);
+ }
+
+ const double distanceToCenter = q.distance(c);
+ if (distanceToCenter == 0) {
+ return p0; // arbitrary
+ }
+
+ CoordinateXY i;
+ i.x = c.x + (q.x - c.x) * (r / distanceToCenter);
+ i.y = c.y + (q.y - c.y) * (r / distanceToCenter);
+
+ if (i.equals2D(p0) || i.equals2D(p2) || containsPointOnCircle(c, p0, p2, isCCW, i)) {
+ return i;
+ }
+
+ if (q.distance(p0) < q.distance(p2)) {
+ return p0;
+ }
+
+ return p2;
+}
+
+double
+CircularArcs::distanceArcSegment(const CoordinateXY& c, double r,
+ const CoordinateXY& p0, const CoordinateXY& p2, bool isCCW,
+ const CoordinateXY& q0, const CoordinateXY& q1)
+{
+ const auto pts = closestPointsArcSegment(c, r, p0, p2, isCCW, q0, q1);
+ return pts[0].distance(pts[1]);
+}
+
+std::array<CoordinateXY, 2>
+CircularArcs::closestPointsArcSegment(const CoordinateXY& c, double r,
+ const CoordinateXY& p0, const CoordinateXY& p2, bool isCCW,
+ const CoordinateXY& q0, const CoordinateXY& q1)
+{
+ if (!std::isfinite(r)) {
+ return geom::LineSegment::closestPoints(p0, p2, q0, q1);
+ }
+
+ // Calculate distance between line and circle center
+ const CoordinateXY d = geom::LineSegment::closestPoint(q0, q1, c);
+ const double dist_cd = d.distance(c);
+
+ if (dist_cd <= r) {
+ // Segment is at least partially within circle; check to see if it intersects
+ CoordinateXY isect0, isect1;
+ auto n = circleIntersectsSegment(c, r, q0, q1, isect0, isect1);
+ if (n > 0 && containsPointOnCircle(c, p0, p2, isCCW, isect0)) {
+ return {isect0, isect0};
+ }
+ if (n > 1 && containsPointOnCircle(c, p0, p2, isCCW, isect1)) {
+ return {isect1, isect1};
+ }
+
+ // Closest point is one of the segment endpoints.
+ const auto e0 = closestPointArcPoint(c, r, p0, p2, isCCW, q0);
+ const auto e1 = closestPointArcPoint(c, r, p0, p2, isCCW, q1);
+ if (e0.distance(q0) < e1.distance(q1)) {
+ return {e0, q0};
+ }
+
+ return {e1, q1};
+ }
+
+ // Find closest point on circle
+ const CoordinateXY g{c.x + (d.x - c.x) * (r / dist_cd),
+ c.y + (d.y - c.y) * (r / dist_cd)};
+
+ const bool ptInArc = containsPointOnCircle(c, p0, p2, isCCW, g);
+
+ if (ptInArc) {
+ return {g, d};
+ }
+
+ // Closest distance is between the segment and one of the arc endpoints
+ const auto e0 = LineSegment::closestPoint(q0, q1, p0);
+ const auto e1 = LineSegment::closestPoint(q0, q1, p2);
+
+ if (e0.distance(p0) < e1.distance(p2)) {
+ return {e0, p0};
+ }
+
+ return {e1, p2};
+}
+
+double
+CircularArcs::distanceArcArc(const CoordinateXY& ca, double ra,
+ const CoordinateXY& a0, const CoordinateXY& a2, bool aCCW,
+ const CoordinateXY& cb, double rb,
+ const CoordinateXY& b0, const CoordinateXY& b2, bool bCCW)
+{
+ const auto pts = closestPointsArcArc(ca, ra, a0, a2, aCCW, cb, rb, b0, b2, bCCW);
+ return pts[0].distance(pts[1]);
+}
+
+std::array<CoordinateXY, 2>
+CircularArcs::closestPointsArcArc(const CoordinateXY& ca, double ra,
+ const CoordinateXY& a0, const CoordinateXY& a2, bool aCCW,
+ const CoordinateXY& cb, double rb,
+ const CoordinateXY& b0, const CoordinateXY& b2, bool bCCW)
+{
+ if (!std::isfinite(ra)) {
+ if (!std::isfinite(rb)) {
+ return LineSegment::closestPoints(a0, a2, b0, b2);
+ } else {
+ auto pts = closestPointsArcSegment(cb, rb, b0, b2, bCCW, a0, a2);
+ std::swap(pts[0], pts[1]);
+ return pts;
+ }
+ } else if (!std::isfinite(rb)) {
+ return closestPointsArcSegment(ca, ra, a0, a2, aCCW, b0, b2);
+ }
+
+ const double d = ca.distance(cb);
+
+ const bool disjoint = d > (ra + rb);
+ const bool contained = d < std::abs(ra - rb);
+ const bool concentric = (d == 0);
+ const bool cocircular = concentric && (ra == rb);
+ const bool overlapping = !(disjoint || contained || concentric);
+
+ std::vector<CoordinateXY> ptsA{a0, a2};
+ std::vector<CoordinateXY> ptsB{b0, b2};
+
+ if (overlapping || cocircular) {
+ auto isect = arcIntersectionPoint(ca, ra, a0, a2, aCCW, cb, rb, b0, b2, bCCW);
+ if (isect.has_value()) {
+ return { isect.value(), isect.value() };
+ }
+ }
+
+ if (!concentric) {
+ // Add points on A and B that intersect line from center-a to center-b
+ {
+ double theta = getAngle(cb, ca);
+ double at0 = getAngle(a0, ca);
+ double at2 = getAngle(a2, ca);
+
+ if (containsAngle(at0, at2, aCCW, theta)) {
+ ptsA.push_back(createPoint(ca, ra, theta));
+ }
+ if (containsAngle(at0, at2, aCCW, theta + MATH_PI)) {
+ ptsA.push_back(createPoint(ca, ra, theta + MATH_PI));
+ }
+ }
+
+ {
+ double theta = getAngle(ca, cb);
+ double bt0 = getAngle(b0, cb);
+ double bt2 = getAngle(b2, cb);
+
+ if (containsAngle(bt0, bt2, bCCW, theta)) {
+ ptsB.push_back(createPoint(cb, rb, theta));
+ }
+ if (containsAngle(bt0, bt2, bCCW, theta + MATH_PI)) {
+ ptsB.push_back(createPoint(cb, rb, theta + MATH_PI));
+ }
+
+ }
+ }
+
+ // Add points on A that intersect line from CA to B0
+ if (containsPointOnCircle(ca, a0, a2, aCCW, b0)) {
+ ptsA.push_back(CircularArcs::createPoint(ca, ra, getAngle(b0, ca)));
+ }
+
+ // Add points on A that intersect line from CA to B2
+ if (containsPointOnCircle(ca, a0, a2, aCCW, b2)) {
+ ptsA.push_back(CircularArcs::createPoint(ca, ra, getAngle(b2, ca)));
+ }
+
+ // Add points on B that intersect line from CB to A0
+ if (containsPointOnCircle(cb, b0, b2, bCCW, a0)) {
+ ptsB.push_back(CircularArcs::createPoint(cb, rb, getAngle(a0, cb)));
+ }
+
+ // Add points on B that intersect line from CB to A2
+ if (containsPointOnCircle(cb, b0, b2, bCCW, a2)) {
+ ptsB.push_back(CircularArcs::createPoint(cb, rb, getAngle(a2, cb)));
+ }
+
+ double minDist = std::numeric_limits<double>::max();
+
+ std::array<CoordinateXY, 2> ret;
+
+ for (const auto& pa : ptsA) {
+ for (const auto& pb : ptsB) {
+ const double dist = pa.distance(pb);
+ if (dist < minDist) {
+ minDist = dist;
+ ret[0] = pa;
+ ret[1] = pb;
+ }
+ }
+ }
+
+ return ret;
+}
+
+std::optional<CoordinateXY>
+CircularArcs::cocircularArcsIntersectionPoint(const CoordinateXY& center, double radius,
+ const CoordinateXY& a0, const CoordinateXY& a2, bool aCCW,
+ const CoordinateXY& b0, const CoordinateXY& b2, bool bCCW)
+{
+ // The body of this function is adapted and simplified from CircularArcIntersector.
+
+ // Make both inputs counter-clockwise for the purpose of determining intersections
+ const double ap0 = Angle::normalizePositive(getAngle(aCCW ? a0 : a2, center));
+ const double ap1 = Angle::normalizePositive(getAngle(aCCW ? a2 : a0, center));
+ const double bp0 = Angle::normalizePositive(getAngle(bCCW ? b0 : b2, center));
+ const double bp1 = Angle::normalizePositive(getAngle(bCCW ? b2 : b0, center));
+
+ // Check if B starts within A
+ if (Angle::isWithinCCW(bp0, ap0, ap1)) {
+ return createPoint(center, radius, bp0);
+ }
+
+ // Check if B ends within A
+ if (Angle::isWithinCCW(bp1, ap0, ap1)) {
+ return createPoint(center, radius, bp1);
+ }
+
+ // Check if A starts within B
+ if (Angle::isWithinCCW(ap0, bp0 , bp1))
+ {
+ return createPoint(center, radius, ap0);
+ }
+
+ return std::nullopt;
+}
+
+std::optional<CoordinateXY>
+CircularArcs::arcIntersectionPoint(const CoordinateXY& ca, double ra,
+ const CoordinateXY& a0, const CoordinateXY& a2, bool aCCW,
+ const CoordinateXY& cb, double rb,
+ const CoordinateXY& b0, const CoordinateXY& b2, bool bCCW)
+{
+ // The body of this function is adapted and simplified from CircularArcIntersector.
+
+ const auto d = ca.distance(cb);
+
+ if (d > ra + rb) {
+ // Circles are disjoint
+ return std::nullopt;
+ }
+
+ if (d < std::abs(ra - rb)) {
+ // One circle contained within the other; arcs cannot intersect
+ return std::nullopt;
+ }
+
+ // a: the distance from c1 to the "radical line", which connects the two intersection points
+ // Expression rewritten by Herbie, https://herbie.uwplse.org/demo/
+ // const double a = (d*d + r1*r1 - r2*r2) / (2*d);
+ const double a = std::fma(ra-rb, (ra + rb) / (d+d), d*0.5);
+
+ if (a == 0 || (d == 0 && ra == rb)) {
+ return cocircularArcsIntersectionPoint(ca, ra, a0, a2, aCCW, b0, b2, bCCW);
+ }
+
+ // Explicitly check endpoint intersections
+ if (a0.equals2D(b0) || a0.equals2D(b2)) {
+ return a0;
+ }
+ if (a2.equals2D(b0) || a2.equals2D(b2)) {
+ return a2;
+ }
+
+ // Compute interior intersection points.
+ const double dx = cb.x - ca.x;
+ const double dy = cb.y - ca.y;
+
+ // point where a line between the two circle center points intersects
+ // the radical line
+ CoordinateXY p{ca.x + a * dx/d, ca.y+a* dy/d};
+
+ // distance from p to the intersection points
+ const double h = std::sqrt(ra*ra - a*a);
+
+ CoordinateXY isect0{p.x + h* dy/d, p.y - h* dx/d };
+ CoordinateXY isect1{p.x - h* dy/d, p.y + h* dx/d };
+
+ if (containsPointOnCircle(ca, a0, a2, aCCW, isect0) &&
+ containsPointOnCircle(cb, b0, b2, bCCW, isect0)) {
+ return isect0;
+ }
+
+ if (containsPointOnCircle(ca, a0, a2, aCCW, isect1) &&
+ containsPointOnCircle(cb, b0, b2, bCCW, isect1)) {
+ return isect1;
+ }
+
+ return std::nullopt;
+}
+
}
}
diff --git a/src/algorithm/PointLocator.cpp b/src/algorithm/PointLocator.cpp
index 7fa7dd7e3..5545d7655 100644
--- a/src/algorithm/PointLocator.cpp
+++ b/src/algorithm/PointLocator.cpp
@@ -53,8 +53,8 @@ PointLocator::locate(const CoordinateXY& p, const Geometry* geom)
return locate(p, ls_geom);
}
- if (geomTypeId == GEOS_POLYGON) {
- const Polygon* poly_geom = static_cast<const Polygon*>(geom);
+ if (geomTypeId == GEOS_POLYGON || geomTypeId == GEOS_CURVEPOLYGON) {
+ const Surface* poly_geom = static_cast<const Surface*>(geom);
return locate(p, poly_geom);
}
@@ -185,32 +185,24 @@ PointLocator::locate(const CoordinateXY& p, const LineString* l)
/* private */
Location
-PointLocator::locateInPolygonRing(const CoordinateXY& p, const LinearRing* ring)
+PointLocator::locateInPolygonRing(const CoordinateXY& p, const Curve* ring)
{
if(!ring->getEnvelopeInternal()->intersects(p)) {
return Location::EXTERIOR;
}
- const CoordinateSequence* cl = ring->getCoordinatesRO();
-
- if(PointLocation::isOnLine(p, cl)) {
- return Location::BOUNDARY;
- }
- if(PointLocation::isInRing(p, cl)) {
- return Location::INTERIOR;
- }
- return Location::EXTERIOR;
+ return PointLocation::locateInRing(p, *ring);
}
/* private */
Location
-PointLocator::locate(const CoordinateXY& p, const Polygon* poly)
+PointLocator::locate(const CoordinateXY& p, const Surface* poly)
{
if(poly->isEmpty()) {
return Location::EXTERIOR;
}
- const LinearRing* shell = poly->getExteriorRing();
+ const Curve* shell = poly->getExteriorRing();
assert(shell);
Location shellLoc = locateInPolygonRing(p, shell);
@@ -223,7 +215,7 @@ PointLocator::locate(const CoordinateXY& p, const Polygon* poly)
// now test if the point lies in or on the holes
for(std::size_t i = 0, n = poly->getNumInteriorRing(); i < n; ++i) {
- const LinearRing* hole = poly->getInteriorRingN(i);
+ const Curve* hole = poly->getInteriorRingN(i);
Location holeLoc = locateInPolygonRing(p, hole);
if(holeLoc == Location::INTERIOR) {
return Location::EXTERIOR;
diff --git a/src/geom/CircularArc.cpp b/src/geom/CircularArc.cpp
index 5ac8f2980..1bca8155a 100644
--- a/src/geom/CircularArc.cpp
+++ b/src/geom/CircularArc.cpp
@@ -12,10 +12,13 @@
*
**********************************************************************/
+#include <geos/algorithm/CircularArcs.h>
#include <geos/geom/CircularArc.h>
#include <geos/triangulate/quadedge/TrianglePredicate.h>
#include <sstream>
+using geos::algorithm::CircularArcs;
+
namespace geos::geom {
CircularArc::CircularArc() :
@@ -192,29 +195,12 @@ CircularArc::~CircularArc()
}
bool
-CircularArc::containsAngle(double theta) const {
+CircularArc::containsAngle(double theta) const
+{
auto t0 = theta0();
auto t2 = theta2();
- if (theta == t0 || theta == t2) {
- return true;
- }
-
- if (getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE) {
- std::swap(t0, t2);
- }
-
- t2 -= t0;
- theta -= t0;
-
- if (t2 < 0) {
- t2 += 2*MATH_PI;
- }
- if (theta < 0) {
- theta += 2*MATH_PI;
- }
-
- return theta >= t2;
+ return CircularArcs::containsAngle(t0, t2, isCCW(), theta);
}
bool
@@ -237,6 +223,46 @@ CircularArc::containsPoint(const CoordinateXY& q) const
return containsPointOnCircle(q);
}
+CoordinateXY
+CircularArc::closestPoint(const CoordinateXY &q) const
+{
+ return CircularArcs::closestPointArcPoint(getCenter(), getRadius(), p0(), p2(), isCCW(), q);
+}
+
+std::array<CoordinateXY, 2>
+CircularArc::closestPoints(const CoordinateXY& q0, const CoordinateXY& q1) const
+{
+ return CircularArcs::closestPointsArcSegment(getCenter(), getRadius(), p0(), p2(), isCCW(), q0, q1);
+}
+
+std::array<CoordinateXY, 2>
+CircularArc::closestPoints(const CircularArc &other) const
+{
+ return CircularArcs::closestPointsArcArc(getCenter(), getRadius(), p0(), p2(), isCCW(),
+ other.getCenter(), other.getRadius(), other.p0(), other.p2(), other.isCCW());
+}
+
+double
+CircularArc::distance(const CoordinateXY &q) const
+{
+ return CircularArcs::distanceArcPoint(getCenter(), getRadius(), p0(), p2(), isCCW(), q);
+}
+
+double
+CircularArc::distance(const CoordinateXY& q0, const CoordinateXY& q1) const
+{
+ return CircularArcs::distanceArcSegment(getCenter(), getRadius(), p0(), p2(), isCCW(), q0, q1);
+}
+
+double
+CircularArc::distance(const CircularArc &other) const
+{
+ return CircularArcs::distanceArcArc(
+ getCenter(), getRadius(), p0(), p2(), isCCW(),
+ other.getCenter(), other.getRadius(), other.p0(), other.p2(), other.isCCW());
+}
+
+
double
CircularArc::getAngle() const
{
@@ -277,19 +303,22 @@ CircularArc::getArea() const {
}
CoordinateXY
-CircularArc::getDirectionPoint() const {
- return algorithm::CircularArcs::getDirectionPoint(getCenter(), getRadius(), theta0(), getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE);
+CircularArc::getDirectionPoint() const
+{
+ return CircularArcs::getDirectionPoint(getCenter(), getRadius(), theta0(), getOrientation() == algorithm::Orientation::COUNTERCLOCKWISE);
}
Envelope
-CircularArc::getEnvelope() const {
+CircularArc::getEnvelope() const
+{
Envelope env;
- algorithm::CircularArcs::expandEnvelope(env, p0(), p1(), p2());
+ CircularArcs::expandEnvelope(env, p0(), p1(), p2());
return env;
}
double
-CircularArc::getLength() const {
+CircularArc::getLength() const
+{
if (isLinear()) {
return p0().distance(p2());
}
@@ -298,7 +327,8 @@ CircularArc::getLength() const {
}
bool
-CircularArc::isUpwardAtPoint(const CoordinateXY& q) const {
+CircularArc::isUpwardAtPoint(const CoordinateXY& q) const
+{
auto quad = geom::Quadrant::quadrant(getCenter(), q);
bool isUpward;
diff --git a/src/operation/distance/ConnectedElementLocationFilter.cpp b/src/operation/distance/ConnectedElementLocationFilter.cpp
index b117627b1..d62339b53 100644
--- a/src/operation/distance/ConnectedElementLocationFilter.cpp
+++ b/src/operation/distance/ConnectedElementLocationFilter.cpp
@@ -47,10 +47,16 @@ void
ConnectedElementLocationFilter::filter_ro(const Geometry* geom)
{
if (geom->isEmpty()) return;
- if((typeid(*geom) == typeid(Point)) ||
- (typeid(*geom) == typeid(LineString)) ||
- (typeid(*geom) == typeid(LinearRing)) ||
- (typeid(*geom) == typeid(Polygon))) {
+
+ const auto type = geom->getGeometryTypeId();
+
+ if (type == GEOS_POINT ||
+ type == GEOS_LINESTRING ||
+ type == GEOS_LINEARRING ||
+ type == GEOS_POLYGON ||
+ type == GEOS_CIRCULARSTRING ||
+ type == GEOS_COMPOUNDCURVE ||
+ type == GEOS_CURVEPOLYGON) {
locations.emplace_back(geom, 0, *(geom->getCoordinate()));
}
}
@@ -60,10 +66,16 @@ ConnectedElementLocationFilter::filter_rw(Geometry* geom)
{
// empty geometries do not provide a location
if (geom->isEmpty()) return;
- if((typeid(*geom) == typeid(Point)) ||
- (typeid(*geom) == typeid(LineString)) ||
- (typeid(*geom) == typeid(LinearRing)) ||
- (typeid(*geom) == typeid(Polygon))) {
+
+ const auto type = geom->getGeometryTypeId();
+
+ if (type == GEOS_POINT ||
+ type == GEOS_LINESTRING ||
+ type == GEOS_LINEARRING ||
+ type == GEOS_POLYGON ||
+ type == GEOS_CIRCULARSTRING ||
+ type == GEOS_COMPOUNDCURVE ||
+ type == GEOS_CURVEPOLYGON) {
locations.emplace_back(geom, 0, *(geom->getCoordinate()));
}
}
diff --git a/src/operation/distance/DistanceOp.cpp b/src/operation/distance/DistanceOp.cpp
index 9a234b4c2..81907466d 100644
--- a/src/operation/distance/DistanceOp.cpp
+++ b/src/operation/distance/DistanceOp.cpp
@@ -23,6 +23,7 @@
#include <geos/operation/distance/GeometryLocation.h>
#include <geos/operation/distance/ConnectedElementLocationFilter.h>
#include <geos/algorithm/Distance.h>
+#include <geos/geom/CircularString.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/LineString.h>
@@ -30,8 +31,8 @@
#include <geos/geom/Polygon.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/LineSegment.h>
-#include <geos/geom/util/PolygonExtracter.h>
-#include <geos/geom/util/LinearComponentExtracter.h>
+#include <geos/geom/util/SimpleCurveExtracter.h>
+#include <geos/geom/util/SurfaceExtracter.h>
#include <geos/geom/util/PointExtracter.h>
#include <geos/util/IllegalArgumentException.h>
#include <geos/util.h>
@@ -113,9 +114,6 @@ DistanceOp::distance()
{
using geos::util::IllegalArgumentException;
- util::ensureNoCurvedComponents(geom[0]);
- util::ensureNoCurvedComponents(geom[1]);
-
if(geom[0] == nullptr || geom[1] == nullptr) {
throw IllegalArgumentException("null geometries are not supported");
}
@@ -208,14 +206,11 @@ DistanceOp::computeMinDistance()
void
DistanceOp::computeContainmentDistance()
{
- using geom::util::PolygonExtracter;
-
- Polygon::ConstVect polys1;
- PolygonExtracter::getPolygons(*(geom[1]), polys1);
-
+ std::vector<const Surface*> polys1;
+ geom::util::SurfaceExtracter::getSurfaces(*(geom[1]), polys1);
#if GEOS_DEBUG
- std::cerr << "PolygonExtracter found " << polys1.size() << " polygons in geometry 2" << std::endl;
+ std::cerr << "SurfaceExtracter found " << polys1.size() << " polygons in geometry 2" << std::endl;
#endif
// NOTE:
@@ -239,11 +234,11 @@ DistanceOp::computeContainmentDistance()
}
}
- Polygon::ConstVect polys0;
- PolygonExtracter::getPolygons(*(geom[0]), polys0);
+ std::vector<const Surface*> polys0;
+ geom::util::SurfaceExtracter::getSurfaces(*(geom[0]), polys0);
#if GEOS_DEBUG
- std::cerr << "PolygonExtracter found " << polys0.size() << " polygons in geometry 1" << std::endl;
+ std::cerr << "SurfaceExtracter found " << polys0.size() << " polygons in geometry 1" << std::endl;
#endif
@@ -267,7 +262,7 @@ DistanceOp::computeContainmentDistance()
/*private*/
void
DistanceOp::computeInside(std::vector<GeometryLocation> & locs,
- const Polygon::ConstVect& polys,
+ const std::vector<const Surface*>& polys,
std::array<GeometryLocation, 2> & locPtPoly)
{
for(auto& loc : locs) {
@@ -276,7 +271,7 @@ DistanceOp::computeInside(std::vector<GeometryLocation> & locs,
if (Location::EXTERIOR != ptLocator.locate(pt, static_cast<const Geometry*>(poly))) {
minDistance = 0.0;
- locPtPoly[0] = std::move(loc);
+ locPtPoly[0] = loc;
locPtPoly[1] = GeometryLocation(poly, pt);
return;
}
@@ -288,7 +283,7 @@ DistanceOp::computeInside(std::vector<GeometryLocation> & locs,
void
DistanceOp::computeFacetDistance()
{
- using geom::util::LinearComponentExtracter;
+ using geom::util::SimpleCurveExtracter;
using geom::util::PointExtracter;
std::array<GeometryLocation, 2> locGeom;
@@ -297,13 +292,13 @@ DistanceOp::computeFacetDistance()
* Geometries are not wholly inside, so compute distance from lines
* and points of one to lines and points of the other
*/
- LineString::ConstVect lines0;
- LineString::ConstVect lines1;
- LinearComponentExtracter::getLines(*(geom[0]), lines0);
- LinearComponentExtracter::getLines(*(geom[1]), lines1);
+ std::vector<const SimpleCurve*> lines0;
+ std::vector<const SimpleCurve*> lines1;
+ SimpleCurveExtracter::getCurves(*(geom[0]), lines0);
+ SimpleCurveExtracter::getCurves(*(geom[1]), lines1);
#if GEOS_DEBUG
- std::cerr << "LinearComponentExtracter found "
+ std::cerr << "SimpleCurveExtracter found "
<< lines0.size() << " lines in geometry 1 and "
<< lines1.size() << " lines in geometry 2 "
<< std::endl;
@@ -368,17 +363,36 @@ DistanceOp::computeFacetDistance()
/*private*/
void
DistanceOp::computeMinDistanceLines(
- const LineString::ConstVect& lines0,
- const LineString::ConstVect& lines1,
+ const std::vector<const SimpleCurve*>& lines0,
+ const std::vector<const SimpleCurve*>& lines1,
std::array<GeometryLocation, 2> & locGeom)
{
- for(const LineString* line0 : lines0) {
- for(const LineString* line1 : lines1) {
+ for(const SimpleCurve* line0 : lines0) {
+ const bool isCurved0 = line0->getGeometryTypeId() == GEOS_CIRCULARSTRING;
+
+ for(const SimpleCurve* line1 : lines1) {
+ const bool isCurved1 = line1->getGeometryTypeId() == GEOS_CIRCULARSTRING;
if (line0->isEmpty() || line1->isEmpty())
continue;
- computeMinDistance(line0, line1, locGeom);
+ if (isCurved0) {
+ const CircularString* cs0 = detail::down_cast<const CircularString*>(line0);
+ if (isCurved1) {
+ const CircularString* cs1 = detail::down_cast<const CircularString*>(line1);
+ computeMinDistance(cs0, cs1, locGeom);
+ } else {
+ computeMinDistance(cs0, detail::down_cast<const LineString*>(line1), locGeom);
+ }
+ } else {
+ if (isCurved1) {
+ const CircularString* cs1 = detail::down_cast<const CircularString*>(line1);
+ computeMinDistance(detail::down_cast<const LineString*>(line0), cs1, locGeom);
+ } else {
+ computeMinDistance(detail::down_cast<const LineString*>(line0), detail::down_cast<const LineString*>(line1), locGeom);
+ }
+ }
+
if(minDistance <= terminateDistance) {
return;
}
@@ -426,17 +440,23 @@ DistanceOp::computeMinDistancePoints(
/*private*/
void
DistanceOp::computeMinDistanceLinesPoints(
- const LineString::ConstVect& lines,
- const Point::ConstVect& points,
+ const std::vector<const SimpleCurve*>& lines,
+ const std::vector<const Point*>& points,
std::array<GeometryLocation, 2> & locGeom)
{
- for(const LineString* line : lines) {
+ for(const SimpleCurve* sc : lines) {
for(const Point* pt : points) {
-
- if (line->isEmpty() || pt->isEmpty())
+ if (sc->isEmpty() || pt->isEmpty())
continue;
- computeMinDistance(line, pt, locGeom);
+ if (sc->getGeometryTypeId() == GEOS_CIRCULARSTRING) {
+ const CircularString* cs = detail::down_cast<const CircularString*>(sc);
+ computeMinDistance(cs, pt, locGeom);
+ } else {
+ const LineString* line = detail::down_cast<const LineString*>(sc);
+ computeMinDistance(line, pt, locGeom);
+ }
+
if(minDistance <= terminateDistance) {
return;
}
@@ -501,6 +521,113 @@ DistanceOp::computeMinDistance(
}
}
+/*private*/
+void
+DistanceOp::computeMinDistance(
+ const CircularString* cs0,
+ const LineString* line1,
+ std::array<GeometryLocation, 2> & locGeom)
+{
+ const Envelope* lineEnv0 = cs0->getEnvelopeInternal();
+ const Envelope* lineEnv1 = line1->getEnvelopeInternal();
+ if(lineEnv0->distance(*lineEnv1) > minDistance) {
+ return;
+ }
+
+ const CoordinateSequence* coord1 = line1->getCoordinatesRO();
+ std::size_t npts1 = coord1->getSize();
+
+ // brute force approach!
+ for(const auto& arc : cs0->getArcs()) {
+ const Envelope arcEnv = arc.getEnvelope();
+
+ if (arcEnv.distanceSquared(*lineEnv1) > minDistance*minDistance) {
+ continue;
+ }
+
+ for(std::size_t j = 0; j < npts1 - 1; ++j) {
+ const CoordinateXY& q0 = coord1->getAt<CoordinateXY>(j);
+ const CoordinateXY& q1 = coord1->getAt<CoordinateXY>(j+1);
+
+ Envelope segEnv(q0, q1);
+
+ if (arcEnv.distanceSquared(segEnv) > minDistance*minDistance) {
+ continue;
+ }
+
+ const auto closestPts = arc.closestPoints(q0, q1);
+
+ const double dist = closestPts[0].distance(closestPts[1]);
+ if(dist < minDistance) {
+ minDistance = dist;
+
+ locGeom[0] = GeometryLocation(cs0, arc.getCoordinatePosition(), closestPts[0]);
+ locGeom[1] = GeometryLocation(line1, j, closestPts[1]);
+ }
+
+ if(minDistance <= terminateDistance) {
+ return;
+ }
+ }
+ }
+}
+
+/*private*/
+void
+DistanceOp::computeMinDistance(
+ const LineString* line0,
+ const CircularString* cs1,
+ std::array<GeometryLocation, 2> & locGeom)
+{
+ computeMinDistance(cs1, line0, locGeom);
+ std::swap(locGeom[0], locGeom[1]);
+}
+
+/*private*/
+void
+DistanceOp::computeMinDistance(
+ const CircularString* cs0,
+ const CircularString* cs1,
+ std::array<GeometryLocation, 2> & locGeom)
+{
+ const Envelope* geomEnv0 = cs0->getEnvelopeInternal();
+ const Envelope* geomEnv1 = cs1->getEnvelopeInternal();
+
+ if(geomEnv0->distance(*geomEnv1) > minDistance) {
+ return;
+ }
+
+ // brute force approach!
+ for(const auto& arc0 : cs0->getArcs()) {
+ const Envelope arcEnv0 = arc0.getEnvelope();
+
+ if (arcEnv0.distanceSquared(*geomEnv1) > minDistance*minDistance) {
+ continue;
+ }
+
+ for (const auto& arc1 : cs1->getArcs()) {
+ const Envelope arcEnv1 = arc1.getEnvelope();
+
+ if (arcEnv0.distanceSquared(arcEnv1) > minDistance*minDistance) {
+ continue;
+ }
+
+ const auto closestPoints = arc0.closestPoints(arc1);
+
+ const double dist = closestPoints[0].distance(closestPoints[1]);
+ if(dist < minDistance) {
+ minDistance = dist;
+ locGeom[0] = GeometryLocation(cs0, arc0.getCoordinatePosition(), closestPoints[0]);
+ locGeom[1] = GeometryLocation(cs1, arc1.getCoordinatePosition(), closestPoints[1]);
+ }
+
+ if(minDistance <= terminateDistance) {
+ return;
+ }
+ }
+ }
+}
+
/*private*/
void
DistanceOp::computeMinDistance(const LineString* line,
@@ -538,6 +665,36 @@ DistanceOp::computeMinDistance(const LineString* line,
}
}
+/*private*/
+void
+DistanceOp::computeMinDistance(const CircularString* cs,
+ const Point* pt,
+ std::array<GeometryLocation, 2> & locGeom)
+{
+ const Envelope* env0 = cs->getEnvelopeInternal();
+ const Envelope* env1 = pt->getEnvelopeInternal();
+ if(env0->distance(*env1) > minDistance) {
+ return;
+ }
+
+ // brute force approach!
+ const CoordinateXY& coord = *pt->getCoordinate();
+ for (const auto& arc : cs->getArcs()) {
+ CoordinateXY closestPt = arc.closestPoint(coord);
+ const double dist = coord.distance(closestPt);
+
+ if(dist < minDistance) {
+ minDistance = dist;
+ locGeom[0] = GeometryLocation(cs, arc.getCoordinatePosition(), closestPt);
+ locGeom[1] = GeometryLocation(pt, 0, coord);
+ }
+
+ if(minDistance <= terminateDistance) {
+ return;
+ }
+ }
+}
+
/* public static */
bool
DistanceOp::isWithinDistance(const geom::Geometry& g0,
diff --git a/tests/unit/algorithm/CircularArcsTest.cpp b/tests/unit/algorithm/CircularArcsTest.cpp
index d1296a32f..cb1da46c8 100644
--- a/tests/unit/algorithm/CircularArcsTest.cpp
+++ b/tests/unit/algorithm/CircularArcsTest.cpp
@@ -4,6 +4,9 @@
#include <geos/algorithm/CircularArcs.h>
#include <geos/geom/CircularArc.h>
#include <geos/geom/Quadrant.h>
+#include <geos/algorithm/Orientation.h>
+
+#include "utility.h"
using geos::geom::CircularArc;
using geos::geom::CoordinateXY;
@@ -20,7 +23,7 @@ struct test_circulararcs_data {
using XY = CoordinateXY;
void checkEnvelope(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2,
- double xmin, double ymin, double xmax, double ymax)
+ double xmin, double ymin, double xmax, double ymax) const
{
{
Envelope e;
@@ -57,6 +60,29 @@ struct test_circulararcs_data {
}
}
+ static void
+ checkArcIntersectionPoint(const CircularArc& a, const CircularArc& b, const CoordinateXY& expected, double tol = 0)
+ {
+ auto actual = CircularArcs::arcIntersectionPoint(a.getCenter(), a.getRadius(), a.p0(), a.p2(), a.isCCW(),
+ b.getCenter(), b.getRadius(), b.p0(), b.p2(), b.isCCW());
+
+ ensure(actual.has_value());
+
+ ensure_equals_xy(actual.value(), expected, tol);
+ }
+
+ static void
+ checkArcDisjoint(const CircularArc& a, const CircularArc& b)
+ {
+ auto actual = CircularArcs::arcIntersectionPoint(a.getCenter(), a.getRadius(), a.p0(), a.p2(), a.isCCW(),
+ b.getCenter(), b.getRadius(), b.p0(), b.p2(), b.isCCW());
+
+ ensure(!actual.has_value());
+ }
+
+ const CircularArc unit_semi = CircularArc::create(XY{-1, 0},
+ XY{0, 1},
+ XY{1, 0});
};
using group = test_group<test_circulararcs_data>;
@@ -328,5 +354,605 @@ void object::test<18>() {
ensure_equals("point is not expected distance from p0", arc2.p0().distance(dp2), arc2.getRadius(), 1e-5);
}
+template<>
+template<>
+void object::test<20>()
+{
+ set_test_name("pt-arc distance, point within unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{0, 0.5}), 0.5);
+}
+
+template<>
+template<>
+void object::test<21>()
+{
+ set_test_name("pt-arc distance, point outside unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{0, 1.5}), 0.5);
+}
+
+template<>
+template<>
+void object::test<22>()
+{
+ set_test_name("pt-arc distance, point outside unit semicircle, sqrt(2) units from endpoint");
+
+ ensure_equals(unit_semi.distance(XY{0, -1}), std::sqrt(2));
+}
+
+template<>
+template<>
+void object::test<23>()
+{
+ set_test_name("pt-arc distance, point outside unit semicircle, sqrt(2)-1 units from endpoint");
+
+ ensure_equals(unit_semi.distance(XY{1, 1}), std::sqrt(2)-1);
+}
+
+template<>
+template<>
+void object::test<24>()
+{
+ set_test_name("pt-arc distance, point on unit semicircle midpoint");
+
+ ensure_equals(unit_semi.distance(XY{0, 1}), 0);
+}
+
+template<>
+template<>
+void object::test<25>()
+{
+ set_test_name("pt-arc distance, point on unit semicircle endpoint");
+
+ ensure_equals(unit_semi.distance(XY{1, 0}), 0);
+}
+
+template<>
+template<>
+void object::test<26>()
+{
+ set_test_name("pt-arc distance, point on unit semicircle center");
+
+ ensure_equals(unit_semi.distance(XY{0, 0}), 1);
+}
+
+template<>
+template<>
+void object::test<27>()
+{
+ set_test_name("pt-arc distance, point inside closed circle");
+
+ ensure_equals(unit_semi.distance(XY{0, 0.5}), 0.5);
+}
+
+template<>
+template<>
+void object::test<28>()
+{
+ set_test_name("arc-seg distance, horizontal edge above unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{-2, 2}, XY{2, 2}), 1.0);
+}
+
+template<>
+template<>
+void object::test<29>()
+{
+ set_test_name("arc-seg distance, vertical edge to the right of unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{2, -2}, XY{2, 2}), 1.0);
+}
+
+template<>
+template<>
+void object::test<30>()
+{
+ set_test_name("arc-seg distance, vertical edge to the left of unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{-2, -2}, XY{-2, 2}), 1.0);
+}
+
+template<>
+template<>
+void object::test<31>()
+{
+ set_test_name("arc-seg distance, vertical edge within unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{0, 0}, XY{0, 0.5}), 0.5);
+}
+
+template<>
+template<>
+void object::test<32>()
+{
+ set_test_name("arc-seg distance, horizontal edge grazing unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{-2, 1}, XY{2, 1}), 0);
+}
+
+template<>
+template<>
+void object::test<33>()
+{
+ set_test_name("arc-seg distance, horizontal edge upper-right of unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{1, 1}, XY{2, 1}), std::sqrt(2) - 1);
+}
+
+template<>
+template<>
+void object::test<34>()
+{
+ set_test_name("arc-seg distance, edge intersecting unit semicircle");
+
+ ensure_equals(unit_semi.distance(XY{0, 0}, XY{2, 2}), 0);
+}
+
+template<>
+template<>
+void object::test<35>()
+{
+ set_test_name("arc-seg distance, edge upper-left of unit semicircle");
+
+ ensure_equals("distance is incorrect", unit_semi.distance(XY{-1, 1}, XY{-2, 2}), std::sqrt(2) - 1, 1e-6);
+}
+
+template<>
+template<>
+void object::test<36>()
+{
+ set_test_name("arc-arc distance, PostGIS ticket #4326");
+
+ CircularArc a = CircularArc::create(
+ XY{-1, 4},
+ XY{0, 5},
+ XY{1, 4});
+
+ CircularArc b = CircularArc::create(
+ XY{1, 6},
+ XY{6, 1},
+ XY{9, 7});
+
+ ensure_equals("distance is incorrect", a.distance(b), 0.0475666, 1e-6);
+}
+
+template<>
+template<>
+void object::test<37>()
+{
+ set_test_name("arc-arc distance, arc inside unit semicircle");
+
+ CircularArc a = CircularArc::create(
+ XY{0, 0.5},
+ XY{-0.3, 0.5},
+ XY{0, 0.6});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 0.271798, 1e-6);
+}
+
+template<>
+template<>
+void object::test<38>()
+{
+ set_test_name("arc-arc distance, arc inside unit semicircle (2)");
+
+ CircularArc a = CircularArc::create(
+ XY{-0.5, 0.5},
+ XY{-0.4, 0.2},
+ XY{0, 0});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 0.292893, 1e-6);
+}
+
+template<>
+template<>
+void object::test<39>()
+{
+ set_test_name("arc-arc distance, arc above unit semicircle");
+
+ CircularArc a = CircularArc::create(
+ XY{-1, 3},
+ XY{0, 2},
+ XY{1, 3});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 1, 1e-6);
+}
+
+template<>
+template<>
+void object::test<40>()
+{
+ set_test_name("arc-arc distance, arc grazes unit semicircle");
+
+ CircularArc a = CircularArc::create(
+ XY{-1, 2},
+ XY{0, 1},
+ XY{1, 2});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<41>()
+{
+ set_test_name("arc-arc distance, circles intersect but arcs do not");
+
+ CircularArc a = CircularArc::create(
+ XY{-1, 1},
+ XY{0, 2},
+ XY{1, 1});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), std::sqrt(2) - 1, 1e-6);
+}
+
+template<>
+template<>
+void object::test<42>()
+{
+ set_test_name("arc-arc distance, circles and arcs intersect");
+
+ CircularArc a = CircularArc::create(
+ XY{-1, 1},
+ XY{0, 0},
+ XY{1, 1});
+
+ ensure_equals(unit_semi.distance(a), 0);
+}
+
+template<>
+template<>
+void object::test<43>()
+{
+ set_test_name("arc-arc distance, closed circle");
+
+ CircularArc a = CircularArc::create(
+ XY{-2, -0.1},
+ XY{1.5, -0.1},
+ XY{-2, -0.1});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 0.480742, 1e-6);
+}
+
+template<>
+template<>
+void object::test<44>()
+{
+ set_test_name("arc-arc distance, concentric and fully parallel");
+
+ CircularArc a = CircularArc::create(
+ XY{-2, 0},
+ XY{0, 2},
+ XY{2, 0});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 1.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<45>()
+{
+ set_test_name("arc-arc distance, concentric with A fully included in unit semicircle's range");
+
+ CircularArc a = CircularArc::create(
+ XY{-0.5 / std::sqrt(2), 0.5 / std::sqrt(2)},
+ XY{0, 0.5},
+ XY{0.5 / std::sqrt(2), 0.5 / std::sqrt(2)});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 0.5, 1e-6);
+}
+
+template<>
+template<>
+void object::test<46>()
+{
+ set_test_name("arc-arc distance, concentric with A partially included in unit semicircle's range");
+
+ CircularArc a = CircularArc::create(
+ XY{-0.5 / std::sqrt(2), -0.5 / std::sqrt(2)},
+ XY{-0.5, 0},
+ XY{-0.5 / std::sqrt(2), 0.5 / std::sqrt(2)});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 0.5, 1e-6);
+}
+
+template<>
+template<>
+void object::test<47>()
+{
+ set_test_name("arc-arc distance, concentric without parallel segments");
+
+ CircularArc a = CircularArc::create(
+ XY{-0.5 / std::sqrt(2), -0.5 / std::sqrt(2)},
+ XY{0, -0.5},
+ XY{0.5 / std::sqrt(2), -0.5 / std::sqrt(2)});
+
+ ensure_equals("distance is incorrect", unit_semi.distance(a), 0.736813, 1e-6);
+}
+
+template<>
+template<>
+void object::test<48>()
+{
+ set_test_name("arc-arc distance, arcs are the same");
+
+ ensure_equals(unit_semi.distance(unit_semi), 0);
+}
+
+template<>
+template<>
+void object::test<49>()
+{
+ set_test_name("arc-arc distance, different orientations (1)");
+
+ CircularArc a = CircularArc::create(XY{-22, 0}, XY{-17, -5}, XY{-12, 0});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 2.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<50>()
+{
+ set_test_name("arc-arc distance, different orientations (2)");
+
+ CircularArc a = CircularArc::create(XY{-19, 0}, XY{-14, -5}, XY{-9, 0});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 1.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<51>()
+{
+ set_test_name("arc-arc distance, different orientations (3)");
+
+ CircularArc a = CircularArc::create(XY{-9, 0}, XY{-4, -5}, XY{1, 0});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 1.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<52>()
+{
+ set_test_name("arc-arc distance, different orientations (4)");
+
+ CircularArc a = CircularArc::create(XY{-1, 0}, XY{4, -5}, XY{9, 0});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 1.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<53>()
+{
+ set_test_name("arc-arc distance, different orientations (5)");
+
+ CircularArc a = CircularArc::create(XY{1, 0}, XY{6, -5}, XY{11, 0});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 1.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<54>()
+{
+ set_test_name("arc-arc distance, different orientations (6)");
+
+ CircularArc a = CircularArc::create(XY{11, 0}, XY{16, -5}, XY{21, 0});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 1.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<55>()
+{
+ set_test_name("arc-arc distance, different orientations (7)");
+
+ CircularArc a = CircularArc::create(XY{-15, -6}, XY{-10, -1}, XY{-5, -6});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 1.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<56>()
+{
+ set_test_name("arc-arc distance, different orientations (8)");
+
+ CircularArc a = CircularArc::create(XY{-5, 0}, XY{0, 5}, XY{5, 0});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 5.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<57>()
+{
+ set_test_name("arc-arc distance, different orientations (9)");
+
+ CircularArc a = CircularArc::create(XY{-5, 0}, XY{0, -5}, XY{5, 0});
+ CircularArc b = CircularArc::create(XY{-10, 0}, XY{0, 10}, XY{10, 0});
+
+ ensure_equals("distance is incorrect", a.distance(b), 5.0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<58>()
+{
+ set_test_name("arc-seg distance, closest point is segment interior");
+
+ XY p0{1, -1};
+ XY p1{2, 0};
+
+ ensure_equals("distance is incorrect", unit_semi.distance(p0, p1), std::sqrt(2)/2, 1e-6);
+}
+
+template<>
+template<>
+void object::test<59>()
+{
+ set_test_name("arc-seg distance, closest distance is between segment interior and arc interior");
+
+ CircularArc a = CircularArc::create(XY{-5, 0}, XY{0, 5}, XY{4, -3});
+
+ XY p0{4, -5};
+ XY p1{7, 0};
+
+ ensure_equals("distance is incorrect", a.distance(p0, p1), 1.0024502, 1e-6);
+}
+
+template<>
+template<>
+void object::test<60>()
+{
+ set_test_name("arc-seg distance (vertical segment tangent to arc)");
+
+ CircularArc a = CircularArc::create(XY{5, 1}, XY{5, 3}, XY{5, 2}, 1, geos::algorithm::Orientation::CLOCKWISE);
+
+ XY q0{4, 4};
+ XY q1{4, 0};
+
+ ensure_equals("distance is incorrect", a.distance(q0, q1), 0);
+}
+
+template<>
+template<>
+void object::test<61>()
+{
+ set_test_name("intersection point: arcs from disjoint circles");
+
+ CircularArc a = CircularArc::create(XY{-1, 0}, XY{0, 1}, XY{1, 0});
+ CircularArc b = CircularArc::create(XY{5, 0}, XY{6, 1}, XY{7, 0});
+
+ checkArcDisjoint(a, b);
+}
+
+template<>
+template<>
+void object::test<62>()
+{
+ set_test_name("intersection point: arcs from contained circles");
+
+ CircularArc a = CircularArc::create(XY{-1, 0}, XY{0, 1}, XY{1, 0});
+ CircularArc b = CircularArc::create(XY{-5, 0}, XY{0, 5}, XY{5, 0});
+
+ checkArcDisjoint(a, b);
+}
+
+template<>
+template<>
+void object::test<63>()
+{
+ set_test_name("intersection point: arcs with endpoint intersections");
+
+ CircularArc a = CircularArc::create(XY{-1, 0}, XY{0, 1}, XY{1, 0});
+ CircularArc b = CircularArc::create(XY{1, 0}, XY{2, -1}, XY{3, 0});
+ CircularArc c = CircularArc::create(XY{3, 0}, XY{4, 1}, XY{5, 0});
+
+ checkArcIntersectionPoint(a, b, XY{1, 0});
+ checkArcIntersectionPoint(b, a, XY{1, 0});
+
+ checkArcIntersectionPoint(b, c, XY{3, 0});
+ checkArcIntersectionPoint(c, b, XY{3, 0});
+
+ checkArcDisjoint(a, c);
+}
+
+template<>
+template<>
+void object::test<64>()
+{
+ set_test_name("intersection point: cocircular arcs with endpoint intersection");
+
+ CircularArc a = CircularArc::create(XY{-5, 0}, XY{0, 5}, XY{5, 0});
+ CircularArc b = CircularArc::create(XY{-5, 0}, XY{0, -5}, XY{ 4, -3});
+
+ checkArcIntersectionPoint(a, b, XY{-5, 0}, 1e-8);
+ checkArcIntersectionPoint(b, a, XY{-5, 0}, 1e-8);
+}
+
+template<>
+template<>
+void object::test<65>()
+{
+ set_test_name("intersection point: cocircular arcs with overlap");
+
+ CircularArc a = CircularArc::create(XY{-5, 0}, XY{0, 5}, XY{5, 0});
+ CircularArc b = CircularArc::create({0, -5}, XY{ 4, -3}, XY{4, 3});
+
+ checkArcIntersectionPoint(a, b, XY{4, 3}, 1e-8);
+ checkArcIntersectionPoint(b, a, XY{5, 0}, 1e-8);
+}
+
+template<>
+template<>
+void object::test<66>()
+{
+ set_test_name("intersection point: disjoint cocircular arcs");
+
+ CircularArc a = CircularArc::create(XY{-4, 3}, XY{0, 5}, XY{4, 3});
+ CircularArc b = CircularArc::create({-5, 0}, XY{ 0, -5}, XY{5, 0});
+
+ checkArcDisjoint(a, b);
+ checkArcDisjoint(b, a);
+}
+
+template<>
+template<>
+void object::test<67>()
+{
+ set_test_name("closest points: two collapsed arcs");
+
+ CircularArc a = CircularArc::create(XY{-10, 0}, XY{-9, 0}, XY{-8, 0}); // horizontal line
+ CircularArc b = CircularArc::create({0, -5}, XY{ 0, 0}, XY{0, 5}); // vertical line
+
+ auto pts = a.closestPoints(b);
+ ensure_equals_xy(pts[0], XY{-8, 0});
+ ensure_equals_xy(pts[1], XY{0, 0});
+}
+
+template<>
+template<>
+void object::test<68>()
+{
+ set_test_name("closest points: arc and collapsed arc");
+
+ CircularArc a = CircularArc::create(XY{-10, -10}, XY{-9, -10}, XY{10, -10}); // horizontal line
+ CircularArc b = CircularArc::create({-5, 5}, XY{ -4, 2}, XY{5, 5});
+
+ auto pts = a.closestPoints(b);
+ ensure_equals_xy(pts[0], XY{0, -10});
+ ensure_equals_xy(pts[1], XY{0, 0});
+
+ pts = b.closestPoints(a);
+ ensure_equals_xy(pts[0], XY{0, 0});
+ ensure_equals_xy(pts[1], XY{0, -10});
+}
+
+template<>
+template<>
+void object::test<69>()
+{
+ set_test_name("distance: overlapping cocircular arcs");
+
+ CircularArc a = CircularArc::create(XY{-5, 0}, XY{0, 5}, XY{5, 0});
+ CircularArc b = CircularArc::create({-4, 3}, XY{0, 5}, XY{ 4, 3});
+
+ ensure_equals(a.distance(b), 0);
+ ensure_equals(b.distance(a), 0);
+}
}
diff --git a/tests/unit/capi/GEOSDistanceTest.cpp b/tests/unit/capi/GEOSDistanceTest.cpp
index 42a895a62..a7e2f3bf4 100644
--- a/tests/unit/capi/GEOSDistanceTest.cpp
+++ b/tests/unit/capi/GEOSDistanceTest.cpp
@@ -171,15 +171,16 @@ template<>
template<>
void object::test<6>()
{
- geom1_ = fromWKT("CIRCULARSTRING (0 0, 1 1, 2 0)");
- geom2_ = fromWKT("LINESTRING (1 1.0001, 2 1)");
+ set_test_name("curved inputs");
- ensure(geom1_);
- ensure(geom2_);
+ geom1_ = fromWKT("CIRCULARSTRING (0 0, 1 1, 2 0)");
+ geom2_ = fromWKT("LINESTRING (1 2, 2 2)");
double dist;
int ret = GEOSDistance(geom1_, geom2_, &dist);
- ensure_equals("curved geometry not supported", ret, 0);
+ ensure_equals(ret, 1);;
+
+ ensure_equals(dist, 1.0);
}
template<>
diff --git a/tests/unit/capi/GEOSDistanceWithinTest.cpp b/tests/unit/capi/GEOSDistanceWithinTest.cpp
index 074c1b59b..0b164d199 100644
--- a/tests/unit/capi/GEOSDistanceWithinTest.cpp
+++ b/tests/unit/capi/GEOSDistanceWithinTest.cpp
@@ -412,14 +412,13 @@ template<>
template<>
void object::test<32>()
{
+ set_test_name("curved inputs");
+
geom1_ = fromWKT("CIRCULARSTRING (0 0, 1 1, 2 0)");
geom2_ = fromWKT("LINESTRING (1 1.0001, 2 1)");
- ensure(geom1_);
- ensure(geom2_);
-
char ret = GEOSDistanceWithin(geom1_, geom2_, 0.1);
- ensure_equals("curved geometry not supported", ret, 2);
+ ensure_equals(ret, 1);
}
template<>
diff --git a/tests/unit/geom/CircularStringTest.cpp b/tests/unit/geom/CircularStringTest.cpp
index 6468696da..24a81086a 100644
--- a/tests/unit/geom/CircularStringTest.cpp
+++ b/tests/unit/geom/CircularStringTest.cpp
@@ -79,7 +79,7 @@ void object::test<2>()
{
// Geometry type functions
ensure_equals("getGeometryType", cs_->getGeometryType(), "CircularString");
- ensure_equals("getGeometryTypdId", cs_->getGeometryTypeId(), geos::geom::GEOS_CIRCULARSTRING);
+ ensure_equals("getGeometryTypeId", cs_->getGeometryTypeId(), geos::geom::GEOS_CIRCULARSTRING);
ensure("isCollection", !cs_->isCollection());
// Geometry size functions
@@ -144,8 +144,8 @@ void object::test<3>()
ensure_THROW(cs_->symDifference(cs_.get()), geos::util::UnsupportedOperationException);
// Distance
- ensure_THROW(cs_->distance(cs_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cs_->isWithinDistance(cs_.get(), 1), geos::util::UnsupportedOperationException);
+ ensure_equals(cs_->distance(cs_.get()), 0);
+ ensure(cs_->isWithinDistance(cs_.get(), 1));
// Valid / Simple
ensure_THROW(cs_->isSimple(), geos::util::UnsupportedOperationException);
diff --git a/tests/unit/geom/CompoundCurveTest.cpp b/tests/unit/geom/CompoundCurveTest.cpp
index c7def3e63..574c41ae6 100644
--- a/tests/unit/geom/CompoundCurveTest.cpp
+++ b/tests/unit/geom/CompoundCurveTest.cpp
@@ -167,8 +167,8 @@ void object::test<3>()
ensure_THROW(cc_->symDifference(cc_.get()), geos::util::UnsupportedOperationException);
// Distance
- ensure_THROW(cc_->distance(cc_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cc_->isWithinDistance(cc_.get(), 1), geos::util::UnsupportedOperationException);
+ ensure_equals(cc_->distance(cc_.get()), 0);
+ ensure(cc_->isWithinDistance(cc_.get(), 1));
// Valid / Simple
ensure_THROW(cc_->isSimple(), geos::util::UnsupportedOperationException);
diff --git a/tests/unit/geom/CurvePolygonTest.cpp b/tests/unit/geom/CurvePolygonTest.cpp
index 0e60aded8..c41dab112 100644
--- a/tests/unit/geom/CurvePolygonTest.cpp
+++ b/tests/unit/geom/CurvePolygonTest.cpp
@@ -172,8 +172,8 @@ void object::test<3>()
ensure_THROW(cp_->symDifference(cp_.get()), geos::util::UnsupportedOperationException);
// Distance
- ensure_THROW(cp_->distance(cp_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(cp_->isWithinDistance(cp_.get(), 1), geos::util::UnsupportedOperationException);
+ ensure_equals(cp_->distance(cp_.get()), 0);
+ ensure(cp_->isWithinDistance(cp_.get(), 1));
// Valid / Simple
ensure_THROW(cp_->isSimple(), geos::util::UnsupportedOperationException);
diff --git a/tests/unit/geom/MultiCurveTest.cpp b/tests/unit/geom/MultiCurveTest.cpp
index 365e77e7b..70b0ada4d 100644
--- a/tests/unit/geom/MultiCurveTest.cpp
+++ b/tests/unit/geom/MultiCurveTest.cpp
@@ -159,8 +159,8 @@ void object::test<3>()
ensure_THROW(mc_->symDifference(mc_.get()), geos::util::UnsupportedOperationException);
// Distance
- ensure_THROW(mc_->distance(mc_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(mc_->isWithinDistance(mc_.get(), 1), geos::util::UnsupportedOperationException);
+ ensure_equals(mc_->distance(mc_.get()), 0);
+ ensure(mc_->isWithinDistance(mc_.get(), 1));
// Valid / Simple
ensure_THROW(mc_->isSimple(), geos::util::UnsupportedOperationException);
diff --git a/tests/unit/geom/MultiSurfaceTest.cpp b/tests/unit/geom/MultiSurfaceTest.cpp
index 1b1d3509a..f12fb9e38 100644
--- a/tests/unit/geom/MultiSurfaceTest.cpp
+++ b/tests/unit/geom/MultiSurfaceTest.cpp
@@ -148,8 +148,8 @@ void object::test<3>()
ensure_THROW(ms_->symDifference(ms_.get()), geos::util::UnsupportedOperationException);
// Distance
- ensure_THROW(ms_->distance(ms_.get()), geos::util::UnsupportedOperationException);
- ensure_THROW(ms_->isWithinDistance(ms_.get(), 1), geos::util::UnsupportedOperationException);
+ ensure_equals(ms_->distance(ms_.get()), 0);
+ ensure(ms_->isWithinDistance(ms_.get(), 1));
// Valid / Simple
ensure_THROW(ms_->isSimple(), geos::util::UnsupportedOperationException);
diff --git a/tests/unit/operation/distance/DistanceOpTest.cpp b/tests/unit/operation/distance/DistanceOpTest.cpp
index bbb5e0850..7c2686773 100644
--- a/tests/unit/operation/distance/DistanceOpTest.cpp
+++ b/tests/unit/operation/distance/DistanceOpTest.cpp
@@ -16,11 +16,18 @@
#include <geos/io/WKTReader.h>
#include <geos/io/WKBReader.h>
#include <geos/geom/CoordinateSequence.h>
+#include <geos/operation/distance/DistanceOp.h>
#include <geos/util.h>
// std
#include <memory>
#include <string>
+#include "utility.h"
+
+using XY = geos::geom::CoordinateXY;
+using geos::operation::distance::DistanceOp;
+using geos::operation::distance::GeometryLocation;
+
namespace tut {
//
// Test Group
@@ -36,6 +43,35 @@ struct test_distanceop_data {
test_distanceop_data()
: wktreader()
{}
+
+ static void checkDistance(const Geometry& g1, const Geometry& g2, double expected, double tol = 0)
+ {
+ ensure_equals("distance is incorrect", g1.distance(&g2), expected, tol);
+ ensure("unexpected isWithinDistance result for distance + tol", DistanceOp::isWithinDistance(g1, g2, expected + tol));
+ if (expected > tol) {
+ ensure("unexpected isWithinDistance result for distance - tol", !DistanceOp::isWithinDistance(g1, g2, expected - tol));
+ }
+
+ ensure_equals("reverse distance is incorrect", g2.distance(&g1), expected, tol);
+ }
+
+ static void checkLocation(const GeometryLocation& loc, const Geometry* component, std::size_t index, const geos::geom::CoordinateXY& c, double tol = 0)
+ {
+ ensure_equals("incorrect component", loc.getGeometryComponent(), component);
+ ensure_equals("incorrect position", loc.getSegmentIndex(), index);
+ ensure_equals_xy(loc.getCoordinate(), c, tol);
+ }
+
+ static void printLocs(const std::array<geos::operation::distance::GeometryLocation, 2>& locs)
+ {
+ if (locs[0].getCoordinate().equals2D(locs[1].getCoordinate())) {
+ std::cout << "POINT (" << locs[0].getCoordinate() << ")" << std::endl;
+
+ } else {
+ std::cout << "LINESTRING (" << locs[0].getCoordinate() << ", " << locs[1].getCoordinate() << ")" << std::endl;
+ }
+ std::cout << locs[0].getCoordinate().distance(locs[1].getCoordinate()) << std::endl;
+ }
};
typedef test_group<test_distanceop_data> group;
@@ -641,6 +677,332 @@ void object::test<27>()
ensure(std::isinf(g2->distance(g1.get())));
}
+template<>
+template<>
+void object::test<28>()
+{
+ set_test_name("CircularString / LineString (disjoint)");
+
+ auto g1 = wktreader.read("CIRCULARSTRING (-7 0, -6 -1, -5 0, 0 5, 4 -3)");
+ auto g2 = wktreader.read("LINESTRING (4 -5, 7 0, 13 0)");
+
+ checkDistance(*g1, *g2, 1.0024502, 1e-6);
+
+ auto locs = DistanceOp(*g1, *g2).nearestLocations();
+
+ checkLocation(locs[0], g1.get(), 2, XY(4.28746, -2.57248), 1e-5);
+ checkLocation(locs[1], g2.get(), 0, XY(5.14706, -3.08824), 1e-5);
+}
+
+template<>
+template<>
+void object::test<29>()
+{
+ set_test_name("CircularString / LineString (intersecting)");
+
+ auto g1 = wktreader.read("CIRCULARSTRING (-6 1, -6 -1, -5 0, 0 5, 4 -3)");
+ auto g2 = wktreader.read("LINESTRING (-4 2, 4 -2, 6 0)");
+
+ checkDistance(*g1, *g2, 0);
+
+ auto locs = DistanceOp(*g1, *g2).nearestLocations();
+
+ checkLocation(locs[0], g1.get(), 2, XY(4.87083, -1.12917), 1e-5);
+ checkLocation(locs[1], g2.get(), 1, XY(4.87083, -1.12917), 1e-5);
+}
+
+template<>
+template<>
+void object::test<30>()
+{
+ set_test_name("CircularString / MultiPoint (disjoint)");
+
+ auto g1 = wktreader.read("CIRCULARSTRING (-7 0, -6 -1, -5 0, 0 5, 4 -3)");
+ auto g2 = wktreader.read("MULTIPOINT (3.5 -5, 0 2, 6 4)");
+
+ checkDistance(*g1, *g2, 2.061553, 1e-6);
+
+ auto locs = DistanceOp(*g1, *g2).nearestLocations();
+
+ checkLocation(locs[0], g1.get(), 2, XY(4, -3), 1e-5);
+ checkLocation(locs[1], g2->getGeometryN(0), 0, XY(3.5, -5), 1e-5);
+}
+
+template<>
+template<>
+void object::test<31>()
+{
+ set_test_name("CircularString / MultiPoint (intersecting)");
+
+ auto g1 = wktreader.read("CIRCULARSTRING (-7 0, -6 -1, -5 0, 0 5, 4 -3)");
+ auto g2 = wktreader.read("MULTIPOINT (-6 15, -6 -1)");
+
+ checkDistance(*g1, *g2, 0);
+
+ auto locs = DistanceOp(*g1, *g2).nearestLocations();
+
+ checkLocation(locs[0], g1.get(), 0u, XY{-6, -1}, 1e-5);
+ checkLocation(locs[1], g2->getGeometryN(1), 0u, XY{-6, -1}, 1e-5);
+}
+
+template<>
+template<>
+void object::test<32>()
+{
+ set_test_name("CircularString / CircularString (disjoint)");
+
+ auto g1 = wktreader.read("CIRCULARSTRING (-5 0, 4 3, 4 -3, 3 0, 0 2, -2 0, -5 0)");
+ auto g2 = wktreader.read("CIRCULARSTRING (0 -2, -1 -1.5, -2 -2, 1 -1, 2 -2)");
+
+ checkDistance(*g1, *g2, 1.33846, 1e-5);
+
+ auto locs = DistanceOp(*g1, *g2).nearestLocations();
+
+ checkLocation(locs[0], g1.get(), 4, XY{-1.69468, 0.147266}, 1e-5);
+ checkLocation(locs[1], g2.get(), 2, XY{-1.06012, -1.03121}, 1e-5);
+}
+
+template<>
+template<>
+void object::test<33>()
+{
+ set_test_name("CircularString / CircularString (intersecting)");
+
+ auto g1 = wktreader.read("CIRCULARSTRING (-5 0, 4 3, 4 -3, 3 0, 0 2, -2 0, -5 0)");
+ auto g2 = wktreader.read("CIRCULARSTRING (0 -2, -1 -2.5, -2 -2, 1 -1, 4 -2)");
+
+ checkDistance(*g1, *g2, 0);
+
+ auto locs = DistanceOp(*g1, *g2).nearestLocations();
+
+ checkLocation(locs[0], g1.get(), 2, XY{3.86825, -1.9045}, 1e-5);
+ checkLocation(locs[1], g2.get(), 2, XY{3.86825, -1.9045}, 1e-5);
+}
+
+template<>
+template<>
+void object::test<34>()
+{
+ set_test_name("CircularString / CompoundCurve (disjoint)");
+
+ auto g1 = wktreader.read("CIRCULARSTRING (-5 0, 4 3, 4 -3, 3 0, 0 2, -2 0, -5 0)");
+ auto g2 = wktreader.read<geos::geom::CompoundCurve>("COMPOUNDCURVE(CIRCULARSTRING (0 -2, -1 -2.5, -2 -2), (-2 -2, -10 0, -5 5))");
+
+ checkDistance(*g1, *g2, 1.16409, 1e-5);
+
+ auto locs = DistanceOp(*g1, *g2).nearestLocations();
+
+ checkLocation(locs[0], g1.get(), 4, XY{-4.42355, -0.194193}, 1e-5);
+ checkLocation(locs[1], g2->getCurveN(1), 0, XY{-4.70588, -1.32353}, 1e-5);
+}
+
+template<>
+template<>
+void object::test<35>()
+{
+ set_test_name("MultiPoint partially in CurvePolygon");
+
+ auto g1 = wktreader.read("CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING(-5 0, 0 5, 5 0), (5 0, -5 0)))");
+ auto g2 = wktreader.read("MULTIPOINT (5 5, 1 1)");
+
+ checkDistance(*g1, *g2, 0);
+}
+
+template<>
+template<>
+void object::test<36>()
+{
+ set_test_name("CurvePolygon contained in CurvePolygon");
+
+ auto g1 = wktreader.read("CURVEPOLYGON (COMPOUNDCURVE (CIRCULARSTRING(-5 0, 0 5, 5 0), (5 0, 0 -5, -5 0)))");
+ auto g2 = wktreader.read("CURVEPOLYGON (CIRCULARSTRING (-1 0, 0 1, 1 0, 0 -1, -1 0))");
+
+ checkDistance(*g1, *g2, 0);
+}
+
+template<>
+template<>
+void object::test<40>()
+{
+ set_test_name("CurvePolygon and point, PostGIS issue #5989");
+
+ auto g1 = wktreader.read("CURVEPOLYGON(COMPOUNDCURVE((129296 142584,94722 100435,91618 97138,57306 60686,26874 28357,13059 34228,14572 65506,14593 65948,14616 66389),CIRCULARSTRING(14616 66389,17955 101124,24417 135415,24655 136418,24895 137421),(24895 137421,25472 139809,19354 141285,0 0,148000 142000,129296 142584)))");
+ auto g2 = wktreader.read("POINT(19925 112376)");
+
+ checkDistance(*g1, *g2, 199.655, 0.001);
+}
+
+template<>
+template<>
+void object::test<41>()
+{
+ set_test_name("CircularString and Point; PostGIS");
+
+ auto g1 = wktreader.read("CIRCULARSTRING(-1 0, 0 1, 1 0)");
+ auto g2 = wktreader.read("POINT (0 0)");
+
+ checkDistance(*g1, *g2, 1, 1e-6);
+}
+
+template<>
+template<>
+void object::test<42>()
+{
+ set_test_name("CircularString and Point (2); PostGIS");
+
+ auto g1 = wktreader.read("CIRCULARSTRING(-3 0, -2 0, -1 0, 0 1, 1 0)");
+ auto g2 = wktreader.read("POINT (0 0)");
+
+ checkDistance(*g1, *g2, 1, 1e-6);
+}
+
+template<>
+template<>
+void object::test<43>()
+{
+ set_test_name("CircularString and CircularString; PostGIS");
+
+ auto g1 = wktreader.read("CIRCULARSTRING(-1 0, 0 1, 1 0)");
+ auto g2 = wktreader.read( "CIRCULARSTRING(0 0, 1 -1, 2 0)");
+
+ checkDistance(*g1, *g2, 1, 1e-6);
+}
+
+template<>
+template<>
+void object::test<44>()
+{
+ set_test_name("CurvePolygon and Point; PostGIS");
+
+ auto g1 = wktreader.read( "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))");
+
+ auto p1 = wktreader.read("POINT(3 14)");
+ auto p2 = wktreader.read("POINT(3 8)");
+ auto p3 = wktreader.read("POINT(6 5)");
+ auto p4 = wktreader.read("POINT(6 4)");
+
+ checkDistance(*g1, *p1, 1, 1e-6);
+ checkDistance(*g1, *p2, 0, 1e-6);
+ checkDistance(*g1, *p3, 1, 1e-6);
+ checkDistance(*g1, *p4, 0, 1e-6);
+}
+
+template<>
+template<>
+void object::test<45>()
+{
+ set_test_name("CurvePolygon and LineString; PostGIS");
+
+ auto g1 = wktreader.read( "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))");
+
+ auto ls1 = wktreader.read("LINESTRING(0 0, 50 0)");
+ auto ls2 = wktreader.read( "LINESTRING(6 0, 10 7)");
+ auto ls3 = wktreader.read( "LINESTRING(4 4, 4 8)");
+ auto ls4 = wktreader.read( "LINESTRING(4 7, 5 6, 6 7)");
+ auto ls5 = wktreader.read( "LINESTRING(10 0, 10 2, 10 0)");
+
+ checkDistance(*g1, *ls1, 0.917484, 1e-6);
+ checkDistance(*g1, *ls2, 0, 1e-6);
+ checkDistance(*g1, *ls3, 0, 1e-6);
+ checkDistance(*g1, *ls4, 0.585786, 1e-6);
+ checkDistance(*g1, *ls5, 1.52913, 1e-6);
+}
+
+template<>
+template<>
+void object::test<46>()
+{
+ set_test_name("CurvePolygon and Polygon; PostGIS");
+
+ auto g1 = wktreader.read( "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))");
+
+ auto p1 = wktreader.read( "POLYGON((10 4, 10 8, 13 8, 13 4, 10 4))");
+ auto p2 = wktreader.read( "POLYGON((9 4, 9 8, 12 8, 12 4, 9 4))");
+ auto p3 = wktreader.read( "POLYGON((1 4, 1 8, 4 8, 4 4, 1 4))");
+
+ checkDistance(*g1, *p1, 0.58415, 1e-6);
+ checkDistance(*g1, *p2, 0);
+ checkDistance(*g1, *p3, 0);
+}
+
+template<>
+template<>
+void object::test<47>()
+{
+ set_test_name("CurvePolygon and CurvePolygon; PostGIS");
+
+ auto g1 = wktreader.read( "CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(1 6, 6 1, 9 7),(9 7, 3 13, 1 6)),COMPOUNDCURVE((3 6, 5 4, 7 4, 7 6),CIRCULARSTRING(7 6,5 8,3 6)))");
+
+ auto cp1 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(-1 4, 0 5, 1 4, 0 3, -1 4))");
+ auto cp2 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(1 4, 2 5, 3 4, 2 3, 1 4))");
+
+ checkDistance(*g1, *cp1, 0.0475666, 1e-6);
+ checkDistance(*g1, *cp2, 0);
+}
+
+template<>
+template<>
+void object::test<48>()
+{
+ set_test_name("MultiSurface and CurvePolygon; PostGIS");
+
+ auto g1 = wktreader.read( "MULTISURFACE(POLYGON((0 0,0 4,4 4,4 0,0 0)),CURVEPOLYGON(CIRCULARSTRING(8 2,10 4,12 2,10 0,8 2)))");
+
+ auto cp1 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(5 2,6 3,7 2,6 1,5 2))");
+ auto cp2 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(4 2,5 3,6 2,5 1,4 2))");
+ auto cp3 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(5 3,6 2,5 1,4 2,5 3))");
+ auto cp4 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(4.5 3,5.5 2,4.5 1,3.5 2,4.5 3))");
+ auto cp5 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(5.5 3,6.5 2,5.5 1,4.5 2,5.5 3))");
+ auto cp6 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(10 3,11 2,10 1,9 2,10 3))");
+ auto cp7 = wktreader.read( "CURVEPOLYGON(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3))");
+ auto cp8 = wktreader.read("CURVEPOLYGON(CIRCULARSTRING(5 7,6 8,7 7,6 6,5 7))");
+
+ checkDistance(*g1, *cp1, 1, 1e-6);
+ checkDistance(*g1, *cp2, 0);
+ checkDistance(*g1, *cp3, 0);
+ checkDistance(*g1, *cp4, 0);
+ checkDistance(*g1, *cp5, 0.5, 1e-6);
+ checkDistance(*g1, *cp6, 0);
+ checkDistance(*g1, *cp7, 0);
+ checkDistance(*g1, *cp8, 2.605551, 1e-6);
+}
+
+template<>
+template<>
+void object::test<49>()
+{
+ set_test_name("MultiCurve and LineString; PostGIS");
+
+ auto g1 = wktreader.read("LINESTRING(0.5 1,0.5 3)");
+ auto g2 = wktreader.read( "MULTICURVE(CIRCULARSTRING(2 3,3 2,2 1,1 2,2 3),(0 0, 0 5))");
+
+ checkDistance(*g1, *g2, 0.5, 1e-6);
+}
+
+template<>
+template<>
+void object::test<50>()
+{
+ set_test_name("CurvePolygon and Point; PostGIS ticket #4326");
+
+ auto g1 = wktreader.read("CURVEPOLYGON(CIRCULARSTRING(7874821 8715927,8907663 8715927,8844683 7750316,7937800 7750316,7874821 8715927))");
+ auto g2 = wktreader.read("POINT(5433865 8243495)");
+
+ checkDistance(*g1, *g2, 2271704.2698450615, 1e-6);
+}
+
+template<>
+template<>
+void object::test<51>()
+{
+ set_test_name("CurvePolygon and Point; PostGIS ticket #5989");
+
+ auto g1 = wktreader.read("CURVEPOLYGON(COMPOUNDCURVE(CIRCULARSTRING(0 0, -1 5, 0 10), (0 10, -10 10, -10 0, 0 0)))");
+ auto g2 = wktreader.read("POINT(-0.5 5)");
+
+ checkDistance(*g1, *g2, 0.5, 1e-6);
+}
+
// TODO: finish the tests by adding:
// LINESTRING - *all*
commit 8d42241e9c557ef9ca6bc97ecc71dd636ea82973
Author: Daniel Baston <dbaston at gmail.com>
Date: Mon Apr 20 17:08:47 2026 -0400
Add SurfaceExtracter
diff --git a/include/geos/geom/util/SurfaceExtracter.h b/include/geos/geom/util/SurfaceExtracter.h
new file mode 100644
index 000000000..1c9922099
--- /dev/null
+++ b/include/geos/geom/util/SurfaceExtracter.h
@@ -0,0 +1,60 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2026 ISciences LLC
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/export.h>
+#include <geos/geom/GeometryFilter.h>
+#include <vector>
+
+namespace geos::geom {
+class Surface;
+}
+
+namespace geos::geom::util {
+
+/**
+ * Extracts all the 2-dimensional (Surface) components from a Geometry.
+ */
+class GEOS_DLL SurfaceExtracter: public GeometryFilter {
+
+public:
+
+ /**
+ * Pushes the Surface components from a single geometry into
+ * the provided vector.
+ */
+ static void getSurfaces(const Geometry& geom, std::vector<const Surface*>& ret);
+
+ /**
+ * Constructs a SurfaceExtracter with a list in which
+ * to store Surfaces found.
+ */
+ SurfaceExtracter(std::vector<const Surface*>& newComps);
+
+ void filter_rw(Geometry* geom) override;
+
+ void filter_ro(const Geometry* geom) override;
+
+private:
+
+ std::vector<const Surface*>& comps;
+
+ // Declare type as noncopyable
+ SurfaceExtracter(const SurfaceExtracter& other) = delete;
+ SurfaceExtracter& operator=(const SurfaceExtracter& rhs) = delete;
+};
+
+} // namespace geos.geom.util
+
diff --git a/src/geom/util/SurfaceExtracter.cpp b/src/geom/util/SurfaceExtracter.cpp
new file mode 100644
index 000000000..9106d3ed5
--- /dev/null
+++ b/src/geom/util/SurfaceExtracter.cpp
@@ -0,0 +1,54 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2026 ISciences LLC
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/geom/util/SurfaceExtracter.h>
+
+#include <geos/geom/Surface.h>
+
+#include <vector>
+
+namespace geos::geom::util {
+
+void
+SurfaceExtracter::getSurfaces(const Geometry& geom, std::vector<const Surface*>& ret)
+{
+ if (!geom.hasDimension(Dimension::A)) {
+ return;
+ }
+
+ SurfaceExtracter pe(ret);
+ geom.apply_ro(&pe);
+}
+
+SurfaceExtracter::SurfaceExtracter(std::vector<const Surface*>& newComps)
+ :
+ comps(newComps)
+{}
+
+void
+SurfaceExtracter::filter_rw(Geometry* geom)
+{
+ if(const Surface* p = dynamic_cast<const Surface*>(geom)) {
+ comps.push_back(p);
+ }
+}
+
+void
+SurfaceExtracter::filter_ro(const Geometry* geom)
+{
+ if(const Surface* p = dynamic_cast<const Surface*>(geom)) {
+ comps.push_back(p);
+ }
+}
+}
commit e5beda552e1c949d234a4c1e45457ee0fd032e04
Author: Daniel Baston <dbaston at gmail.com>
Date: Mon Apr 20 15:07:20 2026 -0400
Add SimpleCurveExtracter
diff --git a/include/geos/geom/util/SimpleCurveExtracter.h b/include/geos/geom/util/SimpleCurveExtracter.h
new file mode 100644
index 000000000..a34c04407
--- /dev/null
+++ b/include/geos/geom/util/SimpleCurveExtracter.h
@@ -0,0 +1,52 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2026 ISciences LLC
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#pragma once
+
+#include <geos/export.h>
+
+#include <vector>
+#include <geos/geom/GeometryComponentFilter.h>
+
+namespace geos::geom {
+class SimpleCurve;
+}
+
+namespace geos::geom::util {
+
+/**
+ * Extracts all SimpleCurve components from a Geometry
+ */
+class GEOS_DLL SimpleCurveExtracter: public GeometryComponentFilter {
+
+public:
+
+ static void getCurves(const Geometry& geom, std::vector<const SimpleCurve*>& ret);
+
+ explicit SimpleCurveExtracter(std::vector<const SimpleCurve*>& newComps);
+
+ void filter_ro(const Geometry* geom) override;
+
+private:
+
+ std::vector<const SimpleCurve*>& comps;
+
+ // Declare type as noncopyable
+ SimpleCurveExtracter(const SimpleCurveExtracter& other) = delete;
+ SimpleCurveExtracter& operator=(const SimpleCurveExtracter& rhs) = delete;
+
+};
+
+} // namespace geos.geom.util
+
diff --git a/src/geom/util/SimpleCurveExtracter.cpp b/src/geom/util/SimpleCurveExtracter.cpp
new file mode 100644
index 000000000..1671b60c2
--- /dev/null
+++ b/src/geom/util/SimpleCurveExtracter.cpp
@@ -0,0 +1,66 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2026 ISciences LLC
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/geom/CompoundCurve.h>
+#include <geos/geom/GeometryComponentFilter.h>
+#include <geos/geom/SimpleCurve.h>
+#include <geos/geom/util/SimpleCurveExtracter.h>
+#include <geos/util.h>
+
+#include <vector>
+
+namespace geos {
+namespace geom { // geos.geom
+namespace util { // geos.geom.util
+
+SimpleCurveExtracter::SimpleCurveExtracter(std::vector<const SimpleCurve*>& newComps)
+ :
+ comps(newComps)
+{}
+
+void
+SimpleCurveExtracter::getCurves(const Geometry& geom, std::vector<const SimpleCurve*>& ret)
+{
+ if (geom.getDimension() == Dimension::P) {
+ return;
+ }
+
+ SimpleCurveExtracter lce(ret);
+ geom.apply_ro(&lce);
+}
+
+void
+SimpleCurveExtracter::filter_ro(const Geometry* geom)
+{
+ if (geom->isEmpty()) {
+ return;
+ }
+
+ const auto typ = geom->getGeometryTypeId();
+
+ if (typ == GEOS_LINEARRING || typ == GEOS_LINESTRING || typ == GEOS_CIRCULARSTRING) {
+ comps.push_back(detail::down_cast<const SimpleCurve*>(geom));
+ }
+
+ if (typ == GEOS_COMPOUNDCURVE) {
+ const CompoundCurve* cc = detail::down_cast<const CompoundCurve*>(geom);
+ for (std::size_t i = 0; i < cc->getNumCurves(); i++) {
+ comps.push_back(cc->getCurveN(i));
+ }
+ }
+}
+
+}
+}
+}
commit 2ef81880b9d0c9032902960f350e58591661f7bd
Author: Daniel Baston <dbaston at gmail.com>
Date: Mon Apr 20 11:17:14 2026 -0400
LineSegment: Add static impl of closestPoints
diff --git a/include/geos/geom/LineSegment.h b/include/geos/geom/LineSegment.h
index f6ef755ac..2f4b78427 100644
--- a/include/geos/geom/LineSegment.h
+++ b/include/geos/geom/LineSegment.h
@@ -417,6 +417,8 @@ public:
*/
double projectionFactor(const CoordinateXY& p) const;
+ static double projectionFactor(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q);
+
/** \brief
* Computes the fraction of distance (in <tt>[0.0, 1.0]</tt>)
* that the projection of a point occurs along this line segment.
@@ -446,6 +448,8 @@ public:
CoordinateXY project(const CoordinateXY& p) const;
+ static CoordinateXY project(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q);
+
/** \brief
* Project a line segment onto this line segment and return the resulting
* line segment.
@@ -471,6 +475,9 @@ public:
///
void closestPoint(const CoordinateXY& p, CoordinateXY& ret) const;
+ static CoordinateXY closestPoint(const CoordinateXY& p0, const CoordinateXY& p1,
+ const CoordinateXY& q);
+
/** \brief
* Returns <code>true</code> if <code>other</code> is
* topologically equal to this LineSegment (e.g. irrespective
@@ -496,6 +503,11 @@ public:
return closestPoints(*line);
}
+ static std::array<CoordinateXY, 2> closestPoints(const CoordinateXY& p0,
+ const CoordinateXY& p1,
+ const CoordinateXY& q0,
+ const CoordinateXY& q1);
+
/**
* Computes an intersection point between two segments,
* if there is one.
@@ -592,6 +604,8 @@ public:
private:
+ static CoordinateXY project(const CoordinateXY& p0, const CoordinateXY& p1, double factor);
+
void project(double factor, CoordinateXY& ret) const;
};
diff --git a/src/geom/LineSegment.cpp b/src/geom/LineSegment.cpp
index c352fc2cf..fd5db79da 100644
--- a/src/geom/LineSegment.cpp
+++ b/src/geom/LineSegment.cpp
@@ -45,14 +45,20 @@ LineSegment::reverse()
std::swap(p0, p1);
}
-/*public*/
double
LineSegment::projectionFactor(const CoordinateXY& p) const
{
- if(p == p0) {
+ return projectionFactor(p0, p1, p);
+}
+
+/*public*/
+double
+LineSegment::projectionFactor(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q)
+{
+ if(q == p0) {
return 0.0;
}
- if(p == p1) {
+ if(q == p1) {
return 1.0;
}
if(p0 == p1) {
@@ -72,7 +78,7 @@ LineSegment::projectionFactor(const CoordinateXY& p) const
double dx = p1.x - p0.x;
double dy = p1.y - p0.y;
double len2 = dx * dx + dy * dy;
- double r = ((p.x - p0.x) * dx + (p.y - p0.y) * dy) / len2;
+ double r = ((q.x - p0.x) * dx + (q.y - p0.y) * dy) / len2;
return r;
}
@@ -103,27 +109,38 @@ LineSegment::project(const Coordinate& p, Coordinate& ret) const
}
CoordinateXY
-LineSegment::project(const CoordinateXY& p) const
+LineSegment::project(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q)
{
- if(p == p0 || p == p1) {
- return p;
+ if(q == p0 || q == p1) {
+ return q;
}
- double r = projectionFactor(p);
+ double r = projectionFactor(p0, p1, q);
double x = p0.x + r * (p1.x - p0.x);
double y = p0.y + r * (p1.y - p0.y);
return CoordinateXY(x, y);
}
-
+CoordinateXY
+LineSegment::project(const CoordinateXY& p) const
+{
+ return project(p0, p1, p);
+}
/*private*/
+CoordinateXY
+LineSegment::project(const CoordinateXY& p0, const CoordinateXY& p1, double factor)
+{
+ if(factor == 1.0) {
+ return p1;
+ }
+
+ return CoordinateXY(p0.x + factor * (p1.x - p0.x), p0.y + factor * (p1.y - p0.y));
+}
+
void
LineSegment::project(double factor, CoordinateXY& ret) const
{
- if( factor == 1.0 )
- ret = p1;
- else
- ret = CoordinateXY(p0.x + factor * (p1.x - p0.x), p0.y + factor * (p1.y - p0.y));
+ ret = project(p0, p1, factor);
}
bool
@@ -156,6 +173,23 @@ LineSegment::project(const LineSegment& seg, LineSegment& ret) const
return true;
}
+CoordinateXY
+LineSegment::closestPoint(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& q)
+{
+ double factor = projectionFactor(p0, p1, q);
+ if(factor > 0 && factor < 1) {
+ return project(p0, p1, factor);
+ }
+
+ double dist0 = p0.distance(q);
+ double dist1 = p1.distance(q);
+ if(dist0 < dist1) {
+ return p0;
+ }
+
+ return p1;
+}
+
//Coordinate*
void
LineSegment::closestPoint(const CoordinateXY& p, CoordinateXY& ret) const
@@ -252,6 +286,56 @@ LineSegment::closestPoints(const LineSegment& line)
return closestPt;
}
+std::array<CoordinateXY, 2>
+LineSegment::closestPoints(const CoordinateXY& p0, const CoordinateXY& p1,
+ const CoordinateXY& q0, const CoordinateXY& q1)
+{
+ // test for intersection
+ algorithm::LineIntersector li;
+ li.computeIntersection(p0, p1, q0, q1);
+ if(li.hasIntersection()) {
+ return {CoordinateXY{li.getIntersection(0)},
+ CoordinateXY{li.getIntersection(0)}};
+ }
+
+ /*
+ * if no intersection closest pair contains at least one endpoint.
+ * Test each endpoint in turn.
+ */
+ std::array<CoordinateXY, 2> closestPt;
+
+ CoordinateXY close00 = closestPoint(p0, p1, q0);
+ double minDistance = close00.distance(q0);
+
+ closestPt[0] = close00;
+ closestPt[1] = q0;
+
+ CoordinateXY close01 = closestPoint(p0, p1, q1);
+ double dist = close01.distance(q1);
+ if(dist < minDistance) {
+ minDistance = dist;
+ closestPt[0] = close01;
+ closestPt[1] = q1;
+ }
+
+ CoordinateXY close10 = closestPoint(q0, q1, p0);
+ dist = close10.distance(p0);
+ if(dist < minDistance) {
+ minDistance = dist;
+ closestPt[0] = p0;
+ closestPt[1] = close10;
+ }
+
+ CoordinateXY close11 = closestPoint(q0, q1, p1);
+ dist = close11.distance(p1);
+ if(dist < minDistance) {
+ closestPt[0] = p1;
+ closestPt[1] = close11;
+ }
+
+ return closestPt;
+}
+
Coordinate
LineSegment::intersection(const LineSegment& line) const
{
diff --git a/src/operation/distance/DistanceOp.cpp b/src/operation/distance/DistanceOp.cpp
index fa25f5be2..9a234b4c2 100644
--- a/src/operation/distance/DistanceOp.cpp
+++ b/src/operation/distance/DistanceOp.cpp
@@ -489,12 +489,7 @@ DistanceOp::computeMinDistance(
if(dist < minDistance) {
minDistance = dist;
- // TODO avoid copy from constructing segs, maybe
- // by making a static closestPoints that takes four
- // coordinate references
- LineSegment seg0{Coordinate(p00), Coordinate(p01)};
- LineSegment seg1{Coordinate(p10), Coordinate(p11)};
- auto closestPt = seg0.closestPoints(seg1);
+ auto closestPt = LineSegment::closestPoints(p00, p01, p10, p11);
locGeom[0] = GeometryLocation(line0, i, closestPt[0]);
locGeom[1] = GeometryLocation(line1, j, closestPt[1]);
@@ -529,12 +524,10 @@ DistanceOp::computeMinDistance(const LineString* line,
if(dist < minDistance) {
minDistance = dist;
- // TODO avoid copy from constructing segs, maybe
- // by making a static closestPoints that takes three
- // coordinate references
- LineSegment seg{Coordinate(coord0->getAt<CoordinateXY>(i)), Coordinate(coord0->getAt<CoordinateXY>(i + 1))};
- Coordinate segClosestPoint;
- seg.closestPoint(*coord, segClosestPoint);
+ const auto& p0 = coord0->getAt<CoordinateXY>(i);
+ const auto& p1 = coord0->getAt<CoordinateXY>(i+1);
+
+ CoordinateXY segClosestPoint = LineSegment::closestPoint(p0, p1, *coord);
locGeom[0] = GeometryLocation(line, i, segClosestPoint);
locGeom[1] = GeometryLocation(pt, 0, *coord);
-----------------------------------------------------------------------
Summary of changes:
include/geos/algorithm/CircularArcIntersector.h | 13 -
include/geos/algorithm/CircularArcs.h | 194 ++++++-
include/geos/algorithm/PointLocator.h | 6 +-
include/geos/geom/CircularArc.h | 8 +
include/geos/geom/LineSegment.h | 14 +
include/geos/geom/util/SimpleCurveExtracter.h | 52 ++
include/geos/geom/util/SurfaceExtracter.h | 60 ++
include/geos/operation/distance/DistanceOp.h | 27 +-
src/algorithm/CircularArcIntersector.cpp | 67 +--
src/algorithm/CircularArcs.cpp | 423 +++++++++++++-
src/algorithm/PointLocator.cpp | 22 +-
src/geom/CircularArc.cpp | 82 ++-
src/geom/LineSegment.cpp | 110 +++-
src/geom/util/SimpleCurveExtracter.cpp | 66 +++
src/geom/util/SurfaceExtracter.cpp | 54 ++
.../distance/ConnectedElementLocationFilter.cpp | 28 +-
src/operation/distance/DistanceOp.cpp | 240 ++++++--
tests/unit/algorithm/CircularArcsTest.cpp | 628 ++++++++++++++++++++-
tests/unit/capi/GEOSDistanceTest.cpp | 11 +-
tests/unit/capi/GEOSDistanceWithinTest.cpp | 7 +-
tests/unit/geom/CircularStringTest.cpp | 6 +-
tests/unit/geom/CompoundCurveTest.cpp | 4 +-
tests/unit/geom/CurvePolygonTest.cpp | 4 +-
tests/unit/geom/MultiCurveTest.cpp | 4 +-
tests/unit/geom/MultiSurfaceTest.cpp | 4 +-
tests/unit/operation/distance/DistanceOpTest.cpp | 362 ++++++++++++
26 files changed, 2279 insertions(+), 217 deletions(-)
create mode 100644 include/geos/geom/util/SimpleCurveExtracter.h
create mode 100644 include/geos/geom/util/SurfaceExtracter.h
create mode 100644 src/geom/util/SimpleCurveExtracter.cpp
create mode 100644 src/geom/util/SurfaceExtracter.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list