[geos-commits] [SCM] GEOS branch main updated. b4aa3dff231341dddc67ef74e246ea32744d7756
git at osgeo.org
git at osgeo.org
Thu Aug 15 19:32:55 PDT 2024
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 b4aa3dff231341dddc67ef74e246ea32744d7756 (commit)
via c354044a6d15d02f58bd4739c313aab7c38e559c (commit)
via f093b649c0c87808c5de38f19d8447e3019a021a (commit)
via d639723a75618fb6c95b9b89ed860fd2bd4a9dd3 (commit)
via b7b28ddd50a7b21ce34397f925e7bbe3e605b2d9 (commit)
via 4f4c4d422ca9d2ddc1631e05c2562a536edaf7ca (commit)
via 8e2567e488e5bd1c200194b22e18f0d83cb885bf (commit)
from 6c9b94a7d85451c2701dc44bba7e16022b32b12c (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 b4aa3dff231341dddc67ef74e246ea32744d7756
Author: Daniel Baston <dbaston at gmail.com>
Date: Thu Aug 15 09:45:35 2024 -0400
CircularArc: Add comment from @rouault on optimizations
diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h
index 699ebd5f1..283eaf3b0 100644
--- a/include/geos/geom/CircularArc.h
+++ b/include/geos/geom/CircularArc.h
@@ -92,6 +92,11 @@ public:
return 2*MATH_PI;
}
+ /// Even Rouault:
+ /// potential optimization?: using crossproduct(p0 - center, p2 - center) = radius * radius * sin(angle)
+ /// could yield the result by just doing a single asin(), instead of 2 atan2()
+ /// actually one should also likely compute dotproduct(p0 - center, p2 - center) = radius * radius * cos(angle),
+ /// and thus angle = atan2(crossproduct(p0 - center, p2 - center) , dotproduct(p0 - center, p2 - center) )
auto t0 = theta0();
auto t2 = theta2();
commit c354044a6d15d02f58bd4739c313aab7c38e559c
Author: Daniel Baston <dbaston at gmail.com>
Date: Wed Jul 17 21:20:14 2024 -0400
Geometry::intersects: support CurvePolygon/Point arguments
diff --git a/src/geom/Geometry.cpp b/src/geom/Geometry.cpp
index 7d5f17c6f..433933587 100644
--- a/src/geom/Geometry.cpp
+++ b/src/geom/Geometry.cpp
@@ -39,6 +39,7 @@
#include <geos/algorithm/InteriorPointLine.h>
#include <geos/algorithm/InteriorPointArea.h>
#include <geos/algorithm/ConvexHull.h>
+#include <geos/algorithm/locate/SimplePointInAreaLocator.h>
#include <geos/geom/prep/PreparedGeometryFactory.h>
#include <geos/operation/intersection/Rectangle.h>
#include <geos/operation/intersection/RectangleIntersection.h>
@@ -322,10 +323,19 @@ Geometry::intersects(const Geometry* g) const
return predicate::RectangleIntersects::intersects(*p, *this);
}
+ auto typ = getGeometryTypeId();
+ if (typ == GEOS_CURVEPOLYGON && g->getGeometryTypeId() == GEOS_POINT) {
+ auto loc = locate::SimplePointInAreaLocator::locatePointInSurface(*g->getCoordinate(), *detail::down_cast<const Surface*>(this));
+ return loc != Location::EXTERIOR;
+ } else if (typ == GEOS_POINT && g->getGeometryTypeId() == GEOS_CURVEPOLYGON) {
+ auto loc = locate::SimplePointInAreaLocator::locatePointInSurface(*getCoordinate(), *detail::down_cast<const Surface*>(g));
+ return loc != Location::EXTERIOR;
+ }
+
#if USE_RELATENG
return operation::relateng::RelateNG::intersects(this, g);
#else
- if (getGeometryTypeId() == GEOS_GEOMETRYCOLLECTION) {
+ if (typ == GEOS_GEOMETRYCOLLECTION) {
auto im = relate(g);
bool res = im->isIntersects();
return res;
diff --git a/tests/unit/capi/GEOSIntersectsTest.cpp b/tests/unit/capi/GEOSIntersectsTest.cpp
index 395d2698c..b86d8c282 100644
--- a/tests/unit/capi/GEOSIntersectsTest.cpp
+++ b/tests/unit/capi/GEOSIntersectsTest.cpp
@@ -234,5 +234,30 @@ void object::test<11>()
ensure_equals("curved geometry not supported", GEOSIntersects(geom2_, geom1_), 2);
}
+// test PIP special case for CurvePolygon
+template<>
+template<>
+void object::test<12>()
+{
+ geom1_ = fromWKT("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))");
+ geom2_ = fromWKT("POINT (0.1556955 0.5355459)");
+
+ // PostGIS would return false here because geom2 is inside geom1
+ // but outside the linearized form of geom1
+ ensure_equals(GEOSIntersects(geom1_, geom2_), 1);
+ ensure_equals(GEOSIntersects(geom2_, geom1_), 1);
+}
+
+template<>
+template<>
+void object::test<13>()
+{
+ geom1_ = fromWKT("CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))");
+ geom2_ = fromWKT("POINT EMPTY");
+
+ ensure_equals(GEOSIntersects(geom1_, geom2_), 0);
+ ensure_equals(GEOSIntersects(geom2_, geom1_), 0);
+}
+
} // namespace tut
commit f093b649c0c87808c5de38f19d8447e3019a021a
Author: Daniel Baston <dbaston at gmail.com>
Date: Thu Jun 20 20:18:30 2024 -0400
SimplePointInAreaLocator: Support CurvePolygon
diff --git a/include/geos/algorithm/locate/IndexedPointInAreaLocator.h b/include/geos/algorithm/locate/IndexedPointInAreaLocator.h
index 44ddf88a1..ac327fe8d 100644
--- a/include/geos/algorithm/locate/IndexedPointInAreaLocator.h
+++ b/include/geos/algorithm/locate/IndexedPointInAreaLocator.h
@@ -60,7 +60,7 @@ private:
// p1 follows p0 in a CoordinateSequence, we know that the address
// of p1 is 16, 24, or 32 bytes greater than the address of p0.
// By packing this offset into the least significant bits of p0,
- // we can retrieve both p0 and p1 while only using 8 byytes.
+ // we can retrieve both p0 and p1 while only using 8 bytes.
std::size_t os = static_cast<std::size_t>(reinterpret_cast<const double*>(p1) - reinterpret_cast<const double*>(p0)) - 2u;
m_p0 = reinterpret_cast<std::size_t>(p0) | os;
diff --git a/include/geos/algorithm/locate/SimplePointInAreaLocator.h b/include/geos/algorithm/locate/SimplePointInAreaLocator.h
index 1149c9cdd..be77ff9d6 100644
--- a/include/geos/algorithm/locate/SimplePointInAreaLocator.h
+++ b/include/geos/algorithm/locate/SimplePointInAreaLocator.h
@@ -22,7 +22,7 @@ namespace geos {
namespace geom {
class Geometry;
class Coordinate;
-class Polygon;
+class Surface;
}
}
@@ -51,7 +51,7 @@ public:
const geom::Geometry* geom);
/** \brief
- * Determines the Location of a point in a [Polygon](@ref geom::Polygon).
+ * Determines the Location of a point in a [Surface](@ref geom::Surface).
*
* The return value is one of:
*
@@ -69,8 +69,8 @@ public:
* @param poly the geometry to test
* @return the Location of the point in the polygon
*/
- static geom::Location locatePointInPolygon(const geom::CoordinateXY& p,
- const geom::Polygon* poly);
+ static geom::Location locatePointInSurface(const geom::CoordinateXY& p,
+ const geom::Surface& poly);
/** \brief
* Determines whether a point is contained in a [Geometry](@ref geom::Geometry),
diff --git a/src/algorithm/locate/IndexedPointInAreaLocator.cpp b/src/algorithm/locate/IndexedPointInAreaLocator.cpp
index 0382c5deb..ce9ac78e4 100644
--- a/src/algorithm/locate/IndexedPointInAreaLocator.cpp
+++ b/src/algorithm/locate/IndexedPointInAreaLocator.cpp
@@ -30,7 +30,6 @@
#include <geos/index/ItemVisitor.h>
#include <algorithm>
-#include <typeinfo>
using geos::geom::CoordinateXY;
diff --git a/src/algorithm/locate/SimplePointInAreaLocator.cpp b/src/algorithm/locate/SimplePointInAreaLocator.cpp
index 5f1550a83..0a7b8e343 100644
--- a/src/algorithm/locate/SimplePointInAreaLocator.cpp
+++ b/src/algorithm/locate/SimplePointInAreaLocator.cpp
@@ -23,9 +23,6 @@
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/LineString.h>
-#include <typeinfo>
-#include <cassert>
-
using namespace geos::geom;
namespace geos {
@@ -63,9 +60,10 @@ SimplePointInAreaLocator::locateInGeometry(const CoordinateXY& p, const Geometry
}
if (geom->getNumGeometries() == 1) {
- if (geom->getGeometryTypeId() == GEOS_POLYGON) {
- auto poly = static_cast<const Polygon*>(geom);
- return locatePointInPolygon(p, poly);
+ auto typ = geom->getGeometryTypeId();
+ if (typ == GEOS_POLYGON || typ == GEOS_CURVEPOLYGON) {
+ auto surface = static_cast<const Surface*>(geom);
+ return locatePointInSurface(p, *surface);
}
}
for (std::size_t i = 0; i < geom->getNumGeometries(); i++) {
@@ -80,28 +78,25 @@ SimplePointInAreaLocator::locateInGeometry(const CoordinateXY& p, const Geometry
}
geom::Location
-SimplePointInAreaLocator::locatePointInPolygon(const CoordinateXY& p, const Polygon* poly)
+SimplePointInAreaLocator::locatePointInSurface(const CoordinateXY& p, const Surface& surface)
{
- if(poly->isEmpty()) {
+ if(surface.isEmpty()) {
return Location::EXTERIOR;
}
- if(!poly->getEnvelopeInternal()->contains(p)) {
+ if(!surface.getEnvelopeInternal()->contains(p)) {
return Location::EXTERIOR;
}
- const LineString* shell = poly->getExteriorRing();
- const CoordinateSequence* cl;
- cl = shell->getCoordinatesRO();
- Location shellLoc = PointLocation::locateInRing(p, *cl);
+ const Curve& shell = *surface.getExteriorRing();
+ Location shellLoc = PointLocation::locateInRing(p, shell);
if(shellLoc != Location::INTERIOR) {
return shellLoc;
}
// now test if the point lies in or on the holes
- for(std::size_t i = 0, n = poly->getNumInteriorRing(); i < n; i++) {
- const LineString* hole = poly->getInteriorRingN(i);
- if(hole->getEnvelopeInternal()->contains(p)) {
- cl = hole->getCoordinatesRO();
- Location holeLoc = RayCrossingCounter::locatePointInRing(p, *cl);
+ for(std::size_t i = 0; i < surface.getNumInteriorRing(); i++) {
+ const Curve& hole = *surface.getInteriorRingN(i);
+ if(hole.getEnvelopeInternal()->contains(p)) {
+ Location holeLoc = RayCrossingCounter::locatePointInRing(p, hole);
if(holeLoc == Location::BOUNDARY) {
return Location::BOUNDARY;
}
diff --git a/src/operation/predicate/RectangleIntersects.cpp b/src/operation/predicate/RectangleIntersects.cpp
index a2feea7e6..c31fe8ce9 100644
--- a/src/operation/predicate/RectangleIntersects.cpp
+++ b/src/operation/predicate/RectangleIntersects.cpp
@@ -166,7 +166,7 @@ protected:
// check rect point in poly (rect is known not to
// touch polygon at this point)
- if(SimplePointInAreaLocator::locatePointInPolygon(rectPt, poly) != geom::Location::EXTERIOR) {
+ if(SimplePointInAreaLocator::locatePointInSurface(rectPt, *poly) != geom::Location::EXTERIOR) {
containsPointVar = true;
return;
}
diff --git a/tests/unit/algorithm/LocatePointInRingTest.cpp b/tests/unit/algorithm/LocatePointInRingTest.cpp
index 605fbad6f..0dfacdff2 100644
--- a/tests/unit/algorithm/LocatePointInRingTest.cpp
+++ b/tests/unit/algorithm/LocatePointInRingTest.cpp
@@ -15,7 +15,6 @@
#include <geos/geom/LineString.h>
#include <geos/geom/Coordinate.h>
// std
-#include <sstream>
#include <string>
#include <memory>
diff --git a/tests/unit/algorithm/locate/SimplePointInAreaLocatorTest.cpp b/tests/unit/algorithm/locate/SimplePointInAreaLocatorTest.cpp
new file mode 100644
index 000000000..0330c1fed
--- /dev/null
+++ b/tests/unit/algorithm/locate/SimplePointInAreaLocatorTest.cpp
@@ -0,0 +1,43 @@
+#include <tut/tut.hpp>
+
+#include <geos/algorithm/locate/SimplePointInAreaLocator.h>
+#include <geos/geom/Location.h>
+#include <geos/io/WKTReader.h>
+
+using geos::geom::CoordinateXY;
+using geos::geom::Location;
+
+namespace tut {
+
+struct test_simplepointinarealocator_data {
+ geos::io::WKTReader reader;
+
+ static void checkLocation(const geos::geom::Geometry&g,
+ const geos::geom::CoordinateXY& pt,
+ geos::geom::Location loc) {
+ geos::algorithm::locate::SimplePointInAreaLocator locator(g);
+ ensure_equals(locator.locate(&pt), loc);
+ }
+
+};
+
+typedef test_group<test_simplepointinarealocator_data> group;
+typedef group::object object;
+
+group test_simplepointinarealocator_group("geos::algorithm::locate::SimplePointInAreaLocator");
+
+template<>
+template<>
+void object::test<1>()
+{
+ auto g = reader.read("MULTISURFACE( CURVEPOLYGON( CIRCULARSTRING( 0 0, 4 0, 4 4, 0 4, 0 0), (1 1, 3 3, 3 1, 1 1)), ((10 10, 14 12, 11 10, 10 10), (11 11, 11.5 11, 11 11.5, 11 11)))");
+
+ checkLocation(*g, {-2, 2}, Location::EXTERIOR); // outside envelope
+ checkLocation(*g, {0, -0.25}, Location::EXTERIOR); // inside envelope, outside shell
+ checkLocation(*g, {1, 1}, Location::BOUNDARY); // vertex of hole
+ checkLocation(*g, {2, 1}, Location::BOUNDARY); // boundary of hole
+ checkLocation(*g, {2.5, 1.5}, Location::EXTERIOR); // inside hole
+
+}
+
+}
commit d639723a75618fb6c95b9b89ed860fd2bd4a9dd3
Author: Daniel Baston <dbaston at gmail.com>
Date: Thu Jun 20 19:52:51 2024 -0400
CircularArc: Make method names more GEOSy
diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h
index 1a87583a3..699ebd5f1 100644
--- a/include/geos/geom/CircularArc.h
+++ b/include/geos/geom/CircularArc.h
@@ -31,20 +31,23 @@ public:
using CoordinateXY = geom::CoordinateXY;
- CircularArc(const CoordinateXY& q0, const CoordinateXY& q1, const CoordinateXY& q2) : p0(q0), p1(q1), p2(q2) {}
+ CircularArc(const CoordinateXY& q0, const CoordinateXY& q1, const CoordinateXY& q2)
+ : p0(q0)
+ , p1(q1)
+ , p2(q2)
+ , m_center_known(false)
+ , m_radius_known(false)
+ , m_orientation_known(false)
+ {}
const CoordinateXY& p0;
const CoordinateXY& p1;
const CoordinateXY& p2;
- bool isZeroLength() const {
- return center().equals2D(p0) || center().equals2D(p1);
- }
-
- bool isCollinear() const {
- return std::isnan(center().x);
- }
-
+ /// Return the orientation of the arc as one of:
+ /// - algorithm::Orientation::CLOCKWISE,
+ /// - algorithm::Orientation::COUNTERCLOCKWISE
+ /// - algorithm::Orientation::COLLINEAR
int orientation() const {
if (!m_orientation_known) {
m_orientation = algorithm::Orientation::index(p0, p1, p2);
@@ -53,7 +56,8 @@ public:
return m_orientation;
}
- const CoordinateXY& center() const {
+ /// Return the center point of the circle associated with this arc
+ const CoordinateXY& getCenter() const {
if (!m_center_known) {
m_center = algorithm::CircularArcs::getCenter(p0, p1, p2);
m_center_known = true;
@@ -62,24 +66,28 @@ public:
return m_center;
}
- double radius() const {
+ /// Return the radius of the circle associated with this arc
+ double getRadius() const {
if (!m_radius_known) {
- m_radius = center().distance(p0);
+ m_radius = getCenter().distance(p0);
m_radius_known = true;
}
return m_radius;
}
- bool isLinear() const {
- return std::isnan(radius());
- }
-
+ /// Return whether this arc forms a complete circle
bool isCircle() const {
return p0.equals(p2);
}
- double angle() const {
+ /// Returns whether this arc forms a straight line (p0, p1, and p2 are collinear)
+ bool isLinear() const {
+ return std::isnan(getRadius());
+ }
+
+ /// Return the inner angle of the sector associated with this arc
+ double getAngle() const {
if (isCircle()) {
return 2*MATH_PI;
}
@@ -100,37 +108,41 @@ public:
return diff;
}
- double length() const {
+ /// Return the length of the arc
+ double getLength() const {
if (isLinear()) {
return p0.distance(p2);
}
- return angle()*radius();
+ return getAngle()*getRadius();
}
+ /// Return the area enclosed by the arc p0-p1-p2 and the line segment p2-p0
double getArea() const {
if (isLinear()) {
return 0;
}
- auto R = radius();
- auto theta = angle();
+ auto R = getRadius();
+ auto theta = getAngle();
return R*R/2*(theta - std::sin(theta));
}
+ /// Return the angle of p0
double theta0() const {
- return std::atan2(p0.y - center().y, p0.x - center().x);
+ return std::atan2(p0.y - getCenter().y, p0.x - getCenter().x);
}
+ /// Return the angle of p2
double theta2() const {
- return std::atan2(p2.y - center().y, p2.x - center().x);
+ return std::atan2(p2.y - getCenter().y, p2.x - getCenter().x);
}
/// Check to see if a coordinate lies on the arc
/// Only the angle is checked, so it is assumed that the point lies on
/// the circle of which this arc is a part.
bool containsPointOnCircle(const CoordinateXY& q) const {
- double theta = std::atan2(q.y - center().y, q.x - center().x);
+ double theta = std::atan2(q.y - getCenter().y, q.x - getCenter().x);
return containsAngle(theta);
}
@@ -141,7 +153,7 @@ public:
return true;
}
- auto dist = std::abs(q.distance(center()) - radius());
+ auto dist = std::abs(q.distance(getCenter()) - getRadius());
if (dist > 1e-8) {
return false;
@@ -184,7 +196,7 @@ public:
/// at the location of a specified point. The point is assumed to
/// be on the arc.
bool isUpwardAtPoint(const CoordinateXY& q) const {
- auto quad = geom::Quadrant::quadrant(center(), q);
+ auto quad = geom::Quadrant::quadrant(getCenter(), q);
bool isUpward;
if (orientation() == algorithm::Orientation::CLOCKWISE) {
diff --git a/src/algorithm/CircularArcs.cpp b/src/algorithm/CircularArcs.cpp
index 0dd89633e..949f18509 100644
--- a/src/algorithm/CircularArcs.cpp
+++ b/src/algorithm/CircularArcs.cpp
@@ -13,12 +13,9 @@
**********************************************************************/
#include <geos/algorithm/CircularArcs.h>
-
-#include "geos/algorithm/Angle.h"
-#include "geos/geom/Envelope.h"
-#include "geos/geom/Quadrant.h"
-
-#include <array>
+#include <geos/algorithm/Orientation.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Quadrant.h>
using geos::geom::CoordinateXY;
diff --git a/src/algorithm/RayCrossingCounter.cpp b/src/algorithm/RayCrossingCounter.cpp
index e7b442b21..9e25deba8 100644
--- a/src/algorithm/RayCrossingCounter.cpp
+++ b/src/algorithm/RayCrossingCounter.cpp
@@ -211,8 +211,8 @@ RayCrossingCounter::shouldCountCrossing(const geom::CircularArc& arc, const geom
/// some Coordinates in the returned array will be equal to CoordinateXY::getNull().
std::array<geom::CoordinateXY, 2>
RayCrossingCounter::pointsIntersectingHorizontalRay(const geom::CircularArc& arc, const geom::CoordinateXY& origin) {
- auto c = arc.center();
- auto R = arc.radius();
+ const auto& c = arc.getCenter();
+ const auto& R = arc.getRadius();
auto dx = std::sqrt(R*R - std::pow(origin.y - c.y, 2) );
@@ -262,7 +262,7 @@ RayCrossingCounter::countArc(const CoordinateXY& p1,
geom::CircularArc arc(p1, p2, p3);
// If the arc is degenerate, process it is two line segments
- if (arc.isCollinear()) {
+ if (arc.isLinear()) {
countSegment(p1, p2);
countSegment(p2, p3);
return;
diff --git a/src/geom/CircularString.cpp b/src/geom/CircularString.cpp
index 7a758f7a1..b666c8eb1 100644
--- a/src/geom/CircularString.cpp
+++ b/src/geom/CircularString.cpp
@@ -61,7 +61,7 @@ CircularString::getLength() const
double tot = 0;
for (std::size_t i = 2; i < coords.size(); i += 2) {
- auto len = CircularArc(coords[i-2], coords[i-1], coords[i]).length();
+ auto len = CircularArc(coords[i-2], coords[i-1], coords[i]).getLength();
tot += len;
}
return tot;
diff --git a/tests/unit/geom/CircularArcTest.cpp b/tests/unit/geom/CircularArcTest.cpp
index a5d3c5772..a1f730091 100644
--- a/tests/unit/geom/CircularArcTest.cpp
+++ b/tests/unit/geom/CircularArcTest.cpp
@@ -16,18 +16,18 @@ struct test_circulararc_data {
void checkAngle(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2, double expected) {
CircularArc arc(p0, p1, p2);
- ensure_equals(p0.toString() + " / " + p1.toString() + " / " + p2.toString(), arc.angle(), expected, eps);
+ ensure_equals(p0.toString() + " / " + p1.toString() + " / " + p2.toString(), arc.getAngle(), expected, eps);
CircularArc rev(p2, p1, p0);
- ensure_equals(p2.toString() + " / " + p1.toString() + " / " + p0.toString(), rev.angle(), expected, eps);
+ ensure_equals(p2.toString() + " / " + p1.toString() + " / " + p0.toString(), rev.getAngle(), expected, eps);
}
void checkLength(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2, double expected) {
CircularArc arc(p0, p1, p2);
- ensure_equals(p0.toString() + " / " + p1.toString() + " / " + p2.toString(), arc.length(), expected, eps);
+ ensure_equals(p0.toString() + " / " + p1.toString() + " / " + p2.toString(), arc.getLength(), expected, eps);
CircularArc rev(p2, p1, p0);
- ensure_equals(p2.toString() + " / " + p1.toString() + " / " + p0.toString(), rev.length(), expected, eps);
+ ensure_equals(p2.toString() + " / " + p1.toString() + " / " + p0.toString(), rev.getLength(), expected, eps);
}
};
commit b7b28ddd50a7b21ce34397f925e7bbe3e605b2d9
Author: Daniel Baston <dbaston at gmail.com>
Date: Thu Jun 20 17:31:24 2024 -0400
Curved geometries: Implement getArea
diff --git a/include/geos/algorithm/Area.h b/include/geos/algorithm/Area.h
index f53938c1a..c3c799f87 100644
--- a/include/geos/algorithm/Area.h
+++ b/include/geos/algorithm/Area.h
@@ -23,12 +23,19 @@
#include <geos/geom/CoordinateSequence.h>
namespace geos {
+
+namespace geom {
+class Curve;
+}
+
namespace algorithm { // geos::algorithm
class GEOS_DLL Area {
public:
+ static double ofClosedCurve(const geom::Curve& ring);
+
/**
* Computes the area for a ring.
*
diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h
index 69ca66746..1a87583a3 100644
--- a/include/geos/geom/CircularArc.h
+++ b/include/geos/geom/CircularArc.h
@@ -108,6 +108,16 @@ public:
return angle()*radius();
}
+ double getArea() const {
+ if (isLinear()) {
+ return 0;
+ }
+
+ auto R = radius();
+ auto theta = angle();
+ return R*R/2*(theta - std::sin(theta));
+ }
+
double theta0() const {
return std::atan2(p0.y - center().y, p0.x - center().x);
}
diff --git a/include/geos/geom/CircularString.h b/include/geos/geom/CircularString.h
index fd43bc2b9..85dc3014e 100644
--- a/include/geos/geom/CircularString.h
+++ b/include/geos/geom/CircularString.h
@@ -41,6 +41,10 @@ public:
return true;
}
+ bool isCurved() const override {
+ return true;
+ }
+
std::unique_ptr<CircularString> reverse() const
{
return std::unique_ptr<CircularString>(reverseImpl());
diff --git a/include/geos/geom/LineString.h b/include/geos/geom/LineString.h
index faf021f20..478970874 100644
--- a/include/geos/geom/LineString.h
+++ b/include/geos/geom/LineString.h
@@ -92,6 +92,10 @@ public:
double getLength() const override;
+ bool isCurved() const override {
+ return false;
+ }
+
/**
* Creates a LineString whose coordinates are in the reverse
* order of this object's
diff --git a/include/geos/geom/SimpleCurve.h b/include/geos/geom/SimpleCurve.h
index a41ea0ee5..1bf674f69 100644
--- a/include/geos/geom/SimpleCurve.h
+++ b/include/geos/geom/SimpleCurve.h
@@ -94,6 +94,8 @@ public:
virtual bool isCoordinate(CoordinateXY& pt) const;
+ virtual bool isCurved() const = 0;
+
bool isEmpty() const override;
/** \brief
diff --git a/src/algorithm/Area.cpp b/src/algorithm/Area.cpp
index e51954dbc..3d2375c95 100644
--- a/src/algorithm/Area.cpp
+++ b/src/algorithm/Area.cpp
@@ -20,6 +20,10 @@
#include <vector>
#include <geos/algorithm/Area.h>
+#include <geos/geom/CircularArc.h>
+#include <geos/geom/Curve.h>
+#include <geos/geom/SimpleCurve.h>
+#include <geos/util/IllegalArgumentException.h>
using geos::geom::CoordinateXY;
@@ -93,7 +97,58 @@ Area::ofRingSigned(const geom::CoordinateSequence* ring)
return sum / 2.0;
}
+double
+Area::ofClosedCurve(const geom::Curve& ring) {
+ if (!ring.isClosed()) {
+ throw util::IllegalArgumentException("Argument is not closed");
+ }
+ double sum = 0;
+
+ for (std::size_t i = 0; i < ring.getNumCurves(); i++) {
+ const geom::SimpleCurve& section = *ring.getCurveN(i);
+
+ if (section.isEmpty()) {
+ continue;
+ }
+
+ const geom::CoordinateSequence& coords = *section.getCoordinatesRO();
+
+ if (section.isCurved()) {
+ for (std::size_t j = 2; j < coords.size(); j += 2) {
+ const CoordinateXY& p0 = coords.getAt<CoordinateXY>(j-2);
+ const CoordinateXY& p1 = coords.getAt<CoordinateXY>(j-1);
+ const CoordinateXY& p2 = coords.getAt<CoordinateXY>(j);
+
+ double triangleArea = 0.5*(p0.x*p2.y - p2.x*p0.y);
+ sum += triangleArea;
+
+ geom::CircularArc arc(p0, p1, p2);
+ if (arc.isLinear()) {
+ continue;
+ }
+
+ double circularSegmentArea = arc.getArea();
+
+ if (algorithm::Orientation::index(p0, p2, p1) == algorithm::Orientation::CLOCKWISE) {
+ sum += circularSegmentArea;
+ } else {
+ sum -= circularSegmentArea;
+ }
+ }
+ } else {
+ for (std::size_t j = 1; j < coords.size(); j++) {
+ const CoordinateXY& p0 = coords.getAt<CoordinateXY>(j-1);
+ const CoordinateXY& p1 = coords.getAt<CoordinateXY>(j);
+
+ double triangleArea = 0.5*(p0.x*p1.y - p1.x*p0.y);
+ sum += triangleArea;
+ }
+ }
+ }
+
+ return std::abs(sum);
+}
} // namespace geos.algorithm
} //namespace geos
diff --git a/src/geom/CurvePolygon.cpp b/src/geom/CurvePolygon.cpp
index 74a5e943b..c5213602c 100644
--- a/src/geom/CurvePolygon.cpp
+++ b/src/geom/CurvePolygon.cpp
@@ -12,6 +12,7 @@
*
**********************************************************************/
+#include <geos/algorithm/Area.h>
#include <geos/geom/Curve.h>
#include <geos/geom/CurvePolygon.h>
#include <geos/geom/CoordinateSequence.h>
@@ -53,7 +54,11 @@ namespace geom {
}
double CurvePolygon::getArea() const {
- throw util::UnsupportedOperationException();
+ double sum = algorithm::Area::ofClosedCurve(*shell);
+ for (const auto& hole : holes) {
+ sum -= algorithm::Area::ofClosedCurve(*hole);
+ }
+ return sum;
}
bool CurvePolygon::hasCurvedComponents() const {
diff --git a/tests/unit/algorithm/AreaTest.cpp b/tests/unit/algorithm/AreaTest.cpp
index 1a81ff08b..f40de5733 100644
--- a/tests/unit/algorithm/AreaTest.cpp
+++ b/tests/unit/algorithm/AreaTest.cpp
@@ -13,12 +13,12 @@
#include <geos/geom/PrecisionModel.h>
#include <geos/io/WKTReader.h>
// std
-#include <sstream>
#include <string>
#include <memory>
using namespace geos;
using namespace geos::geom;
+using geos::algorithm::Area;
namespace tut {
//
@@ -33,7 +33,7 @@ struct test_area_data {
geos::io::WKTReader reader_;
test_area_data():
geom_(nullptr),
- pm_(1),
+ pm_(),
factory_(GeometryFactory::create(&pm_, 0)), reader_(factory_.get())
{
assert(nullptr == geom_);
@@ -48,19 +48,23 @@ struct test_area_data {
void
checkAreaOfRing(std::string wkt, double expectedArea)
{
- std::unique_ptr<Geometry> lineGeom(reader_.read(wkt));
- std::unique_ptr<LineString> line(dynamic_cast<LineString*>(lineGeom.release()));
- ensure(nullptr != line.get());
- const CoordinateSequence* ringSeq = line->getCoordinatesRO();
+ auto ringGeom = reader_.read<Curve>(wkt);
- std::vector<Coordinate> ringCoords;
- ringSeq->toVector(ringCoords);
+ if (const LineString* line = dynamic_cast<const LineString*>(ringGeom.get())) {
+ const CoordinateSequence* ringSeq = line->getCoordinatesRO();
- double actual1 = algorithm::Area::ofRing(ringCoords);
- double actual2 = algorithm::Area::ofRing(ringSeq);
+ std::vector<Coordinate> ringCoords;
+ ringSeq->toVector(ringCoords);
- ensure_equals(actual1, expectedArea);
- ensure_equals(actual2, expectedArea);
+ double actual1 = algorithm::Area::ofRing(ringCoords);
+ double actual2 = algorithm::Area::ofRing(ringSeq);
+
+ ensure_equals(actual1, expectedArea);
+ ensure_equals(actual2, expectedArea);
+ }
+
+ double actual3 = algorithm::Area::ofClosedCurve(*ringGeom);
+ ensure_equals("Area::ofClosedCurve", actual3, expectedArea, 1e-6);
}
void
@@ -115,6 +119,41 @@ void object::test<3>
checkAreaOfRingSigned("LINESTRING (100 200, 100 100, 200 100, 200 200, 100 200)", -10000.0);
}
+template<>
+template<>
+void object::test<4>
+()
+{
+ checkAreaOfRing("CIRCULARSTRING (0 0, 2 2, 4 0, 2 -2, 0 0)", 4*MATH_PI);
+}
+
+template<>
+template<>
+void object::test<5>
+()
+{
+ checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 2, 4 0), (4 0, 0 0))", 2*MATH_PI);
+}
+
+template<>
+template<>
+void object::test<6>
+()
+{
+ // expected area from PostGIS after ST_CurveToLine(geom, 1e-13, 1)
+ checkAreaOfRing("CIRCULARSTRING (0 0, 2 2, 4 0, 2 1, 0 0)", 3.48759);
+}
+
+template<>
+template<>
+void object::test<7>
+()
+{
+ // expected area from PostGIS after ST_CurveToLine(geom, 1e-13, 1)
+ checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 0, 2 1, 2 3, 4 3, 4 5, 1 4, 0.5 0.8, 0 0))", 11.243342);
+ checkAreaOfRing("COMPOUNDCURVE (CIRCULARSTRING (0 0, 2 0, 2 1, 2 3, 4 3), (4 3, 4 5, 1 4, 0 0))", 9.321903);
+}
+
} // namespace tut
diff --git a/tests/unit/geom/CircularArcTest.cpp b/tests/unit/geom/CircularArcTest.cpp
index 00f7bbcbd..a5d3c5772 100644
--- a/tests/unit/geom/CircularArcTest.cpp
+++ b/tests/unit/geom/CircularArcTest.cpp
@@ -73,4 +73,19 @@ void object::test<2>()
checkLength({1.6, 0.4}, {1.6, 0.5}, {1.7, 1}, 0.6122445326877711);
}
+
+// test getArea()
+template<>
+template<>
+void object::test<3>()
+{
+ ensure_equals("half circle, R=2", CircularArc({-2, 0}, {0, 2}, {2, 0}).getArea(), MATH_PI*2);
+
+ ensure_equals("full circle, R=3", CircularArc({-3, 0}, {3, 0}, {-3, 0}).getArea(), MATH_PI*3*3);
+
+ ensure_equals("3/4, mouth up, R=2", CircularArc({-std::sqrt(2), std::sqrt(2)}, {0, -2}, {std::sqrt(2), std::sqrt(2)}).getArea(), MATH_PI*4 - 2*(MATH_PI/2-1), 1e-8);
+
+ ensure_equals("1/4, pointing up, R=2", CircularArc({-std::sqrt(2), std::sqrt(2)}, {0, 2}, {std::sqrt(2), std::sqrt(2)}).getArea(), 2*(MATH_PI/2-1), 1e-8);
+}
+
}
diff --git a/tests/unit/geom/CurvePolygonTest.cpp b/tests/unit/geom/CurvePolygonTest.cpp
index 9cafc6516..0ee7d4baa 100644
--- a/tests/unit/geom/CurvePolygonTest.cpp
+++ b/tests/unit/geom/CurvePolygonTest.cpp
@@ -74,7 +74,7 @@ void object::test<1>()
ensure("getCoordinates", cp->getCoordinates()->isEmpty());
ensure("getCoordinate", cp->getCoordinate() == nullptr);
- ensure_THROW(cp->getArea(), geos::util::UnsupportedOperationException);
+ ensure_equals("getArea", cp->getArea(), 0.0);
ensure_equals("getLength", cp->getLength(), 0.0);
}
@@ -90,7 +90,7 @@ void object::test<2>()
// Geometry size functions
ensure("isEmpty", !cp_->isEmpty());
- ensure_THROW(cp_->getArea(), geos::util::UnsupportedOperationException);
+ ensure_equals("getArea", cp_->getArea(), 9.0526564962674, 1e-8); // expected value from PostGIS with ST_CurveToLine(geom, 1e-13, 1)
ensure_equals("getLength", cp_->getLength(), 19.236489581872586, 1e-8);
ensure_equals("getNumGeometries", cp_->getNumGeometries(), 1u);
ensure_equals("getNumPoints", cp_->getNumPoints(), 14u);
diff --git a/tests/unit/geom/MultiSurfaceTest.cpp b/tests/unit/geom/MultiSurfaceTest.cpp
index 751010c23..1492f2279 100644
--- a/tests/unit/geom/MultiSurfaceTest.cpp
+++ b/tests/unit/geom/MultiSurfaceTest.cpp
@@ -85,7 +85,7 @@ void object::test<2>()
// Geometry size functions
ensure("isEmpty", !ms_->isEmpty());
- ensure_THROW(ms_->getArea(), geos::util::UnsupportedOperationException);
+ ensure_equals("getArea", ms_->getArea(), 4.141592653589132, 1e-6); // expected value from PostGIS with ST_CurveToLine(geom, 1e-13, 1)
ensure_equals("getLength", ms_->getLength(), 10.283185307179586);
ensure_equals("getNumGeometries", ms_->getNumGeometries(), 2u);
ensure_equals("getNumPoints", ms_->getNumPoints(), 10u);
commit 4f4c4d422ca9d2ddc1631e05c2562a536edaf7ca
Author: Daniel Baston <dbaston at gmail.com>
Date: Wed Jun 19 13:17:07 2024 -0400
Curved geometries: Implement getLength()
diff --git a/include/geos/algorithm/CircularArcs.h b/include/geos/algorithm/CircularArcs.h
index 7e3b9b7bb..54f0a9b7d 100644
--- a/include/geos/algorithm/CircularArcs.h
+++ b/include/geos/algorithm/CircularArcs.h
@@ -23,7 +23,6 @@ namespace algorithm {
class GEOS_DLL CircularArcs {
public:
- static constexpr double PI = 3.14159265358979323846;
/// Return the circle center of an arc defined by three points
static geom::CoordinateXY getCenter(const geom::CoordinateXY& p0, const geom::CoordinateXY& p1,
diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h
index 758733738..69ca66746 100644
--- a/include/geos/geom/CircularArc.h
+++ b/include/geos/geom/CircularArc.h
@@ -47,7 +47,7 @@ public:
int orientation() const {
if (!m_orientation_known) {
- m_orientation = algorithm::Orientation::index(center(), p0, p1);
+ m_orientation = algorithm::Orientation::index(p0, p1, p2);
m_orientation_known = true;
}
return m_orientation;
@@ -71,6 +71,43 @@ public:
return m_radius;
}
+ bool isLinear() const {
+ return std::isnan(radius());
+ }
+
+ bool isCircle() const {
+ return p0.equals(p2);
+ }
+
+ double angle() const {
+ if (isCircle()) {
+ return 2*MATH_PI;
+ }
+
+ auto t0 = theta0();
+ auto t2 = theta2();
+
+ if (orientation() == algorithm::Orientation::COUNTERCLOCKWISE) {
+ std::swap(t0, t2);
+ }
+
+ if (t0 < t2) {
+ t0 += 2*MATH_PI;
+ }
+
+ auto diff = t0-t2;
+
+ return diff;
+ }
+
+ double length() const {
+ if (isLinear()) {
+ return p0.distance(p2);
+ }
+
+ return angle()*radius();
+ }
+
double theta0() const {
return std::atan2(p0.y - center().y, p0.x - center().x);
}
@@ -123,11 +160,11 @@ public:
t2 -= t0;
theta -= t0;
- if (t2 < 0){
- t2 += 2*algorithm::CircularArcs::PI;
+ if (t2 < 0) {
+ t2 += 2*MATH_PI;
}
if (theta < 0) {
- theta += 2*algorithm::CircularArcs::PI;
+ theta += 2*MATH_PI;
}
return theta >= t2;
diff --git a/include/geos/geom/CircularString.h b/include/geos/geom/CircularString.h
index 283cb636c..fd43bc2b9 100644
--- a/include/geos/geom/CircularString.h
+++ b/include/geos/geom/CircularString.h
@@ -15,7 +15,6 @@
#pragma once
#include <geos/geom/SimpleCurve.h>
-#include <geos/util/UnsupportedOperationException.h>
namespace geos {
namespace geom {
@@ -35,10 +34,7 @@ public:
GeometryTypeId getGeometryTypeId() const override;
- double getLength() const override
- {
- throw util::UnsupportedOperationException("Cannot calculate length of CircularString");
- }
+ double getLength() const override;
bool hasCurvedComponents() const override
{
diff --git a/src/geom/CircularString.cpp b/src/geom/CircularString.cpp
index f5bb375c8..7a758f7a1 100644
--- a/src/geom/CircularString.cpp
+++ b/src/geom/CircularString.cpp
@@ -12,6 +12,7 @@
*
**********************************************************************/
+#include <geos/geom/CircularArc.h>
#include <geos/geom/CircularString.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/GeometryFactory.h>
@@ -49,6 +50,23 @@ CircularString::getGeometryTypeId() const
return GEOS_CIRCULARSTRING;
}
+double
+CircularString::getLength() const
+{
+ if (isEmpty()) {
+ return 0;
+ }
+
+ const CoordinateSequence& coords = *getCoordinatesRO();
+
+ double tot = 0;
+ for (std::size_t i = 2; i < coords.size(); i += 2) {
+ auto len = CircularArc(coords[i-2], coords[i-1], coords[i]).length();
+ tot += len;
+ }
+ return tot;
+}
+
CircularString*
CircularString::reverseImpl() const
{
diff --git a/tests/unit/capi/GEOSLengthTest.cpp b/tests/unit/capi/GEOSLengthTest.cpp
index 66c2872d1..7be71de71 100644
--- a/tests/unit/capi/GEOSLengthTest.cpp
+++ b/tests/unit/capi/GEOSLengthTest.cpp
@@ -1,6 +1,7 @@
#include <tut/tut.hpp>
// geos
#include <geos_c.h>
+#include <geos/constants.h>
#include "capi_test_utils.h"
@@ -67,7 +68,8 @@ void object::test<4>()
double length = -1;
int ret = GEOSLength(input_, &length);
- ensure_equals("error raised on curved geometry", ret, 0);
+ ensure_equals(ret, 1);
+ ensure_equals(length, geos::MATH_PI);
}
} // namespace tut
diff --git a/tests/unit/geom/CircularArcTest.cpp b/tests/unit/geom/CircularArcTest.cpp
new file mode 100644
index 000000000..00f7bbcbd
--- /dev/null
+++ b/tests/unit/geom/CircularArcTest.cpp
@@ -0,0 +1,76 @@
+#include <tut/tut.hpp>
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CircularArc.h>
+#include <geos/constants.h>
+
+using geos::geom::CoordinateXY;
+using geos::geom::CircularArc;
+using geos::MATH_PI;
+
+namespace tut {
+
+struct test_circulararc_data {
+
+ const double eps = 1e-8;
+
+ void checkAngle(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2, double expected) {
+ CircularArc arc(p0, p1, p2);
+ ensure_equals(p0.toString() + " / " + p1.toString() + " / " + p2.toString(), arc.angle(), expected, eps);
+
+ CircularArc rev(p2, p1, p0);
+ ensure_equals(p2.toString() + " / " + p1.toString() + " / " + p0.toString(), rev.angle(), expected, eps);
+ }
+
+ void checkLength(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2, double expected) {
+ CircularArc arc(p0, p1, p2);
+ ensure_equals(p0.toString() + " / " + p1.toString() + " / " + p2.toString(), arc.length(), expected, eps);
+
+ CircularArc rev(p2, p1, p0);
+ ensure_equals(p2.toString() + " / " + p1.toString() + " / " + p0.toString(), rev.length(), expected, eps);
+ }
+};
+
+using group = test_group<test_circulararc_data>;
+using object = group::object;
+
+group test_circulararc_group("geos::geom::CircularArc");
+
+// test angle() on unit circle
+template<>
+template<>
+void object::test<1>()
+{
+ auto x = std::sqrt(2.0)/2;
+
+ // full circle
+ checkAngle({-1, 0}, {1, 0}, {-1, 0}, 2*MATH_PI);
+
+ // check half-circles
+ checkAngle({-1, 0}, {0, 1}, {1, 0}, MATH_PI); // top
+ checkAngle({-1, 0}, {0, -1}, {1, 0}, MATH_PI); // bottom
+ checkAngle({0, -1}, {-1, 0}, {0, 1}, MATH_PI); // left
+ checkAngle({0, -1}, {1, 0}, {0, 1}, MATH_PI); // right
+
+ // check quadrants
+ checkAngle({-1, 0}, {-x, x}, {0, 1}, MATH_PI/2); // upper left
+ checkAngle({0, 1}, {x, x}, {1, 0}, MATH_PI/2); // upper right
+ checkAngle({0, -1}, {x, -x}, {1, 0}, MATH_PI/2); // lower right
+ checkAngle({0, -1}, {-x, -x}, {-1, 0}, MATH_PI/2); // lower left
+
+ // 3/4
+ checkAngle({-x, x}, {0, -1}, {x, x}, 1.5*MATH_PI); // mouth up
+ checkAngle({-x, -x}, {0, 1}, {x, -x}, 1.5*MATH_PI); // mouth down
+ checkAngle({-x, x}, {1, 0}, {-x, -x}, 1.5*MATH_PI); // mouth left
+ checkAngle({x, -x}, {-1, 0}, {x, x}, 1.5*MATH_PI); // mouth right
+}
+
+// test length()
+template<>
+template<>
+void object::test<2>()
+{
+ checkLength({1.6, 0.4}, {1.6, 0.5}, {1.7, 1}, 0.6122445326877711);
+}
+
+}
diff --git a/tests/unit/geom/CircularStringTest.cpp b/tests/unit/geom/CircularStringTest.cpp
index 8e1eaf4f3..64bbe6524 100644
--- a/tests/unit/geom/CircularStringTest.cpp
+++ b/tests/unit/geom/CircularStringTest.cpp
@@ -58,7 +58,7 @@ void object::test<1>()
ensure(cs->getCoordinate() == nullptr);
ensure_equals(cs->getArea(), 0);
- ensure_THROW(cs_->getLength(), geos::util::UnsupportedOperationException);
+ ensure_equals(cs->getLength(), 0);
}
// Basic Geometry API
@@ -74,7 +74,7 @@ void object::test<2>()
// Geometry size functions
ensure("isEmpty", !cs_->isEmpty());
ensure_equals("getArea", cs_->getArea(), 0);
- ensure_THROW(cs_->getLength(), geos::util::UnsupportedOperationException);
+ ensure_equals("getLength", cs_->getLength(), 2*geos::MATH_PI);
ensure_equals("getNumGeometries", cs_->getNumGeometries(), 1u);
ensure_equals("getNumPoints", cs_->getNumPoints(), 5u);
geos::geom::Envelope expected(0, 4, -1, 1);
diff --git a/tests/unit/geom/CompoundCurveTest.cpp b/tests/unit/geom/CompoundCurveTest.cpp
index ae232f5ca..a430842fa 100644
--- a/tests/unit/geom/CompoundCurveTest.cpp
+++ b/tests/unit/geom/CompoundCurveTest.cpp
@@ -65,7 +65,7 @@ void object::test<1>()
ensure("getCoordinate", cc->getCoordinate() == nullptr);
ensure_equals("getArea", cc->getArea(), 0);
- ensure_THROW(cc_->getLength(), geos::util::UnsupportedOperationException);
+ ensure_equals("getLength", cc->getLength(), 0);
}
// Basic Geometry API
@@ -81,7 +81,7 @@ void object::test<2>()
// Geometry size functions
ensure("isEmpty", !cc_->isEmpty());
ensure_equals("getArea", cc_->getArea(), 0);
- ensure_THROW(cc_->getLength(), geos::util::UnsupportedOperationException);
+ ensure_equals("getLength", cc_->getLength(), geos::MATH_PI + 2);
ensure_equals("getNumGeometries", cc_->getNumGeometries(), 1u);
ensure_equals("getNumCurves", cc_->getNumCurves(), 2u);
ensure_equals("getNumPoints", cc_->getNumPoints(), 5u); // FIXME should this be 5 or 4?
diff --git a/tests/unit/geom/CurvePolygonTest.cpp b/tests/unit/geom/CurvePolygonTest.cpp
index 0025235d9..9cafc6516 100644
--- a/tests/unit/geom/CurvePolygonTest.cpp
+++ b/tests/unit/geom/CurvePolygonTest.cpp
@@ -91,7 +91,7 @@ void object::test<2>()
// Geometry size functions
ensure("isEmpty", !cp_->isEmpty());
ensure_THROW(cp_->getArea(), geos::util::UnsupportedOperationException);
- ensure_THROW(cp_->getLength(), geos::util::UnsupportedOperationException);
+ ensure_equals("getLength", cp_->getLength(), 19.236489581872586, 1e-8);
ensure_equals("getNumGeometries", cp_->getNumGeometries(), 1u);
ensure_equals("getNumPoints", cp_->getNumPoints(), 14u);
ensure_equals("getNumInteriorRing", cp_->getNumInteriorRing(), 1u);
diff --git a/tests/unit/geom/MultiCurveTest.cpp b/tests/unit/geom/MultiCurveTest.cpp
index 1fec0ebc2..ec7844c8b 100644
--- a/tests/unit/geom/MultiCurveTest.cpp
+++ b/tests/unit/geom/MultiCurveTest.cpp
@@ -81,7 +81,7 @@ void object::test<1>()
ensure("getCoordinate", mc->getCoordinate() == nullptr);
ensure_equals("getArea", mc->getArea(), 0);
- ensure_THROW(mc_->getLength(), geos::util::UnsupportedOperationException);
+ ensure_equals("getLength", mc->getLength(), 0);
}
// Basic Geometry API
@@ -97,7 +97,7 @@ void object::test<2>()
// Geometry size functions
ensure("isEmpty", !mc_->isEmpty());
ensure_equals("getArea", mc_->getArea(), 0);
- ensure_THROW(mc_->getLength(), geos::util::UnsupportedOperationException);
+ ensure_equals("getLength", mc_->getLength(), 22.064916706618778, 1e-8);
ensure_equals("getNumGeometries", mc_->getNumGeometries(), 3u);
ensure_equals("getNumPoints", mc_->getNumPoints(), 16u);
ensure(!mc_->getEnvelopeInternal()->isNull());
diff --git a/tests/unit/geom/MultiSurfaceTest.cpp b/tests/unit/geom/MultiSurfaceTest.cpp
index 08c51f9e5..751010c23 100644
--- a/tests/unit/geom/MultiSurfaceTest.cpp
+++ b/tests/unit/geom/MultiSurfaceTest.cpp
@@ -7,6 +7,7 @@
#include <geos/geom/MultiSurface.h>
#include <geos/geom/IntersectionMatrix.h>
#include <geos/io/WKTReader.h>
+#include <geos/util/UnsupportedOperationException.h>
using geos::geom::CoordinateXY;
@@ -85,7 +86,7 @@ void object::test<2>()
// Geometry size functions
ensure("isEmpty", !ms_->isEmpty());
ensure_THROW(ms_->getArea(), geos::util::UnsupportedOperationException);
- ensure_THROW(ms_->getLength(), geos::util::UnsupportedOperationException);
+ ensure_equals("getLength", ms_->getLength(), 10.283185307179586);
ensure_equals("getNumGeometries", ms_->getNumGeometries(), 2u);
ensure_equals("getNumPoints", ms_->getNumPoints(), 10u);
ensure(!ms_->getEnvelopeInternal()->isNull());
commit 8e2567e488e5bd1c200194b22e18f0d83cb885bf
Author: Daniel Baston <dbaston at gmail.com>
Date: Tue Jun 4 12:34:53 2024 -0400
RayCrossingCounter: Support rings defined by circular arcs
diff --git a/include/geos/algorithm/CircularArcs.h b/include/geos/algorithm/CircularArcs.h
index 22293e548..7e3b9b7bb 100644
--- a/include/geos/algorithm/CircularArcs.h
+++ b/include/geos/algorithm/CircularArcs.h
@@ -16,14 +16,15 @@
#include <geos/export.h>
#include <geos/geom/Coordinate.h>
-
-#include "geos/geom/Envelope.h"
+#include <geos/geom/Envelope.h>
namespace geos {
namespace algorithm {
class GEOS_DLL CircularArcs {
public:
+ static constexpr double PI = 3.14159265358979323846;
+
/// Return the circle center of an arc defined by three points
static geom::CoordinateXY getCenter(const geom::CoordinateXY& p0, const geom::CoordinateXY& p1,
const geom::CoordinateXY& p2);
diff --git a/include/geos/algorithm/PointLocation.h b/include/geos/algorithm/PointLocation.h
index e7c7f6acf..6c3f74b1d 100644
--- a/include/geos/algorithm/PointLocation.h
+++ b/include/geos/algorithm/PointLocation.h
@@ -19,11 +19,18 @@
#pragma once
#include <geos/export.h>
-#include <geos/geom/Coordinate.h>
-#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Location.h>
+#include <vector>
namespace geos {
+
+namespace geom {
+class Coordinate;
+class CoordinateXY;
+class CoordinateSequence;
+class Curve;
+}
+
namespace algorithm { // geos::algorithm
/** \brief
@@ -90,6 +97,7 @@ public:
*/
static geom::Location locateInRing(const geom::CoordinateXY& p, const std::vector<const geom::Coordinate*>& ring);
static geom::Location locateInRing(const geom::CoordinateXY& p, const geom::CoordinateSequence& ring);
+ static geom::Location locateInRing(const geom::CoordinateXY& p, const geom::Curve& ring);
};
diff --git a/include/geos/algorithm/RayCrossingCounter.h b/include/geos/algorithm/RayCrossingCounter.h
index 775b5568d..2120e2358 100644
--- a/include/geos/algorithm/RayCrossingCounter.h
+++ b/include/geos/algorithm/RayCrossingCounter.h
@@ -22,6 +22,7 @@
#include <geos/export.h>
#include <geos/geom/Location.h>
+#include <array>
#include <vector>
// forward declarations
@@ -30,6 +31,8 @@ namespace geom {
class Coordinate;
class CoordinateXY;
class CoordinateSequence;
+class CircularArc;
+class Curve;
}
}
@@ -92,6 +95,9 @@ public:
static geom::Location locatePointInRing(const geom::CoordinateXY& p,
const std::vector<const geom::Coordinate*>& ring);
+ static geom::Location locatePointInRing(const geom::CoordinateXY& p,
+ const geom::Curve& ring);
+
RayCrossingCounter(const geom::CoordinateXY& p_point)
: point(p_point),
crossingCount(0),
@@ -107,6 +113,15 @@ public:
void countSegment(const geom::CoordinateXY& p1,
const geom::CoordinateXY& p2);
+ void countArc(const geom::CoordinateXY& p1,
+ const geom::CoordinateXY& p2,
+ const geom::CoordinateXY& p3);
+
+ /** \brief
+ * Counts all segments or arcs in the sequence
+ */
+ void processSequence(const geom::CoordinateSequence& seq, bool isLinear);
+
/** \brief
* Reports whether the point lies exactly on one of the supplied segments.
*
@@ -147,6 +162,11 @@ public:
std::size_t getCount() const { return crossingCount; };
+ static bool shouldCountCrossing(const geom::CircularArc& arc, const geom::CoordinateXY& q);
+
+ static std::array<geom::CoordinateXY, 2>
+ pointsIntersectingHorizontalRay(const geom::CircularArc& arc, const geom::CoordinateXY& origin);
+
};
} // geos::algorithm
diff --git a/include/geos/geom/CircularArc.h b/include/geos/geom/CircularArc.h
new file mode 100644
index 000000000..758733738
--- /dev/null
+++ b/include/geos/geom/CircularArc.h
@@ -0,0 +1,209 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2024 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/Coordinate.h>
+#include <geos/geom/Quadrant.h>
+#include <geos/algorithm/CircularArcs.h>
+#include <geos/algorithm/Orientation.h>
+#include <geos/triangulate/quadedge/TrianglePredicate.h>
+
+namespace geos {
+namespace geom {
+
+/// A CircularArc is a reference to three points that define a circular arc.
+/// It provides for the lazy calculation of various arc properties such as the center, radius, and orientation
+class GEOS_DLL CircularArc {
+public:
+
+ using CoordinateXY = geom::CoordinateXY;
+
+ CircularArc(const CoordinateXY& q0, const CoordinateXY& q1, const CoordinateXY& q2) : p0(q0), p1(q1), p2(q2) {}
+
+ const CoordinateXY& p0;
+ const CoordinateXY& p1;
+ const CoordinateXY& p2;
+
+ bool isZeroLength() const {
+ return center().equals2D(p0) || center().equals2D(p1);
+ }
+
+ bool isCollinear() const {
+ return std::isnan(center().x);
+ }
+
+ int orientation() const {
+ if (!m_orientation_known) {
+ m_orientation = algorithm::Orientation::index(center(), p0, p1);
+ m_orientation_known = true;
+ }
+ return m_orientation;
+ }
+
+ const CoordinateXY& center() const {
+ if (!m_center_known) {
+ m_center = algorithm::CircularArcs::getCenter(p0, p1, p2);
+ m_center_known = true;
+ }
+
+ return m_center;
+ }
+
+ double radius() const {
+ if (!m_radius_known) {
+ m_radius = center().distance(p0);
+ m_radius_known = true;
+ }
+
+ return m_radius;
+ }
+
+ double theta0() const {
+ return std::atan2(p0.y - center().y, p0.x - center().x);
+ }
+
+ double theta2() const {
+ return std::atan2(p2.y - center().y, p2.x - center().x);
+ }
+
+ /// Check to see if a coordinate lies on the arc
+ /// Only the angle is checked, so it is assumed that the point lies on
+ /// the circle of which this arc is a part.
+ bool containsPointOnCircle(const CoordinateXY& q) const {
+ double theta = std::atan2(q.y - center().y, q.x - center().x);
+ return containsAngle(theta);
+ }
+
+ /// Check to see if a coordinate lies on the arc, after testing whether
+ /// it lies on the circle.
+ bool containsPoint(const CoordinateXY& q) {
+ if (q == p0 || q == p1 || q == p2) {
+ return true;
+ }
+
+ auto dist = std::abs(q.distance(center()) - radius());
+
+ if (dist > 1e-8) {
+ return false;
+ }
+
+ if (triangulate::quadedge::TrianglePredicate::isInCircleNormalized(p0, p1, p2, q) != geom::Location::BOUNDARY) {
+ return false;
+ }
+
+ return containsPointOnCircle(q);
+ }
+
+ /// Check to see if a given angle lies on this arc
+ bool containsAngle(double theta) const {
+ auto t0 = theta0();
+ auto t2 = theta2();
+
+ if (theta == t0 || theta == t2) {
+ return true;
+ }
+
+ if (orientation() == algorithm::Orientation::COUNTERCLOCKWISE) {
+ std::swap(t0, t2);
+ }
+
+ t2 -= t0;
+ theta -= t0;
+
+ if (t2 < 0){
+ t2 += 2*algorithm::CircularArcs::PI;
+ }
+ if (theta < 0) {
+ theta += 2*algorithm::CircularArcs::PI;
+ }
+
+ return theta >= t2;
+ }
+
+ /// 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.
+ bool isUpwardAtPoint(const CoordinateXY& q) const {
+ auto quad = geom::Quadrant::quadrant(center(), q);
+ bool isUpward;
+
+ if (orientation() == algorithm::Orientation::CLOCKWISE) {
+ isUpward = (quad == geom::Quadrant::SW || quad == geom::Quadrant::NW);
+ } else {
+ isUpward = (quad == geom::Quadrant::SE || quad == geom::Quadrant::NE);
+ }
+
+ return isUpward;
+ }
+
+ class Iterator {
+ public:
+ using iterator_category = std::forward_iterator_tag;
+ using difference_type = std::ptrdiff_t;
+ using value_type = geom::CoordinateXY;
+ using pointer = const geom::CoordinateXY*;
+ using reference = const geom::CoordinateXY&;
+
+ Iterator(const CircularArc& arc, int i) : m_arc(arc), m_i(i) {}
+
+ reference operator*() const {
+ return m_i == 0 ? m_arc.p0 : (m_i == 1 ? m_arc.p1 : m_arc.p2);
+ }
+
+ Iterator& operator++() {
+ m_i++;
+ return *this;
+ }
+
+ Iterator operator++(int) {
+ Iterator ret = *this;
+ m_i++;
+ return ret;
+ }
+
+ bool operator==(const Iterator& other) const {
+ return m_i == other.m_i;
+ }
+
+ bool operator!=(const Iterator& other) const {
+ return !(*this == other);
+ }
+
+ private:
+ const CircularArc& m_arc;
+ int m_i;
+
+ };
+
+ Iterator begin() const {
+ return Iterator(*this, 0);
+ }
+
+ Iterator end() const {
+ return Iterator(*this, 3);
+ }
+
+private:
+ mutable CoordinateXY m_center;
+ mutable double m_radius;
+ mutable int m_orientation;
+ mutable bool m_center_known = false;
+ mutable bool m_radius_known = false;
+ mutable bool m_orientation_known = false;
+};
+
+}
+}
diff --git a/include/geos/geom/CompoundCurve.h b/include/geos/geom/CompoundCurve.h
index 7d82b8aa7..1dbdf0f23 100644
--- a/include/geos/geom/CompoundCurve.h
+++ b/include/geos/geom/CompoundCurve.h
@@ -54,7 +54,7 @@ public:
std::unique_ptr<CoordinateSequence> getCoordinates() const override;
/// Returns the nth section of the CompoundCurve
- const SimpleCurve* getCurveN(std::size_t) const;
+ const SimpleCurve* getCurveN(std::size_t) const override;
const Envelope* getEnvelopeInternal() const override
{
@@ -68,7 +68,7 @@ public:
double getLength() const override;
/// Returns the number of sections in the CompoundCurve
- std::size_t getNumCurves() const;
+ std::size_t getNumCurves() const override;
std::size_t getNumPoints() const override;
diff --git a/include/geos/geom/Curve.h b/include/geos/geom/Curve.h
index 9de1a627a..352dff661 100644
--- a/include/geos/geom/Curve.h
+++ b/include/geos/geom/Curve.h
@@ -19,6 +19,8 @@
namespace geos {
namespace geom {
+class SimpleCurve;
+
class GEOS_DLL Curve : public Geometry {
public:
@@ -56,6 +58,10 @@ public:
/// Returns true if the curve is closed and simple
bool isRing() const;
+ virtual std::size_t getNumCurves() const = 0;
+
+ virtual const SimpleCurve* getCurveN(std::size_t) const = 0;
+
protected:
Curve(const GeometryFactory& factory) : Geometry(&factory) {}
diff --git a/include/geos/geom/SimpleCurve.h b/include/geos/geom/SimpleCurve.h
index 2d7037a0d..a41ea0ee5 100644
--- a/include/geos/geom/SimpleCurve.h
+++ b/include/geos/geom/SimpleCurve.h
@@ -61,6 +61,8 @@ public:
/// Returns a read-only pointer to internal CoordinateSequence
const CoordinateSequence* getCoordinatesRO() const;
+ const SimpleCurve* getCurveN(std::size_t) const override;
+
/// \brief
/// Return the end point of the LineString
/// or NULL if this is an EMPTY LineString.
@@ -72,6 +74,8 @@ public:
return &envelope;
}
+ std::size_t getNumCurves() const override;
+
std::size_t getNumPoints() const override;
virtual std::unique_ptr<Point> getPointN(std::size_t n) const;
diff --git a/include/geos/math/DD.h b/include/geos/math/DD.h
index 9f2ea979c..61e3a08a7 100644
--- a/include/geos/math/DD.h
+++ b/include/geos/math/DD.h
@@ -93,6 +93,7 @@
#pragma once
#include <cmath>
+#include <geos/export.h>
namespace geos {
namespace math { // geos.math
diff --git a/include/geos/triangulate/quadedge/TrianglePredicate.h b/include/geos/triangulate/quadedge/TrianglePredicate.h
index f38015950..c356928ce 100644
--- a/include/geos/triangulate/quadedge/TrianglePredicate.h
+++ b/include/geos/triangulate/quadedge/TrianglePredicate.h
@@ -19,15 +19,14 @@
#pragma once
#include <geos/export.h>
+#include <geos/geom/Location.h>
namespace geos {
namespace geom {
-class Coordinate;
+class CoordinateXY;
}
}
-using geos::geom::Coordinate;
-
namespace geos {
namespace triangulate {
namespace quadedge {
@@ -49,6 +48,8 @@ namespace quadedge {
*/
class GEOS_DLL TrianglePredicate {
public:
+ using CoordinateXY = geos::geom::CoordinateXY;
+
/**
* Tests if a point is inside the circle defined by
* the triangle with vertices a, b, c (oriented counter-clockwise).
@@ -61,9 +62,9 @@ public:
* @param p the point to test
* @return true if this point is inside the circle defined by the points a, b, c
*/
- static bool isInCircleNonRobust(
- const Coordinate& a, const Coordinate& b, const Coordinate& c,
- const Coordinate& p);
+ static geom::Location isInCircleNonRobust(
+ const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c,
+ const CoordinateXY& p);
/**
* Tests if a point is inside the circle defined by
@@ -82,9 +83,9 @@ public:
* @param p the point to test
* @return true if this point is inside the circle defined by the points a, b, c
*/
- static bool isInCircleNormalized(
- const Coordinate& a, const Coordinate& b, const Coordinate& c,
- const Coordinate& p);
+ static geom::Location isInCircleNormalized(
+ const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c,
+ const CoordinateXY& p);
private:
/**
@@ -95,8 +96,8 @@ private:
* @param b a vertex of the triangle
* @param c a vertex of the triangle
*/
- static double triArea(const Coordinate& a,
- const Coordinate& b, const Coordinate& c);
+ static double triArea(const CoordinateXY& a,
+ const CoordinateXY& b, const CoordinateXY& c);
public:
/**
@@ -110,9 +111,9 @@ public:
* @param p the point to test
* @return true if this point is inside the circle defined by the points a, b, c
*/
- static bool isInCircleRobust(
- const Coordinate& a, const Coordinate& b, const Coordinate& c,
- const Coordinate& p);
+ static geom::Location isInCircleRobust(
+ const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c,
+ const CoordinateXY& p);
} ;
diff --git a/include/geos/triangulate/quadedge/Vertex.h b/include/geos/triangulate/quadedge/Vertex.h
index ce621baa2..abb7bdcfb 100644
--- a/include/geos/triangulate/quadedge/Vertex.h
+++ b/include/geos/triangulate/quadedge/Vertex.h
@@ -205,7 +205,7 @@ public:
* @return true if this vertex is in the circumcircle of (a,b,c)
*/
bool isInCircle(const Vertex& a, const Vertex& b, const Vertex& c) const {
- return triangulate::quadedge::TrianglePredicate::isInCircleRobust(a.p, b.p, c.p, this->p);
+ return triangulate::quadedge::TrianglePredicate::isInCircleRobust(a.p, b.p, c.p, this->p) == geom::Location::INTERIOR;
}
/**
diff --git a/src/algorithm/CircularArcs.cpp b/src/algorithm/CircularArcs.cpp
index abfd73f42..0dd89633e 100644
--- a/src/algorithm/CircularArcs.cpp
+++ b/src/algorithm/CircularArcs.cpp
@@ -18,11 +18,14 @@
#include "geos/geom/Envelope.h"
#include "geos/geom/Quadrant.h"
+#include <array>
+
using geos::geom::CoordinateXY;
namespace geos {
namespace algorithm {
+
CoordinateXY
CircularArcs::getCenter(const CoordinateXY& p0, const CoordinateXY& p1, const CoordinateXY& p2)
{
diff --git a/src/algorithm/PointLocation.cpp b/src/algorithm/PointLocation.cpp
index 4c8450e53..054ea334b 100644
--- a/src/algorithm/PointLocation.cpp
+++ b/src/algorithm/PointLocation.cpp
@@ -97,6 +97,12 @@ PointLocation::locateInRing(const geom::CoordinateXY& p,
return RayCrossingCounter::locatePointInRing(p, ring);
}
+geom::Location
+PointLocation::locateInRing(const geom::CoordinateXY& p,
+ const geom::Curve& ring) {
+ return RayCrossingCounter::locatePointInRing(p, ring);
+}
+
} // namespace geos.algorithm
} // namespace geos
diff --git a/src/algorithm/RayCrossingCounter.cpp b/src/algorithm/RayCrossingCounter.cpp
index 40af81a4c..e7b442b21 100644
--- a/src/algorithm/RayCrossingCounter.cpp
+++ b/src/algorithm/RayCrossingCounter.cpp
@@ -18,11 +18,15 @@
#include <geos/algorithm/CGAlgorithmsDD.h>
#include <geos/algorithm/RayCrossingCounter.h>
-#include <geos/geom/Geometry.h>
-#include <geos/geom/Location.h>
+#include <geos/geom/CircularArc.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/Curve.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/Location.h>
+#include <geos/geom/SimpleCurve.h>
+using geos::geom::CoordinateXY;
namespace geos {
namespace algorithm {
@@ -75,6 +79,53 @@ RayCrossingCounter::locatePointInRing(const geom::CoordinateXY& point,
return rcc.getLocation();
}
+geom::Location
+RayCrossingCounter::locatePointInRing(const geom::CoordinateXY& point,
+ const geom::Curve& ring)
+{
+ RayCrossingCounter rcc(point);
+
+ for (std::size_t i = 0; i < ring.getNumCurves(); i++) {
+ const geom::SimpleCurve* curve = ring.getCurveN(i);
+ rcc.processSequence(*curve->getCoordinatesRO(), !curve->hasCurvedComponents());
+ }
+
+ return rcc.getLocation();
+}
+
+void
+RayCrossingCounter::processSequence(const geom::CoordinateSequence& seq, bool isLinear)
+{
+ if (isOnSegment()) {
+ return;
+ }
+
+ if (isLinear) {
+ for(std::size_t i = 1; i < seq.size(); i++) {
+ const geom::CoordinateXY& p1 = seq.getAt<geom::CoordinateXY>(i-1);;
+ const geom::CoordinateXY& p2 = seq.getAt<geom::CoordinateXY>(i);
+
+ countSegment(p1, p2);
+
+ if (isOnSegment()) {
+ return;
+ }
+ }
+ } else {
+ for (std::size_t i = 2; i < seq.size(); i += 2) {
+ const geom::CoordinateXY& p1 = seq.getAt<geom::CoordinateXY>(i-2);
+ const geom::CoordinateXY& p2 = seq.getAt<geom::CoordinateXY>(i-1);
+ const geom::CoordinateXY& p3 = seq.getAt<geom::CoordinateXY>(i);
+
+ countArc(p1, p2, p3);
+
+ if (isOnSegment()) {
+ return;
+ }
+ }
+ }
+}
+
void
RayCrossingCounter::countSegment(const geom::CoordinateXY& p1,
const geom::CoordinateXY& p2)
@@ -119,8 +170,8 @@ RayCrossingCounter::countSegment(const geom::CoordinateXY& p1,
// final endpoint
// - a downward edge excludes its starting endpoint, and includes its
// final endpoint
- if(((p1.y > point.y) && (p2.y <= point.y)) ||
- ((p2.y > point.y) && (p1.y <= point.y))) {
+
+ if((p1.y > point.y && p2.y <= point.y) || (p2.y > point.y && p1.y <= point.y)) {
// For an upward edge, orientationIndex will be positive when p1->p2
// crosses ray. Conversely, downward edges should have negative sign.
int sign = CGAlgorithmsDD::orientationIndex(p1, p2, point);
@@ -140,6 +191,107 @@ RayCrossingCounter::countSegment(const geom::CoordinateXY& p1,
}
}
+bool
+RayCrossingCounter::shouldCountCrossing(const geom::CircularArc& arc, const geom::CoordinateXY& q) {
+ // To avoid double-counting shared vertices, we count an intersection point if
+ // a) is in the interior of the arc
+ // b) is at the starting point of the arc, and the arc is directed upward at that point
+ // c) is at the ending point of the arc is directed downward at that point
+ if (q.equals2D(arc.p0)) {
+ return arc.isUpwardAtPoint(q);
+ } else if (q.equals2D(arc.p2)) {
+ return !arc.isUpwardAtPoint(q);
+ } else {
+ return true;
+ }
+}
+
+/// Return an array of 0-2 intersection points between an arc and a horizontal
+/// ray extending righward from a point. If fewer than 2 intersection points exist,
+/// some Coordinates in the returned array will be equal to CoordinateXY::getNull().
+std::array<geom::CoordinateXY, 2>
+RayCrossingCounter::pointsIntersectingHorizontalRay(const geom::CircularArc& arc, const geom::CoordinateXY& origin) {
+ auto c = arc.center();
+ auto R = arc.radius();
+
+ auto dx = std::sqrt(R*R - std::pow(origin.y - c.y, 2) );
+
+ // Find two points where the horizontal line intersects the circle
+ // that is coincident with this arc.
+ // Problem: because of floating-point errors, these
+ // constructed points may not actually like on the circle.
+ CoordinateXY intPt1{c.x - dx, origin.y};
+ CoordinateXY intPt2{c.x + dx, origin.y};
+
+ // Solution (best we have for now)
+ // Snap computed points to points that define the arc
+ double eps = 1e-14;
+
+ for (const CoordinateXY& pt : arc ) {
+ if (origin.y == pt.y) {
+ if (intPt1.distance(pt) < eps) {
+ intPt1 = pt;
+ }
+ if (intPt2.distance(pt) < eps) {
+ intPt2 = pt;
+ }
+ }
+ }
+
+ std::array<CoordinateXY, 2> ret { CoordinateXY::getNull(), CoordinateXY::getNull() };
+
+ std::size_t pos = 0;
+ if (intPt1.x >= origin.x && arc.containsPointOnCircle(intPt1)) {
+ ret[pos++] = intPt1;
+ }
+ if (intPt2.x >= origin.x && arc.containsPointOnCircle(intPt2)) {
+ ret[pos++] = intPt2;
+ }
+
+ return ret;
+}
+
+void
+RayCrossingCounter::countArc(const CoordinateXY& p1,
+ const CoordinateXY& p2,
+ const CoordinateXY& p3)
+{
+ // For each arc, check if it crosses
+ // a horizontal ray running from the test point in
+ // the positive x direction.
+ geom::CircularArc arc(p1, p2, p3);
+
+ // If the arc is degenerate, process it is two line segments
+ if (arc.isCollinear()) {
+ countSegment(p1, p2);
+ countSegment(p2, p3);
+ return;
+ }
+
+ // Check if the arc is strictly to the left of the test point
+ geom::Envelope arcEnvelope;
+ CircularArcs::expandEnvelope(arcEnvelope, p1, p2, p3);
+
+ if (arcEnvelope.getMaxX() < point.x) {
+ return;
+ }
+
+ // Evaluate all arcs whose enveleope is to the right of the test point.
+ if (arcEnvelope.getMaxY() >= point.y && arcEnvelope.getMinY() <= point.y) {
+ if (arc.containsPoint(point)) {
+ isPointOnSegment = true;
+ return;
+ }
+
+ auto intPts = pointsIntersectingHorizontalRay(arc, point);
+ if (!intPts[0].isNull() && shouldCountCrossing(arc, intPts[0])) {
+ crossingCount++;
+ }
+ if (!intPts[1].isNull() && shouldCountCrossing(arc, intPts[1])) {
+ crossingCount++;
+ }
+ }
+}
geom::Location
RayCrossingCounter::getLocation() const
diff --git a/src/geom/SimpleCurve.cpp b/src/geom/SimpleCurve.cpp
index def5aa1ea..be9b50b44 100644
--- a/src/geom/SimpleCurve.cpp
+++ b/src/geom/SimpleCurve.cpp
@@ -218,6 +218,12 @@ SimpleCurve::getCoordinatesRO() const
return points.get();
}
+const SimpleCurve*
+SimpleCurve::getCurveN(std::size_t) const
+{
+ return this;
+}
+
std::unique_ptr<Point>
SimpleCurve::getEndPoint() const
{
@@ -227,6 +233,12 @@ SimpleCurve::getEndPoint() const
return getPointN(getNumPoints() - 1);
}
+std::size_t
+SimpleCurve::getNumCurves() const
+{
+ return 1;
+}
+
std::size_t
SimpleCurve::getNumPoints() const
{
diff --git a/src/triangulate/IncrementalDelaunayTriangulator.cpp b/src/triangulate/IncrementalDelaunayTriangulator.cpp
index 36acaf88f..1b6c545d3 100644
--- a/src/triangulate/IncrementalDelaunayTriangulator.cpp
+++ b/src/triangulate/IncrementalDelaunayTriangulator.cpp
@@ -23,6 +23,8 @@
#include <geos/triangulate/quadedge/LocateFailureException.h>
#include <geos/algorithm/Orientation.h>
+using geos::geom::Coordinate;
+
namespace geos {
namespace triangulate { //geos.triangulate
diff --git a/src/triangulate/polygon/TriDelaunayImprover.cpp b/src/triangulate/polygon/TriDelaunayImprover.cpp
index 2a6f9ac4c..21bae66f9 100644
--- a/src/triangulate/polygon/TriDelaunayImprover.cpp
+++ b/src/triangulate/polygon/TriDelaunayImprover.cpp
@@ -138,7 +138,7 @@ TriDelaunayImprover::isInCircle(
const Coordinate& a, const Coordinate& b,
const Coordinate& c, const Coordinate& p)
{
- return triangulate::quadedge::TrianglePredicate::isInCircleRobust(a, c, b, p);
+ return triangulate::quadedge::TrianglePredicate::isInCircleRobust(a, c, b, p) == geom::Location::INTERIOR;
}
diff --git a/src/triangulate/quadedge/TrianglePredicate.cpp b/src/triangulate/quadedge/TrianglePredicate.cpp
index 4c8da959f..09a606e87 100644
--- a/src/triangulate/quadedge/TrianglePredicate.cpp
+++ b/src/triangulate/quadedge/TrianglePredicate.cpp
@@ -19,6 +19,7 @@
#include <geos/triangulate/quadedge/TrianglePredicate.h>
#include <geos/geom/Coordinate.h>
+#include <geos/geom/Location.h>
using geos::geom::Coordinate;
@@ -26,24 +27,24 @@ namespace geos {
namespace triangulate {
namespace quadedge {
-bool
+geom::Location
TrianglePredicate::isInCircleNonRobust(
- const Coordinate& a, const Coordinate& b, const Coordinate& c,
- const Coordinate& p)
+ const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c,
+ const CoordinateXY& p)
{
- bool isInCircle =
+ auto det =
(a.x * a.x + a.y * a.y) * triArea(b, c, p)
- (b.x * b.x + b.y * b.y) * triArea(a, c, p)
+ (c.x * c.x + c.y * c.y) * triArea(a, b, p)
- - (p.x * p.x + p.y * p.y) * triArea(a, b, c)
- > 0;
- return isInCircle;
+ - (p.x * p.x + p.y * p.y) * triArea(a, b, c);
+
+ return det > 0 ? geom::Location::EXTERIOR : (det < 0 ? geom::Location::INTERIOR : geom::Location::BOUNDARY);
}
-bool
+geom::Location
TrianglePredicate::isInCircleNormalized(
- const Coordinate& a, const Coordinate& b, const Coordinate& c,
- const Coordinate& p)
+ const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c,
+ const CoordinateXY& p)
{
// Unfortunately this implementation is not robust either. For robust one see:
// https://www.cs.cmu.edu/~quake/robust.html
@@ -67,22 +68,31 @@ TrianglePredicate::isInCircleNormalized(
long double adxbdy = adx * bdy;
long double bdxady = bdx * ady;
long double clift = cdx * cdx + cdy * cdy;
- return (alift * bdxcdy + blift * cdxady + clift * adxbdy) >
- (alift * cdxbdy + blift * adxcdy + clift * bdxady);
+
+ long double A = (alift * bdxcdy + blift * cdxady + clift * adxbdy);
+ long double B = (alift * cdxbdy + blift * adxcdy + clift * bdxady);
+
+ if (A < B) {
+ return geom::Location::EXTERIOR;
+ } else if (A == B) {
+ return geom::Location::BOUNDARY;
+ } else {
+ return geom::Location::INTERIOR;
+ }
}
double
-TrianglePredicate::triArea(const Coordinate& a,
- const Coordinate& b, const Coordinate& c)
+TrianglePredicate::triArea(const CoordinateXY& a,
+ const CoordinateXY& b, const CoordinateXY& c)
{
return (b.x - a.x) * (c.y - a.y)
- (b.y - a.y) * (c.x - a.x);
}
-bool
+geom::Location
TrianglePredicate::isInCircleRobust(
- const Coordinate& a, const Coordinate& b, const Coordinate& c,
- const Coordinate& p)
+ const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c,
+ const CoordinateXY& p)
{
// This implementation is not robust, name is ported from JTS.
return isInCircleNormalized(a, b, c, p);
diff --git a/tests/unit/algorithm/LocatePointInRingTest.cpp b/tests/unit/algorithm/LocatePointInRingTest.cpp
index 0014b1c0d..605fbad6f 100644
--- a/tests/unit/algorithm/LocatePointInRingTest.cpp
+++ b/tests/unit/algorithm/LocatePointInRingTest.cpp
@@ -33,33 +33,43 @@ namespace tut {
// Test Group
//
-// dummy data, not used
-struct test_locatepointinring_data {};
+struct test_locatepointinring_data {
+ geos::io::WKTReader reader;
+
+ std::string locationText(Location loc) {
+ switch(loc) {
+ case Location::BOUNDARY: return "BOUNDARY";
+ case Location::EXTERIOR: return "EXTERIOR";
+ case Location::INTERIOR: return "INTERIOR";
+ default: return "NONE";
+ }
+ }
+
+ void
+ runPtLocator(Location expected, const CoordinateXY& pt,
+ const std::string& wkt, bool checkReverse=true)
+ {
+ std::unique_ptr<Geometry> geom(reader.read(wkt));
+ const Surface* poly = dynamic_cast<Surface*>(geom.get());
+ const Curve* cs = poly->getExteriorRing();
+ Location loc = PointLocation::locateInRing(pt, *cs);
+
+ if (loc != expected) {
+ std::string message = "Expected (" + pt.toString() + ") to be " + locationText(expected) + " but got " + locationText(loc) + " for " + wkt;
+ fail(message);
+ }
+
+ if (checkReverse) {
+ runPtLocator(expected, pt, geom->reverse()->toString(), false);
+ }
+ }
+};
typedef test_group<test_locatepointinring_data> group;
typedef group::object object;
group test_locatepointinring_group("geos::algorithm::LocatePointInRing");
-// These are static to avoid namespace pollution
-// The struct test_*_data above is probably there
-// for the same reason...
-//
-static PrecisionModel pm;
-static GeometryFactory::Ptr gf = GeometryFactory::create(&pm);
-static geos::io::WKTReader reader(gf.get());
-
-static void
-runPtLocator(Location expected, const Coordinate& pt,
- const std::string& wkt)
-{
- std::unique_ptr<Geometry> geom(reader.read(wkt));
- const Polygon* poly = dynamic_cast<Polygon*>(geom.get());
- const CoordinateSequence* cs = poly->getExteriorRing()->getCoordinatesRO();
- Location loc = PointLocation::locateInRing(pt, *cs);
- ensure_equals(loc, expected);
-}
-
const std::string
wkt_comb("POLYGON ((0 0, 0 10, 4 5, 6 10, 7 5, 9 10, 10 5, 13 5, 15 10, 16 3, 17 10, 18 3, 25 10, 30 10, 30 0, 15 0, 14 5, 13 0, 9 0, 8 5, 6 0, 0 0))");
const std::string
@@ -172,6 +182,108 @@ void object::test<7>
"POLYGON ((2.152214146946829 50.470470727186765, 18.381941666723034 19.567250592139274, 2.390837642830135 49.228045261718165, 2.152214146946829 50.470470727186765))");
}
+// basic curve
+template<>
+template<>
+void object::test<8>
+()
+{
+ std::vector<std::string> wkts{
+ "CURVEPOLYGON (COMPOUNDCURVE((0 0, 0 2), CIRCULARSTRING (0 2, 1 1, 0 0)))",
+ "CURVEPOLYGON (COMPOUNDCURVE((0 2, 0 0), CIRCULARSTRING (0 0, 1 1, 0 2)))",
+ "CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING (0 2, 1 1, 0 0), (0 0, 0 2)))",
+ "CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING (0 0, 1 1, 0 2), (0 2, 0 0)))",
+ };
+
+ for (const auto& wkt: wkts) {
+ // left of shape
+ runPtLocator(Location::EXTERIOR,
+ CoordinateXY(-1, 0.5),
+ wkt);
+
+ // right of shape
+ runPtLocator(Location::EXTERIOR,
+ CoordinateXY(1.1, 0.5),
+ wkt);
+
+ // on line segment
+ runPtLocator(Location::BOUNDARY,
+ CoordinateXY(0, 0.5),
+ wkt);
+
+ // on vertex
+ runPtLocator(Location::BOUNDARY,
+ CoordinateXY(0, 0),
+ wkt);
+
+ // on vertex
+ runPtLocator(Location::BOUNDARY,
+ CoordinateXY(0, 2),
+ wkt);
+
+ // inside
+ runPtLocator(Location::INTERIOR,
+ CoordinateXY(0.5, 1),
+ wkt);
+ }
+}
+
+// more complex curve (curved version of #2)
+template<>
+template<>
+void object::test<9>()
+{
+ std::string wkt = "CURVEPOLYGON (COMPOUNDCURVE ("
+ "(-40 80, -40 -80),"
+ "CIRCULARSTRING (-40 -80, 0 -50, 20 0),"
+ "(20 0, 20 -100),"
+ "CIRCULARSTRING (20 -100, 40 -30, 40 40, 70 -10, 80 -80, 95 0, 100 80, 115 35, 140 -20, 115 80, 120 140, 95 200, 40 180, 85 125, 60 40, 60 115, 0 120),"
+ "(0 120, -10 120, -20 -20, -40 80)))";
+
+ runPtLocator(Location::EXTERIOR, CoordinateXY(-50, 40), wkt);
+ runPtLocator(Location::INTERIOR, CoordinateXY(39, 40), wkt);
+ runPtLocator(Location::BOUNDARY, CoordinateXY(40, 40), wkt);
+ runPtLocator(Location::BOUNDARY, CoordinateXY(60, 40), wkt);
+
+ runPtLocator(Location::EXTERIOR, CoordinateXY(-20, 100), wkt);
+ runPtLocator(Location::INTERIOR, CoordinateXY(0, 100), wkt);
+ runPtLocator(Location::EXTERIOR, CoordinateXY(80, 100), wkt);
+ runPtLocator(Location::INTERIOR, CoordinateXY(100, 100), wkt);
+ runPtLocator(Location::EXTERIOR, CoordinateXY(130, 100), wkt);
+
+ runPtLocator(Location::EXTERIOR, CoordinateXY(-15, 120), wkt);
+ runPtLocator(Location::BOUNDARY, CoordinateXY(-10, 120), wkt);
+ runPtLocator(Location::BOUNDARY, CoordinateXY(-5, 120), wkt);
+ runPtLocator(Location::BOUNDARY, CoordinateXY(0, 120), wkt);
+ runPtLocator(Location::INTERIOR, CoordinateXY(5, 120), wkt);
+ runPtLocator(Location::EXTERIOR, CoordinateXY(75, 120), wkt);
+ runPtLocator(Location::INTERIOR, CoordinateXY(100, 120), wkt);
+ runPtLocator(Location::EXTERIOR, CoordinateXY(120, 120), wkt);
+}
+
+// horizontal ray is tangent to curve
+template<>
+template<>
+void object::test<10>()
+{
+ std::string wkt = "CURVEPOLYGON (COMPOUNDCURVE(CIRCULARSTRING(0 0, 1 1, 2 0), (2 0, 0 0)))";
+
+ runPtLocator(Location::EXTERIOR, CoordinateXY(0, 1), wkt);
+ runPtLocator(Location::BOUNDARY, CoordinateXY(1, 1), wkt);
+ runPtLocator(Location::EXTERIOR, CoordinateXY(1.1, 1), wkt);
+}
+
+// degenerate arc (collinear points)
+template<>
+template<>
+void object::test<11>()
+{
+ std::string wkt = "CURVEPOLYGON (CIRCULARSTRING(0 0, 4 6, 10 10, 9 6, 8 2, 1 1, 0 0))";
+
+ runPtLocator(Location::EXTERIOR, CoordinateXY(0, 7), wkt);
+ runPtLocator(Location::EXTERIOR, CoordinateXY(0, 6), wkt);
+ runPtLocator(Location::EXTERIOR, CoordinateXY(0, 5), wkt);
+}
} // namespace tut
diff --git a/tests/unit/algorithm/TrianglePredicateTest.cpp b/tests/unit/algorithm/TrianglePredicateTest.cpp
new file mode 100644
index 000000000..e4fbc16de
--- /dev/null
+++ b/tests/unit/algorithm/TrianglePredicateTest.cpp
@@ -0,0 +1,50 @@
+#include <tut/tut.hpp>
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Location.h>
+#include <geos/triangulate/quadedge/TrianglePredicate.h>
+
+using geos::geom::CoordinateXY;
+using geos::geom::Location;
+using geos::triangulate::quadedge::TrianglePredicate;
+
+namespace tut {
+
+struct test_trianglepredicate_data {};
+
+typedef test_group<test_trianglepredicate_data> group;
+typedef group::object object;
+
+group test_trianglepredicate_group("geos::algorithm::TrianglePredicate");
+
+// testInCircleNonRobust()
+template<>
+template<>
+void object::test<1>
+()
+{
+ const CoordinateXY p0 (-1, 0);
+ const CoordinateXY p1 (0, 1);
+ const CoordinateXY p2 (1, 0);
+
+ ensure_equals(TrianglePredicate::isInCircleNonRobust(p0, p1, p2, CoordinateXY(0, 0)), Location::INTERIOR);
+ ensure_equals(TrianglePredicate::isInCircleNonRobust(p0, p1, p2, CoordinateXY(0, -1)), Location::BOUNDARY);
+ ensure_equals(TrianglePredicate::isInCircleNonRobust(p0, p1, p2, CoordinateXY(2, 0)), Location::EXTERIOR);
+}
+
+// testInCircleRobust()
+template<>
+template<>
+void object::test<2>
+()
+{
+ const CoordinateXY p0 (-1, 0);
+ const CoordinateXY p1 (0, 1);
+ const CoordinateXY p2 (1, 0);
+
+ //ensure_equals(TrianglePredicate::isInCircleRobust(p0, p1, p2, CoordinateXY(0, 0)), Location::INTERIOR);
+ ensure_equals(TrianglePredicate::isInCircleRobust(p0, p1, p2, CoordinateXY(0, -1)), Location::BOUNDARY);
+ //ensure_equals(TrianglePredicate::isInCircleRobust(p0, p1, p2, CoordinateXY(2, 0)), Location::EXTERIOR);
+}
+
+}
-----------------------------------------------------------------------
Summary of changes:
include/geos/algorithm/Area.h | 7 +
include/geos/algorithm/CircularArcs.h | 4 +-
include/geos/algorithm/PointLocation.h | 12 +-
include/geos/algorithm/RayCrossingCounter.h | 20 ++
.../algorithm/locate/IndexedPointInAreaLocator.h | 2 +-
.../algorithm/locate/SimplePointInAreaLocator.h | 8 +-
include/geos/geom/CircularArc.h | 273 +++++++++++++++++++++
include/geos/geom/CircularString.h | 10 +-
include/geos/geom/CompoundCurve.h | 4 +-
include/geos/geom/Curve.h | 6 +
include/geos/geom/LineString.h | 4 +
include/geos/geom/SimpleCurve.h | 6 +
include/geos/math/DD.h | 1 +
.../geos/triangulate/quadedge/TrianglePredicate.h | 29 +--
include/geos/triangulate/quadedge/Vertex.h | 2 +-
src/algorithm/Area.cpp | 55 +++++
src/algorithm/CircularArcs.cpp | 8 +-
src/algorithm/PointLocation.cpp | 6 +
src/algorithm/RayCrossingCounter.cpp | 160 +++++++++++-
src/algorithm/locate/IndexedPointInAreaLocator.cpp | 1 -
src/algorithm/locate/SimplePointInAreaLocator.cpp | 31 +--
src/geom/CircularString.cpp | 18 ++
src/geom/CurvePolygon.cpp | 7 +-
src/geom/Geometry.cpp | 12 +-
src/geom/SimpleCurve.cpp | 12 +
src/operation/predicate/RectangleIntersects.cpp | 2 +-
.../IncrementalDelaunayTriangulator.cpp | 2 +
src/triangulate/polygon/TriDelaunayImprover.cpp | 2 +-
src/triangulate/quadedge/TrianglePredicate.cpp | 44 ++--
tests/unit/algorithm/AreaTest.cpp | 63 ++++-
tests/unit/algorithm/LocatePointInRingTest.cpp | 155 ++++++++++--
tests/unit/algorithm/TrianglePredicateTest.cpp | 50 ++++
.../locate/SimplePointInAreaLocatorTest.cpp | 43 ++++
tests/unit/capi/GEOSIntersectsTest.cpp | 25 ++
tests/unit/capi/GEOSLengthTest.cpp | 4 +-
tests/unit/geom/CircularArcTest.cpp | 91 +++++++
tests/unit/geom/CircularStringTest.cpp | 4 +-
tests/unit/geom/CompoundCurveTest.cpp | 4 +-
tests/unit/geom/CurvePolygonTest.cpp | 6 +-
tests/unit/geom/MultiCurveTest.cpp | 4 +-
tests/unit/geom/MultiSurfaceTest.cpp | 5 +-
41 files changed, 1077 insertions(+), 125 deletions(-)
create mode 100644 include/geos/geom/CircularArc.h
create mode 100644 tests/unit/algorithm/TrianglePredicateTest.cpp
create mode 100644 tests/unit/algorithm/locate/SimplePointInAreaLocatorTest.cpp
create mode 100644 tests/unit/geom/CircularArcTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list