diff --git a/Triangle.NET/Triangle/Geometry/Contour.cs b/Triangle.NET/Triangle/Geometry/Contour.cs index f3d0ac8..6cd6ae6 100644 --- a/Triangle.NET/Triangle/Geometry/Contour.cs +++ b/Triangle.NET/Triangle/Geometry/Contour.cs @@ -73,28 +73,43 @@ namespace TriangleNet.Geometry return segments; } - public Point FindInteriorPoint() + /// + /// Try to find a point inside the contour. + /// + /// The number of iterations on each segment (default = 8). + /// Threshold for co-linear points (default = 2e-6). + /// Point inside the contour + /// Throws if no point could be found. + /// + /// 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 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). + /// + 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 points) @@ -110,50 +125,81 @@ namespace TriangleNet.Geometry } } - private static Point FindPointInPolygon(List contour) + #region Helper methods + + private static Point FindPointInPolygon(List 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 } }