|
1 | | -# |
2 | | -# Smallest enclosing circle - Library (Python) |
3 | | -# |
4 | | -# Copyright (c) 2020 Project Nayuki |
5 | | -# https://www.nayuki.io/page/smallest-enclosing-circle |
6 | | -# |
7 | | -# This program is free software: you can redistribute it and/or modify |
8 | | -# it under the terms of the GNU Lesser General Public License as published by |
9 | | -# the Free Software Foundation, either version 3 of the License, or |
10 | | -# (at your option) any later version. |
11 | | -# |
12 | | -# This program is distributed in the hope that it will be useful, |
13 | | -# but WITHOUT ANY WARRANTY; without even the implied warranty of |
14 | | -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
15 | | -# GNU Lesser General Public License for more details. |
16 | | -# |
17 | | -# You should have received a copy of the GNU Lesser General Public License |
18 | | -# along with this program (see COPYING.txt and COPYING.LESSER.txt). |
19 | | -# If not, see <http://www.gnu.org/licenses/>. |
20 | | -# |
21 | | - |
22 | | -# fmt: off |
23 | | -# isort: skip_file |
24 | | -import math, random |
| 1 | +""" |
| 2 | +Smallest enclosing cylinder computation. |
| 3 | +
|
| 4 | +Based on the algorithm from: |
| 5 | +https://www.nayuki.io/page/smallest-enclosing-circle |
| 6 | +""" |
| 7 | + |
| 8 | +import math |
| 9 | +from typing import List, Optional, Tuple |
| 10 | + |
25 | 11 | import numpy as np |
26 | 12 |
|
27 | 13 |
|
28 | | -# Data conventions: A point is a pair of floats (x, y). A circle is a triple of floats (center x, center y, radius). |
29 | | - |
30 | | -# Returns the smallest circle that encloses all the given points. Runs in expected O(n) time, randomized. |
31 | | -# Input: A sequence of pairs of floats or ints, e.g. [(0,5), (3.1,-2.7)]. |
32 | | -# Output: A triple of floats representing a circle. |
33 | | -# Note: If 0 points are given, None is returned. If 1 point is given, a circle of radius 0 is returned. |
34 | | -# |
35 | | -# Initially: No boundary points known |
36 | | -def make_circle(points): |
37 | | - # Convert to float and randomize order |
38 | | - shuffled = [(float(x), float(y)) for (x, y) in points] |
39 | | - random.shuffle(shuffled) |
40 | | - |
41 | | - # Progressively add points to circle or recompute circle |
42 | | - c = None |
43 | | - for (i, p) in enumerate(shuffled): |
44 | | - if c is None or not is_in_circle(c, p): |
45 | | - c = _make_circle_one_point(shuffled[ : i + 1], p) |
46 | | - return c |
47 | | - |
48 | | - |
49 | | -# One boundary point known |
50 | | -def _make_circle_one_point(points, p): |
51 | | - c = (p[0], p[1], 0.0) |
52 | | - for (i, q) in enumerate(points): |
53 | | - if not is_in_circle(c, q): |
54 | | - if c[2] == 0.0: |
55 | | - c = make_diameter(p, q) |
| 14 | +def _compute_smallest_circle( |
| 15 | + points: List[Tuple[float, float]] |
| 16 | +) -> Optional[Tuple[float, float, float]]: |
| 17 | + points = [(float(x), float(y)) for x, y in points] |
| 18 | + np.random.shuffle(points) |
| 19 | + |
| 20 | + circle = None |
| 21 | + for i, point in enumerate(points): |
| 22 | + if circle is None or not _point_in_circle(circle, point): |
| 23 | + circle = _compute_circle_with_point(points[: i + 1], point) |
| 24 | + return circle |
| 25 | + |
| 26 | + |
| 27 | +def _compute_circle_with_point( |
| 28 | + points: List[Tuple[float, float]], p: Tuple[float, float] |
| 29 | +) -> Tuple[float, float, float]: |
| 30 | + circle = (p[0], p[1], 0.0) |
| 31 | + for i, q in enumerate(points): |
| 32 | + if not _point_in_circle(circle, q): |
| 33 | + if circle[2] == 0.0: |
| 34 | + circle = _get_circle_from_diameter(p, q) |
56 | 35 | else: |
57 | | - c = _make_circle_two_points(points[ : i + 1], p, q) |
58 | | - return c |
| 36 | + circle = _compute_circle_with_two_points(points[: i + 1], p, q) |
| 37 | + return circle |
59 | 38 |
|
60 | 39 |
|
61 | | -# Two boundary points known |
62 | | -def _make_circle_two_points(points, p, q): |
63 | | - circ = make_diameter(p, q) |
64 | | - left = None |
65 | | - right = None |
66 | | - px, py = p |
67 | | - qx, qy = q |
| 40 | +def _compute_circle_with_two_points( |
| 41 | + points: List[Tuple[float, float]], p: Tuple[float, float], q: Tuple[float, float] |
| 42 | +) -> Tuple[float, float, float]: |
| 43 | + circle = _get_circle_from_diameter(p, q) |
| 44 | + left = right = None |
68 | 45 |
|
69 | | - # For each point not in the two-point circle |
70 | 46 | for r in points: |
71 | | - if is_in_circle(circ, r): |
| 47 | + if _point_in_circle(circle, r): |
72 | 48 | continue |
73 | 49 |
|
74 | | - # Form a circumcircle and classify it on left or right side |
75 | | - cross = _cross_product(px, py, qx, qy, r[0], r[1]) |
76 | | - c = make_circumcircle(p, q, r) |
77 | | - if c is None: |
| 50 | + cross = _compute_cross_product(p[0], p[1], q[0], q[1], r[0], r[1]) |
| 51 | + candidate = _compute_circumcircle(p, q, r) |
| 52 | + |
| 53 | + if candidate is None: |
78 | 54 | continue |
79 | | - 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])): |
80 | | - left = c |
81 | | - 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])): |
82 | | - right = c |
| 55 | + elif cross > 0 and ( |
| 56 | + left is None |
| 57 | + or _compute_cross_product( |
| 58 | + p[0], p[1], q[0], q[1], candidate[0], candidate[1] |
| 59 | + ) |
| 60 | + > _compute_cross_product(p[0], p[1], q[0], q[1], left[0], left[1]) |
| 61 | + ): |
| 62 | + left = candidate |
| 63 | + elif cross < 0 and ( |
| 64 | + right is None |
| 65 | + or _compute_cross_product( |
| 66 | + p[0], p[1], q[0], q[1], candidate[0], candidate[1] |
| 67 | + ) |
| 68 | + < _compute_cross_product(p[0], p[1], q[0], q[1], right[0], right[1]) |
| 69 | + ): |
| 70 | + right = candidate |
83 | 71 |
|
84 | | - # Select which circle to return |
85 | 72 | if left is None and right is None: |
86 | | - return circ |
| 73 | + return circle |
87 | 74 | elif left is None: |
88 | 75 | return right |
89 | 76 | elif right is None: |
90 | 77 | return left |
91 | 78 | else: |
92 | | - return left if (left[2] <= right[2]) else right |
93 | | - |
94 | | - |
95 | | -def make_diameter(a, b): |
96 | | - cx = (a[0] + b[0]) / 2 |
97 | | - cy = (a[1] + b[1]) / 2 |
98 | | - r0 = math.hypot(cx - a[0], cy - a[1]) |
99 | | - r1 = math.hypot(cx - b[0], cy - b[1]) |
100 | | - return (cx, cy, max(r0, r1)) |
101 | | - |
102 | | - |
103 | | -def make_circumcircle(a, b, c): |
104 | | - # Mathematical algorithm from Wikipedia: Circumscribed circle |
105 | | - ox = (min(a[0], b[0], c[0]) + max(a[0], b[0], c[0])) / 2 |
106 | | - oy = (min(a[1], b[1], c[1]) + max(a[1], b[1], c[1])) / 2 |
107 | | - ax = a[0] - ox; ay = a[1] - oy |
108 | | - bx = b[0] - ox; by = b[1] - oy |
109 | | - cx = c[0] - ox; cy = c[1] - oy |
110 | | - d = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2.0 |
111 | | - if d == 0.0: |
| 79 | + return left if left[2] <= right[2] else right |
| 80 | + |
| 81 | + |
| 82 | +def _get_circle_from_diameter( |
| 83 | + a: Tuple[float, float], b: Tuple[float, float] |
| 84 | +) -> Tuple[float, float, float]: |
| 85 | + center_x = (a[0] + b[0]) / 2 |
| 86 | + center_y = (a[1] + b[1]) / 2 |
| 87 | + radius = max( |
| 88 | + math.hypot(center_x - a[0], center_y - a[1]), |
| 89 | + math.hypot(center_x - b[0], center_y - b[1]), |
| 90 | + ) |
| 91 | + return (center_x, center_y, radius) |
| 92 | + |
| 93 | + |
| 94 | +def _compute_circumcircle( |
| 95 | + a: Tuple[float, float], b: Tuple[float, float], c: Tuple[float, float] |
| 96 | +) -> Optional[Tuple[float, float, float]]: |
| 97 | + # Center point |
| 98 | + mid_x = (min(a[0], b[0], c[0]) + max(a[0], b[0], c[0])) / 2 |
| 99 | + mid_y = (min(a[1], b[1], c[1]) + max(a[1], b[1], c[1])) / 2 |
| 100 | + |
| 101 | + ax, ay = a[0] - mid_x, a[1] - mid_y |
| 102 | + bx, by = b[0] - mid_x, b[1] - mid_y |
| 103 | + cx, cy = c[0] - mid_x, c[1] - mid_y |
| 104 | + |
| 105 | + det = (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by)) * 2 |
| 106 | + if det == 0: |
112 | 107 | return None |
113 | | - x = ox + ((ax*ax + ay*ay) * (by - cy) + (bx*bx + by*by) * (cy - ay) + (cx*cx + cy*cy) * (ay - by)) / d |
114 | | - y = oy + ((ax*ax + ay*ay) * (cx - bx) + (bx*bx + by*by) * (ax - cx) + (cx*cx + cy*cy) * (bx - ax)) / d |
115 | | - ra = math.hypot(x - a[0], y - a[1]) |
116 | | - rb = math.hypot(x - b[0], y - b[1]) |
117 | | - rc = math.hypot(x - c[0], y - c[1]) |
118 | | - return (x, y, max(ra, rb, rc)) |
119 | | - |
120 | | - |
121 | | -_MULTIPLICATIVE_EPSILON = 1 + 1e-14 |
122 | 108 |
|
123 | | -def is_in_circle(c, p): |
124 | | - return c is not None and math.hypot(p[0] - c[0], p[1] - c[1]) <= c[2] * _MULTIPLICATIVE_EPSILON |
| 109 | + center_x = ( |
| 110 | + mid_x |
| 111 | + + ( |
| 112 | + (ax * ax + ay * ay) * (by - cy) |
| 113 | + + (bx * bx + by * by) * (cy - ay) |
| 114 | + + (cx * cx + cy * cy) * (ay - by) |
| 115 | + ) |
| 116 | + / det |
| 117 | + ) |
| 118 | + center_y = ( |
| 119 | + mid_y |
| 120 | + + ( |
| 121 | + (ax * ax + ay * ay) * (cx - bx) |
| 122 | + + (bx * bx + by * by) * (ax - cx) |
| 123 | + + (cx * cx + cy * cy) * (bx - ax) |
| 124 | + ) |
| 125 | + / det |
| 126 | + ) |
| 127 | + radius = max(math.hypot(center_x - p[0], center_y - p[1]) for p in (a, b, c)) |
| 128 | + |
| 129 | + return (center_x, center_y, radius) |
| 130 | + |
| 131 | + |
| 132 | +def _point_in_circle( |
| 133 | + circle: Optional[Tuple[float, float, float]], point: Tuple[float, float] |
| 134 | +) -> bool: |
| 135 | + if circle is None: |
| 136 | + return False |
| 137 | + return math.hypot(point[0] - circle[0], point[1] - circle[1]) <= circle[2] * ( |
| 138 | + 1 + 1e-14 |
| 139 | + ) |
| 140 | + |
| 141 | + |
| 142 | +def _compute_cross_product( |
| 143 | + x0: float, y0: float, x1: float, y1: float, x2: float, y2: float |
| 144 | +) -> float: |
| 145 | + return (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0) |
125 | 146 |
|
126 | 147 |
|
127 | | -# Returns twice the signed area of the triangle defined by (x0, y0), (x1, y1), (x2, y2). |
128 | | -def _cross_product(x0, y0, x1, y1, x2, y2): |
129 | | - return (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0) |
| 148 | +def aabc(points: np.ndarray) -> Tuple[float, float, float, float, float]: |
| 149 | + """ |
| 150 | + Compute axis-aligned bounding cylinder for 3D points. |
130 | 151 |
|
| 152 | + Args: |
| 153 | + points: Nx3 array of points |
131 | 154 |
|
132 | | -def aabc(points): |
| 155 | + Returns: |
| 156 | + (center_x, center_y, radius, min_z, max_z) tuple |
| 157 | + """ |
133 | 158 | points = np.asarray(points) |
134 | | - z_min = points[:, 2].min() |
135 | | - z_max = points[:, 2].max() |
136 | | - x, y, r = make_circle(points[:, :2]) |
137 | | - return x, y, r, z_min, z_max |
| 159 | + z_bounds = points[:, 2].min(), points[:, 2].max() |
| 160 | + center_x, center_y, radius = _compute_smallest_circle(points[:, :2]) |
| 161 | + return center_x, center_y, radius, z_bounds[0], z_bounds[1] |
0 commit comments