Carol's Blog
返回

Mesh Deformation Transfer

Deformation Transfer For Triangle Mesh

概述

问题描述

对于两个具有一定视觉相似度的三角网格：原网格 S0S_0 和目标网格 T0T_0，根据原网格已知的变形序列 S={S0, …, S∣S∣}\mathcal S=\{S_0,\,\dots,\, S_{|\mathcal S|}\}，生成目标网格的对应的变形序列T={T0, …,T∣S∣}\mathcal T=\{T_0,\, \dots, T_{|\mathcal S|}\}。

banner

问题假设

  1. S\mathcal S 和 T\mathcal T 中的网格分别具有相同的拓扑结构，但两个集合对应网格之间不要求。
  2. S0S_0 与 T0T_0 应当具有一定的视觉相似度，并且相关点对通过人为标记的方式体现。

基本思想

  1. 根据人为标定的标记点，计算由 S0S_0 到 T0T_0 的三角面片对应关系 M={(s0, t0), …, (s∣M∣, t∣M∣)}M=\{(s_0,\,t_0),\, \dots,\,(s_{\lvert M\rvert},\, t_{\lvert M\rvert})\}。
  2. 对 Si, i∈{1, …, ∣S∣}S_i,\, i\in \{1,\, \dots ,\, |\mathcal S|\} ，由于 SiS_i 与 S0S_0 具有相同的拓扑，可以根据对应关系MM，将SiS_i每个三角面的变换作用到T0T_0，加上一些约束条件，得到变形后的 TiT_i 。

三角面变换

在三维空间中，对于三角面Fi=[v0, v1, v2]F_i = \left[v_0,\, v_1,\, v_2\right] 到另一个三角面Fi~=[v~0, v~1, v~2]\widetilde{F_i} =\left[\widetilde v_0,\, \widetilde v_1,\, \widetilde v_2 \right] 的仿射变换可以分解为线型部分QiQ_i 和非线性部分 di\mathbf d_i。求解该仿射变换需要用四对点到点的关系，对每个三角面引入第四个点：

v3=v0+(v1−v0)×(v2−v0)/∣(v1−v0)×(v2−v0)∣v_3 = v_0 +(v_1 −v_0) \times (v_2 −v_0) \big / \sqrt{|(v_1 −v_0) \times (v_2 −v_0)|}

线性部分Qi=V~iVi−1Q_i = \widetilde V_i {V_i}^{-1}，其中Vi=[v1−v0v2−v0v3−v0]V_i = \left[v_1 - v_0\quad v_2 - v_0 \quad v_3-v_0\right]

对应关系计算

根据标定点，将 S0S_0 在保持拓扑不变的前提下变形为 T0T_0。

根据三角面变换中的定义，将S0S_0中每个三角面的第四个点加到顶点序列之后：

x=(v0, v1, …, vn−1⏟原始顶点, vn, …, vn+m−1⏟新增顶点)x = \left(\underbrace{v_0,\, v_1,\, \dots,\, v_{n-1}}_{\text{原始顶点}},\, \underbrace{v_{n},\, \dots,\, v_{n+m-1}}_{\text{新增顶点}} \right)

其中nn为S0S_0中原始顶点的个数，mm为三角面个数，x∈R3×(n+m)x \in \mathbb R^{3\times (n+m)}。

通过改动xx，来实现S0S_0到T0T_0的变形，具体表现为定义损失函数，最小二乘法搜索最优解：

min⁡x~(wSES, wIEI, wCEC)约束标记点对中原网格顶点与目标网格顶点相同\begin{aligned} \min_{\widetilde x} (w_S E_S,\, w_IE_I,\, w_CE_C) \\ \text{约束标记点对中原网格顶点与目标网格顶点相同} \end{aligned}
损失项表达式备注
平滑性
(smoothness)		ES(x)=∑Qi∑Qj∈adj(Qi)∥Qi−Qj∥F2E_S(x) = \sum_{Q_i}{ \sum_{Q_j \in \text{adj}(Q_i)}{\lVert Q_i -Q_j\rVert_F^2}}对每一个三角面，周围三角面的形变应该尽量与之相似
不变性
(identity)		EI(x)=∑Qi∥Qi−I∥F2E_I(x) = \sum_{Q_i}{\lVert Q_i - \mathbf I \rVert_F^2}每个三角面的形变尽可能小
最近点损失EC(x; c0, …, cn)=∑i∥vi−ci∥F2E_C(x;\, c_0,\,\dots,\,c_n)=\sum_i{\lVert v_i -c_i\rVert^2_F}原网格的每个顶点都应该尽可能贴近与目标网格的最近点

当S0S_0在保持拓扑不变的情况下变形为接近 T0T_0后，计算两个网格三角面之间的对应关系MM。

对S0S_0中的每个三角面，在T0T_0中至少存在一个三角面与之存在最近关系，反之亦然。所以两个网格中三角面的对应关系为多对多，并非映射。

形变迁移

在有原网格与目标网格三角面的对应关系MM, 以及原网格 初始网格S0S_0 与形变网格SiS_i 的每个三角面形变关系后，我们可以直接将S0S_0 到SiS_i的形变迁移到T0T_0 上，得到TiT_i：

min⁡Qi+di∑j=0∣M∣−1∥Ssj−Ttj∥F2约束 Qj vi+dj=Qk vi+dk,∀i,∀j,k∈adj(vi)\begin{aligned} &\min_{Q_i+\mathbf d_i} \sum_{j=0}^{ | M | -1} { \lVert \mathbf S_{s_j}- \mathbf T_{t_j}\rVert_F^2} \\ &\text{约束 } Q_j\, v_i + \mathbf d_j = Q_k \, v_i + \mathbf d_k, \\ &\forall i, \\ &\forall j,k \in \text{adj}(v_i) \end{aligned}

这里的Ssj\mathbf S_{s_j}表述为网格SiS_i中标号为sjs_j的三角形的变换QsjQ_{s_j}，TtjT_{t_j}同理。

对于实际求解，可以将上述对三角面变换的逼近转换到对顶点变换的逼近。

由于三角面变换逼近时，可能会出现边缘撕裂的情况，所以需要添加三角形的邻域约束。而转用顶点表达时，通过同一网格类型拓扑不变的假设，蕴含了变换后的网格不会出现边缘撕裂。

具体的顶点公式推导在后面。

细节推导

三角面变换

以上描述的优化函数大多都是用三角面的形变表示，而优化目标是顶点序列xx，需要进行推导将三角面的形变转换为关于顶点序列 xx 的表达式，即。

Qi=V~i Vi−1→展开顶点序列 xQi=x V^i−1Q_i = \widetilde V_i\, {V_i}^{-1} \quad \xrightarrow[]{\text{展开顶点序列 }x} \quad Q_i= x\,\hat{V}_i^{-1}

计算展开xx

Qi3×3=V~i Vi−1=[v~1i−v~0iv~2i−v~0iv~3i−v~0i] Vi−1=[v~1iv~2iv~3i] Vi−1−[v~0iv~0iv~0i] Vi−1=[v~0v~1…v~n+m−1][…000…100…010…001…]Vi−1−[v~0v~1…v~n+m−1][…111…000…000…000…]Vi−1=[v~0v~1…v~n+m−1][…−1−1−1…100…010…001…]Vi−1=x V^i−1\begin{align*} \underset{3\times 3}{Q_i} &= \widetilde V_i\, {V_i}^{-1} \\[2ex] &= \begin{bmatrix} \widetilde v_1^i - \widetilde v_0^i & \widetilde v_2^i - \widetilde v_0^i & \widetilde v_3^i - \widetilde v_0^i \end{bmatrix} \, {V_i}^{-1} \\[2ex] &= \begin{bmatrix} \widetilde v_1^i & \widetilde v_2^i & \widetilde v_3^i \end{bmatrix}\, {V_i}^{-1} - \begin{bmatrix} \widetilde v_0^i & \widetilde v_0^i & \widetilde v_0^i \end{bmatrix}\, {V_i}^{-1}\\[2ex] &= \begin{bmatrix} \widetilde v_0 & \widetilde v_1 & \dots & \widetilde v_{n+m-1} \end{bmatrix} \begin{bmatrix}&\dots\\0 & 0 & 0\\ & \dots\\1 &0&0\\&\dots&\\0&1&0\\&\dots&\\0&0&1\\&\dots\end{bmatrix} {V_i}^{-1} - \begin{bmatrix} \widetilde v_0 & \widetilde v_1 & \dots & \widetilde v_{n+m-1} \end{bmatrix} \begin{bmatrix}&\dots\\1&1&1\\& \dots\\0&0&0\\&\dots&\\0&0&0\\&\dots&\\0&0&0\\&\dots\end{bmatrix} {V_i}^{-1} \\[2ex] &=\begin{bmatrix} \widetilde v_0 & \widetilde v_1 & \dots & \widetilde v_{n+m-1} \end{bmatrix} \begin{bmatrix}&\dots\\-1&-1&-1\\& \dots\\1&0&0\\&\dots&\\0&1&0\\&\dots&\\0&0&1\\&\dots\end{bmatrix} {V_i}^{-1} \\[2ex] &= x \, \hat{V}_i^{-1} \end{align*}

一般的目标函数

对于一般情况，期望在保持原网格拓扑不变的同时满足每个三角面的变形目标CiC_i，写作：

∑i=0∣M∣−1∥Qi−Ci∥F2=∑i=0∣M∣−1∥x V^i−1−Ci∥F2=∑i=0∣M∣−1∥(x V^i−1−Ci)T∥F2=∑i=0∣M∣−1∥(V^i−1)T xT−CiT∥F2=∥((V^0−1)T(V^1−1)T…(V^∣M∣−1−1)T) xT−(C0TC1T…C∣M∣−1T)∥F2=∥A xT−b∥F2\begin{align*} \sum_{i=0}^{\lvert M \rvert-1}{\big\lVert Q_i -C_i\big\rVert_F^2} &= \sum_{i=0}^{\lvert M \rvert-1}{\big\lVert x\, \hat{V}_i^{-1} - C_i\big\rVert_F^2} \\[2ex] & =\sum_{i=0}^{\lvert M \rvert-1}{ \begin{Vmatrix} \left( x\, \hat{V}_i^{-1} - C_i\right)^T \end{Vmatrix}_F^2 } \\[2ex] &= \sum_{i=0}^{\lvert M \rvert-1}{ \begin{Vmatrix} \left(\hat{V}_i^{-1}\right)^T\, x^T-C_i^T \end{Vmatrix}_F^2 } \\[2ex] &= \begin{Vmatrix} \begin{pmatrix} \left(\hat{V}_0^{-1}\right)^T \\ \left(\hat{V}_1^{-1}\right)^T \\ \dots\\ \left(\hat{V}_{\lvert M \rvert-1}^{-1}\right)^T \end{pmatrix} \,x^T - \begin{pmatrix} C_0^T \\ C_1^T \\ \dots \\ C_{\lvert M \rvert-1}^T \end{pmatrix} \end{Vmatrix} _F^2 \\[2ex] &=\begin{Vmatrix} A\,x^T - b \end{Vmatrix}_F^2 \end{align*}

即目标函数变为 E(x)=∥A xT−b∥F2E (x) = \begin{Vmatrix} A\, x^T - b \end{Vmatrix}_F^2 的形式，使用最小二乘法1，得到

∂E(x)∂x=∂(xTATA x−2bTA x+bT b)∂x=2ATA x−2ATb=0ATA x=ATb\begin{aligned} \frac{\partial E(x)}{\partial x} &= \frac{\partial {\left( x^TA^TA\,x -2b^TA\,x+b^T\,b\right)} }{\partial x} = 2A^TA\,x - 2A^T b = 0 \\[2ex] A^T A\,x &= A^T b \end{aligned}

具体的目标函数

平滑性损失

ES(x)=∑Qi∑Qj∈adj(Qi)∥Qi−Qj∥F2=∥((V^0−1−V^j0−1)T(V^0−1−V^j1−1)T…(V^0−1−V^j∣adj(Qi)∣−1−1)T(V^1−1−V^j0−1)T…) xT−0∥F2=∥AS xT−bS∥F2\begin{align*} E_S(x)&=\sum_{Q_i}{ \sum_{Q_j \in \text{adj}(Q_i)}{\lVert Q_i -Q_j\rVert_F^2}} \\[2ex] &= \begin{Vmatrix} \begin{pmatrix} \left(\hat{V}_0^{-1} - \hat{V}_{j_0}^{-1}\right)^T \\ \left(\hat{V}_0^{-1} - \hat{V}_{j_1}^{-1}\right)^T \\ \dots \\ \left(\hat{V}_0^{-1} - \hat{V}_{j_{\lvert \text{adj}(Q_i) \rvert-1}}^{-1}\right)^T \\ \left(\hat{V}_1^{-1} - \hat{V}_{j_0}^{-1}\right)^T\\ \dots \end{pmatrix} \, x^T - 0 \end{Vmatrix}_F^2 = \begin{Vmatrix} A_S\, x^T - b_S \end{Vmatrix}_F^2 \end{align*}

其中AS∈R3q×(n+m)A_S \in \mathbb R^{3q\times {(n+m)}}， bS∈R3q×3b_S \in \mathbb R^{3q\times 3}，q=∑i∣adj(Qi)∣q = \sum_i {\lvert \text{adj}(Q_i)\rvert}。

不变性损失

EI(x)=∑Qi∥Qi−I∥F2=∥((V^0−1)T(V^1−1)T…(V^m−1−1)T) xT−(II…I)∥=∥AI xT−bI∥F2\begin{align*} E_I(x) &= \sum_{Q_i}{\lVert Q_i - \mathbf I \rVert_F^2}\\[2ex] &= \begin{Vmatrix} \begin{pmatrix} \left(\hat{V}_0^{-1}\right)^T \\ \left(\hat{V}_1^{-1}\right)^T \\ \dots\\ \left(\hat{V}_{m-1}^{-1}\right)^T \end{pmatrix} \,x^T - \begin{pmatrix} \mathbf I \\ \mathbf I \\ \dots \\ \mathbf I \end{pmatrix} \end{Vmatrix} = \begin{Vmatrix} A_I\, x^T - b_I \end{Vmatrix}_F^2 \end{align*}

其中AI∈R3m×(n+m)A_I \in \mathbb R^{3m \times (n+m)}，bI∈R3m×3b_I \in \mathbb R^{3m\times 3}。

最近点损失

EC(x; c0, …, cn)=∑i∥vi−ci∥2=∥I xT−(C0TC1T…Cm+n−1T)∥F2=∥AC xT−bC∥F2\begin{align*} E_C(x;\, c_0,\,\dots,\,c_n)&=\sum_i{\lVert v_i -c_i\rVert^2} \\[2ex] &=\begin{Vmatrix} \mathbf I\, x^T - \begin{pmatrix} C_0^T \\ C_1^T \\ \dots \\ C_{m+n-1}^T \end{pmatrix} \end{Vmatrix}_F^2 = \begin{Vmatrix}A_C \, x^T -b_C \end{Vmatrix}_F^2 \end{align*}

其中AC∈R3m×(n+m)A_C \in \mathbb R^{3m \times (n+m)}，bC∈R3m×3b_C \in \mathbb R^{3m\times 3}。

为了统一使用变量xTx^T，需要有n+mn+m个点的形式，但实际上只需要计算前nn个顶点的最近点。

形变迁移损失

EQ(x)=∑j=0∣M∣−1∥Ssj−Ttj∥F2=∑j=1∣M∣∥SsjT−V^tj−1T xT∥F2=∥−((V^t0−1)T(V^t1−1)T…(V^t∣M∣−1−1)T) xT+(Ss0TSs1T…Ss∣M∣−1T)∥F2=∥AQ xT−bQ∥F2\begin{align*} E_Q(x)&=\sum_{j=0}^{\lvert M\rvert -1} {\lVert \mathbf S_{s_j}- \mathbf T_{t_j}\rVert_F^2}\\[2ex] &=\sum_{j=1}^{\lvert M\rvert}{ \begin{Vmatrix} \mathbf S_{s_j}^T - {\hat{V}_{t_j}^{-1}}^T\, x^T \end{Vmatrix}_F^2 } \\[2ex] &= \begin{Vmatrix} -\begin{pmatrix} \left(\hat{V}_{t_0}^{-1}\right)^T \\ \left(\hat{V}_{t_1}^{-1}\right)^T \\ \dots\\ \left(\hat{V}_{t_{\lvert M \rvert-1}}^{-1}\right)^T \end{pmatrix} \, x^T +\begin{pmatrix} \mathbf S_{s_0}^T \\ \mathbf S_{s_1}^T \\ \dots \\ \mathbf S_{s_{\lvert M\rvert -1}}^T \end{pmatrix} \end{Vmatrix}_F^2\\[2ex] &= \begin{Vmatrix} A_Q\, x^T - b_Q \end{Vmatrix}_F^2 \end{align*}

其中AQ∈R3∣M∣×(n+m)A_Q \in \mathbb R^{3\vert M\rvert \times (n+m)}，bQ∈R3∣M∣×3b_Q \in \mathbb R^{3\vert M\rvert \times 3}，∣M∣\lvert M \rvert是对应关系中元素的个数，满足∣M∣≥m\lvert M\rvert \ge m。

实验

注意事项

对应关系求解: 标记点约束

在求解min⁡xES\min_x E_S、min⁡xEI\min_xE_I和min⁡xEQ\min_x E_Q时，需要约束S0S_0和T0T_0对应的标记点相同：

AT A xT=AT b(Au)T Au (xu)T+(Am)T Am (x~m)T=AT b(Au)T Au (xu)T=AT b−(Am)T Am (x~m)T\begin{align*} A^T\,A\, x^T &= A^T\, b\\[2ex] ({A^u})^T\,A^u\, \left(x^u\right)^T + ({A^m})^T\, A^m\, \left(\widetilde{x}^m\right)^T &= A^T\, b\\[2ex] ({A^u})^T\,A^u\, \left(x^u\right)^T &= A^T\,b - ({A^m})^T\,A^m\, \left(\widetilde{x}^m\right)^T \end{align*}

对于S0S_0，AuA^u、xux^u分别为未标记点的系数矩阵以及顶点序列，AmA^m、xmx^m为已标记点的系数矩阵以及顶点序列，x~m\widetilde{x}^m为T0T_0中已标记点的顶点序列。

再求解出xux^u后再根据约束条件xm=x~mx^m = \widetilde{x}^m 计算出结果 x~\widetilde{x}。

对应关系求解: 目标函数组合与迭代

对应关系计算时，可以将相关目标函数写在一起：

x~=arg⁡min⁡x∥(wSASwIAIwCAC) xT−(wSbSwIbIwCbC)∥F2\widetilde{x} = \underset{x}{\arg \min} \begin{Vmatrix} \begin{pmatrix} w_SA_S \\ w_IA_I\\ w_CA_C \end{pmatrix} \, x^T - \begin{pmatrix} w_S b_S\\ w_I b_I\\ w_C b_C \end{pmatrix} \end{Vmatrix}_F^2

其中wS=1.0w_S=1.0，wI=0.001w_I=0.001，wC=[0, 1 …, 5000]w_C=[0,\,1\, \dots,\, 5000]。

稀疏线性方程组求解2

  1. 直接求解
    • LU分解：Ax=b⟹LUx=bAx=b \Longrightarrow LUx = b
    • Cholesky分解3：Ax=b⟹LTLx=bAx = b \Longrightarrow L^TLx=b
  2. 间接求解
    • Jacobi method / Gauss-Seidel method

优化

  1. 稀疏矩阵的计算：乘法、转置、切片、拼接和索引。
  2. 稀疏方程组求解。
  3. 计算ECE_C时，最近点的计算。
  4. 计算对应关系MM时，最近三角面重心以及法线夹角的计算。

计算流程

计算形变迁移，本质上在优化四个目标函数：ESE_S、EIE_I、ECE_C和EQE_Q。其中系数矩阵几乎都包含V^−1\hat{V}^{-1}的计算。

Footnotes

  1. https://eeweb.engineering.nyu.edu/iselesni/lecture_notes/least_squares/least_squares_SP.pdf ↩

  2. https://zhuanlan.zhihu.com/p/479913328 ↩

  3. https://zhuanlan.zhihu.com/p/112091443 ↩

分享这篇文章：
通过邮件分享这篇文章

分享到微信

微信对普通网页没有开放通用直连分享协议。更稳妥的方式是复制链接、扫码打开，或在支持的设备上调用系统分享。

上一篇
文明与罗马文明
下一篇
曲面参数化概述