-
Notifications
You must be signed in to change notification settings - Fork 1
/
smallestenclosingcircle.py
126 lines (105 loc) · 4.14 KB
/
smallestenclosingcircle.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
#
# Smallest enclosing circle - Library (Python)
#
# 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 <http://www.gnu.org/licenses/>.
#
import math, random
# Data conventions: A point is a pair of floats (x, y). A circle is a triple of floats (center x, center y, radius).
# Returns the smallest circle that encloses all the given points. Runs in expected O(n) time, randomized.
# Input: A sequence of pairs of floats or ints, e.g. [(0,5), (3.1,-2.7)].
# Output: A triple of floats representing a circle.
# Note: If 0 points are given, None is returned. If 1 point is given, a circle of radius 0 is returned.
#
# Initially: No boundary points known
def make_circle(points):
# Convert to float and randomize order
shuffled = [(float(x), float(y)) for (x, y) in points]
random.shuffle(shuffled)
# Progressively add points to circle or recompute circle
c = None
for (i, p) in enumerate(shuffled):
if c is None or not is_in_circle(c, p):
c = _make_circle_one_point(shuffled[ : i + 1], p)
return c
# One boundary point known
def _make_circle_one_point(points, p):
c = (p[0], p[1], 0.0)
for (i, q) in enumerate(points):
if not is_in_circle(c, q):
if c[2] == 0.0:
c = make_diameter(p, q)
else:
c = _make_circle_two_points(points[ : i + 1], p, q)
return c
# Two boundary points known
def _make_circle_two_points(points, p, q):
circ = make_diameter(p, q)
left = None
right = None
px, py = p
qx, qy = q
# For each point not in the two-point circle
for r in points:
if is_in_circle(circ, r):
continue
# Form a circumcircle and classify it on left or right side
cross = _cross_product(px, py, qx, qy, r[0], r[1])
c = make_circumcircle(p, q, r)
if c is None:
continue
elif cross > 0.0 and (left is None or _cross_product(px, py, qx, qy, c[0], c[1]) > _cross_product(px, py, qx, qy, left[0], left[1])):
left = c
elif cross < 0.0 and (right is None or _cross_product(px, py, qx, qy, c[0], c[1]) < _cross_product(px, py, qx, qy, right[0], right[1])):
right = c
# Select which circle to return
if left is None and right is None:
return circ
elif left is None:
return right
elif right is None:
return left
else:
return left if (left[2] <= right[2]) else right
def make_diameter(a, b):
cx = (a[0] + b[0]) / 2
cy = (a[1] + b[1]) / 2
r0 = math.hypot(cx - a[0], cy - a[1])
r1 = math.hypot(cx - b[0], cy - b[1])
return (cx, cy, max(r0, r1))
def make_circumcircle(a, b, c):
# Mathematical algorithm from Wikipedia: Circumscribed circle
ox = (min(a[0], b[0], c[0]) + max(a[0], b[0], c[0])) / 2
oy = (min(a[1], b[1], c[1]) + max(a[1], b[1], c[1])) / 2
ax = a[0] - ox; ay = a[1] - oy
bx = b[0] - ox; by = b[1] - oy
cx = c[0] - ox; cy = c[1] - oy
d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2.0
if d == 0.0:
return None
x = ox + ((ax*ax + ay*ay) * (by - cy) + (bx*bx + by*by) * (cy - ay) + (cx*cx + cy*cy) * (ay - by)) / d
y = oy + ((ax*ax + ay*ay) * (cx - bx) + (bx*bx + by*by) * (ax - cx) + (cx*cx + cy*cy) * (bx - ax)) / d
ra = math.hypot(x - a[0], y - a[1])
rb = math.hypot(x - b[0], y - b[1])
rc = math.hypot(x - c[0], y - c[1])
return (x, y, max(ra, rb, rc))
_MULTIPLICATIVE_EPSILON = 1 + 1e-14
def is_in_circle(c, p):
return c is not None and math.hypot(p[0] - c[0], p[1] - c[1]) <= c[2] * _MULTIPLICATIVE_EPSILON
# Returns twice the signed area of the triangle defined by (x0, y0), (x1, y1), (x2, y2).
def _cross_product(x0, y0, x1, y1, x2, y2):
return (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)