意中基于格网的快速等值线充填算法.docx

1、意中基于格网的快速等值线充填算法意中基于格网的快速等值线充填算法 摘要:本文提出一种基于格网的等值线跟踪，适用于任意边界分割的快速充填算法 ，实现充填的矢量化效果。根据边界线与非封闭等值线间的关系，建立等值线间的拓扑 关系，并以树结构方式存储，以准确快速地实现边界线、非封闭等值线之间的封闭和封 闭等值线间的嵌套。该算法已成功应用于海量多波束地形数据的成图，克服了常用多波 束后处理成图软件栅格充填与等值线之间的失配。 关键词：网格数据；等值线跟踪；拓扑结构；填色 A Fast Algorithm of Color Fill between Contours Based on Grid Data

2、WU Zi-yin, GAO Jin-yao (Key Lab of Submarine GeoSciences,SOA,Hangzhou， Zhejiang,310012) Abstract:Based on contour tracing for grid data,a fast algorithm of color fi ll is proposed to cope with random boundaries surrounding data gaps and to r ealize the vector effect of color fill graph. From the rel

3、ationship between boundary lines and unclosed contours,the topologic structure among contours can be built and stored in a tree structure.Such a topologic structure displ ays the linking and closing between boundary lines and unclosed contours,the nesting among closed polygon in the order of tree ri

4、ngs.It is successfully applied in mapping of huge volume of submarine topographic data observed by multibeam,and it overcomes the inconsistency between color blocks and contou rs produced by mapping software for multibeam postprocessing. Key words:grid data;contour tracing;topologic structure;color

5、fill 1引言 等值线充填，就其实质，是等值线间建立拓扑关系13，也就是建立一种树结 构关系，一直是一个比较棘手的问题。尤其在不规则内外边界，任意等值间距的情况下 ，很难建立等值线间拓扑关系。在实际应用过程中，笔者发现有很多软件不能很好解决 这个问题，例如：GMT,SeaView,Simard,Elac等。本文将阐述在规则格网情况下，在边界 线的基础上，建立等值线间拓扑关系，最终实现任意边界，任意等值间距情况下的等值 线充填。 2算法基本原理 等值线跟踪是在已知格网点的基础上，内插出等值线点，然后跟踪等值点，形成封 闭或非封闭等值线47。非封闭等值线的端点必然落在边界线上，而边界线必然是 封

6、闭的，通过跟踪边界线，再把非封闭等值线端点插入边界线，通过边界线建立等值线 间的拓扑关系，在此基础上跟踪出封闭多边形，实现等值线的填充。该算法分为以下几 大步骤：边界线的跟踪，非封闭等值线端点插入边界线中并排序，非封闭等值线建立拓 扑关系，跟踪封闭多边形并排序，等值线的充填。 2.1边界线的跟踪 要实现边界线的跟踪，首先要找出边界线点。所谓边界线点，就是在某一格网点 的周围至少有一空白格网点或非空的边框点，该点即为边界线点。然后用类似跟踪等值 线的方法即可跟踪出边界线（参见图1）。 图1边界线跟踪 Fig.1Tracing boundaries 笔者用三级链表来存储边界线。链表结构如下： ty

7、pedef struct BOUNDARYS /边界线数目链表 BOUNDARY pBoundary;/ 某条边界线指针 BOUNDARYS *p_Next;/ 下一条边界线数目链表指针 BOUNDARYS（） 。 pBoundaryNULL; p_Next=NULL; 。 BOUNDARYS; typedef struct BOUNDARY/边界线链表 BOUNDARYPOINT *pPoint;/某一边界线点指针 BOUNDARY *p_Next;/ 下一边界线点指针 BOUNDARY() pPoint=NULL; p_Next=NULL; BOUNDARY; BOUNDARY; type

8、def struct BOUNDARYPOINT/边界线点链表 CONTOUR *pContour;/某条等值线指针 boolDirection;/ 等值线方向 BOUNDARY *p_Next;/ 下一边界线点指针 BOUNDARYPOINT() pContour=NULL; p_Next=NULL; BOUNDARYPOINT; 3个链表的关系是: BOUNDARYSBOUNDARYBOUNDARYPOINT; BOUNDARYS链表存储所 有的边界线, BOUNDARY链表存储某一条边界线, BOUNDARYPOINT链表存储非封闭等值线端 点。 2.2非封闭等值线端点插入边界线中并排序

9、 非封闭等值线端点必然落在边界线上，并且是有序排列的，在两个边界线点之间 的等值线端点按等值线值要么从高到低排列，要么从低到高排列（参见图2）。 图2非封闭等值线与边界线关系 Fig.2Relationship between unclosed contours and boundaries 笔者把非封闭等值线端点有序地插入BOUNDARYPOINT（边界线点链表）链表中。 2.3非封闭等值线建立拓扑关系 在非封闭等值线端点插入边界线的基础上，非封闭等值线建立拓扑关系。 为了建立非封闭等值线拓扑关系，必须先建立这样一个等值线树结构： typedef struct CONTOUR/某条等值线链表

10、 floatHeight;/等值线值 double*X, *Y;/ 等值点坐标值 longNum;/ 等值线上的坐标点数目 CONTOURLABEL *pLabel;/等值线标注值链表指针 RECTExtRc;/该等值线的最小扩展框; boolClosed;/等值线是否封闭标识符； CONTOUR *pHeadUp, *pHeadDown, *pTrailUp, *pTrailDown;/与该等值线有拓扑关系的等值线指针 bool HUp,HDown,TUp,TDown;/标识等值线的方向 CONTOUR() Height=-1; X=Y=NULL; Num=-1; pLabel=NULL;

11、Closed=false; p_Next=pHeadUp=pHeadDown=pTrailUp=pTrailDown=NULL; HUp=HDown=TUp=TDown=true; HEIGHTLIST*p_Next;/下一等值线指针 CONTOURS*p_NextContour;/ 下一级等值线指针 CONTOUR; 与某条非封闭等值线有拓扑关系的等值线最多有4条，即两端点各指向两条等值线， 这些等值线也可能指向它本身。（图3的填色效果表明该拓扑关系的可靠性） 图3非封闭等值线与边界线建立拓扑关系 Fig.3Topologic relation between unclosed contou

12、rs and boundaries 在“非封闭等值线端点插入边界线中并排序”中已阐述如何把等值线端点插入边界 线中，通过链表BOUNDARYPOINT（边界线点链表），很容易找到与等值线有拓扑关系的等 值线，并把它的指针写入链表CONTOUR（等值线链表）中。 2.4跟踪封闭多边形并排序 从某一条非封闭等值线出发，沿顺时针或逆时针方向跟踪，必然会返回到这条等值 线， 从而构成封闭多边形，然后把这些封闭多边形按扩展框的大小有序地排列，形成封 闭多边形嵌套系列，即完成封闭区域的跟踪。对于封闭等值线，可做为一封闭区域，在 跟踪好封闭多边形后，一起排序，并加入封闭多边形嵌套链表（图4表明该算法完全可以

13、 建立复杂的拓扑、嵌套关系）。 图4建立等值树并实现充填算法 Fig.4Tree structure of contours built for realizing fill algorithm 为了储存跟踪的封闭多边形，笔者建立了一三重链表结构： typedef struct BLOCKNUM/填充区域系列链表 BLOCKS*pBlocks;/某一填充区域系列指针 BLOCKNUM*p_Next;/下一填充区域系列链表 BLOCKS() pBlocks=NULL; p_Next=NULL; BLOCKNUM； typedef struct BLOCKS/填充区域系列链表 BLOCK*pBlo

14、ck;/某一填充区域指针 BLOCKS*p_Next;/下一填充区域系列链表 BLOCKS() pBlock=NULL; p_Next=NULL; typedef struct BLOCK/填充区域链表 PENhPen;/笔型 BRUSHhBrush;/填充图案 CONTOUR*pContour;/等值线指针 boolDirection;/等值线方向:true:正向；false:反向 BLOCK () pContour=NULL; Direction=true; ; 以上3个链表的关系是：BLOCKBLOCKSBLOCKNUM; BLOCK链表用于存储某一封闭多 边形，BLOCKS链表用于存储

15、某一地形单元的封闭多边形系列（例如：海山，凹地），BL OCKNUM链表用于存储BLOCKS。 2.5等值线的充填 在“跟踪封闭多边形并排序”中已详细阐述如何从非封闭等值线跟踪出封闭多边形 ，并储存在BLOCK（填充区域链表）链表中，同时根据BLOCK链表中CONTOUR（等值线链表 ）链表中等值线值，设置封闭多边形的笔型（PEN）、填充图案（BRUSH）、颜色。建立 了非封闭等值线间的拓扑关系（储存在BLOCK中），同时建立了封闭多边形的嵌套关系（ 储存在BLOCKS中）,再根据笔型、填充图案、颜色，很容易实现等值线的充填。 3和其他算法的比较 有很多软件用栅格填充等值线，即在格网的基础上，

16、把格网再细分成若干小格网， 再按格网的值，以矩形色块的方式实现等值线的填充810，例如：GMT, Simard多 波束系统，Elac多波束系统，SeaView等。这种算法比较简单、快捷，但有它本身的缺陷 ，只适合小比例尺，格网很密的情况下，在大比例尺，格网稀的情况下会出现锯齿状多 边形边界（参见图5）。而本算法在任意比例尺，任意格网情况下都可非常好的实现等值 线的充填。 图5锯齿状多边形边界 Fig.5Zigzag polygon boundaries “Win Surfer”是微机上应用很普遍的一种绘图软件，功能强大，速度快捷；笔者 把该算法的计算速度与“Win Surfer”软件对比了一下

17、，速度不亚于“Win Surfer”软 件。笔者发现，在格网很密，图面变化很复杂的情况下，“Win Surfer”软件充填算法 会出现一些问题（参见图6），而本算法适合于任意情况（参见图7）。 图6“WinSurfer”不能实现充填的范例 Fig.6Graphics produced by “WinSurfer” failed in color fill 图7本算法实现的充填图 Fig.7Map of color fill realized by the fill algorithm proposed in this pap er 4算法的优点及应用前景 4.1设计思路新颖、速度快捷 该算法找

18、到了一种快速建立等值树方法。本算法找到了等值线与边界线间的关系， 通过这种关系，很容易，同时也是很快捷的建立等值线间的拓扑关系。该算法不存在理 论上的缺陷，同时速度不亚于同类商用软件。 4.2适用任意情况 在图面变化复杂，格网密，内外边界多的情况下，有很多软件无法进行等值线充填 ，或者不能很好地进行充填，而本算法适用于任意情况。 4.3应用前景广阔 本算法不但很好地解决了等值线充填问题，而且在GIS应用上有很广阔的应用前景 。该算法准确、快速建立不同等值区域间的拓扑关系，可应用于精确统计特定区域的周 长、面积、体积等。 该算法也可适用于三角网情况下的等值线充填。该算法是在“863”计划海洋领域

19、有 关多波束地形探测技术开发中研制的，它的应用也体现了我们自主开发的多波束后处理 成图软件的特色。 致谢 以上研究是“863”计划海洋领域“820-01-01”课题的研究成果之一。参与该课 题的各位同仁给予了各种帮助，特别感谢金翔龙院士在课题研究过程中给予的具体指导 ，也感谢“126”专项课题的同仁在算法应用、检验中给予的方便。 863计划海洋领域820-01-01课题资助。 作者简介：吴自银，男，28岁，助理研究员。主要从事多波束后处理软件编制GIS研究工 作。 作者单位：国家海洋局海底科学重点实验室,浙江杭州,310012 参考文献： 1王伟，等. GeoStar中图形编辑与拓扑关系的建立

20、J. 武汉测绘科技大学学报 ， 1995，20(5). 2袁修孝，龚健雅. 顾及地形特征线的数字高程模型软件包J. 武汉测绘科技 大学学报，1995，20(5). 3孙玉国. 拓扑空间关系描述与2DT-String空间关系表达D.武汉：武汉测绘科 技大学，1993. 4龚健雅. 顾及地形特征的DEM内插与等高线绘图子系统J. 测绘学报，1990， 15(1). 5 荨 龚健雅. 关于DEM中多面函数内插法几个问题的研究J. 测绘通报，1995(5) . 6龚健雅,苏向辰. 一种快速内插数字地面模型的方法J.测绘通报，1987(2). 7 荨 龚健雅. GIS中矢量栅格一体化数据结构与面向目标数

21、据模型的研究D.武汉 ：武汉测绘科技大学，1992. 8毋河海. 地图数据库系统M. 北京：测绘出版社，1991. 9李德仁，龚健雅，边馥苓. 地理信息系统导论M. 北京：测绘出版社,1993. 10 荨 龚健雅. 整体SIS的数据组织与处理方法M. 武汉：武汉测绘科技大学出版 社，1993.等值线的栅格图形自动生成 王 忠，刘浩吾，柴贺军 (四川大学 水电工程学院) 摘 要：本文针对常规的绘制等值线的方法，提出了一种新的方法，该方法基于计算机图形是由离散点组成，同一数值的所有象素点即构成等值线的原理。由于不需要网格搜索和等值线的游动，使得编制程序简单，绘制的图形精度高，更具有一般性。文中给出

22、了具体的算法，同时利用该方法绘制了某水电工程坝区地形等值线图。 关键词：等值线；程序设计；算法 收稿日期：2000-03-20 作者简介：王忠(1971-)，男，重庆市万县人，现在四川大学攻读岩土工程的博士学位，主要从事岩土工程和计算机应用的研究和开发工作。 等值线是一种形和数的统一，在水利、土木、地质等工程和技术领域内广泛的应用，特别是有限元计算的等值线分析中，它是成果数据表示的重要手段之一。所谓绘制等值线就是对大量离散的、又具有一定规律的几何量值或物理量值，用数学的方法插值并将具有相同量值的点变换成图的过程。常规的等值线绘制通常采用网格法，其绘制的步骤一般为：离散数据网格化；网格点数值化；

23、等值点的计算；等值线的追踪；光滑和标记等值线1。一般来说，这种方法中将要用到的网格是四边形网格2。最近，文献3针对四边形网格的一些问题提出了引入三角形网格的方法。这两种方法的共同之处都是要利用网格及网格上的等值点进行游动追踪，从而使得绘制过程具有如下的缺陷：(1)绘制过程复杂，程序的实现具有一定的难度。这两种方法利用网格结构，首先找出某个四边形网格或三角形网格各边上的等值点，然后对所有网格进行循环，查找等值点，其中涉及许多判断，增加了程序实现的难度。当网格的节点成为某等值点时，必须作为奇点处理，这种情况的处理既降低了图形的精度又增加了绘制的复杂度。(2)绘制的图形的精度不够，并可能在游动追踪时

24、出现相交的情况。上述方法对非网格上的点的处理是用某种曲线拟合的方法。也就是说，该方法进行了两次近似，误差较大。(3)方法不具有一般性，它只能处理具有网格化结构的数据，如果要将某数据转化为网格结构，其中将又要插值，这必将降低图形的精度。针对上述问题，本文提出了绘制等值线的栅格图形法，这里称它为无网格法。该方法不需要网格化的结构，与常规方法相比具有如下的优点：(1)编程简单，容易实现，其中不涉及奇点，也没有等值线的游动追踪，降低了程序设计的复杂性；(2)精度更高，只需一次近似，常规方法则需两次或更多次近似；(3)更具有一般性，没有网格的限制。 1 栅格图形的等值线生成法 栅格图形的绘制有以下特点4

25、：要绘制等值线的区域是有限的，并且是由不连续的点构成的。事实上，栅格图形受到计算机屏幕的限制，人们能够看到的只是几千个或上万个离散的象素组成的图形。比如直线，在计算机上是有限长的，同时是用一个一个的离散的点来表示的，由于人眼的分辨率的局限性，使得它看起来就是连续的。基于上述特点，本文提出了栅格图形的等值线生成法。该方法的基本思想是：利用计算机图形是由离散点组成的特点，只需找出同一等值线上的所有象素点，这些点必然构成该等值线。 这里以矩形山区作等高线为例，来谈谈具体的算法。计算中需要的数据是各个测点的坐标和海拔高度，即(xi,yi,zi)，其中zi代表第i个测点海拔高度，共有M个测点。同时给出将

26、要绘制的等高线的高度，比如从h0开始，每相差h0的高度就绘制一条等高线，共绘制m条等高线。另外已知将要显示的屏幕区域的大小，这里用(StartX,StartY)表示该区域的左上角，(EndX，EndY)表示该区域的右下角。绘制其等高线的算法如下： (1)找出绘制区域左上角和右下角点的xi、yi值，分别用Xmax，Xmin，Ymax,Ymin表示； (2)将坐标(xi,yi)变换成屏幕坐标(SXi,SYi).所需的变换公式为： sxi=xi-Xmin/Xmax-Xmin(EndX-StartX)； syi=yi-Ymin/Ymax-Ymin(EndY-StartY) (1) (3)令istart

27、X,j=StartY; (4)利用求高度的算法(例如距离加权法，最小二乘法等)求出点(i,j)、(i+1,j)和(i,j+1)的高度h1、h2、h3，即图1中的P1、P2和P33个点的高度； (5)检查h1、h2、h3的值，按如下方法确定是否有等高线通过：k=1,h=h0；判断(P1-h)*(P2-h)0是否成立。如果成立，就执行下一步，否则执行；判断P1-h=P2-h是否成立。如果成立，表明有一条等高线通过P1、P2，在这两点打点，并跳到(6)，否则执行下一步；判断P1-hP2-h是否成立。 图1 高程计算示意 如果成立，表明有一条等高线通过P1，就在该处打点，否则就在P2处打点；判断(P1

28、-h)*(P3-h)0是否成立。如果成立，就执行下一步，否则执行；判断P1-hP3-h是否成立。如果成立，在P1,P3这两点打点，并跳到(6)，否则就跳到；判断P1-hP3-h是否成立。如果成立，就在P1处打点，否则就在P3处打点；令k:=k+1，判断km+1是否成立。如果不成立，就执行下一步，否则令h:=h+h0，返回。 (6)令j:j+1，判断j(7)令i:=i+1，判断i(8)结束。 在具体的程序设计中，为了避免重复计算，可以用一个数组来保存所有i+1列P2点的值，同时用另外一个变量来保存P3的值。 从上面的算法可以看出，该方法没有涉及等值线的游动、网格奇点的判断和等值线的连接等问题，从

29、而大大简化了程序的编制，很容易实现，且所绘制的图形不会出现相交线情况。 2 网格化和节点的确定 一个算法的计算时间消耗是人们关注的问题。因*.舅惴谇笕我坏*(i,j)的高度时，需要找出全部对该点的高度有贡献的点。若对整个数组进行搜索，时间消费将是非常之大。为此，引入下面的规则化网格方法来加快速度。首先给出两个概念，即影响域和影响点集。定义如下： 定义1：节点P的影响域O(P)是指该点能够对其它点有所影响的最大的区域。本文中可以指半径为r(预先给定)的闭圆盘，或者为边长为a(预先给定)的正方形。 定义2：节点P的影响点集S(P)是指能够对点P有所影响的点的集合。本文中指元素个数为n(预先给定)的点集，也就是说，能够对点(i,j)的高度有所贡献的已知点的个数只能为n个，这些点一般为到点P最近的n个点。 图2 规则网格化根据上面的定义，在求任一点(i,j)的高度时，只需要找出对该点高度有影响的所有节点，然后根据二维平面拟合的内插算法即