计算机辅助几何设计(7)：细分曲线与曲面

1. 简介

样条曲面的局限性

连续张量积样条曲面的参数域仅在特定的方形网格上定义，因此拓扑受限
将多个参数域装配到一个曲面非常繁琐，很难获得连续性保证
曲线的处理与剪枝很困难

细分方法的优点

提供几何体的粗略表示
获得光滑的表示
通过递归进行实现

细分方法的缺点

没有了参数化的表示

细分方法的基本步骤

细分当前多边形/网格
插入插值点（分裂）
移动点：局部加权平均（均一化）
- 对于所有的点：逼近方法
- 仅对于新生成的点：插值方法

2. 细分曲线

2.1. Corner Cutting Splines

2.1.1. 基本流程：

将每条线段分为两半
顺时针地与其邻接点求平均
重复上述操作

能够收敛至二次B样条曲线

2.1.2. 矩阵表示

记：

$l$阶控制点为：$\pmb p_i^{(l)}$
$l+1$阶分裂点为：$\tilde{\pmb p}_i^{(l+1)}$
$l+1$阶平均控制点：$\pmb p_i^{(l+1)}$

分裂操作的矩阵表示：
$$
\begin{pmatrix}
\vdots\\\tilde x_{2i}^{(l+1)}\\\tilde x_{2i+1}^{(l+1)}\\\vdots
\end{pmatrix}{2n\times 1}=
\begin{pmatrix}
\ddots\\
&1\\
&\frac{1}{2}&\frac{1}{2}\\
&&1\\
&&\frac{1}{2}&\frac{1}{2}\\
&&&&\ddots
\end{pmatrix}{2n\times n}
\begin{pmatrix}
\vdots\\ x_{i}^{(l)}\\ x_{i+1}^{(l)}\\\vdots
\end{pmatrix}{n\times 1}
$$
平均操作的矩阵表示：
$$
\begin{pmatrix}
\vdots\\ x{2i}^{(l+1)}\\ x_{2i+1}^{(l+1)}\\\vdots
\end{pmatrix}{2n\times 1}=
\begin{pmatrix}
\ddots\\
&\frac{1}{2}&\frac{1}{2}\\
&&\frac{1}{2}&\frac{1}{2}\\
&&&&\ddots
\end{pmatrix}{2n\times 2n}
\begin{pmatrix}
\vdots\\\tilde x_{2i}^{(l+1)}\\\tilde x_{2i+1}^{(l+1)}\\\vdots
\end{pmatrix}_{2n\times 1}
$$

2.1.3. Chaikin’s Corner Cutting表示方法

对每条边，插入$\frac{1}{4}$和$\frac{3}{4}$的点，删除旧点，进行重复操作
极限曲线将是$C^1$连续的

对于第$k+1$阶细分，有：
$$
\begin{aligned}
P^{(k+1)}_ {2i}&=\frac{3}{4}P^{(k)}_ i+\frac{1}{4}P^{(k)}_ {i+1}\\\\
P^{(k+1)}_ {2i+1}&=\frac{1}{4}P^{(k)}_ i+\frac{3}{4}P^{(k)}_ {i+1}
\end{aligned}
$$
每次迭代，点数加倍，写成矩阵形式：
$$
\begin{pmatrix}
\vdots\\
p^{k+1}{2i-2}\\
p^{k+1}{2i-2}\\
p^{k+1}{2i-2}\\
p^{k+1}{2i-2}\\
p^{k+1}{2i-2}\\
p^{k+1}{2i-2}\\
\vdots
\end{pmatrix}{2n\times 1}
=\frac{1}{4}
\begin{pmatrix}
\ddots\\
& 3 & 1 & 0 & 0\\
& 1 & 3 & 0 & 0\\
& 0 & 3 & 1 & 0\\
& 0 & 1 & 3 & 0\\
& 0 & 0 & 3 & 1\\
& 0 & 0 & 1 & 3\\
& & & & &\ddots\
\end{pmatrix}{2n\times n}
\begin{pmatrix}
\vdots\\
p^{k+1}{i-1}\\
p^{k+1}{i}\\
p^{k+1}{i+1}\\
p^{k+1}{i+2}\\
\vdots\
\end{pmatrix}_{n\times 1}
$$

2.2. Lane - Riesenfeld Subdivision

2.2.1. 基本流程

各条边加中点
将每条边用中点代替，重复迭代$d$次（$d+1$次B样条细分）

示例：

$d=2$

$$ \begin{aligned} a_1&=\frac{A+B}{2}\\\\ c_1&=\frac{B+C}{2} \end{aligned} $$ $$ \begin{aligned} a_2&=\frac{B+a_1}{2}\\\\ c_2&=\frac{B+c_1}{2} \end{aligned} $$ $$ \begin{aligned} b_1&=\frac{a_2+c_2}{2}\\\\ &=\frac{a_1+2B+c_1}{4}\\\\ &=\frac{A+6B+C}{8} \end{aligned} $$ 即： $$ \begin{pmatrix} a_1\\\\b_1\\\\c_1 \end{pmatrix}= \frac{1}{8}\begin{pmatrix} 4&4&0\\\\ 1&6&1\\\\ 0&4&4 \end{pmatrix} \begin{pmatrix} A\\\\B\\\\C \end{pmatrix} $$

2.2.2. 三次B样条细分的矩阵表示

当$d=2$时，

$$
\begin{pmatrix}\vdots\\
\pmb p_{2i}^{(l+1)}\\
\pmb p_{2i+1}^{(l+1)}\\
\vdots\end{pmatrix}{2n\times1}=
\begin{pmatrix}
\ddots\\
&\frac{1}{8}&\frac{3}{4}&\frac{1}{8}\\
&&\frac{1}{2}&\frac{1}{2}\\
&&\frac{1}{8}&\frac{3}{4}&\frac{1}{8}\\
&&&\frac{1}{2}&\frac{1}{2}\\
&&&&&\ddots
\end{pmatrix}{2n\times n}
\begin{pmatrix}
\vdots\\
\pmb p_{i}^{(l)}\\
\pmb p_{i+1}^{(l)}\\
\vdots
\end{pmatrix}_{n\times 1}
$$

即：

$$
\begin{aligned}
\pmb p_{2i}^{(l+1)}&=\frac{1}{2}\pmb p_i^{(l)}+\frac{1}{2}\pmb p_{i+1}^{(l)}\\\\
\pmb p_{2i+1}^{(l+1)}&=\frac{1}{8}\pmb p_i^{(l)}+\frac{6}{8}\pmb p_{i+1}^{(l)}+\frac{1}{8}\pmb p_{i+2}^{(l)}
\end{aligned}
$$

将平均与分裂操作分开表示：

$$
\begin{pmatrix}
\vdots\\
\pmb p_{2i}^{(l+1)}\\
\pmb p_{2i+1}^{(l+1)}\\
\vdots
\end{pmatrix}{2n\times1}=
\underbrace{
\begin{pmatrix}
\ddots\\
&\frac{1}{4}&\frac{1}{2}&\frac{1}{4}\\
&&\frac{1}{4}&\frac{1}{2}&\frac{1}{4}\\
&&&\frac{1}{4}&\frac{1}{2}&\frac{1}{4}\\
&&&&\frac{1}{4}&\frac{1}{2}&\frac{1}{4}\\
&&&&&&\ddots
\end{pmatrix}{2n\times 2n}}{\mathrm{averaging}}
\underbrace{
\begin{pmatrix}
\ddots\\
&1\\
&&\frac{1}{2}&\frac{1}{2}\\
&&&1\\
&&&\frac{1}{2}&\frac{1}{2}\\
&&&&&\ddots
\end{pmatrix}{2n\times n}}{\mathrm{spliting}}
\begin{pmatrix}
\vdots\\
\pmb p{i}^{(l)}\\
\pmb p_{i+1}^{(l)}\\
\vdots
\end{pmatrix}_{n\times 1}
$$
极限曲线具有$C^2$连续性

核表示为：
$$
h=\frac{1}{8}[\cdots,0,0,1,4,6,4,1,0,0,\cdots]
$$
对于一般情况，核为：
$$
\frac{1}{2^{d+1}}\left(\dots,C_{d+2}^0,\dots,C_{d+2}^{d+2},\dots\right)
$$

2.2.3. 谱收敛分析

谱极限技巧

为了解决以下问题：

为了得到良好的逼近需要进行多次细分
过多的细分可能会导致控制点过剩（类比LOD）
能否直接计算控制点的极限位置？

注意到：

每一个曲线上的点只受一定数量的控制点影响
每个点$\pmb p^{(l+1)}$影响的范围均小于$\pmb p^{(l)}$的影响范围
对于每个邻域，都使用相同的细分矩阵

以三次B样条细分为例，
$$
\left( \begin{array} { c }
{ x _ { - } ^ { ( l + 1 ) } } \\
{ x ^ { ( l + 1 ) } } \\
{ x _ { + } ^ { ( l + 1 ) } } \end{array}
\right)
= \left( \begin{array} { c c c } { \frac { 1 } { 2 } } & { \frac { 1 } { 2 } } & { 0 } \\
{ \frac { 1 } { 8 } } & { \frac { 3 } { 4 } } & { \frac { 1 } { 8 } } \\
{ 0 } & { \frac { 1 } { 2 } } & { \frac { 1 } { 2 } }
\end{array} \right)
\left( \begin{array} { c }
{ x _ { - } ^ { ( l ) } } \\
{ x ^ { ( l ) } } \\
{ x _ { + } ^ { ( l ) } }
\end{array} \right)
$$
其中，

$x_-$为左邻域点
$x$为当前点
$x_+$为右邻域点

对转移矩阵进行相似对角化：
$$
M
=\left( \begin{array} { c c c }
{ \frac { 1 } { 2 } } & { \frac { 1 } { 2 } } & { 0 } \\
{ \frac { 1 } { 8 } } & { \frac { 3 } { 4 } } & { \frac { 1 } { 8 } } \\
{ 0 } & { \frac { 1 } { 2 } } & { \frac { 1 } { 2 } }
\end{array} \right)
=\left( \begin{array} { c c c }
{ 1 } & { -1 } & { -2 } \\
{ 1 } & { 0 } & { 1 } \\
{ 1 } & { 1 } & { -2 }
\end{array} \right)
\left( \begin{array} { c c c }
{ 1 } & { 0 } & { 0 } \\
{ 0 } & { \frac{1}{2} } & { 0 } \\
{ 0 } & { 0 } & { \frac{1}{4} }
\end{array} \right)
\left( \begin{array} { c c c }
{ \frac{1}{6} } & { \frac{2}{3} } & { \frac{1}{6} } \\
{ -\frac{1}{2} } & { 0 } & { \frac{1}{2} } \\
{ -\frac{1}{6} } & { \frac{1}{3} } & { -\frac{1}{6} }
\end{array} \right)\triangleq UDU^{-1}
$$
则求极限：
$$
\begin{aligned}
\begin{pmatrix}
x_-^{(\infty)}\\x^{(\infty)}\\x_+^{(\infty)}
\end{pmatrix}&=
\lim_{l\rightarrow\infty}M^l\begin{pmatrix}
x_-^{(0)}\\x^{(0)}\\x_+^{(0)}
\end{pmatrix}\\
&=\lim_{l\rightarrow\infty}UD^lU^{-1}\begin{pmatrix}
x_-^{(0)}\\x^{(0)}\\x_+^{(0)}
\end{pmatrix}\\
&=U\lim_{l\rightarrow\infty}D^lU^{-1}\begin{pmatrix}
x_-^{(0)}\\x^{(0)}\\x_+^{(0)}
\end{pmatrix}\\
&=U\lim_{l\rightarrow\infty}\begin{pmatrix}
1&0&0\\
0&0&0\\
0&0&0
\end{pmatrix}U^{-1}\begin{pmatrix}
x_-^{(0)}\\x^{(0)}\\x_+^{(0)}
\end{pmatrix}
\end{aligned}
$$
解得：
$$
x^{(\infty)}=\begin{pmatrix}\dfrac{1}{6}&\dfrac{2}{3}&\dfrac{1}{6}\end{pmatrix}\begin{pmatrix}
x_-^{(0)}\\x^{(0)}\\x_+^{(0)}
\end{pmatrix}
$$

收敛的必要条件：最大特征值的绝对值应为1
仿射不变性要求局部矩阵行和为 1，这意味着 $\pmb{1}$ 向量必须是特征值 $1$ 的特征向量

3. 细分曲面

3.1. 张量积B样条细分曲面

从一个四边形网格出发
对每一步细分操作
- 将每个四边形划分成四个（四叉树细分）
- 对顶点进行线性插值
- 使用二维的平均掩膜

3.1.1. 双二次细分曲面

B样条块的矩阵表示：
$$
P(u,v)=
\begin{pmatrix}
1&u&u^2
\end{pmatrix}
MPM^T
\begin{pmatrix}
1\\v\\v^2
\end{pmatrix}
$$
其中，
$$
\begin{aligned}
M&=\frac{1}{2}\begin{pmatrix}
1&1&0\\
-2&2&0\\
1&-2&1
\end{pmatrix}\\\\
P&=\begin{pmatrix}
P_{0,0}&P_{0,1}&P_{0,2}\\
P_{1,0}&P_{1,1}&P_{1,2}\\
P_{2,0}&P_{2,1}&P_{2,2}
\end{pmatrix}
\end{aligned}
$$
考虑$2\times2$块中的一个象限，例如$u,v\in[0,\frac{1}{2}]$，考虑新的曲面块$P’$参数为$u’=\frac{u}{2}$和$v’=\frac{v}{2}$

则有：
$$
\begin{aligned}
P^\prime(u,v)
&=P(\frac{u}{2},\frac{v}{2})\\\\
&=\left[\begin{matrix}
1&\frac{u}{2}&\frac{u^2}{4}
\end{matrix}\right]MPM^\top
\left[\begin{matrix}
1\\\frac{v}{2}\\\frac{v^2}{4}
\end{matrix}\right]\\\\
&\dots\\\\
&=\left[\begin{matrix}
1 & u & u^2
\end{matrix}\right]MP^\prime M^\top \left[\begin{matrix}
1\\v\\v^2
\end{matrix}\right]
\end{aligned}
$$
其中，
$$
P^\prime=SPS^T
$$

$$
S=M^{-1}\left[\begin{matrix}
1 & 0 & 0\\
0 & \frac{1}{2} & 0\\
0 & 0 & \frac{1}{4}
\end{matrix}\right]M=\frac{1}{4}\left[\begin{matrix}
3 & 4 & 0\\
1 & 3 & 0\\
0 & 3 & 1
\end{matrix}\right]
$$

3.1.2. 双三次细分曲面

B样条块的矩阵表示：
$$
P(u,v)=
\begin{pmatrix}
u^3&u^2&u&1
\end{pmatrix}
MPM^T
\begin{pmatrix}
u^3\\u^2\\u\\1
\end{pmatrix}
$$
其中，
$$
M=\frac{1}{6}\begin{pmatrix}
-1&3&-3&1\\
3&-6&3&0\\
-3&0&3&0\\
1&4&1&0
\end{pmatrix}
$$
同样考虑$3\times 3$曲面块的一个象限，例如$u,v\in[0,\frac{1}{2}]$，考虑新的曲面块$P’$参数为$u’=\frac{u}{2}$和$v’=\frac{v}{2}$

能够得到相似的结果：
$$
\begin{aligned}
P’&=SPS^T\\\\
S&=\frac{1}{8}\begin{pmatrix}
4 & 4 & 0 & 0\\
1 & 6 & 1 & 0\\
0 & 4 & 4 & 0\\
0 & 1 & 6 & 1
\end{pmatrix}
\end{aligned}
$$

3.1.3. 细分与平均掩膜Mask

细分掩膜

平均掩膜

3.2. Catmull-Clark Subdivision

3.2.1. 动机

B样条细分曲面存在的问题：

仅适用于内部或常规四边形网格
相比标准的B样条构造并没有灵活性优势
需要一个更好的方法：
- 处理任意拓扑的四边形网格
- 处理边界区域的情况

解决方法：Catmull-Clark细分方法

3.2.2. 算法流程

记：

对于不是四边形的面，称为Non-quad face
所有度不为4的顶点称为奇异点

算法流程：

在每个面中添加一个点，每条边的中点也添加一个点，面上的新顶点连接所有边上的新顶点，如图所示：

可以看到：
- 有几个非四边形面就会多出几个奇异点
- 新多出来的奇异点的度数与原来所在面的边数相等
- 第一次细分后所有面都会变成四边形面，且往后奇异点数目不再增加
进行顶点位置的偏移
- 对于面顶点
  
  $$
  f=\frac{v_1+v_2+v_3+v_4}{4}
  $$
- 对于边顶点
  
  $$
  e=\frac{v_1+v_2+f_1+f_2}{4}
  $$
- 对于原来的顶点
  
  $$
  v=\frac{f_1+f_2+f_3+f_4+2(m_1+m_2+m_3+m_4)+4p}{16}
  $$
  其中，$m$为边中点，$p$为旧顶点