提升小波

前言

经典小波分析是从傅里叶分析的基础上发展起来的，因此它在一定程度上收到傅里叶分析的限制。小波变换和多分辨率分析都是建立在小波的二进伸缩平移的基础上，称之为第一代小波。1996年，Sweldens提出了不依赖傅里叶变换的小波提升方法，提升小波主要具有以下性质。

1. 能够用于构造第一代小波：可以用来提高小波的消失矩，构造具有插值的小波等，可以加深对第一代小波的理解。
2. 能够改进第一代小波的算法性能：完全是原位计算且容易推广到边界。
3. 可用于构造第二代小波：其摆脱了傅里叶变换的限制，同时满足小波变换的好性质，该思想可以用来构造复杂区域上的小波变换。

Z变换下的多分辨率分析

利用Z变换重新认识多分辨率分析的分解重构算法。

之前定义过序列的Z变换 \[ x(z)=\sum_{j\in\mathbb{Z}} x_jz^{-j} \] 定义偶序列对应的Z变换 \[ x_e(z)=\sum_{j\in\mathbb{Z}} x_{2j}z^{-j} \] 则 \[ x_e(z^2)=\sum_{j\in\mathbb{Z}} x_{2j}z^{-2j} \] 同样定义奇序列的Z变换满足 \[ x_o(z^2)=\sum_{j\in\mathbb{Z}} x_{2j+1}z^{-2j} \] 则有 \[ \left( \begin{matrix} x_e(z^2)\\ x_o(z^2) \end{matrix} \right) = \left( \begin{matrix} \frac12&\frac12\\ \frac z2&-\frac z2 \end{matrix} \right) \left( \begin{matrix} x(z)\\ x(-z) \end{matrix} \right) \] 回忆之前的分解重构算法： \[ \begin{aligned} c_{j-1,l}&=2^{-1/2}\sum c_{j,k}\overline{h_{k-2l}}\\ d_{j-1,l}&=2^{-1/2}\sum c_{j,k}\overline{g_{k-2l}}\\ c_{j,l}=2^{-1/2}\sum &c_{j-1,k}h_{l-2k}+2^{-1/2}\sum d_{j-1,k}g_{l-2k} \end{aligned} \] 考虑双尺度系数和小波系数Z变换得到的函数 \[ h(z)=\sum_{j\in\mathbb{Z}} h_jz^{-j},g(z)=\sum_{j\in\mathbb{Z}} g_jz^{-j} \] 由性质可得 \[ g(z)=-z^{-1}\bar{h}(-z^{-1}),\bar g(z)=-z^{-1}h(-z^{-1}) \] 且奇序列偶序列Z变换 \[ h_e(z^2)=\frac{h(z)+h(-z)}{2},h_o(z^2)=\frac{h(z)-h(-z)}{2z^{-1}} \]

\[ g_e(z^2)=\frac{g(z)+g(-z)}{2},g_o(z^2)=\frac{g(z)-g(-z)}{2z^{-1}} \]

分解算法： \[ \begin{aligned} c^{j-1}(z^2)&=\sum_l c_{j-1,l}z^{-2l}\\ &=2^{-1/2}\sum_l(\sum_k c_{j,k}\overline{h_{k-2l}})z^{-2l}\\ &=2^{-1/2}c^j(z)\bar{h}(z^{-1}) \end{aligned} \] 用-z代替z可得 \[ c^{j-1}(z^2)=2^{-1/2}c^j(-z)\bar{h}(-z^{-1}) \] 于是 \[ c^{j-1}(z^2)=2^{-1/2}\frac{c^j(z)\bar{h}(z^{-1})+c^j(-z)\bar{h}(-z^{-1})}{2} \] 同理 \[ d^{j-1}(z^2)=2^{-1/2}\frac{c^j(z)\bar{g}(z^{-1})+c^j(-z)\bar{g}(-z^{-1})}{2} \] 可推出 \[ \left( \begin{matrix} c^{j-1}(z^2)\\ d^{j-1}(z^2) \end{matrix} \right) = 2^{-1/2} \left( \begin{matrix} \bar{h}_e(z^{-2})&\bar{h}_o(z^{-2})\\ \bar{g}_e(z^{-2})&\bar{g}_o(z^{-2}) \end{matrix} \right) \left( \begin{matrix} c^j_e(z^2)\\ c^j_o(z^2) \end{matrix} \right) \] 定义h和g的多相位矩阵为 \[ P(z)= \left( \begin{matrix} h_e(z)&g_e(z)\\ h_o(z)&g_o(z)\\ \end{matrix} \right) \] 则小波分解算法可记为 \[ \left( \begin{matrix} c^{j-1}(z^2)\\ d^{j-1}(z^2) \end{matrix} \right) = 2^{-1/2} \bar{P}(z^{-2})^T \left( \begin{matrix} c^j_e(z^2)\\ c^j_o(z^2) \end{matrix} \right) \] 即 \[ \left( \begin{matrix} c^{j-1}(z)\\ d^{j-1}(z) \end{matrix} \right) = 2^{-1/2} \bar{P}(z^{-1})^T \left( \begin{matrix} c^j_e(z)\\ c^j_o(z) \end{matrix} \right) \] 类似地重构算法： \[ \left( \begin{matrix} c^j_e(z)\\ c^j_o(z) \end{matrix} \right) = 2^{-1/2} P(z) \left( \begin{matrix} c^{j-1}(z)\\ d^{j-1}(z) \end{matrix} \right) \] 完全重构的条件为 \[ P(z)\bar{P}(z^{-1})^T=2I \]

提升小波的构造

可以看出分解重构都只需要用到多相位矩阵\(P(z)\)，简单的想法是将\(P(z)\)分解为一系列简单矩阵的乘积，则只需要作原位计算。

事实上有结论，若\(P(z)\)的行列式等于-2（若多相位矩阵来自于多分辨率分析，则这一点必然成立），

则存在Laurent多项式\(u_i(z),p_i(z)\)和非零常数\(K\)满足 \[ P(z)=\prod_{i=1}^m \left( \begin{matrix} 1&u_i(z)\\ 0&1 \end{matrix} \right) \left( \begin{matrix} 1&0\\ p_i(z)&1 \end{matrix} \right) \left( \begin{matrix} K&0\\ 0&-\frac{2}{K} \end{matrix} \right),p_m(z)=0 \] 下面给出这些Laurent多项式的构造方法

多相位矩阵的因子分解

Laurent多项式

\[ h(z)=\sum_{k=k_a}^{k_b}h_kz^{-k},h_{k_a}\ne 0 ,h_{k_b}\ne 0 \]

称为一个Laurent多项式。次数为 \[ |h(z)|=k_b-k_a \]

Laurent多项式的带余除法可以和普通多项式类似定义，但不唯一，任取其一即可。

构造方法

使用Euclid算法 \[ \begin{aligned} a_{i+1}(z)&=b_i(z)\\ b_{i+1}(z)&=a_i(z)\%b_i(z) \end{aligned} \] 总可以使\(b_n(z)=0\)的\(n\)为偶数且尽量小，则 \[ \left( \begin{matrix} a_n(z)\\ 0 \end{matrix} \right) = \prod_{i=n}^1 \left( \begin{matrix} 0&1\\ 1&-q_i(z) \end{matrix} \right) \left( \begin{matrix} a(z)\\ b(z) \end{matrix} \right) \]

这等价于

\[ \left( \begin{matrix} a(z)\\ b(z) \end{matrix} \right) = \prod_{i=1}^n \left( \begin{matrix} q_i(z)&1\\ 1&0 \end{matrix} \right) \left( \begin{matrix} a_n(z)\\ 0 \end{matrix} \right) \]

从而

\[ \left( \begin{matrix} a(z)\\ b(z) \end{matrix} \right) = \left( \begin{matrix} 1&q_1(z)\\ 0&1 \end{matrix} \right) \left( \begin{matrix} 1&0\\ q_2(z)&1 \end{matrix} \right) \prod_{i=3}^n \left( \begin{matrix} q_i(z)&1\\ 1&0 \end{matrix} \right) \left( \begin{matrix} a_n(z)\\ 0 \end{matrix} \right) \]

对

\[ \left( \begin{matrix} h_e(z)\\ h_o(z)\\ \end{matrix} \right) \]

操作即有

\[ \left( \begin{matrix} h_e(z)\\ h_o(z)\\ \end{matrix} \right) =\prod_{i=1}^{m-1} \left( \begin{matrix} 1&u_i(z)\\ 0&1 \end{matrix} \right) \left( \begin{matrix} 1&0\\ p_i(z)&1 \end{matrix} \right) \left( \begin{matrix} K\\ 0 \end{matrix} \right) \]

记

\[ \left( \begin{matrix} h_e(z)&p_e(z)\\ h_o(z)&p_o(z)\\ \end{matrix} \right) =\prod_{i=1}^{m-1} \left( \begin{matrix} 1&u_i(z)\\ 0&1 \end{matrix} \right) \left( \begin{matrix} 1&0\\ p_i(z)&1 \end{matrix} \right) \left( \begin{matrix} K&0\\ 0&-\frac{2}{K} \end{matrix} \right) \]

结合我们希望得到的形式

\[ P(z) =\prod_{i=1}^{m-1} \left( \begin{matrix} 1&u_i(z)\\ 0&1 \end{matrix} \right) \left( \begin{matrix} 1&0\\ p_i(z)&1 \end{matrix} \right) \left( \begin{matrix} K&0\\ 0&-\frac{2}{K} \end{matrix} \right) \left( \begin{matrix} 1&\frac{2u_m(z)}{K^2}\\ 0&1 \end{matrix} \right) \] 只要取 \[ u_m(z)=-\frac{K^2}{4}(p_o(z)g_e(z)-p_e(z)g_o(z)) \] 即可。

提升算法

\[ P(z)=\prod_{i=1}^m \left( \begin{matrix} 1&u_i(z)\\ 0&1 \end{matrix} \right) \left( \begin{matrix} 1&0\\ p_i(z)&1 \end{matrix} \right) \left( \begin{matrix} K&0\\ 0&-\frac{2}{K} \end{matrix} \right) \]

从而

\[ P^{-1}(z)= \left( \begin{matrix} \frac{1}{K}&0\\ 0&-\frac{K}{2} \end{matrix} \right) \prod_{i=m}^1 \left( \begin{matrix} 1&0\\ -p_i(z)&1 \end{matrix} \right) \left( \begin{matrix} 1&-u_i(z)\\ 0&1 \end{matrix} \right) \]

因此代入到分解重构公式中：

\[ \left( \begin{matrix} c^{j-1}(z)\\ d^{j-1}(z) \end{matrix} \right) = 2^{-1/2} \bar{P}(z^{-1})^T \left( \begin{matrix} c^j_e(z)\\ c^j_o(z) \end{matrix} \right) \]

\[ \left( \begin{matrix} c^j_e(z)\\ c^j_o(z) \end{matrix} \right) = 2^{-1/2} P(z) \left( \begin{matrix} c^{j-1}(z)\\ d^{j-1}(z) \end{matrix} \right) \]

从而分解算法：

1. 先分成奇偶序列（懒小波变换）
2. 保持奇序列，偶序列更新为\(c_e^j(z)-u_i(z)c_o^j(z)\)
3. 保持偶序列，奇序列更新为\(c_o^j(z)-p_i(z)c_e^j(z)\)
4. 重复23，最终做一个放缩变换。

重构算法恰好将之前的步骤全都反过来即可。

以上所有的过程都可原位运算！

提升小波的经典应用——细分曲面的构造

在计算机图形学领域，单一的参数曲面只能表达亏格为1的曲面，一个自然的思路就是利用多片曲面的拼接来表达一个复杂图形。然而边界上的光滑拼接是一个非常难的课题，细分表示就是在这样的背景下产生的。细分曲面的产生思路就类似于多分辨率分析，通过简单的规则不断加入细节就可以得到一个极限曲面，这就是我们想要的光滑曲面。

基于多分辨率分析的细分曲面

对于网格\(M^0\)，假设它的顶点为\(V_i^0\)，则极限曲面\(S(x)\)可以看作定义在\(M^0\)上的参数方程，从而存在\(\phi_i^0(x)\)使得 \[ S(x)=\sum\phi_i^0(x)V_i^0 \] \(\phi_i^0\)在\(V_i^0\)取值1，其它顶点取值0。对于细分后的网格\(M^j\)类似有 \[ S(x)=\sum\phi_i^j(x)V_i^j \] 记 \[ \begin{aligned} \Phi^j(x)&=(\phi_0^j(x),\phi_1^j(x),\cdots)\\ V^j&=(V^j_0,V_1^j,\cdots)\\ T^j(M^0)&=span\{\phi_0^j(x),\phi_1^j(x),\cdots\} \end{aligned} \] 则\(\forall i, \forall j\le k\)有 \[ \phi_i^j(x)\in T^k(M^0) \] 上述规则定义了一个嵌套空间序列 \[ T^0(M^0)\subset T^1(M^0)\subset T^2(M^0)\subset \cdots \]

根据细分规则，有 \[ \Phi^j(x)=\Phi^{j+1}(x)P^j \] 记 \[ S^j(x)=\Phi^j(x)V^j \] 引入内积 \[ \langle f,g\rangle=\sum_{\tau\in F(M^0)}\frac{1}{Area(\tau)}\int_{\tau}f(x)g(x)dx,f,g\in T^j(M^0) \] 对应小波空间即为\(T^j(M^0)\)在\(T^{j+1}(M^0)\)中的补空间。

这样只能整体操作，还有求解线性方程组等操作，时间复杂度非常高。

基于Loop细分的提升小波构造

Loop细分的基网格为三角网格。每次细分，将上一层的每一个三角形分成四个新的三角形。

Loop细分具体规则

1. 对于每条边，找到其所对应的两个三角形，如图

   新的边点e定义为 \[ e=\frac38(V_0+V_1)+\frac18(V_2+V_3) \]

2. 对于原来的每个顶点\(V\)，假设邻居\(V_i,i=1,\cdots, n\)，则更新它的位置为 \[ v=\alpha_nV+\sum_{i=1}^n\beta_nV_i \] 其中 \[ \alpha_n=\frac38+(\frac38+\frac14\cos \frac{2\pi}{n})^2,\beta_n=\frac{1-\alpha_n}{n} \]

Loop细分的提升小波

使用上面的办法，新的顶点位置需要用原来的顶点位置组合更新，则不能原位操作。我们将顶点和边点分别看作奇序列偶序列，希望找到用奇序列更新偶序列，偶序列更新奇序列的办法。

由 \[ e=\frac38(V_0+V_1)+\frac18(V_2+V_3) \]

\[ v=\alpha_nV+\sum_{i=1}^n\beta_nV_i \]

求和得 \[ \sum e_i=\frac{3n}{8}V+\frac58\sum V_i \] 代入得 \[ v=\gamma_nV+\delta_n\sum e_i \] 其中 \[ \gamma_n=\frac85\alpha_n-\frac35,\delta_n=\frac85\beta_n \] 于是Loop细分提升小波重构算法： \[ \begin{cases} e\gets e+\frac38(v_0+v_1)+\frac18(v_2+v_3)\\ v\gets \gamma_n v\\ v_i\gets v_i+\delta_n e \end{cases} \] 分解算法： \[ \begin{cases} v_i\gets v_i-\delta_n e\\ v\gets \frac{1}{\gamma_n} v\\ e\gets e-\frac38(v_0+v_1)-\frac18(v_2+v_3)\\ \end{cases} \] 这样的算法就简单高效得多，但是正交性不佳。

双正交Loop细分小波提升构造方法

找 \[ \psi(x)=\psi_k(x)+\sum_{i=0}^3\omega_i\phi_i(x) \] 使得\(\langle \psi(x),\phi_i(x)\rangle=0\)

设\(\langle \phi_i(x),\phi_j(x)\rangle=a_{i,j},-\langle \psi_k(x),\phi_i(x)\rangle=b_{i}\)，要解\(A\omega=b\)

连续的内积不是很方便计算，所以在实际应用中，内积一般使用的是离散内积。

首先将\(\phi_i,\psi_k\)表示成下一层基函数的线性组合，并假设下一层基函数两两内积相同，从而连续内积即是下一层基函数系数乘积之和。系数见图：

从而可求出\(a_{i,j}\)和\(b_i\)。

注意到上述计算中系数只与顶点的度有关，且大部分情况下度都为6，于是\(w_i\)的值可以预先计算出来。

于是新的Loop细分提升小波算法为 \[ \begin{cases} v_i\gets v_i+\omega_ie\\ e\gets e+\frac38(v_0+v_1)+\frac18(v_2+v_3)\\ v\gets \gamma_n v\\ v_i\gets v_i+\delta_n e \end{cases} \]

\[ \begin{cases} v_i\gets v_i-\delta_n e\\ v\gets \frac{1}{\gamma_n} v\\ e\gets e-\frac38(v_0+v_1)-\frac18(v_2+v_3)\\ v_i\gets v_i-\omega_ie\\ \end{cases} \]