[Classify] Algorithm Implementation??

[? Level ] ★★★

[Abstract] In this article, I use mid point displacement method ,to implement the height generating of the terrain.

[Key Words]fractal、terrain generating、mid point method

?

?

[0引言]

?? 用計算機來摸擬自然景物與現像，如山地層巒疊嶂，大海的波濤洶涌，云彩的萬千變幻，是一個十分吸引人的High Order，要完成這樣的一個模擬：涉及到計算機圖形學、隨機數學、偏微分方程數值解法，以及相關的物理知識等多方面的知識，還要再加點藝術感，這可真算不上是件輕松的活，但是，最后的結果是如此的吸引人，以致于讓很多人對此樂此不疲。

[1原理]

?? 山體地形生成的主要方法有:傅立葉濾波法(Fourier filtering),中點位移法(midpoint displacement)等，其中，中點位移法是標準確的分形幾何的方法，

它是基于Mandelbrot和Van Ness 于1968 年提出的一類一維高斯隨機過程分形布朗運動(Fractional Brownian Motion)而建立的。主要的算法相當好理解，如下：

??????

???????????????????????? ?? 圖1

（1）?? 初始條件，確定初始生成的四個點上的高程值。

Count=0

While Count<Iterate Count

?? 對于當前需細分的所有方格:

????? 令a11,a12,a21,a22為方格的左上，右上，左下，右下的四個點

????? 則可生成新的五個點將當前的方格再細分成四份，生成規則如下：

???????? m1=(a11+a12)/2+r1n

???????? m2=(a12+a22)/2+r2n

???????? m3=(a21+a22)/2+r3n

???????? m4=(a11+a21)/2+r4n

????? 其中，r1n,r2n等為隨機生成的參數，生成規則如下：

r1n=r1n-1*detaT（0，sigma^2）

??????????? r1n-1為第n-1次迭代中的隨機參數，detaT為服從(0,sigma^2)的正態分布的隨機變量。

?? Count<--Count+1

End While

顯然，上述這個過程用遞歸的方式是最容易實現的，如果考慮實時性的要求，

可將遞歸算法用棧改為非遞歸，將會使算法效率得到一定的提高。

[2實現]

?? 服從正態分布隨機變量的算法采用的是乘同余法。

#include <stdio.h>

#include <math.h>

#include <time.h>

//MultipleMethod mid mothed

//Possibility Normal Distribute

?

/*-----Start of generate random num----*/

/*-------------------------------------*/

?

/*Sample Distribution:Norm,use select method.

???? 1.get rand1,rand2;

???? 2.u=2*r1-1,v=2*r2-1,w=u^2+v^2;

???? 3.if w>1 goto 1

???? 4.else? x=u[(-ln(w)/w)]^(1/2);y=v[(-lnw)/w]^(1/2)?????????????

*/

?

float funNorm(float miu,float sigma)

{??

???????? float r1,r2;

???????? float u,v,w;

???????? float x,y;??????

?? ????

???????? do

???????? {

???????? ??????? r1=rand()/(float)32767;

???????? ??????? r2=rand()/(float)32767;

???????? ???

?????????????????? u=2*r1-1;

?????????????????? v=2*r2-1;

??????????????????

?????????????????? w=u*u+v*v;

??????????????????

???????? }while(w>1);

????????

???????? x=u*sqrt(((-log(w))/w));

???????? y=v*sqrt(((-log(w))/w));

????????

???????? return miu+sigma*x;????????????????? //also could return miu+sigma*y;

}

/*-end of generation random num ulti fun.------------------*/

?

#define HILLSIZE 50

#define H??????? 0.5

float gHillHeightMap[HILLSIZE][HILLSIZE];

int gHillHeightMapInt[HILLSIZE][HILLSIZE];

?

void generateFractalHill(int x1,int y1,int x2,int y2,int depth);

?

int main()

{?? ?

???????? int i,j;

???????? FILE * debugFile;

???????? char bDebug;

????????

???????? bDebug=getchar();

????????

???????? if((debugFile=fopen("outputData.txt","w"))==NULL)

???????? ? {

??????????? fprintf(stderr,"open file error!");

??????????? return -1;

????????? } ???

??????? printf("/n");

????????

???? srand((unsigned int)time (NULL));???

???? gHillHeightMap[0][0]=20+50*rand()/(float)32767+funNorm(0,1);

???? gHillHeightMap[0][HILLSIZE-1]=20+50*rand()/(float)32767+funNorm(0,1);

???? gHillHeightMap[HILLSIZE-1][0]=20+50*rand()/(float)32767+funNorm(0,1);

???? gHillHeightMap[HILLSIZE-1][HILLSIZE-1]=20+50*rand()/(float)32767+funNorm(0,1);

????????

???? generateFractalHill(0,0,HILLSIZE-1,HILLSIZE-1,1);

????????

???? //adjust data for useful

???? for(i=0;i<HILLSIZE;i++)

??????? for(j=0;j<HILLSIZE;j++)

?????????? gHillHeightMapInt[i][j]=(int)(10*gHillHeightMap[i][j]);

???????? //output datas

???????? for(i=0;i<HILLSIZE;i++)

????????????? {

???????????????? for(j=0;j<HILLSIZE;j++)

???????????? ??? ?{

???????????????????? fprintf(debugFile,"%d ",gHillHeightMapInt[i][j]);

???????????????????? if(bDebug=='y')

fprintf(stdout,"%d ",gHillHeightMapInt[i][j]);

???????????????? }???

???????????????? fprintf(debugFile,";/r/n");?? //for matlab usage.

???????????????? if(bDebug=='y')

fprintf(stdout,";/r/n");?? //for matlab usage.???

????????????????? ?????? ???

??????????? }??????????????????????????

??????? printf("finished/n");

??????? getchar();

???????? return 0;

}

?????????????????????? //row left,col up,row right,col right ,iter count

void generateFractalHill(int x1,int y1,int x2,int y2,int depth)

{

?? if(x2-x1<=1&&y2-y1<=1) //end condition of recursive.

???????? return;

???? else

???????? {???????????????? gHillHeightMap[x1][(y1+y2)/2]=(gHillHeightMap[x1][y1]+gHillHeightMap[x1][y2])/2+pow(0.5,depth*H)*funNorm(0,1);? //up mid

????????????? ?gHillHeightMap[x2][(y1+y2)/2]=(gHillHeightMap[x2][y1]+gHillHeightMap[x2][y2])/2+pow(0.5,depth*H)*funNorm(0,1);? //down mid????

????????????? ?

????????????? ?gHillHeightMap[(x1+x2)/2][y1]=(gHillHeightMap[x1][y1]+gHillHeightMap[x2][y1])/2+pow(0.5,depth*H)*funNorm(0,1);? //left mid

????????????? ?gHillHeightMap[(x1+x2)/2][y2]=(gHillHeightMap[x1][y2]+gHillHeightMap[x2][y2])/2+pow(0.5,depth*H)*funNorm(0,1);? //right mid

????????????? ? ??????????? ?gHillHeightMap[(x1+x2)/2][(y1+y2)/2]=(gHillHeightMap[x1][y1]+gHillHeightMap[x2][y1]+gHillHeightMap[x1][y2]+gHillHeightMap[x2][y2])/4+pow(0.5,depth*H)*funNorm(0,1);??? //the mid point.

????????????? ?

depth++;

?

//Recursive Part

generateFractalHill(x1,y1,(x1+x2)/2,(y1+y2)/2,depth);//up-left??? martix

generateFractalHill((x1+x2)/2,y1,x2,(y1+y2)/2,depth);//down-left?? martix

generateFractalHill(x1,(y1+y2)/2,(x1+x2)/2,y2,depth);//up-right? martix

generateFractalHill((x1+x2)/2,(y1+y2)/2,x2,y2,depth);//down-right matrix

???????? }

}

[3結果]

?? 實驗結果，在Matlab中的高程成像結果：

???????????????????????????????? ??? 圖2

[參考文獻與網站]

[1] 周麗琨、劉金鵬、陳定方，隨機分形地形生成及其瀏覽，武漢理工在學學報，2001（12）

[2] 聶朝陽，盛建平，基于參數的分形地形可控性研究，計算機工程與科學，2006（1）

?

?