版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/jacke121/article/details/85016402
生成多个互不重叠的不同半径圆
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')