生成每个之间距离最小的3-d随机点？
Posted
技术标签:
【中文标题】生成每个之间距离最小的3-d随机点？【英文标题】：generate 3-d random points with minimum distance between each of them? 【发布时间】：2013-01-24 16:14:41 【问题描述】：那里。 我将使用这个特定字符在 matlab 中生成 10^6 个随机点。 这些点应该在半径为 25 的球体内，它们是 3-D，所以我们有 x、y、z 或 r、theta、phi。 每个点之间有一个最小距离。 首先，我决定生成点然后检查距离，然后省略不具备这些条件的点。但是，它可能会省略很多点。 另一种方法是使用 RSA（随机序列加法），这意味着一个一个地生成点，它们之间的最小距离。例如生成第一个点，然后从点 1 的最小距离中随机生成第二个点。继续直到达到 10^6 个点。 但这需要很多时间，我无法达到 10^6 点，因为为新点搜索适当位置的速度需要很长时间。
我现在正在使用这个程序：
Nmax=10000;
R=25;
P=rand(1,3);
k=1;
while k<Nmax
theta=2*pi*rand(1);
phi=pi*rand(1);
r = R*sqrt(rand(1));
% convert to cartesian
x=r.*sin(theta).*cos(phi);
y=r.*sin(theta).*sin(phi);
z=r.*cos(theta);
P1=[x y z];
r=sqrt((x-0)^2+(y-0)^2+(z-0)^2);
D = pdist2(P1,P,'euclidean');
% euclidean distance
if D>0.146*r^(2/3)
P=[P;P1];
k=k+1;
end
i=i+1;
end
x=P(:,1);y=P(:,2);z=P(:,3); plot3(x,y,z,'.');
如何通过这些条件有效地生成积分？ 谢谢。
【问题讨论】：
Nmax=10000; R=25； P=兰德（1,3）； k=1；而 k我仔细研究了您的算法，并得出结论认为它永远不会起作用 - 至少如果您真的想在该领域获得一百万分，则不会。有一张简单的图片可以解释为什么不这样做 - 这是您需要测试的点数（使用您的 RSA 技术）以获得一个额外的“好”点。如您所见，这仅在几千点处渐近（我针对 200k 点运行了一个稍微快一点的算法来产生这个）：
我不知道您是否曾经尝试过在将它们完美排列时计算出可以放入球体的理论点数，但我开始怀疑这个数字比 1E6 小很多。
可以在here 找到我用来调查此问题的完整代码以及它生成的输出。我从来没有达到我在之前的回答中描述的技术......在你描述的设置中还有太多其他事情发生了。
编辑： 我开始认为，即使有“完美”的安排，也可能无法获得 100 万分。我为自己做了一个简单的模型如下：
假设您从“外壳”(r=25) 开始，并尝试以相等的距离拟合点。如果您将“壳”的面积除以一个“排除盘”（半径为 r_sub_crit）的面积，您会得到该距离处点数的（高）估计值：
numpoints = 4*pi*r^2 / (pi*(0.146 * r^(2/3))^2) ~ 188 * r^(2/3)
下一个“外壳”的半径应小 0.146*r^(2/3) - 但如果您认为这些点的排列非常仔细，您可能会更接近一点.再一次，让我们大方地说贝壳可以比标准更接近 1/sqrt(3)。然后，您可以使用简单的 python 脚本从外壳开始并按照自己的方式进行操作：
import scipy as sc
r = 25
npts = 0
def rc(r):
return 0.146*sc.power(r, 2./3.)
while (r > rc(r)):
morePts = sc.floor(4/(0.146*0.146)*sc.power(r, 2./3.))
npts = npts + morePts
print morePts, ' more points at r = ', r
r = r - rc(r)/sc.sqrt(3)
print 'total number of points fitted in sphere: ', npts
这个的输出是：
1604.0 more points at r = 25
1573.0 more points at r = 24.2793037966
1542.0 more points at r = 23.5725257555
1512.0 more points at r = 22.8795314897
1482.0 more points at r = 22.2001865995
1452.0 more points at r = 21.5343566722
1422.0 more points at r = 20.8819072818
1393.0 more points at r = 20.2427039885
1364.0 more points at r = 19.6166123391
1336.0 more points at r = 19.0034978659
1308.0 more points at r = 18.4032260869
1280.0 more points at r = 17.8156625053
1252.0 more points at r = 17.2406726094
1224.0 more points at r = 16.6781218719
1197.0 more points at r = 16.1278757499
1171.0 more points at r = 15.5897996844
1144.0 more points at r = 15.0637590998
1118.0 more points at r = 14.549619404
1092.0 more points at r = 14.0472459873
1066.0 more points at r = 13.5565042228
1041.0 more points at r = 13.0772594652
1016.0 more points at r = 12.6093770509
991.0 more points at r = 12.1527222975
967.0 more points at r = 11.707160503
943.0 more points at r = 11.2725569457
919.0 more points at r = 10.8487768835
896.0 more points at r = 10.4356855535
872.0 more points at r = 10.0331481711
850.0 more points at r = 9.64102993012
827.0 more points at r = 9.25919600154
805.0 more points at r = 8.88751153329
783.0 more points at r = 8.52584164948
761.0 more points at r = 8.17405144976
740.0 more points at r = 7.83200600865
718.0 more points at r = 7.49957037478
698.0 more points at r = 7.17660957023
677.0 more points at r = 6.86298858965
657.0 more points at r = 6.55857239952
637.0 more points at r = 6.26322593726
618.0 more points at r = 5.97681411037
598.0 more points at r = 5.69920179546
579.0 more points at r = 5.43025383729
561.0 more points at r = 5.16983504778
542.0 more points at r = 4.91781020487
524.0 more points at r = 4.67404405146
506.0 more points at r = 4.43840129415
489.0 more points at r = 4.21074660206
472.0 more points at r = 3.9909446055
455.0 more points at r = 3.77885989456
438.0 more points at r = 3.57435701766
422.0 more points at r = 3.37730048004
406.0 more points at r = 3.1875547421
390.0 more points at r = 3.00498421767
375.0 more points at r = 2.82945327223
360.0 more points at r = 2.66082622092
345.0 more points at r = 2.49896732654
331.0 more points at r = 2.34374079733
316.0 more points at r = 2.19501078464
303.0 more points at r = 2.05264138052
289.0 more points at r = 1.91649661498
276.0 more points at r = 1.78644045325
263.0 more points at r = 1.66233679273
250.0 more points at r = 1.54404945973
238.0 more points at r = 1.43144220603
226.0 more points at r = 1.32437870508
214.0 more points at r = 1.22272254805
203.0 more points at r = 1.1263372394
192.0 more points at r = 1.03508619218
181.0 more points at r = 0.94883272297
170.0 more points at r = 0.867440046252
160.0 more points at r = 0.790771268402
150.0 more points at r = 0.718689381062
140.0 more points at r = 0.65105725389
131.0 more points at r = 0.587737626612
122.0 more points at r = 0.528593100237
113.0 more points at r = 0.473486127367
105.0 more points at r = 0.422279001431
97.0 more points at r = 0.374833844693
89.0 more points at r = 0.331012594847
82.0 more points at r = 0.290676989951
75.0 more points at r = 0.253688551418
68.0 more points at r = 0.219908564725
61.0 more points at r = 0.189198057381
55.0 more points at r = 0.161417773651
49.0 more points at r = 0.136428145311
44.0 more points at r = 0.114089257597
38.0 more points at r = 0.0942608092113
33.0 more points at r = 0.0768020649149
29.0 more points at r = 0.0615717987589
24.0 more points at r = 0.0484282253244
20.0 more points at r = 0.0372289153633
17.0 more points at r = 0.0278306908104
13.0 more points at r = 0.0200894920319
10.0 more points at r = 0.013860207063
8.0 more points at r = 0.00899644813842
5.0 more points at r = 0.00535025545232
total number of points fitted in sphere: 55600.0
这似乎证实了你真的无法达到一百万，无论你如何尝试......
【讨论】：
我想到的另外一件事——点的密度是 r 的函数——当你生成随机点时，你在 r 中均匀采样，你可能不应该这样做。这将使生成稍微更有效率 - 但你仍然面临上述限制......当你查看 3D 分布时，它非常不均匀（你看过我发布的完整代码吗？） 【参考方案2】：您可以做很多事情来改进您的程序 - 算法和代码。
在代码方面，真正让你慢下来的一件事是，你不仅使用for 循环（这很慢），而且在行中
P = [P;P1];
您将元素附加到数组中。每次发生这种情况，Matlab 都需要找到一个新的地方来放置数据，复制过程中的所有点。这很快变得非常缓慢。使用
预分配数组P = zeros(1000000, 3);
跟踪您到目前为止找到的点数 N，并将您的距离计算更改为
D = pdist2(P1, P(1:N, :), 'euclidean');
至少会解决这个问题...
另一个问题是您检查新点与所有先前找到的点 - 因此，当您有 100 个点时，您检查大约 100x100，对于 1000，它是 1000x1000。然后你可以看到这个算法至少是 O(N^3) ......不计算随着密度的增加你会得到更多“未命中”的事实。 N=10E6 的 O(N^3) 算法至少需要 10E18 个周期；如果您有一台 4 GHz 机器，每次比较只有一个时钟周期，您将需要 2.5E9 秒 = 大约 80 年。您可以尝试并行处理，但这只是蛮力 - 谁想要呢？
我建议您考虑将问题分解为更小的部分（字面意思）：例如，如果您将球体分成大约为您的最大距离大小的小框，并且对于每个框，您跟踪什么点在其中，那么您只需要检查“此”框及其直接邻居中的点 - 总共 27 个框。如果您的盒子的宽度为 2.5 毫米，那么您将拥有 100x100x100 = 1M 的盒子。这看起来很多，但是现在您的计算时间将大大减少，因为（在算法结束时）每个盒子平均只有 1 个点...当然，使用您使用的距离标准，您将在中心附近有更多点，但这是一个细节。
您需要的数据结构是一个 100x100x100 的元胞数组，每个元胞都包含到目前为止“在那个元胞中”找到的好点的索引。元胞数组的问题在于它不适合向量化。相反，如果您有内存，则可以将其分配为 10x100x100x100 的 4D 数组，假设每个单元格不超过 10 个点（如果这样做，则必须单独处理；在这里与我合作...） .对尚未找到的点使用-1 的索引
那么你的支票会是这样的：
% initializing:
bigList = zeros(10,102,102,102)-1; % avoid hitting the edge...
NPlist = zeros(102, 102, 102); % track # valid points in each box
bottomcorner = [-25.5, -25.5, -25.5]; % boxes span from -25.5 to +25.5
cellSize = 0.5;
.
% in your loop:
P1= [x, y, z];
cellCoords = ceil(P1/cellSize);
goodFlag = true;
pointsSoFar = bigList(:, cellCoords(1)+(-1:1), cellCoords(2)+(-1:1), cellCoords(3)+(-1:1));
pointsToCheck = find(pointsSoFar>0); % this is where the big gains come...
r=sum(P1.^2);
D = pdist2(P1,P(pointsToCheck, :),'euclidean'); % euclidean distance
if D>0.146*r^(2/3)
P(k,:) = P1;
% check the syntax of this line...
cci = ind2sub([102 102 102], cellCoords(1), cellCoords(2), cellCoords(3));
NP(cci)=NP(cci)+1; % increasing number of points in this box
% you want to handle the case where this > 10!!!
bigList(NP(cci), cci) = k;
k=k+1;
end
....
我不知道你能不能从这里拿走；如果你不能，请在笔记中说明，我可能在这个周末有时间更详细地编写代码。有一些方法可以通过一些矢量化来加快速度，但很快就会变得难以管理。
我认为在空间中随机放置大量点，然后使用上述方法进行巨大的矢量化剔除，可能是要走的路。但我建议先采取一些小步骤......如果你可以让上述工作正常，你可以进一步优化（数组大小等）。
【讨论】：
【参考方案3】：我发现 the reference - “使用三维元胞自动机模拟脑肿瘤生长动力学”，Ansal 等人 (2000)。
我同意这令人费解 - 直到您意识到一件重要的事情。他们在mm 报告他们的结果，但你的代码是用cm 编写的。虽然这看起来微不足道，但“临界半径”的公式 rc = 0.146r^(2/3) 包含一个常数 0.146，即维度 - 维度是 mm^(1/3)，而不是 cm^(1/3)。
当我在我的 python 代码中进行更改以评估可能的格点的数量时，它会跳跃 10 倍。现在他们声称他们正在使用 0.38 的“干扰限制” - 你真的找不到这个数字更多网站。如果您包括该限制，我预测最多可以找到 200k 个点 - 仍然低于他们的 150 万，但不是那么疯狂。
您可以考虑联系作者与他们讨论这个问题吗？如果您想将我包括在对话中，您可以发送电子邮件至：SO（只需两封信）my handle name dot united states。与我在上面发布链接的域相同...
【讨论】：
好的，感谢您的有用帮助。我一定会联系作者，当然会通过电子邮件通知你他们的答案。 还有一件事——我刚刚注意到我用来生成上面曲线的代码是使用 r = 1;因为我们实际上将使用 r = 250（以 mm 为单位），所以情况会发生很大变化。等我有一分钟的时间，我会进行一些计算。 我还意识到我的代码也在 mm 中。你怎么能说它是厘米？ r=25mm。脑中的肿瘤怎么可能是25厘米？以上是关于生成每个之间距离最小的3-d随机点？的主要内容，如果未能解决你的问题，请参考以下文章