细分曲面奇异点处的光滑过渡

Smooth Connection near Singular Points on Subdivision Surfaces

摘要: 首先通过引入形状调节参数c(0≤c≤1)给出了一种形状可调的C-C细分方法，基于该理论，针对细分曲面中的奇异点，本文以奇异点处的2-环作为控制网格，采用循环映射的方法，提出了一种形状可调的G²曲面造型方法，得到Bézier控制点的显式解。与以往的方法相比，生成的曲面不仅在奇异点处达到G²连续，而且解决了曲面设计的可调性问题。本文给出了算法流程和相关数据结构，也给出了相应的实例进行验证。

Abstract: This paper discusses G² smooth connection near singular points on subdivision surfaces. A shape adjustable Catmull-Clark subdivision algorithm is present by introducing subdivision shape ad-justment parameter c(0≤c≤1). Based on this, aimed at the singular points in subdivision surface, based on the 2-ring of the singular points as control mesh, adopting the method of cyclic mapping, this paper proposes a G² surface modeling method which is shape-adjustable. Then we can get the explicit solution of the Bezier control points. Compared with existing methods, the generated surface not only achieves the G² continuous at singular points, but solves the problem of the curved surface design to be adjustable. The algorithm process and data structure are present. Corresponding examples are also given in this paper.

1. 引言

细分曲面与NURBS曲面是现代工业的两大造型技术，细分曲面克服了NURBS曲面处理任意拓扑结构存在的困难，成为有效的建模手段之一。在曲面造型技术中奇异点的出现始终是无法避免的，由于奇异点的存在，最终得到的曲面会出现“褶皱”等现象，严重影响曲面的质量，尤其是在工业产品的设计中，如果奇异点的处理不恰当，会使产品无法投入使用，因此，奇异点的处理技术成为曲面造型的关键点。传统的细分方法的奇异点处仅仅满足 $G^1$ 连续，为了构造出 $G^2$ 连续的曲面，人们进行了很多的探索。2009年，剑桥大学的Thomas J. Cashman [1] 等人采用任意次数的B样条曲线曲面给出了一种对称的非均匀细分算法，并且将该算法推广到曲面，得到了一种奇数次的对称非均匀细分算法，但是该算法在奇异点处没有具体说明，而且只针对奇数阶的情况；后来，Thomas J. Cashman和Malcolm A. Sabin等人 [2] 设计了新的细分格式来处理非均匀细分算法中的奇异点，但是结果并不理想，具有较差的局部性；刘浩 [3] 等人利用流行方法和细分方法构造出 $G^2$ 连续的自由型曲面，实现了流行曲面和Catmull-Clark细分曲面的融合，然而，在构造 $n$ 边域曲面片时没有考虑奇异点的位置信息。本文在Charles Loop [4] 等人研究的基础上，进行了以下改进：(1) 引入了一种形状可调的Catmull-Clark细分方法，将该方法作为C-C细分的前置方法，对曲面进行“混合细分”，使最终的拟合曲面具有一定的可调性；(2) 采用循环映射对奇异点周围的曲面片进行处理时，得到的循环矩阵也含有参数 $c$ ，使每一个曲面片都具有可调性；(3) 给出了双七次Bézier控制点的具体位置以及循环方程组的显式解，得到每一个曲面片的显示表达式，并给出了奇异点处的 $G^2$ 填充结果。

2. 一种形状可调的C-C细分曲面

经典的Catmull-Clark细分基于双三次均匀B-样条，在给定初始控制网格的情况下，生成的结果曲面也就唯一确定，不具有可调性。而在实际应用中，往往大多时候生成的结果曲面不能直接满足要求，需要进行一定的调整。因此，为了增强曲面设计的灵活性和可调性，我们将曲线的可调性 [5] 推广到曲面的可调性，重新设计几何规则和拓扑规则。

对于任意拓扑形状的控制网格，新的细分规则还是由两部分组成：几何点的产生和拓扑结构的建立。

(1) 几何点的产生

给定一个初始网格 ${\Gamma }^{\left(0\right)}$ ，经过 $k$ 次改进的Catmull-Clark细分后得到的控制网格为 ${\Gamma }^{\left(k\right)}$ 。考虑 ${\Gamma }^{\left(k\right)}$ 上度为 $n$ 的顶点 $v^{\left(k\right)}$ 的局部网格，它的邻域由 $n$ 个共边顶点 $e_1^{\left(k\right)},\cdots ,e_n^{\left(k\right)}$ 和 $n$ 个共面顶点 $f_1^{\left(k\right)},\cdots ,f_n^{\left(k\right)}$ 组成，则第 $k+1$ 次细分后产生的新面点、新边点和新顶点为：

新面点(F-点)：用面上所有顶点的均值表示

$f_i^{\left(k+1\right)}=\frac{1}{4}(v^{\left(k\right)}+f_i^{\left(k\right)}+e_i^{\left(k\right)}+e_{i+1}^{(k)})$

新边点(E-点)：用相邻两面上的新面点均值的 $c$ 倍与对应边两端点的均值的 $\left(1-c\right)$ 倍的和来表示

$e_i^{\left(k+1\right)}=\frac{v^{\left(k\right)}+e_i^{\left(k\right)}}{2}c+\frac{f_{i-1}^{\left(k+1\right)}+f_{i+1}^{\left(k+1\right)}}{2}\left(1-c\right)$

新顶点(F-点)：用对应顶点、新面点与其周围一些点的线性组合可以表示为

$v^{(k+1)}=\frac{4cQ}{n}+\frac{8c(1-c)R}{n}+\frac{n-4c(2-c)v^{(k)}}{n}$

其中 $Q$ 表示所有与旧顶点相邻的面上新面点的均值， $R$ 表示所有交于旧顶点的边的中点的均值， $n$ 为对应顶点的度数。同时， $i=1,\cdots ,n$ ，且设 $e_{n+1}^{\left(k\right)}=e_1^{\left(k\right)}$ ， $f_{n+1}^{\left(k\right)}=f_1^{\left(k\right)}$ 。

(2) 拓扑结构的建立

与经典Catmull-Clark的拓扑规则相同，连接每一个新面点与其周围的新边点；连接每一个新顶点与其周围的新边点。

从曲面的细分规则可以推出，给定原控制点 $d_{ij}\left(i,j=0,\cdots ,n\right)$ 则细分后的曲面片可表示为 $P\left(u,v\right)=UCG{C}^{\text{T}}{V}^{\text{T}}$ ，其中 $U=\left[\begin{array}{cccc}1& u& u^2& u^3\end{array}\right]$ ， $V=\left[\begin{array}{cccc}1& v& v^2& v^3\end{array}\right]$ ，相应的参数调节矩阵为

$C=\left[\begin{array}{cccc}\frac{c}{2\left(1+c\right)}& \frac{1}{1+c}& \frac{c}{2\left(1+c\right)}& 0\\ -\frac{c\left(2-c\right)}{1+c}& 0& \frac{c\left(2-c\right)}{1+c}& 0\\ \frac{c\left(5-4c\right)}{2\left(1+c\right)}& \frac{\left(-2c+3\right)\left(c-2\right)}{2\left(1+c\right)}& \frac{4c^2-11c+6}{2\left(1+c\right)}& \frac{c\left(2c-1\right)}{2\left(1+c\right)}\\ \frac{c\left(c-1\right)}{1+c}& \frac{\left(c-1\right)\left(c-2\right)}{1+c}& -\frac{\left(c-1\right)\left(c-2\right)}{1+c}& \frac{c\left(1-c\right)}{1+c}\end{array}\right]$

3. 细分曲面奇异点处的G²连续性

本节通过改进Charles Loop [4] 的理论给出了含有奇异点的 $G^2$ 重建方法。首先，对初始控制网格进行“混合细分”，即先采用形状可调的C-C细分方法进行细分然后再用经典的C-C细分，目的是为了保证曲面片具有B样条边界；其次，在“混合细分”的基础上，以奇异点处的2-环作为控制网格，采用循环映射的方法得到了二阶连续的约束方程组；最后，利用循环矩阵和能量函数优化方法给出了 $G^2$ 填充 $n$ 边洞时，Bézier控制点的显式解，不仅使奇异点周围曲面片之间满足 $G^2$ ，并且和Catmull-Clark细分曲面实现 $G^2$ 拼接，得到整体 $G^2$ 连续的曲面。

3.1. 二阶几何连续条件

定义1：两曲面 $P\left(u,v\right)$ 和 $Q\left(s,r\right)$ 沿公共连接线 $P\left(t\right)=Q\left(t\right)$ 具有 $G^2$ 连续的充要条件 [6] 是两曲面之一可以重新参数化，以至于它们沿公共连接线为 $C^k$ 连续的。如 $Q\left(s,r\right)$ 重新参数化为 $\bar{Q}\left(\bar{s},\bar{r}\right)$ ，使得

$\frac{{\partial }^{i+j}}{\partial u^i\partial v^j}P\left(t\right)=\frac{{\partial }^{i+j}}{\partial s^i\partial r^j}\bar{Q}\left(t\right),i+j=0,1,\cdots ,k.$

根据上述定义可知，给定两个曲面片 $P_i\left(u,v\right)$ 和 $P_{i\text{+}1}\left(u,v\right)$ ： $\left[0,1\right]\times \left[0,1\right]\to R^m$ ，沿公共边界 $C:C\left(v\right)=P_i\left(0,v\right)=P_{i+1}\left(0,v\right)$ 二阶几何连续的充要条件是存在函数 $\left\{\left(p_i\left(v\right),q_i\left(v\right)\right),i=1,2\right\}$ 使得下列表达式沿 $C$ 成立：

$\frac{{\partial }^2}{\partial u_{}^2}{P}_i^{}=p_2\left(v\right)\frac{\partial }{\partial u}{P}_{i+1}+q_2\left(v\right)\frac{\partial }{\partial v}{P}_{i+1}+p_1^2\left(v\right)\frac{{\partial }^2}{\partial u^2}{P}_{i+1}+2p_1\left(v\right)q_1\left(v\right)\frac{{\partial }^2}{\partial u\partial v}{P}_{i+1}+q_1^2\left(v\right)\frac{{\partial }^2}{\partial v^2}{P}_{i+1}.$

3.2. 初始控制网格的“混合细分”

对于任意拓扑形状的散乱点集得到的控制网格，采用本文第1节给出的形状可调的C-C细分算法，对网格进行细分，得到的网格曲面可以表示为

$P\left(u,v\right)=UCH{C}^{\text{T}}{V}^{\text{T}}$

令 $A=CH{C}^{\text{T}}$ ， $H$ 为双三次系数矩阵，则第 $k$ 个曲面片对应的系数矩阵为 $A_k=C_k{H}_k{C}_k^{\text{T}}$ ，相应的参数调节矩阵为

$C_k=\left[\begin{array}{cccc}\frac{c_k}{2\left(1+c_k\right)}& \frac{1}{1+c_k}& \frac{c_k}{2\left(1+c_k\right)}& 0\\ -\frac{c_k\left(2-c_k\right)}{1+c_k}& 0& \frac{c_k\left(2-c_k\right)}{1+c_k}& 0\\ \frac{c_k\left(5-4c_k\right)}{2\left(1+c_k\right)}& \frac{\left(-2c_k+3\right)\left(c_k-2\right)}{2\left(1+c_k\right)}& \frac{4c_k^2-11c_k+6}{2\left(1+c_k\right)}& \frac{c_k\left(2c_k-1\right)}{2\left(1+c_k\right)}\\ \frac{c_k(c_k-1)}{1+c_k}& \frac{\left(c_k-1\right)\left(c_k-2\right)}{1+c_k}& -\frac{\left(c_k-1\right)\left(c_k-2\right)}{1+c_k}& \frac{c_k\left(1-c_k\right)}{1+c_k}\end{array}\right]$

为了保持双三次样条边界，必要时可以在此基础上重新进行Catmull-Clark细分，该过程称为“混合细分”。通过“混合细分”将奇异点从正则网格中隔离出来，构造出奇异点的 $n$ 边洞，对于正则处已经满足二阶连续，下面要做的就是如何在保证二阶几何连续的情况下填充 $n$ 边洞并且与周围的曲面片拼接后曲面整体达到二阶连续。如图1所示。

在图1中，与奇异点直接相邻的边，称为内部边界；与奇异点间接相邻的边称为外部边界。要使最后整张曲面为 $G^2$ 连续，要求曲面片内部边界之间的拼接是 $G^2$ 且外部边界与Catmull-Clark曲面片之间的拼接也满足 $G^2$ 。经过“混合细分”后，构造出奇异点的二环，如下图2所示。

对于 $n$ 边洞的填充，本文构造了双七次B-样条曲面片进行填充，令

$T_i=\left[\begin{array}{cccc}\frac{\partial P_i}{\partial u}& \frac{\partial P_i}{\partial v}& \frac{{\partial }^2{P}_i}{\partial u^2}& \begin{array}{cccc}\frac{{\partial }^2{P}_i}{\partial u\partial v}& \frac{{\partial }^2{P}_i}{\partial v^2}& \frac{{\partial }^3{P}_i}{\partial u^2\partial v}& \begin{array}{cc}\frac{{\partial }^3{P}_i}{\partial u\partial v^2}& \frac{{\partial }^4{P}_i}{\partial u^2\partial v^2}\end{array}\end{array}\end{array}\right]$

则

$T_{i+1}=\left[\begin{array}{cccc}\frac{\partial P_{i+1}}{\partial u}& \frac{\partial P_{i+1}}{\partial v}& \frac{{\partial }^2{P}_{i+1}}{\partial u^2}& \begin{array}{cccc}\frac{{\partial }^2{P}_{i+1}}{\partial u\partial v}& \frac{{\partial }^2{P}_{i+1}}{\partial v^2}& \frac{{\partial }^3{P}_{i+1}}{\partial u^2\partial v}& \begin{array}{cc}\frac{{\partial }^3{P}_{i+1}}{\partial u\partial v^2}& \frac{{\partial }^4{P}_{i+1}}{\partial u^2\partial v^2}\end{array}\end{array}\end{array}\right]$

若 ${\Theta }_i$ 为曲面片 $P_i$ 和 $P_{i+1}$ 之间的参数化矩阵，则由二阶几何连续性定义可以得到

Figure 1. N-sided hole formed by mingle subdivision

图1. 混合细分产生的n边洞

Figure 2. 2-ring around singular vertex on mesh

图2. 混合细分后奇异点的二环

$T_i=T_{i+1}\cdot {\Theta }_i$

那么，具有公共顶点的 $n$ 个曲面片 $P_i$ ( $i=0,\cdots ,n-1$ )，相邻的曲面片之间满足几何连续和循环条件，则 $n$ 个具有公共顶点的曲面片之间满足如下关系

$T_0=T_0\cdot {\Theta }_{n-1}\cdot {\Theta }_{n-2}\cdots {\Theta }_1\cdot {\Theta }_0$

从上式推得

$I={\Theta }_{n-1}\cdot {\Theta }_{n-2}\cdots {\Theta }_1\cdot {\Theta }_0$

要使最后整张曲面为 $G^2$ 连续，要求曲面片内部边界之间的拼接是 $G^2$ 且外部边界与Catmull-Clark曲面片之间的拼接也满足 $G^2$ 。

3.2.1.内部边界的连续性

首先，定义内部边界的相关映射为 ${\phi }_n:\left(u,v\right)\to \left(x,y\right)$ ，令

$\{\begin{array}{l}x=\left(\cos \frac{\text{2π}}{n}\right)\left[\left(1-u\right)v+\sec \frac{\text{2π}}{n}\right]\\ y=\left(\sin \frac{\text{2π}}{n}\right)\left[\left(1-u\right)v+\left(\sec \frac{\text{2π}}{n}\right)uv\right]\end{array}$

$r_n$ 是以原点为中心 $\frac{\text{2π}}{n}$ 为单位的一个逆时针旋转变换

$r_n\left(u,v\right)=\left[\begin{array}{cc}\cos \frac{\text{2π}}{n}& -\sin \frac{\text{2π}}{n}\\ \sin \frac{\text{2π}}{n}& \cos \frac{\text{2π}}{n}\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]$

则可以得到 $P_{i+1}={\phi }_n^{-1}\circ \left(r_n^{-1}\circ {\phi }_n\right)P_i$ 。于是，对于 $n$ 价奇异点周围的曲面片 $P_0,P_1,\cdots ,P_{n-2},P_{n-1}$ 有

$I={\left({\Phi }_n^{-1}\cdot R_n^{-1}\cdot {\Phi }_n\right)}^n={\Phi }_n^{-1}\cdot {\left(R_n^{-1}\right)}^n\cdot {\Phi }_n$ (1)

成立， ${\Phi }_n$ ， $R_n$ 为对应变换 ${\varphi }_n$ ， $r_n$ 的矩阵，并且 ${\left(R_n^{-1}\right)}^n=I$ 。

对于带有内部边界的曲面片 $P_i\left(0,t\right)$ 和 $P_{i+1}\left(t,0\right)$ 达到二阶几何连续，则需要满足下面三个条件

$P_i\left(0,t\right)-P_{i+1}\left(t,0\right)=0$ (2)

$c_n\left(1-t\right)P_{i+1}^u\left(t,0\right)-P_{i+1}^v\left(t,0\right)-c_n\left(t-1\right)P_i{}^v\left(0,t\right)-P_i^u\left(0,t\right)=0$ (3) $\begin{array}{l}2c_n^2\left(1-t\right)\left(P_{i+1}^v\left(t,0\right)-P_i^u\left(0,t\right)\right)+2c_n\left(1-t\right)\left(1+c_n\left(1-t\right)\right)\left(P_{i+1}^{uv}\left(t,0\right)-P_i^{uv}\left(0,t\right)\right)\\ -\left(1+c_n\left(1-t\right)\right)\left(P_{i+1}^{vv}\left(t,0\right)-P_i^{uu}\left(0,t\right)\right)=0\end{array}$ (4)

令双七次张量积曲面片 $P_s\left(u,v\right)={\displaystyle \underset{i=0}{\overset{7}{\sum }}{\displaystyle \underset{j=0}{\overset{7}{\sum }}{p}_{ij}^s{B}_i^7\left(u\right)B_j^7\left(v\right)}}\left(0\le u,v\le 1\right)$ ，则对于 $n$ 边洞的第 $i$ 和 $i+1$ 个曲面片 $P_i\left(u,v\right)$ 和 $P_{i+1}\left(u,v\right)$ 可以得到 $P_i\left(0,t\right)$ 、 $P_{i+1}\left(t,0\right)$ 、 $P_i^u\left(0,t\right)$ 、 $P_i^v\left(0,t\right)$ 、 $P_i^{uv}\left(0,t\right)$ 、 $P_i^{uu}\left(0,t\right)$ 、 $P_{i+1}^u\left(t,0\right)$ 、 $P_{i+1}^v\left(t,0\right)$ 、 $P_{\text{i+1}}^{uv}\left(t,0\right)$ 、 $P_{i+1}^{vv}\left(t,0\right)$ ，将上述结果代入(2)可以得到

$p_{00}^i=p_{07}^{i+1},p_{01}^i=p_{17}^{i+1},p_{02}^i=p_{27}^{i+1},p_{03}^i=p_{37}^{i+1},$

$p_{04}^i=p_{47}^{i+1},p_{05}^i=p_{57}^{i+1},p_{06}^i=p_{77}^{i+1},p_{07}^i=p_{77}^{i+1}.$

代入(3)可得

$c_n\left(-{\alpha }_1-2{\alpha }_2+3{\alpha }_3-7{\alpha }_4+5{\alpha }_5-21{\alpha }_6+7{\alpha }_7-{\alpha }_8\right)$ $+\left({\beta }_1-{\beta }_2+{\beta }_3-{\beta }_4+{\beta }_5-{\beta }_6+{\beta }_7-{\beta }_8\right)=0$ .

$c_n\left(7{\alpha }_1\text{+}13{\alpha }_2-17{\alpha }_3+31{\alpha }_4-19{\alpha }_5+57{\alpha }_6-13{\alpha }_7+{\alpha }_8\right)$ $+\left(-7{\beta }_1+6{\beta }_2-5{\beta }_3+4{\beta }_4-3{\beta }_5+2{\beta }_6-{\beta }_7\right)=0.$ $c_n\left(-21{\alpha }_1-36{\alpha }_2+40{\alpha }_3-54{\alpha }_4+27{\alpha }_5-51{\alpha }_6+6{\alpha }_7\right)$ $+\left(21{\beta }_1-15{\beta }_2+10{\beta }_3-6{\beta }_4+3{\beta }_5-{\beta }_6\right)=0.$ $c_n\left(35{\alpha }_1+55{\alpha }_2-50{\alpha }_3+46{\alpha }_4-17{\alpha }_5+15{\alpha }_6+6{\alpha }_7\right)+\left(-35{\beta }_1+20{\beta }_2-10{\beta }_3+4{\beta }_4-{\beta }_5\right)=0.$ $c_n\left(-35{\alpha }_1-50{\alpha }_2+35{\alpha }_3-19{\alpha }_4+4{\alpha }_5\right)+\left(35{\beta }_1-15{\beta }_2+5{\beta }_3-{\beta }_4\right)=0.$ $c_n\left(21{\alpha }_1+27{\alpha }_2-13{\alpha }_3+3{\alpha }_4\right)+\left(-21{\beta }_1+6{\beta }_2-{\beta }_3\right)=0.$

$c_n\left(-7{\alpha }_1-8{\alpha }_2+2{\alpha }_3\right)+\left(7{\beta }_1-{\beta }_2\right)=0.$

$c_n\left({\alpha }_1+{\alpha }_2\right)-{\beta }_1=0.$

为表示简便，在这里系数 ${\alpha }_i\left(i=1,2,\cdots ,8\right)$ 表示

${\alpha }_1=-7p_{00}^{i+1}-7p_{00}^i,{\alpha }_2=7p_{10}^{i+1}+7p_{01}^i,{\alpha }_3=21p_{20}^{i+1}+21p_{02}^i,{\alpha }_4=35p_{30}^{i+1}+35p_{03}^i,$

${\alpha }_5=35p_{40}^{i+1}+35p_{04}^i,{\alpha }_6=7p_{50}^{i+1}+7p_{05}^i,{\alpha }_7=7p_{60}^{i+1}+7p_{06}^i,{\alpha }_8=7p_{70}^{i+1}+7p_{07}^i.$

对于系数 ${\beta }_i\left(i=1,2,\cdots ,8\right)$ 表示

${\beta }_1=7\left(p_{00}^{i+1}-p_{01}^{i+1}\right)+7\left(p_{10}^i-p_{00}^i\right),{\beta }_2=49\left(p_{10}^{i+1}-p_{11}^{i+1}\right)+49\left(p_{11}^i-p_{01}^i\right),$

${\beta }_3=147\left(p_{20}^{i+1}-p_{21}^{i+1}\right)+147\left(p_{12}^i-p_{02}^i\right),{\beta }_4=245\left(p_{30}^{i+1}-p_{31}^{i+1}\right)+245\left(p_{13}^i-p_{03}^i\right),$

${\beta }_5=245\left(p_{40}^{i+1}-p_{41}^{i+1}\right)+245\left(p_{14}^i-p_{04}^i\right),{\beta }_6=147\left(p_{50}^{i+1}-p_{51}^{i+1}\right)+147\left(p_{15}^i-p_{05}^i\right),$

${\beta }_7=7\left(p_{60}^{i+1}-p_{61}^{i+1}\right)+7\left(p_{16}^i-p_{06}^i\right),{\beta }_8=7\left(p_{70}^{i+1}-p_{71}^{i+1}\right)+7\left(p_{17}^i-p_{07}^i\right).$

以此类推，代入(4)也可以得到8个方程。因此，可以得到24个约束方程，其中含有8个独立方程自动满足循环条件，因此，实际上得到16个约束方程。

3.2.2. 外部边界的连续性

外部边界指的是在奇异点周围，不以奇异点为顶点，并且在拼接过程中与周围曲面片共有的边。在这里可以将外部边界的两个端点分为两类：第一类是与奇异点相邻的点；第二类是与奇异点相对的点。定义映射 ${\psi }_n:\left(u,v\right)\to \left(x,y\right)$ ，其中 $n$ 为对应奇异点的价。对于第一类点，令

$q\left(u,v\right)=\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]$ ， $s\left(u,v\right)=\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]$

对应的链式规则矩阵分别为 $Q$ ， $S$ 为 $n\times n$ 的块对角矩阵，其中，

$Q={\left[\begin{array}{ccccc}1& 0& & & \\ 0& -1& & & \\ & & \ddots & & \\ & & & 1& 0\\ & & & 0& -1\end{array}\right]}_{n\times n}$ ， $S={\left[\begin{array}{ccccc}0& 1& & & \\ 1& 0& & & \\ & & \ddots & & \\ & & & 0& 1\\ & & & 1& 0\end{array}\right]}_{n\times n}$

对于第一类点有

$I=\left(Q\cdot {\psi }_n\cdot S\cdot {\left(R_n^{-1}\cdot {\phi }_n\right)}^{-1}\cdot {\phi }_n\cdot {\psi }_n^{-1}\right)\cdot \left(Q\cdot {\psi }_m\cdot S\cdot {\left(R_m^{-1}\cdot {\phi }_m\right)}^{-1}\cdot {\phi }_m\cdot {\psi }_m^{-1}\right)$

令 ${\psi }_n\left(u,v\right)={\psi }_n\left(v,u\right)$ ，则 $S\cdot {\psi }_n\cdot S={\psi }_n$ ，对于第二类点有

$I={\psi }_k\cdot {\left(S\cdot {\psi }_k\cdot S\right)}^{-1}\cdot {\psi }_l\cdot {\left(S\cdot {\psi }_l\cdot S\right)}^{-1}\cdot {\psi }_m\cdot {\left(S\cdot {\psi }_k\cdot S\right)}^{-1}\cdot {\psi }_n\cdot {\left(S\cdot {\psi }_n\cdot S\right)}^{-1}$

由上述两个式子可以推得

$y\left(u,v\right)=b^4{\left(v\right)}^{\text{T}}\left[\begin{array}{cccc}0& 0& 0& 0\\ \frac{1}{4}& \frac{5c_n^3-3c_n^2-15c_n+18}{12\left(c_n-2\right)\left(2c_n-3\right)}& \frac{c_n^2+2c_n-6}{12\left(c_n-2\right)}& \frac{1}{4}\\ \frac{1}{2}& \frac{c_n+3}{6}& \frac{c_n^2+3c_n-9}{9\left(c_n-2\right)}& \frac{1}{2}\\ \frac{3}{4}& \frac{c_n+9}{12}& \frac{c_n+9}{12}& \frac{3}{4}\\ 1& 1& 1& 1\end{array}\right]b^3(u)$

其中， $c_n=\cos \left(\frac{\text{2π}}{n}\right)$ ，又由对称性可得 $x\left(u,v\right)=y\left(v,u\right)$ 。

在奇异点周围可以得到 $n$ 个“混合细分”的双三次B-样条曲面片，对于外部边界，通过运用双三次节点插入算子，将B样条转化为 $n$ 个Bernstein形式的双三次Bézier曲面片的控制点

$b^k=\frac{1}{36}{Q}^{\text{T}}AQ$ ，

其中 $Q=\left[\begin{array}{cccc}1& 0& 0& 0\\ 4& 4& 2& 1\\ 1& 2& 4& 4\\ 0& 0& 0& 1\end{array}\right]$ ，

$A$ 为“混合细分”矩阵。

对于奇异点处的第 $k$ 个曲面片，有

$A_k=\left[\begin{array}{cccc}1& a_{10}^{k+2}& a_{11}^{k+1}& a_{12}^{k+1}\\ a_{10}^{k-1}& a_{00}^k& a_{10}^{k+1}& a_{20}^{k+1}\\ a_{11}^{k-1}& a_{10}^k& a_{11}^k& a_{12}^k\\ a_{12}^{k-1}& a_{20}^k& a_{20}^k& a_{22}^k\end{array}\right]$

则Bézier曲面片的表达式

$B_i\left(u,v\right)=\frac{1}{36}{b}^3\left(u\right)Q^{\text{T}}\left[\begin{array}{cccc}1& a_{10}^{k+2}& a_{11}^{k+1}& a_{12}^{k+1}\\ a_{10}^{k-1}& a_{00}^k& a_{10}^{k+1}& a_{20}^{k+1}\\ a_{11}^{k-1}& a_{10}^k& a_{11}^k& a_{12}^k\\ a_{12}^{k-1}& a_{20}^k& a_{21}^k& a_{22}^k\end{array}\right]Q{\left(b^3\left(v\right)\right)}^{\text{T}}$

要使外部边界 $P_i\left(1,t\right)$ 满足二阶几何连续只需要满足下面条件即可

$\frac{{\partial }^j}{\partial u^j}{P}_i\left(1,t\right)=\frac{{\partial }^j}{\partial u^j}\left(B_i\circ {\psi }_n\right)\left(1,t\right),j=0,1,2.$

类似于内部边界，可以推得

$P_i\left(1,t\right)$ 、 $B_i\left(1,t\right)$ 、 $\frac{\partial }{\partial u}{P}_i\left(1,t\right)$ 、 $\frac{\partial }{\partial x}{B}_i\left(1,t\right)$ 、 $\frac{\partial }{\partial y}{B}_i\left(1,t\right)$ 、 $\frac{{\partial }^2}{\partial x^2}{B}_i\left(1,t\right)$ 、 $\frac{{\partial }^2}{\partial xy}{B}_i\left(1,t\right)$ 和 $\frac{{\partial }^2}{\partial y^2}{B}_i\left(1,t\right)$ 。同时，由已知条件可以得到

$\frac{\partial }{\partial u}{\psi }_{n,x}\left(1,t\right)=-\frac{1}{2}{\left(1-t\right)}^3-\frac{3\left(c_n+1\right)}{2}t{\left(1-t\right)}^2-\frac{3\left(c_n+1\right)}{2}{t}^2\left(1-t\right)-\frac{1}{2}{t}^3$

$\frac{\partial }{\partial u}{\psi }_{n,y}=\frac{-c_n^2+c_n}{c_n-2}t{\left(1-t\right)}^3+\frac{-2c_n^2+\frac{33}{2}{c}_n-27}{c_n-2}{t}^2{\left(1-t\right)}^2-\frac{3\left(c_n+9\right)}{2}{t}^3\left(1-t\right)$

$\frac{{\partial }^2}{\partial u^2}{\psi }_{n,x}\left(1,t\right)=-9{\left(1-t\right)}^3-3\left(c_n+9\right)t{\left(1-t\right)}^2+\frac{3\left(c_n^2+3c_n-6\right)}{c_n-2}{t}^2\left(1-t\right)-9t^3$

$\frac{{\partial }^2}{\partial u^2}{\psi }_{n,y}\left(1,t\right)=\frac{2c_n^3+2c_n^2}{\left(c_n-2\right)\left(2c_n-3\right)}t{\left(1-t\right)}^3-\frac{2c_n^2}{c_n-2}{t}^2{\left(1-t\right)}^2-3c_n{t}^3\left(1-t\right)$

对于外部边界 $\left(t,1\right)$ 连续性条件，每个方程的化简计算后次数不超过7次，含有8个系数，由对应系数相等，可以得到24个方程，同理，对于外部边界 $\left(1,t\right)$ 也可以得到24个方程。对于n价的奇异点则产生48n个约束方程。对于曲面片 $P_i\left(0,t\right)$ 和 $P_{i+1}\left(t,0\right)$ 的内部边界，由连续性条件可以得到24个约束方程，但是当 $t=1$ 时， $P_i\left(0,1\right)$ 和 $P_{i+1}\left(1,0\right)$ 所确定的9个二阶混合偏导自动满足约束条件，与外部边界的二阶连续性条件重复，因此，对于 $n$ 价的奇异点则产生16n个约束方程。这样，一共得到了64n个方程，对于价为 $n$ 的奇异点，双七次曲面环恰好有64n个未知数，因此可以求得该方程解即所求的控制顶点。

3.3. Bézier控制点的求解

对于简单的公式，可以直接以文本方式输入；对于复杂的公式，可以考虑使用公式编辑器，或者将公式制作成图片后插入文中。编辑公式的过程中要特别注意减号与连字符的区别，前者较长，后者较短。令双七次曲面片控制顶点的系数矩阵为 $C$ ，则对应的第 $i$ 个曲面片的系数矩阵为 $p_i$ ，为了保证曲面片整体的二阶连续性，则需要满足如下条件

$CP=WA$

其中， $C$ 为循环矩阵， $P$ 为双七次曲面片控制顶点，矩阵 $A$ 对应图2二环中第 $i$ 个曲面片的样条顶点。对于循环矩阵 [7] $C=bcirc\left(c_0,c_1,\cdots ,c_{n-1}\right)$ ，它的全部特征根为 ${\lambda }_i={\displaystyle \underset{j=0}{\overset{n}{\sum }}{c}_j{\omega }_i^j}$ ，其中 ${\omega }_j={\text{e}}^{i\frac{2\text{π}}{n+1}j}$ ，并且 $C$ 相似于对角阵 $diag\left({\lambda }_0,\cdots ,{\lambda }_{n-1}\right)$ 。

为了使二次能量函数 [8] 最小，接下来采用Lagrange乘子法 [9] ，可将约束条件转化为下列形式

对于第 $i$ 个块方程，有

$\{\begin{array}{l}{\displaystyle \underset{i=0}{\overset{n-1}{\sum }}{c}_i^{\text{T}}}{\lambda }_i+{\displaystyle \underset{i=0}{\overset{n-1}{\sum }}E{p}_i}=0\hfill \\ {\displaystyle \underset{i=0}{\overset{n-1}{\sum }}{c}_i{p}_i}={\displaystyle \underset{i=0}{\overset{n-1}{\sum }}{\omega }_i{a}_i}\hfill \end{array}$

其中，对于第 $i$ 个方程，只有 $p_i$ 的系数为 $E$ ，其余的均为0。

然后采用傅里叶变换 [10] (DFT)，令 ${\hat{d}}^l=\frac{1}{n}{\displaystyle \underset{k=0}{\overset{n-1}{\sum }}{\bar{E}}^{kl}{d}^k}$ ，其中， $k,l=0,1,\cdots ,n-1$ ， ${\bar{E}}_n^{k_l}=C_n^{k_l}-\sqrt{-1}{S}_n^{k_l}$ ，

Euler系数 $C_n^{k_l}=\cos \frac{2\text{π}{k}_l}{n}$ ， $S_n^{k_l}=\sin \frac{2\text{π}{k}_l}{n}$ 则有

${\hat{\rho }}_j=\frac{1}{n}{\displaystyle \underset{k=0}{\overset{n-1}{\sum }}{\tau }_k}{\rho }_k$ (5)

其中， ${\tau }_k=cos\frac{2\text{π}k}{n}-\sqrt{-1}sin\frac{2\text{π}k}{n}$ 。

将 $p_i,c_i,{\lambda }_i,{\omega }_i,a_i$ 代入(5)式即可得到 ${\hat{p}}_j,{\hat{c}}_j,{\hat{\lambda }}_j,{\hat{\omega }}_j,{\hat{a}}_j$ 。

于是上式可转化为 $\{\begin{array}{l}E{\hat{p}}_j+{\hat{c}}_j^H{\hat{\lambda }}_j=0\hfill \\ {\hat{c}}_j{\hat{p}}_j={\hat{w}}_j{\hat{a}}_j\hfill \end{array}$ 解得 ${\hat{p}}_j=\frac{{\hat{w}}_j{\hat{a}}_j}{{\hat{c}}_j}$ (其中 $j=0,1,\cdots ,n-1$ )，求得 ${\hat{p}}_j$ 后再利用傅里叶逆变换(IDFT)，将 ${\hat{p}}_j,{\hat{c}}_j,{\hat{\lambda }}_j,{\hat{\omega }}_j,{\hat{a}}_j$ 代入(5)式即可反求出相应的双七次控制顶点 $p_i$ 。

4. 曲面造型结果

下面是运用本文的算法对某品牌汽车进行样条重建的实例，对工业设计中汽车车身出现的5-边洞和3-边洞进行局部填充，图3(a)和图4(a)是汽车局部网格图；图3(b)和图4(b)是网格的控制顶点图，为了清楚进行了一些旋转；图3(c)和图4(c)是奇异点附近的网格填充图。整个填充过程如图3和图4所示。

下面是采用本节的方法对汽车局部模型进行曲面造型的结果。该汽车局部散乱点为1452个，“混合细分”后的四边形网格有323个，图中分给出了10个奇异点多变行动填充。图5(a)为“混合细分”得到的四边形控制网格，图5(b)为奇异点的二环结构，图5(c)经过双七次样条拟合后，奇异点周围各曲面片之间满足 $G^2$ 连续，具体过程如图5所示。

(a)3价奇异点局部网格 (b)3奇异点的控制顶点 (c)3价奇异点的拟合

Figure 3. 3-sided hole filling

图3. 3-边洞的填充

(a)5价奇异点局部网格 (b)5奇异点的控制顶点 (c)5价奇异点的拟合

Figure 4. 5-sided hole filling

图4. 5-边洞的填充

(a) 进行“混合细分” (b) 对n边洞进行填充 (c) G²连续曲面

Figure 5. Surface modeling with G² continuity

图5. G²连续曲面造型

5. 小结

本文致力于细分曲面造型中奇异点处的 $G^2$ 连续性研究，给出了奇异点处的 $G^2$ 算法，主要给出了细分曲面中奇异点处的 $G^2$ 重建算法，主要分为三个部分：第一，给出了形状可调的C-C细分算法；第二，具体给出了奇异点处的 $G^2$ 重建方法；第三，给出了奇异点处3边洞和5边洞的填充结果。增加了工业设计中细分曲面的可调性，逐步拓宽了应用范围。

基金项目

国家自然科学基金(51175248/E050603)；南京航空航天大学基本科研业务费资助(NZ2013201)。