[geos-commits] [SCM] GEOS branch master updated. 67d6eafd9a95fdb9890bbcdce6ce85aab0c1dce2

git at osgeo.org git at osgeo.org
Fri Oct 26 09:14:34 PDT 2018


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  67d6eafd9a95fdb9890bbcdce6ce85aab0c1dce2 (commit)
      from  8642d9def6d1df4a14769a74f85b8501f18c527d (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 67d6eafd9a95fdb9890bbcdce6ce85aab0c1dce2
Author: Darafei Praliaskouski <me at komzpa.net>
Date:   Sat Sep 1 16:50:00 2018 +0300

    Fix Delaunay robustness by using longer floating point variable on inCircle

diff --git a/src/triangulate/quadedge/TrianglePredicate.cpp b/src/triangulate/quadedge/TrianglePredicate.cpp
index 51dfcc3..2030ad2 100644
--- a/src/triangulate/quadedge/TrianglePredicate.cpp
+++ b/src/triangulate/quadedge/TrianglePredicate.cpp
@@ -42,22 +42,30 @@ TrianglePredicate::isInCircleNormalized(
 		const Coordinate &a, const Coordinate &b, const Coordinate &c,
 		const Coordinate &p)
 {
-	double adx = a.x - p.x;
-	double ady = a.y - p.y;
-	double bdx = b.x - p.x;
-	double bdy = b.y - p.y;
-	double cdx = c.x - p.x;
-	double cdy = c.y - p.y;
+	// Unfortunately this implementation is not robust either. For robust one see:
+	// https://www.cs.cmu.edu/~quake/robust.html
+	// https://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c
 
-	double abdet = adx * bdy - bdx * ady;
-	double bcdet = bdx * cdy - cdx * bdy;
-	double cadet = cdx * ady - adx * cdy;
-	double alift = adx * adx + ady * ady;
-	double blift = bdx * bdx + bdy * bdy;
-	double clift = cdx * cdx + cdy * cdy;
+	long double adx = a.x - p.x;
+	long double ady = a.y - p.y;
+	long double bdx = b.x - p.x;
+	long double bdy = b.y - p.y;
+	long double cdx = c.x - p.x;
+	long double cdy = c.y - p.y;
 
-	double disc = alift * bcdet + blift * cadet + clift * abdet;
-	return disc > 0;
+	long double bdxcdy = bdx * cdy;
+	long double cdxbdy = cdx * bdy;
+	long double alift = adx * adx + ady * ady;
+
+	long double cdxady = cdx * ady;
+	long double adxcdy = adx * cdy;
+	long double blift = bdx * bdx + bdy * bdy;
+
+	long double adxbdy = adx * bdy;
+	long double bdxady = bdx * ady;
+	long double clift = cdx * cdx + cdy * cdy;
+	return (alift * bdxcdy + blift * cdxady + clift * adxbdy) >
+		(alift * cdxbdy + blift * adxcdy + clift * bdxady);
 }
 
 double
@@ -73,11 +81,9 @@ TrianglePredicate::isInCircleRobust(
 		const Coordinate &a, const Coordinate &b, const Coordinate &c,
 		const Coordinate &p)
 {
-	//checkRobustInCircle(a, b, c, p);
-	//	return isInCircleNonRobust(a, b, c, p);
+	// This implementation is not robust, name is ported from JTS.
 	return isInCircleNormalized(a, b, c, p);
 }
 
 } // namespace geos.geom
 } // namespace geos
-
diff --git a/tests/unit/triangulate/DelaunayTest.cpp b/tests/unit/triangulate/DelaunayTest.cpp
index 21db61e..0467b20 100644
--- a/tests/unit/triangulate/DelaunayTest.cpp
+++ b/tests/unit/triangulate/DelaunayTest.cpp
@@ -128,12 +128,13 @@ namespace tut
 	}
 
 	// 5 - Test Circle
+	// Added a point inside to ensure single possible solution
 	template<>
 	template<>
 	void object::test<5>()
 	{
-		const char * wkt = "POLYGON ((42 30, 41.96 29.61, 41.85 29.23, 41.66 28.89, 41.41 28.59, 41.11 28.34, 40.77 28.15, 40.39 28.04, 40 28, 39.61 28.04, 39.23 28.15, 38.89 28.34, 38.59 28.59, 38.34 28.89, 38.15 29.23, 38.04 29.61, 38 30, 38.04 30.39, 38.15 30.77, 38.34 31.11, 38.59 31.41, 38.89 31.66, 39.23 31.85, 39.61 31.96, 40 32, 40.39 31.96, 40.77 31.85, 41.11 31.66, 41.41 31.41, 41.66 31.11, 41.85 30.77, 41.96 30.39, 42 30))";
-		const char* expectedEdges = "MULTILINESTRING ((41.66 31.11, 41.85 30.77), (41.41 31.41, 41.66 31.11), (41.11 31.66, 41.41 31.41), (40.77 31.85, 41.11 31.66), (40.39 31.96, 40.77 31.85), (40 32, 40.39 31.96), (39.61 31.96, 40 32), (39.23 31.85, 39.61 31.96), (38.89 31.66, 39.23 31.85), (38.59 31.41, 38.89 31.66), (38.34 31.11, 38.59 31.41), (38.15 30.77, 38.34 31.11), (38.04 30.39, 38.15 30.77), (38 30, 38.04 30.39), (38 30, 38.04 29.61), (38.04 29.61, 38.15 29.23), (38.15 29.23, 38.34 28.89), (38.34 28.89, 38.59 28.59), (38.59 28.59, 38.89 28.34), (38.89 28.34, 39.23 28.15), (39.23 28.15, 39.61 28.04), (39.61 28.04, 40 28), (40 28, 40.39 28.04), (40.39 28.04, 40.77 28.15), (40.77 28.15, 41.11 28.34), (41.11 28.34, 41.41 28.59), (41.41 28.59, 41.66 28.89), (41.66 28.89, 41.85 29.23), (41.85 29.23, 41.96 29.61), (41.96 29.61, 42 30), (41.96 30.39, 42 30), (41.85 30.77, 41.96 30.39), (41.66 31.11, 41.96 30.39), (41.41 31.41, 41.96 30.39), (41.41 28.59, 41.96 30.39), (41.41 28.59, 41.
 41 31.41), (38.59 28.59, 41.41 28.59), (38.59 28.59, 41.41 31.41), (38.59 28.59, 38.59 31.41), (38.59 31.41, 41.41 31.41), (38.59 31.41, 39.61 31.96), (39.61 31.96, 41.41 31.41), (39.61 31.96, 40.39 31.96), (40.39 31.96, 41.41 31.41), (40.39 31.96, 41.11 31.66), (38.04 30.39, 38.59 28.59), (38.04 30.39, 38.59 31.41), (38.04 30.39, 38.34 31.11), (38.04 29.61, 38.59 28.59), (38.04 29.61, 38.04 30.39), (39.61 28.04, 41.41 28.59), (38.59 28.59, 39.61 28.04), (38.89 28.34, 39.61 28.04), (40.39 28.04, 41.41 28.59), (39.61 28.04, 40.39 28.04), (41.96 29.61, 41.96 30.39), (41.41 28.59, 41.96 29.61), (41.66 28.89, 41.96 29.61), (40.39 28.04, 41.11 28.34), (38.04 29.61, 38.34 28.89), (38.89 31.66, 39.61 31.96))";
+		const char * wkt = "GEOMETRYCOLLECTION(POLYGON ((42 30, 41.96 29.61, 41.85 29.23, 41.66 28.89, 41.41 28.59, 41.11 28.34, 40.77 28.15, 40.39 28.04, 40 28, 39.61 28.04, 39.23 28.15, 38.89 28.34, 38.59 28.59, 38.34 28.89, 38.15 29.23, 38.04 29.61, 38 30, 38.04 30.39, 38.15 30.77, 38.34 31.11, 38.59 31.41, 38.89 31.66, 39.23 31.85, 39.61 31.96, 40 32, 40.39 31.96, 40.77 31.85, 41.11 31.66, 41.41 31.41, 41.66 31.11, 41.85 30.77, 41.96 30.39, 42 30)), POINT(38.6 30))";
+		const char* expectedEdges = "MULTILINESTRING((41.96 30.39,42 30),(41.96 29.61,42 30),(41.85 30.77,41.96 30.39),(41.85 29.23,41.96 29.61),(41.66 31.11,41.85 30.77),(41.66 28.89,41.85 29.23),(41.41 31.41,41.66 31.11),(41.41 28.59,41.66 28.89),(41.11 31.66,41.41 31.41),(41.11 28.34,41.41 28.59),(40.77 31.85,41.11 31.66),(40.77 28.15,41.11 28.34),(40.39 31.96,40.77 31.85),(40.39 28.04,40.77 28.15),(40 32,40.39 31.96),(40 28,40.39 28.04),(39.61 31.96,40 32),(39.61 28.04,40 28),(39.23 31.85,39.61 31.96),(39.23 28.15,39.61 28.04),(38.89 31.66,39.23 31.85),(38.89 28.34,39.23 28.15),(38.6 30,42 30),(38.6 30,41.96 30.39),(38.6 30,41.96 29.61),(38.6 30,41.85 30.77),(38.6 30,41.85 29.23),(38.6 30,41.66 31.11),(38.6 30,41.66 28.89),(38.6 30,41.41 31.41),(38.6 30,41.41 28.59),(38.6 30,41.11 31.66),(38.6 30,41.11 28.34),(38.6 30,40.77 31.85),(38.6 30,40.77 28.15),(38.6 30,40.39 31.96),(38.6 30,40.39 28.04),(38.6 30,40 32),(38.6 30,40 28),(38.6 30,39.61 31.96),(38.6 30,39.61 28.04),(38.6 30,39.23
  31.85),(38.6 30,39.23 28.15),(38.6 30,38.89 31.66),(38.6 30,38.89 28.34),(38.59 31.41,38.89 31.66),(38.59 31.41,38.6 30),(38.59 28.59,38.89 28.34),(38.59 28.59,38.6 30),(38.34 31.11,38.6 30),(38.34 31.11,38.59 31.41),(38.34 28.89,38.6 30),(38.34 28.89,38.59 28.59),(38.15 30.77,38.6 30),(38.15 30.77,38.34 31.11),(38.15 29.23,38.6 30),(38.15 29.23,38.34 28.89),(38.04 30.39,38.6 30),(38.04 30.39,38.15 30.77),(38.04 29.61,38.6 30),(38.04 29.61,38.15 29.23),(38 30,38.6 30),(38 30,38.04 30.39),(38 30,38.04 29.61))";
 
 		runDelaunay(wkt, false, expectedEdges);
 	}
@@ -172,6 +173,7 @@ namespace tut
 
 		runDelaunay(wkt, false, expectedEdges, 0.001);
 	}
+
 	// 9 - Test for DelaunayTriangulationBuilder::envelope
 	template<>
 	template<>
@@ -189,5 +191,21 @@ namespace tut
 		ensure_equals(env.getHeight() , 107);
 	}
 
+	// 10 - Tolerance robustness
+	template<>
+	template<>
+	void object::test<10>()
+	{
+		const char * wkt = "MULTIPOINT(63.547558624186912368 70.904719023616522122,63.547558624186969212 70.904719023616564755,66.103648384371410884 68.588612471664760051,77.882918707497154287 74.870889977331813725,128.47759065022572145 177.65366864730182783)";
+		const char* expectedEdges = "GEOMETRYCOLLECTION (POLYGON ((63.5475586241869692 70.9047190236165648, 128.4775906502257214 177.6536686473018278, 77.8829187074971543 74.8708899773318137, 63.5475586241869692 70.9047190236165648)), POLYGON ((63.5475586241869692 70.9047190236165648, 77.8829187074971543 74.8708899773318137, 66.1036483843714109 68.5886124716647601, 63.5475586241869692 70.9047190236165648)), POLYGON ((63.5475586241869124 70.9047190236165221, 128.4775906502257214 177.6536686473018278, 63.5475586241869692 70.9047190236165648, 63.5475586241869124 70.9047190236165221)), POLYGON ((63.5475586241869124 70.9047190236165221, 63.5475586241869692 70.9047190236165648, 66.1036483843714109 68.5886124716647601, 63.5475586241869124 70.9047190236165221)))";
+
+		// inCircle predicate can't handle it on machines with 64-bit long double
+		if (sizeof(long double) > sizeof(double))
+		{
+			runDelaunay(wkt, true, expectedEdges, 0.0);
+		}
+	}
+
+
 } // namespace tut
 

-----------------------------------------------------------------------

Summary of changes:
 src/triangulate/quadedge/TrianglePredicate.cpp | 40 +++++++++++++++-----------
 tests/unit/triangulate/DelaunayTest.cpp        | 22 ++++++++++++--
 2 files changed, 43 insertions(+), 19 deletions(-)


hooks/post-receive
-- 
GEOS


More information about the geos-commits mailing list