logo标签列表

在三维空间中以最小的最近邻距离和最大密度随机地给出样本点。

pythonalgorithmnumpyrandomnearest-neighbor
5

我在一个三维空间中有n个点。 我想要随机采样一个子集，使得所有最近邻距离大于r。子集的大小m未知，但是我希望采样点越密集越好，即最大化m。

有类似的问题，但它们都是关于生成点而不是从给定的点中进行采样。
在三维空间中生成具有最小最近邻距离的随机点

生成最小距离的三维随机点？

假设我有300个随机的三维点，

import numpy as np
n = 300
points = np.random.uniform(0, 10, size=(n, 3))

在最大化m的同时，以最小的最近邻距离r=1获取m个点的子集，最快的方法是什么？

你对近似值感兴趣还是结果必须是最优的？ - orlp
1
此问题可以恰当地描述为“在一组欧几里德三维点的单位圆(球)图中查找最大独立集”。 - orlp
在Jallu, R. K.和Das, G. K.的“改进的单位圆图最大独立集算法”中，作者声称该问题（在2D中，意味着3D）是NP难题，并引用Garey，M.，Johnson，D.的书籍“计算机和难解性：NP完全理论指南”作为来源。 - orlp
一个足够快速的近似对我来说已经足够好了，不需要达到最优。我明白在给定300个点的情况下找到全局最优可能需要很长时间。 - Shaun Han
1
@Yatin 没关系。我刚刚发现当最小距离不为1时，David的答案无法工作，所以不需要把第一个更新放回去。我已经在答案部分发布了一个新的更新。 - Shaun Han
显示剩余4条评论
4个回答
2

可能存在一种高效的双目标逼近方案，但为什么要费心呢，因为整数规划在平均情况下非常快速？

import numpy as np

n = 300
points = np.random.uniform(0, 10, size=(n, 3))

from ortools.linear_solver import pywraplp

solver = pywraplp.Solver.CreateSolver("SCIP")
variables = [solver.BoolVar("x[{}]".format(i)) for i in range(n)]
solver.Maximize(sum(variables))
for j, q in enumerate(points):
    for i, p in enumerate(points[:j]):
        if np.linalg.norm(p - q) <= 1:
            solver.Add(variables[i] + variables[j] <= 1)
solver.EnableOutput()
solver.Solve()
print(len([i for (i, variable) in enumerate(variables) if variable.SolutionValue()]))
OP在他们的问题中添加了一些细节（rev [5]（https://stackoverflow.com/revisions/65654970/5）），试图证明他们为什么接受了你的答案。那是对帖子本身的元评论，不应该成为问题的一部分（因此我将其删除）。由于输出来自运行您的代码，您能否将这些细节添加到您的答案中。这样所有相关方都将满意。 - Sabito stands with Ukraine
你能检查一下我的新测试结果吗？由于某些原因，你的代码似乎只在最小距离为1时才能正常工作。 - Shaun Han
1
这不是最优的子集，但使用 KDTree 优化距离计算可以接近最优且不会花费太长时间：
from scipy.spatial import KDTree
import numpy as np

def space_sample(n = 300, low = 0, high = 10, dist = 1):
    points = np.random.uniform(low, high, size=(n, 3))
    k = KDTree(points)
    pairs = np.array(list(k.query_pairs(dist)))
    
    def reduce_pairs(pairs, remove = []):  #iteratively remove the most connected node
        p = pairs[~np.isin(pairs, remove).any(1)]
        if p.size >0:
            count = np.bincount(p.flatten(), minlength = n)
            r = remove + [count.argmax()]
            return reduce_pairs(p, r)
        else:
            return remove
    
    return np.array([p for i, p in enumerate(points) if not(i in reduce_pairs(pairs))])

subset = space_sample()

迭代地删除最连接的节点并不是最优的（保留了约200个300点），但很可能是最接近最优的最快算法（唯一更快的是随机删除）。您可以尝试使用@njitreduce_pairs来加速该部分（如果我以后有时间会尝试）。
我尝试了你的代码。结果还不错，但是从1000个点中抽样需要很长时间。此外，我需要一个能够明确计算距离的代码，因为在我的实际情况中，我需要考虑周期性边界条件。因此，我不能使用KDTree。 - Shaun Han
0

使用30个给定的点测试@David Eisenstat的答案：

from ortools.linear_solver import pywraplp
import numpy as np

def subset_David_Eisenstat(points, r):
    solver = pywraplp.Solver.CreateSolver("SCIP")
    variables = [solver.BoolVar("x[{}]".format(i)) for i in range(len(points))]
    solver.Maximize(sum(variables))
    for j, q in enumerate(points):
        for i, p in enumerate(points[:j]):
            if np.linalg.norm(p - q) <= r:
                solver.Add(variables[i] + variables[j] <= 1)
    solver.EnableOutput()
    solver.Solve()
    indices = [i for (i, variable) in enumerate(variables) if variable.SolutionValue()]
    return points[indices]

points = np.array(
[[ 7.32837882, 12.12765786, 15.01412241],
 [ 8.25164031, 11.14830379, 15.01412241],
 [ 8.21790113, 13.05647987, 13.05647987],
 [ 7.30031002, 13.08276009, 14.05452502],
 [ 9.18056467, 12.0800735 , 13.05183844],
 [ 9.17929647, 11.11270337, 14.03027534],
 [ 7.64737905, 11.48906945, 15.34274827],
 [ 7.01315886, 12.77870699, 14.70301668],
 [ 8.88132414, 10.81243313, 14.68685022],
 [ 7.60617372, 13.39792166, 13.39792166],
 [ 8.85967682, 12.72946394, 12.72946394],
 [ 9.50060705, 11.43361294, 13.37866092],
 [ 8.21790113, 12.08471494, 14.02824481],
 [ 7.32837882, 12.12765786, 16.98587759],
 [ 8.25164031, 11.14830379, 16.98587759],
 [ 7.30031002, 13.08276009, 17.94547498],
 [ 8.21790113, 13.05647987, 18.94352013],
 [ 9.17929647, 11.11270337, 17.96972466],
 [ 9.18056467, 12.0800735 , 18.94816156],
 [ 7.64737905, 11.48906945, 16.65725173],
 [ 7.01315886, 12.77870699, 17.29698332],
 [ 8.88132414, 10.81243313, 17.31314978],
 [ 7.60617372, 13.39792166, 18.60207834],
 [ 8.85967682, 12.72946394, 19.27053606],
 [ 9.50060705, 11.43361294, 18.62133908],
 [ 8.21790113, 12.08471494, 17.97175519],
 [ 7.32837882, 15.01412241, 12.12765786],
 [ 8.25164031, 15.01412241, 11.14830379],
 [ 7.30031002, 14.05452502, 13.08276009],
 [ 9.18056467, 13.05183844, 12.0800735 ],])

当期望的最小距离为1时：

subset1 = subset_David_Eisenstat(points, r=1.)
print(len(subset1))
# Output: 18

检查最小距离：

from scipy.spatial.distance import cdist
dist = cdist(subset1, subset1, metric='euclidean')
# Delete diagonal
res = dist[~np.eye(dist.shape[0],dtype=bool)].reshape(dist.shape[0],-1)
print(np.min(res))
# Output: 1.3285513450926985

将预期的最小距离更改为2：

subset2 = subset_David_Eisenstat(points, r=2.)
print(len(subset2))
# Output: 10

检查最小距离：

from scipy.spatial.distance import cdist
dist = cdist(subset2, subset2, metric='euclidean')
# Delete diagonal
res = dist[~np.eye(dist.shape[0],dtype=bool)].reshape(dist.shape[0],-1)
print(np.min(res))
# Output: 2.0612041004376223
问题出在这一行代码solver.Add(variables[i] + variables[j] <= r)。其中的r应该保持为1--它要求最多只能有一个点出现在所选子集中。 - David Eisenstat
-1

这可能不是很快，但可以迭代3D距离公式，添加到字典中，排序，然后获取ID。

3D距离公式给定点 (x，y，z) 和 (x1，y1，z1)：

m = ((x-x1)^2 + (y-y1)^2)^0.5
d = (m^2 + (z-z1)^2)^0.5

(d 是总距离，^0.5 表示平方根。)

谢谢，但我知道如何计算距离。我更喜欢使用np.linalg.norm :) - Shaun Han
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接