[geos-commits] r3677 - in trunk: include/geos/algorithm src/algorithm tests/unit tests/unit/capi

svn_geos at osgeo.org svn_geos at osgeo.org
Fri Jun 22 08:03:44 PDT 2012


Author: strk
Date: 2012-06-22 08:03:33 -0700 (Fri, 22 Jun 2012)
New Revision: 3677

Added:
   trunk/tests/unit/capi/GEOSGetCentroidTest.cpp
Modified:
   trunk/include/geos/algorithm/CentroidArea.h
   trunk/src/algorithm/CentroidArea.cpp
   trunk/tests/unit/Makefile.am
Log:
Port robustness fix to CentroidArea (#559)

Modified: trunk/include/geos/algorithm/CentroidArea.h
===================================================================
--- trunk/include/geos/algorithm/CentroidArea.h	2012-06-22 15:02:06 UTC (rev 3676)
+++ trunk/include/geos/algorithm/CentroidArea.h	2012-06-22 15:03:33 UTC (rev 3677)
@@ -11,6 +11,10 @@
  * by the Free Software Foundation. 
  * See the COPYING file for more information.
  *
+ **********************************************************************
+ *
+ * Last port: algorithm/CentroidArea.java r612
+ *
  **********************************************************************/
 
 #ifndef GEOS_ALGORITHM_CENTROIDAREA_H
@@ -53,7 +57,8 @@
 	CentroidArea()
 		:
 		basePt(0.0, 0.0),
-		areasum2(0)
+		areasum2(0.0),
+		totalLength(0.0)
 	{}
 
 	~CentroidArea() {}
@@ -74,9 +79,10 @@
 	 */
 	void add(const geom::CoordinateSequence *ring);
 
+	// TODO: deprecate
 	geom::Coordinate* getCentroid() const;
 
-	/// Return false if a centroid couldn't be computed
+	/// Return false if a centroid couldn't be computed ( empty polygon )
 	bool getCentroid(geom::Coordinate& ret) const;
 
 private:
@@ -93,6 +99,10 @@
 	/// partial centroid sum
 	geom::Coordinate cg3;
 
+	// data for linear centroid computation, if needed
+	geom::Coordinate centSum;
+	double totalLength;
+
 	void setBasePoint(const geom::Coordinate &newbasePt);
 
 	void add(const geom::Polygon *poly);
@@ -110,6 +120,17 @@
 	static double area2(const geom::Coordinate &p1, const geom::Coordinate &p2,
 			const geom::Coordinate &p3);
 
+	/**
+	 * Adds the linear segments defined by an array of coordinates
+	 * to the linear centroid accumulators.
+	 *
+	 * This is done in case the polygon(s) have zero-area,
+	 * in which case the linear centroid is computed instead.
+	 *
+	 * @param pts an array of {@link Coordinate}s
+	 */
+	void addLinearSegments(const geom::CoordinateSequence& pts);
+
 };
 
 } // namespace geos::algorithm

Modified: trunk/src/algorithm/CentroidArea.cpp
===================================================================
--- trunk/src/algorithm/CentroidArea.cpp	2012-06-22 15:02:06 UTC (rev 3676)
+++ trunk/src/algorithm/CentroidArea.cpp	2012-06-22 15:03:33 UTC (rev 3677)
@@ -13,6 +13,8 @@
  *
  **********************************************************************
  *
+ * Last port: algorithm/CentroidArea.java r612
+ *
  **********************************************************************/
 
 #include <geos/algorithm/CentroidArea.h>
@@ -56,20 +58,25 @@
 	addShell(ring);
 }
 
+/* TODO: deprecate this */
 Coordinate*
 CentroidArea::getCentroid() const
 {
 	Coordinate *cent = new Coordinate();
-	cent->x = cg3.x/3.0/areasum2;
-	cent->y = cg3.y/3.0/areasum2;
+	getCentroid(*cent); // or return NULL on failure !
 	return cent;
 }
 
 bool
 CentroidArea::getCentroid(Coordinate& ret) const
 {
-	if ( areasum2 == 0.0 ) return false;
-	ret=Coordinate(cg3.x/3.0/areasum2, cg3.y/3.0/areasum2);
+	if ( areasum2 ) {
+		ret = Coordinate(cg3.x/3.0/areasum2, cg3.y/3.0/areasum2);
+	} else if ( totalLength ) {
+		ret = Coordinate(centSum.x/totalLength, centSum.y/totalLength);
+	} else {
+		return false;
+	}
 	return true;
 }
 
@@ -93,11 +100,12 @@
 CentroidArea::addShell(const CoordinateSequence *pts)
 {
 	bool isPositiveArea=!CGAlgorithms::isCCW(pts);
-    std::size_t const n=pts->getSize()-1;
+	std::size_t const n=pts->getSize()-1;
 	for(std::size_t i=0; i<n; ++i)
 	{
 		addTriangle(basePt, pts->getAt(i), pts->getAt(i+1), isPositiveArea);
 	}
+	addLinearSegments(*pts);
 }
 
 void
@@ -109,6 +117,7 @@
 	{
 		addTriangle(basePt, pts->getAt(i), pts->getAt(i+1), isPositiveArea);
 	}
+	addLinearSegments(*pts);
 }
 
 void
@@ -146,5 +155,20 @@
 	return (p2.x-p1.x)*(p3.y-p1.y)-(p3.x-p1.x)*(p2.y-p1.y);
 }
 
+void
+CentroidArea::addLinearSegments(const geom::CoordinateSequence& pts)
+{
+	std::size_t const n = pts.size()-1;
+	for (std::size_t i = 0; i < n; ++i) {
+		double segmentLen = pts[i].distance(pts[i + 1]);
+		totalLength += segmentLen;
+
+		double midx = (pts[i].x + pts[i + 1].x) / 2;
+		centSum.x += segmentLen * midx;
+		double midy = (pts[i].y + pts[i + 1].y) / 2;
+		centSum.y += segmentLen * midy;
+	}
+}
+
 } // namespace geos.algorithm
 } //namespace geos

Modified: trunk/tests/unit/Makefile.am
===================================================================
--- trunk/tests/unit/Makefile.am	2012-06-22 15:02:06 UTC (rev 3676)
+++ trunk/tests/unit/Makefile.am	2012-06-22 15:03:33 UTC (rev 3677)
@@ -104,6 +104,7 @@
 	capi/GEOSCoordSeqTest.cpp \
 	capi/GEOSGeomFromWKBTest.cpp \
 	capi/GEOSGeomToWKTTest.cpp \
+	capi/GEOSGetCentroidTest.cpp \
 	capi/GEOSContainsTest.cpp \
 	capi/GEOSConvexHullTest.cpp \
 	capi/GEOSDistanceTest.cpp \

Added: trunk/tests/unit/capi/GEOSGetCentroidTest.cpp
===================================================================
--- trunk/tests/unit/capi/GEOSGetCentroidTest.cpp	                        (rev 0)
+++ trunk/tests/unit/capi/GEOSGetCentroidTest.cpp	2012-06-22 15:03:33 UTC (rev 3677)
@@ -0,0 +1,155 @@
+// $Id$
+// 
+// Test Suite for C-API GEOSGetCentroid
+
+#include <tut.hpp>
+// geos
+#include <geos_c.h>
+// std
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+
+namespace tut
+{
+    //
+    // Test Group
+    //
+
+    // Common data used in test cases.
+    struct test_capicentroid_data
+    {
+        GEOSGeometry* geom1_;
+        GEOSGeometry* geom2_;
+        GEOSWKTWriter* wktw_;
+        char* wkt_;
+        double area_;
+
+        static void notice(const char *fmt, ...)
+        {
+            std::fprintf( stdout, "NOTICE: ");
+
+            va_list ap;
+            va_start(ap, fmt);
+            std::vfprintf(stdout, fmt, ap);
+            va_end(ap);
+        
+            std::fprintf(stdout, "\n");
+        }
+
+        test_capicentroid_data()
+            : geom1_(0), geom2_(0), wkt_(0)
+        {
+            initGEOS(notice, notice);
+            wktw_ = GEOSWKTWriter_create();
+            GEOSWKTWriter_setTrim(wktw_, 1);
+            GEOSWKTWriter_setRoundingPrecision(wktw_, 8);
+        }       
+
+        ~test_capicentroid_data()
+        {
+            GEOSGeom_destroy(geom1_);
+            GEOSGeom_destroy(geom2_);
+            GEOSWKTWriter_destroy(wktw_);
+            GEOSFree(wkt_);
+            geom1_ = 0;
+            geom2_ = 0;
+            wkt_ = 0;
+            finishGEOS();
+        }
+
+    };
+
+    typedef test_group<test_capicentroid_data> group;
+    typedef group::object object;
+
+    group test_capicentroid_group("capi::GEOSGetCentroid");
+
+    //
+    // Test Cases
+    //
+
+    // Single point
+    template<>
+    template<>
+    void object::test<1>()
+    {
+        geom1_ = GEOSGeomFromWKT("POINT(10 0)");
+
+        ensure( 0 != geom1_ );
+
+        geom2_ = GEOSGetCentroid(geom1_);
+
+        ensure( 0 != geom2_ );
+
+        wkt_ = GEOSWKTWriter_write(wktw_, geom2_);
+
+        ensure_equals(std::string(wkt_), std::string( "POINT (10 0)"));
+
+    }
+
+    // line
+    template<>
+    template<>
+    void object::test<2>()
+    {
+        geom1_ = GEOSGeomFromWKT("LINESTRING(0 0, 10 0)");
+
+        ensure( 0 != geom1_ );
+
+        geom2_ = GEOSGetCentroid(geom1_);
+
+        ensure( 0 != geom2_ );
+
+        wkt_ = GEOSWKTWriter_write(wktw_, geom2_);
+
+        ensure_equals(std::string(wkt_), std::string( "POINT (5 0)"));
+
+    }
+
+    // polygon
+    template<>
+    template<>
+    void object::test<3>()
+    {
+        geom1_ = GEOSGeomFromWKT("POLYGON((0 0, 10 0, 10 10, 0 10, 0 0))");
+
+        ensure( 0 != geom1_ );
+
+        geom2_ = GEOSGetCentroid(geom1_);
+
+        ensure( 0 != geom2_ );
+
+        wkt_ = GEOSWKTWriter_write(wktw_, geom2_);
+
+        ensure_equals(std::string(wkt_), std::string( "POINT (5 5)"));
+
+    }
+
+    // Tiny triangle, see http://trac.osgeo.org/geos/ticket/559
+    template<>
+    template<>
+    void object::test<4>()
+    {
+        geom1_ = GEOSGeomFromWKT(
+"POLYGON(( \
+56.528666666700 25.2101666667, \
+56.529000000000 25.2105000000, \
+56.528833333300 25.2103333333, \
+56.528666666700 25.2101666667))");
+
+        ensure( 0 != geom1_ );
+
+        geom2_ = GEOSGetCentroid(geom1_);
+
+        ensure( 0 != geom2_ );
+
+        wkt_ = GEOSWKTWriter_write(wktw_, geom2_);
+
+        ensure_equals(std::string(wkt_), std::string( "POINT (56.528833 25.210333)" ) );
+
+    }
+
+} // namespace tut
+



More information about the geos-commits mailing list