[geos-commits] [SCM] GEOS branch main updated. 68a5b4d4bfadee1c580fb6702597fe7cfbcc8118
git at osgeo.org
git at osgeo.org
Mon May 29 07:14:36 PDT 2023
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 68a5b4d4bfadee1c580fb6702597fe7cfbcc8118 (commit)
via 523cea8a3afef830cee3b2f9b591bc4dbcd4c98e (commit)
via 2efd5ed41009149db31b896cea98aacf765b004c (commit)
from bfe64bd74a6d79bdf30eed7adc0ecc8f14d210e4 (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 68a5b4d4bfadee1c580fb6702597fe7cfbcc8118
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Mon May 29 07:13:30 2023 -0700
* Port https://github.com/locationtech/jts/pull/981
* References https://github.com/qgis/QGIS/issues/53165
* References https://github.com/libgeos/geos/issues/897
diff --git a/include/geos/operation/buffer/OffsetCurve.h b/include/geos/operation/buffer/OffsetCurve.h
index 1a9f450be..403ca98b7 100644
--- a/include/geos/operation/buffer/OffsetCurve.h
+++ b/include/geos/operation/buffer/OffsetCurve.h
@@ -190,6 +190,12 @@ public:
// Constants
static constexpr int MATCH_DISTANCE_FACTOR = 10000;
+ /**
+ * A QuadSegs minimum value that will prevent generating
+ * unwanted offset curve artifacts near end caps.
+ */
+ static constexpr int MIN_QUADRANT_SEGMENTS = 8;
+
/**
* Creates a new instance for computing an offset curve for a geometry at a given distance.
* with default quadrant segments (BufferParameters::DEFAULT_QUADRANT_SEGMENTS)
@@ -230,7 +236,17 @@ public:
throw util::IllegalArgumentException("OffsetCurve distance must be a finite value");
}
//-- set buffer params, leaving cap style as the default CAP_ROUND
- bufferParams.setQuadrantSegments( bp.getQuadrantSegments());
+
+ /**
+ * Prevent using a very small QuadSegs value, to avoid
+ * offset curve artifacts near the end caps.
+ */
+ int quadSegs = bp.getQuadrantSegments();
+ if (quadSegs < MIN_QUADRANT_SEGMENTS) {
+ quadSegs = MIN_QUADRANT_SEGMENTS;
+ }
+ bufferParams.setQuadrantSegments(quadSegs);
+
bufferParams.setJoinStyle( bp.getJoinStyle());
bufferParams.setMitreLimit( bp.getMitreLimit());
};
diff --git a/tests/unit/capi/GEOSOffsetCurveTest.cpp b/tests/unit/capi/GEOSOffsetCurveTest.cpp
index c2ab6805c..052a5f9df 100644
--- a/tests/unit/capi/GEOSOffsetCurveTest.cpp
+++ b/tests/unit/capi/GEOSOffsetCurveTest.cpp
@@ -104,7 +104,7 @@ void object::test<3>()
{
checkOffset(
"LINESTRING(0 0, 10 0, 10 10)",
- "LINESTRING (12 10, 12 0, 10 -2, 0 -2)",
+ "LINESTRING (0 -2, 10 -2, 10.3901806 -1.9615705, 10.76536686 -1.8477590, 11.11114046 -1.66293922, 11.41421356 -1.41421356, 11.66293922 -1.11114046, 11.84775906 -0.76536686, 11.96157056 -0.3901806, 12 0, 12 10)",
-2, 1, GEOSBUF_JOIN_ROUND, 2,
0.000001
);
@@ -132,7 +132,7 @@ void object::test<5>()
"LINESTRING(33282908 6005055,33282900 6005050,33282892 6005042,33282876 6005007,33282863 6004982,33282866 6004971,33282876 6004975,33282967 6005018,33282999 6005031)",
// Old algorithm
// "LINESTRING EMPTY",
- "LINESTRING (33282939.63869393 6005053.735950421, 33282948.754269257 6005058.043310191, 33282949.873293087 6005058.534544423, 33282982.439409934 6005071.764529393)",
+ "LINESTRING (33282951.601378817 6005059.236579252, 33282982.439409934 6005071.764529393)",
44, 1, GEOSBUF_JOIN_MITRE, 1,
0.000001
);
diff --git a/tests/unit/operation/buffer/OffsetCurveTest.cpp b/tests/unit/operation/buffer/OffsetCurveTest.cpp
index 8ac3a2226..5faea0fe4 100644
--- a/tests/unit/operation/buffer/OffsetCurveTest.cpp
+++ b/tests/unit/operation/buffer/OffsetCurveTest.cpp
@@ -525,7 +525,7 @@ void object::test<38> ()
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)"
+ "LINESTRING (12.928932188134524 27.071067811865476, 42.928932188134524 57.071067811865476, 44.44429766980398 58.314696123025456, 46.1731656763491 59.23879532511287, 48.04909677983872 59.80785280403231, 50 60, 51.95090322016128 59.80785280403231, 53.8268343236509 59.23879532511287, 55.55570233019602 58.314696123025456, 57.071067811865476 57.071067811865476, 87.07106781186548 27.071067811865476)"
);
}
@@ -554,5 +554,32 @@ void object::test<40> ()
}
+// testMinQuadrantSegments
+// See https://github.com/qgis/QGIS/issues/53165
+template<>
+template<>
+void object::test<41> ()
+{
+ checkOffsetCurve(
+ "LINESTRING (553772.0645892698 177770.05079236583, 553780.9235869241 177768.99614978794, 553781.8325485934 177768.41771963477)",
+ -11, 0, BufferParameters::JOIN_MITRE, -1,
+ "LINESTRING (553770.76 177759.13, 553777.54 177758.32)"
+ );
+}
+
+// testMinQuadrantSegments_QGIS
+// See https://github.com/qgis/QGIS/issues/53165#issuecomment-1563214857
+template<>
+template<>
+void object::test<42> ()
+{
+ checkOffsetCurve(
+ "LINESTRING (421 622, 446 625, 449 627)",
+ 133, 0, BufferParameters::JOIN_MITRE, -1,
+ "LINESTRING (405.15 754.05, 416.3 755.39)"
+ );
+}
+
+
} // namespace tut
diff --git a/web/content/specifications/geojson.md b/web/content/specifications/geojson.md
index 502f8ea3e..fb3cee7de 100644
--- a/web/content/specifications/geojson.md
+++ b/web/content/specifications/geojson.md
@@ -14,7 +14,7 @@ Unlike [WKB]({{< ref "wkb" >}}) and [WKT]({{< ref "wkt" >}}), GeoJSON does not r
* [Feature](https://datatracker.ietf.org/doc/html/rfc7946#section-3.2), representation of an object that has a "geometry" and an arbitrary set of other non-geometric "properties".
* [FeatureCollection](https://datatracker.ietf.org/doc/html/rfc7946#section-3.3), representation of a list of Features.
-Since GEOS is almost 100% concerned with geometry operations, there is no GEOS abstraction to translate the non-geometric properties of GeoJSON Feature into, so handling of Feature properties is only available via the C++ [GeoJSON.h](/libgeos/geos/blob/main/include/geos/io/GeoJSON.h) utilities.
+Since GEOS is almost 100% concerned with geometry operations, there is no GEOS abstraction to translate the non-geometric properties of GeoJSON Feature into, so handling of Feature properties is only available via the C++ [GeoJSON.h](https://github.com/libgeos/geos/blob/main/include/geos/io/GeoJSON.h) utilities.
### Writing GeoJSON
commit 523cea8a3afef830cee3b2f9b591bc4dbcd4c98e
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Fri May 26 15:49:06 2023 -0700
Update NEWS for MinimumAreaRectangle
diff --git a/NEWS.md b/NEWS.md
index 881e17287..08cda262c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -23,6 +23,7 @@ xxxx-xx-xx
- Polygonal coverages: CoverageSimplifier (JTS-911, Martin Davis)
- CAPI: GEOSCoverageIsValid, GEOSCoverageSimplifyVW (GH-867, Paul Ramsey)
- CAPI: GEOSGeom_releaseCollection (GH-848)
+ - CAPI: GEOSMinimumRotatedRectangle now uses MinimumAreaRectangle (Paul Ramsey)
- Fixes/Improvements:
- WKTReader: Fix parsing of Z and M flags in WKTReader (#676 and GH-669, Dan Baston)
commit 2efd5ed41009149db31b896cea98aacf765b004c
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Fri May 26 15:46:03 2023 -0700
Add MinimumAreaRectangle, a bounding rectangle algorithm.
Port of https://github.com/locationtech/jts/pull/977
This adds a MinimumAreaRectangle class implementing the
standard "rotating-calipers` algorithm for computing a Minimum-Area Rectangle.
MinimumDiameter.getMinimumRectangle was previously used
for this, but it does not always compute the Minimum-Area
Rectangle. It is kept available, for backwards compatibility.
Also, it computes the Minimum-Width Rectangle, which may be useful.
The CAPI GEOSMinimumRotatedRectangle is now bound to
MinimumAreaRectangle rather than MinimumDiameter.getMinimumRectangle
diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp
index 9941292e5..710004c4e 100644
--- a/capi/geos_ts_c.cpp
+++ b/capi/geos_ts_c.cpp
@@ -20,6 +20,7 @@
#include <geos/algorithm/BoundaryNodeRule.h>
#include <geos/algorithm/MinimumBoundingCircle.h>
#include <geos/algorithm/MinimumDiameter.h>
+#include <geos/algorithm/MinimumAreaRectangle.h>
#include <geos/algorithm/Orientation.h>
#include <geos/algorithm/construct/MaximumInscribedCircle.h>
#include <geos/algorithm/construct/LargestEmptyCircle.h>
@@ -1323,9 +1324,10 @@ extern "C" {
Geometry*
GEOSMinimumRotatedRectangle_r(GEOSContextHandle_t extHandle, const Geometry* g)
{
+ using geos::algorithm::MinimumAreaRectangle;
+
return execute(extHandle, [&]() {
- geos::algorithm::MinimumDiameter m(g);
- auto g3 = m.getMinimumRectangle();
+ auto g3 = MinimumAreaRectangle::getMinimumRectangle(g);
g3->setSRID(g->getSRID());
return g3.release();
});
diff --git a/include/geos/algorithm/Distance.h b/include/geos/algorithm/Distance.h
index b74b87bc0..5d82384d5 100644
--- a/include/geos/algorithm/Distance.h
+++ b/include/geos/algorithm/Distance.h
@@ -49,10 +49,11 @@ public:
* another point of the line (must be different to A)
*/
// formerly distanceLineLine
- static double segmentToSegment(const geom::CoordinateXY& A,
- const geom::CoordinateXY& B,
- const geom::CoordinateXY& C,
- const geom::CoordinateXY& D);
+ static double segmentToSegment(
+ const geom::CoordinateXY& A,
+ const geom::CoordinateXY& B,
+ const geom::CoordinateXY& C,
+ const geom::CoordinateXY& D);
/**
* Computes the distance from a point to a sequence of line segments.
@@ -63,8 +64,9 @@ public:
* a sequence of contiguous line segments defined by their vertices
* @return the minimum distance between the point and the line segments
*/
- static double pointToSegmentString(const geom::CoordinateXY& p,
- const geom::CoordinateSequence* seq);
+ static double pointToSegmentString(
+ const geom::CoordinateXY& p,
+ const geom::CoordinateSequence* seq);
/**
* Computes the distance from a point p to a line segment AB
@@ -80,9 +82,10 @@ public:
* @return the distance from p to line segment AB
*/
// formerly distancePointLine
- static double pointToSegment(const geom::CoordinateXY& p,
- const geom::CoordinateXY& A,
- const geom::CoordinateXY& B);
+ static double pointToSegment(
+ const geom::CoordinateXY& p,
+ const geom::CoordinateXY& A,
+ const geom::CoordinateXY& B);
/**
* Computes the perpendicular distance from a point p to the (infinite) line
@@ -97,10 +100,15 @@ public:
* @return the distance from p to line AB
*/
// formerly distancePointLinePerpendicular
- static double pointToLinePerpendicular(const geom::CoordinateXY& p,
- const geom::CoordinateXY& A,
- const geom::CoordinateXY& B);
+ static double pointToLinePerpendicular(
+ const geom::CoordinateXY& p,
+ const geom::CoordinateXY& A,
+ const geom::CoordinateXY& B);
+ static double pointToLinePerpendicularSigned(
+ const geom::CoordinateXY& p,
+ const geom::CoordinateXY& A,
+ const geom::CoordinateXY& B);
};
} // namespace geos::algorithm
diff --git a/include/geos/algorithm/MinimumAreaRectangle.h b/include/geos/algorithm/MinimumAreaRectangle.h
new file mode 100644
index 000000000..0e340a981
--- /dev/null
+++ b/include/geos/algorithm/MinimumAreaRectangle.h
@@ -0,0 +1,159 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2023 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 <vector>
+#include <memory>
+
+// Forward declarations
+namespace geos {
+ namespace geom {
+ class CoordinateSequence;
+ class CoordinateXY;
+ class Geometry;
+ class GeometryFactory;
+ class LineSegment;
+ class LineString;
+ class Polygon;
+ }
+}
+
+using geos::geom::CoordinateSequence;
+using geos::geom::CoordinateXY;
+using geos::geom::Geometry;
+using geos::geom::GeometryFactory;
+using geos::geom::LineSegment;
+using geos::geom::LineString;
+using geos::geom::Polygon;
+
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+
+
+/**
+ * Computes the minimum-area rectangle enclosing a Geometry.
+ * Unlike the Envelope, the rectangle may not be axis-parallel.
+ *
+ * The first step in the algorithm is computing the convex hull of the Geometry.
+ * If the input Geometry is known to be convex, a hint can be supplied to
+ * avoid this computation.
+ *
+ * In degenerate cases the minimum enclosing geometry
+ * may be a LineString or a Point.
+ *
+ * The minimum-area enclosing rectangle does not necessarily
+ * have the minimum possible width.
+ * Use MinimumDiameter to compute this.
+ *
+ * @see MinimumDiameter
+ * @see ConvexHull
+ *
+ */
+class GEOS_DLL MinimumAreaRectangle {
+
+private:
+
+ // Members
+ const Geometry* m_inputGeom;
+ bool m_isConvex;
+
+ // Methods
+ std::unique_ptr<Geometry> getMinimumRectangle();
+
+ std::unique_ptr<Geometry> computeConvex(const Geometry* convexGeom);
+
+ /**
+ * Computes the minimum-area rectangle for a convex ring of Coordinate.
+ *
+ * This algorithm uses the "dual rotating calipers" technique.
+ * Performance is linear in the number of segments.
+ *
+ * @param ring the convex ring to scan
+ */
+ std::unique_ptr<Polygon> computeConvexRing(const CoordinateSequence* ring);
+
+ std::size_t findFurthestVertex(
+ const CoordinateSequence* pts,
+ const LineSegment& baseSeg,
+ std::size_t startIndex,
+ int orient);
+
+ bool isFurtherOrEqual(double d1, double d2, int orient);
+
+ static double orientedDistance(
+ const LineSegment& seg,
+ const CoordinateXY& p,
+ int orient);
+
+ static std::size_t getNextIndex(
+ const CoordinateSequence* ring,
+ std::size_t index);
+
+ /**
+ * Creates a line of maximum extent from the provided vertices
+ * @param pts the vertices
+ * @param factory the geometry factory
+ * @return the line of maximum extent
+ */
+ static std::unique_ptr<LineString> computeMaximumLine(
+ const CoordinateSequence* pts,
+ const GeometryFactory* factory);
+
+
+public:
+
+ /**
+ * Compute a minimum-area rectangle for a given Geometry.
+ *
+ * @param inputGeom a Geometry
+ */
+ MinimumAreaRectangle(const Geometry* inputGeom)
+ : m_inputGeom(inputGeom)
+ , m_isConvex(false)
+ {};
+
+ /**
+ * Compute a minimum rectangle for a Geometry,
+ * with a hint if the geometry is convex
+ * (e.g. a convex Polygon or LinearRing,
+ * or a two-point LineString, or a Point).
+ *
+ * @param inputGeom a Geometry which is convex
+ * @param isConvex true if the input geometry is convex
+ */
+ MinimumAreaRectangle(const Geometry* inputGeom, bool isConvex)
+ : m_inputGeom(inputGeom)
+ , m_isConvex(isConvex)
+ {};
+
+ /**
+ * Gets the minimum-area rectangular Polygon which encloses the input geometry.
+ * If the convex hull of the input is degenerate (a line or point)
+ * a LineString or Point is returned.
+ *
+ * @param geom the geometry
+ * @return the minimum rectangle enclosing the geometry
+ */
+ static std::unique_ptr<Geometry> getMinimumRectangle(const Geometry* geom);
+
+};
+
+
+} // namespace geos::algorithm
+} // namespace geos
+
diff --git a/include/geos/algorithm/MinimumDiameter.h b/include/geos/algorithm/MinimumDiameter.h
index 956f90d23..ada2d84b7 100644
--- a/include/geos/algorithm/MinimumDiameter.h
+++ b/include/geos/algorithm/MinimumDiameter.h
@@ -3,6 +3,7 @@
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
+ * Copyright (C) 2023 Paul Ramsey <pramsey at cleverelephant.ca>
* Copyright (C) 2005-2006 Refractions Research Inc.
* Copyright (C) 2001-2002 Vivid Solutions Inc.
*
@@ -11,10 +12,6 @@
* by the Free Software Foundation.
* See the COPYING file for more information.
*
- **********************************************************************
- *
- * Last port: algorithm/MinimumDiameter.java r966
- *
**********************************************************************/
#pragma once
@@ -53,12 +50,16 @@ namespace algorithm { // geos::algorithm
*
* This class can also be used to compute a line segment representing
* the minimum diameter, the supporting line segment of the minimum diameter,
- * and a minimum rectangle enclosing the input geometry.
+ * and a minimum-width rectangle of the input geometry.
* This rectangle will have width equal to the minimum diameter, and have
* one side parallel to the supporting segment.
*
- * @see ConvexHull
+ * In degenerate cases the rectangle may be a LineString or a Point.
+ * (Note that this may not be the enclosing rectangle with minimum area;
+ * use MinimumAreaRectangle to compute this.)
*
+ * @see ConvexHull
+ * @see MinimumAreaRectangle
*/
class GEOS_DLL MinimumDiameter {
private:
@@ -75,7 +76,7 @@ private:
void computeWidthConvex(const geom::Geometry* geom);
/**
- * Compute the width information for a ring of {@link geom::Coordinate}s.
+ * Compute the width information for a ring of Coordinate.
* Leaves the width information in the instance variables.
*
* @param pts
@@ -127,7 +128,7 @@ public:
double getLength();
/** \brief
- * Gets the {@link geom::Coordinate} forming one end of the minimum diameter.
+ * Gets the geom::Coordinate forming one end of the minimum diameter.
*
* @return a coordinate forming one end of the minimum diameter
*/
@@ -148,15 +149,17 @@ public:
std::unique_ptr<geom::LineString> getDiameter();
/** \brief
- * Gets the minimum rectangular Polygon which encloses the input geometry.
+ * Gets the rectangular Polygon which encloses the input geometry
+ * and is based on the minimum diameter supporting segment.
*
* The rectangle has width equal to the minimum diameter, and a longer
* length. If the convex hill of the input is degenerate (a line or point)
* a LineString or Point is returned.
- * The minimum rectangle can be used as an extremely generalized
- * representation for the given geometry.
+ * This is not necessarily the rectangle with minimum area.
+ * Use MinimumAreaRectangle to compute this.
*
- * @return the minimum rectangle enclosing the input (or a line or point if degenerate)
+ * @return the the minimum-width rectangle enclosing the geometry
+ * @see MinimumAreaRectangle
*/
std::unique_ptr<geom::Geometry> getMinimumRectangle();
@@ -164,7 +167,8 @@ public:
* Gets the minimum rectangle enclosing a geometry.
*
* @param geom the geometry
- * @return the minimum rectangle enclosing the geometry
+ * @return a rectangle enclosing the input (or a line or point if degenerate)
+ * @see MinimumAreaRectangle
*/
static std::unique_ptr<geom::Geometry> getMinimumRectangle(geom::Geometry* geom);
diff --git a/include/geos/algorithm/Rectangle.h b/include/geos/algorithm/Rectangle.h
new file mode 100644
index 000000000..ef073ee47
--- /dev/null
+++ b/include/geos/algorithm/Rectangle.h
@@ -0,0 +1,96 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2023 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 <geos/geom/LineSegment.h>
+
+// Forward declarations
+namespace geos {
+namespace geom {
+class Coordinate;
+class GeometryFactory;
+class Polygon;
+}
+}
+
+using geos::geom::CoordinateXY;
+using geos::geom::GeometryFactory;
+using geos::geom::LineSegment;
+using geos::geom::Polygon;
+
+namespace geos {
+namespace algorithm {
+
+class GEOS_DLL Rectangle {
+
+public:
+
+ /**
+ * Creates a rectangular {@link Polygon} from a base segment
+ * defining the position and orientation of one side of the rectangle, and
+ * three points defining the locations of the line segments
+ * forming the opposite, left and right sides of the rectangle.
+ * The base segment and side points must be presented so that the
+ * rectangle has CW orientation.
+ *
+ * The rectangle corners are computed as intersections of
+ * lines, which generally cannot produce exact values.
+ * If a rectangle corner is determined to coincide with a side point
+ * the side point value is used to avoid numerical inaccuracy.
+ *
+ * The first side of the constructed rectangle contains the base segment.
+ *
+ * @param baseRightPt the right point of the base segment
+ * @param baseLeftPt the left point of the base segment
+ * @param oppositePt the point defining the opposite side
+ * @param leftSidePt the point defining the left side
+ * @param rightSidePt the point defining the right side
+ * @param factory the geometry factory to use
+ * @return the rectangular polygon
+ */
+ static std::unique_ptr<Polygon>
+ createFromSidePts(
+ const CoordinateXY& baseRightPt,
+ const CoordinateXY& baseLeftPt,
+ const CoordinateXY& oppositePt,
+ const CoordinateXY& leftSidePt,
+ const CoordinateXY& rightSidePt,
+ const GeometryFactory* factory);
+
+private:
+
+ /**
+ * Computes the constant C in the standard line equation Ax + By = C
+ * from A and B and a point on the line.
+ *
+ * @param a the X coefficient
+ * @param b the Y coefficient
+ * @param p a point on the line
+ * @return the constant C
+ */
+ static double
+ computeLineEquationC(double a, double b, const CoordinateXY& p);
+
+ static LineSegment
+ createLineForStandardEquation(double a, double b, double c);
+
+};
+
+
+} // namespace geos::algorithm
+} // namespace geos
+
+
diff --git a/include/geos/geom/LineSegment.h b/include/geos/geom/LineSegment.h
index 21da3212a..714c123d9 100644
--- a/include/geos/geom/LineSegment.h
+++ b/include/geos/geom/LineSegment.h
@@ -188,7 +188,21 @@ public:
return orientationIndex(*seg);
};
-
+ /**
+ * Determines the orientation index of a Coordinate relative to this segment.
+ * The orientation index is as defined in Orientation::index(Coordinate, Coordinate, Coordinate).
+ *
+ * @param p the coordinate to compare
+ *
+ * @return 1 (LEFT) if "p" is to the left of this segment
+ * @return -1 (RIGHT) if "p" is to the right of this segment
+ * @return 0 (COLLINEAR) if "p" is collinear with this segment
+ *
+ */
+ int orientationIndex(const CoordinateXY& p) const
+ {
+ return algorithm::Orientation::index(p0, p1, p);
+ }
/** \brief
* Determines the orientation index of a Coordinate
@@ -255,15 +269,35 @@ public:
return algorithm::Distance::pointToSegment(p, p0, p1);
};
- /** \brief
+ /**
* Computes the perpendicular distance between the (infinite)
* line defined by this line segment and a point.
+ * If the segment has zero length this returns the distance between
+ * the segment and the point.
+ *
+ * @param p the point to compute the distance to
+ * @return the perpendicular distance between the line and point
*/
double distancePerpendicular(const CoordinateXY& p) const
{
+ if (p0.equals2D(p1))
+ return p0.distance(p);
return algorithm::Distance::pointToLinePerpendicular(p, p0, p1);
};
+ /**
+ * Computes the oriented perpendicular distance between the (infinite) line
+ * defined by this line segment and a point.
+ * The oriented distance is positive if the point on the left of the line,
+ * and negative if it is on the right.
+ * If the segment has zero length this returns the distance between
+ * the segment and the point.
+ *
+ * @param p the point to compute the distance to
+ * @return the oriented perpendicular distance between the line and point
+ */
+ double distancePerpendicularOriented(const CoordinateXY& p) const;
+
/** \brief
* Computes the Coordinate that lies a given
* fraction along the line defined by this segment.
@@ -375,6 +409,8 @@ public:
*/
void project(const Coordinate& p, Coordinate& ret) const;
+ CoordinateXY project(const CoordinateXY& p) const;
+
/** \brief
* Project a line segment onto this line segment and return the resulting
* line segment.
diff --git a/src/algorithm/Distance.cpp b/src/algorithm/Distance.cpp
index 278a4002f..ebbba3523 100644
--- a/src/algorithm/Distance.cpp
+++ b/src/algorithm/Distance.cpp
@@ -201,6 +201,29 @@ Distance::pointToSegmentString(const geom::CoordinateXY& p,
}
+/*public static*/
+double
+Distance::pointToLinePerpendicularSigned(
+ const geom::CoordinateXY& p,
+ const geom::CoordinateXY& A,
+ const geom::CoordinateXY& B)
+{
+ // use comp.graphics.algorithms Frequently Asked Questions method
+ /*
+ * (2) s = (Ay-Cy)(Bx-Ax)-(Ax-Cx)(By-Ay)
+ * -----------------------------
+ * L^2
+ *
+ * Then the distance from C to P = |s|*L.
+ */
+ double len2 = (B.x - A.x) * (B.x - A.x) + (B.y - A.y) * (B.y - A.y);
+ double s = ((A.y - p.y) * (B.x - A.x) - (A.x - p.x) * (B.y - A.y))
+ / len2;
+
+ return s * std::sqrt(len2);
+}
+
+
} // namespace geos.algorithm
} //namespace geos
diff --git a/src/algorithm/MinimumAreaRectangle.cpp b/src/algorithm/MinimumAreaRectangle.cpp
new file mode 100644
index 000000000..b3f655fb9
--- /dev/null
+++ b/src/algorithm/MinimumAreaRectangle.cpp
@@ -0,0 +1,269 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2023 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/MinimumAreaRectangle.h>
+#include <geos/algorithm/ConvexHull.h>
+#include <geos/algorithm/Rectangle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LineSegment.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/Point.h>
+#include <geos/geom/Polygon.h>
+#include <geos/util.h>
+#include <geos/util/IllegalArgumentException.h>
+#include <geos/constants.h>
+
+using namespace geos::geom;
+
+
+namespace geos {
+namespace algorithm { // geos.algorithm
+
+
+/* public static */
+std::unique_ptr<Geometry>
+MinimumAreaRectangle::getMinimumRectangle(const Geometry* geom)
+{
+ MinimumAreaRectangle mar(geom);
+ return mar.getMinimumRectangle();
+}
+
+
+/* private */
+std::unique_ptr<Geometry>
+MinimumAreaRectangle::getMinimumRectangle()
+{
+ if (m_inputGeom->isEmpty()) {
+ return m_inputGeom->getFactory()->createPolygon();
+ }
+ if (m_isConvex) {
+ return computeConvex(m_inputGeom);
+ }
+ ConvexHull cvh(m_inputGeom);
+ std::unique_ptr<Geometry> convexGeom = cvh.getConvexHull();
+
+ return computeConvex(convexGeom.get());
+}
+
+
+/* private */
+std::unique_ptr<Geometry>
+MinimumAreaRectangle::computeConvex(const Geometry* convexGeom)
+{
+ const CoordinateSequence* convexHullPts;
+ switch (convexGeom->getGeometryTypeId()) {
+ case GEOS_POLYGON:
+ {
+ auto poly = static_cast<const Polygon*>(convexGeom);
+ convexHullPts = poly->getExteriorRing()->getCoordinatesRO();
+ break;
+ }
+ case GEOS_LINESTRING:
+ {
+ auto line = static_cast<const LineString*>(convexGeom);
+ convexHullPts = line->getCoordinatesRO();
+ break;
+ }
+ case GEOS_POINT:
+ {
+ auto pt = static_cast<const Point*>(convexGeom);
+ convexHullPts = pt->getCoordinatesRO();
+ break;
+ }
+ default:
+ throw util::IllegalArgumentException("computeConvex called with unsupported geometry type");
+ }
+
+ // special cases for lines or points or degenerate rings
+ switch (convexHullPts->size())
+ {
+ case 1:
+ return m_inputGeom->getFactory()->createPoint(convexHullPts->getAt<CoordinateXY>(0));
+ case 2:
+ case 3:
+ return computeMaximumLine(convexHullPts, m_inputGeom->getFactory());
+ default:
+ // TODO: ensure ring is CW
+ return computeConvexRing(convexHullPts);
+ }
+}
+
+
+/* private */
+std::unique_ptr<Polygon>
+MinimumAreaRectangle::computeConvexRing(const CoordinateSequence* ring)
+{
+ // Assert: ring is oriented CW
+
+ double minRectangleArea = DoubleMax;
+ std::size_t minRectangleBaseIndex = NO_COORD_INDEX;
+ std::size_t minRectangleDiamIndex = NO_COORD_INDEX;
+ std::size_t minRectangleLeftIndex = NO_COORD_INDEX;
+ std::size_t minRectangleRightIndex = NO_COORD_INDEX;
+
+ //-- start at vertex after first one
+ std::size_t diameterIndex = 1;
+ std::size_t leftSideIndex = 1;
+ std::size_t rightSideIndex = NO_COORD_INDEX; // initialized once first diameter is found
+
+ LineSegment segBase;
+ LineSegment segDiam;
+ // for each segment, find the next vertex which is at maximum distance
+ for (std::size_t i = 0; i < ring->size() - 1; i++) {
+ segBase.p0 = ring->getAt<CoordinateXY>(i);
+ segBase.p1 = ring->getAt<CoordinateXY>(i + 1);
+ diameterIndex = findFurthestVertex(ring, segBase, diameterIndex, 0);
+
+ const CoordinateXY& diamPt = ring->getAt<CoordinateXY>(diameterIndex);
+ CoordinateXY diamBasePt = segBase.project(diamPt);
+ segDiam.p0 = diamBasePt;
+ segDiam.p1 = diamPt;
+
+ leftSideIndex = findFurthestVertex(ring, segDiam, leftSideIndex, 1);
+
+ //-- init the max right index
+ if (i == 0) {
+ rightSideIndex = diameterIndex;
+ }
+ rightSideIndex = findFurthestVertex(ring, segDiam, rightSideIndex, -1);
+
+ double rectWidth = segDiam.distancePerpendicular(ring->getAt<CoordinateXY>(leftSideIndex))
+ + segDiam.distancePerpendicular(ring->getAt<CoordinateXY>(rightSideIndex));
+ double rectArea = segDiam.getLength() * rectWidth;
+
+ if (rectArea < minRectangleArea) {
+ minRectangleArea = rectArea;
+ minRectangleBaseIndex = i;
+ minRectangleDiamIndex = diameterIndex;
+ minRectangleLeftIndex = leftSideIndex;
+ minRectangleRightIndex = rightSideIndex;
+ }
+ }
+ return Rectangle::createFromSidePts(
+ ring->getAt<CoordinateXY>(minRectangleBaseIndex),
+ ring->getAt<CoordinateXY>(minRectangleBaseIndex + 1),
+ ring->getAt<CoordinateXY>(minRectangleDiamIndex),
+ ring->getAt<CoordinateXY>(minRectangleLeftIndex),
+ ring->getAt<CoordinateXY>(minRectangleRightIndex),
+ m_inputGeom->getFactory());
+}
+
+
+/* private */
+std::size_t
+MinimumAreaRectangle::findFurthestVertex(
+ const CoordinateSequence* pts,
+ const LineSegment& baseSeg,
+ std::size_t startIndex, int orient)
+{
+ double maxDistance = orientedDistance(baseSeg, pts->getAt<CoordinateXY>(startIndex), orient);
+ double nextDistance = maxDistance;
+ std::size_t maxIndex = startIndex;
+ std::size_t nextIndex = maxIndex;
+ //-- rotate "caliper" while distance from base segment is non-decreasing
+ while (isFurtherOrEqual(nextDistance, maxDistance, orient)) {
+ maxDistance = nextDistance;
+ maxIndex = nextIndex;
+
+ nextIndex = getNextIndex(pts, maxIndex);
+ if (nextIndex == startIndex)
+ break;
+ nextDistance = orientedDistance(baseSeg, pts->getAt<CoordinateXY>(nextIndex), orient);
+ }
+ return maxIndex;
+}
+
+
+/* private */
+bool
+MinimumAreaRectangle::isFurtherOrEqual(
+ double d1, double d2,
+ int orient)
+{
+ switch (orient) {
+ case 0: return std::abs(d1) >= std::abs(d2);
+ case 1: return d1 >= d2;
+ case -1: return d1 <= d2;
+ }
+ throw util::IllegalArgumentException("Invalid orientation index");
+}
+
+
+/* private static */
+double
+MinimumAreaRectangle::orientedDistance(
+ const LineSegment& seg,
+ const CoordinateXY& p,
+ int orient)
+{
+ double dist = seg.distancePerpendicularOriented(p);
+ if (orient == 0) {
+ return std::abs(dist);
+ }
+ return dist;
+}
+
+
+/* private static */
+std::size_t
+MinimumAreaRectangle::getNextIndex(
+ const CoordinateSequence* ring,
+ std::size_t index)
+{
+ if (index == NO_COORD_INDEX) index = 0;
+ else index = index + 1;
+
+ if (index >= ring->size() - 1)
+ index = 0;
+ return index;
+}
+
+
+/* private static */
+std::unique_ptr<LineString>
+MinimumAreaRectangle::computeMaximumLine(
+ const CoordinateSequence* pts,
+ const GeometryFactory* factory)
+{
+ //-- find max and min pts for X and Y
+ CoordinateXY ptMinX; ptMinX.setNull();
+ CoordinateXY ptMaxX; ptMaxX.setNull();
+ CoordinateXY ptMinY; ptMinY.setNull();
+ CoordinateXY ptMaxY; ptMaxY.setNull();
+ for (const CoordinateXY& p : pts->items<CoordinateXY>()) {
+ if (ptMinX.isNull() || p.x < ptMinX.x) ptMinX = p;
+ if (ptMaxX.isNull() || p.x > ptMaxX.x) ptMaxX = p;
+ if (ptMinY.isNull() || p.y < ptMinY.y) ptMinY = p;
+ if (ptMaxY.isNull() || p.y > ptMaxY.y) ptMaxY = p;
+ }
+
+ CoordinateXY p0 = ptMinX;
+ CoordinateXY p1 = ptMaxX;
+
+ //-- line is vertical - use Y pts
+ if (p0.x == p1.x) {
+ p0 = ptMinY;
+ p1 = ptMaxY;
+ }
+ CoordinateSequence cs({ p0, p1 });
+ return factory->createLineString(std::move(cs));
+}
+
+
+
+} // namespace geos.algorithm
+} // namespace geos
diff --git a/src/algorithm/MinimumDiameter.cpp b/src/algorithm/MinimumDiameter.cpp
index 879dedb79..813b161f9 100644
--- a/src/algorithm/MinimumDiameter.cpp
+++ b/src/algorithm/MinimumDiameter.cpp
@@ -3,6 +3,7 @@
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
+ * Copyright (C) 2023 Paul Ramsey <pramsey at cleverelephant.ca>
* Copyright (C) 2001-2002 Vivid Solutions Inc.
* Copyright (C) 2005 Refractions Research Inc.
*
@@ -11,10 +12,6 @@
* by the Free Software Foundation.
* See the COPYING file for more information.
*
- **********************************************************************
- *
- * Last port: algorithm/MinimumDiameter.java r966
- *
**********************************************************************/
#include <geos/constants.h>
@@ -59,7 +56,7 @@ namespace algorithm { // geos.algorithm
/**
* Compute a minimum diameter for a giver {@link Geometry}.
*
- * @param geom a Geometry
+ * @param newInputGeom a Geometry
*/
MinimumDiameter::MinimumDiameter(const Geometry* newInputGeom)
{
@@ -78,8 +75,8 @@ MinimumDiameter::MinimumDiameter(const Geometry* newInputGeom)
* (e.g. a convex Polygon or LinearRing,
* or a two-point LineString, or a Point).
*
- * @param geom a Geometry which is convex
- * @param isConvex <code>true</code> if the input geometry is convex
+ * @param newInputGeom a Geometry which is convex
+ * @param newIsConvex <code>true</code> if the input geometry is convex
*/
MinimumDiameter::MinimumDiameter(const Geometry* newInputGeom, const bool newIsConvex)
{
@@ -219,7 +216,7 @@ MinimumDiameter::computeConvexRingMinDiameter(const CoordinateSequence* pts)
unsigned int currMaxIndex = 1;
LineSegment seg;
- // compute the max distance for all segments in the ring, and pick the minimum
+ // for each segment, find a vertex at max distance, and pick the minimum
const std::size_t npts = pts->getSize();
for(std::size_t i = 1; i < npts; ++i) {
seg.p0 = pts->getAt(i - 1);
diff --git a/src/algorithm/Rectangle.cpp b/src/algorithm/Rectangle.cpp
new file mode 100644
index 000000000..7be35a475
--- /dev/null
+++ b/src/algorithm/Rectangle.cpp
@@ -0,0 +1,126 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2023 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/Rectangle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LinearRing.h>
+#include <geos/geom/Polygon.h>
+
+using geos::geom::CoordinateXY;
+using geos::geom::Coordinate;
+using geos::geom::CoordinateSequence;
+using geos::geom::GeometryFactory;
+using geos::geom::LineSegment;
+using geos::geom::LinearRing;
+using geos::geom::Polygon;
+
+
+namespace geos {
+namespace algorithm { // geos.algorithm
+
+
+/* public static */
+std::unique_ptr<Polygon>
+Rectangle::createFromSidePts(
+ const CoordinateXY& baseRightPt,
+ const CoordinateXY& baseLeftPt,
+ const CoordinateXY& oppositePt,
+ const CoordinateXY& leftSidePt,
+ const CoordinateXY& rightSidePt,
+ const GeometryFactory* factory)
+{
+ //-- deltas for the base segment provide slope
+ double dx = baseLeftPt.x - baseRightPt.x;
+ double dy = baseLeftPt.y - baseRightPt.y;
+ // Assert: dx and dy are not both zero
+
+ double baseC = computeLineEquationC(dx, dy, baseRightPt);
+ double oppC = computeLineEquationC(dx, dy, oppositePt);
+ double leftC = computeLineEquationC(-dy, dx, leftSidePt);
+ double rightC = computeLineEquationC(-dy, dx, rightSidePt);
+
+ //-- compute lines along edges of rectangle
+ LineSegment baseLine = createLineForStandardEquation(-dy, dx, baseC);
+ LineSegment oppLine = createLineForStandardEquation(-dy, dx, oppC);
+ LineSegment leftLine = createLineForStandardEquation(-dx, -dy, leftC);
+ LineSegment rightLine = createLineForStandardEquation(-dx, -dy, rightC);
+
+ /**
+ * Corners of rectangle are the intersections of the
+ * base and opposite, and left and right lines.
+ * The rectangle is constructed with CW orientation.
+ * The first side of the constructed rectangle contains the base segment.
+ *
+ * If a corner coincides with a input point
+ * the exact value is used to avoid numerical inaccuracy.
+ */
+ CoordinateXY p0 = rightSidePt.equals2D(baseRightPt) ? baseRightPt
+ : baseLine.lineIntersection(rightLine);
+ CoordinateXY p1 = leftSidePt.equals2D(baseLeftPt) ? baseLeftPt
+ : baseLine.lineIntersection(leftLine);
+ CoordinateXY p2 = leftSidePt.equals2D(oppositePt) ? oppositePt
+ : oppLine.lineIntersection(leftLine);
+ CoordinateXY p3 = rightSidePt.equals2D(oppositePt) ? oppositePt
+ : oppLine.lineIntersection(rightLine);
+
+ CoordinateSequence cs({ p0, p1, p2, p3, p0 });
+ return factory->createPolygon(std::move(cs));
+}
+
+
+
+/* private static */
+double
+Rectangle::computeLineEquationC(double a, double b, const CoordinateXY& p)
+{
+ return a * p.y - b * p.x;
+}
+
+
+/* private static */
+LineSegment
+Rectangle::createLineForStandardEquation(double a, double b, double c)
+{
+ Coordinate p0;
+ Coordinate p1;
+ /*
+ * Line equation is ax + by = c
+ * Slope m = -a/b.
+ * Y-intercept = c/b
+ * X-intercept = c/a
+ *
+ * If slope is low, use constant X values; if high use Y values.
+ * This handles lines that are vertical (b = 0, m = Inf )
+ * and horizontal (a = 0, m = 0).
+ */
+ if (std::abs(b) > std::abs(a)) {
+ //-- abs(m) < 1
+ p0 = Coordinate(0.0, c/b);
+ p1 = Coordinate(1.0, c/b - a/b);
+ }
+ else {
+ //-- abs(m) >= 1
+ p0 = Coordinate(c/a, 0.0);
+ p1 = Coordinate(c/a - b/a, 1.0);
+ }
+ return LineSegment(p0, p1);
+}
+
+
+
+} // namespace geos.algorithm
+} // namespace geos
+
diff --git a/src/geom/LineSegment.cpp b/src/geom/LineSegment.cpp
index 54b85c9ee..c352fc2cf 100644
--- a/src/geom/LineSegment.cpp
+++ b/src/geom/LineSegment.cpp
@@ -96,11 +96,26 @@ LineSegment::project(const Coordinate& p, Coordinate& ret) const
{
if(p == p0 || p == p1) {
ret = p;
+ return;
}
double r = projectionFactor(p);
ret = Coordinate(p0.x + r * (p1.x - p0.x), p0.y + r * (p1.y - p0.y));
}
+CoordinateXY
+LineSegment::project(const CoordinateXY& p) const
+{
+ if(p == p0 || p == p1) {
+ return p;
+ }
+ double r = projectionFactor(p);
+ double x = p0.x + r * (p1.x - p0.x);
+ double y = p0.y + r * (p1.y - p0.y);
+ return CoordinateXY(x, y);
+}
+
+
+
/*private*/
void
LineSegment::project(double factor, CoordinateXY& ret) const
@@ -304,6 +319,17 @@ LineSegment::offset(double offsetDistance)
return ls;
}
+/* public */
+double
+LineSegment::distancePerpendicularOriented(const CoordinateXY& p) const
+{
+ if (p0.equals2D(p1))
+ return p0.distance(p);
+ double dist = distancePerpendicular(p);
+ if (orientationIndex(p) < 0)
+ return -dist;
+ return dist;
+}
/* public */
std::unique_ptr<LineString>
diff --git a/tests/unit/algorithm/MinimumAreaRectangleTest.cpp b/tests/unit/algorithm/MinimumAreaRectangleTest.cpp
new file mode 100644
index 000000000..dad960837
--- /dev/null
+++ b/tests/unit/algorithm/MinimumAreaRectangleTest.cpp
@@ -0,0 +1,182 @@
+//
+// Test Suite for geos::algorithm::MinimumAreaRectangle
+
+#include <tut/tut.hpp>
+// geos
+#include <geos/algorithm/MinimumAreaRectangle.h>
+#include <geos/io/WKTReader.h>
+#include <geos/geom/Geometry.h>
+#include <utility.h>
+// std
+#include <sstream>
+#include <string>
+#include <memory>
+
+using geos::algorithm::MinimumAreaRectangle;
+using geos::io::WKTReader;
+using geos::geom::Geometry;
+
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_minimumarearectangle_data {
+
+ const double TOL = 1e-10;
+
+ WKTReader reader_;
+
+ test_minimumarearectangle_data() {};
+
+ void
+ checkMinRectangle(const std::string& wkt, const std::string& wktExpected)
+ {
+ std::unique_ptr<Geometry> geom = reader_.read(wkt);
+ std::unique_ptr<Geometry> actual = MinimumAreaRectangle::getMinimumRectangle(geom.get());
+ std::unique_ptr<Geometry> expected = reader_.read(wktExpected);
+ ensure_equals_geometry(expected.get(), actual.get(), TOL);
+ }
+
+};
+
+typedef test_group<test_minimumarearectangle_data> group;
+typedef group::object object;
+
+group test_minimumarearectangle_group("geos::algorithm::MinimumAreaRectangle");
+
+//
+// Test Cases
+//
+
+// testEmpty
+template<>
+template<>
+void object::test<1>()
+{
+ checkMinRectangle(
+ "POLYGON EMPTY",
+ "POLYGON EMPTY");
+}
+
+// testLineLengthZero
+template<>
+template<>
+void object::test<2>()
+{
+ checkMinRectangle(
+ "LINESTRING (1 1, 1 1)",
+ "POINT (1 1)");
+}
+
+// testLineHorizontal
+template<>
+template<>
+void object::test<3>()
+{
+ checkMinRectangle(
+ "LINESTRING (1 1, 3 1, 5 1, 7 1)", "LINESTRING (1 1, 7 1)");
+}
+
+// testLineVertical
+template<>
+template<>
+void object::test<4>()
+{
+ checkMinRectangle(
+ "LINESTRING (1 1, 1 4, 1 7, 1 9)",
+ "LINESTRING (1 1, 1 9)");
+ }
+
+// testLineObtuseAngle
+template<>
+template<>
+void object::test<5>()
+{
+ checkMinRectangle(
+ "LINESTRING (1 2, 3 8, 9 8)",
+ "POLYGON ((9 8, 1 2, -1.16 4.88, 6.84 10.88, 9 8))");
+}
+
+// testLineAcuteAngle
+template<>
+template<>
+void object::test<6>()
+{
+ checkMinRectangle(
+ "LINESTRING (5 2, 3 8, 9 8)",
+ "POLYGON ((5 2, 3 8, 8.4 9.8, 10.4 3.8, 5 2))");
+}
+
+// testNotMinDiameter
+template<>
+template<>
+void object::test<7>()
+{
+ checkMinRectangle(
+ "POLYGON ((150 300, 200 300, 300 300, 300 250, 280 120, 210 100, 100 100, 100 240, 150 300))",
+ "POLYGON ((100 100, 100 300, 300 300, 300 100, 100 100))");
+}
+
+// testTriangle
+template<>
+template<>
+void object::test<8>()
+{
+ checkMinRectangle(
+ "POLYGON ((100 100, 200 200, 160 240, 100 100))",
+ "POLYGON ((100 100, 160 240, 208.2758620689651 219.31034482758352, 148.2758620689666 79.31034482758756, 100 100))");
+}
+
+// testConvex
+template<>
+template<>
+void object::test<9>()
+{
+ checkMinRectangle("POLYGON ((3 8, 6 8, 9 5, 7 3, 3 1, 2 4, 3 8))",
+ "POLYGON ((0.2 6.6, 6.6 9.8, 9.4 4.2, 3 1, 0.2 6.6))");
+}
+
+/**
+* Failure case from https://trac.osgeo.org/postgis/ticket/5163
+* @throws Exception
+*/
+// testFlatDiagonal
+template<>
+template<>
+void object::test<10>()
+{
+ // bool error = false;
+ // try {
+ checkMinRectangle(
+ "LINESTRING(-99.48710639268086 34.79029839231914,-99.48370699999998 34.78689899963806,-99.48152167568102 34.784713675318976)",
+ "POLYGON ((-99.48710639268066 34.790298392318675, -99.48710639268066 34.790298392318675, -99.48152167568082 34.78471367531866, -99.48152167568082 34.78471367531866, -99.48710639268066 34.790298392318675))");
+ // } catch (std::exception& e) {
+ // error = true;
+ // // errorMessage = e.what();
+ // }
+ // ensure("caught exception", error == true);
+}
+
+// testBadRectl
+template<>
+template<>
+void object::test<11>()
+{
+ // bool error = false;
+ // try {
+ checkMinRectangle(
+ "POLYGON ((-5.21175 49.944633, -5.77435 50.021367, -5.7997 50.0306, -5.81815 50.0513, -5.82625 50.073567, -5.83085 50.1173, -6.2741 56.758767, -5.93245 57.909, -5.1158 58.644533, -5.07915 58.661733, -3.42575 58.686633, -3.1392 58.6685, -3.12495 58.666233, -1.88745 57.6444, 1.68845 52.715133, 1.7057 52.6829, 1.70915 52.6522, 1.7034 52.585433, 1.3867 51.214033, 1.36595 51.190267, 1.30485 51.121967, 0.96365 50.928567, 0.93025 50.912433, 0.1925 50.7436, -5.21175 49.944633))",
+ "POLYGON ((1.8583607388069103 50.41649058582797, -5.816631979932251 49.904263313964535, -6.395241388167441 58.57389735949991, 1.2797513305717105 59.08612463136336, 1.8583607388069103 50.41649058582797))");
+ // } catch (std::exception& e) {
+ // error = true;
+ // // errorMessage = e.what();
+ // }
+ // ensure("caught exception", error == true);
+}
+
+
+} // namespace tut
+
diff --git a/tests/unit/algorithm/RectangleTest.cpp b/tests/unit/algorithm/RectangleTest.cpp
new file mode 100644
index 000000000..c68810eb8
--- /dev/null
+++ b/tests/unit/algorithm/RectangleTest.cpp
@@ -0,0 +1,104 @@
+//
+// Test Suite for geos::algorithm::Rectangle
+
+#include <tut/tut.hpp>
+// geos
+#include <geos/algorithm/Rectangle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Geometry.h>
+#include <geos/io/WKTReader.h>
+#include <geos/constants.h>
+#include <utility.h>
+// std
+#include <string>
+#include <memory>
+
+using geos::algorithm::Rectangle;
+using geos::geom::Geometry;
+using geos::io::WKTReader;
+
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_rectangle_data {
+
+ const double TOL = 1e-10;
+
+ WKTReader reader_;
+
+ test_rectangle_data() {};
+
+ void
+ checkRectangle(const std::string& wkt, const std::string& wktExpected)
+ {
+
+ std::unique_ptr<Geometry> geom = reader_.read(wkt);
+ const LineString *line = static_cast<LineString*>(geom.get());
+
+ const Coordinate& baseRightPt = line->getCoordinateN(0);
+ const Coordinate& baseLeftPt = line->getCoordinateN(1);
+ const Coordinate& leftSidePt = line->getCoordinateN(2);
+ const Coordinate& oppositePt = line->getCoordinateN(3);
+ const Coordinate& rightSidePt = line->getCoordinateN(4);
+ std::unique_ptr<Polygon> actual = Rectangle::createFromSidePts(
+ baseRightPt, baseLeftPt,
+ oppositePt, leftSidePt,
+ rightSidePt, line->getFactory());
+ std::unique_ptr<Geometry> expected = reader_.read(wktExpected);
+ ensure_equals_geometry(
+ expected.get(),
+ static_cast<Geometry*>(actual.get()),
+ TOL);
+ }
+
+};
+
+typedef test_group<test_rectangle_data> group;
+typedef group::object object;
+
+group test_rectangle_group("geos::algorithm::Rectangle");
+
+//
+// Test Cases
+//
+
+// testOrthogonal()
+template<>
+template<>
+void object::test<1>()
+{
+ checkRectangle(
+ "LINESTRING (9 1, 1 1, 0 5, 7 10, 10 6)",
+ "POLYGON ((0 1, 0 10, 10 10, 10 1, 0 1))"
+ );
+}
+
+// test45()
+template<>
+template<>
+void object::test<2>()
+{
+ checkRectangle(
+ "LINESTRING (10 5, 5 0, 2 1, 2 7, 9 9)",
+ "POLYGON ((-1 4, 6.5 11.5, 11.5 6.5, 4 -1, -1 4))"
+ );
+}
+
+// testCoincidentBaseSides()
+template<>
+template<>
+void object::test<3>()
+{
+ checkRectangle(
+ "LINESTRING (10 5, 7 0, 7 0, 2 7, 10 5)",
+ "POLYGON ((0.2352941176470591 4.0588235294117645, 3.2352941176470598 9.058823529411764, 10 5, 7 0, 0.2352941176470591 4.0588235294117645))"
+ );
+}
+
+
+} // namespace tut
+
diff --git a/tests/unit/capi/GEOSMinimumRotatedRectangleTest.cpp b/tests/unit/capi/GEOSMinimumRotatedRectangleTest.cpp
index 89a1283b5..e8251f466 100644
--- a/tests/unit/capi/GEOSMinimumRotatedRectangleTest.cpp
+++ b/tests/unit/capi/GEOSMinimumRotatedRectangleTest.cpp
@@ -28,10 +28,12 @@ struct test_minimumrotatedrectangle_data : public capitest::utility {
// input
geom1_ = GEOSGeomFromWKT(wkt);
ensure(nullptr != geom1_);
+ GEOSSetSRID(geom1_, 1234);
// result
geom2_ = GEOSMinimumRotatedRectangle(geom1_);
ensure(nullptr != geom2_);
+ ensure("SRID equal", GEOSGetSRID(geom2_) == GEOSGetSRID(geom1_));
// expected
if (expected) {
@@ -93,7 +95,9 @@ template<>
void object::test<5>
()
{
- checkMinRectangle("LINESTRING (1 2, 3 8, 9 6)", "POLYGON ((9 6, 7 10, -1 6, 1 2, 9 6))");
+ checkMinRectangle(
+ "LINESTRING (1 2, 3 8, 9 6)",
+ "POLYGON ((1 2, 3 8, 9 6, 7 0, 1 2))");
}
// Failure case from https://trac.osgeo.org/postgis/ticket/5163
@@ -104,8 +108,20 @@ void object::test<6>
{
checkMinRectangle(
"LINESTRING(-99.48710639268086 34.79029839231914,-99.48370699999998 34.78689899963806,-99.48152167568102 34.784713675318976)",
- "LINESTRING (-99.48710639268086 34.79029839231914, -99.48152167568102 34.784713675318976)"
+ "POLYGON ((-99.48710639 34.79029839, -99.48710639 34.79029839, -99.48152168 34.78471368, -99.48152168 34.78471368, -99.48710639 34.79029839))"
);
}
+// Collection Input
+template<>
+template<>
+void object::test<7>
+()
+{
+ checkMinRectangle(
+ "MULTILINESTRING ((1 2, 3 8, 9 6))",
+ "POLYGON ((1 2, 3 8, 9 6, 7 0, 1 2))");
+}
+
+
} // namespace tut
diff --git a/tests/unit/geom/LineSegmentTest.cpp b/tests/unit/geom/LineSegmentTest.cpp
index c2212b324..31bb2dc7c 100644
--- a/tests/unit/geom/LineSegmentTest.cpp
+++ b/tests/unit/geom/LineSegmentTest.cpp
@@ -29,9 +29,11 @@ struct test_lineseg_data {
geos::geom::LineSegment v1;
double MAX_ABS_ERROR_INTERSECTION = 1e-5;
- void checkLineIntersection(double p1x, double p1y, double p2x, double p2y,
- double q1x, double q1y, double q2x, double q2y,
- double expectedx, double expectedy) {
+ void checkLineIntersection(
+ double p1x, double p1y, double p2x, double p2y,
+ double q1x, double q1y, double q2x, double q2y,
+ double expectedx, double expectedy)
+ {
LineSegment seg1(p1x, p1y, p2x, p2y);
LineSegment seg2(q1x, q1y, q2x, q2y);
@@ -93,6 +95,30 @@ struct test_lineseg_data {
ensure_equals(expectedOrient, orient);
}
+ void checkDistancePerpendicular(
+ double x0, double y0,
+ double x1, double y1,
+ double px, double py,
+ double expected)
+ {
+ LineSegment seg(x0, y0, x1, y1);
+ Coordinate c(px, py);
+ double dist = seg.distancePerpendicular(c);
+ ensure_equals("checkDistancePerpendicular", expected, dist, 0.000001);
+ }
+
+ void checkDistancePerpendicularOriented(
+ double x0, double y0,
+ double x1, double y1,
+ double px, double py,
+ double expected)
+ {
+ LineSegment seg(x0, y0, x1, y1);
+ Coordinate c(px, py);
+ double dist = seg.distancePerpendicularOriented(c);
+ ensure_equals("checkDistancePerpendicularOriented", expected, dist, 0.000001);
+ }
+
test_lineseg_data()
: ph1(0, 2), ph2(10, 2), pv1(0, 0), pv2(0, 10), h1(ph1, ph2), v1(pv1, pv2)
{}
@@ -275,4 +301,36 @@ void object::test<11>()
checkOrientationIndex(seg, 105, 105, 110, 100, -1);
}
+
+// testDistancePerpendicular
+template<>
+template<>
+void object::test<12>()
+{
+ checkDistancePerpendicular(1,1, 1,3, 2,4, 1);
+ checkDistancePerpendicular(1,1, 1,3, 0,4, 1);
+ checkDistancePerpendicular(1,1, 1,3, 1,4, 0);
+ checkDistancePerpendicular(1,1, 2,2, 4,4, 0);
+ //-- zero-length line segment
+ checkDistancePerpendicular(1,1, 1,1, 1,2, 1);
+}
+
+// testDistancePerpendicularOriented
+template<>
+template<>
+void object::test<13>()
+{
+ //-- right of line
+ checkDistancePerpendicularOriented(1,1, 1,3, 2,4, -1);
+ //-- left of line
+ checkDistancePerpendicularOriented(1,1, 1,3, 0,4, 1);
+ //-- on line
+ checkDistancePerpendicularOriented(1,1, 1,3, 1,4, 0);
+ checkDistancePerpendicularOriented(1,1, 2,2, 4,4, 0);
+ //-- zero-length segment
+ checkDistancePerpendicularOriented(1,1, 1,1, 1,2, 1);
+}
+
+
+
} // namespace tut
-----------------------------------------------------------------------
Summary of changes:
NEWS.md | 1 +
capi/geos_ts_c.cpp | 6 +-
include/geos/algorithm/Distance.h | 32 ++-
include/geos/algorithm/MinimumAreaRectangle.h | 159 ++++++++++++
include/geos/algorithm/MinimumDiameter.h | 30 ++-
include/geos/algorithm/Rectangle.h | 96 ++++++++
include/geos/geom/LineSegment.h | 40 ++-
include/geos/operation/buffer/OffsetCurve.h | 18 +-
src/algorithm/Distance.cpp | 23 ++
src/algorithm/MinimumAreaRectangle.cpp | 269 +++++++++++++++++++++
src/algorithm/MinimumDiameter.cpp | 13 +-
src/algorithm/Rectangle.cpp | 126 ++++++++++
src/geom/LineSegment.cpp | 26 ++
tests/unit/algorithm/MinimumAreaRectangleTest.cpp | 182 ++++++++++++++
tests/unit/algorithm/RectangleTest.cpp | 104 ++++++++
.../unit/capi/GEOSMinimumRotatedRectangleTest.cpp | 20 +-
tests/unit/capi/GEOSOffsetCurveTest.cpp | 4 +-
tests/unit/geom/LineSegmentTest.cpp | 64 ++++-
tests/unit/operation/buffer/OffsetCurveTest.cpp | 29 ++-
web/content/specifications/geojson.md | 2 +-
20 files changed, 1197 insertions(+), 47 deletions(-)
create mode 100644 include/geos/algorithm/MinimumAreaRectangle.h
create mode 100644 include/geos/algorithm/Rectangle.h
create mode 100644 src/algorithm/MinimumAreaRectangle.cpp
create mode 100644 src/algorithm/Rectangle.cpp
create mode 100644 tests/unit/algorithm/MinimumAreaRectangleTest.cpp
create mode 100644 tests/unit/algorithm/RectangleTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list