生成多個互不重疊的不同半徑圓
阿新 • • 發佈:2018-12-16
生成多個互不重疊的不同半徑圓
https://www.zhihu.com/question/53012468
固定區域,生成所有圓,統計圓的面積,最後面積最大的就是答案。
import numpy as np from scipy.optimize import minimize from scipy.spatial.distance import pdist import matplotlib.pyplot as plt import os def obj_grad_round(x, radii, pair_index_to_index_pair, repulsion_scale, attraction_scale): n = len(radii) centroids = np.reshape(x, (n, 2)) circles_pdist_sqrd = pdist(centroids, 'sqeuclidean') obj = 0 grad = np.zeros((n, 2)) for pidx, dist_sqrd in enumerate(circles_pdist_sqrd): i, j = pair_index_to_index_pair[pidx] touching_dist_sqrd = (radii[i] + radii[j]) ** 2 grad_inc = 4 * (dist_sqrd - touching_dist_sqrd) * (centroids[i, :] - centroids[j, :]) scale = repulsion_scale if dist_sqrd < touching_dist_sqrd else attraction_scale obj += scale * (dist_sqrd - touching_dist_sqrd) ** 2 grad[i, :] += scale * grad_inc grad[j, :] -= scale * grad_inc return obj, grad.flatten() def obj_grad_square(x, radii, pair_index_to_index_pair, repulsion_scale, attraction_scale): n = len(radii) centroids = np.reshape(x, (n, 2)) circles_pdist_sqrd = pdist(centroids, 'sqeuclidean') obj = 0 grad = np.zeros((n, 2)) for pidx, dist_sqrd in enumerate(circles_pdist_sqrd): i, j = pair_index_to_index_pair[pidx] touching_dist_sqrd = (radii[i] + radii[j]) ** 2 diff = centroids[i, :] - centroids[j, :] if dist_sqrd < touching_dist_sqrd: # repulsion obj += repulsion_scale * (dist_sqrd - touching_dist_sqrd) ** 2 grad_inc = 4 * (dist_sqrd - touching_dist_sqrd) * diff grad[i, :] += repulsion_scale * grad_inc grad[j, :] -= repulsion_scale * grad_inc else: # attraction diff_sqrd = diff ** 2 if max(diff_sqrd) >= touching_dist_sqrd: # do not overlap after projection, needs attraction if diff_sqrd[0] >= diff_sqrd[1]: obj_inc = (diff_sqrd[0] - touching_dist_sqrd) ** 2 grad_inc = np.array([4 * (diff_sqrd[0] - touching_dist_sqrd) * diff[0], 0]) else: obj_inc = (diff_sqrd[1] - touching_dist_sqrd) ** 2 grad_inc = np.array([0, 4 * (diff_sqrd[1] - touching_dist_sqrd) * diff[1]]) obj += attraction_scale * obj_inc grad[i, :] += attraction_scale * grad_inc grad[j, :] -= attraction_scale * grad_inc return obj, grad.flatten() def find_layout(radii, shape='square', repulsion_scale=1e6, attraction_scale=1, attraction_decay=0.3, overlapping_tol=0., seed=None, verbose=False): assert shape == 'round' or shape == 'square' obj_grad = obj_grad_round if shape == 'round' else obj_grad_square n = len(radii) n_pairs = int(n * (n-1) / 2) pair_index_to_index_pair = [None] * n_pairs pidx = 0 for i in range(n-1): for j in range(i+1, n): pair_index_to_index_pair[pidx] = i, j pidx += 1 if seed: np.random.seed(seed) if verbose: options = {'disp': True} history_folder = 'history_%s' % shape if not os.path.exists(history_folder): os.mkdir(history_folder) else: options = None overlapped = True n_trials = 1 x0 = np.random.randn(2 * n) while overlapped: if verbose: print('repulsion scale = %g\tattraction scale = %g' % (repulsion_scale, attraction_scale)) opt_result = minimize(obj_grad, x0, args=(radii, pair_index_to_index_pair, repulsion_scale, attraction_scale), jac=True, options=options) centroids = np.reshape(opt_result.x, (n, 2)) if verbose: if not opt_result.success: print(opt_result.message) draw_layout(centroids, radii, '%s/%02d_attraction=%g.png' % (history_folder, n_trials, attraction_scale)) np.savetxt('%s/%02d_attraction=%g.dat' % (history_folder, n_trials, attraction_scale), centroids, fmt='%g') overlapped = check_overlapping(centroids, radii, overlapping_tol, verbose) if overlapped: attraction_scale *= attraction_decay n_trials += 1 x0 = opt_result.x else: overlapped = False return centroids def check_overlapping(centroids, radii, tol, verbose=False): n = len(radii) for i in range(n-1): for j in range(i+1, n): dij = np.linalg.norm(centroids[i, :] - centroids[j, :]) rhs = radii[i] + radii[j] - tol if dij < rhs: if verbose: print('Circle {i} and circle {j} overlap ({dij} = dist(C_{i}, C_{j}) < r_{i} + r_{j} - tol = {rhs}).'.format(i=i, j=j, dij=dij, rhs=rhs)) return True return False def draw_layout(centroids, radii, not_show_and_save_to=None): fig, ax = plt.subplots() for centroid, radius in zip(centroids, radii): circle = plt.Circle(centroid, radius) ax.add_artist(circle) x = centroids[:, 0] y = centroids[:, 1] xmin = min(x - radii) xmax = max(x + radii) ymin = min(y - radii) ymax = max(y + radii) plt.axis([xmin, xmax, ymin, ymax]) ax.set_aspect('equal') if not_show_and_save_to: plt.savefig(not_show_and_save_to) else: plt.show() plt.close() if __name__ == '__main__': seed = 1 np.random.seed(seed) n_circles = 100 radii = np.random.rand(n_circles) / 10 + 0.1 # [0.1, 0.2) overlapping_tol = 1e-4 shape = 'square' # 'round' or 'square' verbose = False layout = find_layout(radii, shape, overlapping_tol=overlapping_tol, verbose=verbose) draw_layout(layout, radii, 'final_%s.png' % shape) np.savetxt('radii.dat', radii, fmt='%g') np.savetxt('final_%s.dat' % shape, layout, fmt='%g')