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

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 if p2[1] == p1[1]: # horizontal line 

26 return np.inf # Perpendicular bisector is vertical 

27 return -(p2[0] - p1[0]) / (p2[1] - p1[1]) 

28 

29 

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

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

32 

33 Args: 

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

35 

36 Returns: 

37 Circle object containing the center and radius of the circle 

38 """ 

39 num_points = len(matrix) 

40 

41 if num_points == 0: 

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

43 

44 if num_points == 1: 

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

46 

47 if num_points == 2 or num_points == 3: 

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

49 

50 # if num_points == 2: 

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

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

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

54 # 

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

56 

57 # if num_points == 3: 

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

59 # # p1, p2, p3 = matrix 

60 # p = np.array(matrix) 

61 # 

62 # # Midpoints of the sides 

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

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

65 # 

66 # # Slopes of the perpendicular bisectors 

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

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

69 # 

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

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

72 # 

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

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

75 # center_x = mid12[0] 

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

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

78 # center_x = mid23[0] 

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

80 # else: 

81 # # Calculate the intersection of the two perpendicular bisectors 

82 # A, B = m1, -1 

83 # C, D = m2, -1 

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

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

86 # 

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

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

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

90 

91 # try: 

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

93 # except np.linalg.LinAlgError: 

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

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

96 # 

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

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

99 

100 # center_x = (max_x + min_x) / 2 

101 # center_y = (max_y + min_y) / 2 

102 

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

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

105 

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

107 

108 

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

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

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

112 return make_circle_n_points(r) 

113 

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

115 # we know at this stage that n > 0 

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

117 p = points[idx] 

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

119 

120 # Recursively compute the minimum circle without p 

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

122 # R remains empty for now 

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

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

125 # It now calls make_circle_n_points with R = [] 

126 # It returns a circle with radius -inf 

127 # and lands back here where 

128 # n = 1 

129 # p is the final surviving point 

130 # and the circle has radius -inf 

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

132 

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

134 if circle.contains(p): 

135 return circle 

136 

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

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

139 r.append(p) 

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

141 r.pop() 

142 

143 return circle 

144 

145 

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

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

148 

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

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

151 

152 Args: 

153 points: Array or list of 2D points 

154 

155 Returns: 

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

157 """ 

158 if isinstance(points, np.ndarray): 

159 points = list(points) 

160 

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

162 points = points.copy() 

163 # Shuffle the points randomly 

164 secrets.SystemRandom().shuffle(points) 

165 

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