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