From 3e5e56b5ed0b35f463186bdf4f51207123f8a5dd Mon Sep 17 00:00:00 2001 From: Daniel Baston Date: Tue, 4 Jun 2024 12:34:53 -0400 Subject: [PATCH] RayCrossingCounter: Support rings defined by circular arcs --- include/geos/algorithm/CircularArcs.h | 3 +- include/geos/algorithm/PointLocation.h | 12 +- include/geos/algorithm/RayCrossingCounter.h | 20 ++ include/geos/geom/CircularArc.h | 209 ++++++++++++++++++ include/geos/geom/CompoundCurve.h | 4 +- include/geos/geom/Curve.h | 6 + include/geos/geom/SimpleCurve.h | 4 + include/geos/math/DD.h | 1 + .../triangulate/quadedge/TrianglePredicate.h | 25 ++- src/algorithm/CircularArcs.cpp | 3 + src/algorithm/PointLocation.cpp | 6 + src/algorithm/RayCrossingCounter.cpp | 160 +++++++++++++- src/geom/SimpleCurve.cpp | 12 + .../IncrementalDelaunayTriangulator.cpp | 2 + .../quadedge/TrianglePredicate.cpp | 34 ++- .../unit/algorithm/LocatePointInRingTest.cpp | 154 +++++++++++-- .../unit/algorithm/TrianglePredicateTest.cpp | 35 +++ 17 files changed, 635 insertions(+), 55 deletions(-) create mode 100644 include/geos/geom/CircularArc.h create mode 100644 tests/unit/algorithm/TrianglePredicateTest.cpp diff --git a/include/geos/algorithm/CircularArcs.h b/include/geos/algorithm/CircularArcs.h index 22293e5482..b0ec87e073 100644 --- a/include/geos/algorithm/CircularArcs.h +++ b/include/geos/algorithm/CircularArcs.h @@ -16,8 +16,7 @@ #include #include - -#include "geos/geom/Envelope.h" +#include namespace geos { namespace algorithm { diff --git a/include/geos/algorithm/PointLocation.h b/include/geos/algorithm/PointLocation.h index e7c7f6acf7..6c3f74b1d3 100644 --- a/include/geos/algorithm/PointLocation.h +++ b/include/geos/algorithm/PointLocation.h @@ -19,11 +19,18 @@ #pragma once #include -#include -#include #include +#include namespace geos { + +namespace geom { +class Coordinate; +class CoordinateXY; +class CoordinateSequence; +class Curve; +} + namespace algorithm { // geos::algorithm /** \brief @@ -90,6 +97,7 @@ class GEOS_DLL PointLocation { */ static geom::Location locateInRing(const geom::CoordinateXY& p, const std::vector& 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 775b5568d0..2120e2358e 100644 --- a/include/geos/algorithm/RayCrossingCounter.h +++ b/include/geos/algorithm/RayCrossingCounter.h @@ -22,6 +22,7 @@ #include #include +#include #include // forward declarations @@ -30,6 +31,8 @@ namespace geom { class Coordinate; class CoordinateXY; class CoordinateSequence; +class CircularArc; +class Curve; } } @@ -92,6 +95,9 @@ class GEOS_DLL RayCrossingCounter { static geom::Location locatePointInRing(const geom::CoordinateXY& p, const std::vector& 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 @@ class GEOS_DLL RayCrossingCounter { 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 @@ class GEOS_DLL RayCrossingCounter { std::size_t getCount() const { return crossingCount; }; + static bool shouldCountCrossing(const geom::CircularArc& arc, const geom::CoordinateXY& q); + + static std::array + 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 0000000000..c45bdd95da --- /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 +#include +#include +#include +#include +#include + +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*M_PI; + } + if (theta < 0) { + theta += 2*M_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 7d82b8aa7c..1dbdf0f230 100644 --- a/include/geos/geom/CompoundCurve.h +++ b/include/geos/geom/CompoundCurve.h @@ -54,7 +54,7 @@ class GEOS_DLL CompoundCurve : public Curve { std::unique_ptr 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 @@ class GEOS_DLL CompoundCurve : public Curve { 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 9de1a627ae..352dff6617 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 @@ class GEOS_DLL Curve : public Geometry { /// 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 2d7037a0db..a41ea0ee55 100644 --- a/include/geos/geom/SimpleCurve.h +++ b/include/geos/geom/SimpleCurve.h @@ -61,6 +61,8 @@ class GEOS_DLL SimpleCurve : public Curve { /// 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 @@ class GEOS_DLL SimpleCurve : public Curve { return &envelope; } + std::size_t getNumCurves() const override; + std::size_t getNumPoints() const override; virtual std::unique_ptr getPointN(std::size_t n) const; diff --git a/include/geos/math/DD.h b/include/geos/math/DD.h index 9f2ea979c7..61e3a08a72 100644 --- a/include/geos/math/DD.h +++ b/include/geos/math/DD.h @@ -93,6 +93,7 @@ #pragma once #include +#include namespace geos { namespace math { // geos.math diff --git a/include/geos/triangulate/quadedge/TrianglePredicate.h b/include/geos/triangulate/quadedge/TrianglePredicate.h index f380159504..4ba4e39092 100644 --- a/include/geos/triangulate/quadedge/TrianglePredicate.h +++ b/include/geos/triangulate/quadedge/TrianglePredicate.h @@ -19,15 +19,14 @@ #pragma once #include +#include 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). @@ -62,8 +63,8 @@ class GEOS_DLL TrianglePredicate { * @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); + 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 @@ class GEOS_DLL TrianglePredicate { * @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 @@ class GEOS_DLL TrianglePredicate { * @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: /** @@ -111,8 +112,8 @@ class GEOS_DLL TrianglePredicate { * @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); + const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c, + const CoordinateXY& p); } ; diff --git a/src/algorithm/CircularArcs.cpp b/src/algorithm/CircularArcs.cpp index abfd73f42c..0dd89633ee 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 + 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 4c8450e53d..054ea334bc 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 40af81a4cf..e7b442b217 100644 --- a/src/algorithm/RayCrossingCounter.cpp +++ b/src/algorithm/RayCrossingCounter.cpp @@ -18,11 +18,15 @@ #include #include -#include -#include +#include #include #include +#include +#include +#include +#include +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(i-1);; + const geom::CoordinateXY& p2 = seq.getAt(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(i-2); + const geom::CoordinateXY& p2 = seq.getAt(i-1); + const geom::CoordinateXY& p3 = seq.getAt(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 +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 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 def5aa1ea3..be9b50b443 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 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 36acaf88f3..1b6c545d3f 100644 --- a/src/triangulate/IncrementalDelaunayTriangulator.cpp +++ b/src/triangulate/IncrementalDelaunayTriangulator.cpp @@ -23,6 +23,8 @@ #include #include +using geos::geom::Coordinate; + namespace geos { namespace triangulate { //geos.triangulate diff --git a/src/triangulate/quadedge/TrianglePredicate.cpp b/src/triangulate/quadedge/TrianglePredicate.cpp index 4c8da959fb..c7c0eef75e 100644 --- a/src/triangulate/quadedge/TrianglePredicate.cpp +++ b/src/triangulate/quadedge/TrianglePredicate.cpp @@ -19,6 +19,7 @@ #include #include +#include using geos::geom::Coordinate; @@ -28,8 +29,8 @@ namespace quadedge { bool 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 = (a.x * a.x + a.y * a.y) * triArea(b, c, p) @@ -40,10 +41,10 @@ TrianglePredicate::isInCircleNonRobust( return isInCircle; } -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,13 +68,22 @@ 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); @@ -81,11 +91,11 @@ TrianglePredicate::triArea(const Coordinate& a, bool 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); + return isInCircleNormalized(a, b, c, p) == geom::Location::INTERIOR; } } // namespace geos.triangulate.quadedge diff --git a/tests/unit/algorithm/LocatePointInRingTest.cpp b/tests/unit/algorithm/LocatePointInRingTest.cpp index 0014b1c0d8..605fbad6f9 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 geom(reader.read(wkt)); + const Surface* poly = dynamic_cast(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 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 geom(reader.read(wkt)); - const Polygon* poly = dynamic_cast(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 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 0000000000..44469dee09 --- /dev/null +++ b/tests/unit/algorithm/TrianglePredicateTest.cpp @@ -0,0 +1,35 @@ +#include + +#include +#include +#include + +using geos::geom::CoordinateXY; +using geos::geom::Location; +using geos::triangulate::quadedge::TrianglePredicate; + +namespace tut { + +struct test_trianglepredicate_data {}; + +typedef test_group group; +typedef group::object object; + +group test_trianglepredicate_group("geos::algorithm::TrianglePredicate"); + +// testInCircle() +template<> +template<> +void object::test<1> +() +{ + const CoordinateXY p0 (-1, 0); + const CoordinateXY p1 (0, 1); + const CoordinateXY p2 (1, 0); + + ensure_equals(TrianglePredicate::isInCircleNormalized(p0, p1, p2, CoordinateXY(0, 0)), Location::INTERIOR); + ensure_equals(TrianglePredicate::isInCircleNormalized(p0, p1, p2, CoordinateXY(0, -1)), Location::BOUNDARY); + ensure_equals(TrianglePredicate::isInCircleNormalized(p0, p1, p2, CoordinateXY(2, 0)), Location::EXTERIOR); +} + +}