556 lines
14 KiB
C++
556 lines
14 KiB
C++
/********************************************************
|
|
* ██████╗ ██████╗████████╗██╗
|
|
* ██╔════╝ ██╔════╝╚══██╔══╝██║
|
|
* ██║ ███╗██║ ██║ ██║
|
|
* ██║ ██║██║ ██║ ██║
|
|
* ╚██████╔╝╚██████╗ ██║ ███████╗
|
|
* ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝
|
|
* Geophysical Computational Tools & Library (GCTL)
|
|
*
|
|
* Copyright (c) 2022 Yi Zhang (yizhang-geo@zju.edu.cn)
|
|
*
|
|
* GCTL is distributed under a dual licensing scheme. You can redistribute
|
|
* it and/or modify it under the terms of the GNU Lesser General Public
|
|
* License as published by the Free Software Foundation, either version 2
|
|
* of the License, or (at your option) any later version. You should have
|
|
* received a copy of the GNU Lesser General Public License along with this
|
|
* program. If not, see <http://www.gnu.org/licenses/>.
|
|
*
|
|
* If the terms and conditions of the LGPL v.2. would prevent you from using
|
|
* the GCTL, please consider the option to obtain a commercial license for a
|
|
* fee. These licenses are offered by the GCTL's original author. As a rule,
|
|
* licenses are provided "as-is", unlimited in time for a one time fee. Please
|
|
* send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget
|
|
* to include some description of your company and the realm of its activities.
|
|
* Also add information on how to contact you by electronic and paper mail.
|
|
******************************************************/
|
|
|
|
#ifndef _FRACTAL_TOPO_H
|
|
#define _FRACTAL_TOPO_H
|
|
#include <iostream>
|
|
#include <sstream>
|
|
#include <fstream>
|
|
#include <cmath>
|
|
#include <cstring>
|
|
#include <ctime>
|
|
#include <stdlib.h>
|
|
#include <stdio.h>
|
|
#include "ctype.h"
|
|
#ifdef _WINDOWS
|
|
#include "io.h"
|
|
#include "process.h"
|
|
#else
|
|
#include "unistd.h"
|
|
#endif
|
|
|
|
#include "gctl/core.h"
|
|
#include "gctl/utility.h"
|
|
#include "gctl/geometry.h"
|
|
#include "gctl/algorithms.h"
|
|
|
|
using namespace std;
|
|
|
|
struct posi : public gctl::point2dc
|
|
{
|
|
double v;
|
|
};
|
|
|
|
class FracTopo
|
|
{
|
|
public:
|
|
FracTopo();
|
|
~FracTopo();
|
|
|
|
int topo_gener_constant();
|
|
int topo_gener_regular();//生成随机地形
|
|
int topo_gener_random();
|
|
int gauss_filter_regular();
|
|
int gauss_filter_random();
|
|
int output_regular();//输出文件
|
|
int output_random();
|
|
int run(double* range, double* topo, double* para, double* gauss,
|
|
char* filterType, bool bilinear);
|
|
private:
|
|
double xmin,xmax,dx,ymin,ymax,dy;
|
|
double ld,lu,rd,ru;
|
|
int xnum,ynum,order_x,order_y,imax;
|
|
double H;
|
|
double randmax;
|
|
gctl::ellipse_filter_para filter_N;
|
|
gctl::filter_base_e ifarctan;
|
|
|
|
int ntotal;
|
|
int random_num;
|
|
double** topo;
|
|
double** topo_range;
|
|
posi* topo_random;
|
|
};
|
|
|
|
FracTopo::FracTopo()
|
|
{
|
|
topo = nullptr;
|
|
topo_range = nullptr;
|
|
topo_random = nullptr;
|
|
}
|
|
|
|
FracTopo::~FracTopo()
|
|
{
|
|
if (topo!=nullptr)
|
|
{
|
|
for(int j=0; j<ntotal; j++) delete[] topo[j];
|
|
delete[] topo;
|
|
}
|
|
if (topo_range!=nullptr)
|
|
{
|
|
for (int i = 0; i < xnum; i++) delete[] topo_range[i];
|
|
delete[] topo_range;
|
|
}
|
|
if (topo_random!=nullptr)
|
|
{
|
|
delete[] topo_random;
|
|
}
|
|
}
|
|
|
|
int FracTopo::topo_gener_constant()
|
|
{
|
|
xnum = (xmax-xmin)/dx+1;
|
|
ynum = (ymax-ymin)/dy+1;
|
|
|
|
xmax = xmin+(xnum-1)*dx;
|
|
ymax = ymin+(ynum-1)*dy;
|
|
|
|
topo_range = new double* [xnum];
|
|
for(int i=0; i<xnum; i++) topo_range[i] = new double [ynum];
|
|
|
|
double x, y, p1, p2;
|
|
for (int j=0; j<ynum; j++)
|
|
{
|
|
for (int i=0; i<xnum; i++)
|
|
{
|
|
x = xmin+i*dx;
|
|
y = ymin+j*dy;
|
|
p1 = ld + (lu - ld)*(y - ymin)/(ymax - ymin);
|
|
p2 = rd + (ru - rd)*(y - ymin)/(ymax - ymin);
|
|
topo_range[i][j] = p1 + (p2 - p1)*(x - xmin)/(xmax - xmin);;
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int FracTopo::topo_gener_regular()
|
|
{
|
|
int i,j,m,n,R,jmax;
|
|
double temp = 1e10;
|
|
double random;
|
|
|
|
xnum = (xmax-xmin)/dx+1;
|
|
ynum = (ymax-ymin)/dy+1;
|
|
|
|
xmax = xmin+(xnum-1)*dx;
|
|
ymax = ymin+(ynum-1)*dy;
|
|
|
|
order_x = ceil(log(xnum-1)/log(2));
|
|
order_y = ceil(log(ynum-1)/log(2));
|
|
imax = GCTL_MAX(order_x,order_y);
|
|
|
|
ntotal = (int)pow(2.0,imax)+1; //总数据点数为2的imax次方加1
|
|
|
|
topo = new double* [ntotal];//开辟二维动态数组
|
|
for(i=0; i<ntotal; i++) topo[i] = new double [ntotal];
|
|
topo_range = new double* [xnum];
|
|
for(i=0; i<xnum; i++) topo_range[i] = new double [ynum];
|
|
|
|
for (i=0; i<ntotal; i++)//设定地形数据初始值，角点数据必须给出
|
|
{
|
|
for (j=0; j<ntotal; j++)
|
|
{
|
|
if (i==0&&j==0)
|
|
{
|
|
topo[i][j]=ld;//角点初始值；
|
|
}
|
|
else if (i==ntotal-1&&j==0)
|
|
{
|
|
topo[i][j]=lu;//角点初始值；
|
|
}
|
|
else if (i==0&&j==ntotal-1)
|
|
{
|
|
topo[i][j]=rd;//角点初始值；
|
|
}
|
|
else if (i==ntotal-1&&j==ntotal-1)
|
|
{
|
|
topo[i][j]=ru;//角点初始值；
|
|
}
|
|
else
|
|
{
|
|
topo[i][j]=temp;//非角点数值初始化为1e10
|
|
}
|
|
}
|
|
}
|
|
|
|
srand((unsigned) time(nullptr));
|
|
|
|
for (i=1; i<=imax; i++)//开始迭代，生成正方形区域随机地形
|
|
{
|
|
R=int(double(ntotal-1)/pow(2.0,i));
|
|
|
|
jmax=int(pow(4.0,i-1));
|
|
|
|
for (j=1; j<=jmax; j++)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
m=2*R*(j-(ceil((double)j/pow(2.0,i-1))-1) * pow(2.0,i-1))-R;
|
|
n=2*R*(ceil((double)j/pow(2.0,i-1)))-R;
|
|
topo[m][n]=(topo[m-R][n-R]+topo[m+R][n-R]+topo[m-R][n+R]+topo[m+R][n+R])/4+random;
|
|
}
|
|
|
|
for (j=1; j<=jmax; j++)
|
|
{
|
|
m=2*R*(j-(ceil((double)j/pow(2.0,i-1))-1)* pow(2.0,i-1))-R;
|
|
n=2*R*(ceil((double)j/pow(2.0,i-1)))-R;
|
|
|
|
if (topo[m][n-R]==temp)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
if ((n-R)!=0)
|
|
{
|
|
topo[m][n-R]=(topo[m][n-2*R]+topo[m-R][n-R]+topo[m+R][n-R]+topo[m][n])/4+random;
|
|
}
|
|
else
|
|
{
|
|
topo[m][n-R]=(topo[m-R][n-R]+topo[m+R][n-R]+topo[m][n])/3+random;
|
|
}
|
|
}
|
|
|
|
if (topo[m-R][n]==temp)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
if ((m-R)!=0)
|
|
{
|
|
topo[m-R][n]=(topo[m-R][n-R]+topo[m-2*R][n]+topo[m][n]+topo[m-R][n+R])/4+random;
|
|
}
|
|
else
|
|
{
|
|
topo[m-R][n]=(topo[m-R][n-R]+topo[m][n]+topo[m-R][n+R])/3+random;
|
|
}
|
|
}
|
|
|
|
if (topo[m+R][n]==temp)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
if ((m+R)!=(ntotal-1))
|
|
{
|
|
topo[m+R][n]=(topo[m+R][n-R]+topo[m][n]+topo[m+2*R][n]+topo[m+R][n+R])/4+random;
|
|
}
|
|
else
|
|
{
|
|
topo[m+R][n]=(topo[m+R][n-R]+topo[m][n]+topo[m+R][n+R])/3+random;
|
|
}
|
|
}
|
|
|
|
if (topo[m][n+R]==temp)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
if ((n+R)!=(ntotal-1))
|
|
{
|
|
topo[m][n+R]=(topo[m][n]+topo[m-R][n+R]+topo[m+R][n+R]+topo[m][n+2*R])/4+random;
|
|
}
|
|
else
|
|
{
|
|
topo[m][n+R]=(topo[m][n]+topo[m-R][n+R]+topo[m+R][n+R])/3+random;
|
|
}
|
|
}
|
|
}
|
|
|
|
randmax=pow(2.0,-H)*randmax;
|
|
}
|
|
|
|
for (int j=0; j<ntotal; j++)//按预设区域剪裁数组
|
|
{
|
|
for (int i=0; i<ntotal; i++)
|
|
{
|
|
if (i<xnum&&j<ynum)
|
|
{
|
|
topo_range[i][j]=topo[i][j];
|
|
}
|
|
}
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
int FracTopo::topo_gener_random()
|
|
{
|
|
int i,j,m,n,R,jmax;
|
|
double temp = 1e+10;
|
|
double random;
|
|
|
|
dx=dy=ceil(sqrt((xmax-xmin)*(ymax-ymin)/random_num));//按照随机点数预估节点步长
|
|
|
|
xnum = (xmax-xmin)/dx+1;
|
|
ynum = (ymax-ymin)/dy+1;
|
|
|
|
xmax = xmin+(xnum-1)*dx;
|
|
ymax = ymin+(ynum-1)*dy;
|
|
|
|
order_x = ceil(log(xnum)/log(2));
|
|
order_y = ceil(log(ynum)/log(2));
|
|
imax = GCTL_MAX(order_x,order_y);
|
|
|
|
ntotal = (int)pow(2.0,imax)+1; //总数据点数为2的imax次方加1
|
|
|
|
topo = new double* [ntotal];//开辟二维动态数组
|
|
for(i=0; i<ntotal; i++) topo[i] = new double [ntotal];
|
|
topo_random = new posi [random_num];
|
|
|
|
for (i=0; i<ntotal; i++)//设定地形数据初始值，角点数据必须给出
|
|
{
|
|
for (j=0; j<ntotal; j++)
|
|
{
|
|
if (i==0&&j==0)
|
|
{
|
|
topo[i][j]=ld;//角点初始值；
|
|
}
|
|
else if (i==ntotal-1&&j==0)
|
|
{
|
|
topo[i][j]=lu;//角点初始值；
|
|
}
|
|
else if (i==0&&j==ntotal-1)
|
|
{
|
|
topo[i][j]=rd;//角点初始值；
|
|
}
|
|
else if (i==ntotal-1&&j==ntotal-1)
|
|
{
|
|
topo[i][j]=ru;//角点初始值；
|
|
}
|
|
else
|
|
{
|
|
topo[i][j]=temp;//非角点数值初始化为1e10
|
|
}
|
|
}
|
|
}
|
|
|
|
srand((unsigned) time(nullptr));
|
|
|
|
for (i=1; i<=imax; i++)//开始迭代，生成正方形区域随机地形
|
|
{
|
|
R=int(double(ntotal-1)/pow(2.0,i));
|
|
|
|
jmax=int(pow(4.0,i-1));
|
|
|
|
for (j=1; j<=jmax; j++)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
m=2*R*(j-(ceil((double)j/pow(2.0,i-1))-1) * pow(2.0,i-1))-R;
|
|
n=2*R*(ceil((double)j/pow(2.0,i-1)))-R;
|
|
topo[m][n]=(topo[m-R][n-R]+topo[m+R][n-R]+topo[m-R][n+R]+topo[m+R][n+R])/4+random;
|
|
}
|
|
|
|
for (j=1; j<=jmax; j++)
|
|
{
|
|
m=2*R*(j-(ceil((double)j/pow(2.0,i-1))-1)* pow(2.0,i-1))-R;
|
|
n=2*R*(ceil((double)j/pow(2.0,i-1)))-R;
|
|
|
|
if (topo[m][n-R]==temp)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
if ((n-R)!=0)
|
|
{
|
|
topo[m][n-R]=(topo[m][n-2*R]+topo[m-R][n-R]+topo[m+R][n-R]+topo[m][n])/4+random;
|
|
}
|
|
else
|
|
{
|
|
topo[m][n-R]=(topo[m-R][n-R]+topo[m+R][n-R]+topo[m][n])/3+random;
|
|
}
|
|
}
|
|
|
|
if (topo[m-R][n]==temp)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
if ((m-R)!=0)
|
|
{
|
|
topo[m-R][n]=(topo[m-R][n-R]+topo[m-2*R][n]+topo[m][n]+topo[m-R][n+R])/4+random;
|
|
}
|
|
else
|
|
{
|
|
topo[m-R][n]=(topo[m-R][n-R]+topo[m][n]+topo[m-R][n+R])/3+random;
|
|
}
|
|
}
|
|
|
|
if (topo[m+R][n]==temp)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
if ((m+R)!=(ntotal-1))
|
|
{
|
|
topo[m+R][n]=(topo[m+R][n-R]+topo[m][n]+topo[m+2*R][n]+topo[m+R][n+R])/4+random;
|
|
}
|
|
else
|
|
{
|
|
topo[m+R][n]=(topo[m+R][n-R]+topo[m][n]+topo[m+R][n+R])/3+random;
|
|
}
|
|
}
|
|
|
|
if (topo[m][n+R]==temp)
|
|
{
|
|
random=2*randmax*rand()/RAND_MAX-randmax;
|
|
|
|
if ((n+R)!=(ntotal-1))
|
|
{
|
|
topo[m][n+R]=(topo[m][n]+topo[m-R][n+R]+topo[m+R][n+R]+topo[m][n+2*R])/4+random;
|
|
}
|
|
else
|
|
{
|
|
topo[m][n+R]=(topo[m][n]+topo[m-R][n+R]+topo[m+R][n+R])/3+random;
|
|
}
|
|
}
|
|
}
|
|
|
|
randmax=pow(2.0,-H)*randmax;
|
|
}
|
|
|
|
srand((unsigned) time(nullptr));
|
|
int temp_m,temp_n;
|
|
double len_ld,len_lu,len_ru,len_rd;
|
|
for (int i = 0; i < random_num; i++)
|
|
{
|
|
topo_random[i].x = (xmax-xmin)*rand()/RAND_MAX+xmin;
|
|
topo_random[i].y = (ymax-ymin)*rand()/RAND_MAX+ymin;
|
|
temp_m = floor((topo_random[i].x - xmin)/dx);
|
|
temp_n = floor((topo_random[i].y - ymin)/dy);
|
|
len_ld = sqrt(pow(xmin+dx*temp_m - topo_random[i].x,2)+pow(ymin+dy*temp_n - topo_random[i].y,2));
|
|
len_lu = sqrt(pow(xmin+dx*temp_m - topo_random[i].x,2)+pow(ymin+dy*(temp_n+1) - topo_random[i].y,2));
|
|
len_ru = sqrt(pow(xmin+dx*(temp_m+1) - topo_random[i].x,2)+pow(ymin+dy*(temp_n+1) - topo_random[i].y,2));
|
|
len_rd = sqrt(pow(xmin+dx*(temp_m+1) - topo_random[i].x,2)+pow(ymin+dy*temp_n - topo_random[i].y,2));
|
|
topo_random[i].v = (topo[temp_m][temp_n]*len_ld+topo[temp_m+1][temp_n]*len_rd+topo[temp_m+1][temp_n+1]*len_ru+topo[temp_m][temp_n+1]*len_lu)/(len_lu+len_ru+len_rd+len_ld);
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
int FracTopo::output_regular()
|
|
{
|
|
time_t now = time(0);
|
|
char* dt = ctime(&now);
|
|
|
|
cout << "# A fractal terrain data generated by fractopo on " << dt;
|
|
cout << "# For more information please connect zhangyiss@icloud.com" << endl;
|
|
cout << "# Use -h or --help option to see help information" << endl;
|
|
cout << "> Range = " << xmin << "/" << xmax << "/" << ymin << "/" << ymax << endl;
|
|
cout << "> Interval = " << dx << "/" << dy << endl;
|
|
cout << "# x y elevation" << endl;
|
|
for (int j=0; j<ynum; j++)
|
|
{
|
|
for (int i=0; i<xnum; i++)
|
|
{
|
|
cout<<xmin+i*dx<<" "<<ymin+j*dy<<" "<<topo_range[i][j]<<endl;
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int FracTopo::output_random()
|
|
{
|
|
time_t now = time(0);
|
|
char* dt = ctime(&now);
|
|
|
|
cout << "# A fractal terrain data generated by fractopo on " << dt;
|
|
cout << "# For more information please connect zhangyi.cugwuhan@gmail.com" << endl;
|
|
cout << "# Use -h option to see help information" << endl;
|
|
cout << "Count = " << random_num << endl;
|
|
cout << "# x y elevation" << endl;
|
|
for (int i = 0; i < random_num; i++)
|
|
{
|
|
cout<<topo_random[i].x<<" "<<topo_random[i].y<<" "<<topo_random[i].v<<endl;
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int FracTopo::gauss_filter_regular()
|
|
{
|
|
for (int j=0; j<ynum; j++)
|
|
{
|
|
for (int i=0; i<xnum; i++)
|
|
{
|
|
topo_range[i][j] *= gctl::ellipse_filter(xmin+i*dx, ymin+j*dy, filter_N, ifarctan);
|
|
}
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int FracTopo::gauss_filter_random()
|
|
{
|
|
for (int i = 0; i < random_num; i++)
|
|
{
|
|
topo_random[i].v *= gctl::ellipse_filter(topo_random[i].x, topo_random[i].y, filter_N, ifarctan);
|
|
}
|
|
return 0;
|
|
}
|
|
|
|
int FracTopo::run(double* range,double* topo,double* para,double* gauss,char* filterType, bool bilinear)
|
|
{
|
|
if (!strcmp(filterType,"pow"))
|
|
ifarctan = gctl::Eps;
|
|
else if (!strcmp(filterType,"atan"))
|
|
ifarctan = gctl::Atan;
|
|
|
|
if (*(range+5)==1e+30)
|
|
{
|
|
xmin = *range; xmax = *(range+1);
|
|
ymin = *(range+2); ymax = *(range+3);
|
|
random_num = (int) *(range+4);
|
|
ld = *topo; lu = *(topo+1); ru = *(topo+2); rd = *(topo+3);
|
|
H = *para; randmax = *(para+1);
|
|
|
|
if (bilinear) topo_gener_constant();
|
|
else topo_gener_random();
|
|
|
|
if (*gauss!=1e+30)
|
|
{
|
|
filter_N.mu_x = *gauss;
|
|
filter_N.mu_y = *(gauss+1);
|
|
filter_N.sigma = *(gauss+2);
|
|
filter_N.theta = *(gauss+3);
|
|
filter_N.len_x = *(gauss+4);
|
|
filter_N.len_y = *(gauss+5);
|
|
filter_N.magnify = *(gauss+6);
|
|
filter_N.rho = *(gauss+7);
|
|
gauss_filter_random();
|
|
}
|
|
output_random();
|
|
}
|
|
else
|
|
{
|
|
xmin = *range; xmax = *(range+1);
|
|
ymin = *(range+2); ymax = *(range+3);
|
|
dx = *(range+4); dy = *(range+5);
|
|
ld = *topo; lu = *(topo+1); ru = *(topo+2); rd = *(topo+3);
|
|
H = *para; randmax = *(para+1);
|
|
|
|
if (bilinear) topo_gener_constant();
|
|
else topo_gener_regular();
|
|
|
|
if (*gauss!=1e+30)
|
|
{
|
|
filter_N.mu_x = *gauss;
|
|
filter_N.mu_y = *(gauss+1);
|
|
filter_N.sigma = *(gauss+2);
|
|
filter_N.theta = *(gauss+3);
|
|
filter_N.len_x = *(gauss+4);
|
|
filter_N.len_y = *(gauss+5);
|
|
filter_N.magnify = *(gauss+6);
|
|
filter_N.rho = *(gauss+7);
|
|
gauss_filter_regular();
|
|
}
|
|
output_regular();
|
|
}
|
|
return 0;
|
|
}
|
|
#endif |