產生在球面上均勻分布的點

如何產生在球面上均勻分布的點呢？ 這里提供若干種思路。

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