Improve robustness of Contour.FindInteriorPoint

git-svn-id: https://triangle.svn.codeplex.com/svn@77800 0e2699bc-83d4-4a8f-98e7-55e24ab8c7a5

Triangle.NET/Triangle/Geometry/Contour.cs

@@ -73,28 +73,43 @@ namespace TriangleNet.Geometry
             return segments;
         }
 
-        public Point FindInteriorPoint()
+        /// <summary>
+        /// Try to find a point inside the contour.
+        /// </summary>
+        /// <param name="limit">The number of iterations on each segment (default = 8).</param>
+        /// <param name="eps">Threshold for co-linear points (default = 2e-6).</param>
+        /// <returns>Point inside the contour</returns>
+        /// <exception cref="Exception">Throws if no point could be found.</exception>
+        /// <remarks>
+        /// For each corner (index i) of the contour, the 3 points with indices i-1, i and i+1
+        /// are considered and a search on the line through the corner vertex is started (either
+        /// on the bisecting line, or, if <see cref="IPredicates.CounterClockwise"/> is less than
+        /// eps, on the perpendicular line.
+        /// A given number of points will be tested (limit), while the distance to the contour
+        /// boundary will be reduced in each iteration (with a factor of 1/2^i, i = 1 ... limit).
+        /// </remarks>
+        public Point FindInteriorPoint(int limit = 8, double eps = 2e-6)
         {
+            var point = new Point(0.0, 0.0);
+
             if (convex)
             {
                 int count = this.Points.Count;
 
-                var centroid = new Point(0.0, 0.0);
-
                 for (int i = 0; i < count; i++)
                 {
-                    centroid.x += this.Points[i].x;
-                    centroid.y += this.Points[i].y;
+                    point.x += this.Points[i].x;
+                    point.y += this.Points[i].y;
                 }
 
-                // If the hole is convex, use its centroid.
-                centroid.x /= count;
-                centroid.y /= count;
+                // If the contour is convex, use its centroid.
+                point.x /= count;
+                point.y /= count;
 
-                return centroid;
+                return point;
             }
 
-            return FindPointInPolygon(this.Points);
+            return FindPointInPolygon(this.Points, limit, eps);
         }
 
         private void AddPoints(IEnumerable<Vertex> points)
@@ -110,50 +125,81 @@ namespace TriangleNet.Geometry
             }
         }
 
-        private static Point FindPointInPolygon(List<Vertex> contour)
+        #region Helper methods
+
+        private static Point FindPointInPolygon(List<Vertex> contour, int limit, double eps)
         {
             var bounds = new Rectangle();
             bounds.Expand(contour);
 
             int length = contour.Count;
-            int limit = 8;
 
             var test = new Point();
 
-            Point a, b; // Current edge.
-            double cx, cy; // Center of current edge.
-            double dx, dy; // Direction perpendicular to edge.
+            Point a, b, c; // Current corner points.
+
+            double bx, by;
+            double dx, dy;
+            double h;
+
+            var predicates = new RobustPredicates();
+
+            a = contour[0];
+            b = contour[1];
 
             for (int i = 0; i < length; i++)
             {
-                a = contour[i];
-                b = contour[(i + 1) % length];
+                c = contour[(i + 2) % length];
 
-                cx = (a.x + b.x) / 2;
-                cy = (a.y + b.y) / 2;
+                // Corner point.
+                bx = b.X;
+                by = b.Y;
 
-                dx = (b.y - a.y) / 1.374;
-                dy = (a.x - b.x) / 1.374;
+                // NOTE: if we knew the contour points were in counterclockwise order, we
+                // could skip concave corners and search only in one direction.
 
-                for (int j = 1; j <= limit; j++)
+                h = predicates.CounterClockwise(a, b, c);
+
+                if (Math.Abs(h) < eps)
                 {
-                    // Search to the right of the segment.
-                    test.x = cx + dx / j;
-                    test.y = cy + dy / j;
+                    // Points are nearly co-linear. Use perpendicular direction.
+                    dx = (c.Y - a.Y) / 2;
+                    dy = (a.X - c.X) / 2;
+                }
+                else
+                {
+                    // Direction [midpoint(a-c) -> corner point]
+                    dx = (a.X + c.X) / 2 - bx;
+                    dy = (a.Y + c.Y) / 2 - by;
+                }
+
+                // Move around the contour.
+                a = b;
+                b = c;
+
+                h = 1.0;
+
+                for (int j = 0; j < limit; j++)
+                {
+                    // Search in direction.
+                    test.X = bx + dx / h;
+                    test.Y = by + dy / h;
 
                     if (bounds.Contains(test) && IsPointInPolygon(test, contour))
                     {
                         return test;
                     }
 
-                    // Search on the other side of the segment.
-                    test.x = cx - dx / j;
-                    test.y = cy - dy / j;
+                    // Search in opposite direction (see NOTE above).
+                    test.X = bx - dx / h;
+                    test.Y = by - dy / h;
 
                     if (bounds.Contains(test) && IsPointInPolygon(test, contour))
                     {
                         return test;
                     }
+
+                    h = 2.0 * h;
                 }
             }
 
@@ -194,5 +240,7 @@ namespace TriangleNet.Geometry
 
             return inside;
         }
+
+        #endregion
     }
 }