在 3D 空间中生成具有最小最近邻距离的随机点

Question

类似问题：
Generating N random points with certain predefined distance between them

choose n most distant points in R

但它们要么在 matlab 中，要么没有完成所需的任务。

我必须在一个长度为的盒子内创建 N 个点 任意两点之间的距离大于delta。

例如： 假设我在 x、y、z 轴上有一个长度为 10 埃的盒子。
我想在这个盒子里有 200 个随机点，这样最小距离 任意两点之间的距离大于 3 埃。

尝试：

#!python
# -*- coding: utf-8 -*-#
import numpy as np
np.random.seed(100)
np.set_printoptions(2)

box_length = 10
d = box_length
threshold = 6
num_points = 5
x1, y1, z1 = np.random.random(num_points)* box_length, np.random.random(num_points)* box_length, np.random.random(num_points)* box_length
x2, y2, z2 = np.random.random(num_points)* box_length, np.random.random(num_points)* box_length, np.random.random(num_points)* box_length

# print(len(x1))

# just for checking make ponts integers
for i in range(len(x1)):
    x1[i] = int(x1[i])
    x2[i] = int(x2[i])
    y1[i] = int(y1[i])
    y2[i] = int(y2[i])
    z1[i] = int(z1[i])
    z2[i] = int(z2[i])


print(x1)
print(y1)
print(z1)
print("\n")

pt1_lst = []
pt2_lst = []
for i in range(len(x1)):
    a, b, c    = x1[i], y1[i], z1[i]
    a2, b2, c2 = x2[i], y2[i], z2[i]
    dist2      = (a-a2)**2 + (b-b2)**2 + (c-c2)**2

    print("\n")
    print(a,b,c)
    print(a2,b2,c2)
    print(dist2)

    if dist2 > threshold**2:
        pt1 = (a,b,c)
        pt2 = (a2,b2,c2)
        pt1_lst.append(pt1)
        pt2_lst.append(pt2)



print("points")
print(pt1_lst)
print(pt2_lst)

代码问题： 此代码比较 points1 和 points2 的点，但不会在 points1 和 points2 内部进行比较。

可能有更好的算法来解决这个问题，并向那些人致敬 提出解决问题的好主意的人。

谢谢。

附言： 我做了一些研究并尝试找到相关链接，但无法解决 问题。 还有一些相关链接是：

更新::

我在下面尝试了 Stefans 的代码，它适用于 N = 10，但我尝试了 N = 200 并且它使用了非常长的时间（我在 10 分钟后停止了代码）。

有什么有效的方法吗？

我们将不胜感激！

All paths of length L from node n using python
Create Random Points Inside Defined Rectangle with Python

Answer 1

假设我在 x、y、z 轴上有一个长度为 10 埃的盒子。 我想在这个盒子里有 10 个随机点，这样任何两点之间的最小距离都大于 3 埃。

我认为这是可行的，在该框中重复生成十个随机点，直到距离都足够大：

>>> import numpy as np
>>> from itertools import combinations

>>> while True:
        P = np.random.rand(10, 3) * 10
        if all(np.linalg.norm(p - q) > 3
               for p, q in combinations(P, 2)):
            break

>>> P
array([[ 9.02322366,  6.13576854,  3.1745708 ],
       [ 6.48005836,  7.5280536 ,  4.66442095],
       [ 5.78306167,  1.83922896,  9.48337683],
       [ 0.70507032,  0.20737532,  5.31191608],
       [ 3.71977864,  6.40278939,  3.81742814],
       [ 0.03938102,  6.7705456 ,  6.28841217],
       [ 3.27845597,  2.98811665,  4.81792286],
       [ 7.74422021,  9.30027671,  8.1770998 ],
       [ 0.28544716,  0.35155801,  9.77847352],
       [ 4.84536373,  4.21378476,  0.4456017 ]])

大约需要 50 次尝试才能找到一组好的点。这里我试了1000次，20次就好了：

>>> sum(all(np.linalg.norm(p - q) > 3
            for p, q in combinations(np.random.rand(10, 3) * 10, 2))
        for _ in range(1000))
20

Answer 2

这种分布的通用名称是泊松球采样。有已知的 O(n) 可以做到这一点 - 请检查 here

Answer 3

拉丁超立方体将域划分为单元格，随机选择一组单元格，并在每个单元格中放置一个随机点。在一些不幸的情况下，一些点可能会彼此太近，但我们可以检查距离并可能生成一个新的案例。

这是一个基于OpenTURNS 的二维域的实现。检查所有点的距离，但在更复杂的示例中，kd-tree 可能更有效。绘图代码取自here。

"""
Generate random points in a rectangle, with a minimum distance.
"""

# %% Import.

import openturns as ot
import numpy as np
from scipy.spatial.distance import pdist
import matplotlib.pyplot as plt
import matplotlib.animation as animation


# %% Generate points.

# Boundaries.
x_min, x_max = -15, +15
y_min, y_max = -10, +10

# Number of points per case.
n_points = 8

# Number of cases to generate.
n_cases = 50

# Minimum distance.
min_dist = 4

# Uniform distribution.
distribution = ot.ComposedDistribution([ot.Uniform(x_min, x_max),
                                        ot.Uniform(y_min, y_max)])

# Initialize.
ot.RandomGenerator.SetSeed(0)
p_all = np.zeros((n_cases, n_points, 2))
i = 0

# Generate the points.
while i < n_cases:
    # Create a new Latin Hypercube experiment.
    experiment = ot.LHSExperiment(distribution, n_points, False, True)
    # Generate the points.
    dat = np.array(experiment.generate())
    # Check the distance.
    dist = np.min(pdist(dat))
    if dist >= min_dist:
        p_all[i, :, :] = dat
        i += 1
        print('Case {} / {}'.format(i, n_cases))
    else:
        print('Case discarded, min dist = {:.3f}'.format(dist))


# %% Plot.

def init():
    pathcol.set_offsets([[], []])
    return [pathcol]


def update(i, pathcol, points):
    pathcol.set_offsets(points[i])
    return [pathcol]


fig, ax = plt.subplots()
ax.grid()
ax.set_aspect('equal')
ax.set_xlim(1.1 * x_min, 1.1 * x_max)
ax.set_ylim(1.1 * y_min, 1.1 * y_max)
ax.plot([x_min, x_max, x_max, x_min, x_min],
        [y_min, y_min, y_max, y_max, y_min],
        color=(0.3, 0.3, 0.3))

pathcol = ax.scatter([], [])
anim = animation.FuncAnimation(
    fig, update, init_func=init, fargs=(pathcol, p_all), interval=1000, frames=n_cases,
    blit=True, repeat=False)
plt.show()

对于 3D 域，我们可以这样写：

# Boundaries.
z_min, z_max = -20, +20

# Uniform distribution.
distribution = ot.ComposedDistribution([ot.Uniform(x_min, x_max),
                                        ot.Uniform(y_min, y_max),
                                        ot.Uniform(z_min, z_max)])

# Initialize.
p_all = np.zeros((n_cases, n_points, 3))

使用 lhsdesign 和 pdist 在 MATLAB 中进行翻译应该很简单。

正如 Han-Kwang Nienhuys 所指出的，在某些情况下可能存在性能问题，在这种情况下应寻求其他解决方案。

Answer 4

其他几个解决方案在框中取n 随机点，如果任何一对太接近，则丢弃整个集合。这效率不高；找到最坏的违规点并移动它更有效，如下所示：

import numpy as np

def get_sphere_distribution(n, dmin, Ls, maxiter=1e4, allow_wall=True):
    """Get random points in a box with given dimensions and minimum separation.
    
    Parameters:
      
    - n: number of points
    - dmin: minimum distance
    - Ls: dimensions of box, shape (3,) array 
    - maxiter: maximum number of iterations.
    - allow_wall: whether to allow points on wall; 
       (if False: points need to keep distance dmin/2 from the walls.)
        
    Return:
        
    - ps: array (n, 3) of point positions, 
      with 0 <= ps[:, i] < Ls[i]
    - n_iter: number of iterations
    - dratio: average nearest-neighbor distance, divided by dmin.
    
    Note: with a fill density (sphere volume divided by box volume) above about
    0.53, it takes very long. (Random close-packed spheres have a fill density
    of 0.64).
    
    Author: Han-Kwang Nienhuys (2020)
    Copying: BSD, GPL, LGPL, CC-BY, CC-BY-SA
    See Stackoverflow: https://stackoverflow.com/a/62895898/6228891 
    """
    Ls = np.array(Ls).reshape(3)
    if not allow_wall:
        Ls -= dmin
    
    # filling factor; 0.64 is for random close-packed spheres
    # This is an estimate because close packing is complicated near the walls.
    # It doesn't work well for small L/dmin ratios.
    sphere_vol = np.pi/6*dmin**3
    box_vol = np.prod(Ls + 0.5*dmin)
    fill_dens = n*sphere_vol/box_vol
    if fill_dens > 0.64:
        msg = f'Too many to fit in the volume, density {fill_dens:.3g}>0.64'
        raise ValueError(msg)
    
    # initial try   
    ps = np.random.uniform(size=(n, 3)) * Ls
    
    # distance-squared matrix (diagonal is self-distance, don't count)
    dsq = ((ps - ps.reshape(n, 1, 3))**2).sum(axis=2)
    dsq[np.arange(n), np.arange(n)] = np.infty

    for iter_no in range(int(maxiter)):
        # find points that have too close neighbors
        close_counts = np.sum(dsq < dmin**2, axis=1)  # shape (n,)
        n_close = np.count_nonzero(close_counts)
        if n_close == 0:
            break
        
        # Move the one with the largest number of too-close neighbors
        imv = np.argmax(close_counts)
        
        # new positions
        newp = np.random.uniform(size=3)*Ls
        ps[imv]= newp
        
        # update distance matrix
        new_dsq_row = ((ps - newp.reshape(1, 3))**2).sum(axis=-1)
        dsq[imv, :] = dsq[:, imv] = new_dsq_row
        dsq[imv, imv] = np.inf
    else:
        raise RuntimeError(f'Failed after {iter_no+1} iterations.')

    if not allow_wall:
        ps += dmin/2
    
    dratio = (np.sqrt(dsq.min(axis=1))/dmin).mean()
    return ps, iter_no+1, dratio

我尝试了不同的策略来决定移动哪些策略（另请参阅此答案的编辑历史）；似乎每次迭代移动一个坏点比尝试一次移动多个坏点要有效得多。它使填充密度收敛到 0.25 和 0.53 之间存在差异。

dmin=1 和盒子尺寸 LxLxL 的基准测试如下； 3.33 对应于问题中的大小。每个框大小的最高 n 是在 3e+4 次迭代内收敛的最后一个。 d_nn 列是平均最近邻距离。

        L    n  n_iter  d_nn  density
0    3.33   10       9  1.39    0.123
3    3.33   20      40  1.16    0.246
7    3.33   40    5510  1.06    0.493
8    3.33   43    2591  1.03    0.530

9    5.70   10       2  1.73    0.033
14   5.70   60      45  1.22    0.199
16   5.70  150    4331  1.05    0.499
17   5.70  165   25690  1.06    0.549

18  10.00   40       3  1.97    0.030
22  10.00   70      10  1.71    0.053
23  10.00  150      47  1.40    0.113
24  10.00  300     276  1.19    0.225
25  10.00  600    4808  1.07    0.451
27  10.00  720   26418  1.05    0.541

密度值是直径为dmin 的球体的密度（近似于校正壁效应）。随机填充球体的限制约为 0.64；这个算法显然不如把弹珠扔进盒子里。

请注意，问题要求在L=3.33*dmin 框中输入n=200，这将是1.84 的包装密度，并且永远不适合。

作为比较，如果有任何接近的对，则丢弃整个试验解决方案的算法的基准：

        L   n  n_iter  d_nn  density
0    3.33  10      52  1.38    0.123
1    3.33  14     224  1.25    0.172
2    3.33  15   24553  1.15    0.185

3    5.70  10       2  2.02    0.033
4    5.70  20      93  1.42    0.066
5    5.70  29    2089  1.36    0.096
6    5.70  30   10886  1.33    0.100

7   10.00  40      10  2.05    0.030
8   10.00  60    1120  1.79    0.045
9   10.00  68    5521  1.66    0.051
11  10.00  70   28545  1.64    0.053

这个算法需要更多的迭代，而且迭代本身更慢，因为每次迭代都需要重新生成和评估所有点和所有距离。

Answer 5

我通过在网格中生成一个随机点、从随机选项范围中删除该点和相邻点并选择另一个点来解决了类似的问题。