Voronoi APIs added to QuadEdgeSubdivision class, with test added

Contributed by Vishal Tiwari

git-svn-id: http://svn.osgeo.org/geos/trunk@3944 5242fede-7e19-0410-aef8-94bd7d2200fb

include/geos/triangulate/quadedge/QuadEdgeSubdivision.h

@@ -23,8 +23,8 @@
 #include <list>
 #include <stack>
 #include <set>
+#include <vector>
 
-#include <geos/geom/Envelope.h>
 #include <geos/geom/MultiLineString.h>
 #include <geos/triangulate/quadedge/QuadEdgeLocator.h>
 #include <geos/triangulate/quadedge/Vertex.h>
@@ -37,6 +37,8 @@ namespace geom {
 	class GeometryCollection;
 	class GeometryFactory;
 	class Coordinate;
+	class Geometry;
+	class Envelope;
 }
 
 namespace triangulate { //geos.triangulate
@@ -372,7 +374,8 @@ class GEOS_DLL QuadEdgeSubdivision {
 	void getTriangleCoordinates(TriList* triList, bool includeFrame);
 
 private:
-	class TriangleCoordinatesVisitor; 
+	class TriangleCoordinatesVisitor;
+	class TriangleCircumcentreVisitor;
 
 public:
 	/**
@@ -392,7 +395,61 @@ class GEOS_DLL QuadEdgeSubdivision {
 	 * @return a GeometryCollection of triangular Polygons. The caller takes ownership of the returned object.
 	 */
 	std::auto_ptr<geom::GeometryCollection> getTriangles(const geom::GeometryFactory &geomFact);
-
+
+	/**  
+	 * Gets the cells in the Voronoi diagram for this triangulation.
+	 * The cells are returned as a {@link GeometryCollection} of {@link Polygon}s
+	 * The userData of each polygon is set to be the {@link Coordinate}
+	 * of the cell site.  This allows easily associating external 
+	 * data associated with the sites to the cells.
+	 *
+	 * @param geomFact a geometry factory
+	 * @return a GeometryCollection of Polygons
+	 */
+	std::auto_ptr<geom::GeometryCollection> getVoronoiDiagram(const geom::GeometryFactory& geomFact);
+
+	/**  
+	 * Gets a List of {@link Polygon}s for the Voronoi cells 
+	 * of this triangulation.
+	 * The userData of each polygon is set to be the {@link Coordinate}
+	 * of the cell site.  This allows easily associating external 
+	 * data associated with the sites to the cells.
+	 *
+	 * @param geomFact a geometry factory
+	 * @return a List of Polygons
+	 */
+	std::auto_ptr< std::vector<geom::Geometry*> > getVoronoiCellPolygons(const geom::GeometryFactory& geomFact);
+
+	/**
+	 * Gets a collection of {@link QuadEdge}s whose origin
+	 * vertices are a unique set which includes
+	 * all vertices in the subdivision. 
+	 * The frame vertices can be included if required.
+	 * This is useful for algorithms which require traversing the 
+	 * subdivision starting at all vertices.
+	 * Returning a quadedge for each vertex
+	 * is more efficient than 
+	 * the alternative of finding the actual vertices
+	 * using {@link #getVertices} and then locating 
+	 * quadedges attached to them.
+	 *
+	 * @param includeFrame true if the frame vertices should be included
+	 * @return a collection of QuadEdge with the vertices of the subdivision as their origins
+	 */
+	std::auto_ptr<QuadEdgeSubdivision::QuadEdgeList> getVertexUniqueEdges(bool includeFrame);
+
+	/**
+	 * Gets the Voronoi cell around a site specified
+	 * by the origin of a QuadEdge.
+	 * The userData of the polygon is set to be the {@link Coordinate}
+	 * of the site.  This allows attaching external 
+	 * data associated with the site to this cell polygon.
+	 *
+	 * @param qe a quadedge originating at the cell site
+	 * @param geomFact a factory for building the polygon
+	 * @return a polygon indicating the cell extent
+	 */
+	std::auto_ptr<geom::Geometry> getVoronoiCellPolygon(QuadEdge* qe ,const geom::GeometryFactory& geomFact);
 };
 
 } //namespace geos.triangulate.quadedge

src/triangulate/quadedge/QuadEdgeSubdivision.cpp

@@ -18,6 +18,8 @@
 #include <geos/triangulate/quadedge/QuadEdgeSubdivision.h>
 
 #include <vector>
+#include <set>
+#include <iostream>
 
 #include <geos/geom/Polygon.h>
 #include <geos/geom/LineSegment.h>
@@ -26,6 +28,7 @@
 #include <geos/geom/CoordinateArraySequence.h>
 #include <geos/geom/CoordinateSequenceFactory.h>
 #include <geos/geom/CoordinateArraySequenceFactory.h>
+#include <geos/geom/CoordinateList.h>
 #include <geos/geom/GeometryCollection.h>
 #include <geos/geom/GeometryFactory.h>
 #include <geos/util/IllegalArgumentException.h>
@@ -35,9 +38,11 @@
 #include <geos/triangulate/quadedge/LastFoundQuadEdgeLocator.h>
 #include <geos/triangulate/quadedge/LocateFailureException.h>
 #include <geos/triangulate/quadedge/TriangleVisitor.h>
+#include <geos/geom/Triangle.h>
 
-using namespace geos::geom;
 
+using namespace geos::geom;
+using namespace std;
 namespace geos {
 namespace triangulate { //geos.triangulate
 namespace quadedge { //geos.triangulate.quadedge
@@ -389,6 +394,27 @@ QuadEdgeSubdivision::TriangleCoordinatesVisitor : public TriangleVisitor {
     }
 }; 
 
+
+class
+QuadEdgeSubdivision::TriangleCircumcentreVisitor : public TriangleVisitor
+{
+public:
+	void visit(QuadEdge* triEdges[3])
+	{
+		Triangle triangle(triEdges[0]->orig().getCoordinate(),
+				triEdges[1]->orig().getCoordinate(), triEdges[2]->orig().getCoordinate());
+		Coordinate cc;
+		triangle.circumcentre(cc);
+
+		Vertex ccVertex(cc);
+
+		for(int i=0 ; i<3 ; i++){
+			triEdges[i]->rot().setOrig(ccVertex);
+		}
+	}
+};
+
+
 void
 QuadEdgeSubdivision::getTriangleCoordinates(QuadEdgeSubdivision::TriList* triList, bool includeFrame)
 {
@@ -470,6 +496,102 @@ QuadEdgeSubdivision::getTriangles( const GeometryFactory &geomFact)
     return std::auto_ptr<GeometryCollection>(ret);
 }
 
+
+//Methods for VoronoiDiagram
+std::auto_ptr<geom::GeometryCollection> 
+QuadEdgeSubdivision::getVoronoiDiagram(const geom::GeometryFactory& geomFact)
+{
+	std::auto_ptr< std::vector<geom::Geometry*> > vorCells = getVoronoiCellPolygons(geomFact);
+	return std::auto_ptr<GeometryCollection>(geomFact.createGeometryCollection(vorCells.release()));
+}
+
+std::auto_ptr< std::vector<geom::Geometry*> >
+QuadEdgeSubdivision::getVoronoiCellPolygons(const geom::GeometryFactory& geomFact)
+{
+	std::auto_ptr< std::vector<geom::Geometry*> > cells(new std::vector<geom::Geometry*>);
+	TriangleCircumcentreVisitor* tricircumVisitor = new TriangleCircumcentreVisitor();
+	visitTriangles((TriangleVisitor*)tricircumVisitor, true);
+
+	std::auto_ptr<QuadEdgeSubdivision::QuadEdgeList> edges = getVertexUniqueEdges(false);
+
+	for(QuadEdgeSubdivision::QuadEdgeList::iterator it=edges->begin() ; it!=edges->end() ; ++it)
+	{
+		QuadEdge *qe = *it;
+		std::auto_ptr<geom::Geometry> poly = getVoronoiCellPolygon(qe,geomFact);
+
+		cells->push_back(poly.release());
+	}
+	delete tricircumVisitor;
+	return cells;
+}
+std::auto_ptr<geom::Geometry>
+QuadEdgeSubdivision::getVoronoiCellPolygon(QuadEdge* qe ,const geom::GeometryFactory& geomFact)
+{
+	std::vector<Coordinate> cellPts;
+	QuadEdge *startQE = qe;
+	do{
+		Coordinate cc = qe->rot().orig().getCoordinate();
+		cellPts.push_back(cc);
+		qe = &qe->oPrev();
+
+	}while ( qe != startQE);
+
+
+	//CoordList from a vector of Coordinates.
+	geom::CoordinateList coordList(cellPts);
+	//for checking close ring in CoordList class:
+	coordList.closeRing();
+
+	if(coordList.size() < 4)
+	{
+		cout << coordList << endl;
+		coordList.insert(coordList.end(),*(coordList.end()),true);
+	}
+
+	std::auto_ptr<Coordinate::Vect> pts = coordList.toCoordinateArray();
+	std::auto_ptr<geom::Geometry> cellPoly( 
+		geomFact.createPolygon(geomFact.createLinearRing(new geom::CoordinateArraySequence(pts.release())),NULL));
+
+	Vertex v = startQE->orig();
+	Coordinate c(0,0);
+	c = v.getCoordinate();
+	cellPoly->setUserData(reinterpret_cast<void*>(&c));
+	return cellPoly;
+}
+
+std::auto_ptr<QuadEdgeSubdivision::QuadEdgeList>
+QuadEdgeSubdivision::getVertexUniqueEdges(bool includeFrame)
+{
+	std::auto_ptr<QuadEdgeSubdivision::QuadEdgeList> edges(new QuadEdgeList());
+	std::set<Vertex> visitedVertices;
+	for(QuadEdgeSubdivision::QuadEdgeList::iterator it=quadEdges.begin() ; it!=quadEdges.end() ; ++it)
+	{
+		QuadEdge *qe = (QuadEdge*)(*it);
+		Vertex v = qe->orig();
+
+
+		if(visitedVertices.find(v) == visitedVertices.end())	//if v not found
+		{
+			visitedVertices.insert(v);
+			if(includeFrame || ! QuadEdgeSubdivision::isFrameVertex(v))
+			{
+				edges->push_back(qe);
+			}
+		}
+		QuadEdge *qd = &(qe->sym());
+		Vertex vd = qd->orig();
+
+
+		if(visitedVertices.find(vd) == visitedVertices.end()){
+			visitedVertices.insert(vd);
+			if(includeFrame || ! QuadEdgeSubdivision::isFrameVertex(vd)){
+				edges->push_back(qd);
+			}
+		}
+	}
+	return edges;
+}
+
 } //namespace geos.triangulate.quadedge
 } //namespace geos.triangulate
 } //namespace goes

tests/unit/triangulate/quadedge/QuadEdgeSubdivisionTest.cpp

@@ -8,15 +8,19 @@
 #include <geos/triangulate/quadedge/Vertex.h>
 #include <geos/triangulate/quadedge/QuadEdge.h>
 #include <geos/triangulate/quadedge/QuadEdgeSubdivision.h>
+#include <geos/triangulate/DelaunayTriangulationBuilder.h>
 #include <geos/geom/GeometryCollection.h>
 #include <geos/geom/GeometryFactory.h>
+#include <geos/geom/CoordinateSequence.h>
 //#include <geos/io/WKTWriter.h>
+#include <geos/io/WKTReader.h>
 #include <geos/geom/Envelope.h>
 #include <geos/geom/Coordinate.h>
 // std
 #include <stdio.h>
-
+#include <iostream>
 using namespace geos::triangulate::quadedge;
+using namespace geos::triangulate;
 using namespace geos::geom;
 using namespace geos::io;
 
@@ -66,6 +70,42 @@ namespace tut
 		//WKTWriter wkt;
 		//printf("%s\n", wkt.writeFormatted(tris).c_str());
 	}
+	template<>
+	template<>
+	void object::test<2>()
+	{
+		WKTReader reader;
+		Geometry* sites;
+		QuadEdgeSubdivision* subdiv;
+		sites = reader.read("MULTIPOINT ((150 200), (180 270), (275 163))");
+		CoordinateSequence* siteCoords = DelaunayTriangulationBuilder::extractUniqueCoordinates(*sites);
+		Envelope Env = DelaunayTriangulationBuilder::envelope(*siteCoords);
+
+		double expandBy = fmax(Env.getWidth() , Env.getHeight());
+		Env.expandBy(expandBy);
+
+		IncrementalDelaunayTriangulator::VertexList* vertices = DelaunayTriangulationBuilder::toVertices(*siteCoords);
+		subdiv = new quadedge::QuadEdgeSubdivision(Env,0);
+		IncrementalDelaunayTriangulator triangulator(subdiv);
+		triangulator.insertSites(*vertices);
+
+		//Test for getVoronoiDiagram::
+		GeometryFactory geomFact;
+		std::auto_ptr<GeometryCollection> polys = subdiv->getVoronoiDiagram(geomFact);
+		const char *expected_str = "GEOMETRYCOLLECTION (POLYGON ((-5849.974929324658 2268.0517257497568, -4529.9920486948895 2247.139449440667, 221.20588235294116 210.91176470588235, -684.4227119984187 -2848.644297291955, -5849.974929324658 2268.0517257497568)), POLYGON ((212.5 -3774.5, -684.4227119984187 -2848.644297291955, 221.20588235294116 210.91176470588235, 2448.7167655626645 2188.608343256571, 6235.0370264064295 2248.0370264064295, 212.5 -3774.5)), POLYGON ((-4529.9920486948895 2247.139449440667, 2448.7167655626645 2188.608343256571, 221.20588235294116 210.91176470588235, -4529.9920486948895 2247.139449440667)))";
+//		std::cout << polys->toString() << std::endl;
+
+		Geometry *expected = reader.read(expected_str);
+		polys->normalize();
+		expected->normalize();
+		ensure(polys->equalsExact(expected, 1e-7));
+		delete siteCoords;
+		delete sites;
+		delete subdiv;
+		delete vertices;
+		delete expected;
+//		ensure(polys->getCoordinateDimension() == expected->getCoordinateDimension());
+	}
 
 } // namespace tut