寻找重叠圆的新位置

Question

我正在尝试编写一段代码，针对给定的圆列表（list1），能够找到新圆的位置（list2）。list1和list2具有相同的长度，因为对于list1中的每个圆，都必须有一个来自list2的圆。

1. 每对圆（假设是list1中的circle1和list2中的circle2）必须尽可能靠近，
2. list2中的圆不能与list1中的圆重叠，而单个列表中的圆可以相互重叠。

list1是固定的，因此现在我必须找到list2中圆的正确位置。

我编写了这个简单的函数来识别两个圆是否重叠：

def overlap(x1, y1, x2, y2, r1, r2):
    distSq = (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)
    radSumSq = (r1 + r2) * (r1 + r2)
    if (distSq >= radSumSq):
        return False # no overlap
    else:
        return True  #overlap

这是列表1：

包含：

x=[14.11450195 14.14184093 14.15435028 14.16206741 14.16951752 14.17171097
 14.18569565 14.19700241 14.23129082 14.24083233 14.24290752 14.24968338
 14.2518959  14.26536751 14.27209759 14.27612877 14.2904377  14.29187012
 14.29409599 14.29618549 14.30615044 14.31624985 14.3206892  14.3228569
 14.36143875 14.36351967 14.36470699 14.36697292 14.37235737 14.41422081
 14.42583466 14.43226814 14.43319225 14.4437027  14.4557848  14.46592999
 14.47036076 14.47452068 14.47815609 14.52229309 14.53059006 14.53404236
 14.5411644 ] 
y=[-0.35319126 -0.44222349 -0.44763246 -0.35669261 -0.24366629 -0.3998799
 -0.38940558 -0.57744932 -0.45223859 -0.21021004 -0.44250247 -0.45866323
 -0.47203487 -0.51684451 -0.44884869 -0.2018993  -0.40296811 -0.23641759
 -0.18019417 -0.33391538 -0.53565156 -0.45215255 -0.40939832 -0.26936951
 -0.30894437 -0.55504167 -0.47177047 -0.45573688 -0.43100587 -0.5805912
 -0.21770373 -0.199422   -0.17372169 -0.38522363 -0.56950212 -0.56947368
 -0.48770753 -0.24940367 -0.31492445 -0.54263926 -0.53460872 -0.4053807
 -0.43733299]
radius = 0.014

---

可拷贝、可粘贴...

x = [14.11450195,14.14184093,14.15435028,14.16206741,14.16951752,
     14.17171097,14.18569565,14.19700241,14.23129082,14.24083233,
     14.24290752,14.24968338,14.2518959,14.26536751,14.27209759,
     14.27612877,14.2904377,14.29187012,14.29409599,14.29618549,
     14.30615044,14.31624985,14.3206892,14.3228569,14.36143875,
     14.36351967,14.36470699,14.36697292,14.37235737,14.41422081,
     14.42583466,14.43226814,14.43319225,14.4437027,14.4557848,
     14.46592999,14.47036076,14.47452068,14.47815609,14.52229309,
     14.53059006,14.53404236,14.5411644]

y = [-0.35319126,-0.44222349,-0.44763246,-0.35669261,-0.24366629,
     -0.3998799,-0.38940558,-0.57744932,-0.45223859,-0.21021004,
     -0.44250247,-0.45866323,-0.47203487,-0.51684451,-0.44884869,
     -0.2018993,-0.40296811,-0.23641759,-0.18019417,-0.33391538,
     -0.53565156,-0.45215255,-0.40939832,-0.26936951,-0.30894437,
     -0.55504167,-0.47177047,-0.45573688,-0.43100587,-0.5805912,
     -0.21770373,-0.199422,-0.17372169,-0.38522363,-0.56950212,
     -0.56947368,-0.48770753,-0.24940367,-0.31492445,-0.54263926,
     -0.53460872,-0.4053807,-0.43733299]

---

我现在不确定该做什么，我的第一个想法是以列表1的x和y为基础，画出列表2的圆，并进行像 x+c 和 y+c 这样的操作，其中c是一个固定值。然后我可以调用我的重叠函数，如果有重叠，我就可以增加c 的值。 这样我就有了两个 for 循环。现在，我的问题是：

1. 有没有办法避免使用 for 循环？
2. 是否有一种聪明的方法来为列表1中的每个圆找到邻居(列表2中的圆)，而不与列表2中的其他圆重叠？

- Alessandro Peca

---

有没有避免使用for循环的方法？你无法避免迭代，但如果你使用numpy ndarrays，你可以将迭代推入C代码中，这样会更快。 - wwii

你的实际数据集有多大？这两个列表中有多少个圆？一次只处理一个方面，使用小数据集进行优化，然后在一切正常时再进行优化。看起来你可能需要保存每个重叠圆的最近邻信息，然后使用组的几何形状计算一个“安全”向量来移动有问题的圆，这样你就不必迭代随机向量并再次迭代了。你可能需要对这些组使用创意容器。 - wwii

你能提供一个 list2 的例子吗？ - wwii

在我的主题中有这些答案。您已经有了一个带有x和y坐标的list1示例。样本list2类似于list1，具有相同数量的圆和相同的半径。 真实数据集更大，因此我想优化该过程，避免使用for循环。 - Alessandro Peca

但是如果没有对list2的要求进行进一步的说明，就无法回答。 - norok2

Answer 1

这个问题可以看作是一个优化问题，确切地说，是一个具有约束的非线性优化问题。由于优化策略并不总是易于理解，我将尽可能简单地定义问题，并选择一种尽可能通用（但效率较低）且不涉及大量数学的方法。作为剧透：我们将使用scipy库在不到10行代码中阐述问题和最小化过程。

然而，我仍会提供一些提示，让你更深入了解这个问题。

制定问题

作为NLP类问题（非线性规划）制定方案的指南，您可以直接前往原始帖子中的两个要求。

1.每对圆之间必须尽可能接近-＆gt;成本函数的提示 2.圆不能与其他（移动的）圆重叠-＆gt;约束条件的提示

代价函数

让我们从最小化的代价函数制定开始。由于应尽可能少地移动圆（导致尽可能接近的邻域），因此可以选择两个列表中圆之间距离的二次惩罚项作为成本函数：

import scipy.spatial.distance as sd

def cost_function(new_positions, old_positions):
    new_positions = np.reshape(new_positions, (-1, 2))
    return np.trace(sd.cdist(new_positions, old_positions, metric='sqeuclidean'))

为什么选择二次方程？这部分是因为可微性以及随机性原因（将圆视为正态分布的测量误差，最小二乘法则成为极大似然估计）。通过利用代价函数的结构，可以提高优化的效率（消除平方根）。顺便说一句，这个问题与非线性回归有关，其中也使用（非线性）最小二乘法。

现在我们手头有一个代价函数，也有一个评价优化效果的好方法。为了能够比较不同优化策略的解决方案，我们只需将新计算出的位置传递给代价函数即可。

让我们试试吧：例如，让我们使用由Paul Brodersen提出的Voronoi方法计算出的位置。

print(cost_function(new_positions, old_positions))
# prints 0.007999244511697411

如果你问我，这是一个相当不错的价值。考虑到当没有任何位移时，成本函数输出为零，因此这个成本相当接近。现在我们可以尝试通过使用经典优化来超越这个价值！

非线性约束

我们知道，在新集合中，圆不能与其他圆重叠。如果将此转换为约束条件，我们发现距离的下限是半径的2倍，上限则是无穷大。

import scipy.spatial.distance as sd
from scipy.optimize import NonlinearConstraint

def cons_f(x):
    x = np.reshape(x, (-1, 2))
    return sd.pdist(x)

nonlinear_constraint = NonlinearConstraint(cons_f, 2*radius, np.inf, jac='2-point')

在这里，我们通过有限差分（见参数jac='2-point'）来近似雅可比矩阵，使生活变得简单。此时应该说，我们可以通过自己编写一阶和二阶导数的公式而不是使用逼近法来提高效率。但是这留给有兴趣的读者。（因为我们在这里使用的距离计算非常简单，所以并不难。）

另外一个注意事项：您还可以为位置本身设置边界约束，以不超过指定区域。这也可以作为另一个参数使用。（请参见scipy.optimize.Bounds）

现在我们已经准备好了两个要素：代价函数和约束条件。现在让我们最小化整个东西！

from scipy.optimize import minimize

res = minimize(lambda x: cost_function(x, positions), positions.flatten(), method='SLSQP',
               jac="2-point", constraints=[nonlinear_constraint])

正如您所看到的，我们在这里也近似求出了一阶导数。您还可以深入其中，并自行设置导数（分析性地）。

还要注意的是，我们必须始终将参数（指定n个圆的新布局位置的nx2向量）作为平坦向量传递。因此，在代码中可以多次找到重塑操作。

评估、总结和可视化

让我们看看优化结果在成本函数中的表现：

new_positions = np.reshape(res.x, (-1,2))
print(cost_function(new_positions, old_positions))
# prints 0.0010314079483565686

从 Voronoi 方法开始，我们实际上将成本再次降低了 87%！多亏了现代优化策略的强大能力，我们可以在瞬间解决很多问题。

当然，看看移动后的圆形现在是什么样子也是很有趣的： 优化后的圆形

性能：77.1 毫秒 ± 1.17 毫秒

整个代码：

from scipy.optimize import minimize
import scipy.spatial.distance as sd
from scipy.optimize import NonlinearConstraint

# Given by original post
positions = np.array([x, y]).T

def cost_function(new_positions, old_positions):
    new_positions = np.reshape(new_positions, (-1, 2))
    return np.trace(sd.cdist(new_positions, old_positions, metric='sqeuclidean'))

def cons_f(x):
    x = np.reshape(x, (-1, 2))
    return sd.pdist(x)

nonlinear_constraint = NonlinearConstraint(cons_f, 2*radius, np.inf, jac='2-point')
res = minimize(lambda x: cost_function(x, positions), positions.flatten(), method='SLSQP',
               jac="2-point", constraints=[nonlinear_constraint])

- profmori4rty

---

不错！提一下 positions = np.array([x, y]).T 会更好。我想这里的缺点是需要更长的计算时间。 - kilojoules

谢谢提示！我已经添加了。确实，在我的机器上，优化大约需要40秒。 - profmori4rty

顺便说一下，你的方法也非常好。它得出的成本与我的解决方案非常相似：0.0011436334250168178。 - profmori4rty

@kilojoules 我进行了更多的实验并将方法改为SLSQP。这消除了对二阶导数的近似->速度提升570倍（约77毫秒计算时间）。成本保持不变，因为它是最优解。 - profmori4rty

Answer 2

使用NumPy数组，您可以避免使用for循环。

从您的示例设置开始。

import numpy as np

#Using your x and y
c1 = np.array([x,y]).T

# random set of other centers within the same range as c1
c2 = np.random.random((10,2))
np.multiply(c2, c1.max(0)-c1.min(0),out = c2)
np.add(c2, c1.min(0), out=c2)

radius = 0.014
r = radius
min_d = (2*r)*(2*r)

plot_circles(c1,c2)    # see function at end

< p > 从 c1 中每个中心到 c2 中每个中心的距离数组。 < /p>

def dist(c1,c2):
    dx = c1[:,0,None] - c2[:,0]
    dy = c1[:,1,None] - c2[:,1]
    return dx*dx + dy*dy

d = dist(c1,c2)

或者您可以使用scipy.spatial

from scipy.spatial import distance
d = distance.cdist(c1,c2,'sqeuclidean')

创建一个2D的布尔数组，用于表示相交的圆。

intersect = d <= min_d

从两个集合中查找重叠圆的索引。

a,b = np.where(intersect)
plot_circles(c1[a],c2[b])

使用`intersect`或者`a`和`b`索引`c1`、`c2`和`d`，你应该能够获得一组相交的圆，然后找出如何移动`c2`的中心——但我会把这个问题留给另一个问题/答案。如果一个`list2`圆与一个`list1`圆相交，找到两者之间的线并沿着该线移动。如果一个`list2`圆与多个`list1`圆相交，找到两个最靠近的`list1`圆之间的线，并沿垂直于该线的方向移动`list2`圆。你没有提及移动圆的任何限制，所以也许可以随机移动，然后再次找到相交点，但这可能会出现问题。在下面的图像中，大多数红色圆的移动可能很简单，但是蓝色标记的一组可能需要不同的策略。

以下是一些获取组的示例：

>>> for f,g,h in zip(c1[a],c2[b],d[a,b]):
    print(f,g,h)
>>> c1[intersect.any(1)],c2[intersect.any(0)]
>>> for (f,g) in zip(c2,intersect.T):
    if g.any():
            print(f.tolist(),c1[g].tolist())

---

import matplotlib as mpl
from matplotlib import pyplot as plt

def plot_circles(c1,c2):

    bounds = np.array([c1.min(0),c2.min(0),c1.max(0),c2.max(0)])
    xmin, ymin = bounds.min(0)
    xmax, ymax = bounds.max(0)

    circles1 = [mpl.patches.Circle(xy,radius=r,fill=False,edgecolor='g') for xy in c1]
    circles2 = [mpl.patches.Circle(xy,radius=r,fill=False,edgecolor='r') for xy in c2]
    fig = plt.figure()
    ax = fig.add_subplot(111)
    for c in circles2:
        ax.add_artist(c)
    for c in circles1:
        ax.add_artist(c)

    ax.set_xlim(xmin-r,xmax+r)
    ax.set_ylim(ymin-r,ymax+r)

    plt.show()
    plt.close()

- wwii

---

我开始了一个悬赏，希望这算作第二个问题，以找到如何移动c2中心。 - kilojoules

Answer 3

一种解决方法是沿着每个圆之间不想要的间距的梯度进行跟踪，但也许还有更好的方法。这种方法有几个参数需要调整，并且运行需要一些时间。

import matplotlib.pyplot as plt
from scipy.optimize import minimize as mini
import numpy as np
from scipy.optimize import approx_fprime

x = np.array([14.11450195,14.14184093,14.15435028,14.16206741,14.16951752,
     14.17171097,14.18569565,14.19700241,14.23129082,14.24083233,
     14.24290752,14.24968338,14.2518959,14.26536751,14.27209759,
     14.27612877,14.2904377,14.29187012,14.29409599,14.29618549,
     14.30615044,14.31624985,14.3206892,14.3228569,14.36143875,
     14.36351967,14.36470699,14.36697292,14.37235737,14.41422081,
     14.42583466,14.43226814,14.43319225,14.4437027,14.4557848,
     14.46592999,14.47036076,14.47452068,14.47815609,14.52229309,
     14.53059006,14.53404236,14.5411644])

y = np.array([-0.35319126,-0.44222349,-0.44763246,-0.35669261,-0.24366629,
     -0.3998799,-0.38940558,-0.57744932,-0.45223859,-0.21021004,
     -0.44250247,-0.45866323,-0.47203487,-0.51684451,-0.44884869,
     -0.2018993,-0.40296811,-0.23641759,-0.18019417,-0.33391538,
     -0.53565156,-0.45215255,-0.40939832,-0.26936951,-0.30894437,
     -0.55504167,-0.47177047,-0.45573688,-0.43100587,-0.5805912,
     -0.21770373,-0.199422,-0.17372169,-0.38522363,-0.56950212,
     -0.56947368,-0.48770753,-0.24940367,-0.31492445,-0.54263926,
     -0.53460872,-0.4053807,-0.43733299])
radius = 0.014
x0, y0 = (x, y)

def plot_circles(x, y, name='initial'):
   fig, ax = plt.subplots()
   for ii in range(x.size):
      ax.add_patch(plt.Circle((x[ii], y[ii]), radius, color='b', fill=False))
   ax.set_xlim(x.min() - radius, x.max() + radius)
   ax.set_ylim(y.min() - radius, y.max() + radius)
   fig.savefig(name)
   plt.clf()

def spacing(s):
   x, y = np.split(s, 2)
   dX, dY = [np.subtract(*np.meshgrid(xy, xy, indexing='ij')).T
      for xy in [x, y]]
   dXY2 = dX**2 + dY**2
   return np.minimum(dXY2[np.triu_indices(x.size, 1)] - (2 * radius) ** 2, 0).sum() 

plot_circles(x, y)

def spacingJ(s):
   return approx_fprime(s, spacing, 1e-8)

s = np.append(x, y)
for ii in range(50):
   j = spacingJ(s)
   if j.sum() == 0: break
   s += .01 * j 

   x_new, y_new = np.split(s, 2)
   plot_circles(x_new, y_new, 'new%i' % ii)

plot_circles(x_new, y_new, 'new%i' % ii)

https://giphy.com/gifs/x0lWDLZBz5O3gWTbLa

- kilojoules

---

最直接的方法。然而，我认为这些解决方案并不是最优的：例如，假设有三个圆，它们的中心共线，外部两个圆的环在中间环的中心处相接。最优解是将中间环垂直地向左或向右移动到通过外部两个环的中心的线上。您的解决方案是沿着该线移动外部两个环，直到它们不再重叠中心环，导致与最优解相比位移加倍。@wwii已经在他的答案末尾指出了这个问题。 - Paul Brodersen

Answer 4

本答案实现了Lloyd算法的一个变种。基本思想是计算您的点/圆的Voronoi图。这将为每个点分配一个单元格，该单元格是包括该点并具有距离所有其他点最远的中心的区域。

在原始算法中，我们会将每个点移动到其Voronoi单元格的中心。随着时间的推移，这将导致点均匀分布，如这里所示。

在这种变体中，我们只移动重叠另一个点的点。

import numpy as np
import matplotlib.pyplot as plt

from scipy.spatial import Voronoi
from scipy.spatial.distance import cdist


def remove_overlaps(positions, radii, tolerance=1e-6):
    """Use a variation of Lloyds algorithm to move circles apart from each other until none overlap.

    Parameters
    ----------
    positions : array
        The (x, y) coordinates of the circle origins.
    radii : array
        The radii for each circle.
    tolerance : float
        If all circles overlap less than this threshold, the computation stops.
        Higher values leads to faster convergence.

    Returns
    -------
    new_positions : array
        The (x, y) coordinates of the circle origins.

    See also
    --------
    https://en.wikipedia.org/wiki/Lloyd%27s_algorithm

    """

    positions = np.array(positions)
    radii = np.array(radii)

    minimum_distances = radii[np.newaxis, :] + radii[:, np.newaxis]
    minimum_distances[np.diag_indices_from(minimum_distances)] = 0 # ignore distances to self

    # Initialize the first loop.
    distances = cdist(positions, positions)
    displacements = np.max(np.clip(minimum_distances - distances, 0, None), axis=-1)

    while np.any(displacements > tolerance):
        centroids = _get_voronoi_centroids(positions)

        # Compute the direction from each point towards its corresponding Voronoi centroid.
        deltas = centroids - positions
        magnitudes = np.linalg.norm(deltas, axis=-1)
        directions = deltas / magnitudes[:, np.newaxis]

        # Mask NaNs that arise if the magnitude is zero, i.e. the point is already center of the Voronoi cell.
        directions[np.isnan(directions)] = 0

        # Step into the direction of the centroid.
        # Clipping prevents overshooting of the centroid when stepping into the direction of the centroid.
        # We step by half the displacement as the other overlapping point will be moved in approximately the opposite direction.
        positions = positions + np.clip(0.5 * displacements, None, magnitudes)[:, np.newaxis] * directions

        # Initialize next loop.
        distances = cdist(positions, positions)
        displacements = np.max(np.clip(minimum_distances - distances, 0, None), axis=-1)

    return positions


def _get_voronoi_centroids(positions):
    """Construct a Voronoi diagram from the given positions and determine the center of each cell."""
    voronoi = Voronoi(positions)
    centroids = np.zeros_like(positions)
    for ii, idx in enumerate(voronoi.point_region):
        region = [jj for jj in voronoi.regions[idx] if jj != -1] # i.e. ignore points at infinity; TODO: compute correctly clipped regions
        centroids[ii] = np.mean(voronoi.vertices[region], axis=0)
    return centroids


if __name__ == '__main__':

    x = np.array([14.11450195,14.14184093,14.15435028,14.16206741,14.16951752,
                  14.17171097,14.18569565,14.19700241,14.23129082,14.24083233,
                  14.24290752,14.24968338,14.2518959,14.26536751,14.27209759,
                  14.27612877,14.2904377,14.29187012,14.29409599,14.29618549,
                  14.30615044,14.31624985,14.3206892,14.3228569,14.36143875,
                  14.36351967,14.36470699,14.36697292,14.37235737,14.41422081,
                  14.42583466,14.43226814,14.43319225,14.4437027,14.4557848,
                  14.46592999,14.47036076,14.47452068,14.47815609,14.52229309,
                  14.53059006,14.53404236,14.5411644])
    y = np.array([-0.35319126,-0.44222349,-0.44763246,-0.35669261,-0.24366629,
                  -0.3998799,-0.38940558,-0.57744932,-0.45223859,-0.21021004,
                  -0.44250247,-0.45866323,-0.47203487,-0.51684451,-0.44884869,
                  -0.2018993,-0.40296811,-0.23641759,-0.18019417,-0.33391538,
                  -0.53565156,-0.45215255,-0.40939832,-0.26936951,-0.30894437,
                  -0.55504167,-0.47177047,-0.45573688,-0.43100587,-0.5805912,
                  -0.21770373,-0.199422,-0.17372169,-0.38522363,-0.56950212,
                  -0.56947368,-0.48770753,-0.24940367,-0.31492445,-0.54263926,
                  -0.53460872,-0.4053807,-0.43733299])
    radius = 0.014

    positions = np.c_[x, y]
    radii = np.full(len(positions), radius)

    fig, axes = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(14, 7))
    for position, radius in zip(positions, radii):
        axes[0].add_patch(plt.Circle(position, radius, fill=False))
    axes[0].set_xlim(x.min() - radius, x.max() + radius)
    axes[0].set_ylim(y.min() - radius, y.max() + radius)
    axes[0].set_aspect('equal')

    new_positions = remove_overlaps(positions, radii)
    for position, radius in zip(new_positions, radii):
        axes[1].add_patch(plt.Circle(position, radius, fill=False))

    for ax in axes.ravel():
        ax.set_aspect('equal')

    plt.show()

- Paul Brodersen

---

太棒了！您认为这个功能可以扩展到强制圆圈保持在一个框内吗？ - kilojoules

@kilojoules 当然，这很简单。看看我在这里的Lloyds算法的另一个实现或者在这里的stackoverflow。 - Paul Brodersen