单位半球表面上快速均匀分布的随机点

Posted 2023-02-17

【中文标题】单位半球表面上快速均匀分布的随机点

【英文标题】：Fast uniformly distributed random points on the surface of a unit hemisphere

【发布时间】：2011-09-02 07:03:59

【问题描述】：

我正在尝试为蒙特卡洛光线追踪程序在单位球体的表面上生成均匀的随机点。当我说均匀时，我的意思是这些点相对于表面积是均匀分布的。我目前的方法是计算指向正 z 轴和 x-y 平面基点的半球上的均匀随机点。

半球上的随机点表示漫反射灰发射器的热辐射发射方向。

当我使用以下计算时，我得到了正确的结果：

注意：dsfmt* 将返回一个介于 0 和 1 之间的随机数。

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
zenith = asin(sqrt(dsfmt_genrand_close_open(&dsfmtt)));

// Calculate the cartesian point
osRay.c._x = sin(zenith)*cos(azimuthal); 
osRay.c._y = sin(zenith)*sin(azimuthal);
osRay.c._z = cos(zenith);

但是，这非常慢，并且分析表明它占用了很大一部分运行时间。因此，我寻找了一些替代方法：

Marsaglia 1972 拒绝方法

do 
   x1 = 2.0*dsfmt_genrand_open_open(&dsfmtt)-1.0;
   x2 = 2.0*dsfmt_genrand_open_open(&dsfmtt)-1.0;
   S = x1*x1 + x2*x2;
 while(S > 1.0f);


osRay.c._x = 2.0*x1*sqrt(1.0-S);
osRay.c._y = 2.0*x2*sqrt(1.0-S);
osRay.c._z = abs(1.0-2.0*S);

解析笛卡尔坐标计算

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
u = 2*dsfmt_genrand_close_open(&dsfmtt) -1;
w = sqrt(1-u*u);

osRay.c._x = w*cos(azimuthal);
osRay.c._y = w*sin(azimuthal);
osRay.c._z = abs(u);

虽然最后两种方法的运行速度比第一种快几倍，但当我使用它们时，我得到的结果表明它们没有在球体表面生成均匀的随机点，而是给出了有利于赤道的分布。

此外，最后两种方法给出了相同的最终结果，但我确信它们是不正确的，因为我正在与分析解决方案进行比较。

我发现的每个参考资料都表明这些方法确实会产生均匀分布，但是我没有得到正确的结果。

我的实现是否有错误，或者我错过了第二种和第三种方法的基本思想？

【问题讨论】：

您需要什么样的精度？使用您的第一个算法并探索加速三角函数的方法可能是最快的。例如，您可以预先计算正弦和余弦表，然后使用某种快速查找和插值代替全精度计算。不过，我不确定舍入误差会如何影响您的结果。

我曾考虑过使用快速触发函数，但我担心模拟对角度相当敏感。预先计算一个三角值表也会导致缓存未命中的增加。我开始认为第一种方法可能不统一，这就是区别。

您尝试优化第一种方法吗？调用 sin(asin(x)) 等同于使用 x 并可能更改符号。当你知道cos(z)，你可以从等式cos(z)^2 + sin(z)^2 = 1计算sin(z)。

第一种方法可以稍作优化，但速度仍然比不上第二种和第三种方法。

你的第一种方法不统一。请参阅下面的答案。

【参考方案1】：

在单位球体（无论其维度是多少）上生成均匀分布的最简单方法是绘制独立的正态分布并对结果向量进行归一化。

确实，例如在维度3，e^(-x^2/2) e^(-y^2/2) e^(-z^2/2) = e^(-(x^2 + y^2 + z^2)/2) 所以联合分布不随旋转变化。

如果您使用快速正态分布生成器（Ziggurat 或 Ratio-Of-Uniforms）和快速归一化例程（谷歌表示“快速平方根平方根”），这会很快。不需要调用超越函数。

此外，Marsaglia 在半球上并不均匀。由于半球体上 2D 圆盘 点上的对应点不是等距的，因此在赤道附近会有更多点。最后一个似乎是正确的（但是我没有进行计算来确保这一点）。

【讨论】：

绝对同意。这是我在球体上均匀点的首选方法。它还可以推广到任意数量的维度，如果你关心的话......

确实如此。对任意 n 的泛化很有用，因为当 n 变大时会发生有趣的事情，你别无选择。

好吧，我也可以考虑一下，我在研究中确实遇到过它，但是因为我永远不会看到 2 球以上的任何东西，所以我忽略了它。

@TDawg：如果 z 组件的符号错误，您可以简单地翻转它。

【参考方案2】：

如果你取一个单位球体的水平切片，高度为h，它的表面积就是2 pi h。 （这是阿基米德计算球体表面积的方法。）所以z坐标均匀分布在[0,1]中：

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
osRay.c._z = dsfmt_genrand_close_open(&dsfmtt);

xyproj = sqrt(1 - osRay.c._z*osRay.c._z);
osRay.c._x = xyproj*cos(azimuthal); 
osRay.c._y = xyproj*sin(azimuthal);

您还可以通过同时计算cos(azimuthal) 和sin(azimuthal) 来节省一些时间——请参阅this *** question 进行讨论。

编辑添加： 好的，我现在知道这只是您的第三种方法的轻微调整。但它减少了一步。

【参考方案3】：

如果你有一个快速的 RNG，这应该很快：

// RNG::draw() returns a uniformly distributed number between -1 and 1.

void drawSphereSurface(RNG& rng, double& x1, double& x2, double& x3)

    while (true) 
        x1 = rng.draw();
        x2 = rng.draw();
        x3 = rng.draw();
        const double radius = sqrt(x1*x1 + x2*x2 + x3*x3);
        if (radius > 0 && radius < 1) 
            x1 /= radius;
            x2 /= radius;
            x3 /= radius;
            return;
        
       

为了加快速度，您可以将 sqrt 调用移动到 if 块内。

【讨论】：

这在球体上并不统一。它赋予对角线附近的点更多的质量（对应于单位立方体顶点的区域）。

你注意到检查条件radius < 1了吗？

这似乎是 Marsaglia 方法的一种变体，但运行速度要慢一些，但结果相同。

有趣。因此，当半径 > 0 和半径

【参考方案4】：

你试过摆脱asin吗？

azimuthal = 2*PI*dsfmt_genrand_close_open(&dsfmtt);
sin2_zenith = dsfmt_genrand_close_open(&dsfmtt);
sin_zenith = sqrt(sin2_zenith);

// Calculate the cartesian point
osRay.c._x = sin_zenith*cos(azimuthal); 
osRay.c._y = sin_zenith*sin(azimuthal);
osRay.c._z = sqrt(1 - sin2_zenith);

【讨论】：

这给出了你所期望的正确结果（它是使用三角函数识别的操作），并且可以减少运行时间。但是它仍然与第二种方法不同。

【参考方案5】：

我认为您遇到的结果不均匀的问题是因为在极坐标中，圆上的随机点在径向轴上分布不均匀。如果您查看[theta, theta+dtheta]x[r,r+dr] 上的区域，对于固定的theta 和dtheta，r 的不同值的区域会有所不同。直观地说，离中心更远的地方有“更多区域”。因此，您需要缩放随机半径以解决此问题。我没有证据，但缩放比例是r=R*sqrt(rand)，R 是圆的半径，rand 是随机数的开头。

【讨论】：

OP 试图在 3D 球体表面上模拟均匀分布，而不是在 2D 圆的区域上。

@quant_dev：是的，但是当我阅读他的代码时，他是通过使用极坐标然后计算平面上方的高度来做到这一点的。所以除非我读错了，否则我的答案仍然是相关的

我知道面积分数随天顶角而变化，这就是为什么天顶角不能在半球的范围 (0,PI/2) 内随机选择的原因，但我不确定其含义关于建议的功能。我认为他们已经考虑了不断变化的表面积。

【参考方案6】：

第二种和第三种方法实际上会在球体表面上产生均匀分布的随机点，第二种方法 (Marsaglia 1972) 在 Intel Xeon 2.8 GHz 四核上以大约两倍的速度产生最快的运行时间.

正如Alexandre C 所指出的，还有一种使用正态分布的方法，它比我提出的方法更好地扩展到 n 球体。

This link 将为您提供有关选择球体表面上均匀分布的随机点的更多信息。

TonyK 指出的我的初始方法不会产生均匀分布的点，而是在生成随机点时偏向极点。这是我要解决的问题所必需的，但是我只是假设它会生成均匀的随机点。正如Pablo 所建议的那样，可以通过删除 asin() 调用来优化此方法，以减少大约 20% 的运行时间。

【讨论】：

您是否使用单位正态分布变量测试了该方法？我相信它应该与 Marsaglia 的方法具有竞争力，只要您使用快速的正常生成器（制服或 Ziggurat 的比率）。

Alexandre，不，我没有尝试过，因为我发现由于我正在尝试建模的属性，我实际上需要使用不均匀的分布（而是集中在两极）。不过，我确实有另一个领域可以应用它，并会很快对其进行测试。

【参考方案7】：

第一次尝试（错误）

point=[rand(-1,1),rand(-1,1),rand(-1,1)];
len = length_of_vector(point);

已编辑：

怎么样？

while(1)
 point=[rand(-1,1),rand(-1,1),rand(-1,1)];
 len = length_of_vector(point);
 if( len > 1 )
     continue;
 point = point / len
     break

接受度约为 0.4。这意味着您将拒绝 60% 的解决方案。

【讨论】：

哦，我明白了，为什么不呢。因为在长对角线方向上会有更多的点。