/* * Smallest enclosing circle - Library (C++) * * Copyright (c) 2021 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 . */ #include #include #include #include #include "SmallestEnclosingCircle.hpp" using std::size_t; using std::vector; using std::max; using std::min; /*---- Members of struct Point ----*/ Point Point::subtract(const Point &p) const { return Point{x - p.x, y - p.y}; } double Point::distance(const Point &p) const { return std::hypot(x - p.x, y - p.y); } double Point::cross(const Point &p) const { return x * p.y - y * p.x; } /*---- Members of struct Circle ----*/ const Circle Circle::INVALID{Point{0, 0}, -1}; const double Circle::MULTIPLICATIVE_EPSILON = 1 + 1e-14; bool Circle::contains(const Point &p) const { return c.distance(p) <= r * MULTIPLICATIVE_EPSILON; } bool Circle::contains(const vector &ps) const { for (const Point &p : ps) { if (!contains(p)) return false; } return true; } /*---- Smallest enclosing circle algorithm ----*/ static Circle makeSmallestEnclosingCircleOnePoint (const vector &points, size_t end, const Point &p); static Circle makeSmallestEnclosingCircleTwoPoints(const vector &points, size_t end, const Point &p, const Point &q); static std::default_random_engine randGen((std::random_device())()); // Initially: No boundary points known Circle makeSmallestEnclosingCircle(vector points) { // Randomize order std::shuffle(points.begin(), points.end(), randGen); // Progressively add points to circle or recompute circle Circle c = Circle::INVALID; for (size_t i = 0; i < points.size(); i++) { const Point &p = points.at(i); if (c.r < 0 || !c.contains(p)) c = makeSmallestEnclosingCircleOnePoint(points, i + 1, p); } return c; } // One boundary point known static Circle makeSmallestEnclosingCircleOnePoint(const vector &points, size_t end, const Point &p) { Circle c{p, 0}; for (size_t i = 0; i < end; i++) { const Point &q = points.at(i); if (!c.contains(q)) { if (c.r == 0) c = makeDiameter(p, q); else c = makeSmallestEnclosingCircleTwoPoints(points, i + 1, p, q); } } return c; } // Two boundary points known static Circle makeSmallestEnclosingCircleTwoPoints(const vector &points, size_t end, const Point &p, const 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); for (size_t i = 0; i < end; i++) { const Point &r = points.at(i); 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; } Circle makeDiameter(const Point &a, const Point &b) { Point c{(a.x + b.x) / 2, (a.y + b.y) / 2}; return Circle{c, max(c.distance(a), c.distance(b))}; } Circle makeCircumcircle(const Point &a, const Point &b, const Point &c) { // Mathematical algorithm from Wikipedia: Circumscribed circle double ox = (min(min(a.x, b.x), c.x) + max(max(a.x, b.x), c.x)) / 2; double oy = (min(min(a.y, b.y), c.y) + max(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{ox + x, oy + y}; double r = max(max(p.distance(a), p.distance(b)), p.distance(c)); return Circle{p, r}; }