import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from sklearn.datasets import make_blobs
def expand_points_until_touching(points, max_iterations=5000, step_size=0.001):
"""
Expands circles from a set of points until they touch, using an iterative
approach.
Args:
points (np.ndarray): An array of shape (N, 2) containing the (x, y)
coordinates of the points.
max_iterations (int): The maximum number of steps to run the simulation.
step_size (float): The amount to increase the radius in each iteration.
Returns:
np.ndarray: An array of final radii for each point.
"""
num_points = len(points)
radii = np.zeros(num_points)
# A boolean array to track which points can still grow
can_grow = np.ones(num_points, dtype=bool)
for iteration in range(max_iterations):
# Stop the simulation if no more points can grow
if not np.any(can_grow):
print(f"Simulation complete after {iteration} iterations.")
break
# Expand all circles that can still grow
radii[can_grow] += step_size
# Check for collisions between all pairs of circles
for i in range(num_points):
if can_grow[i]:
for j in range(i + 1, num_points):
# Calculate the distance between the centers
distance = np.linalg.norm(points[i] - points[j])
# Check if the sum of radii is greater than or equal to the distance
if radii[i] + radii[j] >= distance:
# Collision detected!
# Calculate the new, precise radii
total_radius = distance
# Distribute the total radius proportionally
# This ensures the circles touch exactly at the collision point
# We use the previous radii as a basis for the ratio
if radii[i] + radii[j] > 0:
ratio = radii[i] / (radii[i] + radii[j])
radii[i] = total_radius * ratio
radii[j] = total_radius * (1 - ratio)
else:
# Handle the case where both radii are 0 at first contact
radii[i] = distance / 2
radii[j] = distance / 2
# Stop both circles from growing further
can_grow[i] = False
can_grow[j] = False
else:
# This 'else' block runs if the inner loop completes without a 'break'
continue
return radii
# Generate some sample data points using make_blobs
# This creates a visually interesting cluster of points
# X, _ = make_blobs(n_samples=50, centers=5, cluster_std=0.8, random_state=123)
# Run the simulation to get the final radii
final_radii = expand_points_until_touching(X)
# Create the plot
fig, ax = plt.subplots(figsize=(10, 10))
# Set plot limits to accommodate the circles
x_min, y_min = X.min(axis=0) - final_radii.max()
x_max, y_max = X.max(axis=0) + final_radii.max()
ax.set_xlim(x_min, x_max)
ax.set_ylim(y_min, y_max)
ax.set_aspect('equal', adjustable='box')
# Plot the original points as small dots
ax.scatter(X[:, 0], X[:, 1], c='gray', s=10, zorder=2)
# Plot the expanded circles
for i, (point, radius) in enumerate(zip(X, final_radii)):
circle = Circle(point, radius, color=f'C{i%10}', alpha=0.6, zorder=1)
ax.add_patch(circle)
ax.set_title("Circle Packing Simulation")
plt.show()