Coverage for src / min_circle / welzl.py: 100%

1"""Welzl's algorithm implementation for finding the minimum enclosing circle. 

2 

3This module provides an implementation of Welzl's randomized algorithm 

4for finding the minimum enclosing circle of a set of points in 2D space. 

5""" 

6 

7import secrets 

8 

9import numpy as np 

10 

11from .cvx import min_circle_cvx 

12from .utils.circle import Circle 

13 

14 

15def perpendicular_slope(p1: np.ndarray, p2: np.ndarray) -> float: 

16 """Calculate the slope of the perpendicular bisector between two points. 

17 

18 Args: 

19 p1: First point as a numpy array [x, y] 

20 p2: Second point as a numpy array [x, y] 

21 

22 Returns: 

23 The slope of the perpendicular bisector line 

24 

25 Examples: 

26 >>> import numpy as np 

27 >>> perpendicular_slope(np.array([0.0, 0.0]), np.array([2.0, 2.0])) 

28 -1.0 

29 >>> perpendicular_slope(np.array([0.0, 0.0]), np.array([1.0, 0.0])) 

30 inf 

31 >>> perpendicular_slope(np.array([0.0, 0.0]), np.array([2.0, 1.0])) 

32 -2.0 

33 """ 

34 if p2[1] == p1[1]: # horizontal line 

35 return float(np.inf) # Perpendicular bisector is vertical 

36 return float(-(p2[0] - p1[0]) / (p2[1] - p1[1])) 

37 

38 

39def make_circle_n_points(matrix: list[np.ndarray]) -> Circle: 

40 """Construct a circle from 0 to 3 points. 

41 

42 Args: 

43 matrix: List of points (up to 3) as numpy arrays 

44 

45 Returns: 

46 Circle object containing the center and radius of the circle 

47 """ 

48 num_points = len(matrix) 

49 

50 if num_points == 0: 

51 return Circle(np.array([0.0, 0.0]), -np.inf) 

52 

53 if num_points == 1: 

54 return Circle(matrix[0], 0) # Single point, radius 0 

55 

56 if num_points == 2 or num_points == 3: 

57 return min_circle_cvx(np.array(matrix), solver="CLARABEL") 

58 

59 msg = f"Expected 0-3 points, got {num_points}" 

60 raise ValueError(msg) 

61 

62 # if num_points == 2: 

63 # # Two points: the center is the midpoint, radius is half the distance 

64 # center = (matrix[0] + matrix[1]) / 2 

65 # radius = np.linalg.norm(matrix[0] - matrix[1]) / 2 

66 # 

67 # return Circle(center=center, radius=radius) 

68 

69 # if num_points == 3: 

70 # # For 3 points: use the circumcenter and circumradius formula 

71 # # p1, p2, p3 = matrix 

72 # p = np.array(matrix) 

73 # 

74 # # Midpoints of the sides 

75 # mid12 = (p[0] + p[1]) * 0.5 

76 # mid23 = (p[1] + p[2]) * 0.5 

77 # 

78 # # Slopes of the perpendicular bisectors 

79 # # Perpendicular slope to line 1-2: m1 

80 # m1 = perpendicular_slope(p[1], p[0]) 

81 # 

82 # # Perpendicular slope to line 2-3: m2 

83 # m2 = perpendicular_slope(p[1], p[2]) 

84 # 

85 # # Use line equations to solve for the intersection (circumcenter) 

86 # if m1 == np.inf: # Line 1-2 is vertical, so we solve for x = mid12[0] 

87 # center_x = mid12[0] 

88 # center_y = m2 * (center_x - mid23[0]) + mid23[1] 

89 # elif m2 == np.inf: # Line 2-3 is vertical, so we solve for x = mid23[0] 

90 # center_x = mid23[0] 

91 # center_y = m1 * (center_x - mid12[0]) + mid12[1] 

92 # else: 

93 # # Calculate the intersection of the two perpendicular bisectors 

94 # A, B = m1, -1 

95 # C, D = m2, -1 

96 # E = m1 * mid12[0] - mid12[1] 

97 # F = m2 * mid23[0] - mid23[1] 

98 # 

99 # # Construct the coefficient matrix and the right-hand side vector 

100 # coeff_matrix = np.array([[A, B], [C, D]]) 

101 # rhs = np.array([E, F]) 

102 

103 # try: 

104 # center_x, center_y = np.linalg.solve(coeff_matrix, rhs) 

105 # except np.linalg.LinAlgError: 

106 # max_x = p[:, 0].max() 

107 # min_x = p[:, 0].min() 

108 # 

109 # max_y = p[:, 1].max() 

110 # min_y = p[:, 1].min() 

111 

112 # center_x = (max_x + min_x) / 2 

113 # center_y = (max_y + min_y) / 2 

114 

115 # center = np.array([center_x, center_y]) 

116 # radius = np.linalg.norm(p - center, axis=1).max() 

117 

118 # return Circle(center=center, radius=radius) 

119 

120 

121def welzl_helper(points: list[np.ndarray], r: list[np.ndarray], n: int) -> Circle: 

122 """Recursive helper function for Welzl's algorithm.""" 

123 if n == 0 or len(r) == 3: 

124 return make_circle_n_points(r) 

125 

126 # Remove a random point by shuffling it to the end 

127 # we know at this stage that n > 0 

128 idx = secrets.SystemRandom().randrange(n) 

129 p = points[idx] 

130 points[idx], points[n - 1] = points[n - 1], points[idx] 

131 

132 # Recursively compute the minimum circle without p 

133 # This is drilling down and open welzl_helper for each individual p! 

134 # R remains empty for now 

135 circle = welzl_helper(points, r, n - 1) 

136 # finally it will arrive at n == 0 and R = [] 

137 # It now calls make_circle_n_points with R = [] 

138 # It returns a circle with radius -inf 

139 # and lands back here where 

140 # n = 1 

141 # p is the final surviving point 

142 # and the circle has radius -inf 

143 # obviously p is not(!) contained in the circle! 

144 

145 # If p is inside the circle, we're done 

146 if circle.contains(p): 

147 return circle 

148 

149 # Otherwise, p must be on the boundary of the minimum enclosing circle 

150 # now we add this final point, points is still empty 

151 r.append(p) 

152 circle = welzl_helper(points, r, n - 1) 

153 r.pop() 

154 

155 return circle 

156 

157 

158def min_circle_welzl(points: np.ndarray | list[np.ndarray]) -> Circle: 

159 """Find the minimum enclosing circle using Welzl's algorithm. 

160 

161 This is the main entry point for Welzl's randomized algorithm. It takes a set of points 

162 and returns the smallest circle that contains all the points. 

163 

164 Args: 

165 points: Array or list of 2D points 

166 

167 Returns: 

168 Circle object containing the center and radius of the minimum enclosing circle 

169 """ 

170 if isinstance(points, np.ndarray): 

171 points = list(points) 

172 

173 # Make a copy of points to avoid modifying the input 

174 points = points.copy() 

175 # Shuffle the points randomly 

176 secrets.SystemRandom().shuffle(points) 

177 

178 return welzl_helper(points, [], len(points))