https://www.cnblogs.com/cofludy/p/5894270.html
如何产生在球面上均匀分布的点呢？ 这里提供若干种思路。
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;
}
效果图是这样滴：
image.png
可以看到在两级的地方点明显的密集。
不要担心，我们还有更漂亮更快速的方式产生球面上的均匀的点。
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°无死角均匀）：
image.png
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.力学方法产生。
这种方法是基于力学原理产生的，球面上的点是主动运动的，他们相互排斥，直到一个最均匀的状态。
思路是： 对球面上的每一个点施加 坐标： 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
网友评论