From 2e83d7ee5e66d50a7d8ac7c169c676c9050e4255 Mon Sep 17 00:00:00 2001
From: wo80 <christian.woltering@tu-dortmund.de>
Date: Wed, 23 Feb 2022 12:42:54 +0100
Subject: [PATCH] Improve contour FindPointInPolygon() reliability (2).

---
 .../Tools/IntersectionHelperTest.cs           | 41 +++++++++++++++++++
 src/Triangle/Geometry/Contour.cs              |  7 ++--
 src/Triangle/Tools/IntersectionHelper.cs      | 32 +++++++++++++++
 3 files changed, 76 insertions(+), 4 deletions(-)
 create mode 100644 src/Triangle.Tests/Tools/IntersectionHelperTest.cs

diff --git a/src/Triangle.Tests/Tools/IntersectionHelperTest.cs b/src/Triangle.Tests/Tools/IntersectionHelperTest.cs
new file mode 100644
index 0000000..d8f2c69
--- /dev/null
+++ b/src/Triangle.Tests/Tools/IntersectionHelperTest.cs
@@ -0,0 +1,41 @@
+using NUnit.Framework;
+using TriangleNet.Geometry;
+using TriangleNet.Tools;
+
+namespace TriangleNet.Tests.Tools
+{
+    public class IntersectionHelperTest
+    {
+        [Test]
+        public void TestIsPointOnSegment()
+        {
+            var a = new Vertex(1.0, 1.0);
+            var b = new Vertex(2.0, 2.0);
+
+            // Test point = segment start point.
+            Assert.IsTrue(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(1.0, 1.0)));
+
+            // Test point = segment end point.
+            Assert.IsTrue(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(2.0, 2.0)));
+
+            // Test point on segment.
+            Assert.IsTrue(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(1.5, 1.5)));
+
+            // Test point collinear, but not on segment.
+            Assert.IsFalse(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(0.0, 0.0)));
+            Assert.IsFalse(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(3.0, 3.0)));
+
+            // Test point not on segment.
+            Assert.IsFalse(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(1.5, 0.5)));
+
+            double eps = 1e-12;
+
+            // Test point collinear near endpoint, but not on segment.
+            Assert.IsFalse(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(2.0 + eps, 2.0 + eps)));
+            Assert.IsFalse(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(2.0 - eps, 2.0 + eps)));
+
+            // Test point collinear near endpoint on segment.
+            Assert.IsTrue(IntersectionHelper.IsPointOnSegment(a, b, new Vertex(2.0 - eps, 2.0 - eps)));
+        }
+    }
+}
diff --git a/src/Triangle/Geometry/Contour.cs b/src/Triangle/Geometry/Contour.cs
index 68cb55d..0617ed6 100644
--- a/src/Triangle/Geometry/Contour.cs
+++ b/src/Triangle/Geometry/Contour.cs
@@ -8,6 +8,7 @@ namespace TriangleNet.Geometry
 {
     using System;
     using System.Collections.Generic;
+    using TriangleNet.Tools;
 
     public class Contour
     {
@@ -245,17 +246,15 @@ namespace TriangleNet.Geometry
         /// <summary>
         /// Work around IsPointInPolygon() failing for points on segments.
         /// </summary>
-        private static bool IsPointOnSegment(Point test, List<Vertex> contour, double esp = 1e-12)
+        private static bool IsPointOnSegment(Point test, List<Vertex> contour, double eps = 1e-12)
         {
-            var p = RobustPredicates.Default;
-
             int count = contour.Count;
 
             int i = count - 1;
 
             for (int j = 0; j < count; j++)
             {
-                if (Math.Abs(p.CounterClockwise(contour[i], test, contour[j])) < esp)
+                if (IntersectionHelper.IsPointOnSegment(contour[i], contour[j], test, eps))
                 {
                     return true;
                 }
diff --git a/src/Triangle/Tools/IntersectionHelper.cs b/src/Triangle/Tools/IntersectionHelper.cs
index 478eafd..bbe663a 100644
--- a/src/Triangle/Tools/IntersectionHelper.cs
+++ b/src/Triangle/Tools/IntersectionHelper.cs
@@ -6,10 +6,42 @@
 
 namespace TriangleNet.Tools
 {
+    using System;
     using TriangleNet.Geometry;
 
     public static class IntersectionHelper
     {
+        /// <summary>
+        /// Check if a given test point lies on a segment.
+        /// </summary>
+        /// <param name="a">The segment start point.</param>
+        /// <param name="b">The segment end point.</param>
+        /// <param name="test">The point to test.</param>
+        /// <param name="eps">Threshold to test collinearity (default = 1e-12).</param>
+        /// <returns></returns>
+        public static bool IsPointOnSegment(Point a, Point b, Point test, double eps = 1e-12)
+        {
+            // The cross product.
+            double cross = (test.Y - a.Y) * (b.X - a.X) - (test.X - a.X) * (b.Y - a.Y);
+
+            // Check if points are collinear.
+            if (Math.Abs(cross) > eps) return false;
+
+            // The dot product (projection of test point onto segment).
+            double dot = (test.X - a.X) * (b.X - a.X) + (test.Y - a.Y) * (b.Y - a.Y);
+
+            // Check if test point is actually between a and b (left of a).
+            if (dot < 0) return false;
+
+            // Length of the segment.
+            double ab = (b.X - a.X) * (b.X - a.X) + (b.Y - a.Y) * (b.Y - a.Y);
+
+            // Check if test point is actually between a and b (right of b).
+            if (dot > ab) return false;
+
+            return true;
+        }
+
         /// <summary>
         /// Compute intersection of two segments.
         /// </summary>