Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #
- # 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
- from numba.pycc import CC
- cc = CC('smallest_circle')
- cc.verbose = True
- # 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
- @cc.export('is_in_circle', 'b1(f8[:], f8[:])')
- def is_in_circle(c, p):
- if (c is not None) and (math.hypot(p[0] - c[0], p[1] - c[1]) <= c[2] * _MULTIPLICATIVE_EPSILON):
- res = True
- else:
- res = False
- return res
- @cc.export('make_circle', '(f8[:,:],)')
- def make_circle(points):
- # points : numpy.ndarray (3,2) float64
- # Convert to float and randomize order
- shuffled = points[np.random.permutation(len(points))]
- # Progressively add points to circle or recompute circle
- first_pass = True
- c = (0,0,0)
- for (i, p) in enumerate(shuffled):
- if first_pass is True or not is_in_circle(c, p):
- first_pass = False
- c = _make_circle_one_point(shuffled[: i + 1], p)
- return c
- # One boundary point known
- @cc.export('make_circle_one_point', '(f8[:,:], f8[:])')
- def _make_circle_one_point(points, p):
- # points : numpy.ndarray (n,2) float64
- # p : numpy.ndarray (2,) float64
- 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
- @cc.export('make_circle_two_points', '(f8[:,:],f8[:],f8[:])')
- def _make_circle_two_points(points, p, q):
- # points : numpy.ndarray (n,2) float64
- # p : numpy.ndarray (2,) float64
- # q : numpy.ndarray (2,) float64
- circ = make_diameter(p, q)
- left = (0, 0, 0)
- left_valid = False
- right = (0, 0, 0)
- right_valid = False
- 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, c_valid = make_circumcircle(p, q, r)
- if c_valid is False:
- continue
- elif cross > 0.0 and (
- left_valid is False or
- _cross_product(px, py, qx, qy, c[0], c[1]) > _cross_product(px, py, qx, qy, left[0], left[1])):
- left = c
- left_valid = True
- elif cross < 0.0 and (right_valid is False or _cross_product(px, py, qx, qy, c[0], c[1]) < _cross_product(px, py, qx, qy, right[0], right[1])):
- right = c
- right_valid = True
- # Select which circle to return
- if left_valid is False and right_valid is False:
- return circ
- elif left_valid is False:
- return right
- elif right_valid is False:
- return left
- else:
- return left if (left[2] <= right[2]) else right
- @cc.export('make_diameter', '(f8[:],f8[:])')
- def make_diameter(a, b):
- # a : numpy.ndarray (2,) float64
- # b : numpy.ndarray (2,) float64
- 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))
- @cc.export('make_circumcircle', '(f8[:], f8[:], f8[:])')
- def make_circumcircle(a, b, c):
- # Mathematical algorithm from Wikipedia: Circumscribed circle
- # a : numpy.ndarray (2,) float64
- # b : numpy.ndarray (2,) float64
- # c : numpy.ndarray (2,) float64
- 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 (0,0,0), False
- 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)), True
- _MULTIPLICATIVE_EPSILON = 1 + 1e-14
- # Returns twice the signed area of the triangle defined by (x0, y0), (x1, y1), (x2, y2).
- @cc.export('cross_product', 'f8(f8,f8,f8,f8,f8,f8)')
- def _cross_product(x0, y0, x1, y1, x2, y2):
- return (x1 - x0) * (y2 - y0) - (y1 - y0) * (x2 - x0)
- if __name__ == '__main__':
- cc.compile()
- import matplotlib.pyplot as plt
- import matplotlib.patches as patches
- import numpy as np
- # Create 30 random points
- points = np.random.rand(30, 2)
- # Compute the smallest enclosing circle
- c = make_circle(points)
- # Plot the points and the enclosing circle
- fig, ax = plt.subplots()
- ax.set_aspect('equal')
- ax.add_patch(patches.Circle(c[:2], c[2], fill=False))
- ax.scatter(points[:, 0], points[:, 1])
- plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement