由不规则离散点集生成TIN与GRID-DEM

由不规则离散点集生成TIN与GRIDDEM现阶段生成数字高程模型（DEM）的方法较多，如以摄影测量得到的像对为数据源跟踪生成等高线及DEM，由机载激光测距仪记录规则点集后生产数据，也可采用传统的地形图扫描后跟踪等高线，记录一连串离散点集，接着运用各类算法进行处理，最后生成不规则三角网（TIN）与规则格网（GRID）DEM的方法。本文主要介绍的就是以等高线（参考图一）和离散点集为数据源，产生TIN与GRIDDEM的技术路线。具体步骤如下：1）跟踪等高线生成离散点集，记录在文本文件中。参考图二和图三。2）读取文本文件中的数据，进行预处理。主要工作是找到XY轴方向上最小最大数值，压缩数据范围，避免数据范围跨度太大或太小，即出现数据分布稠密或稀疏的情况。while(!_demfile.eof()){_demfilepoint3dXYZ[i][0]point3dXYZ[i][1]point3dXYZ[i][2];point3dXYZ[i][2]=point3dXYZ[i][2]/2;//因为XY轴在随后调整，因此相应调Z轴数值if(xMinpoint3dXYZ[i][0])xMin=point3dXYZ[i][0];//得到整个范围的最大与最小数值if(xMaxpoint3dXYZ[i][0])xMax=point3dXYZ[i][0];if(yMinpoint3dXYZ[i][1])yMin=point3dXYZ[i][1];if(yMaxpoint3dXYZ[i][1])yMax=point3dXYZ[i][1];if(zMinpoint3dXYZ[i][2])zMin=point3dXYZ[i][2];if(zMaxpoint3dXYZ[i][2])zMax=point3dXYZ[i][2];i++;}_demfile.close();//文件流读取完毕，关闭_demfile.clear();//文件流清除3）完成TIN中点，线，三角形和网的定义。classAFX_EXT_CLASSTIN_Point{//TIN中的点friendclassTIN;public:intGet_ID(void){return(m_ID);}//ID数值constPOINT3d&Get_Point(void){return(m_Point);}//普通的点doubleGet_X(void){return(m_Point[0]);}//点上的X数值doubleGet_Y(void){return(m_Point[1]);}//点上的Y数值doubleGet_Z(void){return(m_Point[2]);}intGet_Neighbor_Count(void){return(m_nNeighbors);}//邻接点的个数TIN_Point*Get_Neighbor(intiNeighbor){return(iNeighbor=0&&iNeighborm_nNeighbors?m_Neighbors[iNeighbor]:NULL);}//得到某一邻接点intGet_Triangle_Count(void){return(m_nTriangles);}//得到邻接三角形个数TIN_Triangle*Get_Triangle(intiTriangle){return(iTriangle=0&&iTrianglem_nTriangles?m_Triangles[iTriangle]:NULL);}private:TIN_Point(void);//构造函数TIN_Point(intID,POINT3dpptxyz);virtual~TIN_Point(void);intm_ID,m_nNeighbors,m_nTriangles;//本身ID数值，邻接点个数，邻接三角形个数POINT3dm_Point;//自身存储着XY数值TIN_Point**m_Neighbors;//邻接顶点TIN_Triangle**m_Triangles;//邻接三角形bool_Add_Neighbor(TIN_Point*pNeighbor);//增加顶点bool_Add_Triangle(TIN_Triangle*pTriangle);//增加三角形bool_Del_Relations(void);};classAFX_EXT_CLASSTIN_Edge//TIN中的边{friendclassTIN;//申明友元类public://得到边的端点TIN_Point*Get_Point(intiPoint){return(m_Points[iPoint%2]);}private:TIN_Edge(TIN_Point*a,TIN_Point*b);//构造函数virtual~TIN_Edge(void);TIN_Point*m_Points[2];//边的两个端点};4）以找到预处理后离散点集的最大外包三角形作为队列中第一个三角形开始，依次加入点集中的每个顶点，接着判断此新加顶点与已经存在的每个三角形的外接圆是否存在包含关系，如果包含，则改变此关联三角形的边信息后，然后以新加入的顶点和关联三角形的顶点为参数创建新的三角形，最后收尾工作是，以完成的三角化后每个三角形中的顶点为参数，建立TIN中的顶点、边和三角形之间的关系。boolTIN::_Triangulate(TIN_Point**Points,intnPoints,TTIN_Triangle*Triangles,int&nTriangles){inti,j,k,inside,trimax,nedge=0,emax=200,status=0,*complete=NULL;doublexmid,ymid,dmax,xp,yp,x1,y1,x2,y2,x3,y3,xc,yc,r;TTIN_Edge*edges=NULL;//TTIN_Edge:边的结构//-----------------------------------------------------//Allocatememoryforthecompletenesslist,flagforeachtriangle//trimax=4*nPoints;if((complete=(int*)malloc(trimax*sizeof(int)))==NULL)//给即将生成的三角形的标志分配内存{status=1;if(edges)free(edges);//如果不为空，首先释放内存if(complete)free(complete);}//-----------------------------------------------------//Allocatememoryfortheedgelistif((edges=(TTIN_Edge*)malloc(emax*sizeof(TTIN_Edge)))==NULL)//给边的标志分配内存{status=2;if(edges)free(edges);if(complete)free(complete);}//准备找到最大外包三角形作为三角形队列中第一个三角形,并准备不断的添加新顶点_Extent_Update();//更新整个离散顶点的矩形范围dmax=m_Extent.Get_XRange()m_Extent.Get_YRange()?m_Extent.Get_XRange():m_Extent.Get_YRange();//取范围中XY两轴的较大数值xmid=m_Extent.Get_XCenter();//中心点ymid=m_Extent.Get_YCenter();Points[nPoints+0]-m_Point[0]=xmid-20*dmax;//第一个顶点的X数值Points[nPoints+0]-m_Point[1]=ymid-dmax;//第一个顶点的Y数值，此数值取得相当小Points[nPoints+1]-m_Point[0]=xmid;//第二个顶点的X数值Points[nPoints+1]-m_Point[1]=ymid+20*dmax;//第二个顶点的Y数值，居中，但偏高Points[nPoints+2]-m_Point[0]=xmid+20*dmax;//第三个顶点的X数值Points[nPoints+2]-m_Point[1]=ymid-dmax;//第三个顶点的Y数值，第三个顶点与第一个顶点高度相同，在X轴上对称Triangles[0].p1=nPoints+0;//三角形的顶点记录的是顶点数组中的索引Triangles[0].p2=nPoints+1;Triangles[0].p3=nPoints+2;complete[0]=false;nTriangles=1;//初始化三角形个数为1//-----------------------------------------------------//每次增加一个顶点到已存在的TIN格网中for(i=0;inPoints;i++)//对于所有顶点进行循环{xp=Points[i]-m_Point[0];//使用浮点数进行记录yp=Points[i]-m_Point[1];//使用浮点数进行记录nedge=0;//对于每个顶点，初始化边的个数为0//如果顶点位于三角形的外接圆内，则原三角形的三条边被加到一个缓存中，并移走此三角形for(j=0;jnTriangles;j++)//对于所有刚刚建立的三角形进行循环{if(complete[j]){continue;//此三角形的完整性不用破坏，继续下次循环}x1=Points[Triangles[j].p1]-m_Point[0];//三角形的三个顶点y1=Points[Triangles[j].p1]-m_Point[1];x2=Points[Triangles[j].p2]-m_Point[0];y2=Points[Triangles[j].p2]-m_Point[1];x3=Points[Triangles[j].p3]-m_Point[0];y3=Points[Triangles[j].p3]-m_Point[1];inside=_CircumCircle(xp,yp,x1,y1,x2,y2,x3,y3,&xc,&yc,&r);//判断顶点与外接圆的关系if(xc+rxp)//顶点落在外接圆之外{complete[j]=true;//标记此三角形为完整性}if(inside)//顶点落在外接圆里面{//检查是否超出边链表的大小if(nedge+3=emax){emax+=100;//增加边的个数if((edges=(TTIN_Edge*)realloc(edges,emax*sizeof(TTIN_Edge)))==NULL)//增加新的内存分配{status=3;if(edges)free(edges);if(complete)free(complete);}}edges[nedge+0].p1=Triangles[j].p1;//存储边的信息edges[nedge+0].p2=Triangles[j].p2;//顶点的索引：1,2;2,3;3,1edges[nedge+1].p1=Triangles[j].p2;edges[nedge+1].p2=Triangles[j].p3;edges[nedge+2].p1=Triangles[j].p3;edges[nedge+2].p2=Triangles[j].p1;nedge+=3;//对于最外层循环的顶点来说，相关边的个数加上3Triangles[j]=Triangles[nTriangles-1];/将链表中最后一个三角形赋值过来complete[j]=complete[nTriangles-1];nTriangles--;//三角形个数减去j