[geos-commits] [SCM] GEOS branch main updated. 896af228fb27076bdd7f0529d55b3e79e3f10191
git at osgeo.org
git at osgeo.org
Thu Sep 23 14:15:18 PDT 2021
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 896af228fb27076bdd7f0529d55b3e79e3f10191 (commit)
from f2135c8697ba3e7315a6a45871b523fba47d100f (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 896af228fb27076bdd7f0529d55b3e79e3f10191
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Sep 23 14:13:08 2021 -0700
Port https://github.com/locationtech/jts/pull/775, a constrained
delaunay triangulation algorithm.
Adds PolygonTriangulator::triangulate() for a naive ear-clipped
triangulation, and ConstrainedDelaunayTriangulator::triangulate()
for a more aesthetically pleasing polygon triangulation.
Access for the CAPI is via GEOSConstrainedDelaunayTriangulation()
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 13e6548..6257608 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -63,15 +63,6 @@ jobs:
- {
compiler: clang++,
build_type: Release,
- cxxstd: 11,
- arch: 64,
- packages: 'clang',
- cmake: 3.8.*,
- os: ubuntu-16.04
- }
- - {
- compiler: clang++,
- build_type: Release,
cxxstd: 14,
arch: 64,
packages: 'clang',
diff --git a/NEWS b/NEWS
index 264bc79..f73de05 100644
--- a/NEWS
+++ b/NEWS
@@ -15,6 +15,9 @@ Changes in 3.10.0
(Paul Ramsey, Martin Davis)
- CAPI: GEOSWKBWriter_getFlavor, GEOSWKBWriter_setFlavor support
outputting ISO or Extended WKB flavors (#466, Paul Ramsey)
+ - CAPI: GEOSConstrainedDelaunayTriangulation, builds a constrained triangulation
+ of an input Polygon or MultiPolygon, returning a GeometryCollection(Polygon)
+ of the triangles.
- Fixes/Improvements:
- Preserve ordering of lines in overlay results (Martin Davis)
diff --git a/capi/geos_c.cpp b/capi/geos_c.cpp
index 6134356..3c7c8bc 100644
--- a/capi/geos_c.cpp
+++ b/capi/geos_c.cpp
@@ -1625,6 +1625,12 @@ extern "C" {
}
Geometry*
+ GEOSConstrainedDelaunayTriangulation(const Geometry* g)
+ {
+ return GEOSConstrainedDelaunayTriangulation_r(handle, g);
+ }
+
+ Geometry*
GEOSVoronoiDiagram(const Geometry* g, const Geometry* env, double tolerance, int onlyEdges)
{
return GEOSVoronoiDiagram_r(handle, g, env, tolerance, onlyEdges);
diff --git a/capi/geos_c.h.in b/capi/geos_c.h.in
index 0a3cb9d..cdf8469 100644
--- a/capi/geos_c.h.in
+++ b/capi/geos_c.h.in
@@ -941,6 +941,11 @@ extern GEOSGeometry GEOS_DLL * GEOSDelaunayTriangulation_r(
double tolerance,
int onlyEdges);
+/** \see GEOSConstrainedDelaunayTriangulation */
+extern GEOSGeometry GEOS_DLL * GEOSConstrainedDelaunayTriangulation_r(
+ GEOSContextHandle_t handle,
+ const GEOSGeometry *g);
+
/** \see GEOSVoronoiDiagram */
extern GEOSGeometry GEOS_DLL * GEOSVoronoiDiagram_r(
GEOSContextHandle_t extHandle,
@@ -2949,6 +2954,18 @@ extern GEOSGeometry GEOS_DLL * GEOSDelaunayTriangulation(
int onlyEdges);
/**
+* Return a constrained Delaunay triangulation of the vertices of the
+* given polygon(s).
+* For non-polygonal inputs, returns an empty geometry collection.
+*
+* \param g the input geometry whose rings will be used as input
+* \return A newly allocated geometry. NULL on exception.
+* Caller is responsible for freeing with GEOSGeom_destroy().
+*/
+extern GEOSGeometry GEOS_DLL * GEOSConstrainedDelaunayTriangulation(
+ const GEOSGeometry *g);
+
+/**
* Returns the Voronoi polygons of the vertices of the given geometry.
*
* \param g the input geometry whose vertices will be used as sites.
diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp
index f95ca6f..56078eb 100644
--- a/capi/geos_ts_c.cpp
+++ b/capi/geos_ts_c.cpp
@@ -85,6 +85,7 @@
#include <geos/linearref/LengthIndexedLine.h>
#include <geos/triangulate/DelaunayTriangulationBuilder.h>
#include <geos/triangulate/VoronoiDiagramBuilder.h>
+#include <geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h>
#include <geos/util.h>
#include <geos/util/IllegalArgumentException.h>
#include <geos/util/Interrupt.h>
@@ -3707,6 +3708,16 @@ extern "C" {
}
Geometry*
+ GEOSConstrainedDelaunayTriangulation_r(GEOSContextHandle_t extHandle, const Geometry* g1)
+ {
+ using geos::triangulate::polygon::ConstrainedDelaunayTriangulator;
+
+ return execute(extHandle, [&]() -> Geometry* {
+ return ConstrainedDelaunayTriangulator::triangulate(g1).release();
+ });
+ }
+
+ Geometry*
GEOSVoronoiDiagram_r(GEOSContextHandle_t extHandle, const Geometry* g1, const Geometry* env, double tolerance,
int onlyEdges)
{
diff --git a/include/geos/geom/Coordinate.h b/include/geos/geom/Coordinate.h
index e2bde3c..2de7f19 100644
--- a/include/geos/geom/Coordinate.h
+++ b/include/geos/geom/Coordinate.h
@@ -99,6 +99,8 @@ public:
bool equals2D(const Coordinate& other) const;
+ bool equals2D(const Coordinate& other, double tolerance) const;
+
/// 2D only
bool equals(const Coordinate& other) const;
diff --git a/include/geos/geom/Coordinate.inl b/include/geos/geom/Coordinate.inl
index 4389f6f..8af4d3b 100644
--- a/include/geos/geom/Coordinate.inl
+++ b/include/geos/geom/Coordinate.inl
@@ -69,6 +69,18 @@ Coordinate::equals2D(const Coordinate& other) const
}
INLINE bool
+Coordinate::equals2D(const Coordinate& other, double tolerance) const
+{
+ if (std::abs(x - other.x) > tolerance) {
+ return false;
+ }
+ if (std::abs(y - other.y) > tolerance) {
+ return false;
+ }
+ return true;
+}
+
+INLINE bool
Coordinate::equals(const Coordinate& other) const
{
return equals2D(other);
diff --git a/include/geos/geom/CoordinateArraySequence.h b/include/geos/geom/CoordinateArraySequence.h
index 2538098..c4dd1b0 100644
--- a/include/geos/geom/CoordinateArraySequence.h
+++ b/include/geos/geom/CoordinateArraySequence.h
@@ -123,6 +123,8 @@ public:
void expandEnvelope(Envelope& env) const override;
+ void closeRing();
+
std::size_t getDimension() const override;
void apply_rw(const CoordinateFilter* filter) override;
diff --git a/include/geos/geom/Envelope.h b/include/geos/geom/Envelope.h
index ed01723..cd94528 100644
--- a/include/geos/geom/Envelope.h
+++ b/include/geos/geom/Envelope.h
@@ -522,6 +522,11 @@ GEOS_DLL bool operator==(const Envelope& a, const Envelope& b);
GEOS_DLL bool operator!=(const Envelope& a, const Envelope& b);
+/// Strict weak ordering operator for Envelope
+/// This is the C++ equivalent of JTS's compareTo
+GEOS_DLL bool operator< (const Envelope& a, const Envelope& b);
+
+
} // namespace geos::geom
} // namespace geos
diff --git a/include/geos/geom/GeometryFactory.h b/include/geos/geom/GeometryFactory.h
index 23d74c3..3d20f19 100644
--- a/include/geos/geom/GeometryFactory.h
+++ b/include/geos/geom/GeometryFactory.h
@@ -280,6 +280,9 @@ public:
std::unique_ptr<Polygon> createPolygon(std::unique_ptr<LinearRing> && shell,
std::vector<std::unique_ptr<LinearRing>> && holes) const;
+ /// Construct a Polygon from a Coordinate vector, taking ownership of the vector
+ std::unique_ptr<Polygon> createPolygon(std::vector<Coordinate> && coords) const;
+
/// Construct a Polygon with a deep-copy of given arguments
Polygon* createPolygon(const LinearRing& shell,
const std::vector<LinearRing*>& holes) const;
diff --git a/include/geos/geom/Triangle.h b/include/geos/geom/Triangle.h
index 13cec65..cbfbf97 100644
--- a/include/geos/geom/Triangle.h
+++ b/include/geos/geom/Triangle.h
@@ -33,11 +33,10 @@ public:
Coordinate p0, p1, p2;
Triangle(const Coordinate& nP0, const Coordinate& nP1, const Coordinate& nP2)
- :
- p0(nP0),
- p1(nP1),
- p2(nP2)
- {}
+ : p0(nP0)
+ , p1(nP1)
+ , p2(nP2) {}
+
/** \brief
* The inCentre of a triangle is the point which is equidistant
@@ -69,10 +68,71 @@ public:
void circumcentre(Coordinate& resultPoint);
void circumcentreDD(Coordinate& resultPoint);
+ /** Computes the circumcentre of a triangle. */
+ static const Coordinate circumcentre(
+ const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
+
bool isIsoceles();
- /// Computes the circumcentre of a triangle.
- static const Coordinate circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
+ /**
+ * Tests whether a triangle is acute. A triangle is acute if all interior
+ * angles are acute. This is a strict test - right triangles will return
+ * <tt>false</tt>. A triangle which is not acute is either right or obtuse.
+ * <p>
+ * Note: this implementation is not robust for angles very close to 90
+ * degrees.
+ *
+ * @param a a vertex of the triangle
+ * @param b a vertex of the triangle
+ * @param c a vertex of the triangle
+ * @return true if the triangle is acute
+ */
+ static bool isAcute(const Coordinate& a, const Coordinate& b, const Coordinate& c);
+
+ /**
+ * Tests whether a triangle is oriented counter-clockwise.
+ *
+ * @param a a vertex of the triangle
+ * @param b a vertex of the triangle
+ * @param c a vertex of the triangle
+ * @return true if the triangle orientation is counter-clockwise
+ */
+ static bool isCCW(const Coordinate& a, const Coordinate& b, const Coordinate& c);
+
+
+ /**
+ * Tests whether a triangle intersects a point.
+ *
+ * @param a a vertex of the triangle
+ * @param b a vertex of the triangle
+ * @param c a vertex of the triangle
+ * @param p the point to test
+ * @return true if the triangle intersects the point
+ */
+ static bool intersects(const Coordinate& a, const Coordinate& b, const Coordinate& c,
+ const Coordinate& p);
+
+
+ /**
+ * Tests whether a triangle intersects a point.
+ * @param p the point to test
+ * @return true if the triangle intersects the point
+ */
+ bool intersects(const Coordinate& p) { return intersects(p0, p1, p2, p); };
+
+ /**
+ * Tests whether this triangle is oriented counter-clockwise.
+ * @return true if the triangle orientation is counter-clockwise
+ */
+ bool isCCW() { return isCCW(p0, p1, p2); };
+
+ /**
+ * Tests whether this triangle is acute.
+ * @return true if this triangle is acute
+ */
+ bool isAcute() { return isAcute(p0, p1, p2); };
+
+
private:
diff --git a/include/geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h b/include/geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h
new file mode 100644
index 0000000..7d50ace
--- /dev/null
+++ b/include/geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h
@@ -0,0 +1,102 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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
+
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Geometry;
+class GeometryFactory;
+class Polygon;
+}
+namespace triangulate {
+namespace tri {
+class TriList;
+}
+}
+}
+
+using geos::geom::Geometry;
+using geos::geom::GeometryFactory;
+using geos::geom::Polygon;
+using geos::triangulate::tri::TriList;
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/**
+ * Computes the Constrained Delaunay Triangulation of polygons.
+ * The Constrained Delaunay Triangulation of a polygon is a set of triangles
+ * covering the polygon, with the maximum total interior angle over all
+ * possible triangulations. It provides the "best quality" triangulation
+ * of the polygon.
+ * <p>
+ * Holes are supported.
+ */
+class GEOS_DLL ConstrainedDelaunayTriangulator {
+
+private:
+
+ // Members
+ const Geometry* inputGeom;
+ const GeometryFactory* geomFact;
+
+ /**
+ * Computes the triangulation of a single polygon
+ * and returns it as a list of {@link Tri}s.
+ *
+ * @param poly the input polygon
+ * @return list of Tris forming the triangulation
+ */
+ void triangulatePolygon(const Polygon* poly, TriList& triList);
+
+ std::unique_ptr<Geometry> compute();
+
+
+public:
+
+ /**
+ * Constructs a new triangulator.
+ *
+ * @param p_inputGeom the input geometry
+ */
+ ConstrainedDelaunayTriangulator(const Geometry* p_inputGeom)
+ : inputGeom(p_inputGeom)
+ , geomFact(p_inputGeom->getFactory())
+ {}
+
+ /**
+ * Computes the Constrained Delaunay Triangulation of each polygon element in a geometry.
+ *
+ * @param geom the input geometry
+ * @return a GeometryCollection of the computed triangle polygons
+ */
+ static std::unique_ptr<Geometry> triangulate(const Geometry* geom);
+
+
+
+
+};
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/include/geos/triangulate/polygon/PolygonEarClipper.h b/include/geos/triangulate/polygon/PolygonEarClipper.h
new file mode 100644
index 0000000..e6af5c3
--- /dev/null
+++ b/include/geos/triangulate/polygon/PolygonEarClipper.h
@@ -0,0 +1,217 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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/triangulate/polygon/VertexSequencePackedRtree.h>
+
+#include <array>
+#include <memory>
+#include <limits>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+class Polygon;
+class Envelope;
+}
+namespace triangulate {
+namespace tri {
+class TriList;
+}
+}
+}
+
+using geos::geom::Coordinate;
+using geos::geom::Polygon;
+using geos::geom::Envelope;
+using geos::triangulate::tri::TriList;
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/**
+ * Triangulates a polygon using the Ear-Clipping technique.
+ * The polygon is provided as a closed list of contiguous vertices
+ * defining its boundary.
+ * The vertices must have clockwise orientation.
+ *
+ * The polygon boundary must not self-cross,
+ * but may self-touch at points or along an edge.
+ * It may contain repeated points, which are treated as a single vertex.
+ * By default every vertex is triangulated,
+ * including ones which are "flat" (the adjacent segments are collinear).
+ * These can be removed by setting setSkipFlatCorners(boolean)
+ *
+ * The polygon representation does not allow holes.
+ * Polygons with holes can be triangulated by preparing them
+ * with {@link PolygonHoleJoiner}.
+ *
+ * @author Martin Davis
+ *
+ */
+class GEOS_DLL PolygonEarClipper {
+
+private:
+
+ // Members
+
+ static constexpr std::size_t NO_VERTEX_INDEX = std::numeric_limits<std::size_t>::max();
+
+ bool isFlatCornersSkipped = false;
+
+ /**
+ * The polygon vertices are provided in CW orientation.
+ * Thus for convex interior angles
+ * the vertices forming the angle are in CW orientation.
+ */
+ std::vector<Coordinate> vertex;
+ std::vector<std::size_t> vertexNext;
+ std::size_t vertexSize;
+
+ // first available vertex index
+ std::size_t vertexFirst;
+
+ // indices for current corner
+ std::array<std::size_t, 3> cornerIndex;
+
+ /**
+ * Indexing vertices improves ear intersection testing performance a lot.
+ * The polyShell vertices are contiguous, so are suitable for an SPRtree.
+ */
+ VertexSequencePackedRtree vertexCoordIndex;
+
+ // Methods
+
+ std::vector<std::size_t> createNextLinks(std::size_t size) const;
+
+ bool isValidEar(std::size_t cornerIndex, const std::array<Coordinate, 3>& corner);
+
+ /**
+ * Finds another vertex intersecting the corner triangle, if any.
+ * Uses the vertex spatial index for efficiency.
+ *
+ * Also finds any vertex which is a duplicate of the corner apex vertex,
+ * which then requires a full scan of the vertices to confirm ear is valid.
+ * This is usually a rare situation, so has little impact on performance.
+ *
+ * @param cornerIndex the index of the corner apex vertex
+ * @param corner the corner vertices
+ * @return the index of an intersecting or duplicate vertex, or {@link #NO_VERTEX_INDEX} if none
+ */
+ std::size_t findIntersectingVertex(std::size_t cornerIndex, const std::array<Coordinate, 3>& corner) const;
+
+ /**
+ * Scan all vertices in current ring to check if any are duplicates
+ * of the corner apex vertex, and if so whether the corner ear
+ * intersects the adjacent segments and thus is invalid.
+ *
+ * @param cornerIndex the index of the corner apex
+ * @param corner the corner vertices
+ * @return true if the corner ia a valid ear
+ */
+ bool isValidEarScan(std::size_t cornerIndex, const std::array<Coordinate, 3>& corner) const;
+
+ /* private */
+ static Envelope envelope(const std::array<Coordinate, 3>& corner);
+
+ /**
+ * Remove the corner apex vertex and update the candidate corner location.
+ */
+ void removeCorner();
+
+ bool isRemoved(std::size_t vertexIndex) const;
+
+ void initCornerIndex();
+
+ /**
+ * Fetch the corner vertices from the indices.
+ *
+ * @param corner an array for the corner vertices
+ */
+ void fetchCorner(std::array<Coordinate, 3>& cornerVertex) const;
+
+ /**
+ * Move to next corner.
+ */
+ void nextCorner(std::array<Coordinate, 3>& cornerVertex);
+
+ /**
+ * Get the index of the next available shell coordinate starting from the given
+ * index.
+ *
+ * @param index candidate position
+ * @return index of the next available shell coordinate
+ */
+ std::size_t nextIndex(std::size_t index) const;
+
+ bool isConvex(const std::array<Coordinate, 3>& pts) const;
+
+ bool isFlat(const std::array<Coordinate, 3>& pts) const;
+
+ bool hasRepeatedPoint(const std::array<Coordinate, 3>& pts) const;
+
+
+public:
+
+ /**
+ * Creates a new ear-clipper instance.
+ *
+ * @param polyShell the polygon vertices to process
+ */
+ PolygonEarClipper(std::vector<Coordinate>& polyShell);
+
+ /**
+ * Triangulates a polygon via ear-clipping.
+ *
+ * @param polyShell the vertices of the polygon
+ * @param triListResult vector to fill in with the resultant Tri s
+ * @return a list of the Tris
+ */
+ static void triangulate(std::vector<Coordinate>& polyShell, TriList& triListResult);
+
+ /**
+ * Sets whether flat corners formed by collinear adjacent line segments
+ * are included in the triangulation.
+ * Skipping flat corners reduces the number of triangles in the output.
+ * However, it produces a triangulation which does not include
+ * all input vertices. This may be undesirable for downstream processes
+ * (such as computing a Constrained Delaunay Triangulation for
+ * purposes of computing the medial axis).
+ *
+ * The default is to include all vertices in the result triangulation.
+ * This still produces a valid triangulation, with no zero-area triangles.
+ *
+ * Note that repeated vertices are always skipped.
+ *
+ * @param p_isFlatCornersSkipped whether to skip collinear vertices
+ */
+ void setSkipFlatCorners(bool p_isFlatCornersSkipped);
+
+ void compute(TriList& triList);
+
+ std::unique_ptr<Polygon> toGeometry() const;
+
+
+};
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
diff --git a/include/geos/triangulate/polygon/PolygonHoleJoiner.h b/include/geos/triangulate/polygon/PolygonHoleJoiner.h
new file mode 100644
index 0000000..82136ad
--- /dev/null
+++ b/include/geos/triangulate/polygon/PolygonHoleJoiner.h
@@ -0,0 +1,194 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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/geom/Coordinate.h>
+#include <geos/noding/SegmentSetMutualIntersector.h>
+
+#include <unordered_map>
+#include <vector>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Envelope;
+class Geometry;
+class CoordinateSequence;
+class LinearRing;
+}
+namespace triangulate {
+namespace tri {
+class TriList;
+}
+}
+namespace noding {
+class SegmentString;
+}
+}
+
+using geos::geom::Coordinate;
+using geos::geom::CoordinateSequence;
+using geos::geom::Polygon;
+using geos::geom::LinearRing;
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+
+/**
+ * Transforms a polygon with holes into a single self-touching ring
+ * by connecting holes to the exterior shell or to another hole.
+ * The holes are added from the lowest upwards.
+ * As the resulting shell develops, a hole might be added to what was
+ * originally another hole.
+ */
+class GEOS_DLL PolygonHoleJoiner {
+
+private:
+
+ // Members
+
+ static constexpr double EPS = 1.0E-4;
+
+ std::vector<Coordinate> shellCoords;
+
+ // orderedCoords is a copy of shellCoords for sort purposes
+ std::set<Coordinate> orderedCoords;
+
+ // Key: starting end of the cut; Value: list of the other end of the cut
+ std::unordered_map<Coordinate, std::vector<Coordinate>, Coordinate::HashCode> cutMap;
+
+ std::unique_ptr<noding::SegmentSetMutualIntersector> polygonIntersector;
+ const Polygon* inputPolygon;
+
+ // The segstrings allocated in createPolygonIntersector need a
+ // place to hold ownership through the lifecycle of the hole joiner
+ std::vector<std::unique_ptr<noding::SegmentString>> polySegStringStore;
+
+ // Methods
+
+ static std::vector<Coordinate> ringCoordinates(const LinearRing* ring);
+
+ void joinHoles();
+
+ /**
+ * Joins a single hole to the current shellRing.
+ *
+ * @param hole the hole to join
+ */
+ void joinHole(const LinearRing* hole);
+
+ /**
+ * Get the ith shellvertex in shellCoords[] that the current should add after
+ *
+ * @param shellVertex Coordinate of the shell vertex
+ * @param holeVertex Coordinate of the hole vertex
+ * @return the ith shellvertex
+ */
+ std::size_t getShellCoordIndex(const Coordinate& shellVertex, const Coordinate& holeVertex);
+
+ /**
+ * Find the index of the coordinate in ShellCoords ArrayList,
+ * skipping over some number of matches
+ *
+ * @param coord
+ * @return
+ */
+ std::size_t getShellCoordIndexSkip(const Coordinate& coord, std::size_t numSkip);
+
+ /**
+ * Gets a list of shell vertices that could be used to join with the hole.
+ * This list contains only one item if the chosen vertex does not share the same
+ * x value with holeCoord
+ *
+ * @param holeCoord the hole coordinates
+ * @return a list of candidate join vertices
+ */
+ std::vector<Coordinate> getLeftShellVertex(const Coordinate& holeCoord);
+
+ /**
+ * Determine if a line segment between a hole vertex
+ * and a shell vertex lies inside the input polygon.
+ *
+ * @param holeCoord a hole coordinate
+ * @param shellCoord a shell coordinate
+ * @return true if the line lies inside the polygon
+ */
+ bool isJoinable(const Coordinate& holeCoord, const Coordinate& shellCoord) const;
+
+ /**
+ * Tests whether a line segment crosses the polygon boundary.
+ *
+ * @param p0 a vertex
+ * @param p1 a vertex
+ * @return true if the line segment crosses the polygon boundary
+ */
+ bool crossesPolygon(const Coordinate& p0, const Coordinate& p1) const;
+
+ /**
+ * Add hole at proper position in shell coordinate list.
+ * Also adds hole points to ordered coordinates.
+ *
+ * @param shellVertexIndex
+ * @param holeCoords
+ * @param holeVertexIndex
+ */
+ void addHoleToShell(std::size_t shellVertexIndex, const CoordinateSequence* holeCoords, std::size_t holeVertexIndex);
+
+ /**
+ * Sort the holes by minimum X, minimum Y.
+ *
+ * @param poly polygon that contains the holes
+ * @return a list of ordered hole geometry
+ */
+ std::vector<const LinearRing*> sortHoles(const Polygon* poly);
+
+ /**
+ * Gets a list of indices of the leftmost vertices in a ring.
+ *
+ * @param geom the hole ring
+ * @return index of the left most vertex
+ */
+ std::vector<std::size_t> getLeftMostVertex(const LinearRing* ring);
+
+ std::unique_ptr<noding::SegmentSetMutualIntersector> createPolygonIntersector(const Polygon* polygon);
+
+
+public:
+
+ PolygonHoleJoiner(const Polygon* p_inputPolygon);
+
+ static std::vector<Coordinate> join(const Polygon* inputPolygon);
+ static std::unique_ptr<Polygon> joinAsPolygon(const Polygon* inputPolygon);
+
+ /**
+ * Computes the joined ring.
+ *
+ * @return the points in the joined ring
+ */
+ std::vector<Coordinate> compute();
+
+
+};
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/include/geos/triangulate/polygon/PolygonTriangulator.h b/include/geos/triangulate/polygon/PolygonTriangulator.h
new file mode 100644
index 0000000..77271e8
--- /dev/null
+++ b/include/geos/triangulate/polygon/PolygonTriangulator.h
@@ -0,0 +1,107 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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
+
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Geometry;
+class GeometryFactory;
+class Polygon;
+}
+namespace triangulate {
+namespace tri {
+class TriList;
+}
+}
+}
+
+using geos::geom::Geometry;
+using geos::geom::GeometryFactory;
+using geos::geom::Polygon;
+using geos::triangulate::tri::TriList;
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/**
+ * Computes a triangulation of each polygon in a {@link geos::geom::Geometry}.
+ * A polygon triangulation is a non-overlapping set of triangles which
+ * cover the polygon and have the same vertices as the polygon.
+ * The priority is on performance rather than triangulation quality,
+ * so that the output may contain many narrow triangles.
+ *
+ * Holes are handled by joining them to the shell to form a
+ * (self-touching) polygon shell with no holes.
+ * Although invalid, this can be triangulated effectively.
+ *
+ * For better-quality triangulation use geos::triangulate::polygon::ConstrainedDelaunayTriangulator.
+ *
+ * @see geos::triangulate::polygon::ConstrainedDelaunayTriangulator
+ *
+ * @author Martin Davis
+ *
+ */
+class GEOS_DLL PolygonTriangulator {
+
+private:
+
+ // Members
+
+ const Geometry* inputGeom;
+ const GeometryFactory* geomFact;
+
+ std::unique_ptr<Geometry> compute();
+
+ /**
+ * Computes the triangulation of a single polygon
+ *
+ * @return GeometryCollection of triangular polygons
+ */
+ void triangulatePolygon(const Polygon* poly, TriList& triList);
+
+
+public:
+
+ /**
+ * Constructs a new triangulator.
+ *
+ * @param p_inputGeom the input geometry
+ */
+ PolygonTriangulator(const Geometry* p_inputGeom)
+ : inputGeom(p_inputGeom)
+ , geomFact(p_inputGeom->getFactory())
+ {}
+
+ /**
+ * Computes a triangulation of each polygon in a geometry.
+ *
+ * @param geom a geometry containing polygons
+ * @return a GeometryCollection containing the triangle polygons
+ */
+ static std::unique_ptr<Geometry> triangulate(const Geometry* geom);
+
+};
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/include/geos/triangulate/polygon/TriDelaunayImprover.h b/include/geos/triangulate/polygon/TriDelaunayImprover.h
new file mode 100644
index 0000000..eb1e8e0
--- /dev/null
+++ b/include/geos/triangulate/polygon/TriDelaunayImprover.h
@@ -0,0 +1,151 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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/triangulate/tri/TriEdge.h>
+#include <geos/triangulate/tri/Tri.h>
+
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+}
+namespace triangulate {
+namespace tri {
+class TriList;
+}
+}
+}
+
+using geos::geom::Coordinate;
+using geos::triangulate::tri::TriList;
+using geos::triangulate::tri::Tri;
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/**
+ * Improves the quality of a triangulation of Tri via
+ * iterated Delaunay flipping.
+ * This produces the Constrained Delaunay Triangulation
+ * with the constraints being the boundary of the input triangulation.
+ *
+ * @author mdavis
+ *
+ */
+class GEOS_DLL TriDelaunayImprover {
+
+private:
+
+ // Members
+ static constexpr std::size_t MAX_ITERATION = 200;
+ TriList& triList;
+
+ /**
+ * Improves a triangulation by examining pairs of adjacent triangles
+ * (forming a quadrilateral) and testing if flipping the diagonal of
+ * the quadrilateral would produce two new triangles with larger minimum
+ * interior angles.
+ *
+ * @return the number of flips that were made
+ */
+ std::size_t improveScan(TriList& triList);
+
+ /**
+ * Does a flip of the common edge of two Tris if the Delaunay condition is not met.
+ *
+ * @param tri0 a Tri
+ * @param tri1 a Tri
+ * @return true if the triangles were flipped
+ */
+ bool improveNonDelaunay(Tri* tri, TriIndex index);
+
+ /**
+ * Tests if the quadrilateral formed by two adjacent triangles is convex.
+ * opp0-adj0-adj1 and opp1-adj1-adj0 are the triangle corners
+ * and hence are known to be convex.
+ * The quadrilateral is convex if the other corners opp0-adj0-opp1
+ * and opp1-adj1-opp0 have the same orientation (since at least one must be convex).
+ *
+ * @param adj0 adjacent edge vertex 0
+ * @param adj1 adjacent edge vertex 1
+ * @param opp0 corner vertex of triangle 0
+ * @param opp1 corner vertex of triangle 1
+ * @return true if the quadrilateral is convex
+ */
+ static bool isConvex(const Coordinate& adj0, const Coordinate& adj1,
+ const Coordinate& opp0, const Coordinate& opp1);
+
+ /**
+ * Tests if either of a pair of adjacent triangles satisfy the Delaunay condition.
+ * The triangles are opp0-adj0-adj1 and opp1-adj1-adj0.
+ * The Delaunay condition is not met if one opposite vertex
+ * lies is in the circumcircle of the other triangle.
+ *
+ * @param adj0 adjacent edge vertex 0
+ * @param adj1 adjacent edge vertex 1
+ * @param opp0 corner vertex of triangle 0
+ * @param opp1 corner vertex of triangle 1
+ * @return true if the triangles are Delaunay
+ */
+ static bool isDelaunay(const Coordinate& adj0, const Coordinate& adj1,
+ const Coordinate& opp0, const Coordinate& opp1);
+
+ /**
+ * Tests whether a point p is in the circumcircle of a triangle abc
+ * (oriented clockwise).
+ * @param a a vertex of the triangle
+ * @param b a vertex of the triangle
+ * @param c a vertex of the triangle
+ * @param p the point
+ *
+ * @return true if the point is in the circumcircle
+ */
+ static bool
+ isInCircle(const Coordinate& a, const Coordinate& b,
+ const Coordinate& c, const Coordinate& p);
+
+ void improve();
+
+
+public:
+
+ TriDelaunayImprover(TriList& p_triList)
+ : triList(p_triList)
+ {};
+
+ /**
+ * Improves the quality of a triangulation of Tri via
+ * iterated Delaunay flipping.
+ * The Tris are assumed to be linked into a Triangulation
+ * (e.g. via TriangulationBuilder).
+ *
+ * @param triList the list of Tris to flip.
+ */
+ static void improve(TriList& triList);
+
+
+};
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/include/geos/triangulate/polygon/VertexSequencePackedRtree.h b/include/geos/triangulate/polygon/VertexSequencePackedRtree.h
new file mode 100644
index 0000000..616b2ba
--- /dev/null
+++ b/include/geos/triangulate/polygon/VertexSequencePackedRtree.h
@@ -0,0 +1,151 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+class Envelope;
+}
+}
+
+using geos::geom::Coordinate;
+using geos::geom::Envelope;
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/**
+ * A semi-static spatial index for points which occur
+ * in a spatially-coherent sequence.
+ * In particular, this is suitable for indexing the vertices
+ * of a geos::geom::LineString or geos::geom::Polygon ring.
+ *
+ * The index is constructed in a batch fashion on a given sequence of coordinates.
+ * Coordinates can be removed via the remove() method.
+ *
+ * Note that this index queries only the individual points
+ * of the input coordinate sequence,
+ * **not** any line segments which might be lie between them.
+ *
+ * @author Martin Davis
+ *
+ */
+class GEOS_DLL VertexSequencePackedRtree {
+
+private:
+
+
+ /**
+ * Number of items/nodes in a parent node.
+ * Determined empirically. Performance is not too sensitive to this.
+ */
+ static constexpr std::size_t NODE_CAPACITY = 16;
+
+ // Members
+ const std::vector<Coordinate>& items;
+ std::vector<bool> removedItems;
+ std::vector<std::size_t> levelOffset;
+ std::size_t nodeCapacity = NODE_CAPACITY;
+ std::vector<Envelope> bounds;
+
+
+ // Methods
+
+ void build();
+
+ /**
+ * Computes the level offsets.
+ * This is the position in the <tt>bounds</tt> array of each level.
+ *
+ * The levelOffsets array includes a sentinel value of offset[0] = 0.
+ * The top level is always of size 1,
+ * and so also indicates the total number of bounds.
+ *
+ * @return the level offsets
+ */
+ std::vector<std::size_t> computeLevelOffsets();
+
+ static std::size_t ceilDivisor(std::size_t num, std::size_t denom);
+ static std::size_t clampMax(std::size_t x, std::size_t max);
+
+ std::size_t levelNodeCount(std::size_t numNodes);
+
+ std::vector<Envelope> createBounds();
+
+ void fillItemBounds(std::vector<Envelope>& bounds);
+ void fillLevelBounds(std::size_t lvl, std::vector<Envelope>& bounds);
+
+ static Envelope computeNodeEnvelope(const std::vector<Envelope>& bounds,
+ std::size_t start, std::size_t end);
+ static Envelope computeItemEnvelope(const std::vector<Coordinate>& items,
+ std::size_t start, std::size_t end);
+
+ void queryNode(const Envelope& queryEnv,
+ std::size_t level, std::size_t nodeIndex,
+ std::vector<std::size_t>& result) const;
+ void queryNodeRange(const Envelope& queryEnv,
+ std::size_t level, std::size_t nodeStartIndex,
+ std::vector<std::size_t>& result) const;
+ void queryItemRange(const Envelope& queryEnv, std::size_t itemIndex,
+ std::vector<std::size_t>& result) const;
+
+ std::size_t levelSize(std::size_t level) const;
+ bool isNodeEmpty(std::size_t level, std::size_t index);
+ bool isItemsNodeEmpty(std::size_t nodeIndex);
+
+
+public:
+
+ /**
+ * Creates a new tree over the given sequence of coordinates.
+ * The sequence should be spatially coherent to provide query performance.
+ *
+ * @param pts a sequence of points
+ */
+ VertexSequencePackedRtree(const std::vector<Coordinate>& pts);
+
+ std::vector<Envelope> getBounds();
+
+ /**
+ * Removes the input item at the given index from the spatial index.
+ *
+ * @param index the index of the item in the input
+ */
+ void remove(std::size_t index);
+
+ /**
+ * Queries the index to find all items which intersect an extent.
+ * The query result is a list of the indices of input coordinates
+ * which intersect the extent.
+ *
+ * @param queryEnv the query extent
+ * @param result vector to fill with results
+ * @return
+ */
+ void query(const Envelope& queryEnv, std::vector<std::size_t>& result) const;
+
+};
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/include/geos/triangulate/quadedge/TrianglePredicate.h b/include/geos/triangulate/quadedge/TrianglePredicate.h
index a1a7f9f..f380159 100644
--- a/include/geos/triangulate/quadedge/TrianglePredicate.h
+++ b/include/geos/triangulate/quadedge/TrianglePredicate.h
@@ -16,15 +16,21 @@
*
**********************************************************************/
-#ifndef GEOS_TRIANGULATE_QUADEDGE_TRIANGLEPREDICATE_H
-#define GEOS_TRIANGULATE_QUADEDGE_TRIANGLEPREDICATE_H
+#pragma once
#include <geos/export.h>
namespace geos {
-namespace geom { // geos.geom
-
+namespace geom {
class Coordinate;
+}
+}
+
+using geos::geom::Coordinate;
+
+namespace geos {
+namespace triangulate {
+namespace quadedge {
/** \brief
* Algorithms for computing values and predicates
@@ -109,8 +115,9 @@ public:
const Coordinate& p);
} ;
-} // namespace geos.geom
+
+} // namespace geos.triangulate.quadedge
+} // namespace geos.triangulate
} // namespace geos
-#endif //GEOS_TRIANGULATE_QUADEDGE_TRIANGLEPREDICATE_H
diff --git a/include/geos/triangulate/quadedge/Vertex.h b/include/geos/triangulate/quadedge/Vertex.h
index 7be420f..a34f63e 100644
--- a/include/geos/triangulate/quadedge/Vertex.h
+++ b/include/geos/triangulate/quadedge/Vertex.h
@@ -206,7 +206,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 geom::TrianglePredicate::isInCircleRobust(a.p, b.p, c.p, this->p);
+ return triangulate::quadedge::TrianglePredicate::isInCircleRobust(a.p, b.p, c.p, this->p);
}
/**
diff --git a/include/geos/triangulate/tri/Tri.h b/include/geos/triangulate/tri/Tri.h
new file mode 100644
index 0000000..6fe7188
--- /dev/null
+++ b/include/geos/triangulate/tri/Tri.h
@@ -0,0 +1,156 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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/geom/Coordinate.h>
+
+#include <memory>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Geometry;
+class GeometryFactory;
+class Polygon;
+}
+}
+
+using geos::geom::Coordinate;
+using geos::geom::Polygon;
+using geos::geom::GeometryFactory;
+
+typedef int TriIndex;
+
+namespace geos { // geos.
+namespace triangulate { // geos.triangulate
+namespace tri { // geos.triangulate.tri
+
+
+/**
+ * A memory-efficient representation of a triangle in a triangulation.
+ * Contains three vertices, and links to adjacent Tris for each edge.
+ * Tris are constructed independently, and if needed linked
+ * into a triangulation using {@link TriangulationBuilder}.
+ *
+ * @author mdavis
+ *
+ */
+class GEOS_DLL Tri {
+
+private:
+
+ // Members
+ Coordinate p0;
+ Coordinate p1;
+ Coordinate p2;
+
+ /**
+ * triN is the adjacent triangle across the edge pN - pNN.
+ * pNN is the next vertex CW from pN.
+ */
+ Tri* tri0;
+ Tri* tri1;
+ Tri* tri2;
+
+ /**
+ * Replace triOld with triNew
+ *
+ * @param triOld
+ * @param triNew
+ */
+ void replace(Tri* triOld, Tri* triNew);
+
+ /**
+ *
+ * Order: 0: opp0-adj0 edge, 1: opp0-adj1 edge,
+ * 2: opp1-adj0 edge, 3: opp1-adj1 edge
+ *
+ * @param tri
+ * @param index0
+ * @param index1
+ * @return list of adjactent tris
+ */
+ std::vector<Tri*> getAdjacentTris(Tri* tri, TriIndex index0, TriIndex index1);
+
+ void setCoordinates(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
+
+ void flip(Tri* tri, TriIndex index0, TriIndex index1,
+ const Coordinate& adj0, const Coordinate& adj1,
+ const Coordinate& opp0, const Coordinate& opp1);
+
+
+public:
+
+ Tri(const Coordinate& c0, const Coordinate& c1, const Coordinate& c2)
+ : p0(c0)
+ , p1(c1)
+ , p2(c2)
+ , tri0(nullptr)
+ , tri1(nullptr)
+ , tri2(nullptr)
+ {};
+
+ void setAdjacent(Tri* p_tri0, Tri* p_tri1, Tri* p_tri2);
+ void setTri(TriIndex edgeIndex, Tri* tri);
+ void setAdjacent(const Coordinate& pt, Tri* tri);
+
+ /**
+ * Interchanges the vertices of this triangle and a neighbor
+ * so that their common edge
+ * becomes the the other diagonal of the quadrilateral they form.
+ * Neighbour triangle links are modified accordingly.
+ *
+ * @param index the index of the adjacent tri to flip with
+ */
+ void flip(TriIndex index);
+
+ void validate();
+ void validateAdjacent(TriIndex index);
+
+ std::pair<const Coordinate&, const Coordinate&> getEdge(Tri* neighbor) const;
+
+ const Coordinate& getEdgeStart(TriIndex i) const;
+ const Coordinate& getEdgeEnd(TriIndex i) const;
+
+ bool hasCoordinate(const Coordinate& v) const;
+ const Coordinate& getCoordinate(TriIndex i) const;
+
+ TriIndex getIndex(const Coordinate& p) const;
+ TriIndex getIndex(Tri* tri) const;
+
+ Tri* getAdjacent(TriIndex i) const;
+ bool hasAdjacent(TriIndex i) const;
+ bool isAdjacent(Tri* tri) const;
+ int numAdjacent() const;
+
+ static TriIndex next(TriIndex i);
+ static TriIndex prev(TriIndex i);
+ static TriIndex oppVertex(TriIndex edgeIndex);
+ static TriIndex oppEdge(TriIndex vertexIndex);
+ Coordinate midpoint(TriIndex edgeIndex) const;
+
+ std::unique_ptr<Polygon> toPolygon(const GeometryFactory* gf) const;
+
+ friend std::ostream& operator << (std::ostream& os, const Tri&);
+
+};
+
+
+
+
+} // namespace geos.triangulate.tri
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/include/geos/triangulate/tri/TriEdge.h b/include/geos/triangulate/tri/TriEdge.h
new file mode 100644
index 0000000..8058856
--- /dev/null
+++ b/include/geos/triangulate/tri/TriEdge.h
@@ -0,0 +1,80 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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 <iostream>
+
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+}
+}
+
+using geos::geom::Coordinate;
+
+namespace geos { // geos.
+namespace triangulate { // geos.triangulate
+namespace tri { // geos.triangulate.tri
+
+
+/**
+ * Represents an edge in a {@link Tri},
+ * to be used as a key for looking up Tris
+ * while building a triangulation.
+ * The edge value is normalized to allow lookup
+ * of adjacent triangles.
+ *
+ * @author mdavis
+ */
+class GEOS_DLL TriEdge {
+
+private:
+
+ void normalize();
+
+
+public:
+
+ // Members
+ Coordinate p0;
+ Coordinate p1;
+
+ TriEdge(const Coordinate& a, const Coordinate& b)
+ : p0(a)
+ , p1(b)
+ {
+ normalize();
+ }
+
+ struct GEOS_DLL HashCode {
+ std::size_t operator()(const TriEdge& te) const;
+ };
+
+ friend bool operator == (const TriEdge& te0, const TriEdge& te1);
+
+ friend std::ostream& operator << (std::ostream& os, const TriEdge& te);
+
+};
+
+
+
+
+} // namespace geos.triangulate.tri
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/include/geos/triangulate/tri/TriList.h b/include/geos/triangulate/tri/TriList.h
new file mode 100644
index 0000000..f66dccb
--- /dev/null
+++ b/include/geos/triangulate/tri/TriList.h
@@ -0,0 +1,93 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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/triangulate/tri/Tri.h>
+
+#include <geos/export.h>
+#include <iostream>
+#include <deque>
+#include <array>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+class Geometry;
+}
+}
+
+using geos::geom::Coordinate;
+using geos::geom::Geometry;
+
+namespace geos { // geos.
+namespace triangulate { // geos.triangulate
+namespace tri { // geos.triangulate.tri
+
+
+/**
+ * Represents an edge in a {@link Tri},
+ * to be used as a key for looking up Tris
+ * while building a triangulation.
+ * The edge value is normalized to allow lookup
+ * of adjacent triangles.
+ *
+ * @author mdavis
+ */
+class GEOS_DLL TriList {
+
+private:
+
+ // Members
+ std::deque<Tri> triStore;
+ std::vector<Tri*> tris;
+
+ // Methods
+ Tri* create(const Coordinate& c0, const Coordinate& c1, const Coordinate& c2);
+
+
+public:
+
+ TriList() {};
+
+ void add(const Coordinate& c0, const Coordinate& c1, const Coordinate& c2);
+
+ void add(std::array<Coordinate, 3>& corner)
+ {
+ add(corner[0], corner[1], corner[2]);
+ }
+
+ std::unique_ptr<Geometry> toGeometry(const GeometryFactory* geomFact) const;
+
+ friend std::ostream& operator << (std::ostream& os, TriList& te);
+
+ // Support for iterating on TriList
+ typedef std::vector<Tri*>::iterator iterator;
+ typedef std::vector<Tri*>::const_iterator const_iterator;
+ size_t size() const { return tris.size(); }
+ bool empty() const { return tris.empty(); }
+ iterator begin() { return tris.begin(); }
+ iterator end() { return tris.end(); }
+ const_iterator begin() const { return tris.begin(); }
+ const_iterator end() const { return tris.end(); }
+ Tri* operator [] (std::size_t index) { return tris[index]; }
+
+};
+
+
+} // namespace geos.triangulate.tri
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/include/geos/triangulate/tri/TriangulationBuilder.h b/include/geos/triangulate/tri/TriangulationBuilder.h
new file mode 100644
index 0000000..fab1534
--- /dev/null
+++ b/include/geos/triangulate/tri/TriangulationBuilder.h
@@ -0,0 +1,86 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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/triangulate/tri/TriEdge.h>
+
+#include <memory>
+#include <unordered_map>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+}
+namespace triangulate {
+namespace tri {
+class Tri;
+class TriList;
+}
+}
+}
+
+using geos::geom::Coordinate;
+
+namespace geos { // geos.
+namespace triangulate { // geos.triangulate
+namespace tri { // geos.triangulate.tri
+
+
+/**
+ * Builds a triangulation from a set of {@link Tri}s
+ * by linking adjacent Tris.
+ *
+ * @author mdavis
+ *
+ */
+class GEOS_DLL TriangulationBuilder {
+
+private:
+
+ // Members
+ std::unordered_map<TriEdge, Tri*, TriEdge::HashCode> triMap;
+
+ Tri* find(const Coordinate& p0, const Coordinate& p1) const;
+
+ void add(Tri* tri);
+
+ void addAdjacent(Tri* tri, Tri* adj, const Coordinate& p0, const Coordinate& p1);
+
+
+public:
+
+ TriangulationBuilder(TriList& triList);
+
+ /**
+ * Builds the triangulation of a set of {@link Tri}s.
+ *
+ * @param triList the list of Tris
+ */
+ static void build(TriList& triList);
+
+
+};
+
+
+
+
+
+
+
+} // namespace geos.triangulate.tri
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/geom/CoordinateArraySequence.cpp b/src/geom/CoordinateArraySequence.cpp
index 56ef86d..b750f10 100644
--- a/src/geom/CoordinateArraySequence.cpp
+++ b/src/geom/CoordinateArraySequence.cpp
@@ -236,6 +236,14 @@ CoordinateArraySequence::setOrdinate(std::size_t index, std::size_t ordinateInde
}
void
+CoordinateArraySequence::closeRing()
+{
+ if(!isEmpty() && front() != back()) {
+ add(front());
+ }
+}
+
+void
CoordinateArraySequence::apply_rw(const CoordinateFilter* filter)
{
for(auto& coord : vect) {
diff --git a/src/geom/Envelope.cpp b/src/geom/Envelope.cpp
index 17e8fb5..bda4422 100644
--- a/src/geom/Envelope.cpp
+++ b/src/geom/Envelope.cpp
@@ -165,6 +165,40 @@ operator==(const Envelope& a, const Envelope& b)
return a.equals(&b);
}
+bool
+operator< (const Envelope& a, const Envelope& b)
+{
+ /*
+ * Compares two envelopes using lexicographic ordering.
+ * The ordering comparison is based on the usual numerical
+ * comparison between the sequence of ordinates.
+ * Null envelopes are less than all non-null envelopes.
+ */
+ if (a.isNull()) {
+ // null == null
+ if (b.isNull())
+ return false;
+ // null < notnull
+ else
+ return true;
+ }
+ // notnull > null
+ if (b.isNull())
+ return false;
+
+ // compare based on numerical ordering of ordinates
+ if (a.getMinX() < b.getMinX()) return true;
+ if (a.getMinX() > b.getMinX()) return false;
+ if (a.getMinY() < b.getMinY()) return true;
+ if (a.getMinY() > b.getMinY()) return false;
+ if (a.getMaxX() < b.getMaxX()) return true;
+ if (a.getMaxX() > b.getMaxX()) return false;
+ if (a.getMaxY() < b.getMaxY()) return true;
+ if (a.getMaxY() > b.getMaxY()) return false;
+ return false; // == is not strictly <
+}
+
+
/*public*/
size_t
Envelope::hashCode() const
diff --git a/src/geom/GeometryFactory.cpp b/src/geom/GeometryFactory.cpp
index 45f10ea..a6f6c30 100644
--- a/src/geom/GeometryFactory.cpp
+++ b/src/geom/GeometryFactory.cpp
@@ -571,6 +571,19 @@ const
return std::unique_ptr<Polygon>(new Polygon(std::move(shell), *this));
}
+/*public*/
+std::unique_ptr<Polygon>
+GeometryFactory::createPolygon(std::vector<Coordinate> && coords)
+const
+{
+ const geom::CoordinateSequenceFactory* csf = getCoordinateSequenceFactory();
+ std::unique_ptr<geom::CoordinateSequence> cs = csf->create(std::move(coords));
+ std::unique_ptr<geom::LinearRing> lr = createLinearRing(std::move(cs));
+ std::unique_ptr<geom::Polygon> ply = createPolygon(std::move(lr));
+ return ply;
+}
+
+
std::unique_ptr<Polygon>
GeometryFactory::createPolygon(std::unique_ptr<LinearRing> && shell, std::vector<std::unique_ptr<LinearRing>> && holes)
const
diff --git a/src/geom/Triangle.cpp b/src/geom/Triangle.cpp
index 008fc28..70e7f8a 100644
--- a/src/geom/Triangle.cpp
+++ b/src/geom/Triangle.cpp
@@ -15,6 +15,11 @@
#include <geos/geom/Triangle.h>
#include <geos/geom/Coordinate.h>
#include <geos/algorithm/CGAlgorithmsDD.h>
+#include <geos/algorithm/Orientation.h>
+#include <geos/algorithm/Angle.h>
+
+using geos::algorithm::Angle;
+using geos::algorithm::Orientation;
namespace geos {
namespace geom { // geos::geom
@@ -90,5 +95,44 @@ Triangle::det(double m00, double m01, double m10, double m11) const
}
+/* public static */
+bool
+Triangle::isAcute(const Coordinate& a, const Coordinate& b, const Coordinate& c)
+{
+ if (!Angle::isAcute(a, b, c))
+ return false;
+ if (!Angle::isAcute(b, c, a))
+ return false;
+ if (!Angle::isAcute(c, a, b))
+ return false;
+ return true;
+}
+
+
+/* public static */
+bool
+Triangle::isCCW(const Coordinate& a, const Coordinate& b, const Coordinate& c)
+{
+ return Orientation::COUNTERCLOCKWISE == Orientation::index(a, b, c);
+}
+
+
+/* public static */
+bool
+Triangle::intersects(const Coordinate& a, const Coordinate& b, const Coordinate& c, const Coordinate& p)
+{
+ int exteriorIndex = isCCW(a, b, c) ?
+ Orientation::CLOCKWISE : Orientation::COUNTERCLOCKWISE;
+ if (exteriorIndex == Orientation::index(a, b, p))
+ return false;
+ if (exteriorIndex == Orientation::index(b, c, p))
+ return false;
+ if (exteriorIndex == Orientation::index(c, a, p))
+ return false;
+ return true;
+}
+
+
+
} // namespace geos::geom
} // namespace geos
diff --git a/src/triangulate/polygon/ConstrainedDelaunayTriangulator.cpp b/src/triangulate/polygon/ConstrainedDelaunayTriangulator.cpp
new file mode 100644
index 0000000..0d42754
--- /dev/null
+++ b/src/triangulate/polygon/ConstrainedDelaunayTriangulator.cpp
@@ -0,0 +1,95 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+#include <geos/geom/util/PolygonExtracter.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/triangulate/tri/TriList.h>
+#include <geos/triangulate/tri/TriangulationBuilder.h>
+#include <geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h>
+#include <geos/triangulate/polygon/PolygonHoleJoiner.h>
+#include <geos/triangulate/polygon/PolygonEarClipper.h>
+#include <geos/triangulate/polygon/TriDelaunayImprover.h>
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/* public static */
+std::unique_ptr<Geometry>
+ConstrainedDelaunayTriangulator::triangulate(const Geometry* geom)
+{
+ ConstrainedDelaunayTriangulator cdt(geom);
+ return cdt.compute();
+}
+
+
+/* private */
+std::unique_ptr<Geometry>
+ConstrainedDelaunayTriangulator::compute()
+{
+ // short circuit empty input case
+ if(inputGeom->isEmpty()) {
+ auto gf = inputGeom->getFactory();
+ return gf->createGeometryCollection();
+ }
+
+ std::vector<const Polygon*> polys;
+ geom::util::PolygonExtracter::getPolygons(*inputGeom, polys);
+
+ TriList triList;
+ for (const Polygon* poly : polys) {
+ // Skip empty component polygons
+ if (poly->isEmpty())
+ continue;
+ triangulatePolygon(poly, triList);
+ }
+ return triList.toGeometry(geomFact);
+}
+
+
+/**
+* Computes the triangulation of a single polygon
+* and returns it as a list of {@link Tri}s.
+*
+* @param poly the input polygon
+* @return list of Tris forming the triangulation
+*/
+void
+ConstrainedDelaunayTriangulator::triangulatePolygon(const Polygon* poly, TriList& triList)
+{
+ /**
+ * Normalize to ensure that shell and holes have canonical orientation.
+ *
+ * TODO: perhaps better to just correct orientation of rings?
+ */
+ std::unique_ptr<Polygon> polyNorm = poly->clone();
+ polyNorm->normalize();
+
+ std::vector<Coordinate> polyShell = PolygonHoleJoiner::join(polyNorm.get());
+ PolygonEarClipper::triangulate(polyShell, triList);
+ tri::TriangulationBuilder::build(triList);
+ TriDelaunayImprover::improve(triList);
+ return;
+}
+
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/triangulate/polygon/PolygonEarClipper.cpp b/src/triangulate/polygon/PolygonEarClipper.cpp
new file mode 100644
index 0000000..30d6c7e
--- /dev/null
+++ b/src/triangulate/polygon/PolygonEarClipper.cpp
@@ -0,0 +1,368 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/Angle.h>
+#include <geos/algorithm/Orientation.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CoordinateArraySequence.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/Polygon.h>
+#include <geos/geom/Triangle.h>
+#include <geos/triangulate/tri/TriList.h>
+#include <geos/triangulate/polygon/PolygonEarClipper.h>
+#include <geos/util/IllegalStateException.h>
+
+using geos::algorithm::Orientation;
+using geos::algorithm::Angle;
+using geos::triangulate::tri::TriList;
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/* public */
+PolygonEarClipper::PolygonEarClipper(std::vector<Coordinate>& polyShell)
+ : vertex(polyShell)
+ , vertexSize(polyShell.size()-1)
+ , vertexFirst(0)
+ , vertexCoordIndex(polyShell)
+{
+ vertexNext = createNextLinks(vertexSize);
+ initCornerIndex();
+}
+
+
+/* private */
+std::vector<std::size_t>
+PolygonEarClipper::createNextLinks(std::size_t size) const
+{
+ std::vector<std::size_t> next(size);
+ for (std::size_t i = 0; i < size; i++) {
+ next[i] = i + 1;
+ }
+ next[size - 1] = 0;
+ return next;
+}
+
+
+/* public static */
+void
+PolygonEarClipper::triangulate(std::vector<Coordinate>& polyShell, TriList& triListResult)
+{
+ PolygonEarClipper clipper(polyShell);
+ return clipper.compute(triListResult);
+}
+
+
+/* public */
+void
+PolygonEarClipper::setSkipFlatCorners(bool p_isFlatCornersSkipped)
+{
+ isFlatCornersSkipped = p_isFlatCornersSkipped;
+}
+
+
+/* public */
+void
+PolygonEarClipper::compute(TriList& triList)
+{
+ /**
+ * Count scanned corners, to catch infinite loops
+ * (which indicate an algorithm bug)
+ */
+ std::size_t cornerScanCount = 0;
+
+ std::array<Coordinate, 3> corner;
+ fetchCorner(corner);
+
+ /**
+ * Scan continuously around vertex ring,
+ * until all ears have been found.
+ */
+ while (true) {
+ /**
+ * Non-convex corner- remove if flat, or skip
+ * (a concave corner will turn into a convex corner
+ * after enough ears are removed)
+ */
+ if (! isConvex(corner)) {
+ // remove the corner if it is flat or a repeated point
+ bool isCornerRemoved = hasRepeatedPoint(corner)
+ || (isFlatCornersSkipped && isFlat(corner));
+ if (isCornerRemoved) {
+ removeCorner();
+ }
+ cornerScanCount++;
+ if (cornerScanCount > 2 * vertexSize) {
+ throw util::IllegalStateException("Unable to find a convex corner");
+ }
+ }
+ /**
+ * Convex corner - check if it is a valid ear
+ */
+ else if (isValidEar(cornerIndex[1], corner)) {
+ triList.add(corner[0], corner[1], corner[2]);
+ removeCorner();
+ cornerScanCount = 0;
+ }
+ if (cornerScanCount > 2 * vertexSize) {
+ throw util::IllegalStateException("Unable to find a valid ear");
+ }
+
+ //--- done when all corners are processed and removed
+ if (vertexSize < 3) {
+ return;
+ }
+
+ /**
+ * Skip to next corner.
+ * This is done even after an ear is removed,
+ * since that creates fewer skinny triangles.
+ */
+ nextCorner(corner);
+ }
+}
+
+/* private */
+bool
+PolygonEarClipper::isValidEar(std::size_t cornerIdx, const std::array<Coordinate, 3>& corner)
+{
+ std::size_t intApexIndex = findIntersectingVertex(cornerIdx, corner);
+ //--- no intersections found
+ if (intApexIndex == NO_VERTEX_INDEX)
+ return true;
+ //--- check for duplicate corner apex vertex
+ if ( vertex[intApexIndex].equals2D(corner[1]) ) {
+ //--- a duplicate corner vertex requires a full scan
+ return isValidEarScan(cornerIdx, corner);
+ }
+ return false;
+}
+
+
+/* private */
+std::size_t
+PolygonEarClipper::findIntersectingVertex(std::size_t cornerIdx, const std::array<Coordinate, 3>& corner) const
+{
+ Envelope cornerEnv = envelope(corner);
+ std::vector<std::size_t> result;
+ vertexCoordIndex.query(cornerEnv, result);
+
+ std::size_t dupApexIndex = NO_VERTEX_INDEX;
+ //--- check for duplicate vertices
+ for (std::size_t i = 0; i < result.size(); i++) {
+ std::size_t vertIndex = result[i];
+
+ if (vertIndex == cornerIdx ||
+ vertIndex == vertex.size() - 1 ||
+ isRemoved(vertIndex))
+ continue;
+
+ const Coordinate& v = vertex[vertIndex];
+ /**
+ * If another vertex at the corner is found,
+ * need to do a full scan to check the incident segments.
+ * This happens when the polygon ring self-intersects,
+ * usually due to hold joining.
+ * But only report this if no properly intersecting vertex is found,
+ * for efficiency.
+ */
+ if (v.equals2D(corner[1])) {
+ dupApexIndex = vertIndex;
+ }
+ //--- don't need to check other corner vertices
+ else if (v.equals2D(corner[0]) || v.equals2D(corner[2])) {
+ continue;
+ }
+ //--- this is a properly intersecting vertex
+ else if (geom::Triangle::intersects(corner[0], corner[1], corner[2], v))
+ return vertIndex;
+ }
+ if (dupApexIndex != NO_VERTEX_INDEX) {
+ return dupApexIndex;
+ }
+ return NO_VERTEX_INDEX;
+}
+
+
+/* private */
+bool
+PolygonEarClipper::isValidEarScan(std::size_t cornerIdx, const std::array<Coordinate, 3>& corner) const
+{
+ double cornerAngle = Angle::angleBetweenOriented(corner[0], corner[1], corner[2]);
+
+ std::size_t currIndex = nextIndex(vertexFirst);
+ std::size_t prevIndex = vertexFirst;
+ for (std::size_t i = 0; i < vertexSize; i++) {
+ const Coordinate& vPrev = vertex[prevIndex];
+ const Coordinate& v = vertex[currIndex];
+ /**
+ * Because of hole-joining vertices can occur more than once.
+ * If vertex is same as corner[1],
+ * check whether either adjacent edge lies inside the ear corner.
+ * If so the ear is invalid.
+ */
+ if (currIndex != cornerIdx && v.equals2D(corner[1])) {
+ const Coordinate& vNext = vertex[nextIndex(currIndex)];
+
+ //TODO: for robustness use segment orientation instead
+ double aOut = Angle::angleBetweenOriented(corner[0], corner[1], vNext);
+ double aIn = Angle::angleBetweenOriented(corner[0], corner[1], vPrev);
+ if (aOut > 0 && aOut < cornerAngle ) {
+ return false;
+ }
+ if (aIn > 0 && aIn < cornerAngle) {
+ return false;
+ }
+ if (aOut == 0 && aIn == cornerAngle) {
+ return false;
+ }
+ }
+
+ //--- move to next vertex
+ prevIndex = currIndex;
+ currIndex = nextIndex(currIndex);
+ }
+ return true;
+}
+
+
+/* private static */
+Envelope
+PolygonEarClipper::envelope(const std::array<Coordinate, 3>& corner)
+{
+ Envelope cornerEnv(corner[0], corner[1]);
+ cornerEnv.expandToInclude(corner[2]);
+ return cornerEnv;
+}
+
+
+/* private */
+void
+PolygonEarClipper::removeCorner()
+{
+ std::size_t cornerApexIndex = cornerIndex[1];
+ if (vertexFirst == cornerApexIndex) {
+ vertexFirst = vertexNext[cornerApexIndex];
+ }
+ vertexNext[cornerIndex[0]] = vertexNext[cornerApexIndex];
+ vertexCoordIndex.remove(cornerApexIndex);
+ vertexNext[cornerApexIndex] = NO_VERTEX_INDEX;
+ vertexSize--;
+ //-- adjust following corner indexes
+ cornerIndex[1] = nextIndex(cornerIndex[0]);
+ cornerIndex[2] = nextIndex(cornerIndex[1]);
+}
+
+
+/* private */
+bool
+PolygonEarClipper::isRemoved(std::size_t vertexIndex) const
+{
+ return NO_VERTEX_INDEX == vertexNext[vertexIndex];
+}
+
+
+/* private */
+void
+PolygonEarClipper::initCornerIndex()
+{
+ cornerIndex[0] = 0;
+ cornerIndex[1] = 1;
+ cornerIndex[2] = 2;
+}
+
+
+/* private */
+void
+PolygonEarClipper::fetchCorner(std::array<Coordinate, 3>& cornerVertex) const
+{
+ cornerVertex[0] = vertex[cornerIndex[0]];
+ cornerVertex[1] = vertex[cornerIndex[1]];
+ cornerVertex[2] = vertex[cornerIndex[2]];
+}
+
+
+/* private */
+void
+PolygonEarClipper::nextCorner(std::array<Coordinate, 3>& cornerVertex)
+{
+ if ( vertexSize < 3 ) {
+ return;
+ }
+ cornerIndex[0] = nextIndex(cornerIndex[0]);
+ cornerIndex[1] = nextIndex(cornerIndex[0]);
+ cornerIndex[2] = nextIndex(cornerIndex[1]);
+ fetchCorner(cornerVertex);
+}
+
+
+/* private */
+std::size_t
+PolygonEarClipper::nextIndex(std::size_t index) const
+{
+ return vertexNext[index];
+}
+
+
+/* private static */
+bool
+PolygonEarClipper::isConvex(const std::array<Coordinate, 3>& pts) const
+{
+ return Orientation::CLOCKWISE == Orientation::index(pts[0], pts[1], pts[2]);
+}
+
+
+/* private static */
+bool
+PolygonEarClipper::isFlat(const std::array<Coordinate, 3>& pts) const
+{
+ return Orientation::COLLINEAR == Orientation::index(pts[0], pts[1], pts[2]);
+}
+
+
+/* private static */
+bool
+PolygonEarClipper::hasRepeatedPoint(const std::array<Coordinate, 3>& pts) const
+{
+ return pts[1].equals2D(pts[0]) || pts[1].equals2D(pts[2]);
+}
+
+
+/* public */
+std::unique_ptr<Polygon>
+PolygonEarClipper::toGeometry() const
+{
+ auto gf = geom::GeometryFactory::create();
+ std::unique_ptr<geom::CoordinateArraySequence> cs(new geom::CoordinateArraySequence());
+ std::size_t index = vertexFirst;
+ for (std::size_t i = 0; i < vertexSize; i++) {
+ const Coordinate& v = vertex[index];
+ index = nextIndex(index);
+ cs->add(v, true);
+ }
+ cs->closeRing();
+ auto lr = gf->createLinearRing(std::move(cs));
+ return gf->createPolygon(std::move(lr));
+}
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
+
diff --git a/src/triangulate/polygon/PolygonHoleJoiner.cpp b/src/triangulate/polygon/PolygonHoleJoiner.cpp
new file mode 100644
index 0000000..c4afe1d
--- /dev/null
+++ b/src/triangulate/polygon/PolygonHoleJoiner.cpp
@@ -0,0 +1,365 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+#include <geos/algorithm/LineIntersector.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CoordinateSequenceFactory.h>
+#include <geos/geom/LinearRing.h>
+#include <geos/geom/Polygon.h>
+#include <geos/noding/BasicSegmentString.h>
+#include <geos/noding/MCIndexSegmentSetMutualIntersector.h>
+#include <geos/noding/SegmentIntersectionDetector.h>
+#include <geos/noding/SegmentString.h>
+#include <geos/noding/SegmentStringUtil.h>
+#include <geos/util/IllegalStateException.h>
+#include <geos/util/IllegalArgumentException.h>
+
+#include <geos/triangulate/polygon/PolygonHoleJoiner.h>
+
+#include <limits>
+
+
+using geos::geom::GeometryFactory;
+using geos::geom::CoordinateSequence;
+using geos::geom::CoordinateSequenceFactory;
+using geos::noding::SegmentString;
+using geos::noding::SegmentSetMutualIntersector;
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+PolygonHoleJoiner::PolygonHoleJoiner(const Polygon* p_inputPolygon)
+ : polygonIntersector(nullptr)
+ , inputPolygon(p_inputPolygon)
+{
+ polygonIntersector = createPolygonIntersector(p_inputPolygon);
+ if(p_inputPolygon->getNumPoints() < 4)
+ throw util::IllegalArgumentException("Input polygon has too few points");
+}
+
+
+/* public static */
+std::unique_ptr<Polygon>
+PolygonHoleJoiner::joinAsPolygon(const Polygon* inputPolygon)
+{
+ std::vector<Coordinate> coords = join(inputPolygon);
+ const geom::GeometryFactory* gf = inputPolygon->getFactory();
+ return gf->createPolygon(std::move(coords));
+}
+
+
+/* public static */
+std::vector<Coordinate>
+PolygonHoleJoiner::join(const Polygon* inputPolygon)
+{
+ PolygonHoleJoiner joiner(inputPolygon);
+ return joiner.compute();
+}
+
+
+/* public */
+std::vector<Coordinate>
+PolygonHoleJoiner::compute()
+{
+ //--- copy the input polygon shell coords
+ shellCoords = ringCoordinates(inputPolygon->getExteriorRing());
+ if (inputPolygon->getNumInteriorRing() != 0) {
+ joinHoles();
+ }
+ return shellCoords;
+}
+
+
+/* private static */
+std::vector<Coordinate>
+PolygonHoleJoiner::ringCoordinates(const LinearRing* ring)
+{
+ const CoordinateSequence* cs = ring->getCoordinatesRO();
+ std::vector<Coordinate> coords(cs->size());
+ for (std::size_t i = 0; i < cs->size(); i++) {
+ coords[i] = cs->getAt(i);
+ }
+ return coords;
+}
+
+
+/* private */
+void
+PolygonHoleJoiner::joinHoles()
+{
+ orderedCoords.insert(shellCoords.begin(), shellCoords.end());
+ std::vector<const LinearRing*> orderedHoles = sortHoles(inputPolygon);
+ for (std::size_t i = 0; i < orderedHoles.size(); i++) {
+ joinHole(orderedHoles.at(i));
+ }
+}
+
+
+/* private */
+void
+PolygonHoleJoiner::joinHole(const LinearRing* hole)
+{
+ /**
+ * 1) Get a list of HoleVertex Index.
+ * 2) Get a list of ShellVertex.
+ * 3) Get the pair that has the shortest distance between them.
+ * This pair is the endpoints of the cut
+ * 4) The selected ShellVertex may occurs multiple times in
+ * shellCoords[], so find the proper one and add the hole after it.
+ */
+ const CoordinateSequence* holeCoordSeq = hole->getCoordinatesRO();
+ std::vector<std::size_t> holeLeftVerticesIndex = getLeftMostVertex(hole);
+ const Coordinate& holeCoord = holeCoordSeq->getAt(holeLeftVerticesIndex[0]);
+ std::vector<Coordinate> shellCoordsList = getLeftShellVertex(holeCoord);
+ Coordinate shellCoord = shellCoordsList.at(0);
+ std::size_t shortestHoleVertexIndex = 0;
+ //--- pick the shell-hole vertex pair that gives the shortest distance
+ if (std::abs(shellCoord.x - holeCoord.x) < EPS) {
+ double shortest = std::numeric_limits<double>::max();
+ for (std::size_t i = 0; i < holeLeftVerticesIndex.size(); i++) {
+ for (std::size_t j = 0; j < shellCoordsList.size(); j++) {
+ double shellCoordY = shellCoordsList[j].y;
+ double holeLeftY = (holeCoordSeq->getAt(holeLeftVerticesIndex[i])).y;
+ double currLength = std::abs(shellCoordY - holeLeftY);
+ if (currLength < shortest) {
+ shortest = currLength;
+ shortestHoleVertexIndex = i;
+ shellCoord = shellCoordsList[j];
+ }
+ }
+ }
+ }
+ std::size_t shellVertexIndex = getShellCoordIndex(shellCoord,
+ holeCoordSeq->getAt(holeLeftVerticesIndex[shortestHoleVertexIndex]));
+ addHoleToShell(shellVertexIndex, holeCoordSeq, holeLeftVerticesIndex[shortestHoleVertexIndex]);
+}
+
+
+/* private */
+std::size_t
+PolygonHoleJoiner::getShellCoordIndex(const Coordinate& shellVertex, const Coordinate& holeVertex)
+{
+ std::size_t numSkip = 0;
+ std::vector<Coordinate> newValueList;
+ newValueList.emplace_back(holeVertex);
+ auto it = cutMap.find(shellVertex);
+ if (it != cutMap.end()) { // found
+ std::vector<Coordinate>& cutMapCoords = it->second; // value
+ for (const Coordinate& coord : cutMapCoords) {
+ if ( coord.y < holeVertex.y ) {
+ numSkip++;
+ }
+ }
+ cutMapCoords.emplace_back(holeVertex);
+ } else {
+ cutMap[shellVertex] = newValueList;
+ }
+ it = cutMap.find(holeVertex);
+ if (it == cutMap.end()) { // not found
+ cutMap[holeVertex] = newValueList;
+ }
+ return getShellCoordIndexSkip(shellVertex, numSkip);
+}
+
+
+/* private */
+std::size_t
+PolygonHoleJoiner::getShellCoordIndexSkip(const Coordinate& coord, std::size_t numSkip)
+{
+ for (std::size_t i = 0; i < shellCoords.size(); i++) {
+ if (shellCoords[i].equals2D(coord, EPS)) {
+ if (numSkip == 0)
+ return i;
+ numSkip--;
+ }
+ }
+ throw util::IllegalStateException("Vertex is not in shellcoords");
+}
+
+
+/* private */
+std::vector<Coordinate>
+PolygonHoleJoiner::getLeftShellVertex(const Coordinate& holeCoord)
+{
+ std::vector<Coordinate> list;
+ auto it = orderedCoords.upper_bound(holeCoord);
+ // scroll forward until x changes
+ while ((*it).x == holeCoord.x) {
+ it++;
+ }
+ it = orderedCoords.lower_bound(*it);
+ do {
+ it--;
+ } while (!isJoinable(holeCoord, *it) && it != orderedCoords.begin());
+ const Coordinate& closest = *it;
+ list.emplace_back(closest);
+ if ( closest.x != holeCoord.x )
+ return list;
+
+ double chosenX = closest.x;
+ double closestX = closest.x;
+ list.clear();
+ while (chosenX == closestX) {
+ list.emplace_back(closest);
+ it = orderedCoords.lower_bound(closest);
+ // reached the start of the list
+ if (it == orderedCoords.begin())
+ return list;
+ // paper over difference between Java TreeSet.lower()
+ // and std::set::lower_bound()
+ it--;
+ if (it == orderedCoords.begin())
+ return list;
+ closestX = (*it).x;
+ }
+ return list;
+}
+
+
+/* private */
+bool
+PolygonHoleJoiner::isJoinable(const Coordinate& holeCoord, const Coordinate& shellCoord) const
+{
+ /**
+ * Since the line runs between a hole and the shell,
+ * it is inside the polygon if it does not cross the polygon boundary.
+ */
+ bool isJoinable = ! crossesPolygon(holeCoord, shellCoord);
+ /*
+ //--- slow code for testing only
+ LineString join = geomFact.createLineString(new Coordinate[] { holeCoord, shellCoord });
+ bool isJoinableSlow = inputPolygon.covers(join)
+ if (isJoinableSlow != isJoinable) {
+ System.out.println(WKTWriter.toLineString(holeCoord, shellCoord));
+ }
+ //Assert.isTrue(isJoinableSlow == isJoinable);
+ */
+ return isJoinable;
+}
+
+
+/* private */
+bool
+PolygonHoleJoiner::crossesPolygon(const Coordinate& p0, const Coordinate& p1) const
+{
+ // Build up a two-point SegmentString to pass to
+ // the SegmentIntersectionDetector
+ std::vector<Coordinate> coords;
+ coords.emplace_back(p0);
+ coords.emplace_back(p1);
+ auto* csf = inputPolygon->getFactory()->getCoordinateSequenceFactory();
+ auto cs = csf->create(std::move(coords));
+ noding::BasicSegmentString segString(cs.get(), nullptr);
+
+ std::vector<const SegmentString*> segStrings;
+ segStrings.push_back(&segString);
+
+ algorithm::LineIntersector li;
+ noding::SegmentIntersectionDetector segInt(&li);
+ segInt.setFindProper(true);
+ polygonIntersector->setSegmentIntersector(&segInt);
+ polygonIntersector->process(&segStrings);
+ return segInt.hasProperIntersection();
+}
+
+
+/* private */
+void
+PolygonHoleJoiner::addHoleToShell(std::size_t shellVertexIndex,
+ const CoordinateSequence* holeCoords, std::size_t holeVertexIndex)
+{
+ std::vector<Coordinate> newCoords;
+ newCoords.emplace_back(shellCoords[shellVertexIndex]);
+
+ std::size_t nPts = holeCoords->size() - 1;
+ std::size_t i = holeVertexIndex;
+ do {
+ newCoords.emplace_back(holeCoords->getAt(i));
+ i = (i + 1) % nPts;
+ } while (i != holeVertexIndex);
+ newCoords.emplace_back(holeCoords->getAt(holeVertexIndex));
+
+ // Insert newCoords into shellCoords, starting at shellVertextIndex
+ shellCoords.insert(
+ shellCoords.begin() + static_cast<long>(shellVertexIndex),
+ newCoords.begin(),
+ newCoords.end());
+ // Insert all newCoords into orderedCoords
+ orderedCoords.insert(newCoords.begin(), newCoords.end());
+}
+
+
+/* private */
+std::vector<const LinearRing*>
+PolygonHoleJoiner::sortHoles(const Polygon* poly)
+{
+ std::vector<const LinearRing*> holes;
+ for (std::size_t i = 0; i < poly->getNumInteriorRing(); i++) {
+ holes.push_back(poly->getInteriorRingN(i));
+ }
+ // void std::sort( RandomIt first, RandomIt last, Compare comp );
+ std::sort(holes.begin(), holes.end(), [](const LinearRing* a, const LinearRing* b) {
+ return *(a->getEnvelopeInternal()) < *(b->getEnvelopeInternal());
+ });
+
+ return holes;
+}
+
+
+/* private */
+std::vector<std::size_t>
+PolygonHoleJoiner::getLeftMostVertex(const LinearRing* ring)
+{
+ const CoordinateSequence* cs = ring->getCoordinatesRO();
+ std::vector<std::size_t> list;
+ double minX = ring->getEnvelopeInternal()->getMinX();
+ for (std::size_t i = 0; i < cs->size(); i++) {
+ if (std::abs(cs->getAt(i).x - minX) < EPS) {
+ list.push_back(i);
+ }
+ }
+ return list;
+}
+
+/* private */
+std::unique_ptr<SegmentSetMutualIntersector>
+PolygonHoleJoiner::createPolygonIntersector(const Polygon* polygon)
+{
+ std::vector<const SegmentString*> polySegStrings;
+ noding::SegmentStringUtil::extractSegmentStrings(polygon, polySegStrings);
+
+ // Put the ownership on the PolygonHoleJoiner, so that
+ // when the joining is done, the SegmentStrings go away.
+ // TODO: xxxxxx Make extractSegmentStrings() return non-const pointers
+ // so that this cast and others go away
+ for (const SegmentString* ss: polySegStrings)
+ polySegStringStore.emplace_back(const_cast<SegmentString*>(ss));
+
+ std::unique_ptr<SegmentSetMutualIntersector> ssmi(new noding::MCIndexSegmentSetMutualIntersector());
+ ssmi->setBaseSegments(&polySegStrings);
+ return ssmi;
+}
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/triangulate/polygon/PolygonTriangulator.cpp b/src/triangulate/polygon/PolygonTriangulator.cpp
new file mode 100644
index 0000000..731426e
--- /dev/null
+++ b/src/triangulate/polygon/PolygonTriangulator.cpp
@@ -0,0 +1,85 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/util/PolygonExtracter.h>
+#include <geos/triangulate/tri/TriList.h>
+#include <geos/triangulate/polygon/PolygonTriangulator.h>
+#include <geos/triangulate/polygon/PolygonHoleJoiner.h>
+#include <geos/triangulate/polygon/PolygonEarClipper.h>
+
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/* public static */
+std::unique_ptr<Geometry>
+PolygonTriangulator::triangulate(const Geometry* geom)
+{
+ PolygonTriangulator clipper(geom);
+ return clipper.compute();
+}
+
+/* private */
+std::unique_ptr<Geometry>
+PolygonTriangulator::compute()
+{
+ // Short circuit empty case
+ if(inputGeom->isEmpty()) {
+ auto gf = inputGeom->getFactory();
+ return gf->createGeometryCollection();
+ }
+
+ std::vector<const Polygon*> polys;
+ geom::util::PolygonExtracter::getPolygons(*inputGeom, polys);
+
+ TriList triList;
+ for (const Polygon* poly : polys) {
+ // Skip empty component polygons
+ if (poly->isEmpty())
+ continue;
+ triangulatePolygon(poly, triList);
+ }
+ return triList.toGeometry(geomFact);
+}
+
+/* private */
+void
+PolygonTriangulator::triangulatePolygon(const Polygon* poly, TriList& triList)
+{
+ /**
+ * Normalize to ensure that shell and holes have canonical orientation.
+ *
+ * TODO: perhaps better to just correct orientation of rings?
+ */
+ std::unique_ptr<Polygon> polyNorm = poly->clone();
+ polyNorm->normalize();
+
+ std::vector<Coordinate> polyShell = PolygonHoleJoiner::join(polyNorm.get());
+
+ PolygonEarClipper::triangulate(polyShell, triList);
+ //Tri::validate(triList);
+
+ return;
+}
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/triangulate/polygon/TriDelaunayImprover.cpp b/src/triangulate/polygon/TriDelaunayImprover.cpp
new file mode 100644
index 0000000..1e4e588
--- /dev/null
+++ b/src/triangulate/polygon/TriDelaunayImprover.cpp
@@ -0,0 +1,149 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+#include <geos/algorithm/Orientation.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/triangulate/quadedge/TrianglePredicate.h>
+#include <geos/triangulate/tri/Tri.h>
+#include <geos/triangulate/tri/TriList.h>
+#include <geos/triangulate/polygon/TriDelaunayImprover.h>
+
+
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/* public static */
+void
+TriDelaunayImprover::improve(TriList& triList)
+{
+ TriDelaunayImprover improver(triList);
+ improver.improve();
+}
+
+/* private */
+void
+TriDelaunayImprover::improve()
+{
+ for (std::size_t i = 0; i < MAX_ITERATION; i++) {
+ std::size_t improveCount = improveScan(triList);
+ //System.out.println("improve #" + i + " - count = " + improveCount);
+ if (improveCount == 0) {
+ return;
+ }
+ }
+}
+
+/* private */
+std::size_t
+TriDelaunayImprover::improveScan(TriList& tris)
+{
+ std::size_t improveCount = 0;
+ if (tris.size() == 0) return 0; // Fend off infinite loop
+ for (std::size_t i = 0; i < tris.size() - 1; i++) {
+ Tri* tri = tris[i];
+ for (TriIndex j = 0; j < 3; j++) {
+ if (improveNonDelaunay(tri, j)) {
+ improveCount++;
+ }
+ }
+ }
+ return improveCount;
+}
+
+/* private */
+bool
+TriDelaunayImprover::improveNonDelaunay(Tri* tri, TriIndex index)
+{
+ if (tri == nullptr) {
+ return false;
+ }
+ Tri* tri1 = tri->getAdjacent(index);
+ if (tri1 == nullptr) {
+ return false;
+ }
+ //tri0.validate();
+ //tri1.validate();
+
+ TriIndex index1 = tri1->getIndex(tri);
+ const Coordinate& adj0 = tri->getCoordinate(index);
+ const Coordinate& adj1 = tri->getCoordinate(Tri::next(index));
+ const Coordinate& opp0 = tri->getCoordinate(Tri::oppVertex(index));
+ const Coordinate& opp1 = tri1->getCoordinate(Tri::oppVertex(index1));
+
+ /**
+ * The candidate new edge is opp0 - opp1.
+ * Check if it is inside the quadrilateral formed by the two triangles.
+ * This is the case if the quadrilateral is convex.
+ */
+ if (!isConvex(adj0, adj1, opp0, opp1)) {
+ return false;
+ }
+
+ /**
+ * The candidate edge is inside the quadrilateral. Check to see if the flipping
+ * criteria is met. The flipping criteria is to flip if the two triangles are
+ * not Delaunay (i.e. one of the opposite vertices is in the circumcircle of the
+ * other triangle).
+ */
+ if (!isDelaunay(adj0, adj1, opp0, opp1)) {
+ tri->flip(index);
+ return true;
+ }
+ return false;
+}
+
+/* private static */
+bool
+TriDelaunayImprover::isConvex(
+ const Coordinate& adj0, const Coordinate& adj1,
+ const Coordinate& opp0, const Coordinate& opp1)
+{
+ int dir0 = algorithm::Orientation::index(opp0, adj0, opp1);
+ int dir1 = algorithm::Orientation::index(opp1, adj1, opp0);
+ bool isConvex = (dir0 == dir1);
+ return isConvex;
+}
+
+
+/* private static */
+bool
+TriDelaunayImprover::isDelaunay(
+ const Coordinate& adj0, const Coordinate& adj1,
+ const Coordinate& opp0, const Coordinate& opp1)
+{
+ if (isInCircle(adj0, adj1, opp0, opp1)) return false;
+ if (isInCircle(adj1, adj0, opp1, opp0)) return false;
+ return true;
+}
+
+
+/* private static */
+bool
+TriDelaunayImprover::isInCircle(
+ const Coordinate& a, const Coordinate& b,
+ const Coordinate& c, const Coordinate& p)
+{
+ return triangulate::quadedge::TrianglePredicate::isInCircleRobust(a, c, b, p);
+}
+
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/triangulate/polygon/VertexSequencePackedRtree.cpp b/src/triangulate/polygon/VertexSequencePackedRtree.cpp
new file mode 100644
index 0000000..906757f
--- /dev/null
+++ b/src/triangulate/polygon/VertexSequencePackedRtree.cpp
@@ -0,0 +1,327 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Envelope.h>
+#include <geos/triangulate/polygon/VertexSequencePackedRtree.h>
+
+#include <algorithm>
+
+namespace geos {
+namespace triangulate {
+namespace polygon {
+
+
+/**
+* Creates a new tree over the given sequence of coordinates.
+* The sequence should be spatially coherent to provide query performance.
+*
+* @param pts a sequence of points
+*/
+/* public */
+VertexSequencePackedRtree::VertexSequencePackedRtree(const std::vector<Coordinate>& pts)
+ : items(pts)
+ , removedItems(pts.size(), false)
+{
+ build();
+}
+
+/* public */
+std::vector<Envelope>
+VertexSequencePackedRtree::getBounds()
+{
+ std::vector<Envelope> boundsCopy = bounds;
+ // std::vector<Envelope> boundsCopy;
+ // std::copy(bounds.begin(), bounds.end(), boundsCopy.begin());
+ return boundsCopy;
+}
+
+//-- Index Build --------------------------------------------------------------
+
+/* private */
+void
+VertexSequencePackedRtree::build()
+{
+ levelOffset = computeLevelOffsets();
+ bounds = createBounds();
+}
+
+/**
+* Computes the level offsets.
+* This is the position in the <tt>bounds</tt> array of each level.
+*
+* The levelOffsets array includes a sentinel value of offset[0] = 0.
+* The top level is always of size 1,
+* and so also indicates the total number of bounds.
+*
+* @return the level offsets
+*/
+/* private */
+std::vector<std::size_t>
+VertexSequencePackedRtree::computeLevelOffsets()
+{
+ std::vector<std::size_t> offsets;
+ offsets.push_back(0u);
+ std::size_t levelSize = items.size();
+ std::size_t currOffset = 0;
+ do {
+ levelSize = levelNodeCount(levelSize);
+ currOffset += levelSize;
+ offsets.push_back(currOffset);
+ } while (levelSize > 1);
+ return offsets;
+}
+
+/* private static */
+std::size_t
+VertexSequencePackedRtree::ceilDivisor(std::size_t num, std::size_t denom)
+{
+ std::size_t div = num / denom;
+ return div * denom >= num ? div : div + 1;
+}
+
+/* private */
+std::size_t
+VertexSequencePackedRtree::levelNodeCount(std::size_t numNodes)
+{
+ return ceilDivisor(numNodes, nodeCapacity);
+}
+
+/* private */
+std::vector<Envelope>
+VertexSequencePackedRtree::createBounds()
+{
+ std::size_t bndsSize = levelOffset[levelOffset.size() - 1] + 1;
+ std::vector<Envelope> bnds(bndsSize);
+ fillItemBounds(bnds);
+
+ for (std::size_t lvl = 1; lvl < levelOffset.size(); lvl++) {
+ fillLevelBounds(lvl, bnds);
+ }
+ return bnds;
+}
+
+
+/* private */
+void
+VertexSequencePackedRtree::fillItemBounds(std::vector<Envelope>& bnds)
+{
+ std::size_t nodeStart = 0;
+ std::size_t bndIndex = 0;
+ do {
+ std::size_t nodeEnd = clampMax(nodeStart + nodeCapacity, items.size());
+ bnds[bndIndex++] = computeItemEnvelope(items, nodeStart, nodeEnd);
+ nodeStart = nodeEnd;
+ }
+ while (nodeStart < items.size());
+}
+
+
+/* private */
+void
+VertexSequencePackedRtree::fillLevelBounds(std::size_t lvl, std::vector<Envelope>& bnds)
+{
+ std::size_t levelStart = levelOffset[lvl - 1];
+ std::size_t levelEnd = levelOffset[lvl];
+ std::size_t nodeStart = levelStart;
+ std::size_t levelBndIndex = levelOffset[lvl];
+ do {
+ std::size_t nodeEnd = clampMax(nodeStart + nodeCapacity, levelEnd);
+ bnds[levelBndIndex++] = computeNodeEnvelope(bnds, nodeStart, nodeEnd);
+ nodeStart = nodeEnd;
+ }
+ while (nodeStart < levelEnd);
+}
+
+
+/* private static */
+Envelope
+VertexSequencePackedRtree::computeNodeEnvelope(const std::vector<Envelope>& bnds,
+ std::size_t start, std::size_t end)
+{
+ Envelope env;
+ for (std::size_t i = start; i < end; i++) {
+ env.expandToInclude(bnds[i]);
+ }
+ return env;
+}
+
+/* private static */
+Envelope
+VertexSequencePackedRtree::computeItemEnvelope(const std::vector<Coordinate>& items,
+ std::size_t start, std::size_t end)
+{
+ Envelope env;
+ for (std::size_t i = start; i < end; i++) {
+ env.expandToInclude(items[i]);
+ }
+ return env;
+}
+
+//-- Index Query --------------------------------------------------------------
+
+/**
+* Queries the index to find all items which intersect an extent.
+* The query result is a list of the indices of input coordinates
+* which intersect the extent.
+*
+* @param queryEnv the query extent
+* @return an array of the indices of the input coordinates
+*/
+/* public */
+void
+VertexSequencePackedRtree::query(const Envelope& queryEnv, std::vector<std::size_t>& result) const
+{
+ std::size_t level = levelOffset.size() - 1;
+ queryNode(queryEnv, level, 0, result);
+}
+
+
+/* private */
+void
+VertexSequencePackedRtree::queryNode(const Envelope& queryEnv,
+ std::size_t level, std::size_t nodeIndex,
+ std::vector<std::size_t>& result) const
+{
+ std::size_t boundsIndex = levelOffset[level] + nodeIndex;
+ Envelope nodeEnv = bounds[boundsIndex];
+
+ //--- node is empty
+ if (nodeEnv.isNull())
+ return;
+ if (! queryEnv.intersects(nodeEnv))
+ return;
+
+ std::size_t childNodeIndex = nodeIndex * nodeCapacity;
+ if (level == 0) {
+ queryItemRange(queryEnv, childNodeIndex, result);
+ }
+ else {
+ queryNodeRange(queryEnv, level - 1, childNodeIndex, result);
+ }
+}
+
+/* private */
+void
+VertexSequencePackedRtree::queryNodeRange(const Envelope& queryEnv,
+ std::size_t level, std::size_t nodeStartIndex,
+ std::vector<std::size_t>& result) const
+{
+ std::size_t levelMax = levelSize(level);
+ for (std::size_t i = 0; i < nodeCapacity; i++) {
+ std::size_t index = nodeStartIndex + i;
+ if (index >= levelMax)
+ return;
+ queryNode(queryEnv, level, index, result);
+ }
+}
+
+
+/* private */
+std::size_t
+VertexSequencePackedRtree::levelSize(std::size_t level) const
+{
+ return levelOffset[level + 1] - levelOffset[level];
+}
+
+
+/* private */
+void
+VertexSequencePackedRtree::queryItemRange(const Envelope& queryEnv,
+ std::size_t itemIndex, std::vector<std::size_t>& result) const
+{
+ for (std::size_t i = 0; i < nodeCapacity; i++) {
+ std::size_t index = itemIndex + i;
+ if (index >= items.size())
+ return;
+ const Coordinate& p = items[index];
+ bool removed = removedItems[index];
+ if ( (!removed) && queryEnv.contains(p))
+ result.push_back(index);
+ }
+}
+
+//-- Index Modify --------------------------------------------------------------
+
+/**
+* Removes the input item at the given index from the spatial index.
+*
+* @param index the index of the item in the input
+*/
+/* public */
+void
+VertexSequencePackedRtree::remove(std::size_t index)
+{
+ removedItems[index] = true;
+
+ //--- prune the item parent node if all its items are removed
+ std::size_t nodeIndex = index / nodeCapacity;
+ if (! isItemsNodeEmpty(nodeIndex))
+ return;
+
+ bounds[nodeIndex].setToNull();
+
+ if (levelOffset.size() <= 2)
+ return;
+
+ //-- prune the node parent if all children removed
+ std::size_t nodeLevelIndex = nodeIndex / nodeCapacity;
+ if (! isNodeEmpty(1, nodeLevelIndex))
+ return;
+
+ std::size_t nodeIndex1 = levelOffset[1] + nodeLevelIndex;
+ bounds[nodeIndex1].setToNull();
+
+ // TODO: propagate removal up the tree nodes?
+}
+
+/* private */
+bool
+VertexSequencePackedRtree::isNodeEmpty(std::size_t level, std::size_t index)
+{
+ std::size_t start = index * nodeCapacity;
+ std::size_t end = clampMax(start + nodeCapacity, levelOffset[level]);
+ for (std::size_t i = start; i < end; i++) {
+ if (!bounds[i].isNull())
+ return false;
+ }
+ return true;
+}
+
+/* private */
+bool
+VertexSequencePackedRtree::isItemsNodeEmpty(std::size_t nodeIndex)
+{
+ std::size_t start = nodeIndex * nodeCapacity;
+ std::size_t end = clampMax(start + nodeCapacity, items.size());
+ for (std::size_t i = start; i < end; i++) {
+ if (!removedItems[i]) return false;
+ }
+ return true;
+}
+
+/* private static */
+std::size_t
+VertexSequencePackedRtree::clampMax(std::size_t x, std::size_t max)
+{
+ return x > max ? max: x;
+}
+
+
+} // namespace geos.triangulate.polygon
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/triangulate/quadedge/TrianglePredicate.cpp b/src/triangulate/quadedge/TrianglePredicate.cpp
index 4584b61..4c8da95 100644
--- a/src/triangulate/quadedge/TrianglePredicate.cpp
+++ b/src/triangulate/quadedge/TrianglePredicate.cpp
@@ -20,8 +20,11 @@
#include <geos/geom/Coordinate.h>
+using geos::geom::Coordinate;
+
namespace geos {
-namespace geom { // geos.geom
+namespace triangulate {
+namespace quadedge {
bool
TrianglePredicate::isInCircleNonRobust(
@@ -85,5 +88,6 @@ TrianglePredicate::isInCircleRobust(
return isInCircleNormalized(a, b, c, p);
}
-} // namespace geos.geom
+} // namespace geos.triangulate.quadedge
+} // namespace geos.triangulate
} // namespace geos
diff --git a/src/triangulate/tri/Tri.cpp b/src/triangulate/tri/Tri.cpp
new file mode 100644
index 0000000..696be7e
--- /dev/null
+++ b/src/triangulate/tri/Tri.cpp
@@ -0,0 +1,406 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/Orientation.h>
+#include <geos/algorithm/LineIntersector.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CoordinateSequenceFactory.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LinearRing.h>
+#include <geos/geom/Polygon.h>
+#include <geos/triangulate/tri/Tri.h>
+#include <geos/util/IllegalArgumentException.h>
+
+
+using geos::util::IllegalArgumentException;
+using geos::algorithm::Orientation;
+
+namespace geos { // geos
+namespace triangulate { // geos.triangulate
+namespace tri { // geos.triangulate.tri
+
+
+/* public */
+void
+Tri::setAdjacent(Tri* p_tri0, Tri* p_tri1, Tri* p_tri2)
+{
+ tri0 = p_tri0;
+ tri1 = p_tri1;
+ tri2 = p_tri2;
+}
+
+/* public */
+void
+Tri::setTri(TriIndex edgeIndex, Tri* tri)
+{
+ switch (edgeIndex) {
+ case 0: tri0 = tri; return;
+ case 1: tri1 = tri; return;
+ case 2: tri2 = tri; return;
+ }
+ assert(false); // never reach here
+}
+
+/* private */
+void
+Tri::setCoordinates(const Coordinate& np0, const Coordinate& np1, const Coordinate& np2)
+{
+ p0 = np0;
+ p1 = np1;
+ p2 = np2;
+}
+
+/* public */
+void
+Tri::setAdjacent(const Coordinate& pt, Tri* tri)
+{
+ TriIndex index = getIndex(pt);
+ setTri(index, tri);
+ // TODO: validate that tri is adjacent at the edge specified
+}
+
+/**
+* Replace triOld with triNew
+*
+* @param triOld
+* @param triNew
+*/
+/* private */
+void
+Tri::replace(Tri* triOld, Tri* triNew)
+{
+ if ( tri0 != nullptr && tri0 == triOld ) {
+ tri0 = triNew;
+ }
+ else if ( tri1 != nullptr && tri1 == triOld ) {
+ tri1 = triNew;
+ }
+ else if ( tri2 != nullptr && tri2 == triOld ) {
+ tri2 = triNew;
+ }
+}
+
+/* public */
+void
+Tri::flip(TriIndex index)
+{
+ Tri* tri = getAdjacent(index);
+ TriIndex index1 = tri->getIndex(this);
+
+ Coordinate adj0 = getCoordinate(index);
+ Coordinate adj1 = getCoordinate(next(index));
+ Coordinate opp0 = getCoordinate(oppVertex(index));
+ Coordinate opp1 = tri->getCoordinate(oppVertex(index1));
+
+ flip(tri, index, index1, adj0, adj1, opp0, opp1);
+}
+
+/* private */
+void
+Tri::flip(Tri* tri, TriIndex index0, TriIndex index1,
+ const Coordinate& adj0, const Coordinate& adj1,
+ const Coordinate& opp0, const Coordinate& opp1)
+{
+ //System.out.println("Flipping: " + this + " -> " + tri);
+ //validate();
+ //tri.validate();
+
+ setCoordinates(opp1, opp0, adj0);
+ tri->setCoordinates(opp0, opp1, adj1);
+
+ /**
+ * Order: 0: opp0-adj0 edge, 1: opp0-adj1 edge,
+ * 2: opp1-adj0 edge, 3: opp1-adj1 edge
+ */
+ std::vector<Tri*> adjacent = getAdjacentTris(tri, index0, index1);
+ setAdjacent(tri, adjacent[0], adjacent[2]);
+ //--- update the adjacent triangles with new adjacency
+ if ( adjacent[2] != nullptr ) {
+ adjacent[2]->replace(tri, this);
+ }
+ tri->setAdjacent(this, adjacent[3], adjacent[1]);
+ if ( adjacent[1] != nullptr ) {
+ adjacent[1]->replace(this, tri);
+ }
+ //validate();
+ //tri.validate();
+}
+
+/**
+*
+* Order: 0: opp0-adj0 edge, 1: opp0-adj1 edge,
+* 2: opp1-adj0 edge, 3: opp1-adj1 edge
+*
+* @param tri
+* @param index0
+* @param index1
+* @return
+*/
+/* private */
+std::vector<Tri*>
+Tri::getAdjacentTris(Tri* triAdj, TriIndex index, TriIndex indexAdj)
+{
+ std::vector<Tri*> adj(4);
+ adj[0] = getAdjacent(prev(index));
+ adj[1] = getAdjacent(next(index));
+ adj[2] = triAdj->getAdjacent(next(indexAdj));
+ adj[3] = triAdj->getAdjacent(prev(indexAdj));
+ return adj;
+}
+
+/* public */
+void
+Tri::validate()
+{
+ if ( Orientation::CLOCKWISE != Orientation::index(p0, p1, p2) ) {
+ throw IllegalArgumentException("Tri is not oriented correctly");
+ }
+
+ validateAdjacent(0);
+ validateAdjacent(1);
+ validateAdjacent(2);
+}
+
+/* public */
+void
+Tri::validateAdjacent(TriIndex index)
+{
+ Tri* tri = getAdjacent(index);
+ if (tri == nullptr) return;
+
+ assert(isAdjacent(tri));
+ assert(tri->isAdjacent(this));
+
+ // const Coordinate& e0 = getCoordinate(index);
+ // const Coordinate& e1 = getCoordinate(next(index));
+ // TriIndex indexNeighbor = tri->getIndex(this);
+ // const Coordinate& n0 = tri->getCoordinate(indexNeighbor);
+ // const Coordinate& n1 = tri->getCoordinate(next(indexNeighbor));
+ // assert(e0.equals2D(n1)); // Edge coord not equal
+ // assert(e1.equals2D(n0)); // Edge coord not equal
+
+ //--- check that no edges cross
+ algorithm::LineIntersector li;
+ for (TriIndex i = 0; i < 3; i++) {
+ for (TriIndex j = 0; j < 3; j++) {
+ const Coordinate& p00 = getCoordinate(i);
+ const Coordinate& p01 = getCoordinate(next(i));
+ const Coordinate& p10 = tri->getCoordinate(j);
+ const Coordinate& p11 = tri->getCoordinate(next(j));
+ li.computeIntersection(p00, p01, p10, p11);
+ assert(!li.isProper());
+ }
+ }
+}
+
+/* public */
+std::pair<const Coordinate&, const Coordinate&>
+Tri::getEdge(Tri* neighbor) const
+{
+ TriIndex index = getIndex(neighbor);
+ TriIndex nextTri = next(index);
+
+ // const Coordinate& e0 = getCoordinate(index);
+ // const Coordinate& e1 = getCoordinate(nextTri);
+ // assert (neighbor->hasCoordinate(e0));
+ // assert (neighbor->hasCoordinate(e1));
+ // TriIndex iN = neighbor->getIndex(e0);
+ // TriIndex iNPrev = prev(iN);
+ // assert (neighbor->getIndex(e1) == iNPrev);
+
+ std::pair<const Coordinate&, const Coordinate&> edge(getCoordinate(index), getCoordinate(nextTri));
+ return edge;
+}
+
+/* public */
+const Coordinate&
+Tri::getEdgeStart(TriIndex i) const
+{
+ return getCoordinate(i);
+}
+
+/* public */
+const Coordinate&
+Tri::getEdgeEnd(TriIndex i) const
+{
+ return getCoordinate(next(i));
+}
+
+/* public */
+bool
+Tri::hasCoordinate(const Coordinate& v) const
+{
+ if ( p0.equals(v) || p1.equals(v) || p2.equals(v) ) {
+ return true;
+ }
+ return false;
+}
+
+/* public */
+const Coordinate&
+Tri::getCoordinate(TriIndex i) const
+{
+ if ( i == 0 ) {
+ return p0;
+ }
+ if ( i == 1 ) {
+ return p1;
+ }
+ return p2;
+}
+
+/* public */
+TriIndex
+Tri::getIndex(const Coordinate& p) const
+{
+ if ( p0.equals2D(p) )
+ return 0;
+ if ( p1.equals2D(p) )
+ return 1;
+ if ( p2.equals2D(p) )
+ return 2;
+ return -1;
+}
+
+/* public */
+TriIndex
+Tri::getIndex(Tri* tri) const
+{
+ if ( tri0 == tri )
+ return 0;
+ if ( tri1 == tri )
+ return 1;
+ if ( tri2 == tri )
+ return 2;
+ return -1;
+}
+
+/* public */
+Tri*
+Tri::getAdjacent(TriIndex i) const
+{
+ switch(i) {
+ case 0: return tri0;
+ case 1: return tri1;
+ case 2: return tri2;
+ }
+ assert(false); // Never get here
+ return nullptr;
+}
+
+/* public */
+bool
+Tri::hasAdjacent(TriIndex i) const
+{
+ return nullptr != getAdjacent(i);
+}
+
+/* public */
+bool
+Tri::isAdjacent(Tri* tri) const
+{
+ return getIndex(tri) >= 0;
+}
+
+/* public */
+int
+Tri::numAdjacent() const
+{
+ int num = 0;
+ if ( tri0 != nullptr )
+ num++;
+ if ( tri1 != nullptr )
+ num++;
+ if ( tri2 != nullptr )
+ num++;
+ return num;
+}
+
+/* public static */
+TriIndex
+Tri::next(TriIndex i)
+{
+ switch (i) {
+ case 0: return 1;
+ case 1: return 2;
+ case 2: return 0;
+ }
+ return -1;
+}
+
+/* public static */
+TriIndex
+Tri::prev(TriIndex i)
+{
+ switch (i) {
+ case 0: return 2;
+ case 1: return 0;
+ case 2: return 1;
+ }
+ return -1;
+}
+
+/* public static */
+TriIndex
+Tri::oppVertex(TriIndex edgeIndex)
+{
+ return prev(edgeIndex);
+}
+
+/* public static */
+TriIndex
+Tri::oppEdge(TriIndex vertexIndex)
+{
+ return next(vertexIndex);
+}
+
+/* public */
+Coordinate
+Tri::midpoint(TriIndex edgeIndex) const
+{
+ const Coordinate& np0 = getCoordinate(edgeIndex);
+ const Coordinate& np1 = getCoordinate(next(edgeIndex));
+ double midX = (np0.x + np1.x) / 2;
+ double midY = (np0.y + np1.y) / 2;
+ return Coordinate(midX, midY);
+}
+
+/* public */
+std::unique_ptr<geom::Polygon>
+Tri::toPolygon(const geom::GeometryFactory* gf) const
+{
+ std::vector<Coordinate> coords(4);
+ coords[0] = p0; coords[1] = p1;
+ coords[2] = p2; coords[3] = p0;
+
+ return gf->createPolygon(std::move(coords));
+}
+
+std::ostream&
+operator<<(std::ostream& os, const Tri& tri)
+{
+ os << "POLYGON ((";
+ os << tri.p0 << ", ";
+ os << tri.p1 << ", ";
+ os << tri.p2 << ", ";
+ os << tri.p0 << "))";
+ return os;
+}
+
+
+} // namespace geos.triangulate.tri
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/triangulate/tri/TriEdge.cpp b/src/triangulate/tri/TriEdge.cpp
new file mode 100644
index 0000000..5d58674
--- /dev/null
+++ b/src/triangulate/tri/TriEdge.cpp
@@ -0,0 +1,72 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/geom/Coordinate.h>
+#include <geos/triangulate/tri/TriEdge.h>
+#include <geos/export.h>
+
+#include <iostream>
+
+namespace geos { // geos
+namespace triangulate { // geos.triangulate
+namespace tri { // geos.triangulate.tri
+
+
+/* private */
+void
+TriEdge::normalize()
+{
+ if ( p0.compareTo(p1) < 0 ) {
+ Coordinate tmp = p0;
+ p0 = p1;
+ p1 = tmp;
+ }
+}
+
+
+/* public */
+std::size_t
+TriEdge::HashCode::operator()(TriEdge const& te) const
+{
+ geom::Coordinate::HashCode coordHash;
+
+ std::size_t result = 17;
+ result ^= coordHash(te.p0);
+ result ^= coordHash(te.p1);
+ return result;
+}
+
+/* public */
+bool
+operator == (TriEdge const& te0, TriEdge const& te1)
+{
+ return te0.p0 == te1.p0 && te0.p1 == te1.p1;
+}
+
+
+
+std::ostream&
+operator<<(std::ostream& os, const TriEdge& te)
+{
+ os << "LINESTRING (";
+ os << te.p0 << ", ";
+ os << te.p1 << ")";
+ return os;
+}
+
+
+} // namespace geos.triangulate.tri
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/triangulate/tri/TriList.cpp b/src/triangulate/tri/TriList.cpp
new file mode 100644
index 0000000..40f4081
--- /dev/null
+++ b/src/triangulate/tri/TriList.cpp
@@ -0,0 +1,75 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/Polygon.h>
+#include <geos/triangulate/tri/TriList.h>
+#include <geos/export.h>
+
+
+namespace geos { // geos
+namespace triangulate { // geos.triangulate
+namespace tri { // geos.triangulate.tri
+
+
+/* private */
+Tri*
+TriList::create(const Coordinate& c0, const Coordinate& c1, const Coordinate& c2)
+{
+ triStore.emplace_back(c0, c1, c2);
+ Tri* newTri = &triStore.back();
+ return newTri;
+}
+
+
+/* public */
+void
+TriList::add(const Coordinate& c0, const Coordinate& c1, const Coordinate& c2)
+{
+ Tri* newTri = create(c0, c1, c2);
+ tris.push_back(newTri);
+}
+
+/* public */
+std::unique_ptr<Geometry>
+TriList::toGeometry(const geom::GeometryFactory* geomFact) const
+{
+ std::vector<std::unique_ptr<Geometry>> geoms;
+ for (auto* tri: tris) {
+ std::unique_ptr<Geometry> geom = tri->toPolygon(geomFact);
+ geoms.emplace_back(geom.release());
+ }
+ return geomFact->createGeometryCollection(std::move(geoms));
+}
+
+
+std::ostream&
+operator<<(std::ostream& os, TriList& triList)
+{
+ os << "TRILIST ";
+ os << "[" << triList.size() << "] (";
+ for (auto* tri: triList) {
+ os << " " << *tri << "," << std::endl;
+ }
+ os << ")";
+ return os;
+}
+
+
+} // namespace geos.triangulate.tri
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/src/triangulate/tri/TriangulationBuilder.cpp b/src/triangulate/tri/TriangulationBuilder.cpp
new file mode 100644
index 0000000..405a640
--- /dev/null
+++ b/src/triangulate/tri/TriangulationBuilder.cpp
@@ -0,0 +1,103 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************/
+
+
+#include <geos/geom/Coordinate.h>
+#include <geos/triangulate/tri/Tri.h>
+#include <geos/triangulate/tri/TriList.h>
+#include <geos/triangulate/tri/TriangulationBuilder.h>
+
+
+namespace geos { // geos
+namespace triangulate { // geos.triangulate
+namespace tri { // geos.triangulate.tri
+
+
+/* public static */
+void
+TriangulationBuilder::build(TriList& triList)
+{
+ TriangulationBuilder tb(triList);
+}
+
+
+/* private */
+TriangulationBuilder::TriangulationBuilder(TriList& triList)
+{
+ for (Tri* tri : triList) {
+ add(tri);
+ }
+}
+
+/* private */
+Tri*
+TriangulationBuilder::find(const Coordinate& p0, const Coordinate& p1) const
+{
+ const TriEdge e(p0, p1);
+ auto it = triMap.find(e);
+ if (it == triMap.end()) {
+ // not found
+ return nullptr;
+ }
+ else {
+ // found
+ return it->second;
+ }
+}
+
+/* private */
+void
+TriangulationBuilder::add(Tri* tri)
+{
+ const Coordinate& p0 = tri->getCoordinate(0);
+ const Coordinate& p1 = tri->getCoordinate(1);
+ const Coordinate& p2 = tri->getCoordinate(2);
+
+ // get adjacent triangles, if any
+ Tri* n0 = find(p0, p1);
+ Tri* n1 = find(p1, p2);
+ Tri* n2 = find(p2, p0);
+
+ tri->setAdjacent(n0, n1, n2);
+ addAdjacent(tri, n0, p0, p1);
+ addAdjacent(tri, n1, p1, p2);
+ addAdjacent(tri, n2, p2, p0);
+}
+
+/* private */
+void
+TriangulationBuilder::addAdjacent(Tri* tri, Tri* adj, const Coordinate& p0, const Coordinate& p1)
+{
+ /**
+ * If adjacent is null, this tri is first one to be recorded for edge
+ */
+ if (adj == nullptr) {
+
+ triMap.emplace(
+ std::piecewise_construct,
+ std::forward_as_tuple(p0, p1),
+ std::forward_as_tuple(tri));
+
+ // TriEdge e(p0, p1);
+ // triMap.insert(std::make_pair(e, tri));
+ return;
+ }
+ adj->setAdjacent(p1, tri);
+}
+
+
+} // namespace geos.triangulate.tri
+} // namespace geos.triangulate
+} // namespace geos
+
diff --git a/tests/unit/capi/GEOSConstrainedDelaunayTriangulationTest.cpp b/tests/unit/capi/GEOSConstrainedDelaunayTriangulationTest.cpp
new file mode 100644
index 0000000..e428134
--- /dev/null
+++ b/tests/unit/capi/GEOSConstrainedDelaunayTriangulationTest.cpp
@@ -0,0 +1,78 @@
+//
+// Test Suite for C-API GEOSConstrainedDelaunayTriangulation
+
+#include <tut/tut.hpp>
+// geos
+#include <geos_c.h>
+// std
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+#include <memory>
+
+#include "capi_test_utils.h"
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used in test cases.
+struct test_capiconstrainedgeosdelaunaytriangulation_data : public capitest::utility {
+ test_capiconstrainedgeosdelaunaytriangulation_data() {
+ GEOSWKTWriter_setTrim(wktw_, 1);
+ }
+};
+
+typedef test_group<test_capiconstrainedgeosdelaunaytriangulation_data> group;
+typedef group::object object;
+
+group test_capiconstrainedgeosdelaunaytriangulation_group("capi::GEOSConstrainedDelaunayTriangulation");
+
+//
+// Test Cases
+//
+
+// Empty polygon
+template<>
+template<>
+void object::test<1>
+()
+{
+ geom1_ = GEOSGeomFromWKT("POLYGON EMPTY");
+ ensure_equals(GEOSisEmpty(geom1_), 1);
+
+ geom2_ = GEOSConstrainedDelaunayTriangulation(geom1_);
+ ensure (geom2_ != nullptr);
+ ensure_equals(GEOSisEmpty(geom2_), 1);
+ ensure_equals(GEOSGeomTypeId(geom2_), GEOS_GEOMETRYCOLLECTION);
+}
+
+// Single point
+template<>
+template<>
+void object::test<2>
+()
+{
+ geom1_ = GEOSGeomFromWKT("POINT(0 0)");
+ geom2_ = GEOSConstrainedDelaunayTriangulation(geom1_);
+ ensure_equals(GEOSisEmpty(geom2_), 1);
+ ensure_equals(GEOSGeomTypeId(geom2_), GEOS_GEOMETRYCOLLECTION);
+}
+
+// Polygon
+template<>
+template<>
+void object::test<3>
+()
+{
+ geom1_ = GEOSGeomFromWKT("POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10))");
+ geom2_ = GEOSGeomFromWKT("GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 10, 10 10)), POLYGON ((90 90, 20 40, 90 10, 90 90)))");
+ geom3_ = GEOSConstrainedDelaunayTriangulation(geom1_);
+ ensure_geometry_equals(geom2_, geom3_);
+}
+
+
+
+} // namespace tut
+
diff --git a/tests/unit/triangulate/polygon/ConstrainedDelaunayTriangulatorTest.cpp b/tests/unit/triangulate/polygon/ConstrainedDelaunayTriangulatorTest.cpp
new file mode 100644
index 0000000..d181b43
--- /dev/null
+++ b/tests/unit/triangulate/polygon/ConstrainedDelaunayTriangulatorTest.cpp
@@ -0,0 +1,126 @@
+//
+// Test Suite for geos::triangulate::polygon::ConstrainedDelaunayTriangulator
+//
+// tut
+#include <tut/tut.hpp>
+#include <utility.h>
+
+// geos
+#include <geos/geom/Geometry.h>
+#include <geos/io/WKTReader.h>
+#include <geos/io/WKTWriter.h>
+#include <geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h>
+
+// std
+#include <stdio.h>
+
+using geos::triangulate::polygon::ConstrainedDelaunayTriangulator;
+using geos::geom::Geometry;
+
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_constraineddelaunay_data {
+
+ geos::io::WKTReader r;
+
+ test_constraineddelaunay_data() {}
+
+ void checkTri(const std::string& wkt, const std::string& wktExpected)
+ {
+ std::unique_ptr<Geometry> geom = r.read(wkt);
+ std::unique_ptr<Geometry> actual = ConstrainedDelaunayTriangulator::triangulate(geom.get());
+ std::unique_ptr<Geometry> expected = r.read(wktExpected);
+ // std::cout << std::endl << "actual" << std::endl << *actual << std::endl;
+ // std::cout << std::endl << "expected" << std::endl << *expected << std::endl;
+ ensure_equals_geometry(expected.get(), actual.get());
+ }
+
+ /**
+ * Check union of result equals original geom
+ * @param wkt
+ */
+ void checkTri(const std::string& wkt)
+ {
+ std::unique_ptr<Geometry> geom = r.read(wkt);
+ std::unique_ptr<Geometry> actual = ConstrainedDelaunayTriangulator::triangulate(geom.get());
+ std::unique_ptr<Geometry> actualUnion = actual->Union();
+ // std::cout << std::endl << "actual" << std::endl << *actual << std::endl;
+ // std::cout << std::endl << "actualUnion" << std::endl << *actualUnion << std::endl;
+ // std::cout << std::endl << "geom" << std::endl << *geom << std::endl;
+ ensure_equals_geometry(geom.get(), actualUnion.get());
+ }
+
+};
+
+
+typedef test_group<test_constraineddelaunay_data> group;
+typedef group::object object;
+
+group test_constraineddelaunay_group("geos::triangulate::polygon::ConstrainedDelaunayTriangulator");
+
+
+// testQuad
+template<>
+template<>
+void object::test<1>()
+{
+ checkTri(
+ "POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 10, 10 10)), POLYGON ((90 90, 20 40, 90 10, 90 90)))"
+ );
+}
+
+// testPent
+template<>
+template<>
+void object::test<2>()
+{
+ checkTri(
+ "POLYGON ((10 10, 20 40, 90 90, 100 50, 90 10, 10 10))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 10, 10 10)), POLYGON ((90 90, 20 40, 100 50, 90 90)), POLYGON ((100 50, 20 40, 90 10, 100 50)))"
+ );
+}
+
+// testHoleCW
+template<>
+template<>
+void object::test<3>()
+{
+ checkTri(
+ "POLYGON ((10 90, 90 90, 90 20, 10 10, 10 90), (30 70, 80 70, 50 30, 30 70))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 10 90, 30 70, 10 10)), POLYGON ((10 10, 30 70, 50 30, 10 10)), POLYGON ((80 70, 30 70, 90 90, 80 70)), POLYGON ((10 90, 30 70, 90 90, 10 90)), POLYGON ((80 70, 90 90, 90 20, 80 70)), POLYGON ((90 20, 10 10, 50 30, 90 20)), POLYGON ((90 20, 50 30, 80 70, 90 20)))"
+ );
+}
+
+// testMultiPolygon
+template<>
+template<>
+void object::test<4>()
+{
+ checkTri(
+ "MULTIPOLYGON (((10 10, 20 50, 50 50, 40 20, 10 10)), ((20 60, 60 60, 90 20, 90 90, 20 60)), ((10 90, 10 70, 40 70, 50 90, 10 90)))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 50, 40 20, 10 10)), POLYGON ((50 50, 20 50, 40 20, 50 50)), POLYGON ((90 90, 90 20, 60 60, 90 90)), POLYGON ((90 90, 60 60, 20 60, 90 90)), POLYGON ((10 70, 10 90, 40 70, 10 70)), POLYGON ((50 90, 10 90, 40 70, 50 90)))"
+ );
+}
+
+// testFail
+template<>
+template<>
+void object::test<5>()
+{
+ checkTri(
+ "POLYGON ((110 170, 138 272, 145 286, 152 296, 160 307, 303 307, 314 301, 332 287, 343 278, 352 270, 385 99, 374 89, 359 79, 178 89, 167 91, 153 99, 146 107, 173 157, 182 163, 191 170, 199 176, 208 184, 218 194, 226 203, 198 252, 188 247, 182 239, 175 231, 167 223, 161 213, 156 203, 155 198, 110 170))"
+ );
+}
+
+
+
+
+} // namespace tut
+
+
diff --git a/tests/unit/triangulate/polygon/PolygonTriangulatorTest.cpp b/tests/unit/triangulate/polygon/PolygonTriangulatorTest.cpp
new file mode 100644
index 0000000..4c3d0d5
--- /dev/null
+++ b/tests/unit/triangulate/polygon/PolygonTriangulatorTest.cpp
@@ -0,0 +1,186 @@
+//
+// Test Suite for geos::triangulate::polygon::PolygonTriangulator
+//
+// tut
+#include <tut/tut.hpp>
+#include <utility.h>
+
+// geos
+#include <geos/geom/Geometry.h>
+#include <geos/io/WKTReader.h>
+// #include <geos/io/WKTWriter.h>
+#include <geos/triangulate/polygon/PolygonTriangulator.h>
+
+// std
+#include <stdio.h>
+
+using geos::triangulate::polygon::PolygonTriangulator;
+using geos::geom::Geometry;
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_polygontriangulator_data {
+
+ geos::io::WKTReader r;
+
+ test_polygontriangulator_data() {}
+
+ void checkTri(const std::string& wkt, const std::string& wktExpected)
+ {
+ std::unique_ptr<Geometry> geom = r.read(wkt);
+ std::unique_ptr<Geometry> actual = PolygonTriangulator::triangulate(geom.get());
+ std::unique_ptr<Geometry> expected = r.read(wktExpected);
+ ensure_equals_geometry(expected.get(), actual.get());
+ }
+
+ /**
+ * Check union of result equals original geom
+ * @param wkt
+ */
+ void checkTri(const std::string& wkt)
+ {
+ std::unique_ptr<Geometry> geom = r.read(wkt);
+ std::unique_ptr<Geometry> actual = PolygonTriangulator::triangulate(geom.get());
+ std::unique_ptr<Geometry> actualUnion = actual->Union();
+ // std::cout << std::endl << "actual" << std::endl << *actual << std::endl;
+ // std::cout << std::endl << "actualUnion" << std::endl << *actualUnion << std::endl;
+ // std::cout << std::endl << "geom" << std::endl << *geom << std::endl;
+ ensure_equals_geometry(geom.get(), actualUnion.get());
+ }
+
+};
+
+
+typedef test_group<test_polygontriangulator_data> group;
+typedef group::object object;
+
+group test_polygontriangulator_group("geos::triangulate::polygon::PolygonTriangulator");
+
+
+// testQuad
+template<>
+template<>
+void object::test<1>()
+{
+ checkTri(
+ "POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 90, 10 10)), POLYGON ((90 90, 90 10, 10 10, 90 90)))"
+ );
+}
+
+// testPent
+template<>
+template<>
+void object::test<2>()
+{
+ checkTri(
+ "POLYGON ((10 10, 20 40, 90 90, 100 50, 90 10, 10 10))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 90, 10 10)), POLYGON ((90 90, 100 50, 90 10, 90 90)), POLYGON ((90 10, 10 10, 90 90, 90 10)))"
+ );
+}
+
+// testHoleCW
+template<>
+template<>
+void object::test<3>()
+{
+ checkTri(
+ "POLYGON ((10 90, 90 90, 90 20, 10 10, 10 90), (30 70, 80 70, 50 30, 30 70))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 10 90, 30 70, 10 10)), POLYGON ((80 70, 30 70, 10 90, 80 70)), POLYGON ((10 10, 30 70, 50 30, 10 10)), POLYGON ((80 70, 10 90, 90 90, 80 70)), POLYGON ((90 20, 10 10, 50 30, 90 20)), POLYGON ((80 70, 90 90, 90 20, 80 70)), POLYGON ((90 20, 50 30, 80 70, 90 20)))"
+ );
+}
+
+// testTouchingHoles
+template<>
+template<>
+void object::test<4>()
+{
+ checkTri(
+ "POLYGON ((10 90, 90 90, 90 10, 10 10, 10 90), (20 80, 50 70, 30 30, 20 80), (70 20, 50 70, 80 80, 70 20))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 10 90, 20 80, 10 10)), POLYGON ((50 70, 20 80, 10 90, 50 70)), POLYGON ((10 10, 20 80, 30 30, 10 10)), POLYGON ((30 30, 50 70, 70 20, 30 30)), POLYGON ((80 80, 50 70, 10 90, 80 80)), POLYGON ((90 10, 10 10, 30 30, 90 10)), POLYGON ((80 80, 10 90, 90 90, 80 80)), POLYGON ((90 10, 30 30, 70 20, 90 10)), POLYGON ((70 20, 80 80, 90 90, 70 20)), POLYGON ((90 90, 90 10, 70 20, 90 90)))"
+ );
+}
+
+// testRepeatedPoints
+template<>
+template<>
+void object::test<5>()
+{
+ checkTri(
+ "POLYGON ((71 195, 178 335, 178 335, 239 185, 380 210, 290 60, 110 70, 71 195))",
+ "GEOMETRYCOLLECTION (POLYGON ((71 195, 178 335, 239 185, 71 195)), POLYGON ((71 195, 239 185, 290 60, 71 195)), POLYGON ((71 195, 290 60, 110 70, 71 195)), POLYGON ((239 185, 380 210, 290 60, 239 185)))"
+ );
+}
+
+// testMultiPolygon
+template<>
+template<>
+void object::test<6>()
+{
+ checkTri(
+ "MULTIPOLYGON (((10 10, 20 50, 50 50, 40 20, 10 10)), ((20 60, 60 60, 90 20, 90 90, 20 60)), ((10 90, 10 70, 40 70, 50 90, 10 90)))",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 50, 50 50, 10 10)), POLYGON ((50 50, 40 20, 10 10, 50 50)), POLYGON ((90 90, 90 20, 60 60, 90 90)), POLYGON ((60 60, 20 60, 90 90, 60 60)), POLYGON ((10 70, 10 90, 50 90, 10 70)), POLYGON ((50 90, 40 70, 10 70, 50 90)))"
+ );
+}
+
+// testCeeShape
+template<>
+template<>
+void object::test<7>()
+{
+ checkTri(
+ "POLYGON ((110 170, 138 272, 145 286, 152 296, 160 307, 303 307, 314 301, 332 287, 343 278, 352 270, 385 99, 374 89, 359 79, 178 89, 167 91, 153 99, 146 107, 173 157, 182 163, 191 170, 199 176, 208 184, 218 194, 226 203, 198 252, 188 247, 182 239, 175 231, 167 223, 161 213, 156 203, 155 198, 110 170))"
+ );
+}
+
+// Linestring input
+template<>
+template<>
+void object::test<8>()
+{
+ checkTri(
+ "LINESTRING (110 170, 138 272, 145 286, 152 296, 160 307, 303 307, 314 301, 332 287, 343 278, 352 270, 385 99, 374 89, 359 79, 178 89, 167 91, 153 99, 146 107, 173 157, 182 163, 191 170, 199 176, 208 184, 218 194, 226 203, 198 252, 188 247, 182 239, 175 231, 167 223, 161 213, 156 203, 155 198, 110 170)",
+ "GEOMETRYCOLLECTION EMPTY"
+ );
+}
+
+// Empty input
+template<>
+template<>
+void object::test<9>()
+{
+ checkTri(
+ "POLYGON EMPTY",
+ "GEOMETRYCOLLECTION EMPTY"
+ );
+}
+
+// Empty input in collection
+template<>
+template<>
+void object::test<10>()
+{
+ checkTri(
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10)), POLYGON EMPTY)",
+ "GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 90, 10 10)), POLYGON ((90 90, 90 10, 10 10, 90 90)))"
+ );
+}
+
+// Empty mistyped input
+template<>
+template<>
+void object::test<11>()
+{
+ checkTri(
+ "POINT EMPTY",
+ "GEOMETRYCOLLECTION EMPTY"
+ );
+}
+
+} // namespace tut
+
+
diff --git a/tests/unit/triangulate/polygon/VertexSequencePackedRtreeTest.cpp b/tests/unit/triangulate/polygon/VertexSequencePackedRtreeTest.cpp
new file mode 100644
index 0000000..1b0b10c
--- /dev/null
+++ b/tests/unit/triangulate/polygon/VertexSequencePackedRtreeTest.cpp
@@ -0,0 +1,175 @@
+//
+// Test Suite for geos::triangulate::polygon::VertexSequencePackedRtree
+//
+// tut
+#include <tut/tut.hpp>
+
+// geos
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/Point.h>
+#include <geos/io/WKTReader.h>
+#include <geos/io/WKTWriter.h>
+#include <geos/triangulate/polygon/VertexSequencePackedRtree.h>
+
+// std
+#include <stdio.h>
+
+using geos::triangulate::polygon::VertexSequencePackedRtree;
+using geos::geom::Point;
+// using geos::geom::Envelope;
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_packedrtree_data {
+
+ geos::io::WKTReader r;
+ std::vector<Coordinate> coords;
+ std::vector<std::size_t> expected;
+
+ test_packedrtree_data() {}
+
+ std::unique_ptr<VertexSequencePackedRtree>
+ createSPRtree(std::string& multipointwkt)
+ {
+ coords.clear();
+ auto geom_a = r.read(multipointwkt);
+ for (std::size_t i = 0; i < geom_a->getNumGeometries(); i++) {
+ auto pt = dynamic_cast<const Point*>(geom_a->getGeometryN(i));
+ coords.emplace_back(pt->getX(), pt->getY());
+ }
+ std::unique_ptr<VertexSequencePackedRtree> vspr(new VertexSequencePackedRtree(coords));
+ return vspr;
+ }
+
+ Envelope
+ getEnvelope(double minx, double miny, double maxx, double maxy)
+ {
+ return Envelope(minx, maxx, miny, maxy);
+ }
+
+ bool
+ isEqualResult(std::vector<std::size_t>& expectedIds,
+ std::vector<std::size_t>& resultIds)
+ {
+ for (std::size_t i = 0; i < resultIds.size(); i++) {
+ if (expectedIds[i] != resultIds[i])
+ return false;
+ }
+ return true;
+ }
+
+ void
+ checkQuery(std::string& wkt,
+ const Envelope& queryEnv,
+ std::vector<std::size_t> expectedIds)
+ {
+ auto tree = createSPRtree(wkt);
+ std::vector<std::size_t> resultIds;
+ tree->query(queryEnv, resultIds);
+ ensure("result size differs from expected", expected.size() == resultIds.size());
+ ensure("result values differ from expected", isEqualResult(expectedIds, resultIds));
+ }
+
+};
+
+
+typedef test_group<test_packedrtree_data> group;
+typedef group::object object;
+
+group test_packedrtree_group("geos::triangulate::polygon::VertexSequencePackedRtree");
+
+
+// test1
+template<>
+template<>
+void object::test<1>
+()
+{
+ std::string wkt("MULTIPOINT((1 1))");
+ expected = {0};
+ checkQuery(wkt, getEnvelope(1,1, 4,4), expected);
+}
+
+// test2
+template<>
+template<>
+void object::test<2>
+()
+{
+ std::string wkt("MULTIPOINT((0 0), (1 1))");
+ expected = {1};
+ checkQuery(wkt, getEnvelope(1,1, 4,4), expected);
+}
+
+
+// test3
+template<>
+template<>
+void object::test<3>
+()
+{
+ std::string wkt("MULTIPOINT((0 0), (1 1), (2 2), (3 3), (4 4), (5 5))");
+
+ expected = {2,3,4};
+ checkQuery(wkt, getEnvelope(2,2, 4,4), expected);
+
+ expected = {0};
+ checkQuery(wkt, getEnvelope(0,0, 0,0), expected);
+}
+
+// test10
+template<>
+template<>
+void object::test<4>
+()
+{
+ std::string wkt("MULTIPOINT((0 0), (1 1), (2 2), (3 3), (4 4), (5 5), (6 6), (7 7), (8 8), (9 9), (10 10))");
+
+ expected = {2,3,4};
+ checkQuery(wkt, getEnvelope(2,2, 4,4), expected);
+
+ expected = {7,8};
+ checkQuery(wkt, getEnvelope(7,7, 8,8), expected);
+
+ expected = {0};
+ checkQuery(wkt, getEnvelope(0,0, 0,0), expected);
+}
+
+// test6WithDups
+template<>
+template<>
+void object::test<5>
+()
+{
+ std::string wkt("MULTIPOINT((0 0), (1 1), (2 2), (3 3), (4 4), (5 5), (4 4), (3 3), (2 2), (1 1), (0 0))");
+
+ expected = {2,3,4,6,7,8};
+ checkQuery(wkt, getEnvelope(2,2, 4,4), expected);
+
+ expected = {0, 10};
+ checkQuery(wkt, getEnvelope(0,0, 0,0), expected);
+}
+
+
+// test10
+template<>
+template<>
+void object::test<6>
+()
+{
+ std::string wkt("MULTIPOINT((0 0), (1 1), (2 2), (3 3), (4 4), (5 5), (6 6), (7 7), (8 8), (9 9), (10 10), (11 11), (12 12), (13 13), (14 14), (15 15), (16 16), (17 17), (18 18), (17 17), (16 16))");
+
+ expected = {2,3,4};
+ checkQuery(wkt, getEnvelope(2,2, 4,4), expected);
+}
+
+
+
+} // namespace tut
+
+
-----------------------------------------------------------------------
Summary of changes:
.github/workflows/ci.yml | 9 -
NEWS | 3 +
capi/geos_c.cpp | 6 +
capi/geos_c.h.in | 17 +
capi/geos_ts_c.cpp | 11 +
include/geos/geom/Coordinate.h | 2 +
include/geos/geom/Coordinate.inl | 12 +
include/geos/geom/CoordinateArraySequence.h | 2 +
include/geos/geom/Envelope.h | 5 +
include/geos/geom/GeometryFactory.h | 3 +
include/geos/geom/Triangle.h | 74 +++-
.../polygon/ConstrainedDelaunayTriangulator.h | 102 ++++++
.../geos/triangulate/polygon/PolygonEarClipper.h | 217 +++++++++++
.../geos/triangulate/polygon/PolygonHoleJoiner.h | 194 ++++++++++
.../geos/triangulate/polygon/PolygonTriangulator.h | 107 ++++++
.../geos/triangulate/polygon/TriDelaunayImprover.h | 151 ++++++++
.../polygon/VertexSequencePackedRtree.h | 151 ++++++++
.../geos/triangulate/quadedge/TrianglePredicate.h | 19 +-
include/geos/triangulate/quadedge/Vertex.h | 2 +-
include/geos/triangulate/tri/Tri.h | 156 ++++++++
include/geos/triangulate/tri/TriEdge.h | 80 ++++
include/geos/triangulate/tri/TriList.h | 93 +++++
.../geos/triangulate/tri/TriangulationBuilder.h | 86 +++++
src/geom/CoordinateArraySequence.cpp | 8 +
src/geom/Envelope.cpp | 34 ++
src/geom/GeometryFactory.cpp | 13 +
src/geom/Triangle.cpp | 44 +++
.../polygon/ConstrainedDelaunayTriangulator.cpp | 95 +++++
src/triangulate/polygon/PolygonEarClipper.cpp | 368 +++++++++++++++++++
src/triangulate/polygon/PolygonHoleJoiner.cpp | 365 ++++++++++++++++++
src/triangulate/polygon/PolygonTriangulator.cpp | 85 +++++
src/triangulate/polygon/TriDelaunayImprover.cpp | 149 ++++++++
.../polygon/VertexSequencePackedRtree.cpp | 327 +++++++++++++++++
src/triangulate/quadedge/TrianglePredicate.cpp | 8 +-
src/triangulate/tri/Tri.cpp | 406 +++++++++++++++++++++
src/triangulate/tri/TriEdge.cpp | 72 ++++
src/triangulate/tri/TriList.cpp | 75 ++++
src/triangulate/tri/TriangulationBuilder.cpp | 103 ++++++
.../GEOSConstrainedDelaunayTriangulationTest.cpp | 78 ++++
.../ConstrainedDelaunayTriangulatorTest.cpp | 126 +++++++
.../polygon/PolygonTriangulatorTest.cpp | 186 ++++++++++
.../polygon/VertexSequencePackedRtreeTest.cpp | 175 +++++++++
42 files changed, 4194 insertions(+), 25 deletions(-)
create mode 100644 include/geos/triangulate/polygon/ConstrainedDelaunayTriangulator.h
create mode 100644 include/geos/triangulate/polygon/PolygonEarClipper.h
create mode 100644 include/geos/triangulate/polygon/PolygonHoleJoiner.h
create mode 100644 include/geos/triangulate/polygon/PolygonTriangulator.h
create mode 100644 include/geos/triangulate/polygon/TriDelaunayImprover.h
create mode 100644 include/geos/triangulate/polygon/VertexSequencePackedRtree.h
create mode 100644 include/geos/triangulate/tri/Tri.h
create mode 100644 include/geos/triangulate/tri/TriEdge.h
create mode 100644 include/geos/triangulate/tri/TriList.h
create mode 100644 include/geos/triangulate/tri/TriangulationBuilder.h
create mode 100644 src/triangulate/polygon/ConstrainedDelaunayTriangulator.cpp
create mode 100644 src/triangulate/polygon/PolygonEarClipper.cpp
create mode 100644 src/triangulate/polygon/PolygonHoleJoiner.cpp
create mode 100644 src/triangulate/polygon/PolygonTriangulator.cpp
create mode 100644 src/triangulate/polygon/TriDelaunayImprover.cpp
create mode 100644 src/triangulate/polygon/VertexSequencePackedRtree.cpp
create mode 100644 src/triangulate/tri/Tri.cpp
create mode 100644 src/triangulate/tri/TriEdge.cpp
create mode 100644 src/triangulate/tri/TriList.cpp
create mode 100644 src/triangulate/tri/TriangulationBuilder.cpp
create mode 100644 tests/unit/capi/GEOSConstrainedDelaunayTriangulationTest.cpp
create mode 100644 tests/unit/triangulate/polygon/ConstrainedDelaunayTriangulatorTest.cpp
create mode 100644 tests/unit/triangulate/polygon/PolygonTriangulatorTest.cpp
create mode 100644 tests/unit/triangulate/polygon/VertexSequencePackedRtreeTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list