Use tighter bounds and faster filters in orientation and isInCircle p…

…redicates.
libgeos · Sep 14, 2024 · 4143db7 · 4143db7
1 parent 701c868
commit 4143db7
Show file tree

Hide file tree

Showing 2 changed files with 59 additions and 45 deletions.
diff --git a/include/geos/algorithm/CGAlgorithmsDD.h b/include/geos/algorithm/CGAlgorithmsDD.h
@@ -20,6 +20,7 @@
 
 #include <geos/export.h>
 #include <geos/math/DD.h>
+#include <cmath>
 
 // Forward declarations
 namespace geos {
@@ -92,41 +93,14 @@ class GEOS_DLL CGAlgorithmsDD {
         double pbx, double pby,
         double pcx, double pcy)
     {
-        /**
-         * A value which is safely greater than the relative round-off
-         * error in double-precision numbers
-         */
-        double constexpr DP_SAFE_EPSILON =  1e-15;
-
-        double detsum;
         double const detleft = (pax - pcx) * (pby - pcy);
         double const detright = (pay - pcy) * (pbx - pcx);
         double const det = detleft - detright;
-
-        if(detleft > 0.0) {
-            if(detright <= 0.0) {
-                return orientation(det);
-            }
-            else {
-                detsum = detleft + detright;
-            }
-        }
-        else if(detleft < 0.0) {
-            if(detright >= 0.0) {
-                return orientation(det);
-            }
-            else {
-                detsum = -detleft - detright;
-            }
-        }
-        else {
-            return orientation(det);
-        }
-
-        double const errbound = DP_SAFE_EPSILON * detsum;
-        if((det >= errbound) || (-det >= errbound)) {
-            return orientation(det);
-        }
+        // Coefficient due to https://doi.org/10.1007/s10543-015-0574-9
+        double const error = std::abs(detleft + detright)
+                             * 3.3306690621773714e-16;
+        if (std::abs(det) >= error)
+            return (det > 0) - (det < 0);
         return CGAlgorithmsDD::FAILURE;
     };
 

diff --git a/src/triangulate/quadedge/TrianglePredicate.cpp b/src/triangulate/quadedge/TrianglePredicate.cpp
@@ -29,16 +29,30 @@ namespace quadedge {
 
 geom::Location
 TrianglePredicate::isInCircleNonRobust(
-    const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c,
-    const CoordinateXY& p)
+    const CoordinateXY& p, const CoordinateXY& q, const CoordinateXY& r,
+    const CoordinateXY& t)
 {
-    auto det =
-        (a.x * a.x + a.y * a.y) * triArea(b, c, p)
-        - (b.x * b.x + b.y * b.y) * triArea(a, c, p)
-        + (c.x * c.x + c.y * c.y) * triArea(a, b, p)
-        - (p.x * p.x + p.y * p.y) * triArea(a, b, c);
-
-    return det > 0 ? geom::Location::EXTERIOR : (det < 0 ? geom::Location::INTERIOR : geom::Location::BOUNDARY);
+    auto qpx = q.x - p.x;
+    auto qpy = q.y - p.y;
+    auto rpx = r.x - p.x;
+    auto rpy = r.y - p.y;
+    auto tpx = t.x - p.x;
+    auto tpy = t.y - p.y;
+    auto tqx = t.x - q.x;
+    auto tqy = t.y - q.y;
+    auto rqx = r.x - q.x;
+    auto rqy = r.y - q.y;
+    auto qpxtpy = qpx * tpy;
+    auto qpytpx = qpy * tpx;
+    auto tpxtqx = tpx * tqx;
+    auto tpytqy = tpy * tqy;
+    auto qpxrpy = qpx * rpy;
+    auto qpyrpx = qpy * rpx;
+    auto rpxrqx = rpx * rqx;
+    auto rpyrqy = rpy * rqy;
+    auto det = (qpxtpy - qpytpx) * (rpxrqx + rpyrqy)
+               - (qpxrpy - qpyrpx) * (tpxtqx + tpytqy);
+    return static_cast<geom::Location>( (det > 0) - (det < 0) + 1 );
 }
 
 geom::Location
@@ -91,11 +105,37 @@ TrianglePredicate::triArea(const CoordinateXY& a,
 
 geom::Location
 TrianglePredicate::isInCircleRobust(
-    const CoordinateXY& a, const CoordinateXY& b, const CoordinateXY& c,
-    const CoordinateXY& p)
+    const CoordinateXY& q, const CoordinateXY& p, const CoordinateXY& r,
+    const CoordinateXY& t)
 {
-    // This implementation is not robust, name is ported from JTS.
-    return isInCircleNormalized(a, b, c, p);
+    // This implementation is not robust and defaults to BOUNDARY in case of
+    // uncertainty.
+    auto qpx = q.x - p.x;
+    auto qpy = q.y - p.y;
+    auto rpx = r.x - p.x;
+    auto rpy = r.y - p.y;
+    auto tpx = t.x - p.x;
+    auto tpy = t.y - p.y;
+    auto tqx = t.x - q.x;
+    auto tqy = t.y - q.y;
+    auto rqx = r.x - q.x;
+    auto rqy = r.y - q.y;
+    auto qpxtpy = qpx * tpy;
+    auto qpytpx = qpy * tpx;
+    auto tpxtqx = tpx * tqx;
+    auto tpytqy = tpy * tqy;
+    auto qpxrpy = qpx * rpy;
+    auto qpyrpx = qpy * rpx;
+    auto rpxrqx = rpx * rqx;
+    auto rpyrqy = rpy * rqy;
+    auto det = (qpxtpy - qpytpx) * (rpxrqx + rpyrqy)
+               - (qpxrpy - qpyrpx) * (tpxtqx + tpytqy);
+    auto deterror = ((std::abs(qpxtpy) + std::abs(qpytpx))
+                    * (std::abs(rpxrqx) + std::abs(rpyrqy))
+                    + (std::abs(qpxrpy) + std::abs(qpyrpx))
+                    * (std::abs(tpxtqx) + std::abs(tpytqy)))
+                    * 9.99200719823023e-16;
+    return static_cast<geom::Location>( (det > deterror) - (det < -deterror) + 1 );
 }
 
 } // namespace geos.triangulate.quadedge