#ifndef _TRACKLINE_H #define _TRACKLINE_H #include "sysDefine.h" class trackLine { public: trackLine(){} ~trackLine(){} int getInputNode(char*); int initBox(char*); int initLine(char*); int interLine(); int outLine(char*); private: node2dArray inputNode; node2dArray lineNode; int xnum,ynum; double dx,dy; double xmin,xmax,ymin,ymax; _2iArray boxIndex; node2dArray boxNode; }; int trackLine::getInputNode(char* para) { char order[1024] = "0,1,2"; char filename[1024] = "NULL"; int orders[3]; double temp_d; _1dArray tempRow; string temp_str; stringstream temp_ss; //首先尝试分离文件名和+d命令 if (2 != sscanf(para,"%[^+]+d%s",filename,order)) { //分离失败 直接将para解释为文件名 strcpy(filename,para); } ifstream datain; if(open_infile(datain,filename)) return -1; if (3 != sscanf(order,"%d,%d,%d",&orders[0],&orders[1],&orders[2])) { cout << BOLDRED << "error ==> " << RESET << "wrong order parameters for the file: " << filename << endl; return -1; } node2d tempNode; while(getline(datain,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); } //读入指定位置的数据 tempNode.x = tempRow.at(orders[0]); tempNode.y = tempRow.at(orders[1]); tempNode.val = tempRow.at(orders[2]); tempNode.id = inputNode.size(); inputNode.push_back(tempNode); } datain.close(); return 0; } int trackLine::initBox(char* para) { //定位二维点位的中心位置 xmin = ymin = BDL_MAX; xmax = ymax = BDL_MIN; 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; } double cx,cy; cx = 0.5*(xmin+xmax); cy = 0.5*(ymin+ymax); //确定盒子的边长 double pointArea; /*这里我们发现使用任意两点间的最短距离作为自动网格化大小并不能取得很好的效果 经常可能会造成插值点无值的状况 因此我们这里将使用平均点位面积估算网格大小 效果出众 而且这样做仅需要很少的运算*/ if (2 != sscanf(para,"%lf/%lf",&dx,&dy)) { pointArea = (xmax-xmin)*(ymax-ymin)/inputNode.size(); dx = dy = sqrt(pointArea); } //确定盒子在x和y方向的半盒子数 向上取整 int halfNum_x = ceil((xmax-cx)/dx); int halfNum_y = ceil((ymax-cy)/dy); //重新确定盒子的范围 向外推出一层 xmin = cx - halfNum_x*dx - dx; xmax = cx + halfNum_x*dx + dx; ymin = cy - halfNum_y*dy - dy; ymax = cy + halfNum_y*dy + dy; //确定盒子在x和y方向的盒子数 xnum = int ((xmax-xmin)/dx); ynum = int ((ymax-ymin)/dy); //初始化盒子上的节点 node2d tempNode; 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.val = BDL_MAX; tempNode.id = boxNode.size(); boxNode.push_back(tempNode); } } //将输入点位放置在盒子中 int temp_m,temp_n; boxIndex.resize(xnum*ynum); 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); boxIndex.at(temp_n*xnum+temp_m).push_back(inputNode.at(i).id); } //查看每个盒子节点周围四个盒子中的输入点 插值确定盒子节点上的值 double tempRes,oneDist,totalDist; int rectIndex; int nodeIndex; for (int j = 0; j < ynum-1; j++) { for (int k = 0; k < xnum-1; k++) { tempRes = 0; totalDist = 0; nodeIndex = (j+1)*(xnum+1)+k+1; //查看四个格子内的顶点 并利用距离平方反比插值 for (int n = 0; n < 2; n++) { for (int p = 0; p < 2; p++) { rectIndex = (j+n)*xnum+k+p; if (!boxIndex.at(rectIndex).empty()) { for (int b = 0; b < boxIndex.at(rectIndex).size(); b++) { totalDist += 1.0/(1e-30 +pow(inputNode.at(boxIndex.at(rectIndex).at(b)).x - boxNode.at(nodeIndex).x,2) +pow(inputNode.at(boxIndex.at(rectIndex).at(b)).y - boxNode.at(nodeIndex).y,2)); } } } } if (totalDist == 0) { boxNode.at(nodeIndex).val = BDL_MAX; continue; } //计算插值 for (int n = 0; n < 2; n++) { for (int p = 0; p < 2; p++) { rectIndex = (j+n)*xnum+k+p; if (!boxIndex.at(rectIndex).empty()) { for (int b = 0; b < boxIndex.at(rectIndex).size(); b++) { oneDist = 1.0/(1e-30 +pow(inputNode.at(boxIndex.at(rectIndex).at(b)).x - boxNode.at(nodeIndex).x,2) +pow(inputNode.at(boxIndex.at(rectIndex).at(b)).y - boxNode.at(nodeIndex).y,2)); tempRes += inputNode.at(boxIndex.at(rectIndex).at(b)).val*oneDist/totalDist; } } } } boxNode.at(nodeIndex).val = tempRes; } } return 0; } int trackLine::initLine(char* para) { char order[1024] = "0,1"; char filename[1024] = "NULL"; int orders[2]; double temp_d; _1dArray tempRow; string temp_str; stringstream temp_ss; node2d tempNode; double length; double xinterval; node2d n1,n2,snode,enode; //首先尝试从参数初始化待插值点 if (5 == sscanf(para,"%lf/%lf/%lf/%lf/%lf",&n1.x,&n1.y,&n2.x,&n2.y,&xinterval)) { snode = n1; enode = n2; //计算水平方向剖面长度 double plen = sqrt(pow(enode.x-snode.x,2)+pow(enode.y-snode.y,2)); 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.px = length; tempNode.id = lineNode.size(); lineNode.push_back(tempNode); length += xinterval; } } //首先尝试分离文件名和+d命令 else if (2 == sscanf(para,"%[^+]+d%s",filename,order)) { ifstream datain; if(open_infile(datain,filename)) return -1; if (2 != sscanf(order,"%d,%d",&orders[0],&orders[1])) { cout << BOLDRED << "error ==> " << RESET << "wrong order parameters for the file: " << filename << endl; return -1; } while(getline(datain,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); } //读入指定位置的数据 tempNode.x = tempRow.at(orders[0]); tempNode.y = tempRow.at(orders[1]); tempNode.id = lineNode.size(); lineNode.push_back(tempNode); } datain.close(); //计算待插值点在剖面上的累积长度 length = 0; lineNode.at(0).px = 0.0; for (int i = 1; i < lineNode.size(); i++) { length += node2dDistance(lineNode.at(i-1),lineNode.at(i)); lineNode.at(i).px = length; } } //分离失败 直接将para解释为文件名 else { strcpy(filename,para); ifstream datain; if(open_infile(datain,filename)) return -1; while(getline(datain,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); } //读入指定位置的数据 tempNode.x = tempRow.at(orders[0]); tempNode.y = tempRow.at(orders[1]); tempNode.id = lineNode.size(); lineNode.push_back(tempNode); } datain.close(); //计算待插值点在剖面上的累积长度 length = 0; lineNode.at(0).px = 0.0; for (int i = 1; i < lineNode.size(); i++) { length += node2dDistance(lineNode.at(i-1),lineNode.at(i)); lineNode.at(i).px = length; } } return 0; } int trackLine::interLine() { int temp_m,temp_n; for (int i = 0; i < lineNode.size(); i++) { //检查待插值点的位置是否在盒子内 if (lineNode.at(i).x < xmin + dx || lineNode.at(i).x > xmax - dx || lineNode.at(i).y < ymin + dy || lineNode.at(i).y > ymax - dy) { lineNode.at(i).val = BDL_MAX; continue; } temp_m = floor((lineNode.at(i).x - xmin)/dx); temp_n = floor((lineNode.at(i).y - ymin)/dy); lineNode.at(i).val = interRectDist(xmin+temp_m*dx,ymin+temp_n*dy,dx,dy, lineNode.at(i).x,lineNode.at(i).y, boxNode.at(temp_n*(xnum+1)+temp_m).val, boxNode.at(temp_n*(xnum+1)+temp_m+1).val, boxNode.at((temp_n+1)*(xnum+1)+temp_m+1).val, boxNode.at((temp_n+1)*(xnum+1)+temp_m).val); } return 0; } int trackLine::outLine(char* filename) { time_t now = time(0); char* dt = ctime(&now); ofstream lineOut; if (open_outfile(lineOut,filename)) return -1; lineOut << "# This file is generated by trackLine on " << dt; lineOut << "# For more information please contact the author via: zhangyi.cugwuhan@gmail.com" << endl; lineOut << "# XonLine XonPlane YonPlane Val" << endl; for (int i = 0; i < lineNode.size(); i++) { if (lineNode.at(i).val == BDL_MAX) { lineOut << lineNode.at(i).px << " " << lineNode.at(i).x << " " << lineNode.at(i).y << " NaN" << endl; } else lineOut << lineNode.at(i).px << " " << lineNode.at(i).x << " " << lineNode.at(i).y << " " << lineNode.at(i).val << endl; } lineOut.close(); return 0; } #endif