栅格地图矢量化是地理信息系统数据采集的重要方式，可分为自动矢量化和手动矢量化两大类[1, 2, 3]。线状要素是地图中最重要的组成部分,而面状要素用图者最关心的是其外围线[4]，利用外围线可以方便地生成面状要素[5, 6]。而点状符号及文字则不适合自动矢量化[7, 8, 9, 10]。自动矢量化方法目前已有多种，可分为基于细化算法的方法和非细化的方法两类[9, 10, 11, 12, 13, 14]。2005年前《测绘学报》上发表过较多关于自动矢量化的文章，在一定程度上推动了自动矢量化理论的发展[13, 15]，近年来相关研究主要提出了非细化的方法、Zhang快速细化算法、矢量化质量控制理论等[11, 14]。采用步长搜索框的非细化算法在强调矢量化速度的情况下矢量化精度通常不能满足要求，基于“边缘夹比”的非细化方法则费时较长且极易受栅格线体边缘的影响，使矢量线体呈锯齿状，不符合原栅格地图线划所表示的位置和形状；基于细化算法的方法则先对栅格线体进行细化求得骨架线，然后自动追踪骨架线生成矢量线，再对矢量线进行压缩，该过程费时较长但精度可取。传统的细化算法计算速度慢、易受线体边缘影响。近年来出现的Zhang快速细化算法的改进算法较为优秀，其细化结果受线体边缘影响小[16]。但这些算法的细化结果仅包括背景像素及线体像素两种，使得自动追踪矢量化过程中将像素行列号转换为矢量坐标以及矢量线节点压缩计算量很大，并且压缩后的矢量线可能偏离栅格线划。为此，本文设计基于扫描细化算法的自动矢量化方法，该方法速度较快、精度较高、在交点处不相互打断、后期处理量较小。
1 扫描细化算法的原理扫描细化算法是根据栅格线状符号的形态特征及导数法则设计的直接搜寻特征像素的算法，该算法将线体细化像素分为普通细化像素(标记为2)和特征像素两类，特征像素又分为端点(标记为5)、节点(标记为3)和交点(标记为4)3种，这加快细化计算速度，并减少自动追踪矢量化时坐标转换和节点压缩。
1.1 栅格线体的形态特征及导数原则模拟地图经扫描得到的栅格图像中所有线划符号均是以相邻的线体像素组合而成的，这些像素具有一定的规律。如图 1(a)，为使地图中线划符号的拐弯能被分辨出来，至少在拐弯处的两内边界上留下1倍线宽为背景像素，这样该拐弯处的中心位置的横跨像素个数将达到3倍线宽，那么在保证2倍线宽的情况下，线体可继续保留该方向序贯的趋势。
|
| 图 1 栅格线划像素的形态特征 Fig. 1 Morphology of the line in raster image |
当确定两个邻近的特征像素后，线体的走向仍然有沿着该两个特征像素所指方向的推进趋势，这就是导数原则。如图 1(b)，圆圈表示刚获取的两个特征像素，虚线表示继续前推的趋势。
1.2 算法原理通过已获得的特征像素按导数法则推进，如果跳出线体边界或推进长度达到2倍线宽则通过垂直于推进方向的方向进行截断来获取该处的中心像素作为特征像素，然后将此特征像素与上一特征像素利用直线反连，如此往复即可细化线体。
1.3 搜寻起始特征像素(1) 按图 2(a)进行推扫，记录线体该处的横向宽度Ln1,并得出该处的横向中心像素(圆圈处)，以及跳出线体的位置为STXHP。
(2) 按图 2(b)在该中心像素处向上及向下推扫，获取该处纵向的上行像素数a和下行像素数b,纵向宽度Ln2=a+b,并获取纵向上的中心像素。
(3) 按图 2(c)获得该像素左右的像素数a和b，该处的横向宽度Ln3=a+b，并求得该处的横向中心像素(菱形处)。
判断：如果Ln3≤1.5Ln2，将第(1)步确定的中心像素标记为5及STP，将第(3)步确定的中心像素标记为3及SCDP，拟定线宽Ln=Ln3；否则将第(2)步确定的中心像素标记为5及STP，拟定线宽Ln=Ln2。在STP处横向向右偏移0.5Ln记录为TP(待定点),如图 2(e)，获取该方向上的中心像素(三角形处),标记为3及SCDP，之后将STP与SCDP之间的像素标记为2。
凡此过程中遇到已标记为细化像素则直接搜寻下一线体的起始像素。
|
| 图 2 搜寻起始特征像素 Fig. 2 Seeking the start feature pixel and end feature pixel |
如图 3(a)已知上一特征像素BP(圆圈处)及当前特征像素NP(菱形处)，按照BP—NP方向前推，直到跳出线体或递推达到2倍线宽时止，标记该像素为TP(图 3(b)中三角形)。之后在TP处按照截断方向(JDL)推扫获取中心像素(图 3(b)中圆形)。并将BP置为NP，NP置为刚获取的中心特征像素；如果在递推过程中发现递推像素已细化则进行交点判断处理。之后将NP和BP之间的像素标记为2(此步称为反连)。
|
| 图 3 待定点及中心像素的获取 Fig. 3 Seeking the undetermined pixel and the center pixel |
(1) 前推方式可以用BP—NP的近似导数来确定，设纵向像素差及横向像素差分别为dx及dy，近似导数P按式(1)计算。若P>1则令Hya=1，Zya=P前推时先按纵向推进，至Zya时横移1像素并纵移1像素，偏移方向与dx及dy符号相同，如此往复; 若P=1则令Hya=Zya=1，每次横向偏移1像素，同时纵向偏移1像素;若P<1则令Zya=1,Hya=P,先横向偏移，至Hya时再纵向偏移1像素，如此往复。
说明：由于预处理后线体的平均宽度小于10像素，所以2倍平均线宽小于20像素，那么设定前推限制为30像素时，可使推进到二倍线宽时不会产生偏移。
(2) 反连方式可采用计算机图形学中绘制直线算法，如DDA算法。
(3) 截断方向按式(2)计算。如图 4(b)表示截断方向JDL的实际方向
式中，“&&”表示“逻辑与”运算。
|
| 图 4 截断方向以及线尾判断 Fig. 4 Scanning direction and the judgment of the end of line |
线体本端细化完毕的第1种条件是相交，第2种条件是悬挂线尾。如图 4(c)，前推时若前方像素是已细化像素，则遇到了相交，此时在该细化像素周围的0.5Ln范围内搜寻特征像素，若搜寻到则将该特征像素标记为4；如果没有搜寻到特征像素，则直接将之标记为4，再反连细化。
第2种条件如图 4(d)，前推长度小于0.3Ln，并且其左右推扫像素数a+b<2Ln。此时将NP标记为5(线头线尾像素)，之后反连细化；如果a+b≥2Ln,则从NP处向a与b中较大的一方偏移0.5Ln再继续前推细化。
2 扫描细化和自动矢量化算法 2.1 扫描细化算法(1) 从上到下从左到右搜寻线体的起始特征像素STP、SCDP以及线宽Ln，置XH=1(XH=3时跳出本线的细化进程)，令BP为STP，NP为SCDP。
(2) 依据NP和BP计算前推方式及JDL，从NP反连至BP。
(3) 在NP处前推判断是否为交叉点及线结尾，如果既非结尾亦非交叉点，则继续前推直到获取到截断像素TP，然后进行第(4)步。如果是结尾交叉则判断是否为闭合线结尾，若是则标记结尾点为5然后执行第(7)步；如果是悬挂线结尾则将结尾点标记为5，XH++，然后执行第(6)步；如果是一般交叉点，则将之标记为4，XH++，然后执行第(6)步。
(4) 在TP处根据JDL获取中心像素，并将之标记为3，并记录为NP，将原NP记录为BP。
(5) 重复执行第(2)步至第(4)步直到本线体扫描完毕，执行第(7)步。
(6) 如果XH=3则执行第(7)步，如果XH=2则将STP记录为NP,并置其为3，再将SCDP标记为BP,计算推扫方式及截断方式，然后执行第(3)步。
(7) 从STXHP开始从左至右搜索线体的起始特征像素并将第一象特征元标记为BP和STP，第二特征像素标记为NP和SCDP,获取线宽Ln，置XH=1，获取新的STXHP，然后执行第(2)步，如果搜索到最后一行最后一列仍然没有进入线体，则栅格地图细化完毕。
2.2 自动追踪矢量化算法(1) 搜寻端点像素，如果已搜寻到图像的最后一像素则执行第(4)步；否则将之标记为NP，线号增1，将该像素转为矢量点，将点号存入当前线中。
(2) 从NP处搜寻8邻域中的非BP的下一细化像素，标记为NP，将原NP标为BP。
(3) 若刚搜寻的像素是非端点的特征像素则将该像素转为矢量点，将点号存入当前线中，然后重复第(2)步，若为端点像素则将该像素转为矢量点，将点号存入当前线中，然后执行第(1)步，若为交点，则作为普特征点处理，将交点存入交点表，其下一像素的搜索按从上一特征像素指向交点的方向为起始搜索，保证在交点处不断开。
(4) 从交点表中调出交点作为起始点，追踪方法同上，仅一端追踪完毕后还要从本交点反向追踪，所有交点追踪完毕后输出矢量数据文件。
此过程中凡被标为NP的像素除交点像素外均置为0,且在判断下一细化像素时，从BP—NP的方向开始同时沿顺时针及逆时针逐像素判断。
3 试验及分析 3.1 试 验利用C#语言编写程序实现该算法，采用该程序对不同类型的扫描地图进行自动矢量化后统计其位置精度。图 5为扫描细化算法的试验(局部)；图 6为对彩色线划道路网进行色彩提取自动矢量化及后处理的过程。之后通过对4幅不同等高线疏密程度的扫描地形图按每天工作8 h采取如下4种方法来矢量化：扫描细化法自动矢量化、MapGIS半自动矢量化、MapGIS全自动矢量化及精细的手动矢量化，分别统计出用时量及位置误差。
|
| 图 5 扫描细化结果 Fig. 5 The result of scan-thinning |
|
| 图 6 彩色区划地图矢量化 Fig. 6 Auto-vectorization of the colored-area-graphics |
图 5试验中分辨率为300dpi，则1像素相当于纸质地图上0.085mm，统计如表 1，可见细化结果符合《数字线划地形图、数字高程模型质量要求》(GB/T 1794.1—2000)的精度要求。
| 点类 | 最小误差 | 最大误差 | 中误差 | 允许误差/mm | |||
| /像素 | /mm | /像素 | /mm | /像素 | /mm | ||
| 特征点 | 0 | 0 | 0.75 | 0.064 | 0.3 | 0.026 | 0.2 |
| 普通细化点 | 0.1 | 0.008 5 | 2 | 0.17 | 1.2 | 0.102 | 0.2 |
通过图 5(b)、(c)两图可直观地看出不进行光滑时矢量线大致居于栅格线体的中心，光滑后整个线体更加趋于栅格线体中心。 3.2.2 时间复杂度分析
同传统算法相同，扫描细化算法具有O(N)的起始追踪点寻找时间。令K为线体像素数(远小于N)，T为平均线宽，那么线体中心像素数则为K/T。由于每向前搜索2倍线宽就进行横向截断扫描1倍线宽，那么需要O(((2+1)/2)K/T)≈O(1.5K/T)搜寻特征像素的时间以及O(K/T)的反连时间，所以扫描细化算法的总时间为O(N+2.5K/T)；传统算法中最快的边缘跟踪剥皮算法所需总时间为O(N+ K)，所以，当T≥3pixel时扫描细化算法效率高于边缘跟踪剥皮算法；而Zhang快速细化算法需利用4个3×3模板进行n轮比对，直到没有可被删除的像素为止，然后进行一轮冗余细化像素的剔除，所以总时间为O((4n+1)N),可见该算法效率并不高。
针对扫描细化算法的自动追踪矢量化算法坐标转换计算量及冗余坐标点压缩计算量远小于传统算法的计算量。
对试验中采用4种不同方法进行矢量化的用时及误差进行统计如表 2。 可见，全手工方法用时约为扫描细化算法的8倍，MapGIS半自动矢量化方法则约为4倍，MapGIS全自动矢量化及后处理方法为9倍。从精度上说，MapGIS的半自动及全自动矢量化均不及扫描细化法及手工精细矢量化法。
| 图幅 | 矢量化时间/h | 平均线体偏移(纸图)/mm | ||||||
| SC | MAN | MA | MN | SC | MAN | MA | MN | |
| Ⅰ | 7 | 40 | 90 | 84 | 0.05 | 0.19 | 0.18 | 0.05 |
| Ⅱ | 12 | 52 | 104 | 107 | 0.08 | 0.16 | 0.19 | 0.07 |
| Ⅲ | 18 | 74 | 131 | 120 | 0.1 | 0.18 | 0.18 | 0.12 |
| Ⅳ | 23 | 80 | 150 | 146 | 0.07 | 0.18 | 0.18 | 0.09 |
| 注：SC为扫描细化法;MAN为MapGIS半自动矢量化;MA为MapGIS全自动矢量化;MN为手工法。 | ||||||||
可见扫描细化算法矢量化精度高、效率高，可大大减少矢量化工作的成本，提高矢量化成果的质量。并且能直接得出栅格线划符号的中心特征像素，并将之单独标识，缩减自动矢量化时追踪线体的坐标转换时间，且自动矢量化出来的线体较少出现锯齿状波动，减少自动矢量化的后处理时间。
| [1] | GONG Jianya. The Basic of Geographic Information System[M]. Beijing: Science Press, 2001:33-41.(龚健雅.地理信息系统基础[M].北京:科学出版社,2001:33-41.) |
| [2] | ZHANG Donghui, HE Zhengwei, YANG Bin. Study and Implementation of Raster Image Automatic Vector System[J].Computer Engineering and Applications, 2010,46(10):171-174.(张东辉,何正伟,杨斌.栅格图像自动矢量化系统的研究与实现[J].计算机工程与应用,2010,46(10):171-174.) |
| [3] | KHOTANZAD A, ZINK E. Contour Line and Geographic Feature Extraction from USGS Color Topographical Paper Maps[J].PAMI, 2003, 25(1):18-31. |
| [4] | ZHAI Renjian, WU Fang, ZHU Li, et al. Structured Representation of Curve Shape[J].Acta Geodaetica et Cartographica Sinica,2009,38(2):175-182.(翟仁健,武芳,朱丽,等.曲线形态的结构化表达[J].测绘学报,2009,38(2):175-182.) |
| [5] | DENG Min, MA Hangying. The Hierarchical Representation of Topological Relations between a Line and an Area[J]. Acta Geodaetica et Cartographica Sinica, 2008,37(4):507-520.(邓敏,马杭英.线与面目标拓扑关系的层次表达方法[J].测绘学报,2008,37(4):507-520.) |
| [6] | GUO Qingsheng, DU Xiaochu, LIU Hao. Research on Quantitative Representation and Abstraction of Topological Relation between Two Regions[J] . Acta Geodaetica et Cartographica Sinica, 2005, 34 (2): 123-128. (郭庆胜, 杜晓初, 刘浩. 空间拓扑关系定量描述与抽象方法研究[J]. 测绘学报, 2005, 34 (2): 123-128. ) |
| [7] | ZHANG Xin, CHEN Guoxiong, ZHONG Ershun. Line Object Pick-up Based on Improved Grid Thinning Algorithm[J].Geo-information Science, 2007,9(3):25-27.(张欣,陈国雄,钟耳顺.优化栅格细化算法的现状地物提取[J].地球信息科学,2007,9(3):25-27.) |
| [8] | WANG Yong, LI Chaokui. An Excellent Vectorization Algorithm and Its Programming for Raster Data of Linear Feature[J].Engineering of Surveying and Mapping, 2009,18(5):44-50.(王勇,李朝奎.线性要素栅格数据矢量化快速简便算法设计与实现[J].测绘工程,2009,18(5):44-50.) |
| [9] | LI Yanying, GUO Jingjun. Research on Integral Vectorization Method for the Scanned Image of Large-scale Topomap[J]. Bulletin of Surveying and Mapping, 2000(8):20-21. (李岩影, 过静珺.大比例尺地形图扫描图像的整体矢量化方法研究[J].测绘通报, 2000(8):20-21.) |
| [10] | CHEN Zhengguang, WU Yushu, WANG Yufang. A New Vectorization Method of Contour[J].Computer Engineering and Applications,2004(3):84-86. (陈争光,吴裕树,王玉芳.一种新型的地形等高线线性体矢量化方法[J].计算机工程与应用,2004(3):84-86.) |
| [11] | FU Qinghua, NI Shaoxiang, GUO Jian, et al. Vectorization of Raster Data and Solution to Relevant Problems[J]. Geo-Information Science,2004, 22(2):85-89. (扶卿华,倪绍祥,郭剑,等.栅格数据矢量化及其相关问题的解决方法[J].地球信息科学,2004, 22(2):85-89.) |
| [12] | FU Yong, QIN Zhiyuan, GONG Jingjing.Research Actuality and Thoughts on the Automatic Extraction of Linear Objects from Remote Sensing Images[J]. Journal of Geodesy and Geodynamics, 2002(2):112-116.(傅勇, 秦志远,龚菁菁.遥感影像上线状目标自动提取研究现状与思考[J].大地测量与地球动力学,2002(2):112-116.) |
| [13] | WU Huayi, GONG Jianya, LI Deren. Non-boundary Run-length Encoding System for Raster and Its Relevant Algorithms[J]. Acta Geodaetica et Cartographica Sinica,1998,27(1): 63-68.(吴华意, 龚健雅, 李德仁.无边界游程编码及其矢栅直接相互转换算法[J].测绘学报,1998,27(1): 63-68.) |
| [14] | JIA Hongling, WANG Zuoyong, CUI Haiyan, et al. On the Quality Control and Evaluation Method of Scanned Vector Data[J]. Bulletin of Surveying and Mapping,2006(4):69-71. (贾红玲,王作勇,崔海燕,等. 扫描矢量化的数据质量控制与评价方法浅析[J].测绘通报,2006(4):69-71.) |
| [15] | XIE Shunping,DU Jinkang, WANG Lachun et al. Approach of Vectorization for GIS Raster Data Based on Run-length Encoding System[J]. Acta Geodaetica et Cartographica Sinica, 2004,33 (4 ): 323-327. (谢顺平,都金康,王腊春,等.基于游程编码的GIS栅格数据矢量化方法[J].测绘学报,2004,33 (4 ): 323-327.) |
| [16] | YANG Ronghao, YU Daijun. Solutions about Some Minutial Problems in Automatic Vectorizing of Scanned Map[J]. Surveying and Mapping of Sichuan, 2008,31(3):127-130.(杨容浩,余代俊.扫描地图自动矢量化若干细节问题的解决方法[J].四川测绘,2008,31(3):127-130.) |
| [17] | YANG Yun, ZHU Changqing, SUN Qun. Vectorization of Linear Features on Color Scanned Topographic Maps[J]. Journal of Computer-aided Design and Computer Graphics, 2009,21(4):533-541.(杨云,朱长青,孙群.彩色扫描地图中线目标的矢量化方法[J]. 计算机辅助设计与图形学学报,2009,21(4):533-541.) |