logo标签列表

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

c++mathrandomgeometry
22

我正在尝试在单位球面上生成均匀随机点以用于蒙特卡罗光线追踪程序。当我说均匀时，我的意思是这些点在表面积方面是均匀分布的。我的当前方法是计算指向正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);

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

马萨利亚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);

这两种方法比第一种运行速度快几倍，但使用它们时，我得到的结果表明它们并没有在球面上生成均匀随机点，而是产生了偏向赤道的分布。

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

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

我的实现中是否存在错误，或者我是否遗漏了第二和第三种方法中的基本思想？

1
你需要什么样的精度？使用你的第一个算法并探索加速三角函数可能是最快的方法。例如，你可以预先计算正弦和余弦表，然后使用某种快速查找和插值而不是完全精确的计算。但我不确定舍入误差会如何影响你的结果。 - Daniel Pryden
我曾考虑使用快速三角函数，但我担心模拟相当敏感于角度。预计算三角值表也会增加缓存未命中率。我开始认为第一种方法可能不是均匀的，这就是区别所在。 - cubiclewar
1
你尝试过优化第一种方法吗？调用sin(asin(x))等价于使用可能有符号变化的x。当您知道cos(z)时，可以根据等式cos(z)^2 + sin(z)^2 = 1计算sin(z)。 - quant_dev
你的第一个方法不规范。请看下面我的回答。 - TonyK
https://github.com/ampl/gsl/blob/644e768630841bd085cb7121085a688c4ff424d0/randist/sphere.c#L66 - zwcloud
显示剩余5条评论
7个回答
16

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

实际上，例如在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圆盘上的对应点<->半球上的点不是等距的，因此您会在赤道附近拥有更多点。最后一个似乎是正确的（但我没有进行计算以确保这一点）。

绝对同意。这是我在球面上生成均匀点的首选方法。如果你关心的话，它也可以推广到任意维度... - Ian Ross
确实。推广到任意n是有用的，而且当n变大时会发生有趣的事情，你别无选择。 - Alexandre C.
好的，我也可以考虑这个。在我的研究中，我确实遇到过它，但由于我永远不会查看高于2球体的任何内容，所以我将其忽略了。 - cubiclewar
@TDawg：如果z分量符号有误，你可以简单地翻转它。 - Alexandre C.
4
如果你在高度为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(方位角)和sin(方位角)来节省一些时间 - 可以参考这个stackoverflow问题进行讨论。 编辑添加：好的，我现在看到这只是您第三种方法的微调。但是它省去了一步。
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块内部。
1
这在球体上不是均匀的。它会给靠近对角线的点更多的质量（对应于单位立方体的顶点所在的区域）。 - Alexandre C.
3
你注意到了对条件 radius < 1 的检查吗？ - quant_dev
1
这似乎是Marsaglia方法的一个变体，但运行速度要慢得多，并且给出相同的结果。 - cubiclewar
有趣。因此，当半径> 0且半径<1时，我们将立方体点云截断为球形点云。天才！这很有效！ - cmarangu
2

你尝试过摆脱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);
这会得到正确的结果，正如您所期望的那样（它是使用三角恒等式进行操作），并且可以显著减少运行时间。但是，它仍然不如第二种方法。 - cubiclewar
1

我认为你在非均匀结果方面遇到的问题是由于在极坐标中，圆上的随机点在径向轴上分布不均匀。如果您查看 [theta，theta+dtheta]x[r，r+dr] 上的区域，对于固定的 theta 和 dtheta，不同的 r 值将具有不同的面积。 直观地说，越靠近中心，"面积" 就会更多。 因此，您需要对随机半径进行缩放以解决这个问题。我没有找到证据，但缩放应为r=R*sqrt（rand），其中R 是圆的半径，rand 是指随机数。

OP 正在尝试在 3D 球体表面模拟均匀分布，而不是在 2D 圆形区域上。 - quant_dev
@quant_dev：是的，但是根据我阅读他的代码，他是通过使用极坐标然后计算平面上方的高度来实现的。所以除非我理解有误，否则我的答案仍然相关。 - carlpett
我明白该区域的分数会随着天顶角的变化而改变，这就是为什么在半球范围内不能随意选择天顶角（0，PI / 2），但我不确定对所提出的函数的影响。我以为它们已经考虑到了表面积的变化。 - cubiclewar
1

第二种和第三种方法实际上可以在球面上产生均匀分布的随机点，其中第二种方法（Marsaglia 1972）在Intel Xeon 2.8 GHz四核处理器上的运行速度大约是其他方法的两倍。

正如Alexandre C所指出的那样，还有一种使用正态分布的方法，比我提出的方法更适用于n维球体。

这个链接将为您提供有关在球面上选择均匀分布的随机点的进一步信息。

正如TonyK所指出的那样，我的初始方法并不会产生均匀分布的点，而是在生成随机点时会偏向极点。尽管这符合我试图解决的问题要求，但我只是假设它会生成均匀随机点。如Pablo所建议的那样，通过删除asin()调用可以优化此方法，从而将运行时间缩短约20%。

你测试过使用单位正态分布变量的方法吗？我相信，只要你使用快速正态生成器（比如均匀比率法或者Ziggurat算法），它应该可以与Marsaglia的方法相媲美。 - Alexandre C.
Alexandre，不，我还没有尝试过，因为我发现由于我正在尝试建模的属性，我实际上需要使用一个不均匀的分布（而是在极点处集中）。不过我有另一个可以应用它的领域，很快就会测试它。 - cubiclewar
-1

第一次尝试（错误）

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%的解决方案。

哦，我明白了，为什么不呢？因为你在长对角线方向上会得到更多的分数。 - Luka Rahne
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接