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

---
 src/Triangle.Tests/Geomerty/ContourTest.cs | 40 ++++++++++++++++++++++
 src/Triangle/Geometry/Contour.cs           | 29 ++++++++++++++--
 2 files changed, 67 insertions(+), 2 deletions(-)

diff --git a/src/Triangle.Tests/Geomerty/ContourTest.cs b/src/Triangle.Tests/Geomerty/ContourTest.cs
index 97b429e..c7a41cf 100644
--- a/src/Triangle.Tests/Geomerty/ContourTest.cs
+++ b/src/Triangle.Tests/Geomerty/ContourTest.cs
@@ -1,5 +1,7 @@
 using NUnit.Framework;
 using TriangleNet.Geometry;
+using System;
+using System.Collections.Generic;
 
 namespace TriangleNet.Tests.Geometry
 {
@@ -23,5 +25,43 @@ namespace TriangleNet.Tests.Geometry
 
             Assert.IsTrue(p.X > 0d && p.X < 1d && p.Y > 0d && p.Y < 1d);
         }
+
+        [Test]
+        public void TestFindInteriorPointL()
+        {
+            // L-shaped contour (FindPointInPolygon() produces a test candidate
+            // which lies exactly on a segment where IsPointInPolygon() returns
+            // true, so IsPointOnSegment() is actually needed here).
+            var points = new List<Vertex>()
+            {
+                new Vertex(3, 1),
+                new Vertex(1, 1),
+                new Vertex(1, 3),
+                new Vertex(2, 3),
+                new Vertex(2, 2),
+                new Vertex(3, 2)
+            };
+
+            var contour = new Contour(points);
+
+            var poly = new Polygon(6);
+
+            poly.Add(contour, true);
+
+            var h = poly.Holes[0];
+            var p = RobustPredicates.Default;
+
+            int count = points.Count;
+            int i = count - 1;
+
+            for (int j = 0; j < count; j++)
+            {
+                double ccw = p.CounterClockwise(points[i], h, points[j]);
+
+                Assert.Greater(Math.Abs(ccw), 1e-12);
+
+                i = j;
+            }
+        }
     }
 }
diff --git a/src/Triangle/Geometry/Contour.cs b/src/Triangle/Geometry/Contour.cs
index 68ab9af..68cb55d 100644
--- a/src/Triangle/Geometry/Contour.cs
+++ b/src/Triangle/Geometry/Contour.cs
@@ -185,7 +185,7 @@ namespace TriangleNet.Geometry
                     test.x = bx + dx * h;
                     test.y = by + dy * h;
 
-                    if (bounds.Contains(test) && IsPointInPolygon(test, contour))
+                    if (bounds.Contains(test) && IsPointInPolygon(test, contour) && !IsPointOnSegment(test, contour))
                     {
                         return test;
                     }
@@ -194,7 +194,7 @@ namespace TriangleNet.Geometry
                     test.x = bx - dx * h;
                     test.y = by - dy * h;
 
-                    if (bounds.Contains(test) && IsPointInPolygon(test, contour))
+                    if (bounds.Contains(test) && IsPointInPolygon(test, contour) && !IsPointOnSegment(test, contour))
                     {
                         return test;
                     }
@@ -241,6 +241,31 @@ namespace TriangleNet.Geometry
             return inside;
         }
 
+
+        /// <summary>
+        /// Work around IsPointInPolygon() failing for points on segments.
+        /// </summary>
+        private static bool IsPointOnSegment(Point test, List<Vertex> contour, double esp = 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)
+                {
+                    return true;
+                }
+
+                i = j;
+            }
+
+            return false;
+        }
+
         #endregion
     }
 }