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

Fast uniformly distributed random points on the surface of a unit hemisphere

本文关键字：分布 随机 表面上 单位      更新时间：2023-10-16

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

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

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

注意:dsfmt* is将返回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);

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

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

我找到的每一个参考文献都表明，这些方法确实产生了均匀分布，但我没有得到正确的结果。

是否有一个错误在我的实现或我错过了一个基本的想法在第二个和第三个方法?

在单位球(无论其维数是多少)上生成均匀分布的最简单方法是绘制独立的正态分布并将结果向量归一化。

事实上，例如在维度3中，e^(-x^2/2) e^(-y^2/2) e^(-z^2/2) = e^(-(x^2 + y^2 + z^2)/2)所以联合分布通过旋转是不变的

如果你使用快速正态分布生成器(Ziggurat或ratio - of - uniform)和快速归一化程序(谷歌"快速反平方根")，这是快速的。不需要超越函数调用。

同时，Marsaglia在半球上也是不均匀的。你会在赤道附近有更多的点，因为二维圆盘上的对应点<->半球体上的点不是等距的。最后一个似乎是正确的(但是我没有做计算来确保这一点)。

如果你取高度为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)也可以节省一些时间——参见这个stackoverflow问题进行讨论。

编辑添加:好，我现在看到，这只是你的第三种方法的一个小调整。但它减少了一个步骤。

如果你有一个快速的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块内。

你试过摆脱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);

我认为您遇到的非均匀结果的问题是因为在极坐标中，圆上的随机点在径向轴上不均匀分布。如果你看[theta, theta+dtheta]x[r,r+dr]的面积，对于固定的theta和dtheta, r的不同值的面积会有所不同。从直觉上看，离中心更远的地方有"更多的区域"。因此，您需要缩放随机半径来考虑这一点。我没有得到证据，但缩放是r=R*sqrt(rand), R是圆的半径，rand开始是随机数。

第二和第三种方法实际上在球体表面上产生均匀分布的随机点，第二种方法(Marsaglia 1972)产生最快的运行时间，大约是英特尔至强2.8 GHz四核处理器速度的两倍。

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

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

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

第一次尝试(错误)

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

相关文章：