数学杂志  2021, Vol. 41 Issue (3): 247-256   PDF    
扩展功能
加入收藏夹
复制引文信息
加入引用管理器
Email Alert
RSS
本文作者相关文章
冯峰
蒋维
一种带约束限制的三次B样条曲线矢量数据压缩算法
冯峰, 蒋维    
武汉大学数学与统计学院, 湖北 武汉 430072
摘要:为了便于大型矢量数据高效的检索分析,存储和传输,事先对矢量数据进行压缩是极为必要的.本文基于B样条良好的局部性和光滑性,利用带约束条件限制的三次B样条拟合方法对曲线矢量数据进行压缩.为了验证所提出算法的高效性,本文给出了9种不同的曲线矢量数据压缩算例,并同时与传统的Douglas-Peucker矢量压缩算法进行对比.数值算例表明,本文所提出的曲线矢量数据压缩算法明显优于传统的Douglas-Peucker压缩算法.该算法不仅能够保证曲线整体的二阶光滑性,还能够显著地降低数据的压缩率,因而具有广泛的应用前景(例如自动驾驶).
关键词曲线矢量数据压缩    三次B样条    整体C2-连续    Douglas-Peucker算法    约束条件    
A CUBIC B-SPLINE-BASED VECTOR DATA COMPRESSION ALGORITHM WITH BOUNDARY CONSTRAINTS
FENG Feng, JIANG Wei    
School of Mathematics and Statistics, Wuhan University, Wuhan 430072, China
Abstract: In order to efficiently retrieve, analyze, store and transmit large amount of vector data, it is extremely necessary to compress these vector data in advance. Based on elegant properties of the B-spline (e.g., locality and smoothness), we propose a cubic B-spline-based algorithm to compress the vector data with boundary constraints. The proposed cubic B-spline vector data compression algorithm is tested on nine examples with curve vector data. We also compare numerical results produced by the proposed algorithm with these of the classical Douglas-Peucker compression algorithm. Numerical results show that the proposed cubic B-spline-based vector compression algorithm not only can significantly reduce the compression rate, but also can produce highly accurate compression curve with C2-smoothness. Therefore, the algorithm has many important potential applications (e.g., automatic drive).
Keywords: curve vector data compression     cubic B-spline     C2-smoothness     Douglas-Peucker algorithm     boundary constraints    
1 引言

伴随着定位技术的高速发展和地理空间数据采集能力的不断增强, 地理空间数据的体量正变得越来越大, 因此对地理空间数据进行高效率的压缩, 以便后续的储存和传输显得尤为重要. 一般而言, 地理空间数据通常可以分为栅格数据和矢量数据两种数据结构[1]. 其中, 矢量数据因具有拓扑关系容易表达和冗余度较低等优点, 占据了地理空间数据的绝大部分. 曲线矢量数据压缩是地理空间数据压缩中一类很重要也是很常见的子问题, 例如: 在电子地图中的道路数据就可以被认为是一类曲线矢量数据, 它可以由一连串的采样于真实道路的坐标点来表示. 但是, 在大多数情况下, 电子地图中的道路矢量数据信息量需要保证其准确性以描述其位置和相应的拓扑关系, 因此往往数据量会很大, 势必会带来数据的冗余. 然而在智能交通, 无人驾驶等领域, 对高精度导航地图系统中的道路数据进行高效地检索、分析、存储和传输是尤为关键的. 矢量数据压缩能够极大地降低道路数据存储和传输的开销, 也会明显提高相应的数据处理能力, 因此针对矢量数据的压缩算法研究是十分必要的.

针对曲线矢量数据压缩, 传统的压缩算法[1-3]包括: 垂距限值法, 角度限值法, Douglas-Peucker算法[4], Li-OpenShaw简化算法[5], 预测压缩算法[6], 小波法[7, 8]等. 其中, 最经典的算法是Douglas-Peucker算法, 该算法可以看作一种从整体到局部的曲线矢量数据压缩后保留点集数据的方法. Douglas-Peucker算法可以递归式的通过编程实现, 并可以达到$ \mathcal{O}(n\log n) $的时间复杂度, 其基本思路为: 分别找到每段曲线中距离曲线首尾端点连线最大的中间点, 比较该最大距离与预先给定的误差阈值, 如果最大距离小于给定阈值则将该直线作为这段曲线的近似, 否则以此作为新的端点将该段曲线分为两段, 并对每段曲线重复上述过程, 直到所有曲线段处理完毕. Douglas-Peucker算法简单并且易于实现, 但是该方法却也存在一些问题, 例如选择不同的首尾点执行算法会得到不同的压缩结果, 压缩后的的曲线(即折线) 在节点处没有光滑性. 另外, 最主要的不足在于该算法没有考虑曲线矢量数据的拓扑关系, 如果没有进行必要的预处理的话, 压缩结果可能会改变原始数据的拓扑结构.

矢量数据的压缩本质上是由于原始数据信息存在一定的冗余, 因此要对数据进行压缩和结构的化简. 矢量数据压缩的核心在于在不扰乱数据拓扑关系的前提下, 找到一个尽可能少的数据集合能够在一定的精度范围内反映原始数据. 准确地来讲, 我们将矢量数据用一系列有序矢量点集[9]表示,

$ \begin{equation} \mathcal{P} = \left\{ {\bf{P}}_{1}, {\bf{P}}_{2}, {\bf{P}}_{3}, \ldots, {\bf{P}}_{|\mathcal{P}|-1}, {\bf{P}}_{|\mathcal{P}|} \right\}, \end{equation} $ (1.1)

矢量数据压缩的目标则是从该点集$ \mathcal{P} $中抽取或构造一个具有代表性的子集

$ \begin{equation} \mathcal{P}^{\prime} = \left\{ {\bf{P}}_{1}^{\prime}, {\bf{P}}_{2}^{\prime}, {\bf{P}}_{3}^{\prime}, \ldots, {\bf{P}}_{|\mathcal{P^\prime}|-1}^{\prime}, {\bf{P}}_{|\mathcal{P^\prime}|}^{\prime} \right\}, \end{equation} $ (1.2)

用来表示该矢量数据, 并且需要在保证拓扑结构等关系不改变的情况下在数据量上有尽可能小的压缩比. 我们将压缩比定义如下,

$ \begin{equation} \alpha = \frac{|\mathcal{P^\prime}|}{|\mathcal{P}|}, \end{equation} $ (1.3)

其中, $ |\mathcal{P}| $为原始矢量数据量的大小, $ |\mathcal{P^\prime}| $为压缩后的数据量大小. 通常, 压缩后的数据是指从原始矢量数据提取出的一个子集合, 这是因为大部分压缩算法是基于分析原始矢量数据并保留一部分点集. 但是, 压缩数据的作用是能够代表原始矢量数据, 因此压缩数据可以为任意点集, 只要其经由某种方式重构后能够尽可能的近似原始矢量数据.

为了克服传统Douglas-Peucker算法无法保证压缩后的曲线具有整体光滑性的缺陷, 本文采用三次B样条曲线拟合方法来构造相应的数据压缩算法. B样条是一类特殊的样条函数, 它可以描述为一组分段多项式基函数的线性组合, 其良好的局部性和光滑性使得基于B样条的拟合技术广泛应用于计算机辅助几何设计(CAD) [10, 11, 13], 信号处理[14, 15]和图像处理[16]等领域. 一般来说, 由于B样条具有良好的性质, B样条曲线仅需较少的参数就可以精确的表示非常复杂的形状, 因此利用B样条对矢量数据进行拟合可以仅使用较少的参数满足较高的误差精度要求. 本文将利用B样条拟合技术对曲线矢量数据进行压缩, 同时满足压缩数据过程中对首尾端点的约束要求.

2 B样条曲线简介

B样条曲线是一组分段多项式基函数的线性组合, 而这些基函数确定了曲线的连续性. 具体来说, 一条平面上的B样条曲线$ \bf{X}: \mathbb{R} \rightarrow \mathbb{R}^2 $, i.e., $ \rho \mapsto \bf{X}(\rho) $可以表示为[10]

$ \begin{equation} {\bf{X}}(\rho) = \sum\limits_{i = 0}^n {\bf{P}}_i N_{i, k}(\rho), \end{equation} $ (2.1)

其中系数$ {\bf{P}}_i = \left( x_{i}, y_{i} \right)^T \in \mathbb{R}^2, i = 0, 1, \ldots, n $, 是该B样条曲线的第$ i $个控制点, $ N_{i, k}(\rho) \in \mathbb{R}[\rho] $为定义在实数域$ \mathbb{R} $上的第$ i $$ k $次B样条基函数.

B样条基函数有许多等价定义, 而Cox-de Boor递归公式是应用最为广泛的一种标准方法. Cox-de Boor递归公式将一个高次B样条基函数表示为一些较低次B样条基函数的线性组合, 即对于所有正整数$ k>0 $来说, $ N_{i, k}(\rho) $为两个$ (k-1) $次B样条基函数的线性组合. 我们将第$ i $$ k $次B样条基函数的Cox-de Boor递归公式定义为:

$ \begin{equation} {N_{i, k}(\rho) = \frac{\rho-\rho_{i}}{\rho_{i+k}-\rho_{i}} N_{i, k-1}(\rho)+\frac{\rho_{i+k+1}-\rho}{\rho_{i+k+1}-\rho_{i+1}} N_{i+1, k-1}(\rho)}, \end{equation} $ (2.2)

其中$ i = 0, 1, 2, \ldots, n $, $ k \geq 1 $, 并且我们规定$ 0/0 $$ 0 $. 特别地, 我们有

$ \begin{equation} N_{i, 0}(\rho) = \left\{ \begin{array}{ll} {1, } & {\rho_{i} \leq \rho \leq \rho_{i+1}, } \\ {0, } & {\rm{otherwise}.} \end{array} \right. \end{equation} $ (2.3)

其中$ U = \{\rho_0, \cdots, \rho_{m}\} $被称为节点向量, $ \rho_{i} $称为节点, 节点个数必须满足$ m = n+k+1 $. $ U $是一个非递减序列, 即$ \rho_i\leq \rho_{i+1}, i = 0, \ldots, m-1 $. 图 1-3分别展示了次数为$ 1 $, $ 2 $$ 3 $次的B样条基函数, 并同时给出构成其函数的两个较低次的B样条基函数. 例如, 在图 3中, 三次B样条基函数$ N_{i, 3} $由两个较低次B样条基函数$ N_{i, 2} $$ N_{i+1, 2} $线性组合而成.

图 1 一次B样条基函数

图 2 二次B样条基函数

图 3 三次B样条基函数

从递推式可以看出, B样条基函数有许多有用的性质, 我们将一些常用且重要的性质[10-13] 罗列如下:

$ \bullet $  非负性和局部支撑性: 对于任意正整数$ k \geq 0 $, $ N_{i, k}(\rho) $为定义在节点向量上的一个非负多项式, 并且我们有: 当$ \rho \notin \left[ \rho_i, \rho_{i+k+1} \right] $时, $ N_{i, k}(\rho) = 0 $.

$ \bullet $  规范性: 对于任意正整数$ k \geq 0 $, $ \sum_i N_{i, k}(\rho) \equiv 1 $.

$ \bullet $  连续可微性: 对于任意正整数$ k \geq 1 $, B样条基函数在重复度为$ r $的节点处是$ k-r $次可微的, 即$ N_{i, k}(\rho) \in C^{k-r}(\mathbb{R}) $; 而在节点内部, 基函数$ N_{i, k}(\rho) $是无限次可微的, 即$ N_{i, k}(\rho) \in C^{\infty}(\mathbb{R} \setminus U) $.

$ \bullet $  线性无关性: 一组B样条基函数$ \left\{ N_{i, k}(\rho), N_{i+1, k}(\rho), \ldots, N_{i+k, k}(\rho) \right \} $在区间$ \left[ \rho_{i}, \rho_{i+k+1} \right] $是线性无关的.

利用B样条基函数的连续可微性, B样条曲线$ \bf{X}(\rho) $$ r $次微分可以通过B样条基函数的$ r $次微分直接计算得到, 即

$ \begin{equation} {\bf{X}}^{(r)}(\rho) = \sum\limits_{i = 0}^n {\bf{P}}_i N_{i, k}^{(r)}(\rho). \end{equation} $ (2.4)

大部分矢量数据是采样于具有连续性的曲线, 例如电子地图中的道路数据主要由一系列的线段, 圆弧和缓和曲线等图形要素组成, 并且它们的连接处具有一定的光滑性(一般而言, 它们所代表的图形要素组合需要满足直到二阶的可微性). 对于这类矢量数据, 数据压缩算法不仅要保证在压缩率尽可能小(即压缩后的点尽可能少)的情况下能够尽可能好地保证其误差精度要求, 还需要满足其隐含的光滑性要求. 在本文中, 我们将给出基于带约束条件的三次B样条曲线拟合矢量压缩算法. 由于B样条具有良好的性质, 只需较少的参数就可以重构出较为复杂的形状, 因此将B样条应用于此类矢量数据压缩可以很好地满足实际需求.

3 三次B样条曲线拟合

一般来说, 一组给定曲线矢量数据是从某些具有一定光滑性的曲线中采样得到的一系列样本数据点$ \{{\bf{Q}}_k\}_{k = 0}^K\subset \mathbb{R}^2 $. 在给定控制点数量$ n $的条件下, 三次B样条曲线拟合技术即意味着需要找到一个阶数为$ 3 $的B样条曲线$ {\bf{X}}(\rho) = \sum_{i = 0}^n {\bf{P}}_i N_{i, 3}(\rho) $满足所有数据点到样条曲线的误差和最小, 即

$ \begin{equation} \min \ \sum\limits_{k = 0}^{K}\left|{\bf{Q}}_{k}-{\bf{X}}\left(\overline{\rho}_{k}\right)\right|^{2}, \end{equation} $ (3.1)

$ \overline{\rho}_{k} $为预先给定的数据点参数值, 这里我们利用积累弦长法确定, 即有

$ \begin{equation} \overline{\rho}_k = \overline{\rho}_{k-1}+\frac{\left|{\bf{Q}}_k-{\bf{Q}}_{k-1}\right|}{\sum\limits_{i = 1}^{K}\left|{\bf{Q}}_i-{\bf{Q}}_{i-1}\right|} , \ k = 1, \ldots, K-1, \end{equation} $ (3.2)

其中$ \overline{\rho}_0 = 0 $, $ \overline{\rho}_{K} = 1 $. 另外, 在不引起误会的前提下, 在下面我们将阶数为$ 3 $的B样条基函数$ N_{i, 3}(\rho) $简记为$ N_i(\rho) $.

由B样条曲线的定义可以知道[13], 一条三次B样条曲线$ {\bf{X}}(\rho) = \sum_{i = 0}^n {\bf{P}}_i N_i(\rho) $由其控制点$ {\bf{P}}_i $和定义其基函数$ N_i(\rho) $的节点向量$ U = \{\rho_0, \cdots, \rho_{m}\} $唯一决定. 为了能够利用尽可能少的信息来表示曲线矢量数据, 我们预先给定相应的节点向量, 在此我们选择准均匀分布的节点向量, 即定义域两端节点具有$ 4 $阶重复度, 其余节点为均匀分布: $ \rho_{-3} = \rho_{-2} = \rho_{-1} = \rho_0 = 0, \rho_1 = h, \ldots, \rho_{n-1} = (n-1)h, \rho_n = \rho_{n+1} = \rho_{n+2} = \rho_{n+3}( = nh) = 1 $, 其中$ h = 1/n $. 需要注意的是, 如此构造的三次B样条曲线的两个首尾端点将与其对应的首尾控制点重合, 即:

$ \begin{equation} {\bf{X}}(0) = {\bf{P}}_{0}, \qquad {\bf{X}}(1) = {\bf{P}}_{n}. \end{equation} $ (3.3)

此时, 当采用上述准均匀分布的节点向量, 并使用积累弦长法对数据点进行参数化时, 上述极小化问题(8) 便可以转化为求解未知控制点的线性最小二乘问题[10, 13]. 应用标准的最小二乘技术, 在给定控制点个数$ n $的前提下, 上述极小化问题(3.1) 可以转化为关于下列矩阵迹的极小问题:

$ \begin{equation} \min\limits_{\bf{P}} \ {\textbf{tr}}\bigl[({\bf{Q}}-{\bf{N}}{\bf{P}})^T({\bf{Q}}-{\bf{N}}{\bf{P}})\bigr], \end{equation} $ (3.4)

其中, $ \textbf{tr}(\cdot) $表示求矩阵的迹, $ \bf{N} $$ (K+1)\times(n+1) $型矩阵, 即

$ \begin{equation} \bf{N} = \begin{bmatrix} {N_{0}\left(\overline{\rho}_{0}\right)} & {N_{1}\left(\overline{\rho}_{0}\right)} & {\cdots} & {N_{n}\left(\overline{\rho}_{0}\right)} \\ {N_{0}\left(\overline{\rho}_{1}\right)} & {N_{1}\left(\overline{\rho}_{1}\right)} & {\cdots} & {N_{n}\left(\overline{\rho}_{1}\right)} \\ {\vdots} & {\vdots} & {\ddots} & {\vdots} \\ {N_{0}\left(\overline{\rho}_{K}\right)} & {N_{1}\left(\overline{\rho}_{K}\right)} & {\dots} & {N_{n}\left(\overline{\rho}_{K}\right)} \end{bmatrix}, \end{equation} $ (3.5)

$ {\bf{Q}} = \left[{{\bf{Q}}_{0}} , \ {{\bf{Q}}_{1}} , \ {\dots} , \ {{\bf{Q}}_{K}}\right]^T $$ (K+1) \times 2 $型二维数据向量集, $ {\bf{P}} = \left[{{\bf{P}}_{0}} , \ {{\bf{P}}_{1}} , \ {\dots} , \ {{\bf{P}}_{n}}\right]^T $$ (n+1) \times 2 $型未知控制点二维向量集. 对极小化问题(11) 中的目标函数求关于控制点向量$ \bf{P} $的偏导数, 并将其置为$ \bf{0} $向量, 可以得到一个以控制点向量为未知量的线性代数方程组:

$ \begin{equation} \left({\bf{N}}^{T}{\bf{N}}\right) {\bf{P}} = {\bf{N}}^T {\bf{Q}}. \end{equation} $ (3.6)

并且, 由Schoenberg-Whitney定理[10]可推知, $ \bf{N}^{T}\bf{N} $是非奇异的. 于是, 当利用积累弦长法对数据点进行参数化, 并采用准均匀分布的节点向量构造相应的三次B样条基函数时, 在给定三次B样条控制点数目条件下, 数据点的三次B样条曲线拟合所得到的结果是唯一的:

$ \begin{equation} {\bf{P}} = \left({\bf{N}}^{T}{\bf{N}}\right)^{-1} {\bf{N}}^T {\bf{Q}}. \end{equation} $ (3.7)

这也就意味着对于任意曲线矢量数据来说, 在给定控制点数目前提下, 应用三次B样条拟合技术对其进行曲线拟合可以得到唯一的控制点序列.

4 一种带约束的三次B样条曲线矢量数据压缩算法

在上一节我们考虑了将曲线矢量数据作为一个整体应用最小二乘三次B样条拟合技术进行数据拟合. 然而, 在实际应用情况下, 对于大部分曲线矢量数据压缩来说, 其首尾端点需要满足必要的约束条件. 本节中, 我们将基于一种带约束的三次B样条曲线拟合算法, 提出一种新的带约束三次B样条拟合压缩算法.

与上节类似, 带约束三次B样条最小二乘曲线拟合问题可表述为: 给定一组曲线矢量数据集$ \{{\bf{Q}}_k\}_{k = 0}^K\subset \mathbb{R}^2 $以及控制点个数$ n $, 需要找到一个阶数为$ 3 $的B样条曲线$ {\bf{X}}(\rho) = \sum_{i = 0}^n {\bf{P}}_i N_{i}(\rho) $, 使得所有数据点到样条曲线的误差和最小, 即式(3.1), 同时该样条曲线在首尾端点满足给定的具体值: $ {\bf{X}}(0) $, $ {\bf{X}}(1) $, $ {\bf{X}}^{\prime\prime}(0) $, $ {\bf{X}}^{\prime\prime}(1) $.

与上一节相同, 我们同样采用准均匀分布的节点向量, 并使用积累弦长法对数据点进行参数化. 利用式(2.4), 我们可以得出如下等式:

$ \begin{equation} {\bf{X}}^{\prime\prime}(0) = 6(n-2)^2{\bf{P}}_{0}-9(n-2)^2{\bf{P}}_{1}+3(n-2)^2{\bf{P}}_{2}, \end{equation} $ (4.1)
$ \begin{equation} {\bf{X}}^{\prime\prime}(1) = 6(n-2)^2{\bf{P}}_{n}-9(n-2)^2{\bf{P}}_{n-1}+3(n-2)^2{\bf{P}}_{n-2}. \end{equation} $ (4.2)

结合约束条件(3.3), 该问题可以转化为带等式约束的二次规划问题[10], 写成向量形式为:

$ \begin{equation} \begin{array}{l} \min\limits_{\bf{P}} \ \textbf{tr}\bigl[({\bf{Q}}-{\bf{N}}{\bf{P}})^T({\bf{Q}}-{\bf{N}}{\bf{P}})\bigr], \\ \rm{subject to} \quad {\bf{M}}{ \bf{P}} = {\bf{T}}, \end{array} \end{equation} $ (4.3)

其中, $ {\bf{N}} $$ (K+1) \times (n+1) $矩阵, $ {\bf{Q}} = \left[{{\bf{Q}}_{0}} , \ {{\bf{Q}}_{1}} , \ {\dots} , \ {{\bf{Q}}_{K}}\right]^T $$ (K+1) \times 2 $二维数据向量, $ {\bf{P}} = \left[{{\bf{P}}_{0}} , \ {{\bf{P}}_{1}} , \ {\dots} , \ {{\bf{P}}_{n}}\right]^T $$ (n+1) \times 2 $未知控制点向量, 而$ {\bf{M}} $$ {\bf{T}} $分别为等式约束线性方程组的$ 4 \times (n+1) $系数矩阵和$ 4 \times 2 $维系数向量:

$ \begin{equation} \bf{M} = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 0 & 1 \\ 6(n-2)^2 & -9(n-2)^2 & 3(n-2)^2 & \cdots & 0 & 0 & 0 \\ 0 & 0 & 0 & \cdots & 3(n-2)^2 & -9(n-2)^2 & 6(n-2)^2 \\ \end{bmatrix}, \end{equation} $ (4.4)
$ \begin{equation} {\bf{T}} = \left[{\bf{X}}(0), \ {\bf{X}}(1), \ {\bf{X}}^{\prime\prime}(0), \ {\bf{X}}^{\prime\prime}(1)\right]^T. \end{equation} $ (4.5)

利用Lagrange乘子法, 令$ {\bf{A}} = \left[\lambda_{ij}\right] $$ 4\times 2 $维Lagrange乘子向量, 则带等式约束的二次规划问题(4.3) 可以转化为如下最小化问题:

$ \begin{equation} \min\limits_{\bf{P}}\, \textbf{tr}\bigl[ \left({\bf{Q}}-{\bf{N}}{\bf{P}}\right) ^T \left({\bf{Q}}-{\bf{N}}{\bf{P}}\right)+{\bf{A}}^{T}({\bf{M}} {\bf{P}}-{\bf{T}})\bigr]. \end{equation} $ (4.6)

分别对$ \bf{A} $$ \bf{P} $求导并将结果取成$ \bf{0} $向量, 得极小解满足:

$ \begin{equation} \left\{ \begin{array}{lll} ({\bf{N}}^T{\bf{N}}){\bf{P}}-{\bf{N}}^T{\bf{Q}}+{\bf{M}}^T{\bf{A}} & = &{\bf{0}}, \\ \bf{M} \bf{P}-\bf{T} & = &{\bf{0}}. \end{array} \right. \end{equation} $ (4.7)

由此可以得到分块矩阵形式如下

$ \begin{equation} \begin{bmatrix}{\bf{N}^{T}\bf{N}} & {\bf{M}^{T}} \\ {\bf{M}} & {\bf{0}}\end{bmatrix} \begin{bmatrix}{\bf{P}} \\ {\bf{A}}\end{bmatrix} = \begin{bmatrix}{\bf{N}^{T} \bf{Q}} \\ {\bf{T}}\end{bmatrix}. \end{equation} $ (4.8)

$ {\bf{N}}^{T}{\bf{N}} $$ {\bf{M}}\left({\bf{N}}^{T} {\bf{N}}\right)^{-1} {\bf{M}}^{T} $非奇异, 则有唯一解:

$ \begin{equation} {\bf{P}} = \left({\bf{N}}^{T} {\bf{N}}\right)^{-1} {\bf{N}}^{T} {\bf{Q}}-\left({\bf{N}}^{T} {\bf{N}}\right)^{-1} {\bf{M}}^{T} {\bf{A}}, \end{equation} $ (4.9)
$ \begin{equation} {\bf{A}} = ({\bf{M}}\left({\bf{N}}^{T}{\bf{N}}\right)^{-1} {\bf{M}}^{T})^{-1}({\bf{M}}\left({\bf{N}}^{T} {\bf{N}}\right)^{-1} {\bf{N}}^{T} {\bf{Q}}-{\bf{T}}). \end{equation} $ (4.10)

同样由Schoenberg-Whitney定理[10]可推知, 当采用本文所述的准均匀分布的节点向量, 并使用积累弦长法对数据点进行参数化时, $ {\bf{N}}^{T}{\bf{N}} $是非奇异的, 又由于$ \bf{M} $矩阵是行满秩的, 于是$ {\bf{M}}\left({\bf{N}}^{T} {\bf{N}}\right)^{-1} {\bf{M}}^{T} $也是非奇异的. 因此, 在给定三次B样条控制点数目条件下, 对数据点的带约束的三次B样条曲线拟合是唯一的. 图 4给出了利用带约束条件限制的三次B样条拟合示例, 其中控制点数目分别选为$ 8 $$ 16 $, 给定边界条件$ {\bf{X}}(0) = {\bf{Q}}_0 $, $ {\bf{X}}(1) = {\bf{Q}}_K $, $ {\bf{X}}^{\prime\prime}(0) = \bf{0} $, $ {\bf{X}}^{\prime\prime}(1) = \bf{0} $.

图 4 带约束条件限制的三次B样条曲线拟合, 左图: 8个控制点; 右图: 16个控制点.

基于上述带约束的三次B样条曲线拟合方法, 我们很容易提出一种新的三次B样条曲线矢量数据压缩算法. 该算法主要步骤如下:

(1) 利用积累弦长法(3.2) 将曲线矢量数据参数化$ \{\left({\bf{Q}}_k, \ \overline{\rho}_{k}\right)\}_{k = 0}^K $.

(2) 给定误差阈值$ \varepsilon $, 选取三次B样条初始控制点数目$ n $ (这里的$ n $值初始选取不应过大)求解带等式约束的二次规划问题, 由式(4.9) 得到相应的控制点序列, 其中三次B样条的节点向量设定为准均匀分布.

(3) 计算所有数据点到拟合的三次样条曲线的误差, 如果数据点到样条曲线的最大误差小于给定误差阈值$ \varepsilon $:

$ \begin{equation} \max \limits_{k = 0, \ldots , K}\left|{\bf{Q}}_{k}-{\bf{X}}\left(\overline{\rho}_{k}\right)\right| < \varepsilon, \end{equation} $ (4.11)

则将该控制点序列作为该矢量数据的压缩结果; 否则逐渐增大控制点数量, 重复第(2) 步.

在实际计算中, 上述算法可以与二分法结合, 从而快速地找到既满足误差精度要求, 同时又使得控制点(即压缩存储点)的数量尽可能少的最优存储曲线. 我们将结合二分法的三次B样条曲线矢量数据压缩算法详细步骤在下列算法1中给出.

5 数值结果比较

为了测试本文所提出的带约束条件限制三次B样条曲线矢量数据压缩算法, 同时与传统的Douglas-Peucker算法[4]进行对比, 本节给出相应的数值算例并比较其数据压缩率结果. 我们将三次B样条曲线矢量数据压缩算法的压缩率定义为

$ \begin{equation} \alpha_{\rm{bsp}} = \frac{\rm{三次\; B\; 样条曲线中控制点数目}}{\rm{曲线矢量数据中矢量点数目}}, \end{equation} $ (5.1)

而Douglas-Peucker算法的压缩率定义为

$ \begin{equation} \alpha_{\rm{dp}} = \frac{\rm{压缩后剩余矢量点数目}}{\rm{曲线矢量数据中矢量点数目}}. \end{equation} $ (5.2)

图 5展示了$ 9 $种不同的曲线矢量数据, 它们都是从随机选取的连续曲线中等距采样得到的含有$ 1000 $个点的矢量数据集.

图 5 一些曲线矢量数据的例子

下面分别利用本文所提出的三次B样条曲线矢量数据压缩算法和Douglas-Peucker算法对这些曲线矢量数据进行压缩, 表 1表 2分别给出了在不同误差阈值$ \varepsilon $下利用Douglas-Peucker算法和本文提出的算法对这些曲线矢量数据进行压缩的压缩率结果. 从表 1表 2中可以看出, 对于图 5中的9种不同的曲线矢量数据, 在三种给定的误差阈值下, 本文所提出的三次B样条曲线矢量数据压缩算法对曲线矢量数据进行压缩的压缩率结果均明显优于经典的Douglas-Peucker算法.

表 1 Douglas-Peucker算法的压缩率αdp数值结果

表 2 基于三次样条压缩算法的压缩率αbsp的数值结果
6 结论

基于三次B样条良好的局部性和光滑性, 本文将三次B样条应用于曲线矢量数据压缩. 通过给定矢量数据端点边界条件限制, 矢量数据压缩的问题首先被转化为一类带约束条件的三次B样条最小二乘曲线拟合问题, 进而通过求解带等式约束的二次规划问题得到三次B样条曲线控制点作为压缩结果. 由数值实验结果可以看出, 本文所提出的三次B样条矢量数据算法, 相对于传统的Douglas-Peucker矢量数据压缩算法, 在给定相同的误差阈值条件下, 能够明显的降低矢量数据压缩的压缩率. 同时, 由于继承了B样条曲线的性质, 由本算法压缩后的曲线还整体地满足$ C^2 $连续性. 除此之外, 本文所提出的算法由于额外增加了端点边界条件的限制, 更能够体现出原始矢量数据的特征. 需要强调的是, 本文所提出的三次B样条压缩算法可以很容易地推广到多维(例如三维矢量数据) 情况, 并且由于B样条良好的局部性, 针对大数据量的情况, 三次B样条压缩算法也可以利用并行计算来显著提高算法的效率.

参考文献
[1] 温永宁, 等. 矢量空间数据渐进传输研究进展[J]. 地理与地理信息科学, 2011, 27(6): 6–12.
[2] 张振鑫, 张维, 刘嫔, 等. 矢量地图数据简化研究进展[J]. 测绘工程, 2016, 25(6): 10–14. DOI:10.3969/j.issn.1006-7949.2016.06.003
[3] 杨建宇, 杨崇俊, 明冬萍, 等. WebGIS系统中矢量数据的压缩与化简方法综述[J]. 计算机工程与应用, 2005, 40(32): 36–38.
[4] Douglas D, Peucker T. Algorithms for the reduction of the number of points required to represent a line or its caricature[J]. The Canadian Cartographer, 1973, 10(2): 112–122. DOI:10.3138/FM57-6770-U75U-7727
[5] Li Z, Openshaw S. Algorithms for automated line generalization based on a natural principle of objective generalization[J]. International Journal of Geographical Information Systems, 1992, 6(5): 373–389. DOI:10.1080/02693799208901921
[6] 黄培之. 具有预测功能的曲线矢量数据压缩方法[J]. 测绘学报, 1995, 024(4): 316–320.
[7] 马伯宁, 冷志光, 汤晓安, 等. 具有误差修正的线矢量数据小波变换[J]. 计算机辅助设计与图形学学报, 2011(11): 1825–1829.
[8] 吴凡. 基于小波分析的线状特征数据无级表达[J]. 武汉大学学报(信息科学版), 2004, 29(6): 488–491.
[9] 单玉香. 矢量数据压缩模型与算法的研究[D]. 太原理工大学, 2004.
[10] Piegl L, Tiller W. The NURBS book[M]. Berlin: Springer, 1997.
[11] De Boor C. A practical guide to splines, volume 27[M]. New York: Springer-verlag, 1978.
[12] Kincaid D, Cheney W. Numerical analysis: mathematics of scientific computing, volume 2[M]. American Mathematical Society, 2009.
[13] 施法中. 计算机辅助几何设计与非均匀有理B样条[M]. 高等教育出版社, 2001.
[14] Unser M, Aldroubi A, Eden M. B-spline signal processing. i. theory[J]. IEEE Transactions on Signal Processing, 1993, 41(2): 821–833. DOI:10.1109/78.193220
[15] Unser M, Aldroubi A, Eden M. B-spline signal processing. ii. efficiency design and applications[J]. IEEE Transactions on Signal Processing, 1993, 41(2): 834–848. DOI:10.1109/78.193221
[16] Unser M. Splines: A perfect fit for signal and image processing[J]. IEEE Signal Processing Magazine, 1999, 16(6): 22–38. DOI:10.1109/79.799930