Toric曲面的渐进迭代逼近

Progressive Iterative Approximation of Toric Surfaces

摘要: 渐进迭代逼近(PIA)是一种直观有效的数据拟合方法。当给定数据点的参数域为不规则的凸多边形时，需要对参数域剖分来用多片曲面拟合，然后考虑相邻曲面片的拼接。Toric曲面是Bézier曲面的推广，它的参数域可以调整为任意凸多边形。使用Toric曲面做渐进迭代逼近即可以保留渐进迭代逼近的优点，也可以整体对数据点进行拟合，无需考虑曲面的重构与拼接。本篇文章定义了一种对凸多边形上的参数点进行字典排序的方法。并实现了一种用Toric曲面做渐进迭代逼近的算法。我们还用具体的数值例子证明方法有效。

Abstract: Progressive iterative approximation (PIA) is an intuitive and effective data fitting method. When the parameter domain of a given data point is an irregular convex polygon, the parameter domain needs to be partitioned to be fitted by a multi-piece surface. Then we consider the stitching of adja-cent surface patches. Toric surfaces are a generalization of Bézier surfaces whose parametric do-main can be adjusted to any convex polygon. Using Toric surface for progressive iterative approxi-mation can not only retain the advantages of progressive iterative approximation, but also fit the data points as a whole, without considering the reconstruction and splicing of the surface. This pa-per defines a lexicographic method for sorting the points of a convex polygon. A progressive itera-tive approximation algorithm using Toric surfaces is also implemented. We also use specific nu-merical examples to prove that the method is effective.

1. 引言

数据拟合是处理科学和工程问题中大量采样和实验数据的常用方法。它通过将离散的数据拟合为连续的代数表示来反映数据的基本趋势。渐进迭代逼近(Progressive Iterative Approximation, PIA)是一种有效的数据拟合方法，它使用几何迭代法动态求解，具有明显的几何意义。渐进迭代逼近方法被广泛应用于几何设计、数据拟合、网格和样条体的生成、逆向工程以及图像处理等领域。它的优点包括直观操作、简单算法和强适应性。

曲线曲面渐进迭代逼近性质的研究自上世纪70年代开始。齐东旭等 [1] 提出了均匀三次B样条曲线盈亏修正算法，并发现这种曲线具有PIA性质。Lin等 [2] 证明非均的三次B样条曲线曲面也具有PIA特性。2005年，Lin等 [3] 证明了在非退化配置矩阵的情况下，曲线(包括张量积曲面)具有渐进迭代逼近性质。张莉 [4] 将PIA的适用范围推广到三角域上，给出了三角片面Bézier曲面具有PIA性质的条件。胡倩倩 [5] 等人对PIA算法进行了改进实现了三角B-B曲面加速逼近的算法，使算法更加高效。Liu [6] 等人研究了张量积Bézier曲面的PIA性质。近年来PIA算法在各种曲线曲面逼近上得到了广泛的应用。季康松 [7] 等给PIA加入了法相约束，使隐式曲线能够更好地拟合散乱数据点及其几何特征。吴硕琳 [8] 等在PIA算法的基础上提出了非均匀三次B样条曲线Hermite插值算法，并证明了该算法是收敛的。周晨 [9] 等实现了用非均匀3次B样条拟合曲线的PIA性质实现了一种矢量地图曲线化简方法，展现出PIA算法在实际生活中有广泛的应用。

Toric曲面是张量积Bézier曲面和三角Bézier曲面片的推广，具有许多与Bézier曲面相类似的良好性质。2002年Ktasauskas [10] 基于Toric簇和Torie理想，定义了一类与整数点集下的有理多边形参数曲面Toric曲面。Toric曲面与Bézier曲面相比，在曲面拼接 [11] 、曲面补洞 [12] 、过渡曲面 [13] 等应用中，Toric方法需要更少的曲面片，整体性更好。汪涵 [14] 等在Toric-Bézier曲线上加入了伸缩因子实现曲线的自由变形，提高了曲线的可调性和预见性。孙兰银等 [15] 利用最小二乘法实现了Toric曲面拟合数据点。2020年，Li [16] 等定义了一种特殊的Toric曲线并将其命名为GT-Bézier曲线，这种曲线是Bézier曲线的拓展。Yu [17] 发现了一种特殊的Toric曲线，并证明了这种曲线具有PIA性。

本文在现有渐进迭代逼近方法的研究基础上，研究Toric曲面的渐进迭代逼近算法。得出当普半径小于1时，Toric曲面同样可以渐进迭代逼近数据点。定义了一种字典排序的方法，将选取的参数点按字典排序，将Toric曲面做渐进迭代逼近使其逐渐靠近这些数据点，最终得到了拟合这些数据点的Toric曲面。本文的第一章介绍了Toric曲面的定义和字典排序方法为下文做准备。第二章讲解了Toric曲面做渐进迭代逼近的过程。第三章给出了数值实例可以直观地看出拟合效果。相较于别的曲面渐进迭代逼近，Toric曲面的参数域选取更加灵活，可以是任意凸多边形。对于不规则形状的参数域，可以直接用Toric曲面做整体拟合，无需考虑曲面的重构与拼接。

2. 预备知识

定义1. (Toric-Bernstein基函数)：设在平面中，坐标为整数的有限个点构成的集合为 $A\subset {\mathbb{Z}}^2$ 。设该点集A的凸包 $conv\left(A\right)$ 为 $\Delta A$ 。令n边形 $\Delta A$ 的第r条边的边界方程为 $l_r\left(u,v\right):a_{r}u+b_{r}v+c_r=0$ ， $r=1,2,\cdots ,n$ 。这里的 $a_r,b_r,c_r$ 互素且向量 $\langle a_r,b_r\rangle$ 指向 $\Delta A$ 的内部。那么对于集合A中的每一个整数格点 $\left(i,j\right)\in A$ ，定义Toric-Bernstein基函数为：

$B_{i,j}\left(u,v\right)=C_{i,j}{l}_1{\left(u,v\right)}^{l_1\left(i,j\right)}{l}_2{\left(u,v\right)}^{l_2\left(i,j\right)}\cdots l_n{\left(u,v\right)}^{l_n\left(i,j\right)}\left(u,v\right)\in \Delta A$ (1)

其中 $C_{i,j}>0$ 为自定义系数。

定义2. (Toric曲面)：对于有限整数点集A，定义在 $\Delta A$ 上的Toric Bézier曲面为 $F_{A,\omega ,B}\left(u,v\right)$ 为：

$F_{A,\omega ,B}\left(u,v\right)=\frac{{{\displaystyle \sum }}_{\left(i,j\right)\in A}{\omega }_{i,j}{p}_{i,j}{B}_{i,j}\left(u,v\right)}{{{\displaystyle \sum }}_{\left(i,j\right)\in A}{\omega }_{i,j}{B}_{i,j}\left(u,v\right)}\left(u,v\right)\in \Delta A$ (2)

其中 $B_{i,j}\left(u,v\right)$ 是Toric-Bernstein基函数， $p_{i,j}\in {ℝ}^3$ 为对应整数格点 $\left(i,j\right)\in A$ 的控制顶点， ${\omega }_{i,j}$ 为控制顶点的权因子。

因为Toric曲面定义在 $\Delta A$ 上， $\Delta A$ 可能是一个不规则图形。为了描述方便需要给这些整数点集进行排序，下面给出字典排序的定义：

定义3. (字典排序)：设平面 ${\mathbb{Z}}^2$ 中有两个坐标为整数的点 $\left(i,j\right)$ 和 $\left(i^{\prime },j^{\prime }\right)$ ，若 $\left(i,j\right)>\left(i^{\prime },j^{\prime }\right)$ 成立，当且仅当 $i>i^{\prime }$ 或( $i=i^{\prime }$ 且 $j>j^{\prime }$ )。

设点集A有n个坐标为整数的点，可以用字典排序的方法从小到大将点集A重新排列为 $A=\left\{a_1,a_2,\cdots ,a_n\right\}$ 。那么Toric-Bernstein基函数可以写成，

$B_{a_k}\left(u,v\right)=C_{a_k}{l}_1{\left(u,v\right)}^{l_1\left(a_k\right)}{l}_2{\left(u,v\right)}^{l_2\left(a_k\right)}\cdots l_n{\left(u,v\right)}^{l_n\left(a_k\right)}\left(u,v\right)\in \Delta A,k=1,2,\cdots ,n$

其中 $C_{a_k}$ 为自定义系数。同理，Toric曲面可以写成，

$F_{A,\omega ,B}\left(u,v\right)=\frac{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{p}_{a_k}{B}_{a_k}\left(u,v\right)}{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{B}_{a_k}\left(u,v\right)}\left(u,v\right)\in \Delta A$

3. Toric曲面渐进迭代逼近

已知空间 ${ℝ}^3$ 中有n个数据点的集合 $\left\{p_{i,j}|\left(i,j\right)\in A\right\}$ ，每一个数据点 $p_{i,j}$ 对应的参数值为 $\left(i,j\right)$ ，它们组成平面有限点集A。将点集A中的点进行从小到大的字典排序可以将A写成点列 $A=\left\{a_1,a_2,\cdots ,a_n\right\}$ ，这样集合 $\left\{p_{i,j}|\left(i,j\right)\in A\right\}$ 可以按照角标的字典排序变为 ${\left\{p_{a_k}\right\}}_{k=1}^n$ 。同样的也对点集A中的每一顶点的参数值进行字典排序，设排序后的参数值为集合 ${\left\{t_{a_k}\right\}}_{k=1}^n$ 。通过点集A可以得到凸包 $\Delta A$ ，并得到 $\Delta A$ 对应的边界方程以及A中任意一点对应的Toric-Bernstein 基函数。通过预先自定义好的Toric-Bernstein基系数 $C_{a_k}$ ，以及每个参数点对应的权因子 ${\omega }_{a_k}$ 。设定迭代的初始控制顶点为 ${\left\{p_{a_k}^0=p_{a_k}\right\}}_{k=1}^n$ 。这样我们可以得出初始迭代Toric曲面：

$F_{A,\omega ,B}^0\left(u,v\right)=\frac{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{p}_{a_k}^0{B}_{a_k}\left(u,v\right)}{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{B}_{a_k}\left(u,v\right)}\left(u,v\right)\in \Delta A$

之后计算每个控制顶点需要调整的向量 ${\left\{{\delta }_i\right\}}_{i=0}^n$

${\delta }_{a_k}^0=p_{a_k}-F_{A,\omega ,B}^0\left(t_{a_k}\right)k=1,2,\cdots ,n$

我们可以得到第二次迭代Toric曲面的控制顶点

$p_{a_k}^1=p_{a_k}^0+{\delta }_{a_k}^{0}k=1,2,\cdots ,n$

最终得到新的曲面 $F_{A,\omega ,B}^1\left(u,v\right)$

$F_{A,\omega ,B}^1\left(u,v\right)=\frac{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{p}_{a_k}^1{B}_{a_k}\left(u,v\right)}{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{B}_{a_k}\left(u,v\right)}\left(i,j\right)\in A$

同理，我们可以通过曲面 $F_{A,\omega ,B}^m\left(u,v\right)$ 得到曲面 $F_{A,\omega ,B}^{i+1}\left(u,v\right)$

${\delta }_{a_k}^m=p_{a_k}-F_{A,\omega ,B}^m\left(t_{a_k}\right)k=1,2,\cdots ,n$

$p_{a_k}^{m+1}=p_{a_k}^m+{\delta }_{a_k}^{m}k=1,2,\cdots ,n$

$F_{A,\omega ,B}^{m+1}\left(u,v\right)=\frac{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{p}_{a_k}^{k+1}{B}_{a_k}\left(u,v\right)}{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{B}_{a_k}\left(u,v\right)}\left(u,v\right)\in \Delta A$

如果 $\underset{m\to \infty }{\lim }{F}_{A,\omega ,B}^m\left(a_k\right)=p_{a_k},k=1,2,\cdots ,n$ 成立。我们可以称初始曲面 $F_{A,\omega ,B}^m\left(a_k\right)$ 有PIA性。下面我们分析一下渐进迭代逼近收敛的条件。

${\delta }_{a_k}^{m+1}={\delta }_{a_k}^m-\underset{i=1}{\overset{n}{{\displaystyle \sum }}}{\delta }_{a_i}^m\frac{{\omega }_{a_i}{B}_{a_i}\left(t_{a_k}\right)}{{{\displaystyle \sum }}_{j=1}^n{\omega }_{a_j}{B}_{a_j}\left(t_{a_k}\right)}k=1,2,\cdots ,n;m=0,1,\cdots$

上式写成矩阵的形式为

${\left[{\delta }_{a_1}^{m+1},{\delta }_{a_2}^{m+1},\cdots ,{\delta }_{a_n}^{m+1}\right]}^{\text{T}}=\left(I-B\right){\left[{\delta }_{a_1}^m,{\delta }_{a_2}^m,\cdots ,{\delta }_{a_n}^m\right]}^{\text{T}}m=0,1,\cdots$

这里的I是指 $\left(n+1\right)\times \left(n+1\right)$ 的单位矩阵，矩阵B定义为

$B:=\left[\begin{array}{cccc}\frac{{\omega }_{a_1}{B}_{a_1}\left(t_{a_1}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_1}\right)}& \frac{{\omega }_{a_2}{B}_{a_2}\left(t_{a_1}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_1}\right)}& \cdots & \frac{{\omega }_{a_n}{B}_{a_n}\left(t_{a_1}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_1}\right)}\\ \frac{{\omega }_{a_1}{B}_{a_1}\left(t_{a_2}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_2}\right)}& \frac{{\omega }_{a_2}{B}_{a_2}\left(t_{a_2}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_2}\right)}& \cdots & \frac{{\omega }_{a_n}{B}_{a_n}\left(t_{a_2}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_2}\right)}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{{\omega }_{a_1}{B}_{a_1}\left(t_{a_n}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_n}\right)}& \frac{{\omega }_{a_2}{B}_{a_2}\left(t_{a_n}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_n}\right)}& \cdots & \frac{{\omega }_{a_n}{B}_{a_n}\left(t_{a_n}\right)}{{\displaystyle \underset{j=1}{\overset{n}{\sum }}{\omega }_{a_j}}{B}_{a_j}\left(t_{a_n}\right)}\end{array}\right]$ (3)

当矩阵 $\left(I-B\right)$ 的普半径 $\rho \left(I-B\right)<1$ 时渐进迭代收敛，当我们要使用Toric Bézier曲面做渐进迭代逼近时，需要提前计算出 $\rho \left(I-B\right)$ 以确保收敛。

Toric曲面的渐进迭代逼近算法：

输入：空间 ${ℝ}^3$ 中的n个数据点 ${\left\{p_{a_k}\right\}}_{k=1}^n$ ，每个数据点对应的权因子 ${\left\{{\omega }_{a_k}\right\}}_{k=1}^n$ 和迭代次数 $\hat{m}$ 。

输出：对数据点 ${\left\{p_{a_k}\right\}}_{k=1}^n$ 的m次渐进迭代逼近后的Toric Bézier曲面。

Step 1：根据数据点集 ${\left\{p_{a_k}\right\}}_{k=1}^n$ 的确定点集A以及 $\Delta A$ ，并由此得出Toric-Bernstein基函数 ${\left\{B_{a_k}\left(u,v\right)\right\}}_{k=1}^n$ ，代入到矩阵(3)中求出普半径 $\rho \left(I-B\right)$

Step 2：当普半径 $\rho \left(I-B\right)<1$ 时：

Step 2.1：令迭代次数 $m=0$ 并计算初始Toric Bézier曲面控制顶点 ${\left\{p_{a_k}^0=p_{a_k}\right\}}_{k=1}^n$

Step 2.2：计算误差 ${\delta }_{a_k}^m=p_{a_k}-F_{A,\omega ,B}^m\left(a_k\right)$ ，和新的控制顶点 $p_{a_k}^{m+1}=p_{a_k}^m+{\delta }_{a_k}^m$

Step 2.3：若 $m<\hat{m}$ ，令 $m=m+1$ 回到第Step 2.2；若 $m=\hat{m}$ ，输出Toric Bézier曲面

$F_{A,\omega ,B}^{m+1}\left(u,v\right)=\frac{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{p}_{a_k}^{k+1}{B}_{a_k}\left(u,v\right)}{{{\displaystyle \sum }}_{k=1}^n{\omega }_{a_k}{B}_{a_k}\left(u,v\right)}\left(u,v\right)\in \Delta A$

Step 3：当普半径 $\rho \left(I-B\right)\ge 1$ 时：重新调整权因子 ${\left\{{\omega }_{a_k}\right\}}_{k=1}^n$ 或取样点 ${\left\{p_{a_k}\right\}}_{k=1}^n$

4. 数值实例

例1：设空间 ${ℝ}^3$ 中有8个数据点按词典顺序为 $p_{0,0}=\left(0,0,1\right)$ ， $p_{0,1}=\left(0,1,2\right)$ ， $p_{0,2}=\left(0,2,1\right)$ ， $p_{1,0}=\left(1,0,2\right)$ ， $p_{1,1}=\left(1,1,3\right)$ ， $p_{1,2}=\left(1,2,1\right)$ ， $p_{2,0}=\left(2,0,1\right)$ ， $p_{2,1}=\left(2,1,1\right)$ ，视每一点的权因子都为1。可得 $A=\left\{a_1,a_2,\cdots ,a_8\right\}=\left\{\left(0,0\right),\left(0,1\right),\cdots ,\left(2,1\right)\right\}$ 。 $\Delta A$ 如图1所示。计算出toric-Bernstein基函数 ${\left\{B_{a_k}\left(u,v\right)\right\}}_{k=1}^n$ ，并求出普半径 $\rho \left(I-B\right)=0.875$ 证明在此基函数下Toric Bézier曲面渐进迭代逼近收敛。逼近效果如图2~4所示：

Figure 1. $\mathcal{A}$ and $\Delta \mathcal{A}$

图1. $\mathcal{A}$ 和 $\Delta \mathcal{A}$

Figure 2. The initial Toric surface

图2. 初始Toric曲面

Figure 3. The surface after 1 iteration

图3. 迭代1次后的曲面

Figure 4. The surface after 50 iteration

图4. 迭代50次后的曲面

例2：如图5所示的凸图形 $A=\left\{a_1,a_2,\cdots ,a_{19}\right\}=\left\{\left(0,0\right),\left(0,1\right),\cdots ,\left(4,4\right)\right\}$ 中取19个数据采集点，并在函数上 $f\left(u,v\right)=u^2+v^2$ 上找到每个数据采集点对应的数值大小如图6，并将他们按照字典排序，不失一般性可设每个拟合点对应的权因子都为1。求出普半径 $\rho \left(I-B\right)=0.966$ 证明迭代是收敛的。迭代次数与渐进逼近效果如图7所示。

Figure 5. $\mathcal{A}$ and $\Delta \mathcal{A}$

图5. $\mathcal{A}$ 和 $\Delta \mathcal{A}$

Figure 6. $f\left(u,v\right)$ and data point ${\left\{p_{a_k}\right\}}_{k=1}^n$

图6. $f\left(u,v\right)$ 和数据点 ${\left\{p_{a_k}\right\}}_{k=1}^n$

Figure 7. Iteration effect

图7. 迭代效果

例3：如(图8)所示的凸图形中取72个数据点以及每个点所对应权重大小，适当的调节对应点权重的

大小可以有更好的拟合效果。并用这些点来拟合平面 $y=\sin \left(\frac{u}{2}\right)\text{con}\left(\frac{v}{2}\right)$ (如图9)中的一部分，其拟合效

果如图10所示。

Figure 8. Sampling points and weight distribution

图8. 取样点和权重分布

Figure 9. $f\left(u,v\right)$ and data point ${\left\{p_{a_k}\right\}}_{k=1}^n$

图9. $f\left(u,v\right)$ 和数据点 ${\left\{p_{a_k}\right\}}_{k=1}^n$

Figure 10. Iteration effect

图10. 迭代效果

5. 总结与展望

本篇文章定义了一种新的对凸多边形上的参数点进行字典排序的方法，方便对不规则凸多边形上的点进行排序与编号。并实现了一种用Toric曲面做渐进迭代逼近的算法。该方法可以直接逼近不规则的凸多边形域上的数据点，而无需考虑曲面的重构与拼接。相较于张量积Bézier曲面 [6] 和三角片面Bézier [4] [5] 曲面的渐进迭代逼近算法，本文的算法更加灵活，可以为任意凸多边形，所以适用场景更多。

总结研究成果，发现仍有以下问题可以进一步探讨：

1) Toric曲面的渐进迭代逼近收敛性是受toric-Bernstein基函数(公式1)中的自定义系数 $C_{i,j}$ 和Toric Bézier曲面(公式2)的权因子 ${\omega }_{i,j}$ 影响的。给出Toric曲面渐进迭代逼近收敛的约束条件是我们进一步考虑的问题。

2) 矩阵B (矩阵3)为全正矩阵是渐进迭代收敛充分条件。如何在一定约束条件下证明(矩阵3)为完全正矩阵也是之后值得研究的问题之一。

NOTES

^*通讯作者。