/* * Smallest enclosing circle - Library (C#) * * Copyright (c) 2020 Project Nayuki * https://www.nayuki.io/page/smallest-enclosing-circle * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published by * the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with this program (see COPYING.txt and COPYING.LESSER.txt). * If not, see . */ using System; using System.Collections.Generic; public sealed class SmallestEnclosingCircle { /* * Returns the smallest circle that encloses all the given points. Runs in expected O(n) time, randomized. * Note: If 0 points are given, a circle of radius -1 is returned. If 1 point is given, a circle of radius 0 is returned. */ // Initially: No boundary points known public static Circle MakeCircle(IList points) { // Clone list to preserve the caller's data, do Durstenfeld shuffle List shuffled = new List(points); Random rand = new Random(); for (int i = shuffled.Count - 1; i > 0; i--) { int j = rand.Next(i + 1); Point temp = shuffled[i]; shuffled[i] = shuffled[j]; shuffled[j] = temp; } // Progressively add points to circle or recompute circle Circle c = Circle.INVALID; for (int i = 0; i < shuffled.Count; i++) { Point p = shuffled[i]; if (c.r < 0 || !c.Contains(p)) c = MakeCircleOnePoint(shuffled.GetRange(0, i + 1), p); } return c; } // One boundary point known private static Circle MakeCircleOnePoint(List points, Point p) { Circle c = new Circle(p, 0); for (int i = 0; i < points.Count; i++) { Point q = points[i]; if (!c.Contains(q)) { if (c.r == 0) c = MakeDiameter(p, q); else c = MakeCircleTwoPoints(points.GetRange(0, i + 1), p, q); } } return c; } // Two boundary points known private static Circle MakeCircleTwoPoints(List points, Point p, Point q) { Circle circ = MakeDiameter(p, q); Circle left = Circle.INVALID; Circle right = Circle.INVALID; // For each point not in the two-point circle Point pq = q.Subtract(p); foreach (Point r in points) { if (circ.Contains(r)) continue; // Form a circumcircle and classify it on left or right side double cross = pq.Cross(r.Subtract(p)); Circle c = MakeCircumcircle(p, q, r); if (c.r < 0) continue; else if (cross > 0 && (left.r < 0 || pq.Cross(c.c.Subtract(p)) > pq.Cross(left.c.Subtract(p)))) left = c; else if (cross < 0 && (right.r < 0 || pq.Cross(c.c.Subtract(p)) < pq.Cross(right.c.Subtract(p)))) right = c; } // Select which circle to return if (left.r < 0 && right.r < 0) return circ; else if (left.r < 0) return right; else if (right.r < 0) return left; else return left.r <= right.r ? left : right; } public static Circle MakeDiameter(Point a, Point b) { Point c = new Point((a.x + b.x) / 2, (a.y + b.y) / 2); return new Circle(c, Math.Max(c.Distance(a), c.Distance(b))); } public static Circle MakeCircumcircle(Point a, Point b, Point c) { // Mathematical algorithm from Wikipedia: Circumscribed circle double ox = (Math.Min(Math.Min(a.x, b.x), c.x) + Math.Max(Math.Max(a.x, b.x), c.x)) / 2; double oy = (Math.Min(Math.Min(a.y, b.y), c.y) + Math.Max(Math.Max(a.y, b.y), c.y)) / 2; double ax = a.x - ox, ay = a.y - oy; double bx = b.x - ox, by = b.y - oy; double cx = c.x - ox, cy = c.y - oy; double d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2; if (d == 0) return Circle.INVALID; double x = ((ax*ax + ay*ay) * (by - cy) + (bx*bx + by*by) * (cy - ay) + (cx*cx + cy*cy) * (ay - by)) / d; double y = ((ax*ax + ay*ay) * (cx - bx) + (bx*bx + by*by) * (ax - cx) + (cx*cx + cy*cy) * (bx - ax)) / d; Point p = new Point(ox + x, oy + y); double r = Math.Max(Math.Max(p.Distance(a), p.Distance(b)), p.Distance(c)); return new Circle(p, r); } } public struct Circle { public static readonly Circle INVALID = new Circle(new Point(0, 0), -1); private const double MULTIPLICATIVE_EPSILON = 1 + 1e-14; public Point c; // Center public double r; // Radius public Circle(Point c, double r) { this.c = c; this.r = r; } public bool Contains(Point p) { return c.Distance(p) <= r * MULTIPLICATIVE_EPSILON; } public bool Contains(ICollection ps) { foreach (Point p in ps) { if (!Contains(p)) return false; } return true; } } public struct Point { public double x; public double y; public Point(double x, double y) { this.x = x; this.y = y; } public Point Subtract(Point p) { return new Point(x - p.x, y - p.y); } public double Distance(Point p) { double dx = x - p.x; double dy = y - p.y; return Math.Sqrt(dx * dx + dy * dy); } // Signed area / determinant thing public double Cross(Point p) { return x * p.y - y * p.x; } }