# Copyright (c) 2013 The Chromium OS Authors. All rights reserved. # Use of this source code is governed by a BSD-style license that can be # found in the LICENSE file. """minicircle: calculating the minimal enclosing circle given a set of points Reference: [1] Emo Welzl. Smallest enclosing disks (balls and ellipsoids). http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.46.1450 [2] Circumscribed circle. http://en.wikipedia.org/wiki/Circumscribed_circle - get_two_farthest_clusters(): Classify the points into two farthest clusters - get_radii_of_two_minimal_enclosing_circles(): Get the radii of the two minimal enclosing circles """ import copy from sets import Set from elements import Point, Circle def _mini_circle_2_points(p1, p2): """Derive the mini circle with p1 and p2 composing the diameter. This also handles the special case when p1 and p2 are located at the same coordinate. """ center = Point((p1.x + p2.x) * 0.5, (p1.y + p2.y) * 0.5) radius = center.distance(p1) return Circle(center, radius) def _mini_circle_3_points(A, B, C): """Derive the mini circle enclosing arbitrary three points, A, B, C. @param A: a point (possibly a vertex of a triangle) @param B: a point (possibly a vertex of a triangle) @param C: a point (possibly a vertex of a triangle) """ # Check if this is an obtuse triangle or a right triangle, # including the special cases # (1) the 3 points are on the same line # (2) any 2 points are located at the same coordinate # (3) all 3 points are located at the same coordinate a = B.distance(C) b = C.distance(A) c = A.distance(B) if (a ** 2 >= b ** 2 + c ** 2): return _mini_circle_2_points(B, C) elif (b ** 2 >= c ** 2 + a ** 2): return _mini_circle_2_points(C, A) elif (c ** 2 >= a ** 2 + b ** 2): return _mini_circle_2_points(A, B) # It is an acute triangle. Refer to Reference [2]. D = 2 * (A.x * (B.y - C.y) + B.x *(C.y - A.y) + C.x * (A.y - B.y)) x = ((A.x ** 2 + A.y ** 2) * (B.y - C.y) + (B.x ** 2 + B.y ** 2) * (C.y - A.y) + (C.x ** 2 + C.y ** 2) * (A.y - B.y)) / D y = ((A.x ** 2 + A.y ** 2) * (C.x - B.x) + (B.x ** 2 + B.y ** 2) * (A.x - C.x) + (C.x ** 2 + C.y ** 2) * (B.x - A.x)) / D center = Point(x, y) radius = center.distance(A) return Circle(center, radius) def _b_minicircle0(R): """build minicircle0: build the mini circle with an empty P and has R on the boundary. @param R: boundary points, a set of points which should be on the boundary of the circle to be built """ if len(R) == 0: return Circle(None, None) if len(R) == 1: return Circle(R.pop(), 0) elif len(R) == 2: p1 = R.pop() p2 = R.pop() return _mini_circle_2_points(p1, p2) else: return _mini_circle_3_points(*list(R)) def _b_minicircle(P, R): """build minicircle: build the mini circle enclosing P and has R on the boundary. @param P: a set of points that should be enclosed in the circle to be built @param R: boundary points, a set of points which should be on the boundary of the circle to be built """ P = copy.deepcopy(P) R = copy.deepcopy(R) if len(P) == 0 or len(R) == 3: D = _b_minicircle0(R) else: p = P.pop() D = _b_minicircle(P, R) if (not D) or (p not in D): R.add(p) D = _b_minicircle(P, R) return D def _make_Set_of_Points(points): """Convert the points to a set of Point objects. @param points: could be a list/set of pairs_of_ints/Point_objects. """ return Set([p if isinstance(p, Point) else Point(*p) for p in points]) def minicircle(points): """Get the minimal enclosing circle of the points. @param points: a list of points which should be enclosed in the circle to be built """ P = _make_Set_of_Points(points) return _b_minicircle(P, Set()) if P else None