在球面上随机均匀分布点的算法

在球面上随机均匀分布点的算法

文/Sobereva(2) 2015-Dec-12


写个程序时遇到一个问题,需要得到在一个单位球面上随机分布一批点的坐标,想了想办法,第一个办法是选取一个起始点(1,0,0),然后依次绕着X、Y、Z轴随机旋转一定的角度,Fortran代码如下
x=1
y=0
z=0
!Step 1: Rotate about Z axis
call RANDOM_NUMBER(rotval)
rotval=rotval*2*pi
xtmp=cos(rotval)*x-sin(rotval)*y
ytmp=sin(rotval)*x+cos(rotval)*y
x=xtmp
y=ytmp
!Step 2: Rotate about Y axis
call RANDOM_NUMBER(rotval)
rotval=rotval*2*pi
xtmp=cos(rotval)*x-sin(rotval)*z
ztmp=sin(rotval)*x+cos(rotval)*z
x=xtmp
z=ztmp
!Step 3: Rotate about X axis
call RANDOM_NUMBER(rotval)
rotval=rotval*2*pi
ytmp=cos(rotval)*y-sin(rotval)*z
ztmp=sin(rotval)*y+cos(rotval)*z
y=ytmp
z=ztmp

点不是很多的话,看起来还蛮均匀,能用,但点数多了就会发现点的分布并不均匀,100000个点时如下所示,可以隐约看到有一个环状分布比较密集。



出现此问题的原因是,虽然第一步得到了一个均匀随机的环状分布,但是第二步的时候在球的两极已经有点的聚集现象了,第三步时把这个问题弱化了,但是聚集的两极经过绕着X轴转一圈后,还是略微形成了环状密集分布区域。


感觉不够完美,又在网上找找有没有什么其它方法,看到一个做法,对应的Fortran代码:
call RANDOM_NUMBER(rotval)
theta = pi*rotval
call RANDOM_NUMBER(rotval)
phi = 2*pi*rotval
x=sin(theta)*cos(phi)
y=sin(theta)*sin(phi)
z=cos(theta)
试了下,发现很坑爹,点在两极严重聚集。给出这样算法的人真是不负责任,10000个点时分布如下所示




于是搜英文资料,发现一个页面很好,http://mathworld.wolfram.com/SpherePointPicking.html,专门介绍了获得球面上随机分布点的算法,没想到还专门有算法解决这个问题。其中最简单实用的Marsaglia的方法,非常好。过程是,在(-1,1)区间内随机取两个值x1和x2,若这两个值的平方和大于等于1则重新选取。然后球面上的点的x、y、z的坐标用x1和x2就可以很简单得到,Fortran代码如下
do while(.true.)
    call RANDOM_NUMBER(x1)
    call RANDOM_NUMBER(x2)
    x1=2*(x1-0.5D0)
    x2=2*(x2-0.5D0)
    if (x1**2+x2**2<1) exit
end do
x=2*x1*dsqrt(1-x1**2-x2**2)
y=2*x2*dsqrt(1-x1**2-x2**2)
z=1-2*(x1**2+x2**2)

10000个点的时候分布如下,非常均匀完美

已有 5 条评论

  1. xchen

    这个问题我在做分子随机转动时候遇到过,如果按照欧拉角转动,生成三个角的均匀分布,就会导致第二种情况出现的问题,大部分分布在两极,如果随机生成均匀分布的是cos(theta),可以避免上述问题,但是算法复杂。

    wolfram上讨论的是Euler–Rodrigues formula方法,不做转动矩阵的话用不到。

  2. Junjie Li

    我也玩过这个~ 是给每个点一个电荷,搞个简单的动力学演化求个收敛,还挺凑合的

  3. Anonymous

    有这么麻烦么?对x,y,z独立随机生成,然后对坐标归一化到1/sqrt(x^2+y^2+z^2), 这有逻辑缺陷么?

    1. Anonymous

      当然,必须是高斯分布。。。不能是均匀分布

    2. sobereva

      Marsaglia丝毫不麻烦,用你的做法只会更麻烦。你的方法是否可行一试便知。

评论已关闭