From 0a32f150493afe818c126235be324ad6aa4bfecd Mon Sep 17 00:00:00 2001 From: Paul Ramsey Date: Thu, 16 Dec 2021 15:26:24 -0800 Subject: [PATCH] Port JTS https://github.com/locationtech/jts/pull/810 New OffsetCurve implementation. Closes #530 --- capi/geos_ts_c.cpp | 11 +- include/geos/geom/GeometryFactory.h | 3 + include/geos/geom/LineSegment.h | 16 + include/geos/geom/LineString.h | 3 + include/geos/geom/util/GeometryMapper.h | 115 ++++++ ...veSetBuilder.h => BufferCurveSetBuilder.h} | 34 +- include/geos/operation/buffer/BufferOp.h | 5 + include/geos/operation/buffer/OffsetCurve.h | 271 ++++++++++++++ .../operation/buffer/OffsetCurveBuilder.h | 5 + .../operation/buffer/OffsetSegmentGenerator.h | 32 +- .../geos/operation/buffer/SegmentMCIndex.h | 56 +++ include/geos/operation/union/OverlapUnion.h | 2 +- include/geos/operation/union/UnionStrategy.h | 2 +- src/geom/GeometryFactory.cpp | 9 + src/geom/LineSegment.cpp | 12 + src/geom/LineString.cpp | 9 + src/geom/util/GeometryMapper.cpp | 102 +++++ src/operation/buffer/BufferBuilder.cpp | 10 +- ...tBuilder.cpp => BufferCurveSetBuilder.cpp} | 48 +-- src/operation/buffer/BufferOp.cpp | 14 +- src/operation/buffer/OffsetCurve.cpp | 350 ++++++++++++++++++ src/operation/buffer/OffsetCurveBuilder.cpp | 36 +- src/operation/buffer/SegmentMCIndex.cpp | 62 ++++ tests/unit/capi/GEOSOffsetCurveTest.cpp | 282 +++++++------- tests/unit/geom/LineSegmentTest.cpp | 67 ++++ tests/unit/geom/util/GeometryMapperTest.cpp | 118 ++++++ .../unit/operation/buffer/OffsetCurveTest.cpp | 303 +++++++++++++++ 27 files changed, 1758 insertions(+), 219 deletions(-) create mode 100644 include/geos/geom/util/GeometryMapper.h rename include/geos/operation/buffer/{OffsetCurveSetBuilder.h => BufferCurveSetBuilder.h} (91%) create mode 100644 include/geos/operation/buffer/OffsetCurve.h create mode 100644 include/geos/operation/buffer/SegmentMCIndex.h create mode 100644 src/geom/util/GeometryMapper.cpp rename src/operation/buffer/{OffsetCurveSetBuilder.cpp => BufferCurveSetBuilder.cpp} (90%) create mode 100644 src/operation/buffer/OffsetCurve.cpp create mode 100644 src/operation/buffer/SegmentMCIndex.cpp create mode 100644 tests/unit/geom/util/GeometryMapperTest.cpp create mode 100644 tests/unit/operation/buffer/OffsetCurveTest.cpp diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index bcc2913e5e..09f8018603 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -62,6 +62,7 @@ #include #include #include +#include #include #include #include @@ -170,6 +171,7 @@ using geos::algorithm::distance::DiscreteHausdorffDistance; using geos::operation::buffer::BufferBuilder; using geos::operation::buffer::BufferParameters; +using geos::operation::buffer::OffsetCurve; using geos::operation::distance::IndexedFacetDistance; using geos::operation::geounion::CascadedPolygonUnion; using geos::operation::overlayng::OverlayNG; @@ -1166,13 +1168,8 @@ extern "C" { ); bp.setMitreLimit(mitreLimit); - bool isLeftSide = true; - if(width < 0) { - isLeftSide = false; - width = -width; - } - BufferBuilder bufBuilder(bp); - std::unique_ptr g3 = bufBuilder.bufferLineSingleSided(g1, width, isLeftSide); + OffsetCurve oc(*g1, width, bp); + std::unique_ptr g3 = oc.getCurve(); g3->setSRID(g1->getSRID()); return g3.release(); }); diff --git a/include/geos/geom/GeometryFactory.h b/include/geos/geom/GeometryFactory.h index 3d20f190ef..64906392b5 100644 --- a/include/geos/geom/GeometryFactory.h +++ b/include/geos/geom/GeometryFactory.h @@ -299,6 +299,9 @@ class GEOS_DLL GeometryFactory { std::unique_ptr createLineString( std::unique_ptr && coordinates) const; + std::unique_ptr createLineString( + std::vector && coordinates) const; + /// Construct a LineString with a deep-copy of given argument LineString* createLineString( const CoordinateSequence& coordinates) const; diff --git a/include/geos/geom/LineSegment.h b/include/geos/geom/LineSegment.h index 1d78f7e1f8..95b64a5a75 100644 --- a/include/geos/geom/LineSegment.h +++ b/include/geos/geom/LineSegment.h @@ -216,6 +216,22 @@ class GEOS_DLL LineSegment { double offsetDistance, Coordinate& ret) const; + + /** + * Computes the {@link LineSegment} that is offset from + * the segment by a given distance. + * The computed segment is offset to the left of the line if the offset distance is + * positive, to the right if negative. + * + * @param offsetDistance the distance the point is offset from the segment + * (positive is to the left, negative is to the right) + * @return a line segment offset by the specified distance + * + * @throws IllegalStateException if the segment has zero length + */ + LineSegment offset(double offsetDistance); + + /** \brief * Compute the projection factor for the projection of the point p * onto this LineSegment. diff --git a/include/geos/geom/LineString.h b/include/geos/geom/LineString.h index 6639b1c091..e5240d3414 100644 --- a/include/geos/geom/LineString.h +++ b/include/geos/geom/LineString.h @@ -208,6 +208,9 @@ class GEOS_DLL LineString: public Geometry { LineString(CoordinateSequence::Ptr && pts, const GeometryFactory& newFactory); + LineString(std::vector && pts, + const GeometryFactory& newFactory); + LineString* cloneImpl() const override { return new LineString(*this); } LineString* reverseImpl() const override; diff --git a/include/geos/geom/util/GeometryMapper.h b/include/geos/geom/util/GeometryMapper.h new file mode 100644 index 0000000000..121b6d50c7 --- /dev/null +++ b/include/geos/geom/util/GeometryMapper.h @@ -0,0 +1,115 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2021 Paul Ramsey + * + * 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. + * + ********************************************************************** + * + * Last port: geom/util/GeometryMapper.java + * + **********************************************************************/ + +#pragma once + +#include +#include +#include + +#include + +// Forward declarations +namespace geos { +namespace geom { +class Geometry; +class GeometryCollection; +class GeometryFactory; +} +} + +namespace geos { +namespace geom { // geos.geom +namespace util { // geos.geom.util + + +/** + * Methods to map various collections + * of {@link Geometry}s + * via defined mapping functions. + * + * @author Martin Davis + * + */ +class GEOS_DLL GeometryMapper { + +public: + + /** + * An interface for geometry functions that map a geometry input to a geometry output. + * The output may be nullptr if there is no valid output value for + * the given input value. + */ + typedef std::function(const Geometry&)> mapOp; + + /** + * Maps the members of a {@link Geometry} + * (which may be atomic or composite) + * into another Geometry of most specific type. + * null results are skipped. + * In the case of hierarchical {@link GeometryCollection}s, + * only the first level of members are mapped. + * + * @param geom the input atomic or composite geometry + * @param op the mapping operation + * @return a result collection or geometry of most specific type + */ + static std::unique_ptr map( + const Geometry& geom, + mapOp op); + + /** + * Maps the atomic elements of a {@link Geometry} + * (which may be atomic or composite) + * using a mapOp mapping operation + * into an atomic Geometry or a flat collection + * of the most specific type. + * null and empty values returned from the mapping operation + * are discarded. + * + * @param geom the geometry to map + * @param emptyDim the dimension of empty geometry to create + * @param op the mapping operation + * @return the mapped result + */ + static std::unique_ptr flatMap( + const Geometry& geom, + int emptyDim, + mapOp op); + + + +private: + + static void flatMap( + const Geometry& geom, + mapOp op, + std::vector>& mapped); + + static void addFlat( + std::unique_ptr& geom, + std::vector>& geomList); + + + +}; + +} // namespace geos.geom.util +} // namespace geos.geom +} // namespace geos + diff --git a/include/geos/operation/buffer/OffsetCurveSetBuilder.h b/include/geos/operation/buffer/BufferCurveSetBuilder.h similarity index 91% rename from include/geos/operation/buffer/OffsetCurveSetBuilder.h rename to include/geos/operation/buffer/BufferCurveSetBuilder.h index dbaa5aa823..8858f1c62f 100644 --- a/include/geos/operation/buffer/OffsetCurveSetBuilder.h +++ b/include/geos/operation/buffer/BufferCurveSetBuilder.h @@ -13,15 +13,15 @@ * ********************************************************************** * - * Last port: operation/buffer/OffsetCurveSetBuilder.java r378 (JTS-1.12) + * Last port: operation/buffer/BufferCurveSetBuilder.java r378 (JTS-1.12) * **********************************************************************/ -#ifndef GEOS_OP_BUFFER_OFFSETCURVESETBUILDER_H -#define GEOS_OP_BUFFER_OFFSETCURVESETBUILDER_H +#pragma once #include #include +#include #include @@ -35,6 +35,7 @@ namespace geos { namespace geom { class Geometry; class CoordinateSequence; +class PrecisionModel; class GeometryCollection; class Point; class LineString; @@ -50,6 +51,7 @@ class SegmentString; namespace operation { namespace buffer { class OffsetCurveBuilder; +class BufferParameters; } } } @@ -59,7 +61,7 @@ namespace operation { // geos.operation namespace buffer { // geos.operation.buffer /** - * \class OffsetCurveSetBuilder + * \class BufferCurveSetBuilder * * \brief * Creates all the raw offset curves for a buffer of a Geometry. @@ -68,7 +70,7 @@ namespace buffer { // geos.operation.buffer * final buffer area. * */ -class GEOS_DLL OffsetCurveSetBuilder { +class GEOS_DLL BufferCurveSetBuilder { private: @@ -80,7 +82,7 @@ class GEOS_DLL OffsetCurveSetBuilder { std::vector newLabels; const geom::Geometry& inputGeom; double distance; - OffsetCurveBuilder& curveBuilder; + OffsetCurveBuilder curveBuilder; /// The raw offset curves computed. /// This class holds ownership of std::vector elements. @@ -210,8 +212,8 @@ class GEOS_DLL OffsetCurveSetBuilder { double bufferDistance); // Declare type as noncopyable - OffsetCurveSetBuilder(const OffsetCurveSetBuilder& other) = delete; - OffsetCurveSetBuilder& operator=(const OffsetCurveSetBuilder& rhs) = delete; + BufferCurveSetBuilder(const BufferCurveSetBuilder& other) = delete; + BufferCurveSetBuilder& operator=(const BufferCurveSetBuilder& rhs) = delete; /** * Computes orientation of a ring using a signed-area orientation test. @@ -231,11 +233,20 @@ class GEOS_DLL OffsetCurveSetBuilder { public: /// Constructor - OffsetCurveSetBuilder(const geom::Geometry& newInputGeom, - double newDistance, OffsetCurveBuilder& newCurveBuilder); + BufferCurveSetBuilder( + const geom::Geometry& newInputGeom, + double newDistance, + const geom::PrecisionModel* newPm, + const BufferParameters& newBufParams) + : inputGeom(newInputGeom) + , distance(newDistance) + , curveBuilder(newPm, newBufParams) + , curveList() + , isInvertOrientation(false) + {}; /// Destructor - ~OffsetCurveSetBuilder(); + ~BufferCurveSetBuilder(); /** \brief * Computes the set of raw offset curves for the buffer. @@ -279,4 +290,3 @@ class GEOS_DLL OffsetCurveSetBuilder { #pragma warning(pop) #endif -#endif // ndef GEOS_OP_BUFFER_OFFSETCURVESETBUILDER_H diff --git a/include/geos/operation/buffer/BufferOp.h b/include/geos/operation/buffer/BufferOp.h index 4280c20e60..5cd99719d4 100644 --- a/include/geos/operation/buffer/BufferOp.h +++ b/include/geos/operation/buffer/BufferOp.h @@ -171,6 +171,11 @@ class GEOS_DLL BufferOp { int quadrantSegments = BufferParameters::DEFAULT_QUADRANT_SEGMENTS, int endCapStyle = BufferParameters::CAP_ROUND); + static std::unique_ptr bufferOp( + const geom::Geometry* g, + double distance, + BufferParameters& bufParms); + /** \brief * Initializes a buffer computation for the given geometry. * diff --git a/include/geos/operation/buffer/OffsetCurve.h b/include/geos/operation/buffer/OffsetCurve.h new file mode 100644 index 0000000000..a624157b9d --- /dev/null +++ b/include/geos/operation/buffer/OffsetCurve.h @@ -0,0 +1,271 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2021 Paul Ramsey + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU Lesser General Public Licence as published + * by the Free Software Foundation. + * See the COPYING file for more information. + * + **********************************************************************/ + +#pragma once + +#include +#include +#include +#include + + +#ifdef _MSC_VER +#pragma warning(push) +#pragma warning(disable: 4251) // warning C4251: needs to have dll-interface to be used by clients of class +#endif + +// Forward declarations +namespace geos { +namespace geom { +class Geometry; +class LineString; +class LinearRing; +class Polygon; +class CoordinateSequence; +class Coordinate; +} +namespace operation { +namespace buffer { +class SegmentMCIndex; +} +} +namespace index { +namespace chain { +class MonotoneChain; +} +} +} + +using geos::geom::Geometry; +using geos::geom::GeometryFactory; +using geos::geom::LineString; +using geos::geom::LinearRing; +using geos::geom::Polygon; +using geos::geom::CoordinateSequence; +using geos::geom::Coordinate; + +namespace geos { +namespace operation { +namespace buffer { + +/** + * Computes an offset curve from a geometry. + * The offset curve is a linear geometry which is offset a specified distance + * from the input. + * If the offset distance is positive the curve lies on the left side of the input; + * if it is negative the curve is on the right side. + * + * The offset curve of a line is a LineString which + * The offset curve of a Point is an empty LineString. + * The offset curve of a Polygon is the boundary of the polygon buffer (which + * may be a MultiLineString. + * For a collection the output is a MultiLineString of the element offset curves. + * + * The offset curve is computed as a single contiguous section of the geometry buffer boundary. + * In some geometric situations this definition is ill-defined. + * This algorithm provides a "best-effort" interpretation. + * In particular: + * + * * For self-intersecting lines, the buffer boundary includes + * offset lines for both left and right sides of the input line. + * Only a single contiguous portion on the specified side is returned. + * * If the offset corresponds to buffer holes, only the largest hole is used. + * + * Offset curves support setting the number of quadrant segments, + * the join style, and the mitre limit (if applicable) via + * the BufferParameters. + * + * @author Martin Davis + */ +class GEOS_DLL OffsetCurve { + + +private: + + // Constants + static constexpr int NEARNESS_FACTOR = 10000; + + // Members + const Geometry& inputGeom; + double distance; + BufferParameters bufferParams; + double matchDistance; + const GeometryFactory* geomFactory; + + // Methods + + std::unique_ptr computeCurve(const LineString& lineGeom, double distance); + + std::unique_ptr offsetSegment(const CoordinateSequence* pts, double distance); + + static std::unique_ptr getBufferOriented(const LineString& geom, double distance, BufferParameters& bufParms); + + /** + * Extracts the largest polygon by area from a geometry. + * Used here to avoid issues with non-robust buffer results which have spurious extra polygons. + * + * @param geom a geometry + * @return the polygon element of largest area + */ + static std::unique_ptr extractMaxAreaPolygon(const Geometry& geom); + + static std::unique_ptr extractLongestHole(const Polygon& poly); + + std::unique_ptr computeCurve( + const CoordinateSequence* bufferPts, + std::vector& rawOffsetList); + + int markMatchingSegments(const Coordinate& p0, const Coordinate& p1, + SegmentMCIndex& segIndex, const CoordinateSequence* bufferPts, + std::vector& isInCurve); + + static double subsegmentMatchFrac(const Coordinate& p0, const Coordinate& p1, + const Coordinate& seg0, const Coordinate& seg1, double matchDistance); + + /** + * Extracts a section of a ring of coordinates, starting at a given index, + * and keeping coordinates which are flagged as being required. + * + * @param ring the ring of points + * @param startIndex the index of the start coordinate + * @param isExtracted flag indicating if coordinate is to be extracted + * @return + */ + static void extractSection(const CoordinateSequence* ring, int iStartIndex, + std::vector& isExtracted, std::vector& extractedPoints); + + static std::size_t next(std::size_t i, std::size_t size); + + + /* private */ + class MatchCurveSegmentAction : public index::chain::MonotoneChainSelectAction + { + + private: + + const Coordinate& p0; + const Coordinate& p1; + const CoordinateSequence* bufferPts; + double matchDistance; + std::vector& isInCurve; + double minFrac = -1; + int minCurveIndex = -1; + + public: + + MatchCurveSegmentAction( + const Coordinate& p_p0, const Coordinate& p_p1, + const CoordinateSequence* p_bufferPts, double p_matchDistance, + std::vector& p_isInCurve) + : p0(p_p0) + , p1(p_p1) + , bufferPts(p_bufferPts) + , matchDistance(p_matchDistance) + , isInCurve(p_isInCurve) + , minFrac(-1) + , minCurveIndex(-1) + {}; + + void select(const index::chain::MonotoneChain& mc, std::size_t segIndex) override; + void select(const geom::LineSegment& seg) override { (void)seg; return; }; + + int getMinCurveIndex() { return minCurveIndex; } + }; + + +public: + + /** + * Creates a new instance for computing an offset curve for a geometryat a given distance. + * with default quadrant segments BufferParameters::DEFAULT_QUADRANT_SEGMENTS + * and join style BufferParameters::JOIN_STYLE. + * + * @param geom the geometry to offset + * @param dist the offset distance (positive = left, negative = right) + * + * \see BufferParameters + */ + OffsetCurve(const Geometry& geom, double dist) + : inputGeom(geom) + , distance(dist) + , matchDistance(std::abs(dist)/NEARNESS_FACTOR) + , geomFactory(geom.getFactory()) + {}; + + /** + * Creates a new instance for computing an offset curve for a geometry at a given distance. + * allowing the quadrant segments and join style and mitre limit to be set + * via BufferParameters. + * + * @param geom + * @param dist + * @param bp + */ + OffsetCurve(const Geometry& geom, double dist, BufferParameters& bp) + : inputGeom(geom) + , distance(dist) + , bufferParams(bp) + , matchDistance(std::abs(dist)/NEARNESS_FACTOR) + , geomFactory(geom.getFactory()) + {}; + + /** + * Computes the offset curve of a geometry at a given distance, + * and for a specified quadrant segments, join style and mitre limit. + * + * @param geom a geometry + * @param dist the offset distance (positive = left, negative = right) + * @param quadSegs the quadrant segments (-1 for default) + * @param joinStyle the join style (-1 for default) + * @param mitreLimit the mitre limit (-1 for default) + * @return the offset curve + */ + static std::unique_ptr getCurve( + const Geometry& geom, + double dist, int quadSegs, BufferParameters::JoinStyle joinStyle, double mitreLimit); + + static std::unique_ptr getCurve(const Geometry& geom, double dist); + std::unique_ptr getCurve(); + + /** + * Gets the raw offset line. + * The quadrant segments and join style and mitre limit to be set + * via BufferParameters. + * + * The raw offset line may contain loops and other artifacts which are + * not present in the true offset curve. + * The raw offset line is matched to the buffer ring (which is clean) + * to extract the offset curve. + * + * @param geom the linestring to offset + * @param dist the offset distance + * @param bufParams the buffer parameters to use + * @param lineList the vector to populate with the return value + * @return the raw offset line + */ + static void rawOffset(const LineString& geom, double dist, BufferParameters& bufParams, std::vector& lineList); + static void rawOffset(const LineString& geom, double dist, std::vector& lineList); + +}; + + +} // namespace geos::operation::buffer +} // namespace geos::operation +} // namespace geos + +#ifdef _MSC_VER +#pragma warning(pop) +#endif + + diff --git a/include/geos/operation/buffer/OffsetCurveBuilder.h b/include/geos/operation/buffer/OffsetCurveBuilder.h index ca62587c67..e4448eb19f 100644 --- a/include/geos/operation/buffer/OffsetCurveBuilder.h +++ b/include/geos/operation/buffer/OffsetCurveBuilder.h @@ -154,6 +154,10 @@ class GEOS_DLL OffsetCurveBuilder { double distance, std::vector& lineList); + void getOffsetCurve(const geom::CoordinateSequence* inputPts, + double p_distance, + std::vector& lineList); + private: double distance; @@ -196,6 +200,7 @@ class GEOS_DLL OffsetCurveBuilder { OffsetSegmentGenerator& segGen); + // Declare type as noncopyable OffsetCurveBuilder(const OffsetCurveBuilder& other) = delete; OffsetCurveBuilder& operator=(const OffsetCurveBuilder& rhs) = delete; diff --git a/include/geos/operation/buffer/OffsetSegmentGenerator.h b/include/geos/operation/buffer/OffsetSegmentGenerator.h index 56e1d8c445..c4cab7721d 100644 --- a/include/geos/operation/buffer/OffsetSegmentGenerator.h +++ b/include/geos/operation/buffer/OffsetSegmentGenerator.h @@ -147,6 +147,22 @@ class GEOS_DLL OffsetSegmentGenerator { segList.addPts(pts, isForward); } + /** \brief + * Compute an offset segment for an input segment on a given + * side and at a given distance. + * + * The offset points are computed in full double precision, + * for accuracy. + * + * @param seg the segment to offset + * @param side the side of the segment the offset lies on + * @param distance the offset distance + * @param offset the points computed for the offset segment + */ + static void computeOffsetSegment(const geom::LineSegment& seg, + int side, double distance, + geom::LineSegment& offset); + private: /** @@ -302,22 +318,6 @@ class GEOS_DLL OffsetSegmentGenerator { /// void addInsideTurn(int orientation, bool addStartPoint); - /** \brief - * Compute an offset segment for an input segment on a given - * side and at a given distance. - * - * The offset points are computed in full double precision, - * for accuracy. - * - * @param seg the segment to offset - * @param side the side of the segment the offset lies on - * @param distance the offset distance - * @param offset the points computed for the offset segment - */ - void computeOffsetSegment(const geom::LineSegment& seg, - int side, double distance, - geom::LineSegment& offset); - /** * Adds points for a circular fillet around a reflex corner. * diff --git a/include/geos/operation/buffer/SegmentMCIndex.h b/include/geos/operation/buffer/SegmentMCIndex.h new file mode 100644 index 0000000000..a6ba27a79e --- /dev/null +++ b/include/geos/operation/buffer/SegmentMCIndex.h @@ -0,0 +1,56 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2021 Paul Ramsey + * + * This is free software; you can redistribute and/or modify it under + * the terms of the GNU Lesser General Public Licence as published + * by the Free Software Foundation. + * See the COPYING file for more information. + * + **********************************************************************/ + +#pragma once + +#include +#include + + +// Forward declarations +namespace geos { +namespace geom { +class CoordinateSequence; +class Envelope; +} +} + +using geos::geom::CoordinateSequence; +using geos::geom::Envelope; +using namespace geos::index; + +namespace geos { +namespace operation { +namespace buffer { + +class GEOS_DLL SegmentMCIndex { + +private: + + strtree::TemplateSTRtree index; + std::vector segChains; + + void buildIndex(const CoordinateSequence* segs); + +public: + + SegmentMCIndex(const CoordinateSequence* segs); + + void query(const Envelope* env, chain::MonotoneChainSelectAction& action); +}; + + +} // geos.operation.buffer +} // geos.operation +} // geos diff --git a/include/geos/operation/union/OverlapUnion.h b/include/geos/operation/union/OverlapUnion.h index 383b513720..171a930aab 100644 --- a/include/geos/operation/union/OverlapUnion.h +++ b/include/geos/operation/union/OverlapUnion.h @@ -1,7 +1,7 @@ /********************************************************************** * * GEOS - Geometry Engine Open Source - * http://trac.osgeo.org/geos + * http://geos.osgeo.org * * Copyright (C) 2011 Sandro Santilli * diff --git a/include/geos/operation/union/UnionStrategy.h b/include/geos/operation/union/UnionStrategy.h index 31955655f6..5b560d98c8 100644 --- a/include/geos/operation/union/UnionStrategy.h +++ b/include/geos/operation/union/UnionStrategy.h @@ -1,7 +1,7 @@ /********************************************************************** * * GEOS - Geometry Engine Open Source - * http://trac.osgeo.org/geos + * http://geos.osgeo.org * * Copyright (C) 2020 Paul Ramsey * diff --git a/src/geom/GeometryFactory.cpp b/src/geom/GeometryFactory.cpp index a6f6c30d81..25f7d0bde7 100644 --- a/src/geom/GeometryFactory.cpp +++ b/src/geom/GeometryFactory.cpp @@ -643,6 +643,15 @@ const return std::unique_ptr(new LineString(std::move(newCoords), *this)); } +/*public*/ +std::unique_ptr +GeometryFactory::createLineString(std::vector && newCoords) +const +{ + // Can't use make_unique with protected constructor + return std::unique_ptr(new LineString(std::move(newCoords), *this)); +} + /*public*/ LineString* GeometryFactory::createLineString(const CoordinateSequence& fromCoords) diff --git a/src/geom/LineSegment.cpp b/src/geom/LineSegment.cpp index ee475cd7a5..9d90223eec 100644 --- a/src/geom/LineSegment.cpp +++ b/src/geom/LineSegment.cpp @@ -302,6 +302,18 @@ LineSegment::pointAlongOffset(double segmentLengthFraction, ret = Coordinate(offsetx, offsety); } +/* public */ +LineSegment +LineSegment::offset(double offsetDistance) +{ + Coordinate offset0, offset1; + pointAlongOffset(0, offsetDistance, offset0); + pointAlongOffset(1, offsetDistance, offset1); + LineSegment ls(offset0, offset1); + return ls; +} + + /* public */ std::unique_ptr LineSegment::toGeometry(const GeometryFactory& gf) const diff --git a/src/geom/LineString.cpp b/src/geom/LineString.cpp index 7abed82fa5..38cacadaaf 100644 --- a/src/geom/LineString.cpp +++ b/src/geom/LineString.cpp @@ -106,6 +106,15 @@ LineString::LineString(CoordinateSequence::Ptr && newCoords, validateConstruction(); } +/*public*/ +LineString::LineString(std::vector && newCoords, + const GeometryFactory& factory) + : + Geometry(&factory), + points(new CoordinateArraySequence(std::move(newCoords))) +{ + validateConstruction(); +} std::unique_ptr LineString::getCoordinates() const diff --git a/src/geom/util/GeometryMapper.cpp b/src/geom/util/GeometryMapper.cpp new file mode 100644 index 0000000000..32787ec020 --- /dev/null +++ b/src/geom/util/GeometryMapper.cpp @@ -0,0 +1,102 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2021 Paul Ramsey + * + * 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. + * + ********************************************************************** + * + * Last port: geom/util/GeometryMapper.java + * + **********************************************************************/ + +#include +#include +#include +#include + +namespace geos { +namespace geom { // geos.geom +namespace util { // geos.geom.util + + +/* public static */ +std::unique_ptr +GeometryMapper::map(const Geometry& geom, mapOp op) +{ + std::vector> mapped; + for (std::size_t i = 0; i < geom.getNumGeometries(); i++) { + const Geometry* subgeom = geom.getGeometryN(i); + std::unique_ptr g = op(*subgeom); + if (g != nullptr) + mapped.push_back(std::move(g)); + } + return geom.getFactory()->buildGeometry(std::move(mapped)); +} + + + +/* public static */ +std::unique_ptr +GeometryMapper::flatMap(const Geometry& geom, int emptyDim, mapOp op) +{ + std::vector> mapped; + flatMap(geom, op, mapped); + + if (mapped.empty()) { + return geom.getFactory()->createEmpty(emptyDim); + } + if (mapped.size() == 1) { + std::unique_ptr rg(mapped.at(0).release()); + return rg; + } + return geom.getFactory()->buildGeometry(std::move(mapped)); +} + +/* private static */ +void +GeometryMapper::flatMap(const Geometry& geom, mapOp op, std::vector>& mapped) +{ + for (std::size_t i = 0; i < geom.getNumGeometries(); i++) { + const Geometry* g = geom.getGeometryN(i); + if (g->isCollection()) { + flatMap(*g, op, mapped); + } + else { + std::unique_ptr res = op(*g); + if (res != nullptr && ! res->isEmpty()) { + addFlat(res, mapped); + } + } + } +} + +/* private static */ +void +GeometryMapper::addFlat(std::unique_ptr& geom, std::vector>& geomList) +{ + if (geom->isEmpty()) return; + if (geom->isCollection()) { + auto col = static_cast(geom.get()); + std::vector> geoms = col->releaseGeometries(); + for (auto& subgeom: geoms) { + addFlat(subgeom, geomList); + } + } + else { + geomList.emplace_back(geom.release()); + } +} + + + +} // namespace geos.geom.util +} // namespace geos.geom +} // namespace geos + diff --git a/src/operation/buffer/BufferBuilder.cpp b/src/operation/buffer/BufferBuilder.cpp index a9420349b6..b03b65e552 100644 --- a/src/operation/buffer/BufferBuilder.cpp +++ b/src/operation/buffer/BufferBuilder.cpp @@ -29,7 +29,7 @@ #include #include #include -#include +#include #include #include #include @@ -390,10 +390,8 @@ BufferBuilder::buffer(const Geometry* g, double distance) { // This scope is here to force release of resources owned by - // OffsetCurveSetBuilder when we're doing with it - - OffsetCurveBuilder curveBuilder(precisionModel, bufParams); - OffsetCurveSetBuilder curveSetBuilder(*g, distance, curveBuilder); + // BufferCurveSetBuilder when we're doing with it + BufferCurveSetBuilder curveSetBuilder(*g, distance, precisionModel, bufParams); curveSetBuilder.setInvertOrientation(isInvertOrientation); GEOS_CHECK_FOR_INTERRUPTS(); @@ -401,7 +399,7 @@ BufferBuilder::buffer(const Geometry* g, double distance) std::vector& bufferSegStrList = curveSetBuilder.getCurves(); #if GEOS_DEBUG - std::cerr << "OffsetCurveSetBuilder got " << bufferSegStrList.size() + std::cerr << "BufferCurveSetBuilder got " << bufferSegStrList.size() << " curves" << std::endl; #endif // short-circuit test diff --git a/src/operation/buffer/OffsetCurveSetBuilder.cpp b/src/operation/buffer/BufferCurveSetBuilder.cpp similarity index 90% rename from src/operation/buffer/OffsetCurveSetBuilder.cpp rename to src/operation/buffer/BufferCurveSetBuilder.cpp index 1b49d217c0..d26859bd67 100644 --- a/src/operation/buffer/OffsetCurveSetBuilder.cpp +++ b/src/operation/buffer/BufferCurveSetBuilder.cpp @@ -14,7 +14,7 @@ * ********************************************************************** * - * Last port: operation/buffer/OffsetCurveSetBuilder.java r378 (JTS-1.12) + * Last port: operation/buffer/BufferCurveSetBuilder.java r378 (JTS-1.12) * **********************************************************************/ @@ -23,7 +23,7 @@ #include #include #include -#include +#include #include #include #include @@ -62,17 +62,9 @@ namespace geos { namespace operation { // geos.operation namespace buffer { // geos.operation.buffer -OffsetCurveSetBuilder::OffsetCurveSetBuilder(const Geometry& newInputGeom, - double newDistance, OffsetCurveBuilder& newCurveBuilder): - inputGeom(newInputGeom), - distance(newDistance), - curveBuilder(newCurveBuilder), - curveList(), - isInvertOrientation(false) -{ -} -OffsetCurveSetBuilder::~OffsetCurveSetBuilder() + +BufferCurveSetBuilder::~BufferCurveSetBuilder() { for(std::size_t i = 0, n = curveList.size(); i < n; ++i) { SegmentString* ss = curveList[i]; @@ -85,7 +77,7 @@ OffsetCurveSetBuilder::~OffsetCurveSetBuilder() /* public */ std::vector& -OffsetCurveSetBuilder::getCurves() +BufferCurveSetBuilder::getCurves() { add(inputGeom); return curveList; @@ -93,7 +85,7 @@ OffsetCurveSetBuilder::getCurves() /*public*/ void -OffsetCurveSetBuilder::addCurves(const std::vector& lineList, +BufferCurveSetBuilder::addCurves(const std::vector& lineList, geom::Location leftLoc, geom::Location rightLoc) { for(std::size_t i = 0, n = lineList.size(); i < n; ++i) { @@ -104,7 +96,7 @@ OffsetCurveSetBuilder::addCurves(const std::vector& lineLis /*private*/ void -OffsetCurveSetBuilder::addCurve(CoordinateSequence* coord, +BufferCurveSetBuilder::addCurve(CoordinateSequence* coord, geom::Location leftLoc, geom::Location rightLoc) { #if GEOS_DEBUG @@ -134,7 +126,7 @@ OffsetCurveSetBuilder::addCurve(CoordinateSequence* coord, /*private*/ void -OffsetCurveSetBuilder::add(const Geometry& g) +BufferCurveSetBuilder::add(const Geometry& g) { if(g.isEmpty()) { #if GEOS_DEBUG @@ -173,7 +165,7 @@ OffsetCurveSetBuilder::add(const Geometry& g) /*private*/ void -OffsetCurveSetBuilder::addCollection(const GeometryCollection* gc) +BufferCurveSetBuilder::addCollection(const GeometryCollection* gc) { for(std::size_t i = 0, n = gc->getNumGeometries(); i < n; i++) { const Geometry* g = gc->getGeometryN(i); @@ -183,7 +175,7 @@ OffsetCurveSetBuilder::addCollection(const GeometryCollection* gc) /*private*/ void -OffsetCurveSetBuilder::addPoint(const Point* p) +BufferCurveSetBuilder::addPoint(const Point* p) { // a zero or negative width buffer of a point is empty if(distance <= 0.0) { @@ -201,7 +193,7 @@ OffsetCurveSetBuilder::addPoint(const Point* p) /*private*/ void -OffsetCurveSetBuilder::addLineString(const LineString* line) +BufferCurveSetBuilder::addLineString(const LineString* line) { if (curveBuilder.isLineOffsetEmpty(distance)) { return; @@ -231,7 +223,7 @@ OffsetCurveSetBuilder::addLineString(const LineString* line) /*private*/ void -OffsetCurveSetBuilder::addPolygon(const Polygon* p) +BufferCurveSetBuilder::addPolygon(const Polygon* p) { double offsetDistance = distance; @@ -294,7 +286,7 @@ OffsetCurveSetBuilder::addPolygon(const Polygon* p) /* private */ void -OffsetCurveSetBuilder::addRingBothSides(const CoordinateSequence* coord, double p_distance) +BufferCurveSetBuilder::addRingBothSides(const CoordinateSequence* coord, double p_distance) { addRingSide(coord, p_distance, Position::LEFT, @@ -309,7 +301,7 @@ OffsetCurveSetBuilder::addRingBothSides(const CoordinateSequence* coord, double /* private */ void -OffsetCurveSetBuilder::addRingSide(const CoordinateSequence* coord, +BufferCurveSetBuilder::addRingSide(const CoordinateSequence* coord, double offsetDistance, int side, geom::Location cwLeftLoc, geom::Location cwRightLoc) { @@ -322,7 +314,7 @@ OffsetCurveSetBuilder::addRingSide(const CoordinateSequence* coord, Location leftLoc = cwLeftLoc; Location rightLoc = cwRightLoc; #if GEOS_DEBUG - std::cerr << "OffsetCurveSetBuilder::addPolygonRing: CCW: " << Orientation::isCCW(coord) << std::endl; + std::cerr << "BufferCurveSetBuilder::addPolygonRing: CCW: " << Orientation::isCCW(coord) << std::endl; #endif bool isCCW = isRingCCW(coord); if (coord->size() >= LinearRing::MINIMUM_VALID_SIZE && isCCW) @@ -355,7 +347,7 @@ OffsetCurveSetBuilder::addRingSide(const CoordinateSequence* coord, /* private static*/ bool -OffsetCurveSetBuilder::isRingCurveInverted( +BufferCurveSetBuilder::isRingCurveInverted( const CoordinateSequence* inputPts, double dist, const CoordinateSequence* curvePts) { @@ -396,7 +388,7 @@ OffsetCurveSetBuilder::isRingCurveInverted( */ /* private static */ double -OffsetCurveSetBuilder::maxDistance(const CoordinateSequence* pts, const CoordinateSequence* line) { +BufferCurveSetBuilder::maxDistance(const CoordinateSequence* pts, const CoordinateSequence* line) { double maxDistance = 0; for (std::size_t i = 0; i < pts->size(); i++) { const Coordinate& p = pts->getAt(i); @@ -410,7 +402,7 @@ OffsetCurveSetBuilder::maxDistance(const CoordinateSequence* pts, const Coordin /*private*/ bool -OffsetCurveSetBuilder::isErodedCompletely(const LinearRing* ring, +BufferCurveSetBuilder::isErodedCompletely(const LinearRing* ring, double bufferDistance) { const CoordinateSequence* ringCoord = ring->getCoordinatesRO(); @@ -436,7 +428,7 @@ OffsetCurveSetBuilder::isErodedCompletely(const LinearRing* ring, /*private*/ bool -OffsetCurveSetBuilder::isTriangleErodedCompletely( +BufferCurveSetBuilder::isTriangleErodedCompletely( const CoordinateSequence* triangleCoord, double bufferDistance) { Triangle tri(triangleCoord->getAt(0), triangleCoord->getAt(1), triangleCoord->getAt(2)); @@ -451,7 +443,7 @@ OffsetCurveSetBuilder::isTriangleErodedCompletely( /*private*/ bool -OffsetCurveSetBuilder::isRingCCW(const CoordinateSequence* coords) const +BufferCurveSetBuilder::isRingCCW(const CoordinateSequence* coords) const { bool isCCW = algorithm::Orientation::isCCWArea(coords); //--- invert orientation if required diff --git a/src/operation/buffer/BufferOp.cpp b/src/operation/buffer/BufferOp.cpp index 95dd9e2ccb..7eeee2351f 100644 --- a/src/operation/buffer/BufferOp.cpp +++ b/src/operation/buffer/BufferOp.cpp @@ -87,16 +87,26 @@ BufferOp::precisionScaleFactor(const Geometry* g, /*public static*/ std::unique_ptr -BufferOp::bufferOp(const Geometry* g, double distance, +BufferOp::bufferOp(const Geometry* g, double dist, int quadrantSegments, int nEndCapStyle) { BufferOp bufOp(g); bufOp.setQuadrantSegments(quadrantSegments); bufOp.setEndCapStyle(nEndCapStyle); - return bufOp.getResultGeometry(distance); + return bufOp.getResultGeometry(dist); } +/*public static*/ +std::unique_ptr +BufferOp::bufferOp(const geom::Geometry* g, double dist, + BufferParameters& bufParms) +{ + BufferOp bufOp(g, bufParms); + return bufOp.getResultGeometry(dist); +} + + /*public*/ std::unique_ptr BufferOp::getResultGeometry(double nDistance) diff --git a/src/operation/buffer/OffsetCurve.cpp b/src/operation/buffer/OffsetCurve.cpp new file mode 100644 index 0000000000..d1389ae078 --- /dev/null +++ b/src/operation/buffer/OffsetCurve.cpp @@ -0,0 +1,350 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2021 Paul Ramsey + * + * 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 +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +using namespace geos::index::chain; +using namespace geos::geom; +using geos::geom::util::GeometryMapper; + +namespace geos { +namespace operation { +namespace buffer { + +/* public static */ +std::unique_ptr +OffsetCurve::getCurve(const Geometry& geom, double distance) +{ + OffsetCurve oc(geom, distance); + return oc.getCurve(); +} + + +/* public static */ +std::unique_ptr +OffsetCurve::getCurve(const Geometry& geom, + double dist, + int quadSegs, + BufferParameters::JoinStyle joinStyle, + double mitreLimit) +{ + BufferParameters bufParms; + if (quadSegs >= 0) bufParms.setQuadrantSegments(quadSegs); + if (joinStyle >= 0) bufParms.setJoinStyle(joinStyle); + if (mitreLimit >= 0) bufParms.setMitreLimit(mitreLimit); + OffsetCurve oc(geom, dist, bufParms); + return oc.getCurve(); +} + + +/* public */ +std::unique_ptr +OffsetCurve::getCurve() +{ + GeometryMapper::mapOp GetCurveMapOp = [this](const Geometry& geom)->std::unique_ptr { + + if (geom.getGeometryTypeId() == GEOS_POINT) return nullptr; + if (geom.getGeometryTypeId() == GEOS_POLYGON) { + auto boundary = geom.buffer(distance)->getBoundary(); + + if (boundary->getGeometryTypeId() == GEOS_LINEARRING) { + const LinearRing& ring = static_cast(geom); + auto ringCs = ring.getCoordinatesRO(); + std::unique_ptr ls(geom.getFactory()->createLineString(*ringCs)); + return ls; + } + return boundary; + } + return computeCurve(static_cast(geom), distance); + }; + + return GeometryMapper::flatMap(inputGeom, 1, GetCurveMapOp); +} + + +/* public static */ +void +OffsetCurve::rawOffset(const LineString& geom, double dist, + std::vector& lineList) +{ + BufferParameters bp; + rawOffset(geom, dist, bp, lineList); + return; +} + + +/* public static */ +void +OffsetCurve::rawOffset(const LineString& geom, double distance, BufferParameters& bufParams, + std::vector& lineList) +{ + OffsetCurveBuilder ocb(geom.getFactory()->getPrecisionModel(), bufParams); + ocb.getOffsetCurve(geom.getCoordinatesRO(), distance, lineList); + return; +} + + +/* private */ +std::unique_ptr +OffsetCurve::computeCurve(const LineString& lineGeom, double p_distance) +{ + //-- first handle special/simple cases + if (lineGeom.getNumPoints() < 2 || lineGeom.getLength() == 0.0) { + return geomFactory->createLineString(); + } + if (lineGeom.getNumPoints() == 2) { + return offsetSegment(lineGeom.getCoordinatesRO(), p_distance); + } + + std::vector rawOffsetLines; + rawOffset(lineGeom, p_distance, bufferParams, rawOffsetLines); + if (rawOffsetLines.empty() || rawOffsetLines[0]->size() == 0) { + for (auto* cs: rawOffsetLines) + delete cs; + return geomFactory->createLineString(); + } + /** + * Note: If the raw offset curve has no + * narrow concave angles or self-intersections it could be returned as is. + * However, this is likely to be a less frequent situation, + * and testing indicates little performance advantage, + * so not doing this. + */ + + std::unique_ptr bufferPoly = getBufferOriented(lineGeom, p_distance, bufferParams); + + //-- first try matching shell to raw curve + const CoordinateSequence* shell = bufferPoly->getExteriorRing()->getCoordinatesRO(); + std::unique_ptr offsetCurve = computeCurve(shell, rawOffsetLines); + if (! offsetCurve->isEmpty() || bufferPoly->getNumInteriorRing() == 0) { + for (auto* cs: rawOffsetLines) + delete cs; + return offsetCurve; + } + + //-- if shell didn't work, try matching to largest hole + auto longestHole = extractLongestHole(*bufferPoly); + const CoordinateSequence* holePts = longestHole ? longestHole->getCoordinatesRO() : nullptr; + offsetCurve = computeCurve(holePts, rawOffsetLines); + for (auto* cs: rawOffsetLines) + delete cs; + return offsetCurve; +} + + /* private */ +std::unique_ptr +OffsetCurve::offsetSegment(const CoordinateSequence* pts, double p_distance) +{ + LineSegment ls(pts->getAt(0), pts->getAt(1)); + LineSegment offsetSeg = ls.offset(p_distance); + std::vector coords = { offsetSeg.p0, offsetSeg.p1 }; + return geomFactory->createLineString(std::move(coords)); +} + +/* private static */ +std::unique_ptr +OffsetCurve::getBufferOriented(const LineString& geom, double p_distance, BufferParameters& bufParms) +{ + std::unique_ptr buffer = BufferOp::bufferOp(&geom, std::abs(p_distance), bufParms); + std::unique_ptr bufferPoly = extractMaxAreaPolygon(*buffer); + //-- for negative distances (Right of input) reverse buffer direction to match offset curve + if (p_distance < 0) { + bufferPoly = bufferPoly->reverse(); + } + return bufferPoly; +} + + +/* private static */ +std::unique_ptr +OffsetCurve::extractMaxAreaPolygon(const Geometry& geom) +{ + if (geom.getNumGeometries() == 1) { + const Polygon& poly = static_cast(geom); + poly.clone(); + } + + double maxArea = 0; + const Polygon* maxPoly = nullptr; + for (std::size_t i = 0; i < geom.getNumGeometries(); i++) { + const Polygon* poly = static_cast(geom.getGeometryN(i)); + double area = poly->getArea(); + if (maxPoly == nullptr || area > maxArea) { + maxPoly = poly; + maxArea = area; + } + } + return maxPoly->clone(); +} + +/* private static */ +std::unique_ptr +OffsetCurve::extractLongestHole(const Polygon& poly) +{ + const LinearRing* largestHole = nullptr; + double maxLen = -1; + for (std::size_t i = 0; i < poly.getNumInteriorRing(); i++) { + const LinearRing* hole = poly.getInteriorRingN(i); + double len = hole->getLength(); + if (len > maxLen) { + largestHole = hole; + maxLen = len; + } + } + return largestHole ? largestHole->clone() : nullptr; +} + +/* private */ +std::unique_ptr +OffsetCurve::computeCurve(const CoordinateSequence* bufferPts, std::vector& rawOffsetList) +{ + std::vector isInCurve; + isInCurve.resize(bufferPts->size() - 1, false); + + SegmentMCIndex segIndex(bufferPts); + + int curveStart = -1; + CoordinateSequence* cs = rawOffsetList[0]; + for (std::size_t i = 0; i < cs->size() - 1; i++) { + int index = markMatchingSegments( + cs->getAt(i), cs->getAt(i+1), + segIndex, bufferPts, isInCurve); + if (curveStart < 0) { + curveStart = index; + } + } + std::vector curvePts; + extractSection(bufferPts, curveStart, isInCurve, curvePts); + return geomFactory->createLineString(std::move(curvePts)); +} + +void +OffsetCurve::MatchCurveSegmentAction::select(const MonotoneChain& mc, std::size_t segIndex) +{ + (void)mc; // Quiet unused variable warning + + /** + * A curveRingPt segment may match all or only a portion of a single raw segment. + * There may be multiple curve ring segs that match along the raw segment. + * The one closest to the segment start is recorded as the offset curve start. + */ + double frac = subsegmentMatchFrac(bufferPts->getAt(segIndex), bufferPts->getAt(segIndex+1), p0, p1, matchDistance); + //-- no match + if (frac < 0) return; + + isInCurve[segIndex] = true; + + //-- record lowest index + if (minFrac < 0 || frac < minFrac) { + minFrac = frac; + minCurveIndex = static_cast(segIndex); + } +} + +/* private */ +int +OffsetCurve::markMatchingSegments( + const Coordinate& p0, const Coordinate& p1, + SegmentMCIndex& segIndex, const CoordinateSequence* bufferPts, + std::vector& isInCurve) +{ + Envelope matchEnv(p0, p1); + matchEnv.expandBy(matchDistance); + MatchCurveSegmentAction action(p0, p1, bufferPts, matchDistance, isInCurve); + segIndex.query(&matchEnv, action); + return action.getMinCurveIndex(); +} + + +/* private static */ +double +OffsetCurve::subsegmentMatchFrac(const Coordinate& p0, const Coordinate& p1, + const Coordinate& seg0, const Coordinate& seg1, double matchDistance) +{ + if (matchDistance < algorithm::Distance::pointToSegment(p0, seg0, seg1)) + return -1; + if (matchDistance < algorithm::Distance::pointToSegment(p1, seg0, seg1)) + return -1; + //-- matched - determine position as fraction + LineSegment seg(seg0, seg1); + return seg.segmentFraction(p0); +} + + +/* private static */ +void +OffsetCurve::extractSection(const CoordinateSequence* ring, int iStartIndex, + std::vector& isExtracted, std::vector& extractedPoints) +{ + if (iStartIndex < 0) + return; + + CoordinateList coordList; + std::size_t startIndex = static_cast(iStartIndex); + std::size_t i = startIndex; + do { + coordList.insert(coordList.end(), ring->getAt(i), false); + if (! isExtracted[i]) { + break; + } + i = next(i, ring->size() - 1); + } while (i != startIndex); + //-- handle case where every segment is extracted + if (isExtracted[i]) { + coordList.insert(coordList.end(), ring->getAt(i), false); + } + + //-- if only one point found return empty LineString + if (coordList.size() == 1) + return; + + std::copy(coordList.begin(), coordList.end(), + std::back_inserter(extractedPoints)); + + return; +} + +/* private static */ +std::size_t +OffsetCurve::next(std::size_t i, std::size_t size) { + i += 1; + return (i < size) ? i : 0; +} + + + + +} // namespace geos.operation.buffer +} // namespace geos.operation +} // namespace geos + diff --git a/src/operation/buffer/OffsetCurveBuilder.cpp b/src/operation/buffer/OffsetCurveBuilder.cpp index 4058c8d611..c250084731 100644 --- a/src/operation/buffer/OffsetCurveBuilder.cpp +++ b/src/operation/buffer/OffsetCurveBuilder.cpp @@ -81,6 +81,40 @@ OffsetCurveBuilder::getLineCurve(const CoordinateSequence* inputPts, segGen->getCoordinates(lineList); } +/* public */ +void +OffsetCurveBuilder::getOffsetCurve( + const CoordinateSequence* inputPts, + double p_distance, + std::vector& lineList) +{ + distance = p_distance; + + // a zero width offset curve is empty + if (p_distance == 0.0) return; + bool isRightSide = p_distance < 0.0; + + double posDistance = std::abs(p_distance); + std::unique_ptr segGen = getSegGen(posDistance); + if (inputPts->size() <= 1) { + computePointCurve(inputPts->getAt(0), *segGen); + } + else { + computeSingleSidedBufferCurve(*inputPts, isRightSide, *segGen); + } + + segGen->getCoordinates(lineList); + + // for right side line is traversed in reverse direction, so have to reverse generated line + if (isRightSide) { + for (auto* cs: lineList) { + CoordinateSequence::reverse(cs); + } + } + return; +} + + /* private */ void OffsetCurveBuilder::computePointCurve(const Coordinate& pt, @@ -277,7 +311,7 @@ OffsetCurveBuilder::computeSingleSidedBufferCurve( const CoordinateSequence& inputPts, bool isRightSide, OffsetSegmentGenerator& segGen) { - double distTol = simplifyTolerance(distance); + double distTol = simplifyTolerance(std::abs(distance)); if(isRightSide) { diff --git a/src/operation/buffer/SegmentMCIndex.cpp b/src/operation/buffer/SegmentMCIndex.cpp new file mode 100644 index 0000000000..9441f5a9e0 --- /dev/null +++ b/src/operation/buffer/SegmentMCIndex.cpp @@ -0,0 +1,62 @@ +/********************************************************************** + * + * GEOS - Geometry Engine Open Source + * http://geos.osgeo.org + * + * Copyright (C) 2021 Paul Ramsey + * + * 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 +#include +#include +#include +#include +#include + +using geos::geom::CoordinateSequence; +using geos::geom::Envelope; +using namespace geos::index; + +namespace geos { +namespace operation { +namespace buffer { + +/* public */ +SegmentMCIndex::SegmentMCIndex(const CoordinateSequence* segs) +{ + buildIndex(segs); +} + +/* private */ +void +SegmentMCIndex::buildIndex(const CoordinateSequence* segs) +{ + chain::MonotoneChainBuilder::getChains(segs, nullptr, segChains); + for (chain::MonotoneChain& mc : segChains) { + index.insert(&(mc.getEnvelope()), &mc); + } + return; +} + +/* public */ +void +SegmentMCIndex::query(const Envelope* env, chain::MonotoneChainSelectAction& action) +{ + index.query(*env, [&env, &action](const chain::MonotoneChain* indexChain) { + indexChain->select(*env, action); + return true; + } + ); +} + + +} // namespace geos.operation.buffer +} // namespace geos.operation +} // namespace geos + diff --git a/tests/unit/capi/GEOSOffsetCurveTest.cpp b/tests/unit/capi/GEOSOffsetCurveTest.cpp index 6603221082..9c4651b87a 100644 --- a/tests/unit/capi/GEOSOffsetCurveTest.cpp +++ b/tests/unit/capi/GEOSOffsetCurveTest.cpp @@ -22,6 +22,44 @@ struct test_capioffsetcurve_data : public capitest::utility { test_capioffsetcurve_data() { GEOSWKTWriter_setTrim(wktw_, 1); } + + void debugOutput (const char* label, GEOSGeometry* geom) + { + if (geom) { + char* wkt_c = GEOSWKTWriter_write(wktw_, geom); + printf("%s: %s\n", label, wkt_c); + if (wkt_c) GEOSFree(wkt_c); + } + else { + printf("%s: NULL\n", label); + } + } + + void checkOffset( + const char* wkt, const char* expected, + double width, int quadSegs, int joinStyle, double mitreLimit, + double checkTolerance) + { + static int debug = 0; + // input + geom1_ = GEOSGeomFromWKT(wkt); + ensure(nullptr != geom1_); + if (debug) debugOutput("Input", geom1_); + + // result + geom2_ = GEOSOffsetCurve(geom1_, width, quadSegs, joinStyle, mitreLimit); + ensure(nullptr != geom2_); + if (debug) debugOutput("Result", geom2_); + + // expected + if (expected) { + geom3_ = GEOSGeomFromWKT(expected); + if (debug) debugOutput("Expected", geom3_); + ensure(nullptr != geom3_); + ensure_geometry_equals(geom2_, geom3_, checkTolerance); + } + } + }; typedef test_group group; @@ -36,222 +74,176 @@ group test_capioffsetcurve_group("capi::GEOSOffsetCurve"); // Straight, left template<> template<> -void object::test<1> -() +void object::test<1>() { - geom1_ = GEOSGeomFromWKT("LINESTRING(0 0, 10 0)"); - - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, 2, 0, GEOSBUF_JOIN_ROUND, 2); - - ensure(nullptr != geom2_); - - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - - ensure_equals(std::string(wkt_), std::string("LINESTRING (0 2, 10 2)")); - + checkOffset( + "LINESTRING(0 0, 10 0)", + "LINESTRING (0 2, 10 2)", + 2, 0, GEOSBUF_JOIN_ROUND, 2, + 0.000001 + ); } // Straight, right template<> template<> -void object::test<2> -() +void object::test<2>() { - geom1_ = GEOSGeomFromWKT("LINESTRING(0 0, 10 0)"); - - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, -2, 0, GEOSBUF_JOIN_ROUND, 2); - - ensure(nullptr != geom2_); - - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - - ensure_equals(std::string(wkt_), std::string( - "LINESTRING (10 -2, 0 -2)" - )); + checkOffset( + "LINESTRING(0 0, 10 0)", + "LINESTRING (10 -2, 0 -2)", + -2, 0, GEOSBUF_JOIN_ROUND, 2, + 0.000001 + ); } // Outside curve template<> template<> -void object::test<3> -() +void object::test<3>() { - geom1_ = GEOSGeomFromWKT("LINESTRING(0 0, 10 0, 10 10)"); - - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, -2, 1, GEOSBUF_JOIN_ROUND, 2); - - ensure(nullptr != geom2_); - - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - - ensure_equals(std::string(wkt_), std::string( - "LINESTRING (12 10, 12 0, 10 -2, 0 -2)" - )); + checkOffset( + "LINESTRING(0 0, 10 0, 10 10)", + "LINESTRING (12 10, 12 0, 10 -2, 0 -2)", + -2, 1, GEOSBUF_JOIN_ROUND, 2, + 0.000001 + ); } // Inside curve template<> template<> -void object::test<4> -() +void object::test<4>() { - geom1_ = GEOSGeomFromWKT("LINESTRING(0 0, 10 0, 10 10)"); - - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, 2, 1, GEOSBUF_JOIN_ROUND, 2); - - ensure(nullptr != geom2_); - - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - - ensure_equals(std::string(wkt_), std::string( - "LINESTRING (0 2, 8 2, 8 10)" - )); + checkOffset( + "LINESTRING(0 0, 10 0, 10 10)", + "LINESTRING (0 2, 8 2, 8 10)", + 2, 1, GEOSBUF_JOIN_ROUND, 2, + 0.000001 + ); } // See http://trac.osgeo.org/postgis/ticket/413 template<> template<> -void object::test<5> -() +void object::test<5>() { - geom1_ = GEOSGeomFromWKT("LINESTRING(33282908 6005055,33282900 6005050,33282892 6005042,33282876 6005007,33282863 6004982,33282866 6004971,33282876 6004975,33282967 6005018,33282999 6005031)"); - - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, 44, 1, GEOSBUF_JOIN_MITRE, 1); - - ensure(nullptr != geom2_); - - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - - ensure_equals(std::string(wkt_), std::string( - "LINESTRING EMPTY" - )); + checkOffset( + "LINESTRING(33282908 6005055,33282900 6005050,33282892 6005042,33282876 6005007,33282863 6004982,33282866 6004971,33282876 6004975,33282967 6005018,33282999 6005031)", + // Old algorithm + // "LINESTRING EMPTY", + "LINESTRING (33282915.828957077 6005042.473668674, 33282949.298243735 6005058.282106612, 33282949.3293186 6005058.295747999, 33282982.439409934 6005071.764529393)", + 44, 1, GEOSBUF_JOIN_MITRE, 1, + 0.000001 + ); } // 0 distance // See http://trac.osgeo.org/postgis/ticket/454 template<> template<> -void object::test<6> -() +void object::test<6>() { - geom1_ = GEOSGeomFromWKT("LINESTRING(0 0, 10 0)"); - - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, 0, 0, GEOSBUF_JOIN_ROUND, 2); + checkOffset( + "LINESTRING(0 0, 10 0)", + "LINESTRING (0 0, 10 0)", + 0, 0, GEOSBUF_JOIN_ROUND, 2, + 0.000001 + ); +} - ensure(nullptr != geom2_); +// left-side and right-side curve +// See http://trac.osgeo.org/postgis/ticket/633 +template<> +template<> +void object::test<7>() +{ + const char* wkt0 = "LINESTRING (" + "665.7317504882812500 133.0762634277343700," + "1774.4752197265625000 19.9391822814941410," + "756.2413940429687500 466.8306579589843700," + "626.1337890625000000 1898.0147705078125000," + "433.8007202148437500 404.6052856445312500)"; - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); + double width = 57.164000837203; - ensure_equals(std::string(wkt_), std::string( - "LINESTRING (0 0, 10 0)" - )); + // left-sided + checkOffset( + wkt0, + NULL, + width, 8, GEOSBUF_JOIN_MITRE, 5.57, + 0.000001 + ); + + ensure(GEOSGeomGetNumPoints(geom2_) >= GEOSGeomGetNumPoints(geom1_)); } // left-side and right-side curve // See http://trac.osgeo.org/postgis/ticket/633 template<> template<> -void object::test<7> -() +void object::test<8>() { - std::string wkt0("LINESTRING (" + const char* wkt0 = "LINESTRING (" "665.7317504882812500 133.0762634277343700," "1774.4752197265625000 19.9391822814941410," "756.2413940429687500 466.8306579589843700," "626.1337890625000000 1898.0147705078125000," - "433.8007202148437500 404.6052856445312500)"); - - geom1_ = GEOSGeomFromWKT(wkt0.c_str()); - ensure(nullptr != geom1_); + "433.8007202148437500 404.6052856445312500)"; double width = 57.164000837203; - // left-sided - { - geom2_ = GEOSOffsetCurve(geom1_, width, 8, GEOSBUF_JOIN_MITRE, 5.57); - ensure(nullptr != geom2_); - // likely, 5 >= 5 - ensure(GEOSGeomGetNumPoints(geom2_) >= GEOSGeomGetNumPoints(geom1_)); - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - GEOSGeom_destroy(geom2_); - GEOSFree(wkt_); - //ensure_equals(std::string(wkt_), ...); - } - // right-sided - { - width = -width; - geom2_ = GEOSOffsetCurve(geom1_, width, 8, GEOSBUF_JOIN_MITRE, 5.57); - ensure(nullptr != geom2_); - // likely, 5 >= 7 - ensure(GEOSGeomGetNumPoints(geom2_) >= GEOSGeomGetNumPoints(geom1_)); - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - //ensure_equals(std::string(wkt_), ...); - } + checkOffset( + wkt0, + NULL, + -1 * width, 8, GEOSBUF_JOIN_MITRE, 5.57, + 0.000001 + ); + + ensure(GEOSGeomGetNumPoints(geom2_) >= GEOSGeomGetNumPoints(geom1_)); } // Test duplicated inner vertex in input // See http://trac.osgeo.org/postgis/ticket/602 template<> template<> -void object::test<8> -() +void object::test<9>() { - double width = -1; - - geom1_ = GEOSGeomFromWKT("LINESTRING(0 0,0 10,0 10,10 10)"); - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, width, 8, GEOSBUF_JOIN_ROUND, 0); - ensure("Unexpected exception", nullptr != geom2_); - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - ensure_equals(std::string(wkt_), "LINESTRING (10 9, 1 9, 1 0)"); + checkOffset( + "LINESTRING(0 0,0 10,0 10,10 10)", + "LINESTRING (10 9, 1 9, 1 0)", + -1, 8, GEOSBUF_JOIN_ROUND, 0, + 0.000001 + ); } // Test duplicated final vertex in input // See http://trac.osgeo.org/postgis/ticket/602 template<> template<> -void object::test<9> -() +void object::test<10>() { - double width = -1; - - geom1_ = GEOSGeomFromWKT("LINESTRING(0 0,0 10,0 10)"); - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, width, 8, GEOSBUF_JOIN_ROUND, 0); - ensure("Unexpected exception", nullptr != geom2_); - wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - ensure_equals(std::string(wkt_), "LINESTRING (1 10, 1 0)"); + checkOffset( + "LINESTRING(0 0,0 10,0 10)", + "LINESTRING (1 10, 1 0)", + -1, 8, GEOSBUF_JOIN_ROUND, 0, + 0.000001 + ); } // Test only duplicated vertex in input // See http://trac.osgeo.org/postgis/ticket/602 template<> template<> -void object::test<10> -() +void object::test<11>() { - double width = -1; - - geom1_ = GEOSGeomFromWKT("LINESTRING(0 10,0 10,0 10)"); - ensure(nullptr != geom1_); - - geom2_ = GEOSOffsetCurve(geom1_, width, 8, GEOSBUF_JOIN_ROUND, 0); - ensure("Missing expected exception", nullptr == geom2_); + checkOffset( + "LINESTRING(0 10, 0 10, 0 10)", + "LINESTRING EMPTY", + -1, 8, GEOSBUF_JOIN_ROUND, 0, + 0.000001 + ); } } // namespace tut diff --git a/tests/unit/geom/LineSegmentTest.cpp b/tests/unit/geom/LineSegmentTest.cpp index eb95aae8dc..d930e28cb7 100644 --- a/tests/unit/geom/LineSegmentTest.cpp +++ b/tests/unit/geom/LineSegmentTest.cpp @@ -2,9 +2,12 @@ // Test Suite for geos::geom::LineSegment class. #include +#include + // geos #include #include + // std #include @@ -39,6 +42,36 @@ struct test_lineseg_data { ensure("checkLineIntersection", dist <= MAX_ABS_ERROR_INTERSECTION); } + void checkOffsetPoint( + double x0, double y0, + double x1, double y1, + double segFrac, double offset, + double expectedX, double expectedY) + { + LineSegment seg(x0, y0, x1, y1); + Coordinate actual; + seg.pointAlongOffset(segFrac, offset, actual); + Coordinate expected(expectedX, expectedY); + ensure_equals(actual.x, expected.x); + ensure_equals(actual.y, expected.y); + } + + void checkOffsetLine( + double x0, double y0, + double x1, double y1, + double offset, + double expectedX0, double expectedY0, + double expectedX1, double expectedY1) + { + LineSegment seg(x0, y0, x1, y1); + LineSegment actual = seg.offset(offset); + + Coordinate expected0(expectedX0, expectedY0); + Coordinate expected1(expectedX1, expectedY1); + ensure_equals_xyz(actual.p0, expected0); + ensure_equals_xyz(actual.p1, expected1); + } + test_lineseg_data() : ph1(0, 2), ph2(10, 2), pv1(0, 0), pv2(0, 10), h1(ph1, ph2), v1(pv1, pv2) {} @@ -153,5 +186,39 @@ void object::test<7> 35613477.772841461, 4257160.5339209242 ); } +// testOffsetLine +template<> +template<> +void object::test<8>() +{ + const double ROOT2 = std::sqrt(2.0); + checkOffsetLine(0, 0, 10, 10, 0, 0, 0, 10, 10 ); + checkOffsetLine(0, 0, 10, 10, ROOT2, -1, 1, 9, 11 ); + checkOffsetLine(0, 0, 10, 10, -ROOT2, 1, -1, 11, 9); +} + +// testOffsetPoint +template<> +template<> +void object::test<9>() +{ + double ROOT2 = std::sqrt(2.0); + checkOffsetPoint(0, 0, 10, 10, 0.0, ROOT2, -1, 1); + checkOffsetPoint(0, 0, 10, 10, 0.0, -ROOT2, 1, -1); + + checkOffsetPoint(0, 0, 10, 10, 1.0, ROOT2, 9, 11); + checkOffsetPoint(0, 0, 10, 10, 0.5, ROOT2, 4, 6); + + checkOffsetPoint(0, 0, 10, 10, 0.5, -ROOT2, 6, 4); + checkOffsetPoint(0, 0, 10, 10, 0.5, -ROOT2, 6, 4); + + checkOffsetPoint(0, 0, 10, 10, 2.0, ROOT2, 19, 21); + checkOffsetPoint(0, 0, 10, 10, 2.0, -ROOT2, 21, 19); + + checkOffsetPoint(0, 0, 10, 10, 2.0, 5 * ROOT2, 15, 25); + checkOffsetPoint(0, 0, 10, 10, -2.0, 5 * ROOT2, -25, -15); +} + + } // namespace tut diff --git a/tests/unit/geom/util/GeometryMapperTest.cpp b/tests/unit/geom/util/GeometryMapperTest.cpp new file mode 100644 index 0000000000..9bff66e859 --- /dev/null +++ b/tests/unit/geom/util/GeometryMapperTest.cpp @@ -0,0 +1,118 @@ +// +// Test Suite for geos::geom::util::GeometryCombiner class. + +// tut +#include +// geos +#include +#include + +#include +// std +#include + +namespace tut { + +using namespace geos::geom; +using geos::geom::util::GeometryMapper; + + +// Common data used by tests +struct test_geometrymapper_data { + + geos::io::WKTReader wktreader_; + + GeometryMapper::mapOp KEEP_LINE = [](const Geometry& geom)->std::unique_ptr { + if (geom.getGeometryTypeId() == GEOS_POINT) { + return geom.getFactory()->createEmpty(1); + } + if (geom.getGeometryTypeId() == GEOS_LINESTRING) { + return geom.clone(); + } + return nullptr; + }; + + GeometryMapper::mapOp BOUNDARY = [](const Geometry& geom)->std::unique_ptr { + return geom.getBoundary(); + }; + + void checkFlatMap(const std::string& wkt, int dim, GeometryMapper::mapOp op, const std::string& wktExpected) + { + auto geom = wktreader_.read(wkt); + auto actual = GeometryMapper::flatMap(*geom, dim, op); + auto expected = wktreader_.read(wktExpected); + ensure_equals_geometry(actual.get(), expected.get()); + } + +}; + +typedef test_group group; +typedef group::object object; + +group test_geometrymapper_group("geos::geom::util::GeometryMapper"); + + +// testFlatMapInputEmpty +template<> +template<> +void object::test<1>() +{ + checkFlatMap( + "GEOMETRYCOLLECTION(POINT EMPTY, LINESTRING EMPTY)", + 1, KEEP_LINE, "LINESTRING EMPTY"); +} + +// testFlatMapInputMulti +template<> +template<> +void object::test<2>() +{ + checkFlatMap( + "GEOMETRYCOLLECTION( MULTILINESTRING((0 0, 1 1), (1 1, 2 2)), LINESTRING(2 2, 3 3))", + 1, KEEP_LINE, "MULTILINESTRING ((0 0, 1 1), (1 1, 2 2), (2 2, 3 3))"); +} + +// testFlatMapResultEmpty +template<> +template<> +void object::test<3>() +{ + checkFlatMap( + "GEOMETRYCOLLECTION( LINESTRING(0 0, 1 1), LINESTRING(1 1, 2 2))", + 1, KEEP_LINE, "MULTILINESTRING((0 0, 1 1), (1 1, 2 2))"); + + checkFlatMap( + "GEOMETRYCOLLECTION( POINT(0 0), POINT(0 0), LINESTRING(0 0, 1 1))", + 1, KEEP_LINE, "LINESTRING(0 0, 1 1)"); + + checkFlatMap( + "MULTIPOINT((0 0), (1 1))", + 1, KEEP_LINE, "LINESTRING EMPTY"); +} + +// testFlatMapResultNull +template<> +template<> +void object::test<4>() +{ + checkFlatMap( + "GEOMETRYCOLLECTION( POINT(0 0), LINESTRING(0 0, 1 1), POLYGON ((1 1, 1 2, 2 1, 1 1)))", + 1, KEEP_LINE, "LINESTRING(0 0, 1 1)"); +} + +// testFlatMapBoundary +template<> +template<> +void object::test<5>() +{ + checkFlatMap( + "GEOMETRYCOLLECTION( POINT(0 0), LINESTRING(0 0, 1 1), POLYGON ((1 1, 1 2, 2 1, 1 1)))", + 0, BOUNDARY, "GEOMETRYCOLLECTION (POINT (0 0), POINT (1 1), LINESTRING (1 1, 1 2, 2 1, 1 1))"); + + checkFlatMap( + "LINESTRING EMPTY", + 0, BOUNDARY, "POINT EMPTY"); +} + + +} // namespace tut diff --git a/tests/unit/operation/buffer/OffsetCurveTest.cpp b/tests/unit/operation/buffer/OffsetCurveTest.cpp new file mode 100644 index 0000000000..8a6ff0217c --- /dev/null +++ b/tests/unit/operation/buffer/OffsetCurveTest.cpp @@ -0,0 +1,303 @@ +// +// Test Suite for geos::operation::buffer::BufferOp class. + +// tut +#include +#include + +// geos +#include +#include +#include +#include +// #include + +// std +#include +#include + +namespace tut { +// +// Test Group +// +using geos::operation::buffer::OffsetCurve; +using geos::operation::buffer::BufferParameters; + +// Common data used by tests +struct test_offsetcurve_data { + geos::io::WKTReader wktreader; + + test_offsetcurve_data() {}; + + void checkOffsetCurve(const std::string& wkt, double distance, + int quadSegs, int joinStyle, double mitreLimit, + const std::string& wktExpected) + { + checkOffsetCurve(wkt, distance, quadSegs, joinStyle, mitreLimit, wktExpected, 0.05); + } + + void checkOffsetCurve(const std::string& wkt, double distance, + int quadSegs, int joinStyle, double mitreLimit, + const std::string& wktExpected, double tolerance) + { + BufferParameters::JoinStyle js = static_cast(joinStyle); + std::unique_ptr geom = wktreader.read(wkt); + std::unique_ptr result = OffsetCurve::getCurve(*geom, distance, quadSegs, js, mitreLimit); + std::unique_ptr expected = wktreader.read(wktExpected); + ensure_equals_geometry(result.get(), expected.get(), tolerance); + } + + void checkOffsetCurve(const std::string& wkt, double distance, const std::string& wktExpected) + { + checkOffsetCurve(wkt, distance, wktExpected, 0.05); + } + + void checkOffsetCurve(const std::string& wkt, double distance, const std::string& wktExpected, double tolerance) + { + std::unique_ptr geom = wktreader.read(wkt); + std::unique_ptr result = OffsetCurve::getCurve(*geom, distance); + std::unique_ptr expected = wktreader.read(wktExpected); + + // geos::io::WKTWriter wktwriter; + // wktwriter.setRoundingPrecision(2); + // std::cout << std::endl; + // std::cout << "Expected: " << wktwriter.write(expected.get()) << std::endl; + // std::cout << " Result: " << wktwriter.write(result.get()) << std::endl; + + ensure_equals_geometry(result.get(), expected.get(), tolerance); + } + +}; + +typedef test_group group; +typedef group::object object; + +group test_offsetcurve_group("geos::operation::buffer::OffsetCurve"); + + +// testPoint +template<> +template<> +void object::test<1> () +{ + checkOffsetCurve( + "POINT (0 0)", 1, + "LINESTRING EMPTY" + ); +} + +// testEmpty +template<> +template<> +void object::test<2> () +{ + checkOffsetCurve( + "LINESTRING EMPTY", 1, + "LINESTRING EMPTY" + ); +} + +// testZeroLenLine +template<> +template<> +void object::test<3> () +{ + checkOffsetCurve( + "LINESTRING (1 1, 1 1)", 1, + "LINESTRING EMPTY" + ); +} + +// testSegment1Short +template<> +template<> +void object::test<4> () +{ + checkOffsetCurve( + "LINESTRING (2 2, 2 2.0000001)", 1, + "LINESTRING (1 2, 1 2.0000001)", + 0.00000001 + ); +} + +// testSegment1 +template<> +template<> +void object::test<5> () +{ + checkOffsetCurve( + "LINESTRING (0 0, 9 9)", 1, + "LINESTRING (-0.71 0.71, 8.29 9.71)" + ); +} + +// testSegments2 +template<> +template<> +void object::test<6> () +{ + checkOffsetCurve( + "LINESTRING (0 0, 9 9, 25 0, 30 15)", 1, + "LINESTRING (-0.71 0.71, 8.29 9.71, 8.44 9.83, 8.6 9.92, 8.77 9.97, 8.96 10, 9.14 9.99, 9.32 9.95, 9.49 9.87, 24.43 1.47, 29.05 15.32)" + ); +} + +// testZigzagOneEndCurved4 +template<> +template<> +void object::test<7> () +{ + checkOffsetCurve( + "LINESTRING (1 3, 6 3, 4 5, 9 5)", 1, + "LINESTRING (1 4, 3.59 4, 3.29 4.29, 3.17 4.44, 3.08 4.62, 3.02 4.8, 3 5, 3.02 5.2, 3.08 5.38, 3.17 5.56, 3.29 5.71, 3.44 5.83, 3.62 5.92, 3.8 5.98, 4 6, 9 6)" + ); +} + +// testEmptyResult +template<> +template<> +void object::test<8> () +{ + checkOffsetCurve( + "LINESTRING (3 5, 5 7, 7 5)", -4, + "LINESTRING EMPTY" + ); +} + +// testSelfCross +template<> +template<> +void object::test<9> () +{ + checkOffsetCurve( + "LINESTRING (50 90, 50 10, 90 50, 10 50)", 10, + "LINESTRING (60 90, 60 60)" ); +} + +// testSelfCrossNeg +template<> +template<> +void object::test<10> () +{ + checkOffsetCurve( + "LINESTRING (50 90, 50 10, 90 50, 10 50)", -10, + "LINESTRING (40 90, 40 60, 10 60)" ); +} + +// testRing +template<> +template<> +void object::test<11> () +{ + checkOffsetCurve( + "LINESTRING (10 10, 50 90, 90 10, 10 10)", -10, + "LINESTRING (26.18 20, 50 67.63, 73.81 20, 26.18 20)" ); +} + +// testClosedCurve +template<> +template<> +void object::test<12> () +{ + checkOffsetCurve( + "LINESTRING (30 70, 80 80, 50 10, 10 80, 60 70)", 10, + "LINESTRING (45 83.2, 78.04 89.81, 80 90, 81.96 89.81, 83.85 89.23, 85.59 88.29, 87.11 87.04, 88.35 85.5, 89.27 83.76, 89.82 81.87, 90 79.9, 89.79 77.94, 89.19 76.06, 59.19 6.06, 58.22 4.3, 56.91 2.77, 55.32 1.53, 53.52 0.64, 51.57 0.12, 49.56 0.01, 47.57 0.3, 45.68 0.98, 43.96 2.03, 42.49 3.4, 41.32 5.04, 1.32 75.04, 0.53 76.77, 0.09 78.63, 0.01 80.53, 0.29 82.41, 0.93 84.2, 1.89 85.85, 3.14 87.28, 4.65 88.45, 6.34 89.31, 8.17 89.83, 10.07 90, 11.96 89.81, 45 83.2)" + ); +} + +// testMultiLine +template<> +template<> +void object::test<13> () +{ + checkOffsetCurve( + "MULTILINESTRING ((20 30, 60 10, 80 60), (40 50, 80 30))", 10, + "MULTILINESTRING ((24.47 38.94, 54.75 23.8, 70.72 63.71), (44.47 58.94, 84.47 38.94))" + ); +} + +// testPolygon +template<> +template<> +void object::test<14> () +{ + checkOffsetCurve( + "POLYGON ((100 200, 200 100, 100 100, 100 200))", 10, + "LINESTRING (90 200, 90.19 201.95, 90.76 203.83, 91.69 205.56, 92.93 207.07, 94.44 208.31, 96.17 209.24, 98.05 209.81, 100 210, 101.95 209.81, 103.83 209.24, 105.56 208.31, 107.07 207.07, 207.07 107.07, 208.31 105.56, 209.24 103.83, 209.81 101.95, 210 100, 209.81 98.05, 209.24 96.17, 208.31 94.44, 207.07 92.93, 205.56 91.69, 203.83 90.76, 201.95 90.19, 200 90, 100 90, 98.05 90.19, 96.17 90.76, 94.44 91.69, 92.93 92.93, 91.69 94.44, 90.76 96.17, 90.19 98.05, 90 100, 90 200)" + ); +} + +// testPolygon +template<> +template<> +void object::test<15> () +{ + checkOffsetCurve( + "POLYGON ((100 200, 200 100, 100 100, 100 200))", -10, + "LINESTRING (110 175.86, 175.86 110, 110 110, 110 175.86)" + ); +} + +// testPolygonWithHole +template<> +template<> +void object::test<16> () +{ + checkOffsetCurve( + "POLYGON ((20 80, 80 80, 80 20, 20 20, 20 80), (30 70, 70 70, 70 30, 30 30, 30 70))", 10, + "MULTILINESTRING ((10 80, 10.19 81.95, 10.76 83.83, 11.69 85.56, 12.93 87.07, 14.44 88.31, 16.17 89.24, 18.05 89.81, 20 90, 80 90, 81.95 89.81, 83.83 89.24, 85.56 88.31, 87.07 87.07, 88.31 85.56, 89.24 83.83, 89.81 81.95, 90 80, 90 20, 89.81 18.05, 89.24 16.17, 88.31 14.44, 87.07 12.93, 85.56 11.69, 83.83 10.76, 81.95 10.19, 80 10, 20 10, 18.05 10.19, 16.17 10.76, 14.44 11.69, 12.93 12.93, 11.69 14.44, 10.76 16.17, 10.19 18.05, 10 20, 10 80), (40 60, 40 40, 60 40, 60 60, 40 60))" + ); +} + +// testPolygonWithHole +template<> +template<> +void object::test<17> () +{ + checkOffsetCurve( + "POLYGON ((20 80, 80 80, 80 20, 20 20, 20 80), (30 70, 70 70, 70 30, 30 30, 30 70))", -10, + "LINESTRING EMPTY" + ); +} + + //--------------------------------------- + +// testQuadSegs +template<> +template<> +void object::test<18> () +{ + checkOffsetCurve( + "LINESTRING (20 20, 50 50, 80 20)", + 10, 2, -1, -1, + "LINESTRING (12.93 27.07, 42.93 57.07, 50 60, 57.07 57.07, 87.07 27.07)" + ); +} + +// testJoinBevel +template<> +template<> +void object::test<19> () +{ + checkOffsetCurve( + "LINESTRING (20 20, 50 50, 80 20)", + 10, -1, BufferParameters::JOIN_BEVEL, -1, + "LINESTRING (12.93 27.07, 42.93 57.07, 57.07 57.07, 87.07 27.07)" + ); +} + +// testJoinMitre +template<> +template<> +void object::test<20> () +{ + checkOffsetCurve( + "LINESTRING (20 20, 50 50, 80 20)", + 10, -1, BufferParameters::JOIN_MITRE, -1, + "LINESTRING (12.93 27.07, 50 64.14, 87.07 27.07)" + ); +} + + + +} // namespace tut