/******************************************************** * ██████╗ ██████╗████████╗██╗ * ██╔════╝ ██╔════╝╚══██╔══╝██║ * ██║ ███╗██║ ██║ ██║ * ██║ ██║██║ ██║ ██║ * ╚██████╔╝╚██████╗ ██║ ███████╗ * ╚═════╝ ╚═════╝ ╚═╝ ╚══════╝ * 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 . * * 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 _CUTPROFILE_H #define _CUTPROFILE_H #include "gctl/core.h" #include "gctl/utility.h" #include "gctl/algorithms.h" struct node { int id; double px,pz; //切割剖面上的二维横纵坐标值 /*我们在建立盒子的过程中并不专注输入点是否在直角坐标或是球坐标下但在球坐标下必须要指定网格化参数*/ double x,y,z; //在球坐标数据插值中 x y值对应经纬度 z值对应深度或者相对于参考球的高程 double val; //目前只允许每个点附带一个数据值 node() { id = -1; x = y = z = val = GCTL_BDL_MAX; px = pz = GCTL_BDL_MAX; } void sph_init(double lon,double lat,double ele) { px = lon; pz = lat; x = lon; y = lat; z = ele; } void info() { std::cout << id << " " << x << " " << y << " " << z << " " << px << " " << pz << " " << val << std::endl; } }; typedef std::vector nodeArray; node c2s_node(node n,double r,double R) { node m; double rad = gctl::ellipse_radius_2d(R, r, n.y*GCTL_Pi/180.0); // 这里还有疑问应该没问题 m.x = rad*sin((0.5 - n.y/180.0)*GCTL_Pi)*cos((2.0 + n.x/180.0)*GCTL_Pi); m.y = rad*sin((0.5 - n.y/180.0)*GCTL_Pi)*sin((2.0 + n.x/180.0)*GCTL_Pi); m.z = rad*cos((0.5 - n.y/180.0)*GCTL_Pi); return m; } double interCubeDistSph(double r,double R,double x0,double y0,double z0,double dx,double dy,double dz,double x,double y,double z, double d0,double d1,double d2,double d3,double d4,double d5,double d6,double d7) { node tempNode,tempNode2; double res = GCTL_BDL_MAX; double total_dist = 0; double dist[8] = {0,0,0,0,0,0,0,0}; double val[8]; val[0] = d0; val[1] = d1; val[2] = d2; val[3] = d3; val[4] = d4; val[5] = d5; val[6] = d6; val[7] = d7; //检查八个角点值取有效角点的值进行插值如果8个角点均为BDL_MAX 则返回BDL_MAX //八个角点中至少有一个为有效值初始化res为0 注意我们现在这个算法的主要疑点在于可能将角点值广播到整个块体范围内 for (int i = 0; i < 8; i++) { if (val[i] != GCTL_BDL_MAX) { res = 0; break; } } //计算八个角点与待插值点的距离信息 tempNode.x = x-x0; tempNode.y = y-y0; tempNode.z = z-z0; tempNode2 = c2s_node(tempNode,r,R); dist[0] = 1.0/(1e-30+tempNode.x*tempNode.x+tempNode.y*tempNode.y+tempNode.z*tempNode.z); tempNode.x = x-dx-x0; tempNode.y = y-y0; tempNode.z = z-z0; tempNode2 = c2s_node(tempNode,r,R); dist[1] = 1.0/(1e-30+tempNode.x*tempNode.x+tempNode.y*tempNode.y+tempNode.z*tempNode.z); tempNode.x = x-dx-x0; tempNode.y = y-dy-y0; tempNode.z = z-z0; tempNode2 = c2s_node(tempNode,r,R); dist[2] = 1.0/(1e-30+tempNode.x*tempNode.x+tempNode.y*tempNode.y+tempNode.z*tempNode.z); tempNode.x = x-x0; tempNode.y = y-dy-y0; tempNode.z = z-z0; tempNode2 = c2s_node(tempNode,r,R); dist[3] = 1.0/(1e-30+tempNode.x*tempNode.x+tempNode.y*tempNode.y+tempNode.z*tempNode.z); tempNode.x = x-x0; tempNode.y = y-y0; tempNode.z = z-dz-z0; tempNode2 = c2s_node(tempNode,r,R); dist[4] = 1.0/(1e-30+tempNode.x*tempNode.x+tempNode.y*tempNode.y+tempNode.z*tempNode.z); tempNode.x = x-dx-x0; tempNode.y = y-y0; tempNode.z = z-dz-z0; tempNode2 = c2s_node(tempNode,r,R); dist[5] = 1.0/(1e-30+tempNode.x*tempNode.x+tempNode.y*tempNode.y+tempNode.z*tempNode.z); tempNode.x = x-dx-x0; tempNode.y = y-dy-y0; tempNode.z = z-dz-z0; tempNode2 = c2s_node(tempNode,r,R); dist[6] = 1.0/(1e-30+tempNode.x*tempNode.x+tempNode.y*tempNode.y+tempNode.z*tempNode.z); tempNode.x = x-x0; tempNode.y = y-dy-y0; tempNode.z = z-dz-z0; tempNode2 = c2s_node(tempNode,r,R); dist[7] = 1.0/(1e-30+tempNode.x*tempNode.x+tempNode.y*tempNode.y+tempNode.z*tempNode.z); for (int i = 0; i < 8; i++) { if (val[i] != GCTL_BDL_MAX) { total_dist += dist[i]; } } for (int i = 0; i < 8; i++) { if (val[i] != GCTL_BDL_MAX) { res += val[i]*dist[i]/total_dist; } } return res; } class cutProfile { private: nodeArray inputNode; //输入点数组 nodeArray profileNode; //待插值点数组 int xnum,ynum,znum; double dx,dy,dz; double xmin,xmax,ymin,ymax,zmin,zmax; std::vector > boxIndex; nodeArray boxNode; double refr; double refR; int orders[4]; int maxorder; public: cutProfile(){} ~cutProfile(){} int setorder(std::string); //设置读取列顺序 int setref(std::string); //设置球坐标下的参考半径 int getInputNode(std::string); //输入模型数据点深度最好用负数表示 /*初始化一个规则网络把输入点装进入这里我们有两种生成这个规则网络的方式 1 auto自动生成适用于自由分布的三维空间数据点云 2 para参数生成适用于有一定分布规律的三维空间数据点云所需的参数为三个方向的网络尺寸*/ int initBox(std::string); /*初始化一个平面切面平行于z轴由两个点确定需要指定剖面横纵向剖分间隔*/ int initVecProfile(std::string); /*初始化一个球面上的点*/ int initSphProfile(std::string); /*插值剖面值*/ int interProfile(); /*输出剖面输出剖面二维坐标与实际三维坐标*/ int outProfile(std::string, int); }; int cutProfile::setorder(std::string para) { if (4 != sscanf(para.c_str(),"%d,%d,%d,%d",&orders[0],&orders[1],&orders[2],&orders[3])) { throw gctl::runtime_error("Invalid column indexes."); } else { maxorder = -1; for (int i = 0; i < 4; i++) { if(maxorder < orders[i]) maxorder = orders[i]; } } return 0; } int cutProfile::setref(std::string para) { if (para == "WGS84") { refr = GCTL_WGS84_PoleRadius; refR = GCTL_WGS84_EquatorRadius; } else if (para == "EarthR") { refr = refR = GCTL_Earth_Radius; } else if (2 != sscanf(para.c_str(), "%lf/%lf", &refr, &refR)) { throw gctl::runtime_error("Invalid reference ellipsoid."); } return 0; } int cutProfile::getInputNode(std::string filename) { double temp_d; std::vector tempRow; std::string temp_str; std::stringstream temp_ss; std::ifstream modelin; gctl::open_infile(modelin,filename); node tempNode; while(getline(modelin,temp_str)) { if(*(temp_str.begin()) == '#') continue; temp_ss.clear(); temp_ss << temp_str; if(!tempRow.empty()) tempRow.clear(); while (temp_ss >> temp_d) { tempRow.push_back(temp_d); } if(tempRow.size() < maxorder+1) { throw gctl::runtime_error("Invalid input file format at: "+temp_str); } else { tempNode.x = tempRow.at(orders[0]); tempNode.y = tempRow.at(orders[1]); tempNode.z = tempRow.at(orders[2]); tempNode.val = tempRow.at(orders[3]); tempNode.id = inputNode.size(); inputNode.push_back(tempNode); } } modelin.close(); return 0; } int cutProfile::initBox(std::string para) { //确定点云的范围 xmin = ymin = zmin = 1e+30; xmax = ymax = zmax = -1e+30; for (int i = 0; i < inputNode.size(); i++) { if (inputNode.at(i).x < xmin) xmin = inputNode.at(i).x; if (inputNode.at(i).x > xmax) xmax = inputNode.at(i).x; if (inputNode.at(i).y < ymin) ymin = inputNode.at(i).y; if (inputNode.at(i).y > ymax) ymax = inputNode.at(i).y; if (inputNode.at(i).z < zmin) zmin = inputNode.at(i).z; if (inputNode.at(i).z > zmax) zmax = inputNode.at(i).z; } //确定点云的中心位置 double cx,cy,cz; cx = 0.5*(xmin+xmax); cy = 0.5*(ymin+ymax); cz = 0.5*(zmin+zmax); //若有可用的三方向块体尺寸则使用输入的尺寸 //否则自动计算点云内任意两点间的最短距离作为三方向上的块体尺寸 /* 这里我们发现使用任意两点间的最短距离作为自动网格化大小并不能取得很好的效果经常可能会造成插值点无值的状况因此我们这里将使用平均点位体积估算网格大小效果出众而且这样做仅需要很少的运算注意在处理球坐标点时需要指定三方向网格尺寸不能自动确定 */ double pointSpace; if (3 != sscanf(para.c_str(),"%lf/%lf/%lf",&dx,&dy,&dz)) { pointSpace = (zmax-zmin)*(ymax-ymin)*(xmax-xmin)/inputNode.size(); dx = dy = dz = gctl::sqrtn(pointSpace, 3); } //从中心位置向三个方向推进确定各个方向的块体数 int halfNum_x = ceil((xmax-cx)/dx); int halfNum_y = ceil((ymax-cy)/dy); int halfNum_z = ceil((zmax-cz)/dz); //更新点云范围注意这里我们把盒子向外延伸一层以适应算法的要求 xmin = cx - halfNum_x*dx - dx; xmax = cx + halfNum_x*dx + dx; ymin = cy - halfNum_y*dy - dy; ymax = cy + halfNum_y*dy + dy; zmin = cz - halfNum_z*dz - dz; zmax = cz + halfNum_z*dz + dz; //计算盒子在三方向的个数 xnum = (xmax-xmin)/dx; ynum = (ymax-ymin)/dy; znum = (zmax-zmin)/dz; //初始化boxNode node tempNode; for (int i = 0; i < znum+1; i++) { for (int j = 0; j < ynum+1; j++) { for (int k = 0; k < xnum+1; k++) { tempNode.x = xmin + k*dx; tempNode.y = ymin + j*dy; tempNode.z = zmin + i*dz; tempNode.val = GCTL_BDL_MAX; tempNode.id = boxNode.size(); boxNode.push_back(tempNode); } } } //初始化盒子内点云列表数据将点云索引放入盒子内 int temp_m,temp_n,temp_k; boxIndex.resize(xnum*ynum*znum); for (int i = 0; i < inputNode.size(); i++) { //以块体的左上角位置为准 temp_m = floor((inputNode.at(i).x - xmin)/dx); temp_n = floor((inputNode.at(i).y - ymin)/dy); temp_k = floor((inputNode.at(i).z - zmin)/dz); boxIndex.at(temp_k*xnum*ynum+temp_n*xnum+temp_m).push_back(inputNode.at(i).id); } //插值计算盒子节点上的数据值注意最外层的节点为外扩点所以数据值直接设为BDL_MAX node tempBoxNode; double tempRes,oneDist,totalDist; int cubeIndex; int nodeIndex; for (int i = 0; i < znum-1; i++) { for (int j = 0; j < ynum-1; j++) { for (int k = 0; k < xnum-1; k++) { tempRes = 0; totalDist = 0; nodeIndex = (i+1)*(xnum+1)*(ynum+1)+(j+1)*(xnum+1)+k+1; //查看八个格子内的顶点并利用距离平方反比插值注意这里我们会使用球坐标转换程序只是我们默认的参考半径是0 for (int m = 0; m < 2; m++) { for (int n = 0; n < 2; n++) { for (int p = 0; p < 2; p++) { cubeIndex = (i+m)*xnum*ynum+(j+n)*xnum+k+p; if (!boxIndex.at(cubeIndex).empty()) { for (int b = 0; b < boxIndex.at(cubeIndex).size(); b++) { if (refr == 0 && refR == 0) // 参考球为0 表示输入点在直角坐标系下 { tempNode = inputNode.at(boxIndex.at(cubeIndex).at(b)); tempBoxNode = boxNode.at(nodeIndex); } else { tempNode = c2s_node(inputNode.at(boxIndex.at(cubeIndex).at(b)),refr,refR); tempBoxNode = c2s_node(boxNode.at(nodeIndex),refr,refR); } totalDist += 1.0/(1e-30+pow(tempNode.x - tempBoxNode.x,2)+pow(tempNode.y - tempBoxNode.y,2)+pow(tempNode.z - tempBoxNode.z,2)); } } } } } //如果盒子内没有顶点则初始化boxNode值为BDL_MAX if (totalDist == 0) { boxNode.at(nodeIndex).val = GCTL_BDL_MAX; continue; } //计算插值 for (int m = 0; m < 2; m++) { for (int n = 0; n < 2; n++) { for (int p = 0; p < 2; p++) { cubeIndex = (i+m)*xnum*ynum+(j+n)*xnum+k+p; if (!boxIndex.at(cubeIndex).empty()) { for (int b = 0; b < boxIndex.at(cubeIndex).size(); b++) { if (refr == 0 && refR == 0) // 参考球为0 表示输入点在直角坐标系下 { tempNode = inputNode.at(boxIndex.at(cubeIndex).at(b)); tempBoxNode = boxNode.at(nodeIndex); } else { tempNode = c2s_node(inputNode.at(boxIndex.at(cubeIndex).at(b)),refr,refR); tempBoxNode = c2s_node(boxNode.at(nodeIndex),refr,refR); } oneDist = 1.0/(1e-30+pow(tempNode.x - tempBoxNode.x,2)+pow(tempNode.y - tempBoxNode.y,2)+pow(tempNode.z - tempBoxNode.z,2)); tempRes += inputNode.at(boxIndex.at(cubeIndex).at(b)).val*oneDist/totalDist; } } } } } boxNode.at(nodeIndex).val = tempRes; } } } return 0; } //node n1,node n2,double xinterval,double zinterval int cutProfile::initVecProfile(std::string para) { //注意n2的深度值应该大于n1的深度 double xinterval,zinterval; node n1,n2,snode,enode; if (8 != sscanf(para.c_str(),"%lf/%lf/%lf/%lf/%lf/%lf/%lf/%lf",&n1.x,&n1.y,&n1.z,&n2.x,&n2.y,&n2.z,&xinterval,&zinterval)) { throw gctl::runtime_error("Invalid profile parameter."); } snode = n1; enode = n2; snode.z = n1.zn2.z?n1.z:n2.z; //计算水平和纵向方向剖面长度 double plen = sqrt(pow(enode.x-snode.x,2)+pow(enode.y-snode.y,2)); double pdep = enode.z - snode.z; node tempNode; double depth = snode.z; double length; while(depth <= enode.z) { length = 0; while(length <= plen) { tempNode.x = snode.x+(enode.x-snode.x)*length/plen; tempNode.y = snode.y+(enode.y-snode.y)*length/plen; tempNode.z = depth; tempNode.px = length; tempNode.pz = depth; tempNode.id = profileNode.size(); profileNode.push_back(tempNode); length += xinterval; } depth += zinterval; } return 0; } //lonmin lonmax latmin latmax dlon dlat ele int cutProfile::initSphProfile(std::string para) { node tempNode; double lon,lat; double lonmin,lonmax,latmin,latmax,dlon,dlat,eleva; if (7 != sscanf(para.c_str(),"%lf/%lf/%lf/%lf/%lf/%lf/%lf",&lonmin,&lonmax,&latmin,&latmax,&dlon,&dlat,&eleva)) { throw gctl::runtime_error("Invalid spherical profile parameter."); } lat = latmin; while(lat <= latmax) { lon = lonmin; while (lon <= lonmax) { tempNode.sph_init(lon,lat,eleva); tempNode.id = profileNode.size(); profileNode.push_back(tempNode); lon += dlon; } lat += dlat; } return 0; } int cutProfile::interProfile() { int temp_m,temp_n,temp_k; for (int i = 0; i < profileNode.size(); i++) { //检查待插值点的位置是否在盒子内 if (profileNode.at(i).x < xmin + dx || profileNode.at(i).x > xmax - dx || profileNode.at(i).y < ymin + dy || profileNode.at(i).y > ymax - dy || profileNode.at(i).z < zmin + dz || profileNode.at(i).z > zmax - dz) { profileNode.at(i).val = GCTL_BDL_MAX; continue; } temp_m = floor((profileNode.at(i).x - xmin)/dx); temp_n = floor((profileNode.at(i).y - ymin)/dy); temp_k = floor((profileNode.at(i).z - zmin)/dz); if (refr == 0 && refR == 0) { profileNode.at(i).val = gctl::cube_interpolate(xmin+temp_m*dx, ymin+temp_n*dy, zmin+temp_k*dz, dx,dy,dz, profileNode.at(i).x, profileNode.at(i).y, profileNode.at(i).z, boxNode.at(temp_k*(xnum+1)*(ynum+1)+temp_n*(xnum+1)+temp_m).val, boxNode.at(temp_k*(xnum+1)*(ynum+1)+temp_n*(xnum+1)+temp_m+1).val, boxNode.at(temp_k*(xnum+1)*(ynum+1)+(temp_n+1)*(xnum+1)+temp_m+1).val, boxNode.at(temp_k*(xnum+1)*(ynum+1)+(temp_n+1)*(xnum+1)+temp_m).val, boxNode.at((temp_k+1)*(xnum+1)*(ynum+1)+temp_n*(xnum+1)+temp_m).val, boxNode.at((temp_k+1)*(xnum+1)*(ynum+1)+temp_n*(xnum+1)+temp_m+1).val, boxNode.at((temp_k+1)*(xnum+1)*(ynum+1)+(temp_n+1)*(xnum+1)+temp_m+1).val, boxNode.at((temp_k+1)*(xnum+1)*(ynum+1)+(temp_n+1)*(xnum+1)+temp_m).val); } else { profileNode.at(i).val = interCubeDistSph(refr,refR,xmin+temp_m*dx, ymin+temp_n*dy, zmin+temp_k*dz, dx,dy,dz, profileNode.at(i).x, profileNode.at(i).y, profileNode.at(i).z, boxNode.at(temp_k*(xnum+1)*(ynum+1)+temp_n*(xnum+1)+temp_m).val, boxNode.at(temp_k*(xnum+1)*(ynum+1)+temp_n*(xnum+1)+temp_m+1).val, boxNode.at(temp_k*(xnum+1)*(ynum+1)+(temp_n+1)*(xnum+1)+temp_m+1).val, boxNode.at(temp_k*(xnum+1)*(ynum+1)+(temp_n+1)*(xnum+1)+temp_m).val, boxNode.at((temp_k+1)*(xnum+1)*(ynum+1)+temp_n*(xnum+1)+temp_m).val, boxNode.at((temp_k+1)*(xnum+1)*(ynum+1)+temp_n*(xnum+1)+temp_m+1).val, boxNode.at((temp_k+1)*(xnum+1)*(ynum+1)+(temp_n+1)*(xnum+1)+temp_m+1).val, boxNode.at((temp_k+1)*(xnum+1)*(ynum+1)+(temp_n+1)*(xnum+1)+temp_m).val); } } return 0; } int cutProfile::outProfile(std::string filename, int type) { double temprad; std::ofstream profileOut; gctl::open_outfile(profileOut,filename); if (type == 0) { profileOut << "# ProfileX ProfileZ X Y Z Val" << std::endl; for (int i = 0; i < profileNode.size(); i++) { if (profileNode.at(i).val == GCTL_BDL_MAX) { profileOut << profileNode.at(i).px << " " << profileNode.at(i).pz << " " << profileNode.at(i).x << " " << profileNode.at(i).y << " " << profileNode.at(i).z << " NaN" << std::endl; } else profileOut << profileNode.at(i).px << " " << profileNode.at(i).pz << " " << profileNode.at(i).x << " " << profileNode.at(i).y << " " << profileNode.at(i).z << " " << profileNode.at(i).val << std::endl; } } else if (type == 1) { profileOut << "# ProfileDeg ProfileEle Lon Lat Radius Val" << std::endl; for (int i = 0; i < profileNode.size(); i++) { temprad = gctl::ellipse_radius_2d(refR, refr, profileNode.at(i).y*GCTL_Pi/180.0) + profileNode.at(i).z; if (profileNode.at(i).val == GCTL_BDL_MAX) { profileOut << profileNode.at(i).px << " " << profileNode.at(i).pz << " " << profileNode.at(i).x << " " << profileNode.at(i).y << " " << temprad << " NaN" << std::endl; } else profileOut << profileNode.at(i).px << " " << profileNode.at(i).pz << " " << profileNode.at(i).x << " " << profileNode.at(i).y << " " << temprad << " " << profileNode.at(i).val << std::endl; } } else if (type == 2) { profileOut << "#Lon Lat Ele Radius Val" << std::endl; for (int i = 0; i < profileNode.size(); i++) { temprad = gctl::ellipse_radius_2d(refR, refr, profileNode.at(i).y*GCTL_Pi/180.0) + profileNode.at(i).z; if (profileNode.at(i).val == GCTL_BDL_MAX) { profileOut << profileNode.at(i).x << " " << profileNode.at(i).y << " " << profileNode.at(i).z << " " << temprad << " NaN" << std::endl; } else profileOut << profileNode.at(i).x << " " << profileNode.at(i).y << " " << profileNode.at(i).z << " " << temprad << " " << profileNode.at(i).val << std::endl; } } profileOut.close(); return 0; } #endif