产生在球面上均匀分布的点

如何产生在球面上均匀分布的点呢？ 这里提供若干种思路。

1. 生成两个随机分布的角度，并且按照球坐标的形式产生。

缺点： 产生的点是按照角度均匀分布的， 但是注意这并不是球面上均匀分布的点，后面可以看到两极的地方角度明显密集。并且，由于在计算过程中有大量的三角函数的计算，程序的效率不高。 

 要求不高时可以采用这种方法。

思路是，采用10807 方法产生一组两个随机数，分别作为角度使用。将产生的随机数输出到plt 文件中，使用tecplot绘图。（tecplot是非常好用的CFD后处理可视化的软件，强烈推荐使用）。关于10807 方法产生随机数，有空就另外再开一篇帖子。

下面附上 C++  代码：

#include <iostream>
#include<cmath>
#include<fstream>

using namespace std;

class cRandom
{
    public:
        cRandom(int x,double y):seed(x),random(y){};
        cRandom():seed(0),random(0){};

        int seed;
        double random;
};

cRandom my_random(int z)
// 16807 way to create random numbers
// z is the seed number, num is the total random number to create
{
    //z(n+1)=(a*z(n)+b) mod m
    //describe m=a*q+r to avoid that the number is large than the computer can bear
    const int m=pow(2,31)-1;
    const int a=16807;
    const int q=127773;
    const int r=2836;

    int temp=a*(z%q)-r*(z/q);

    if(temp<0)
    {
        temp=m+temp;
    }
    //z is the seed number
    z=temp;
    double t = z*1.0/m;

    cRandom cr;
    cr.random=t;
    cr.seed=z;

    return cr;
}

int main()
{
    cout << "Hello world!" << endl;
    const int num = pow(10,5);
    int z1 = 20;
    int z2=112;

    ofstream fcout;
    fcout.open("result.plt");
    fcout<<" title=random  "<<endl;
    fcout<<"  VARIABLES = "X","Y ","Z " "<<endl;
    fcout<<"zone I= "<<num<<",datapacking=POINT"<<endl;
    cRandom sita(z1,0.0);
    cRandom pesi(z2,0.0);
    const double pi = 3.141592;

    for(int i=0;i!=num;++i)
    {
        sita=my_random(pesi.seed);
        pesi=my_random(sita.seed);

        double x=1.0*sin(pi*sita.random)*cos(2*pi*pesi.random);
        double y=1.0*sin(pi*sita.random)*sin(2*pi*pesi.random);
        double z=1.0*cos(pi*sita.random);

        fcout<<x<<" "<<y<<" " <<z<<endl;
    }


    fcout.close();
    fcout<<"this program is finished"<<endl;

    return 0;
}

 效果图是这样滴：

             

可以看到在两级的地方点明显的密集。

 不要担心，我们还有更漂亮更快速的方式产生球面上的均匀的点。

2. 三维球面上的Marsaglia 方法

这是一种基于变换抽样的方法，产生的结果非常漂亮。原理我会单开一篇帖子讲述，这里仅仅给出实行的方法。

 step1：
　　　　随机抽样产生一对均匀分布的随机数 u ，v   ；这里u，v 在[-1,1] 范围内
　step2 ：
　　　　计算  r^2 = u^2+v^2; 
　　　　如果 r^2 > 1 则重新抽样，直到满足   r^2 < 1  .
  step3 ：
　　　　计算
　　
　　　　x=2*u*sqrt(1-r2);

　　　　y=2*v*sqrt(1-r2);

　　　　z=1-2*r2;　　

 好了，基本原理非常简单，直接上代码：

#include <iostream>
#include<cmath>
#include<fstream>

using namespace std;

class cRandom
{
    public:
        cRandom(int x,double y):seed(x),random(y){};
        cRandom():seed(0),random(0){};

        int seed;
        double random;
};

cRandom my_random(int z)
// 16807 way to create random numbers
// z is the seed number, num is the total random number to create
{
    //z(n+1)=(a*z(n)+b) mod m
    //describe m=a*q+r to avoid that the number is large than the computer can bear
    const int m=pow(2,31)-1;
    const int a=16807;
    const int q=127773;
    const int r=2836;

    int temp=a*(z%q)-r*(z/q);

    if(temp<0)
    {
        temp=m+temp;
    }
    //z is the seed number
    z=temp;
    double t = z*1.0/m;

    cRandom cr;
    cr.random=t;
    cr.seed=z;

    return cr;
}

int main()
{
    cout << "Hello world!" << endl;
    const int num = pow(10,5);
    int z1 = 20;
    int z2=112;

    ofstream fcout;
    fcout.open("result.plt");
    fcout<<" title=random  "<<endl;
    fcout<<"  VARIABLES = "X","Y ","Z " "<<endl;
    fcout<<"zone I= "<<num/2<<",datapacking=POINT"<<endl;
    cRandom sita(z1,0.0);
    cRandom pesi(z2,0.0);


    for(int i=0;i!=num;++i)
    {
        sita=my_random(pesi.seed);
        pesi=my_random(sita.seed);

        double u=2*sita.random-1.0;
        double v=2*pesi.random-1.0;

        double r2=pow(u,2)+pow(v,2);
        if(r2<1)
        {
           double x=2*u*sqrt(1-r2);
           double y=2*v*sqrt(1-r2);
           double z=1-2*r2;

           fcout<<x<<" "<<y<<" " <<z<<endl;
        }
    }

    fcout.close();
    fcout<<"this program is finished"<<endl;

    return 0;
}

效果图是这样的（360°无死角均匀）：

 

3. 直接抽样法产生球面上均匀分布的点

这种方法是使用直接抽样法，产生球面上随机分布的点，与变换抽样法不同。

#include <iostream>
#include <fstream>
#include<cmath>
using namespace std;

double is_independence(const int N);
    
void my_random(ofstream & fcout,const int num,int z)
// 16807 way to create random numbers
// z is the seed number, num is the total random number to create
{
    //z(n+1)=(a*z(n)+b) mod m
    //describe m=a*q+r to avoid that the number is large than the computer can bear
    const int m=pow(2,31)-1;
    const int a=16807;
    const int q=127773;
    const int r=2836;

    for (int i=0;i!=num;++i)
    {
        int temp=a*(z%q)-r*(z/q);
    
        if(temp<0)
        {
            temp=m+temp;
        }
        //z is the seed number
        z=temp;
        double t = z*1.0/m;
         
        fcout<<t<<endl;
    }
}

void my_randomz(ofstream & fcout,const int num,int z)
// another way to create random numbers
{
    const int m=100000001;
    const int a=329;
    const int q=303951;
    const int r=122;

    for (int i=0;i!=num;++i)
    {
        int temp=a*(z%q)-r*(z/q);

        if(temp<0)
        {
            temp=m+temp;
        }
        z=temp;
        double t = z*1.0/m;

        fcout<<t<<endl;
    }
}

double is_uniform(const int N,const int K)
//judge whether the random is uniform
{

    ifstream fcin;
    fcin.open("re2.dat");
    double esp(0);
    double total(0);
    for (int i = 0; i != N; ++ i)
    {
        double temp;
        fcin>>temp;
        total += pow(temp,K);
    }
    esp = total/N-1.0-1/(K+1);
    fcin.close();    
    return esp;
}

double is_independence(const int N)
//judge whether the randoms are independence
{
    ifstream fcin;
    fcin.open("re2.dat");    // open the file
    if(!fcin)
    {
        cout<<"error! :open file error!"<<endl;
        return -1.0;
    }

    double esp(0);
    double xSq(0);
    double x(0);
    double xy(0);
    
    double temp(0);      //start to input the random number
    fcin>>temp;
    double oldTemp(0);

    for(int i=0;i!=N;++i)
    {
        oldTemp=temp;

        x+=temp;
        xSq+=temp*temp;

        fcin>>temp;
        xy +=temp*oldTemp;
    }
    
    fcin.close();

    double cResult(0);
    cResult=(xy-x*x)/(xSq-x*x);
    return cResult;
}

int main()
{
    cout<<"hello,world!"<<endl;
    int num=pow(10,7);  //the total number of randoms
    int z=45;     //the seed of the program
/*
    // get random number using 10807 and another way
    //save the result in two files

    ofstream fcout;
    fcout.open("re1.dat");
    my_random(fcout,num,z); //call 16807 way to create random number
    fcout.close();
    fcout.open("re2.dat");
    my_randomz(fcout,num,z); //call 16807 way to create random number
    fcout.close();
*/
    
    //const int N = pow(10,3);  // the number of the random
    int N(0);
    cout<<"inter the number of the checked number:"<<endl;
    cin >> N;

    const int K = 10;  //the section number divide
    double esp(0);
    esp=is_uniform(N,K);
    cout<<"uniform number is "<<"esp="<<esp<<endl;

    esp=is_independence(N);
    cout<<"indenpendence number is "<<"esp="<<esp<<endl;

    return 0;
}

 

4.力学方法产生。

这种方法是基于力学原理产生的，球面上的点是主动运动的，他们相互排斥，直到一个最均匀的状态。

这里是传送门： http://zhidao.baidu.com/link?url=hN6nj60BC0RvMZorGtT6VToPH2T9cXcC26MwL2G94e_nFSUPXUstInmXxwydCx-0PI3l4jlKnCv2ZvyLrVsA4oe3b-3PldLS2V2D-hcH0OO

思路是： 对球面上的每一个点施加 坐标： x,y,z; 施加速度，vx,vy,vz ; 计算其余点对当前点的作用力，求解一个微分方程组。初始状态随机分布，每一个计算步更新当前粒子的运动速度和位移。

传送门的matlab代码：

function []=main()
    N=12; %点数量
    a=rand(N,1)*2*pi;  %根据随机求面均匀分布，先生成一个初始状态
    b=asin(rand(N,1)*2-1);
    r0=[cos(a).*cos(b),sin(a).*cos(b),sin(b)];
    v0=zeros(size(r0));
    G=1e-2;%斥力常数，试验这个值比较不错
 
   for ii=1:200%模拟200步，一般已经收敛，其实可以在之下退出
         [rn,vn]=countnext(r0,v0,G);%更新状态
          r0=rn;v0=vn;
    end

    plot3(rn(:,1),rn(:,2),rn(:,3),'.');hold on;%画结果
    [xx,yy,zz]=sphere(50); 
    h2=surf(xx,yy,zz); %画一个单位球做参考
    set(h2,'edgecolor','none','facecolor','r','facealpha',0.7);
    axis equal;
    axis([-1 1 -1 1 -1 1]);
    hold off;

 end

function [rn vn]=countnext(r,v,G) %更新状态的函数
%r存放每点的x，y，z数据，v存放每点的速度数据
    num=size(r,1);
    dd=zeros(3,num,num); %各点间的矢量差
    
    for m=1:num-1
          for n=m+1:num
              dd(:,m,n)=(r(m,:)-r(n,:))';
              dd(:,n,m)=-dd(:,m,n);
          end
     end 

      L=sqrt(sum(dd.^2,1));%各点间的距离
      L(L<1e-2)=1e-2; %距离过小限定
      F=sum(dd./repmat(L.^3,[3 1 1]),3)';%计算合力

      Fr=r.*repmat(dot(F,r,2),[1 3]); %计算合力径向分量
      Fv=F-Fr; %切向分量

       rn=r+v;  %更新坐标
       rn=rn./repmat(sqrt(sum(rn.^2,2)),[1 3]);
       vn=v+G*Fv;%跟新速度
end