带参数的三重两步插值细分方法

Parameterized Ternary Two-Step Interpolatory Subdivision Scheme

摘要: 两步插值细分格式是一种新型式的细分格式，它具有在偶数次加细时插值控制网格的特点。这种特点使得它在计算机辅助几何设计中具有潜在的吸引力。本文给出了一种新的三重两步插值曲线细分格式。证明了通过调节参数其极限曲线可以达到$C^2$ 连续。其次，本文还分析了该细分格式的Hӧlder指数，并且给出了其基函数图像。数值实例表明，该两步插值细分方法可以使极限曲线较好的接近初始控制多边形并保证插值，该方法是合理有效的。

Abstract: The two-step interpolatory subdivision scheme is a novel type of subdivision method characterized by interpolating the control mesh at even refinement steps. Such a feature renders it potentially attractive in computer-aided geometric design (CAGD). In this paper, we propose a new ternary two-step interpolatory curve subdivision scheme. It is proved that the limit curve may attain $C^2$ -continuity by adjusting parameters. Furthermore, the Hölder exponent of the scheme is analyzed, and the graph of the basic limit function is provided. Numerical examples demonstrate that the proposed two-step interpolatory subdivision method yields limit curves that closely approximate the initial control polygon while ensuring interpolation, thus verifying the method’s reasonableness and effectiveness.

1. 引言

在计算机辅助设计领域之中，细分方法是生成光滑曲线和曲面的重要工具。细分方法的基本原理是给出一组初始控制点集，通过不断地重复应用加细规则来生成新的控制点集。随着迭代次数的增加，控制点集中的控制点会逐渐逼近形成一个光滑的曲线或曲面。

对于细分方法可以从多个角度进行分类。从拓扑规则来看，细分格式可以分为基本型细分[1]和对偶型细分[2]两类；根据细分过程中生成的极限曲线是不是通过初始控制点，细分方法还可以以生成曲线的特点划分为逼近型细分[3] [4]和插值型细分[5] [6]。通过逼近型细分方法生成的几何形状不一定经过所有初始控制点，而插值型细分方法则确保生成的几何形状能够精确地通过所有给定的控制点，从而为精确的几何建模提供了更为有力的支持。

细分方法的发展历程可以追溯至20世纪70年代，Chaikin算法[7]通过迭代割角生成二次B样条曲线，开启了细分方法的先河。1978年，Doo-Sabin [8]和Catmull-Clark [9]分别提出针对张量积曲面和任意拓扑网格的细分算法，奠定了逼近型细分的基础。然而，这些方法生成的曲线曲面并不通过初始控制点，无法满足插值的需求。插值型细分的突破出现在1987年，Dyn [10]等人提出四点插值细分方案，使细分后的曲线严格通过初始控制点，实现了几何精确性与光滑性的平衡。这一里程碑式的工作进程为多步插值细分方法奠定了理论基础。

多步插值细分格式及其可细化函数在计算机辅助几何设计(CAGD)、数值偏微分方程和逼近理论中都有应用。多步插值细分方法的发展源于计算机辅助几何设计和数值分析领域对高精度、高光滑性曲线曲面生成的需求。基于以上，且鉴于多步插值的优质性质，本文建立了一种新的三重两步插值曲线细分方法。证明了通过调节参数其极限曲线可以达到$C^2$ 连续。其次，还分析了该格式的Hӧlder指数并给出了其基函数。数值实例表明，该两步插值细分方法可以使极限曲线较好的接近初始控制多边形并保证插值，该方法是合理有效的。

2. 基础知识

定义2.1 设初始控制点集为$P^0=\left\{p_j^0\in {ℝ}^d,j\in \mathbb{Z}\right\}$ ，通过对$P^0$ 持续添加新点以生成细分曲线。设$P^k=\left\{p_j^k\in {ℝ}^d,j\in \mathbb{Z}\right\}$ 为第$k$ 次细分后生成的新控制点集，其中$k\ge 0$ 且$k\in \mathbb{Z}$ 。则m重细分格式的加细过程可表示为

$p_i^{k+1}={\displaystyle {\sum }_{j\in Z}{a}_{i-mj}{p}_j^k},i\in \mathbb{Z},$

其中$a=\left\{a_i,i\in \mathbb{Z}\right\}$ 被称为掩模。记$P^{k+1}=S_a{P}^k$ ，$S_a$ 称为细分算子，可简记为S。

定义2.2 设细分掩模$a=\left\{a_i,i\in \mathbb{Z}\right\}\in l_0\left(\mathbb{Z}\right)$ ，$l_0\left(\mathbb{Z}\right)$ 是一个定义在$\mathbb{Z}$ 上的有限的支撑序列集合，其对应的Laurent多项式则记为$a\left(z\right)={\displaystyle {\sum }_{i\in Z}{a}_i{z}^i}$ ，称之为细分掩模a对应的生成函数。如果掩模满足

$a\left(c_a-k\right)=a\left(k\right),\forall k\in \mathbb{Z},c_a\in \mathbb{Z}.$

则可以称该掩模是关于$\frac{c_a}{2}$ 对称的。当$c_a=1$ 时，对应的细分格式称为对偶型细分格式。当$c_a=0$ 时，对应的细分格式则称为基本型细分格式。

对于满足${\displaystyle {\sum }_{i\in Z}a\left(i\right)}=m$ 的序列，通过傅里叶变换可以得到$\widehat{\phi }\left(\xi \right)={\displaystyle {\prod }_{j=1}^{\infty }\frac{1}{m}a}\left({\text{e}}^{-i{m}^{-j}\xi }\right),\xi \in ℝ$。其中傅里叶变换的定义为$\widehat{f}\left(\xi \right)={\displaystyle \int f\left(x\right){\text{e}}^{-ix\xi }\text{d}x}$，$\widehat{\phi }\left(0\right)=1$ ，$\phi$ 是m重可加细函数，满足加细方程 $\phi ={\displaystyle {\sum }_{k\in Z}a\left(k\right)\phi \left(m-k\right)}$ 。

定义2.3 如果掩模a满足

$a\left(z\right)={\left(1+z+z^2+z^3+\cdots +z^{m-1}\right)}^{J}b\left(z\right),b\in l_0\left(Z\right),J\in N_0$ .

则称该掩模a具有$J$ 阶求和规则。为方便起见，定义$\text{sr}\left(a,M\right)=J$ ，其中$J$ 是满足上述条件的最大整数。

引理2.1 [11]若m重细分格式$S$ 的n阶差分格式为$S_n$ ，则存在一个该细分格式$S$ 的一阶差分格式的细分格式$S_1$ ，使得

$d{P}^k=S_{1}d{P}^{k-1}$ ,

其中，$P^k=S^k{P}^0$ ，$d{P}^k=\left\{{\left(d{P}^k\right)}_i=m^k\left(p_{i+1}^k-p_i^k\right),i\in Z\right\}$ 。

当且仅当对任意初始控制网格$f^0$ ，$\frac{1}{m}{S}_1$ 必须收敛到零函数时，细分格式$S$ 才是一致收敛的，此时，$S_n$ 的生成多项式为

$a^{\left(n\right)}\left(z\right)={\displaystyle {\sum }_{i\in Z}{a}_i^{\left(n\right)}{z}^i}={\left(\frac{mz}{1+z+\cdots +z^{m-1}}\right)}^{n}a\left(z\right)$ .

引理2.2 [12]如果$m$ 重细分格式$S$ 的掩模a和差分格式$S_j\left(j=1,2,\cdots ,n+1\right)$ 均存在，并且$\frac{1}{m}{S}_{n+1}$ 对于任意的初始控制网格$f^0$ 都收敛到零函数，那么m重细分格式$S$ 是$C^n$ 连续的。

为证明细分格式$S$ 是$C^n$ 连续的，只需要证明$n+1$ 阶差分格式$S_{n+1}$ 存在，并且存在一个正整数$L$ ，使得${\Vert {\left(\frac{1}{m}{S}_{n+1}\right)}^L\Vert }_{\infty }<1$ 。其中

${\Vert {\left(\frac{1}{m}{S}_{n+1}\right)}^L\Vert }_{\infty }=\max \left\{{\displaystyle {\sum }_{j\in Z}\left|a_{i+m^{L}j}^{\left[n+1,L\right]}\right|},i=0,1,\cdots ,m^L-1\right\},$

$a^{\left[n+1,L\right]}\left(z\right)={\displaystyle {\sum }_{i,j}{a}_{i+m^{L}j}^{\left[n+1,L\right]}{z}^{i+m^{L}j}}=\frac{1}{m^L}{\displaystyle {\prod }_{j=0}^{L-1}{a}^{\left(n+1\right)}\left(z^{mj}\right)}.$

3. 单参数三重两步插值细分

3.1. 单参数三重两步插值细分格式

给定初始的控制点集$P^0=\left\{p_i^0\in {ℝ}^d,i\in \mathbb{Z}\right\}$ ，现设第k次细分后的控制点集为$P^k=\left\{p_i^k\in {ℝ}^d,i\in \mathbb{Z}\right\}$ 。根据定义2.1可定义细分规则为$p_i^{k+1}={\displaystyle {\sum }_{j\in Z}{a}_{i-3j}{p}_j^k},i\in \mathbb{Z}$ 。本文给出新的单参数三重两步插值细分格式：

$\{\begin{array}{l}{p}_{3i}^{k+1}=\frac{-3-4b_{-1}}{216}{p}_{i-2}^k+\left(\frac{19}{144}+\frac{b_{-1}}{27}\right)p_{i-1}^k+\frac{137}{144}{p}_i^k-\left(\frac{11}{144}+\frac{b_{-1}}{27}\right)p_{i+1}^k+\frac{3+8b_{-1}}{432}{p}_{i+2}^k,\\ p_{3i+1}^{k+1}=\frac{3+8b_{-1}}{432}{p}_{i-2}^k-\left(\frac{11}{144}+\frac{b_{-1}}{27}\right)p_{i-1}^k+\frac{137}{144}{p}_i^k+\left(\frac{19}{144}+\frac{b_{-1}}{27}\right)p_{i+1}^k+\frac{-3-4b_{-1}}{216}{p}_{i+2}^k,\\ p_{3i+2}^{k+1}=-\frac{1}{16}{p}_{i-1}^k+\frac{9}{16}{p}_i^k+\frac{9}{16}{p}_{i+1}^k-\frac{1}{16}{p}_{i+2}^k.\end{array}$

该掩模带一个参数$b_{-1}$ ，满足$a_i=a_{1-i},i=1,2,\cdots ,7$ ，对应的生成函数为

$\begin{array}{l}a\left(z\right)=z^{-6}\left(z+1\right){\left(z^2+z+1\right)}^3(\left(8b_{-1}+3\right)z^6-2\left(20b_{-1}+9\right)z^5+2\left(44b_{-1}+9\right)z^4\\ -2\left(56b_{-1}-9\right)z^3+2\left(44b_{-1}+9\right)z^2-2\left(20b_{-1}+9\right)z+8b_{-1}+3)/432.\end{array}$ (3.1)

两步插值细分的具体步骤是对初始控制点集第一步进行格式(3.1)细分后，第二步再次进行细分格式(3.1)的细分。将格式(3.1)的两次加细作用看作一次细分并做平移，可得到一个九重的细分格式，其生成函数为

$\begin{array}{c}{z}^{-2}a\left(z\right)a\left(z^3\right)=z^{-26}{\left(z^2+z+1\right)}^3{\left(z^6+z^3+1\right)}^3(\left(8b_{-1}+3\right)z^7-\left(32b_{-1}+15\right)z^6+48b_{-1}{z}^5\\ +\left(36-24b_{-1}\right)z^4+\left(36-24b_{-1}\right)z^3+48b_{-1}{z}^2-\left(32b_{-1}+15\right)z+8b_{-1}\\ +3)(\left(8b_{-1}+3\right)z^{21}-\left(32b_{-1}+15\right)z^{18}+48b_{-1}{z}^{15}+\left(36-24b_{-1}\right)z^{12}\\ +\left(36-24b_{-1}\right)z^9+48b_{-1}{z}^6-\left(32b_{-1}+15\right)z^3+8b_{-1}+3)/186624.\end{array}$ (3.2)

计算该九重细分格式(3.2)的加细规则为$p_{9i}^{k+1}=p_i^k$ ，因此该九重细分格式(3.2)是一个插值型细分格式，格式(3.1)是一个两步插值型细分格式。由于细分格式(3.2)是细分格式(3.1)的两步格式，因此两者的整数阶光滑度、Hӧlder指数和基函数均一致。

3.2. 连续性分析

下面定理给出细分格式(3.1)和(3.2)的整数阶的连续性分析。

定理3.1 当参数$b_{-1}$ 满足$-10.2816<b_{-1}<10.372$ 时，细分格式(3.1)和格式(3.2)$C^0$ 连续；当参数$b_{-1}$ 满足$-2.34233<b_{-1}<2.7915$ 时，细分格式(3.1)和格式(3.2)$C^1$ 连续；当参数$b_{-1}$ 满足 $-0.469119<b_{-1}<0.483101$ 时，细分格式(3.1)和格式(3.2)$C^2$ 连续。

证明 由引理2.2知细分格式(3.1)若$C^n$ 连续，则有正整数$L$ 使得${\Vert {\left(\frac{1}{m}{S}_{n+1}\right)}^L\Vert }_{\infty }<1$ 。为了统一，下面分析中对三重细分格式(3.1)采用引理2.2中$m=3,L=2$ 的情形，其恰好对应九重细分格式(3.2)采用引理2.1中$m=9,L=1$ 的情形。两种情形下的细分格式光滑性情况一致。

现对细分格式(3.1)的一阶差分细分序列$S_1$ 的${\left(\frac{1}{3}{S}_1\right)}^2$ 做分组求和，可得范数为

$\begin{array}{l}{\Vert {\left(\frac{1}{3}{S}_1\right)}^2\Vert }_{\infty }=\text{max}\{({\left|8b_{-1}+3\right|}^2+\left|-128b_{-1}^2+2952b_{-1}+1053\right|+2(\left|-32b_{-1}^2+1308b_{-1}+747\right|\\ +\left|-32b_{-1}^2+1692b_{-1}+4923\right|+6\left|-16b_{-1}^2+178b_{-1}+1119\right|)/186624,\\ (4\left|63-52b_{-1}\right|+20\left|4b_{-1}+9\right|+\left|-16b_{-1}^2+130b_{-1}+57\right|\end{array}$

$\begin{array}{l}+2\left|16b_{-1}^2-122b_{-1}+1053\right|+\left|16b_{-1}^2+14b_{-1}+3\right|)/15552,\\ \left(7776\left|\frac{4b_{-1}}{243}+\frac{1}{36}\right|+\left|-32b_{-1}^2+46b_{-1}+1053\right|+\left|32b_{-1}^2+82b_{-1}+27\right|\right)/7776,\\ (\left|-64b_{-1}^2+144b_{-1}+63\right|+2\left|-32b_{-1}^2-516b_{-1}+99\right|+12\left|-16b_{-1}^2-158b_{-1}+1959\right|\\ +2\left|32b_{-1}^2+900b_{-1}-1683\right|+\left|128b_{-1}^2+1272b_{-1}+459\right|)/186624,\\ (12\left|32b_{-1}^2+178b_{-1}+1539\right|+2\left|64b_{-1}^2-204b_{-1}+2511\right|+2\left|64b_{-1}^2+564b_{-1}+1215\right|\\ +\left|128b_{-1}^2+120b_{-1}+27\right|+\left|256b_{-1}^2+480b_{-1}+297\right|)/186624\}.\end{array}$

令不等式${\Vert {\left(\frac{1}{3}{S}_1\right)}^2\Vert }_{\infty }<1$ 可得$-10.2816<b_{-1}<\mathrm{10.372.}$ 此时细分格式(3.1)和格式(3.2)$C^0$ 连续，细分格式收敛。

同理对细分格式(3.1)的二阶差分细分序列$S_2$ 的${\left(\frac{1}{3}{S}_2\right)}^2$ 分组求和，可得范数为

$\begin{array}{l}{\Vert {\left(\frac{1}{3}{S}_2\right)}^2\Vert }_{\infty }=\max \{\left({\left|8b_{-1}+3\right|}^2+48\left|9-8b_{-1}\right|+\left|-64b_{-1}^2-432b_{-1}+1575\right|\right)/10368,\\ \left|-320b_{-1}^2+672b_{-1}+369\right|/20736+\left|-\frac{5b_{-1}^2}{324}+\frac{b_{-1}}{72}+\frac{25}{2304}\right|+\left|-\frac{b_{-1}^2}{81}-\frac{b_{-1}}{864}+\frac{1}{768}\right|\\ +\left|\frac{b_{-1}^2}{81}-\frac{5b_{-1}}{288}+\frac{79}{768}\right|,(\left|-8b_{-1}^2+8b_{-1}+5\right|+\left|-8b_{-1}^2+16b_{-1}+11\right|+\left|8b_{-1}^2-40b_{-1}+11\right|\\ +\left|8b_{-1}^2+16b_{-1}+5\right|)/288,(\left|-32b_{-1}^2+206b_{-1}+63\right|+\left|16b_{-1}^2+14b_{-1}+3\right|\\ +\left|16b_{-1}^2+110b_{-1}+339\right|+\left|32b_{-1}^2-110b_{-1}+81\right|)/1728,(\left|-16b_{-1}^2-94b_{-1}+339\right|\\ +\left|-16b_{-1}^2+2b_{-1}+3\right|+\left|32b_{-1}^2+146b_{-1}+63\right|+\left|32b_{-1}^2+242b_{-1}-81\right|)/1728\}.\end{array}$

令不等式${\Vert {\left(\frac{1}{3}{S}_2\right)}^2\Vert }_{\infty }<1$ 可得$-2.34233<b_{-1}<2.7915$ 。此时细分格式(3.1)和格式(3.2)$C^1$ 连续。

对细分格式(3.1)的三阶差分细分序列$S_3$ 的${\left(\frac{1}{3}{S}_3\right)}^2$ 分组求和，可得范数为

$\begin{array}{l}{\Vert {\left(\frac{1}{3}{S}_3\right)}^2\Vert }_{\infty }=\text{max}\{\frac{1}{2}\left(\left|\left(3-2b_{-1}\right)b_{-1}\right|+\left|b_{-1}\left(2b_{-1}+1\right)\right|\right),\frac{1}{48}(\left|\left(15-16b_{-1}\right)b_{-1}\right|+12\left|b_{-1}\left(2b_{-1}+3\right)\right|\\ +\left|b_{-1}\left(8b_{-1}+3\right)\right|),\frac{{\left|8b_{-1}+3\right|}^2}{2304}+\left|-\frac{5b_{-1}^2}{36}+\frac{121b_{-1}}{96}+\frac{37}{256}\right|+\left|\frac{b_{-1}^2}{9}+\frac{3b_{-1}}{32}+\frac{5}{256}\right|\\ +\left|\frac{5b_{-1}^2}{18}-\frac{5b_{-1}}{16}+\frac{65}{256}\right|,\frac{1}{64}\left|-32b_{-1}^2-62b_{-1}+39\right|+\left|\frac{7}{256}-\frac{7b_{-1}^2}{36}\right|+\left|\frac{11b_{-1}^2}{36}+\frac{61b_{-1}}{96}+\frac{59}{256}\right|,\\ \left|-\frac{7b_{-1}^2}{18}+\frac{3b_{-1}}{8}+\frac{65}{256}\right|+\left|\frac{13b_{-1}^2}{36}+\frac{49b_{-1}}{96}+\frac{37}{256}\right|+\frac{3}{64}\left|16b_{-1}^2-22b_{-1}+7\right|\}.\end{array}$

令不等式${\Vert {\left(\frac{1}{3}{S}_3\right)}^2\Vert }_{\infty }<1$ 可得$-0.469119<b_{-1}<0.483101$ 。此时细分格式(3.1)和格式(3.2)$C^2$ 连续。

3.3. Hӧlder指数

本节给出了参数在生成函数$C^2$ 连续范围内的极限基函数的Hӧlder指数图像。相对于定理3.1给出的整数阶光滑性，Hӧlder指数可对光滑度作出更为精确而细致的刻画。本文采用Matlab ttoolbox [13]计算细分格式的Hӧlder指数。图1给出了Hӧlder指数随参数$b_{-1}$ 从−1.5到1.5按步长0.1计算的结果，呈现出了先增大后减小的趋势。最开始当参数$b_{-1}$ 增大时，Hӧlder指数也相应地增大；当$b_{-1}=0.1$ 时，Hӧlder指数达到最高值2.327734；当参数$b_{-1}$ 超过0.1之后，Hӧlder指数随着参数的增大而减小。

Figure 1. The Hölder exponent curve of the subdivision scheme (3.1) and (3.2)

图1. 细分格式(3.1)和(3.2)的Hӧlder指数曲线

图2给出了当$b_{-1}=-0.2,-0.1,0,0.1,0.2$ 时细分格式(3.1)、(3.2)所对应的基函数图象，可以看出细分格式具有良好的光滑性。

Figure 2. The limit basic function of the subdivision scheme (3.1) and (3.2) when $b_{-1}=-0.2,-0.1,0,1,2$

图2. 当$b_{-1}=-0.2,-0.1,0,1,2$ 时细分格式(3.1)和(3.2)的基函数

3.4. 数值实例

图3系统选择参数$b_{-1}$ 的取值，给出细分格式(3.1)在$b_{-1}=-0.469$ 、$b_{-1}=0.483$ ($C^2$ 连续区间端点值)、$b_{-1}=0$ 、$b_{-1}=0.1$ (Hӧlder指数最大值)以及随机点$b_{-1}=-0.2$ 时曲线细分一次和细分六次的效果图。图3(a)，图3(c)，图3(e)，图3(g)，图3(i)展示了初始控制网格(以蓝色表示)及其经过一次细分后的结果(以红色表示)。图3(b)，图3(d)，图3(f)，图3(h)，图3(j)的图像则展现了初始控制点集经过六次细分加细后的控制网格。由于本文采用了三重细分格式，控制顶点的数量呈现出快速增长的趋势。六次细分后的控制点变得相对密集，曲线变得更加光滑。参数$b_{-1}$ 在$C^2$ 连续区间内细分后的曲线都具有较好的连续性。同时对比$b_{-1}=-0.2$ 和$b_{-1}=0.1$ 两组图像，可以看出$b_{-1}=-0.2$ 的曲线更加光滑。这说明前面的连续性分析以及

(a) (b)

(c) (d)

(e) (f)

(g) (h)

(i) (j)

Figure 3. The subdivision effect diagram of the two-step interpolatory subdivision scheme (3.1); (a) 1 subdivision refinement $\left(b_{-1}=-0.469\right)$ ; (b) 6 subdivision refinements $\left(b_{-1}=-0.469\right)$ ; (c) 1 subdivision refinement $\left(b_{-1}=0.483\right)$ ; (d) 6 subdivision refinements $\left(b_{-1}=0.483\right)$ ; (e) 1 subdivision refinement $\left(b_{-1}=0\right)$ ; (f) 6 subdivision refinements $\left(b_{-1}=0\right)$ ; (g) 1 subdivision refinement $\left(b_{-1}=-0.2\right)$ ; (h) 6 subdivision refinements $\left(b_{-1}=-0.2\right)$ ; (i) 1 subdivision refinement $\left(b_{-1}=0.1\right)$ ; (j) 6 subdivision refinements $\left(b_{-1}=0.1\right)$

图3. 两步插值细分格式(3.1)的细分效果图；(a) 细分$\left(b_{-1}=-0.469\right)$ 1次；(b) 细分$\left(b_{-1}=-0.469\right)$ 6次；(c) 细分$\left(b_{-1}=0.483\right)$ 1次；(d) 细分$\left(b_{-1}=0.483\right)$ 6次；(e) 细分$\left(b_{-1}=0\right)$ 1次；(f) 细分$\left(b_{-1}=0\right)$ 6次；(g) 细分$\left(b_{-1}=-0.2\right)$ 1次；(h) 细分$\left(b_{-1}=-0.2\right)$ 6次；(i) 细分$\left(b_{-1}=0.1\right)$ 1次；(j) 细分$\left(b_{-1}=0.1\right)$ 6次

Hӧlder指数分析是正确的。此外，观察细分后的图像，可以发现原始控制网格在视觉上完全落在极限曲线上。这表明，细分格式(3.1)可以达到较好的插值效果。

NOTES

^*通讯作者。