1# Copyright (c) 2013 The Chromium OS Authors. All rights reserved. 2# Use of this source code is governed by a BSD-style license that can be 3# found in the LICENSE file. 4 5"""minicircle: calculating the minimal enclosing circle given a set of points 6 7 Reference: 8 [1] Emo Welzl. Smallest enclosing disks (balls and ellipsoids). 9 http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.46.1450 10 [2] Circumscribed circle. http://en.wikipedia.org/wiki/Circumscribed_circle 11 12 - get_two_farthest_clusters(): Classify the points into two farthest clusters 13 14 - get_radii_of_two_minimal_enclosing_circles(): Get the radii of the 15 two minimal enclosing circles 16 17""" 18 19 20import copy 21 22from sets import Set 23 24from elements import Point, Circle 25 26 27def _mini_circle_2_points(p1, p2): 28 """Derive the mini circle with p1 and p2 composing the diameter. 29 This also handles the special case when p1 and p2 are located at 30 the same coordinate. 31 """ 32 center = Point((p1.x + p2.x) * 0.5, (p1.y + p2.y) * 0.5) 33 radius = center.distance(p1) 34 return Circle(center, radius) 35 36 37def _mini_circle_3_points(A, B, C): 38 """Derive the mini circle enclosing arbitrary three points, A, B, C. 39 40 @param A: a point (possibly a vertex of a triangle) 41 @param B: a point (possibly a vertex of a triangle) 42 @param C: a point (possibly a vertex of a triangle) 43 """ 44 # Check if this is an obtuse triangle or a right triangle, 45 # including the special cases 46 # (1) the 3 points are on the same line 47 # (2) any 2 points are located at the same coordinate 48 # (3) all 3 points are located at the same coordinate 49 a = B.distance(C) 50 b = C.distance(A) 51 c = A.distance(B) 52 if (a ** 2 >= b ** 2 + c ** 2): 53 return _mini_circle_2_points(B, C) 54 elif (b ** 2 >= c ** 2 + a ** 2): 55 return _mini_circle_2_points(C, A) 56 elif (c ** 2 >= a ** 2 + b ** 2): 57 return _mini_circle_2_points(A, B) 58 59 # It is an acute triangle. Refer to Reference [2]. 60 D = 2 * (A.x * (B.y - C.y) + B.x *(C.y - A.y) + C.x * (A.y - B.y)) 61 x = ((A.x ** 2 + A.y ** 2) * (B.y - C.y) + 62 (B.x ** 2 + B.y ** 2) * (C.y - A.y) + 63 (C.x ** 2 + C.y ** 2) * (A.y - B.y)) / D 64 y = ((A.x ** 2 + A.y ** 2) * (C.x - B.x) + 65 (B.x ** 2 + B.y ** 2) * (A.x - C.x) + 66 (C.x ** 2 + C.y ** 2) * (B.x - A.x)) / D 67 68 center = Point(x, y) 69 radius = center.distance(A) 70 return Circle(center, radius) 71 72 73def _b_minicircle0(R): 74 """build minicircle0: build the mini circle with an empty P and has R 75 on the boundary. 76 77 @param R: boundary points, a set of points which should be on the boundary 78 of the circle to be built 79 """ 80 if len(R) == 0: 81 return Circle(None, None) 82 if len(R) == 1: 83 return Circle(R.pop(), 0) 84 elif len(R) == 2: 85 p1 = R.pop() 86 p2 = R.pop() 87 return _mini_circle_2_points(p1, p2) 88 else: 89 return _mini_circle_3_points(*list(R)) 90 91 92def _b_minicircle(P, R): 93 """build minicircle: build the mini circle enclosing P and has R on the 94 boundary. 95 96 @param P: a set of points that should be enclosed in the circle to be built 97 @param R: boundary points, a set of points which should be on the boundary 98 of the circle to be built 99 """ 100 P = copy.deepcopy(P) 101 R = copy.deepcopy(R) 102 if len(P) == 0 or len(R) == 3: 103 D = _b_minicircle0(R) 104 else: 105 p = P.pop() 106 D = _b_minicircle(P, R) 107 if (not D) or (p not in D): 108 R.add(p) 109 D = _b_minicircle(P, R) 110 return D 111 112 113def _make_Set_of_Points(points): 114 """Convert the points to a set of Point objects. 115 116 @param points: could be a list/set of pairs_of_ints/Point_objects. 117 """ 118 return Set([p if isinstance(p, Point) else Point(*p) for p in points]) 119 120 121def minicircle(points): 122 """Get the minimal enclosing circle of the points. 123 124 @param points: a list of points which should be enclosed in the circle to be 125 built 126 """ 127 P = _make_Set_of_Points(points) 128 return _b_minicircle(P, Set()) if P else None 129