矢量和栅格是两种最基本的空间数据模型，在应用中各具优势。依据一定准则，将矢量数据离散化到对应尺度格网上的过程称为离散化或格网化，这是矢量数据与栅格数据同构处理的重要环节。矢量数据模型将现实世界的实体抽象为点、线、面几何要素的组合^[1]，因此几何要素的离散化是矢量离散化的基础。其中，线和面要素的离散化首先需要确定线的格网路径或面的格网边界，核心问题都是线的离散化，即离散线生成^[2-3]。

在计算机图形学及地理信息系统领域，已研发多种离散线生成算法，具有代表性的有数字微分法^[4]、中点画线法^[5]和Bresenham算法^[6]。Wu等^[7]基于Bresenham算法提出了双步直线生成算法；在此基础上，贾银亮等^[8]提出了六步直线生成算法；黄斌茂等^[9]提出了基于自适应步长的直线生成算法。这些算法多是针对矩形格网提出，而在很多应用中，除了矩形，还有三角形、六边形格网，适用于后两类格网的离散线生成算法并不多见。尽管Vince^[10]建立了多维格网的离散线数学模型并给出相应算法，但其仅适用于可由单一泰森多边形（Voronoi）单元（如正方形、六边形、立方体等）平移得到的格网系统。对于在地形建模、空间分析、模拟仿真中广泛应用的三角形格网，该算法并不适用。单元非一致相邻且方向不一致的特点导致三角形格网上的离散线问题比较复杂^[11]，相关研究较少。Freeman^[12]提出了基于顶点的三角形格网离散线生成算法，该算法根据单元顶点到矢量线的距离筛选单元，导致部分离散线单元偏离矢量线。Nagy^[13]提出了基于相邻序列的三角形格网离散线生成算法，该算法仅以格网单元的相邻关系及路径最短为约束条件，无法确保路径唯一。Sarkar等^[14]提出了沿三角形单元边界的最短路径搜寻算法，但该算法忽略了单元中心与矢量线的关系。

本文根据三角形格网的结构特点，利用格网间的弱对偶关系，借助六边形格网离散线等效建立了三角形格网离散线数学模型，并降维求解，尝试从理论层面简化问题；并设计实现了平面三角形格网离散线转化生成算法。对比实验证明了本文算法的效果更优。

1 对偶及弱对偶定义

在数学中，对偶定义为存在于概念、理论及数学结构之间的一对一的转化关系^[15]。格网的对偶是指不同格网之间特有的位置关系，连接相邻格网单元的中心并将其作为顶点，可得到相应的对偶格网^[16]。三角形格网与六边形格网存在对偶关系，设三角形格网单元的边长为l_t，则相应对偶六边形格网单元的边长$ {l_h} = \frac{1}{{\sqrt 3 }}{l_t}$，且两者相互垂直，如图 1所示。

作为格网单元计算、定位的参考点^[17]，通常以单元中心代替格网单元。因此，为了简化格网转换算法，应尽量使三角形格网与六边形格网单元中心重合。显然，图 1所示的对偶关系无法满足该需求。将六边形单元边长调整为${l_h}{\rm{'}} = \frac{1}{3}{l_t}$，且平行于三角形单元的边，则中心不在三角形顶点处的六边形单元与三角形单元中心重合，如图 2所示。本文称该关系为弱对偶，可确保两类格网单元中心最大限度重合，便于格网转换。

2 三角形格网离散线模型

基于三角形格网与六边形格网之间的弱对偶关系，可实现六边形格网与三角形格网之间离散线模型的转化，简化三角形格网离散线模型的求解过程。

2.1 六边形格网离散线模型

在六边形格网中，由某一单元中心指向其相邻单元中心的一组向量称为方向向量。非共线的两个方向向量线性无关，则由起点单元到终点单元的离散线可由任意两个非共线方向向量的有序排列表示。在实际应用中，选择与矢量线ab夹角最小的两个方向向量生成离散线，称为最优方向向量，记作$V' = \left\{ {\mathit{\boldsymbol{v}}_0^{'}, \mathit{\boldsymbol{v}}_1^{'}} \right\}$，如图 3所示。

令$ W' = \left\{ {\mathit{\boldsymbol{w}}_1^{'}, \mathit{\boldsymbol{w}}_2^{'} \cdots \mathit{\boldsymbol{w}}_N^{'}} \right\}$，$ \mathit{\boldsymbol{w}}_k^{'}$为$V' $中任意元素，离散线单元可表示为$\mathit{\boldsymbol{u}}_k^{'} = \mathit{\boldsymbol{a}} + \mathit{\boldsymbol{w}}_1^{'} + \mathit{\boldsymbol{w}}_2^{'} + \cdots + \mathit{\boldsymbol{w}}_k^{'} $，则以矢量线ab的端点所在单元为起点与终点的离散线为$\mathit{\boldsymbol{a}}, \mathit{\boldsymbol{u}}_1^{'}, \mathit{\boldsymbol{u}}_2^{'} \cdots \mathit{\boldsymbol{u}}_N^{'}\left( \mathit{\boldsymbol{b}} \right)$。仅给定终点与起点，无法唯一确定离散线，因此给出如下限制条件：（1）离散线单元数目最小；（2）在满足（1）的所有离散线中，筛选出的离散线应使$ \mathop {{\rm{max}}}\limits_{1 \le k \le N} d\left( {\mathit{\boldsymbol{u}}_k^{'}, \mathit{\boldsymbol{ab}}} \right)$最小, $d\left( {\mathit{\boldsymbol{u}}_k^{'}, \mathit{\boldsymbol{ab}}} \right)$为$ \mathit{\boldsymbol{u}}_k^{'}$到ab的距离。满足限制条件的离散线是唯一的，且是ab对应的拟合精度最高的离散线，如图 4所示。

2.2 弱对偶三角形格网的离散线模型

如图 4所示的六边形格网中隐含对应的弱对偶三角形格网，即六边形格网离散线与相应的三角形格网离散线之间存在转化关系。如图 5所示，依据弱对偶关系，六边形格网离散线单元可分为两类，第一类中心与三角形单元中心重合，如蓝色点所示；第二类中心位于三角形单元顶点，如绿色点所示。令${\mathit{\boldsymbol{u}}_k} $和$ \mathit{\boldsymbol{u}}_k^{'}$分别表示以点k为中心的三角形离散线单元和六边形离散线单元。对于第一类情况，可直接将$ \mathit{\boldsymbol{u}}_k^{'}$的中心作为${\mathit{\boldsymbol{u}}_k} $的中心；对于第二类情况，需进一步讨论离散线单元中心的对应关系。

如图 6（a）所示，六边形格网离散线单元u'_A的中心位于三角形顶点，其相邻六边形格网离散线单元的中心（如图中蓝色点所示）对应两个三角形单元，当这两个三角形单元方向相同时，令u'_B为u'_A前一单元，向量$\mathit{\boldsymbol{v}}_0^{'}$、$ \mathit{\boldsymbol{v}}_1^{'}$为ab对应的最优方向向量，则$\mathit{\boldsymbol{u}}_A^{'}$对应的三角形单元u_c的中心可通过将$\mathit{\boldsymbol{u}}_B^{'}$的中心平移$\mathit{\boldsymbol{v}}_1^{'}$得到。

如图 6（b）所示，$\mathit{\boldsymbol{u}}_A^{'}$前后两个相邻六边形格网离散线单元的中心（如图中蓝色点所示）所对应的三角形单元方向相反。在此情况下，${\mathit{\boldsymbol{u}}_D}$、${\mathit{\boldsymbol{u}}_E}$、${\mathit{\boldsymbol{u}}_F}$、${\mathit{\boldsymbol{u}}_G}$为$\mathit{\boldsymbol{u}}_A^{'}$可能对应的三角形格网离散线单元，由$\mathit{\boldsymbol{u}}_A^{'}$和$\mathit{\boldsymbol{u}}_B^{'}$的中心平移$\mathit{\boldsymbol{v}}_0^{'}$，可分别得到${\mathit{\boldsymbol{u}}_D}$和${\mathit{\boldsymbol{u}}_E}$的中心，$\mathit{\boldsymbol{u}}_A^{'}$和$\mathit{\boldsymbol{u}}_C^{'}$的中心平移$ - \mathit{\boldsymbol{v}}_0^{'}$可分别得到${\mathit{\boldsymbol{u}}_F}$和${\mathit{\boldsymbol{u}}_G}$的中心。进一步根据$\mathit{\boldsymbol{u}}_A^{'}$、${\mathit{\boldsymbol{u}}_D}$、${\mathit{\boldsymbol{u}}_E}$的中心与矢量线的位置关系，可判断出$\mathit{\boldsymbol{u}}_A^{'}$对应的三角形单元，若$\mathit{\boldsymbol{u}}_A^{'}$、${\mathit{\boldsymbol{u}}_D}$、${\mathit{\boldsymbol{u}}_E}$的中心位于矢量线异侧，则$\mathit{\boldsymbol{u}}_A^{'}$对应的三角形单元为${\mathit{\boldsymbol{u}}_D}$、${\mathit{\boldsymbol{u}}_E}$；反之，则为${\mathit{\boldsymbol{u}}_F}$、${\mathit{\boldsymbol{u}}_G}$。

利用以上对应关系，以六边形格网离散线模型为基础，通过简单判断，可得到矢量线ab对应的三角形格网离散线模型，实现两类离散线模型的转化。

3 离散线模型的降维求解 3.1 六边形格网离散线模型一维等价形式

如§2.1所述，格网中的理想离散线需满足限制条件。确定最优方向向量后，可保证离散线单元数目最小。为了得到垂直距离$d\left( {\mathit{\boldsymbol{u}}_k^{'}, \mathit{\boldsymbol{ab}}} \right)$，过起点a做垂直于ab的直线H，令p_H为H上的垂直投影算子, 则最优方向向量集合$V' = \left\{ {\mathit{\boldsymbol{v}}_0^{'}, \mathit{\boldsymbol{v}}_1^{'}} \right\}$在H上的垂直投影, 构成一维直线上的向量集合$V = \left\{ {{\mathit{\boldsymbol{v}}_0} = {p_H}\left( {\mathit{\boldsymbol{v}}_0^{'}} \right), {\mathit{\boldsymbol{v}}_1} = {p_H}\left( {\mathit{\boldsymbol{v}}_1^{'}} \right)} \right\}$。令$W = \left\{ {{\mathit{\boldsymbol{w}}_1}, {\mathit{\boldsymbol{w}}_2} \cdots {\mathit{\boldsymbol{w}}_N}} \right\}{\rm{}}, {\mathit{\boldsymbol{w}}_k}$为V中任意元素，定义${\hat u_k} = {\mathit{\boldsymbol{w}}_1} + {\mathit{\boldsymbol{w}}_2} + \cdots + {\mathit{\boldsymbol{w}}_k}$, 连接$0,{\mathit{\boldsymbol{\hat u}}_1},{\mathit{\boldsymbol{\hat u}}_2} \cdots {\mathit{\boldsymbol{\hat u}}_N}$可形成一维直线上的路径，若${\hat u_N} = 0$, 则对应闭合路径$P:\left( {0, {{\mathit{\boldsymbol{\hat u}}}_1}, {{\mathit{\boldsymbol{\hat u}}}_2} \cdots {{\mathit{\boldsymbol{\hat u}}}_N}\left( 0 \right)} \right)$，该闭合路径偏离原点的最远距离可由$L\left( P \right) = \mathop {{\rm{max}}}\limits_{1 \le k \le N} \left\| {{{\mathit{\boldsymbol{\hat u}}}_k}} \right\|$（$\left\| {{{\mathit{\boldsymbol{\hat u}}}_k}} \right\|$为${\mathit{\boldsymbol{\hat u}}_k}$到原点的距离）表示。由于ab垂直于H，则a、b在H上的投影重合于一点，因此，六边形格网离散线$\mathit{\boldsymbol{a}}, \mathit{\boldsymbol{u}}_1^{'}, \mathit{\boldsymbol{u}}_2^{'} \cdots \mathit{\boldsymbol{u}}_N^{'}\left( \mathit{\boldsymbol{b}} \right)$与一维直线上的闭合路径$P:\left( {0, {{\mathit{\boldsymbol{\hat u}}}_1}, {{\mathit{\boldsymbol{\hat u}}}_2} \cdots {{\mathit{\boldsymbol{\hat u}}}_N}\left( 0 \right)} \right)$对应，且有${\mathit{\boldsymbol{\hat u}}_k} = {p_H}\left( {{{\mathit{\boldsymbol{\hat u}}}_k}} \right)$，$d\left( {\mathit{\boldsymbol{u}}_k^{'},\mathit{\boldsymbol{ab}}} \right) = \left\| {{{\hat u}_k}} \right\|$，如图 7所示。因此，求解二维平面内使$\mathop {{\rm{max}}}\limits_{1 \le k \le N} d\left( {\mathit{\boldsymbol{u}}_k^{'},\mathit{\boldsymbol{ab}}} \right)$ 达到最短的理想离散线等价于求解一维直线上min(L(P))对应的闭合路径^[10]。

3.2 三角形格网离散线求解算法

根据§3.1的分析，求解闭合路径的关键是使L(P)最小，采取“贪心”算法，每一步所选择的向量都保证当前路径最接近原点可达到目的。在六边形格网中给定矢量线ab，利用夹角最小原则确定最优方向向量集合V'，进一步将V'投影到直线H可得集合V。集合V'与V的元素对应，利用“贪心”算法求出每一步中集合V中的元素，由V'对应的最优方向向量可得弱对偶六边形格网中的理想离散线，再根据§2.2中离散线模型的转化关系，可得出相应三角形格网中的离散线。该算法实现了一维闭合路径→二维弱对偶六边形格网离散线→二维三角形格网离散线的求解过程，逐步简化三角形格网离散线生成问题，并大幅提升运算效率。算法具体流程如图 8所示。

4 对比实验及分析

为了验证本文算法的效果，将其分别与Freeman算法^[12]和全路径算法^[2]进行对比实验。3种算法均采用C#语言实现，除核心步骤外，其余操作完全相同，最大限度地保证结果客观、公平。全部代码均编译为Release版本，并在一台SAMSUNG 450R5U笔记本（硬件配置：Intel Core i5-3 230 MB CPU @ 2.60 GHz，8 GB RAM；操作系统：Windows 7 x64旗舰版；开发工具：Visual Studio 2012）上运行。实验数据来源于原国家测绘地理信息局1:25万基础地理信息数据库，包含全国34个省（直辖市、自治区）级行政区界。将每个省界视为封闭的任意多边形，对多边形的每一条边分别采用本文算法、Freeman算法和全路径算法进行处理，生成相应的离散线，其部分结果如图 9、图 10所示。图 9为细节效果对比，其中三角形单元边长约为3 m。图 10为整体效果对比，其中三角形单元边长约为3 km。

由图 9中的细节效果对比可以看出，本文算法与Freeman算法均可在三角形格网内生成离散线，但细节上存在明显差异。在倾斜、水平直线及转角处，Freeman算法生成的部分离散线单元并非直线经过的单元，离散线单元与直线的偏离较大，无法精确拟合直线。相比之下，本文算法在各种情况下均能选出直线经过的三角形单元，拟合效果理想。这种离散线与直线的精确拟合有利于离散化后的采样与分析。由图 10整体效果对比可看出，相比于Freeman算法，本文算法生成的离散线具有更加良好的连续性和平滑性。同时，当格网单元边长较大时，Freeman算法会出现直线仅沿格网单元边缘行进的情况，无法确保直线与离散线的精确拟合。

两种算法原理上的差别造成了以上效果差异。本文算法以格网单元中心到直线的垂直距离为依据筛选离散线单元，可确保离散线单元中心贴近直线。而Freeman算法依据单元顶点到直线的垂直距离选取离散线单元的方式则忽略了单元中心与直线的关系，容易造成格网单元偏离直线。

全路径算法涉及的大量坐标按距离排序的过程十分耗时，而本文算法仅在确定最优方向向量时涉及少量三角函数运算，后续过程经降维处理仅涉及加减运算，计算量远小于全路径算法，因而效率更高。对比实验表明，全路径算法与本文算法效果相同，但计算效率（时间比）差异较大。如表 1、图 11所示（单元边长约为3 m），本文算法的计算效率平均约为全路径算法的9.32倍。

5 结语

本文根据三角形格网与六边形格网之间的弱对偶关系，以及二维六边形格网离散线与一维闭合路径之间的等价关系，在三角形格网离散线、弱对偶六边形格网离散线、一维闭合路径之间建立转化关系，逐步降低问题的复杂度。

相比于现有算法，本文算法效果理想，且在效率上优势明显。算法涉及的弱对偶和降维思想为解决格网离散化问题提供了新的思路，对进一步解决格网问题具有重要的辅助作用。同时，该思想能够扩展到多维的情况，对于高维度问题的研究也具有重要意义。