基于分形随机生成月面数字高程模型

一些前置知识或说明

分形：通常被定义为“一个粗糙或零碎的几何形状，可以分成数个部分，且每一部分都（至少近似地）是整体缩小后的形状”，即具有自相似的性质。

通俗的说，分形即是某一个几何形状，它可以分成若干部分，每一部分放大后和整个图形类似，并且每一部分被分开的子图形也具有此种性质。在本文的应用中只需大概理解即可，不必深究。

Diamond-square algorithm：

即钻石形-正方形算法。算法过程如图所示。用此算法来生成随机而又自然的高程数据。根据分形理论分为若干子问题。每一子问题处理一个正方形内的部分高程数据。其过程如下

对于每个子问题：

1.初始状态为正方形四角高程数据已知，随机范围[-s,s]已知。即图1

2.计算中心点数据，即钻石的尖点（上下左右四个钻石）。计算方法为算出四角平均值之后再加上在[-s,s]内随机的数据。如图2

3.计算钻石顶面的中心点，方法同2。如图3

4.此子问题解决，将此问题分解为四个子问题，即四个小正方形，并把s置为0.5s后重复上述过程。如图4，5

初始化可根据问题具体内容给出四角的高程值以及光滑度，即s的取值。

月面石块以及撞击坑数量的数学模型 ：根据统计数据建立出月球平缓月海内石块数量nor与石块直径d_r的函数关系如下：

\[ lg(nor)=-2.589\times lg(d_r)-3.3955 \]

平缓月海内撞击坑数量noc与撞击坑直径d_c的函数关系如下：

\[ lg(noc)=-2\times lg(d_c)-1\quad ,d_c\leq40\\ lg(noc)=-3\times lg(d_c)+0.62\quad ,d_c>40 \]

月面撞击坑形状的数学模型

上述函数数据来源于：张伍, 党兆龙, 贾阳, 等. 月面数字地形构造方法研究[J]. 航天器环境工程, 2008, 25(4):301-306.感兴趣的读者请自行查阅

生成过程

1. 生成月面地形

   1.1根据Diamond-square algorithm随机填充高程点，生成月面的自然起伏。

2. 生成撞击坑，给出最小撞击坑直径

   2.1 由最小直径，根据上述数学模型计算出撞击坑数量noc

   2.2 在[1,noc]范围内随机生成数量random_n

   2.3 根据random_n反向计算对应直径d_c

   2.4 随机生成撞击坑中心点x,y，并随机此撞击坑类型（新鲜，年轻，年老）

   2.5 根据撞击坑形状的数学模型计算并填充高程点

   2.6 循环上述过程noc次

3. 生成石块

   3.1 方法同上，石块形状模型取标准抛物线

4. 输出高程数据

生成结果

采用此方法生成1km x 1km ，分辨率为0.5m的月面高程模型，输出2000000个xyz点对

生成结果在global mapper中预览如下：

可以看出效果还是可以接受的。

代码

#include<bits/stdc++.h>

using namespace std;
const int MAX=2e3+3;
const double min_rock_d=1.5;
const double min_crater_d=10;
const int length=2e3;
const double resolution=0.5;
const double smoothness=25;
const double h1_d[3][2]= {0.23,0.17,0.11,0.25,0.19,0.13}; //h1_d[i][0] for lower_limit, h1_d[i][1] for upper_limit
const double h2_d[3][2]= {0.022,0.016,0.008,0.06,0.045,0.03}; // ditto


double dem[MAX][MAX];
bool vis[MAX][MAX];
struct node
{
    int xs,xt,ys,yt;
    double s;
    node(int _xs,int _xt,int _ys,int _yt,double _s):xs(_xs),xt(_xt),ys(_ys),yt(_yt),s(_s){}
};

inline double sq(double x)
{
    return x*x;
}

inline double rock(double d,double h,double x,double y)
{
    double res=h-4.0*h*(sq(x)+sq(y))/sq(d);
    if(res<0)return 0.0;
    return res;
}

inline double crater(double d,double h1,double h2,double x,double y)
{
    x=sqrt(sq(x)+sq(y));
    double res;
    double rd=rand()%2+1;
    if(x+rd<=d/2)
        res=4.0*h1/sq(d)*sq(x)+h2-h1;
    else
        res=-4.0*h2/sq(d)*sq(x)+4.0*h2/d*(x);
    if(x+rd>d/2&&res<0)
        return 0.0;
    return res;
}

inline int c_d_n(double d)
{
    if(d<=40.0)
        return int(pow(10.0,-2.0*log10(d)-1.0)*length*length/4);
    else
        return int(pow(10.0,-3.0*log10(d)+0.62)*length*length/4);
}
inline double c_n_d(int n)
{
    double nn=double(n)/sq(double(length/2));
    if(n<65)
        return pow(10.0,(log10(nn)-0.62)/-3.0);
    else
        return pow(10.0,(log10(nn)+1.0)/-2.0);
}
inline int r_d_n(double d)
{
    return int(pow(10.0,-2.589*log10(d)-3.3955)*length*length/4);
}
inline double r_n_d(int n)
{
    double nn=double(n)/sq(double(length/2));
    return pow(10.0,(log10(nn)+3.3955)/-2.589);
}

void fill_rock(double d,double h,int xs,int xt,int ys,int yt,int xm,int ym)
{
    if(xs<0)xs=0;
    if(xt>=length)xt=length-1;
    if(ys<0)ys=0;
    if(yt>=length)yt=length-1;
    for(int i=xs; i<=xt; i++)
        for(int j=ys; j<=yt; j++)
        {
            double xx=double(i-xm)*resolution,yy=double(j-ym)*resolution;
            dem[i][j]+=rock(d,h,xx,yy);
        }
}
void fill_crater(double d,double h1,double h2,int xs,int xt,int ys,int yt,int xm,int ym)
{
    if(xs<0)xs=0;
    if(xt>=length)xt=length-1;
    if(ys<0)ys=0;
    if(yt>=length)yt=length-1;
    for(int i=xs; i<=xt; i++)
        for(int j=ys; j<=yt; j++)
        {
            double xx=double(i-xm)*resolution,yy=double(j-ym)*resolution;
            dem[i][j]+=crater(d,h1,h2,xx,yy);
            vis[i][j]=1;
        }
}

void terrain_genaration(void)
{
    srand(time(0));
    queue<node>Q;
    Q.push(node(0,length-1,0,length-1,smoothness));
    dem[0][0]=dem[0][length-1]=dem[length-1][0]=dem[length-1][length-1]=10.0;
    while(!Q.empty())
    {
        node cur=Q.front();
        Q.pop();
        if(cur.xt-cur.xs<2)continue;
        int xs=cur.xs,xt=cur.xt,ys=cur.ys,yt=cur.yt;
        int mx=(xs+xt)>>1,my=(ys+yt)>>1;
        double rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        double avg=(dem[xs][ys]+dem[xs][yt]+dem[xt][ys]+dem[xt][yt])/4.0;
        dem[mx][my]=avg+rd;
        rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        dem[xs][my]=avg+rd;
        rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        dem[xt][my]=avg+rd;
        rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        dem[mx][ys]=avg+rd;
        rd=double(rand()%int(max(cur.s*2000.0,1.0)))/1000.0-cur.s;
        dem[mx][yt]=avg+rd;
        double nxts=cur.s/2.0;
        Q.push(node(xs,mx,ys,my,nxts));
        Q.push(node(mx,xt,ys,my,nxts));
        Q.push(node(xs,mx,my,yt,nxts));
        Q.push(node(mx,xt,my,yt,nxts));
    }
}

void Creater(int noc,int a,int b)
{
    double minn=1e9+7;
    for(int i=0;i<a;i++)
    {
        int num=rand()%noc+b;
        double d=c_n_d(num);
        minn=min(minn,d);
        int t=rand()%3;
        double h1=h1_d[t][0]+(h1_d[t][1]-h1_d[t][0])*double(rand()%10+1)/10.0;
        h1*=d;
        double h2=h2_d[t][0]+(h2_d[t][1]-h2_d[t][0])*double(rand()%10+1)/10.0;
        h2*=d;
        int r=int(d/resolution);
        int x=rand()%length,y=rand()%length;
        while(vis[x][y])
        {
            x=rand()%length;
            y=rand()%length;
        }
        fill_crater(d,h1,h2,x-r,x+r,y-r,y+r,x,y);
    }
}

int main()
{
    srand(time(0));
    terrain_genaration();
    int noc=c_d_n(min_crater_d);
    cout<<"noc="<<noc<<endl;
    Creater(noc,noc*0.3,noc*2.5);
    Creater(noc,noc*0.2,noc*0.5);
    Creater(noc,noc*0.15,2);
    int nor=r_d_n(min_rock_d);
    cout<<"nor="<<nor<<endl;
    for(int i=0; i<nor; i++)
    {
        int num=rand()%nor+1;
        double d=r_n_d(num);
        double h=double(rand()%50+50)/100.0*d;
        int r=int(d/resolution)/2;
        int x=rand()%length,y=rand()%length;
        fill_rock(d,h,x-r,x+r,y-r,y+r,x,y);
    }
    FILE *fpwrite=fopen("maptest.txt","w");
    for(int i=0;i<500;i++)
        for(int j=0;j<500;j++)
        {
            int si=i*4,ti=i*4+4,sj=j*4,tj=j*4+4;
            for(int k=si;k<ti;k++)
                for(int l=sj;l<tj;l++)
                    fprintf(fpwrite,"%.5f %.5f %.5f\n",double(k)*0.5,double(l)*0.5,dem[k][l]);
        }
}

交流

如果本文有什么错误或漏洞诚邀读者指正。如果有什么问题也欢迎与我交流。

邮箱：SChLiu907@gmail.com