Circle Packing

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()