[geos-commits] [SCM] GEOS branch master updated. 02ddad983f2a72bba5c25bac69d8c5f0ce01b276
git at osgeo.org
git at osgeo.org
Thu Jan 24 14:43:34 PST 2019
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, master has been updated
via 02ddad983f2a72bba5c25bac69d8c5f0ce01b276 (commit)
via 5dee78d28de214fc46833615d069123b0ca8a393 (commit)
via 6283586e79e81682d754bc27d8d7d39e48e427c3 (commit)
from 397f07e0f819120b2c5b4de521c5ccd217224718 (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 02ddad983f2a72bba5c25bac69d8c5f0ce01b276
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Jan 24 14:42:56 2019 -0800
MinimumBoundingCircle passes tests
diff --git a/src/algorithm/MinimumBoundingCircle.cpp b/src/algorithm/MinimumBoundingCircle.cpp
index 233b019..f178ae0 100644
--- a/src/algorithm/MinimumBoundingCircle.cpp
+++ b/src/algorithm/MinimumBoundingCircle.cpp
@@ -173,7 +173,7 @@ MinimumBoundingCircle::computeCirclePoints()
return;
}
if (input->getNumPoints() == 1) {
- extremalPts[0] = *(input->getCoordinate());
+ extremalPts.push_back(*(input->getCoordinate()));
return;
}
@@ -305,140 +305,5 @@ MinimumBoundingCircle::pointWithMinAngleWithSegment(std::vector<Coordinate>& pts
}
-#if 0
-
-
-///////
-
-/*public*/
-InteriorPointArea::InteriorPointArea(const Geometry *g)
-{
- foundInterior=false;
- maxWidth=0.0;
- factory=g->getFactory();
- add(g);
-}
-
-/*public*/
-InteriorPointArea::~InteriorPointArea()
-{
-}
-
-/*public*/
-bool
-InteriorPointArea::getInteriorPoint(Coordinate& ret) const
-{
- if ( ! foundInterior ) return false;
-
- ret=interiorPoint;
- return true;
-}
-
-/*public*/
-void
-InteriorPointArea::add(const Geometry *geom)
-{
- const Polygon *poly = dynamic_cast<const Polygon*>(geom);
- if ( poly ) {
- addPolygon(geom);
- return;
- }
-
- const GeometryCollection *gc = dynamic_cast<const GeometryCollection*>(geom);
- if ( gc )
- {
- for(std::size_t i=0, n=gc->getNumGeometries(); i<n; i++) {
- add(gc->getGeometryN(i));
- }
- }
-}
-
-/*private*/
-void
-InteriorPointArea::addPolygon(const Geometry *geometry)
-{
- if (geometry->isEmpty()) return;
-
- Coordinate intPt;
- double width;
-
- unique_ptr<LineString> bisector ( horizontalBisector(geometry) );
- if ( bisector->getLength() == 0.0 ) {
- width = 0;
- intPt = bisector->getCoordinateN(0);
- }
- else {
- unique_ptr<Geometry> intersections ( bisector->intersection(geometry) );
- const Geometry *widestIntersection = widestGeometry(intersections.get());
- const Envelope *env = widestIntersection->getEnvelopeInternal();
- width=env->getWidth();
- env->centre(intPt);
- }
- if (!foundInterior || width>maxWidth) {
- interiorPoint = intPt;
- maxWidth = width;
- foundInterior=true;
- }
-}
-
-//@return if geometry is a collection, the widest sub-geometry; otherwise,
-//the geometry itself
-const Geometry*
-InteriorPointArea::widestGeometry(const Geometry *geometry)
-{
- const GeometryCollection *gc = dynamic_cast<const GeometryCollection*>(geometry);
- if ( gc ) {
- return widestGeometry(gc);
- } else {
- return geometry;
- }
-}
-
-const Geometry*
-InteriorPointArea::widestGeometry(const GeometryCollection* gc) {
- if (gc->isEmpty()) {
- return gc;
- }
- const Geometry* p_widestGeometry=gc->getGeometryN(0);
-
- // scan remaining geom components to see if any are wider
- for(std::size_t i=1, n=gc->getNumGeometries(); i<n; i++) // start at 1
- {
- const Envelope *env1(gc->getGeometryN(i)->getEnvelopeInternal());
- const Envelope *env2(p_widestGeometry->getEnvelopeInternal());
- if (env1->getWidth()>env2->getWidth()) {
- p_widestGeometry=gc->getGeometryN(i);
- }
- }
- return p_widestGeometry;
-}
-
-/* private */
-LineString*
-InteriorPointArea::horizontalBisector(const Geometry *geometry)
-{
- const Envelope *envelope=geometry->getEnvelopeInternal();
-
- /**
- * Original algorithm. Fails when geometry contains a horizontal
- * segment at the Y midpoint.
- */
- // Assert: for areas, minx <> maxx
- //double avgY=avg(envelope->getMinY(),envelope->getMaxY());
-
- double bisectY = SafeBisectorFinder::getBisectorY(*dynamic_cast<const Polygon *>(geometry));
- vector<Coordinate>*cv=new vector<Coordinate>(2);
- (*cv)[0].x = envelope->getMinX();
- (*cv)[0].y = bisectY;
- (*cv)[1].x = envelope->getMaxX();
- (*cv)[1].y = bisectY;
-
- CoordinateSequence *cl = factory->getCoordinateSequenceFactory()->create(cv);
-
- LineString *ret = factory->createLineString(cl);
- return ret;
-}
-#endif
-
} // namespace geos.algorithm
} // namespace geos
diff --git a/tests/unit/algorithm/MinimumBoundingCircleTest.cpp b/tests/unit/algorithm/MinimumBoundingCircleTest.cpp
index 34bbe96..bcc4185 100644
--- a/tests/unit/algorithm/MinimumBoundingCircleTest.cpp
+++ b/tests/unit/algorithm/MinimumBoundingCircleTest.cpp
@@ -61,16 +61,18 @@ namespace tut
std::unique_ptr<Geometry> actual(geomFact->createMultiPoint(exPts));
double actualRadius = mbc.getRadius();
geos::geom::Coordinate actualCentre = mbc.getCentre();
- //cout << "Centre = " << actualCentre << " Radius = " << actualRadius;
geomOut.reset(reader.read(wktOut));
- bool isEqual = geomOut->equals(geom.get());
+ bool isEqual = actual->equals(geomOut.get());
+
// need this hack because apparently equals does not work for MULTIPOINT EMPTY
if (geomOut->isEmpty() && geom->isEmpty())
isEqual = true;
- // if (!isEqual) {
- // System.out.println("Actual = " + actual + ", Expected = " + expected);
- // }
+ if (!isEqual) {
+ std::cout << "Centre = " << actualCentre << " Radius = " << actualRadius << " isEqual = " << isEqual << std::endl;
+ std::cout << "Actual:" << std::endl << actual->toString() << std::endl;
+ std::cout << "Expected:" << std::endl << geomOut->toString() << std::endl;
+ }
ensure(isEqual);
if (centreOut.isNull()) {
@@ -118,5 +120,79 @@ namespace tut
}
+ template<>
+ template<>
+ void object::test<2>()
+ {
+ Coordinate c(15, 15);
+ doMinimumBoundingCircleTest(
+ "MULTIPOINT ((10 10), (20 20))",
+ "MULTIPOINT ((10 10), (20 20))",
+ c,
+ 7.0710678118654755);
+ }
+
+ template<>
+ template<>
+ void object::test<3>()
+ {
+ Coordinate c(20, 20);
+ doMinimumBoundingCircleTest(
+ "MULTIPOINT ((10 10), (20 20), (30 30))",
+ "MULTIPOINT ((10 10), (30 30))",
+ c,
+ 14.142135623730951);
+ }
+
+ template<>
+ template<>
+ void object::test<4>()
+ {
+ Coordinate c(15, 15);
+ doMinimumBoundingCircleTest(
+ "MULTIPOINT ((10 10), (20 20), (10 20))",
+ "MULTIPOINT ((10 10), (20 20), (10 20))",
+ c,
+ 7.0710678118654755);
+ }
+
+ template<>
+ template<>
+ void object::test<5>()
+ {
+ Coordinate c(150, 100);
+ doMinimumBoundingCircleTest(
+ "POLYGON ((100 100, 200 100, 150 90, 100 100))",
+ "MULTIPOINT ((100 100), (200 100))",
+ c,
+ 50);
+ }
+
+ template<>
+ template<>
+ void object::test<6>()
+ {
+ Coordinate c(15, 15);
+ doMinimumBoundingCircleTest(
+ "MULTIPOINT ((10 10), (20 20), (10 20), (15 19))",
+ "MULTIPOINT ((10 10), (20 20), (10 20))",
+ c,
+ 7.0710678118654755);
+ }
+
+ template<>
+ template<>
+ void object::test<7>()
+ {
+ Coordinate c(26284.84180271327, 65267.114509082545);
+ doMinimumBoundingCircleTest(
+ "POLYGON ((26426 65078, 26531 65242, 26096 65427, 26075 65136, 26426 65078))",
+ "MULTIPOINT ((26531 65242), (26075 65136), (26096 65427))",
+ c,
+ 247.4360455914027);
+ }
+
+
+
} // namespace tut
commit 5dee78d28de214fc46833615d069123b0ca8a393
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Jan 24 14:17:48 2019 -0800
Draft tests for MinimumBoundingCircle
diff --git a/include/geos/algorithm/MinimumBoundingCircle.h b/include/geos/algorithm/MinimumBoundingCircle.h
index 104c16b..6ba76d2 100644
--- a/include/geos/algorithm/MinimumBoundingCircle.h
+++ b/include/geos/algorithm/MinimumBoundingCircle.h
@@ -68,7 +68,7 @@ public:
centre.setNull();
}
- ~MinimumBoundingCircle();
+ ~MinimumBoundingCircle() {};
/**
* Gets a geometry which represents the Minimum Bounding Circle.
diff --git a/src/algorithm/MinimumBoundingCircle.cpp b/src/algorithm/MinimumBoundingCircle.cpp
index 9c6316b..233b019 100644
--- a/src/algorithm/MinimumBoundingCircle.cpp
+++ b/src/algorithm/MinimumBoundingCircle.cpp
@@ -156,7 +156,7 @@ MinimumBoundingCircle::computeCentre()
void
MinimumBoundingCircle::compute()
{
- if (extremalPts.size() > 0) return;
+ if (!extremalPts.empty()) return;
computeCirclePoints();
computeCentre();
diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am
index 46dada0..186dd34 100644
--- a/tests/unit/Makefile.am
+++ b/tests/unit/Makefile.am
@@ -51,6 +51,7 @@ geos_unit_SOURCES = \
algorithm/InteriorPointAreaTest.cpp \
algorithm/LengthTest.cpp \
algorithm/LocatePointInRingTest.cpp \
+ algorithm/MinimumBoundingCircleTest.cpp \
algorithm/MinimumDiameterTest.cpp \
algorithm/OrientationIndexFailureTest.cpp \
algorithm/PointLocatorTest.cpp \
diff --git a/tests/unit/algorithm/MinimumBoundingCircleTest.cpp b/tests/unit/algorithm/MinimumBoundingCircleTest.cpp
new file mode 100644
index 0000000..34bbe96
--- /dev/null
+++ b/tests/unit/algorithm/MinimumBoundingCircleTest.cpp
@@ -0,0 +1,122 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2011 Sandro Santilli <strk at kbt.io>
+ *
+ * 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.
+ *
+ **********************************************************************/
+
+//
+// Test Suite for geos::algorithm::MinimumBoundingCircle
+
+#include <tut/tut.hpp>
+// geos
+#include <geos/algorithm/MinimumBoundingCircle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/PrecisionModel.h>
+#include <geos/io/WKTReader.h>
+
+// std
+#include <sstream>
+#include <string>
+#include <memory>
+
+namespace tut
+{
+ //
+ // Test Group
+ //
+
+ // dummy data, not used
+ struct test_minimumboundingcircle_data {
+ typedef geos::geom::Geometry Geometry;
+ typedef geos::geom::Coordinate Coordinate;
+ typedef geos::geom::PrecisionModel PrecisionModel;
+ typedef geos::geom::GeometryFactory GeometryFactory;
+ typedef geos::algorithm::MinimumBoundingCircle MinimumBoundingCircle;
+
+ geos::io::WKTReader reader;
+ std::unique_ptr<Geometry> geom;
+ std::unique_ptr<Geometry> geomOut;
+ GeometryFactory::Ptr geomFact = GeometryFactory::create();
+
+ test_minimumboundingcircle_data()
+ {}
+
+ void doMinimumBoundingCircleTest(std::string wktIn, std::string wktOut,
+ geos::geom::Coordinate ¢reOut, double radiusOut)
+ {
+
+ geom.reset(reader.read(wktIn));
+ MinimumBoundingCircle mbc(geom.get());
+ std::vector<Coordinate> exPts = mbc.getExtremalPoints();
+ std::unique_ptr<Geometry> actual(geomFact->createMultiPoint(exPts));
+ double actualRadius = mbc.getRadius();
+ geos::geom::Coordinate actualCentre = mbc.getCentre();
+ //cout << "Centre = " << actualCentre << " Radius = " << actualRadius;
+
+ geomOut.reset(reader.read(wktOut));
+ bool isEqual = geomOut->equals(geom.get());
+ // need this hack because apparently equals does not work for MULTIPOINT EMPTY
+ if (geomOut->isEmpty() && geom->isEmpty())
+ isEqual = true;
+ // if (!isEqual) {
+ // System.out.println("Actual = " + actual + ", Expected = " + expected);
+ // }
+ ensure(isEqual);
+
+ if (centreOut.isNull()) {
+ ensure(centreOut.distance(actualCentre) < 0.0001);
+ }
+ if (radiusOut >= 0) {
+ ensure(fabs(radiusOut - actualRadius) < 0.0001);
+ }
+ }
+
+ void doMinimumBoundingCircleTest(std::string wktIn, std::string wktOut)
+ {
+ geos::geom::Coordinate c;
+ c.setNull();
+ doMinimumBoundingCircleTest(wktIn, wktOut, c, -1.0);
+ }
+
+ };
+
+ typedef test_group<test_minimumboundingcircle_data> group;
+ typedef group::object object;
+
+ group test_minimumboundingcircle_group("geos::algorithm::MinimumBoundingCircle");
+
+
+
+
+
+
+ //
+ // Test Cases
+ //
+
+
+ template<>
+ template<>
+ void object::test<1>()
+ {
+ Coordinate c(10, 10);
+ doMinimumBoundingCircleTest(
+ "POINT (10 10)",
+ "POINT (10 10)",
+ c,
+ 0);
+
+ }
+
+} // namespace tut
+
commit 6283586e79e81682d754bc27d8d7d39e48e427c3
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Thu Jan 24 11:56:34 2019 -0800
First draft of MinimumBoundingCircle port
diff --git a/include/geos/algorithm/MinimumBoundingCircle.h b/include/geos/algorithm/MinimumBoundingCircle.h
new file mode 100644
index 0000000..104c16b
--- /dev/null
+++ b/include/geos/algorithm/MinimumBoundingCircle.h
@@ -0,0 +1,139 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2019 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.
+ *
+ **********************************************************************
+ *
+ * Last port: algorithm/MinimumBoundingCircle.java 2019-01-23
+ *
+ **********************************************************************/
+
+#ifndef GEOS_ALGORITHM_MINIMUMBOUNDINGCIRCLE_H
+#define GEOS_ALGORITHM_MINIMUMBOUNDINGCIRCLE_H
+
+#include <geos/export.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/Point.h>
+#include <geos/geom/Triangle.h>
+
+#include <vector>
+
+// Forward declarations
+// namespace geos {
+// namespace geom {
+// class GeometryCollection;
+// }
+// }
+
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+
+class GEOS_DLL MinimumBoundingCircle {
+
+private:
+
+ // member variables
+ const geom::Geometry *input;
+ std::vector<geom::Coordinate> extremalPts;
+ geom::Coordinate centre;
+ double radius;
+
+ void computeCentre();
+ void compute();
+ void computeCirclePoints();
+ geom::Coordinate lowestPoint(std::vector<geom::Coordinate>& pts);
+ geom::Coordinate pointWitMinAngleWithX(std::vector<geom::Coordinate>& pts, geom::Coordinate& P);
+ geom::Coordinate pointWithMinAngleWithSegment(std::vector<geom::Coordinate>& pts,
+ geom::Coordinate& P, geom::Coordinate& Q);
+
+
+public:
+
+ MinimumBoundingCircle(const geom::Geometry *geom):
+ input(nullptr),
+ radius(0.0)
+ {
+ input = geom;
+ centre.setNull();
+ }
+
+ ~MinimumBoundingCircle();
+
+ /**
+ * Gets a geometry which represents the Minimum Bounding Circle.
+ * If the input is degenerate (empty or a single unique point),
+ * this method will return an empty geometry or a single Point geometry.
+ * Otherwise, a Polygon will be returned which approximates the
+ * Minimum Bounding Circle.
+ * (Note that because the computed polygon is only an approximation,
+ * it may not precisely contain all the input points.)
+ *
+ * @return a Geometry representing the Minimum Bounding Circle.
+ */
+ geom::Geometry* getCircle();
+
+ /**
+ * Gets a geometry representing a line between the two farthest points
+ * in the input.
+ * These points will be two of the extremal points of the Minimum Bounding Circle.
+ * They also lie on the convex hull of the input.
+ *
+ * @return a LineString between the two farthest points of the input
+ * @return a empty LineString if the input is empty
+ * @return a Point if the input is a point
+ */
+ geom::Geometry* getFarthestPoints();
+
+ /**
+ * Gets a geometry representing the diameter of the computed Minimum Bounding
+ * Circle.
+ *
+ * @return the diameter LineString of the Minimum Bounding Circle
+ * @return a empty LineString if the input is empty
+ * @return a Point if the input is a point
+ */
+ geom::Geometry* getDiameter();
+
+ /**
+ * Gets the extremal points which define the computed Minimum Bounding Circle.
+ * There may be zero, one, two or three of these points,
+ * depending on the number of points in the input
+ * and the geometry of those points.
+ *
+ * @return the points defining the Minimum Bounding Circle
+ */
+ std::vector<geom::Coordinate> getExtremalPoints();
+
+ /**
+ * Gets the centre point of the computed Minimum Bounding Circle.
+ *
+ * @return the centre point of the Minimum Bounding Circle
+ * @return null if the input is empty
+ */
+ geom::Coordinate getCentre();
+
+ /**
+ * Gets the radius of the computed Minimum Bounding Circle.
+ *
+ * @return the radius of the Minimum Bounding Circle
+ */
+ double getRadius();
+
+};
+
+} // namespace geos::algorithm
+} // namespace geos
+
+#endif // GEOS_ALGORITHM_MINIMUMBOUNDINGCIRCLE_H
+
diff --git a/include/geos/geom/Triangle.h b/include/geos/geom/Triangle.h
index e80bb95..1a78ec8 100644
--- a/include/geos/geom/Triangle.h
+++ b/include/geos/geom/Triangle.h
@@ -66,6 +66,8 @@ public:
*/
void circumcentre(Coordinate& resultPoint);
+ static const Coordinate circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2);
+
private:
/**
diff --git a/src/algorithm/Makefile.am b/src/algorithm/Makefile.am
index f48c0bf..2579281 100644
--- a/src/algorithm/Makefile.am
+++ b/src/algorithm/Makefile.am
@@ -23,6 +23,7 @@ libalgorithm_la_SOURCES = \
InteriorPointPoint.cpp \
LineIntersector.cpp \
Length.cpp \
+ MinimumBoundingCircle.cpp \
MinimumDiameter.cpp \
NotRepresentableException.cpp \
Orientation.cpp \
diff --git a/src/algorithm/MinimumBoundingCircle.cpp b/src/algorithm/MinimumBoundingCircle.cpp
new file mode 100644
index 0000000..9c6316b
--- /dev/null
+++ b/src/algorithm/MinimumBoundingCircle.cpp
@@ -0,0 +1,444 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2013 Sandro Santilli <strk at kbt.io>
+ * Copyright (C) 2005-2006 Refractions Research Inc.
+ * Copyright (C) 2001-2002 Vivid Solutions Inc.
+ *
+ * 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: algorithm/MinimumBoundingCircle.java 2019-01-24
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/Angle.h>
+#include <geos/algorithm/MinimumBoundingCircle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequenceFactory.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryCollection.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/Polygon.h>
+#include <geos/geom/Triangle.h>
+#include <geos/util/GEOSException.h>
+
+#include <math.h> // sqrt
+#include <memory> // for unique_ptr
+#include <typeinfo>
+#include <vector>
+
+using namespace geos::geom;
+
+namespace geos {
+namespace algorithm { // geos.algorithm
+
+/*public*/
+Geometry*
+MinimumBoundingCircle::getCircle()
+{
+ //TODO: ensure the output circle contains the extermal points.
+ //TODO: or maybe even ensure that the returned geometry contains ALL the input points?
+
+ compute();
+ if (centre.isNull())
+ return input->getFactory()->createPolygon();
+ Point *centrePoint = input->getFactory()->createPoint(centre);
+ if (radius == 0.0)
+ return centrePoint;
+ return centrePoint->buffer(radius);
+}
+
+/*public*/
+Geometry*
+MinimumBoundingCircle::getFarthestPoints()
+{
+ compute();
+ switch (extremalPts.size()) {
+ case 0:
+ return input->getFactory()->createLineString();
+ case 1:
+ return input->getFactory()->createPoint(centre);
+ }
+
+ size_t dims = input->getDimension();
+ size_t len = 2;
+ CoordinateSequence *cs = input->getFactory()->getCoordinateSequenceFactory()->create(len, dims);
+ cs->add(extremalPts[0], true);
+ cs->add(extremalPts[extremalPts.size() - 1], true);
+ return input->getFactory()->createLineString(cs);
+}
+
+/*public*/
+Geometry*
+MinimumBoundingCircle::getDiameter()
+{
+ compute();
+ switch (extremalPts.size()) {
+ case 0:
+ return input->getFactory()->createLineString();
+ case 1:
+ return input->getFactory()->createPoint(centre);
+ }
+ size_t dims = input->getDimension();
+ size_t len = 2;
+ CoordinateSequence *cs = input->getFactory()->getCoordinateSequenceFactory()->create(len, dims);
+ // TODO: handle case of 3 extremal points, by computing a line from one of
+ // them through the centre point with len = 2*radius
+ cs->add(extremalPts[0], true);
+ cs->add(extremalPts[1], true);
+ return input->getFactory()->createLineString(cs);
+}
+
+/*public*/
+std::vector<Coordinate>
+MinimumBoundingCircle::getExtremalPoints()
+{
+ compute();
+ return extremalPts;
+}
+
+/*public*/
+Coordinate
+MinimumBoundingCircle::getCentre()
+{
+ compute();
+ return centre;
+}
+
+/*public*/
+double
+MinimumBoundingCircle::getRadius()
+{
+ compute();
+ return radius;
+}
+
+/*private*/
+void
+MinimumBoundingCircle::computeCentre()
+{
+ switch (extremalPts.size()) {
+ case 0: {
+ centre.setNull();
+ break;
+ }
+ case 1: {
+ centre = extremalPts[0];
+ break;
+ }
+ case 2: {
+ double xAvg = (extremalPts[0].x + extremalPts[1].x) / 2.0;
+ double yAvg = (extremalPts[0].y + extremalPts[1].y) / 2.0;
+ Coordinate c(xAvg, yAvg);
+ centre = c;
+ break;
+ }
+ case 3: {
+ centre = Triangle::circumcentre(extremalPts[0], extremalPts[1], extremalPts[2]);
+ break;
+ }
+ default: {
+ util::GEOSException("Logic failure in MinimumBoundingCircle algorithm!");\
+ }
+ }
+}
+
+/*private*/
+void
+MinimumBoundingCircle::compute()
+{
+ if (extremalPts.size() > 0) return;
+
+ computeCirclePoints();
+ computeCentre();
+ if (!centre.isNull())
+ radius = centre.distance(extremalPts[0]);
+}
+
+/*private*/
+void
+MinimumBoundingCircle::computeCirclePoints()
+{
+ // handle degenerate or trivial cases
+ if (input->isEmpty()) {
+ return;
+ }
+ if (input->getNumPoints() == 1) {
+ extremalPts[0] = *(input->getCoordinate());
+ return;
+ }
+
+ /**
+ * The problem is simplified by reducing to the convex hull.
+ * Computing the convex hull also has the useful effect of eliminating duplicate points
+ */
+ Geometry *convexHull = input->convexHull();
+
+ CoordinateSequence *cs = convexHull->getCoordinates();
+ std::vector<Coordinate> pts;
+ cs->toVector(pts);
+
+ // strip duplicate final point, if any
+ if (pts.front().equals2D(pts.back())) {
+ pts.pop_back();
+ }
+
+ /**
+ * Optimization for the trivial case where the CH has fewer than 3 points
+ */
+ if (pts.size() <= 2) {
+ extremalPts = pts;
+ return;
+ }
+
+ // find a point P with minimum Y ordinate
+ Coordinate P = lowestPoint(pts);
+
+ // find a point Q such that the angle that PQ makes with the x-axis is minimal
+ Coordinate Q = pointWitMinAngleWithX(pts, P);
+
+ /**
+ * Iterate over the remaining points to find
+ * a pair or triplet of points which determine the minimal circle.
+ * By the design of the algorithm,
+ * at most <tt>pts.length</tt> iterations are required to terminate
+ * with a correct result.
+ */
+ for (auto p: pts) {
+ Coordinate R = pointWithMinAngleWithSegment(pts, P, Q);
+
+ // if PRQ is obtuse, then MBC is determined by P and Q
+ if (algorithm::Angle::isObtuse(P, R, Q)) {
+ extremalPts.push_back(P);
+ extremalPts.push_back(Q);
+ return;
+ }
+ // if RPQ is obtuse, update baseline and iterate
+ if (algorithm::Angle::isObtuse(R, P, Q)) {
+ P = R;
+ continue;
+ }
+ // if RQP is obtuse, update baseline and iterate
+ if (algorithm::Angle::isObtuse(R, Q, P)) {
+ Q = R;
+ continue;
+ }
+ // otherwise all angles are acute, and the MBC is determined by the triangle PQR
+ extremalPts.push_back(P);
+ extremalPts.push_back(Q);
+ extremalPts.push_back(R);
+ return;
+ }
+ // never get here
+ util::GEOSException("Logic failure in MinimumBoundingCircle algorithm!");
+}
+
+/*private*/
+Coordinate
+MinimumBoundingCircle::lowestPoint(std::vector<Coordinate>& pts)
+{
+ Coordinate min = pts[0];
+ for (auto pt: pts) {
+ if (pt.y < min.y)
+ min = pt;
+ }
+ return min;
+}
+
+
+/*private*/
+Coordinate
+MinimumBoundingCircle::pointWitMinAngleWithX(std::vector<Coordinate>& pts, Coordinate& P)
+{
+ double minSin = std::numeric_limits<double>::max();
+ Coordinate minAngPt;
+ minAngPt.setNull();
+ for (auto p: pts) {
+
+ if (p == P) continue;
+
+ /**
+ * The sin of the angle is a simpler proxy for the angle itself
+ */
+ double dx = p.x - P.x;
+ double dy = p.y - P.y;
+ if (dy < 0) dy = -dy;
+ double len = sqrt(dx * dx + dy * dy);
+ double sin = dy / len;
+
+ if (sin < minSin) {
+ minSin = sin;
+ minAngPt = p;
+ }
+ }
+ return minAngPt;
+}
+
+
+/*private*/
+Coordinate
+MinimumBoundingCircle::pointWithMinAngleWithSegment(std::vector<Coordinate>& pts, Coordinate& P, Coordinate& Q)
+{
+ double minAng = std::numeric_limits<double>::max();
+ Coordinate minAngPt;
+ minAngPt.setNull();
+ for (auto p: pts) {
+ if (p == P) continue;
+ if (p == Q) continue;
+
+ double ang = Angle::angleBetween(P, p, Q);
+ if (ang < minAng) {
+ minAng = ang;
+ minAngPt = p;
+ }
+ }
+ return minAngPt;
+}
+
+
+#if 0
+
+
+///////
+
+/*public*/
+InteriorPointArea::InteriorPointArea(const Geometry *g)
+{
+ foundInterior=false;
+ maxWidth=0.0;
+ factory=g->getFactory();
+ add(g);
+}
+
+/*public*/
+InteriorPointArea::~InteriorPointArea()
+{
+}
+
+/*public*/
+bool
+InteriorPointArea::getInteriorPoint(Coordinate& ret) const
+{
+ if ( ! foundInterior ) return false;
+
+ ret=interiorPoint;
+ return true;
+}
+
+/*public*/
+void
+InteriorPointArea::add(const Geometry *geom)
+{
+ const Polygon *poly = dynamic_cast<const Polygon*>(geom);
+ if ( poly ) {
+ addPolygon(geom);
+ return;
+ }
+
+ const GeometryCollection *gc = dynamic_cast<const GeometryCollection*>(geom);
+ if ( gc )
+ {
+ for(std::size_t i=0, n=gc->getNumGeometries(); i<n; i++) {
+ add(gc->getGeometryN(i));
+ }
+ }
+}
+
+/*private*/
+void
+InteriorPointArea::addPolygon(const Geometry *geometry)
+{
+ if (geometry->isEmpty()) return;
+
+ Coordinate intPt;
+ double width;
+
+ unique_ptr<LineString> bisector ( horizontalBisector(geometry) );
+ if ( bisector->getLength() == 0.0 ) {
+ width = 0;
+ intPt = bisector->getCoordinateN(0);
+ }
+ else {
+ unique_ptr<Geometry> intersections ( bisector->intersection(geometry) );
+ const Geometry *widestIntersection = widestGeometry(intersections.get());
+ const Envelope *env = widestIntersection->getEnvelopeInternal();
+ width=env->getWidth();
+ env->centre(intPt);
+ }
+ if (!foundInterior || width>maxWidth) {
+ interiorPoint = intPt;
+ maxWidth = width;
+ foundInterior=true;
+ }
+}
+
+//@return if geometry is a collection, the widest sub-geometry; otherwise,
+//the geometry itself
+const Geometry*
+InteriorPointArea::widestGeometry(const Geometry *geometry)
+{
+ const GeometryCollection *gc = dynamic_cast<const GeometryCollection*>(geometry);
+ if ( gc ) {
+ return widestGeometry(gc);
+ } else {
+ return geometry;
+ }
+}
+
+const Geometry*
+InteriorPointArea::widestGeometry(const GeometryCollection* gc) {
+ if (gc->isEmpty()) {
+ return gc;
+ }
+ const Geometry* p_widestGeometry=gc->getGeometryN(0);
+
+ // scan remaining geom components to see if any are wider
+ for(std::size_t i=1, n=gc->getNumGeometries(); i<n; i++) // start at 1
+ {
+ const Envelope *env1(gc->getGeometryN(i)->getEnvelopeInternal());
+ const Envelope *env2(p_widestGeometry->getEnvelopeInternal());
+ if (env1->getWidth()>env2->getWidth()) {
+ p_widestGeometry=gc->getGeometryN(i);
+ }
+ }
+ return p_widestGeometry;
+}
+
+/* private */
+LineString*
+InteriorPointArea::horizontalBisector(const Geometry *geometry)
+{
+ const Envelope *envelope=geometry->getEnvelopeInternal();
+
+ /**
+ * Original algorithm. Fails when geometry contains a horizontal
+ * segment at the Y midpoint.
+ */
+ // Assert: for areas, minx <> maxx
+ //double avgY=avg(envelope->getMinY(),envelope->getMaxY());
+
+ double bisectY = SafeBisectorFinder::getBisectorY(*dynamic_cast<const Polygon *>(geometry));
+ vector<Coordinate>*cv=new vector<Coordinate>(2);
+ (*cv)[0].x = envelope->getMinX();
+ (*cv)[0].y = bisectY;
+ (*cv)[1].x = envelope->getMaxX();
+ (*cv)[1].y = bisectY;
+
+ CoordinateSequence *cl = factory->getCoordinateSequenceFactory()->create(cv);
+
+ LineString *ret = factory->createLineString(cl);
+ return ret;
+}
+#endif
+
+} // namespace geos.algorithm
+} // namespace geos
diff --git a/src/geom/Triangle.cpp b/src/geom/Triangle.cpp
index fda964c..4a9685c 100644
--- a/src/geom/Triangle.cpp
+++ b/src/geom/Triangle.cpp
@@ -52,6 +52,17 @@ void Triangle::circumcentre(Coordinate& result)
result = Coordinate(ccx,ccy);
}
+/* public static */
+const Coordinate
+Triangle::circumcentre(const Coordinate& p0, const Coordinate& p1, const Coordinate& p2)
+{
+ Triangle t(p0, p1, p2);
+ Coordinate c;
+ t.circumcentre(c);
+ return c;
+}
+
+
double Triangle::det(double m00 , double m01 , double m10 , double m11) const
{
return m00 * m11 - m01 * m10;
-----------------------------------------------------------------------
Summary of changes:
include/geos/algorithm/MinimumBoundingCircle.h | 139 +++++++++
include/geos/geom/Triangle.h | 2 +
src/algorithm/Makefile.am | 1 +
src/algorithm/MinimumBoundingCircle.cpp | 309 +++++++++++++++++++++
src/geom/Triangle.cpp | 11 +
tests/unit/Makefile.am | 1 +
tests/unit/algorithm/MinimumBoundingCircleTest.cpp | 198 +++++++++++++
7 files changed, 661 insertions(+)
create mode 100644 include/geos/algorithm/MinimumBoundingCircle.h
create mode 100644 src/algorithm/MinimumBoundingCircle.cpp
create mode 100644 tests/unit/algorithm/MinimumBoundingCircleTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list