找到重叠圆圈的新位置

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]

---

现在我不确定我必须做什么，我的第一个想法是绘制 list2 的圆圈，从列表一中获取 x 和 y 并执行类似 和 的操作x+c，y+c其中c是一个固定值。然后我可以调用重叠函数，如果存在重叠，我可以增加该c值。这样我就有了2个for循环。现在，我的问题是：

1. 有办法避免for循环吗？
2. 是否有一个聪明的解决方案来为 list1 中的每个圆找到邻居（list2 中的圆）（不与 list2 中的其他圆重叠）？

Answer 1

使用 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

c1从每个中心到每个中心的距离数组c2

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

为相交的圆创建一个二维布尔数组。

intersect = d <= min_d

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

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

使用intersectora和b来索引 c1、c2 和 d，您应该能够获得相交圆组，然后找出如何移动 c2 中心 - 但我将把它留给另一个问题/答案。如果一个list2圆与一个list1圆相交 - 找到两个圆之间的线并沿着该线移动。如果一个list2圆与多个圆相交-沿着垂直于该圆的线list1找到两个closestlist1 litst2` 圆之间的线。circles and move the您没有提到移动圆圈的任何限制，因此可能会随机移动，然后再次找到相交，但这可能会出现问题。在下图中，弄清楚如何移动大部分红色圆圈可能很简单，但蓝色圆圈中的组可能需要不同的策略。

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

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

Answer 2

这个问题很好地可以看作是一个优化问题。更准确地说，是一个带有约束的非线性优化问题。\n由于优化策略并不总是那么容易理解，因此我将尽可能简单地定义问题，并选择一种尽可能通用（但效率较低）的方法，并且不涉及很多数学。剧透一下：我们将使用 scipy 库用不到 10 行代码来表述问题和最小化过程。

\n

然而，我仍然会提供一些提示，告诉你在哪里可以更脏。

\n

制定问题

\n

作为 NLP 类问题（非线性编程）的制定指南，您可以直接转到原始帖子中的两个要求。

\n

\n
1. 每对圆必须尽可能靠近 -> 成本函数提示
\n
2. 圆圈不得与其他（移动的）圆圈重叠 -> 约束提示
\n

\n

成本函数

\n

让我们从要最小化的成本函数的公式开始。\n由于圆应该尽可能少地移动（产生最接近的可能邻域），因此两个列表的圆之间的距离的二次惩罚项可以选择成本函数：

\n

import scipy.spatial.distance as sd\n\ndef cost_function(new_positions, old_positions):\n    new_positions = np.reshape(new_positions, (-1, 2))\n    return np.trace(sd.cdist(new_positions, old_positions, metric=\'sqeuclidean\'))\n

\n

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

\n

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

\n

让我们尝试一下：例如，让我们使用 Voronoi 方法（由 Paul Brodersen）计算的位置。

\n

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

\n

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

\n

非线性约束

\n

我们知道圆圈不能与新集合中的其他圆圈重叠。如果我们将其转化为约束，我们会发现距离的下限是半径的 2 倍，而上限就是无穷大。

\n

import scipy.spatial.distance as sd\nfrom scipy.optimize import NonlinearConstraint\n\ndef cons_f(x):\n    x = np.reshape(x, (-1, 2))\n    return sd.pdist(x)\n\nnonlinear_constraint = NonlinearConstraint(cons_f, 2*radius, np.inf, jac=\'2-point\')\n

\n

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

\n

另请注意：您还可以为位置本身设置边界约束，使其不超过指定区域。然后可以将其用作另一个参数。（看scipy.optimize.Bounds）

\n

最小化约束下的成本函数

\n

现在我们已经有了成本函数和约束这两个要素。现在让我们最小化整个事情！

\n

from scipy.optimize import minimize\n\nres = minimize(lambda x: cost_function(x, positions), positions.flatten(), method=\'SLSQP\',\n               jac="2-point", constraints=[nonlinear_constraint])\n

\n

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

\n

另请注意，我们必须始终将参数（指定 n 个圆的新布局的位置的 nx2 向量）作为平面向量传递。因此，在代码中可以多次发现重塑。

\n

评估、总结和可视化

\n

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

\n

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

\n

从 Voronoi 方法开始，我们实际上又降低了 87% 的成本！借助现代优化策略的力量，我们可以立即解决很多问题。

\n

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

\n

表现：77.1 ms \xc2\xb1 1.17 ms

\n

整个代码：

\n

from scipy.optimize import minimize\nimport scipy.spatial.distance as sd\nfrom scipy.optimize import NonlinearConstraint\n\n# Given by original post\npositions = np.array([x, y]).T\n\ndef cost_function(new_positions, old_positions):\n    new_positions = np.reshape(new_positions, (-1, 2))\n    return np.trace(sd.cdist(new_positions, old_positions, metric=\'sqeuclidean\'))\n\ndef cons_f(x):\n    x = np.reshape(x, (-1, 2))\n    return sd.pdist(x)\n\nnonlinear_constraint = NonlinearConstraint(cons_f, 2*radius, np.inf, jac=\'2-point\')\nres = minimize(lambda x: cost_function(x, positions), positions.flatten(), method=\'SLSQP\',\n               jac="2-point", constraints=[nonlinear_constraint])\n

\n