""" Circle packing algorithm(s) """ import numpy as np import planar as pl def fill(contour, rmin=0., rmax=1.e99, min_sep=0., nmaxpf=None, tolerance=1e4, tol_ta=1, distribution="uniform", distr_params=None): # In case it's not an array: try: _ = contour.size except: contour = np.array(contour) # Polygon pts = [tuple(item) for item in contour.tolist()] polygon = pl.Polygon(pts) center = (contour[:-1, 0].mean(), contour[:-1, 1].mean()) diameter = 2 * np.hypot(contour[:-1, 0] - center[0], contour[:-1, 1] - center[1]).max() radius = 0.5 * diameter # List of fibers rlist = [] positions = [] tryanother = 0 nonstop = True # Cell counter cc = 0 # minsep = min_sep * 0.5 while nonstop: # Random radius if distribution == "uniform": r = np.random.uniform(rmin, rmax, 1) elif distribution == "gamma": valid = False while not valid: # Work with diameters first, then back to radius mean = 2. * distr_params["mean"] shape = distr_params["shape"] scale = mean / shape r = 0.5 * np.random.gamma(shape, scale, 1) valid = (r >= rmin) & (r <= rmax) allowed_r = r + min_sep # Go trying different positions k = 0 while (k < tolerance): # Random position inside the polygon correct = False while not correct: # Try a random position dx = np.random.uniform(-radius, -radius + diameter) dy = np.random.uniform(-radius, -radius + diameter) x = center[0] + dx y = center[1] + dy # 1: Check if the center is inside the polygon correct = polygon.contains_point((x, y)) # 2: Check that the whole fiber is inside the polygon # For this, I actually check that all the points of the # polygon are farther than the fiber radius cx, cy = contour[:, 0], contour[:, 1] correct = correct & (np.hypot(cx - x, cy - y) > allowed_r).all() k += 1 # Check for clashes with other existing sub-units clash = False for i, xy in enumerate(positions): xc, yc = xy if (np.hypot(x - xc, y - yc) < allowed_r + rlist[i]): # Clash. This sub-unit doesn't fit here clash = True break if clash: k += 1 # If I tried too many times with a sub-unit and it # didn't fit, stop if k >= tolerance: tryanother += 1 if tryanother >= tol_ta: nonstop = False # Else, keep trying else: # No clash at all. Good position xy = (x, y) positions.append(xy) rlist.append(r) cc += 1 k = tolerance if cc >= nmaxpf: nonstop = False positions, rlist = np.array(positions), np.array(rlist) r = rlist[..., 0] x, y = positions.T return x, y, r