As-Rigid-As-Possible Shape Interpolation(ASAP)



0 尽量刚体的平面形状插值

代码也实现了,把理解的过程写下来。

  • 第一是试试写md和MJ;
  • 第二是记录一下看文章的思路,公式怎么推导的,怎么转换成代码的。

1 Introduction

As-Rigid-As-Possible(ASAP) Shape Interpolation 首先说说文章做什么:给定两个平面网格(一般是同构的),平面形状插值就是计算两个网格之间的变换,如何从一个网格平滑地变换到另一个网格。例如下面论文中给出的变换,输入都是开始和结束两个网格,中间过程都是通过插值得到的。

我们很自由地会想到线性插值——很简单,将网格中对应点按照连线每次移动一个百分比的距离,就实现了最简单的差值。但是效果很差:比如下面a)的大象与长颈鹿的变换,大象鼻子就在变换的过程发生了萎缩,我们更希望得到的是b)这种带有弧线的变换,看起来更自然。本文就是说明如何实现曲线变换啦。

网格变换示例
图1 网格变换示例

2 Transforming Shapes

2.1 单个三角形的变换

首先从单个三角形的变换开始理解,再推广到一个三角形网格上

2.1.1 变换方程

定义源三角形顶点$P=(\vec{p_1},\vec{p_2},\vec{p_3})$和目标三角形$Q=(\vec{q_1},\vec{q_2},\vec{q_3})$,两个三角形的顶点序号相对应。使得$P$变形为$Q$的仿射变换用矩阵$A$表示: $$ A\vec{p_i}+\vec{l}= \begin{pmatrix} a_1 & a_2 \\ a_3 & a_4 \end{pmatrix} \vec{p_i}+ \begin{pmatrix} l_x \\ l_y \end{pmatrix} =\vec{q_i},\quad i\in\lbrace1,2,3\rbrace\label{1}\tag{1} $$

  • 矩阵$A$的作用是将三角形$P$进行一个仿射变换
  • 矩阵$l$的作用是将三角形移动到新的位置

我们不考虑$l$矩阵,因为该矩阵对形状是没有贡献的,只对位置有变化作用。我们最关心的是对旋转、放缩有关系的$A$矩阵

2.1.2 $A$矩阵

假设三角形的中间插值形状为$V(t)=(\vec{v_1}(t),\vec{v_2}(t),\vec{v_3}(t))$,定义$V(t)=A(t)P$,这里需要对$A$矩阵进行一个插值,不然就直接得到了$Q$结果。那么要怎么设计$A$矩阵的插值函数才是有意义的?$A(t)$应该要使得变形理想,作者就列出下面的一些变换要求:

  • 随着$t$的变换必须是对称的
  • 旋转角和放缩比应该是线性的
  • 三角形在旋转过程不能翻转
  • 点的变换路径必须是简单的

作者认为,最基础的想法就是对$A$矩阵分解,使之含有旋转和正放缩因子部分。那么最自然的就是SVD分解: $$ A=R_{\alpha}DR_{\beta}=R_{\alpha} \begin{pmatrix} s_x & 0 \\ 0 & s_y \end{pmatrix} R_{\beta}, \quad s_x,s_y\gt0\label{2}\tag{2} $$ 随后,在作者尝试多种试验后,他们发现Shoemake提出了一种基于SVD的分解表现最好,将上述矩阵分成两部分,一个是单一的旋转阵,另一个是对称阵: $$ A=R_{\alpha}DR_{\beta}=R_{\alpha}(R_{\beta}R_{\beta}^T)DR_{\beta}=\\ (R_{\alpha}R_{\beta})(R_{\beta}^{T}DR_{\beta})= R_{\gamma}S=R_{\gamma} \begin{pmatrix} s_x & s_h \\ s_h & s_y \end{pmatrix}\label{3}\tag{3} $$ 得到两种分解以后就是对公式$\ref{2}$、$\ref{3}$进行插值做实验,它们的插值公式如下: $$ 对于\ref{2}:A_{\alpha,\beta}=R_{t\alpha}((1-t)I+tD)R_{t\beta}, \quad t\in [0,1]\label{4}\tag{4} $$ $$ 对于\ref{3}:A_{\gamma}=R_{t\gamma}((1-t)I+tS), \quad t\in [0,1]\label{5}\tag{5} $$

  • 这里要解释一下$R_{t\gamma}$的含义: $$ R_{\gamma}= \begin{pmatrix} a & b \\ b & a \end{pmatrix}= \begin{pmatrix} cos\theta & sin\theta \\ sin\theta & cos\theta \end{pmatrix} $$ $$ R_{t\gamma}= \begin{pmatrix} cos(t\theta) & sin(t\theta) \\ sin(t\theta) & cos(t\theta) \end{pmatrix}, \quad t\in [0,1]\label{6}\tag{6} $$
  • 即,我们得到$R_{\gamma}$矩阵后,要根据矩阵其中一个值反求出$\theta$然后对$\theta$进行线性插值,才能得到$R_{t\gamma}$

2.1.3 单个三角形变形实验

最简单的就是线性变换了$A(t)=(1-t)I+tA$,然后我们把$\ref{2}$和$\ref{3}$结果进行对比如下:

图2 4种变换结果

图中红色的三角形是源或目标三角形,黑色线条是对应顶点的移动轨迹。

  • 图2a就是线性插值的结果;作者首先实验的是该结果,当然这个结果是很差的;
  • 然后通过SVD进行了图2b的实验,使用的是公式$\ref{4}$。但是从轨迹上来看,旋转角度过大,超过了$\pi$,然后作者在$\ref{4}$的基础上进行了一个减$2\pi$的处理,得到了图2c的结果,不过这个旋转轨迹更加诡异了,不可行。
  • 最终采用公式$\ref{5}$的结果得到了一个比较好的旋转轨迹,图2d。

2.2 三角网格

2.2.1 能量函数

  • 那么对一个三角网格中的每一个三角形都做一个这种变换就好了吗?这肯定不是的,因为假设两个三角形共享一个顶点,对于这个共享点,一个三角形顶点的轨迹一般不会适用于另一个三角形。也就是说,假设我对这两个三角网格分别计算插值过程,然后显示出来,那个共享点就很可能会在插值过程中分离。

来明确一些符号:现在不是考虑一个三角形了,而是考虑三角网格$\tau =\lbrace T_{\lbrace i,j,k \rbrace}\rbrace$。每一个源三角形$P_{\lbrace i,j,k \rbrace}=(\vec{_{p_{i}}},\vec{_{p_{j}}},\vec{_{p_{k}}})$对应一个目标三角形$Q_{\lbrace i,j,k \rbrace}=(\vec{_{q_{i}}},\vec{_{q_{j}}},\vec{_{q_{k}}})$,像2.1中一样,对每一个三角形计算公式$\ref{5}$的映射$A_{\lbrace i,j,k \rbrace}(t)$。大多数的顶点对应于不止一个三角形,如上所说,对所有三角形使用理想的映射$A_{\lbrace i,j,k \rbrace}(t)$,那么网格在变形过程一般都不会形成一个完整的网格(也就是许多分离的三角形)。 令 $$ V(t)=(\vec{v_{1}}(t),...,\vec{v_{n}}(t)),\quad t\in[0,1]\label{7}\tag{7} $$ 为每个顶点(顶点总数$n$)的理想路径,当然$V(0)$就是源三角网格,$V(1)$就是目标三角网格。 定义一个源三角形$P_{\lbrace i,j,k \rbrace}(t)$到其插值过程$\vec{v_{i}}(t),\vec{v_{j}}(t),\vec{v_{k}}(t)$的实际变换仿射矩阵为$B_{\lbrace i,j,k \rbrace}(t)$,即一个实际三角形变换为: $$ B_{\lbrace i,j,k \rbrace}(t)\vec{p_{f}}+\vec{l}=\vec{v_f}(t),\quad f\in\lbrace i,j,k \rbrace\label{8}\tag{8} $$

  • 公式$\ref{8}$与$\ref{1}$的意义是一样的,只是$A$是理想变换矩阵,$B$是实际变化呢矩阵。那么要怎么计算这个实际变换矩阵?

文章给出的方法就是计算所有$A-B$的二范数($\begin{Vmatrix}\cdot\end{Vmatrix}$)的和,因为$B$是待定的,那么必然有一个解,使得这个这个和最小,即求下式的最小值: $$ E_{V(t)}=\sum_{\lbrace i,j,k \rbrace\in \tau} \begin{Vmatrix} A_{\lbrace i,j,k \rbrace}(t)-B_{\lbrace i,j,k \rbrace}(t) \end{Vmatrix}^2\label{9}\tag{9} $$ $$ E_{V(t)}=u^T \begin{pmatrix} c & G^T \\ G & H \end{pmatrix}u,\\ 其中u^T=(1,v_{2x}(t),v_{2y}(t),...v_{nx}(t),v_{ny}(t))\label{10}\tag{10} $$

  • 下面重点说明怎么从公式$\ref{9}$变成公式$\ref{10}$:根据公式$\ref{8}$,可以将$B$表示为中间点的函数:

$$\begin{align}&B_{\lbrace i,j,k \rbrace}(t)\vec{p_{f}}+\vec{l}=\vec{v_f}(t)\\&\Rightarrow\begin{pmatrix}b_1 & b_2\\b_3 & b_4\end{pmatrix}\begin{pmatrix}p_{1x} & p_{2x} & p_{3x}\\p_{1y} & p_{2y} & p_{3y}\end{pmatrix}+\begin{pmatrix}l_x & l_x & l_x\\l_y&l_y& l_y\end{pmatrix}=\begin{pmatrix}v_{1x} & v_{2x} & v_{3x}\\v_{1y} & v_{2y} & v_{3y}\end{pmatrix}\\&\Rightarrow\begin{pmatrix}p_{1x} & p_{1y} & 1 \\p_{2x} & p_{2y} & 1 \\p_{3y} & p_{3y} & 1\end{pmatrix}\begin{pmatrix}b_1 & b_2 \\b_3 & b_4 \\l_x & l_y\end{pmatrix}=\begin{pmatrix}v_{1x} & v_{1y} \\v_{2x} & v_{2y} \\v_{3y} & v_{3y}\end{pmatrix}\\&\Rightarrow x=A^{-1}B=\begin{pmatrix}p_{1x} & p_{1y} & 1 \\p_{2x} & p_{2y} & 1 \\p_{3y} & p_{3y} & 1\end{pmatrix}^{-1}B=\begin{pmatrix}h_1 & h_2 & h_3 \\h_4 & h_5 & h_6 \\h_7 & h_8 & h_9\end{pmatrix}B\\&\Rightarrow\begin{pmatrix}b_1 & b_2 \\b_3 & b_4 \\l_x & l_y\end{pmatrix}=\begin{pmatrix}h_1 & h_2 & h_3 \\h_4 & h_5 & h_6 \\h_7 & h_8 & h_9\end{pmatrix}\begin{pmatrix}v_{1x} & v_{1y} \\v_{2x} & v_{2y} \\v_{3y} & v_{3y}\end{pmatrix}\\&\Rightarrow\begin{cases}b_1=h_1v_{1x}+h_2v_{2x}+h_3v_{3x} \\b_2=h_1v_{1y}+h_2v_{2y}+h_3v_{3y} \\b_3=h_4v_{1x}+h_5v_{2x}+h_6v_{3x} \\b_4=h_4v_{1y}+h_5v_{2y}+h_6v_{3y} \end{cases}\end{align}$$ 得到了$B$的以中间点为变量的表达式,对于一个具体的时间$t$,$A$可由$\ref{5}$得到。 则可以计算$\ref{9}$,其中 $$ \begin{Vmatrix} A_{\lbrace i,j,k \rbrace}(t)-B_{\lbrace i,j,k \rbrace}(t) \end{Vmatrix}^2= \begin{Vmatrix} a_1-b_1 & a_2-b_2 \\ a_3-b_3 & a_4-b_4 \end{Vmatrix}^2 =\sum_{i=1}^{4}{(a_{i}-b_{i})^2}\label{11}\tag{11} $$

  • 上式包含6个未知数$v_{1x...3x}\,v_{1y...3y}$,我们假定$(p_{1x},p_{1y})$到$(q_{1x},q_{1y})$的移动是线性插值,则根据具体的$t$就可以得到$v_{1x}、v_{1y}$接着就可以求得剩下四个未知数,这也就是为什么$u^T$的第一项是1的原因。
  • 公式$\ref{9}$其实就是一个关于$v$的二次多项式,那么肯定就能写成一个二次型$\ref{10}$。
  • 从编程的角度来说,如何把二次多项式写成矩阵的形式?对所有的三角形计算出二次型的各个系数,例如我们将$(a_1-b_1)^2$展开,那么必然会产生$cv_{2x}v_{3x}$一项,那么这一项怎么填写进矩阵?以$(x+y)^2$为例: $$ (x+y)^2=x^2+xy+yx+y^2= \begin{pmatrix} x & y \end{pmatrix} \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} $$ $$ \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix}理解为 \begin{array}{c|lcr} & \text{x} & \text{y} \\ \hline x & 1 & 1 \\ y & 1 & 1 \end{array} $$ 平方项的系数加在对角线上;含有两个未知数的项,其系数折半加在矩阵对称位置。
  • 因为一个点例如$v5$可能会被多个三角形所包含,所以在矩阵的某些位置就会被加上许多的值,这个矩阵填写完成后,我们就可以进行下一步的计算。 实际上,我也是这样写的:方法比较笨拙,但是结果还算正确
    for (size_t i = 0; i != src_faces.size(); ++i){
    //……
    M(0, 0) += Ar(0, 0)*Ar(0, 0) + Ar(1, 0)*Ar(1, 0);
    M(F(pm[0] * 2 - 1), F(pm[0] * 2 - 1)) += (A_n(0, 0)*A_n(0, 0) + A_n(1, 0)*A_n(1, 0))*(pm[0] == 0 ? u0.at(0)*u0.at(0) : 1);//square??
    //……
    M(F(pm[0] * 2 - 1), 0) += -(Ar(0, 0)*A_n(0, 0) + Ar(1, 0)*A_n(1, 0))*(pm[0] == 0 ? u0.at(0) : 1);
    M(0, F(pm[0] * 2 - 1)) += -(Ar(0, 0)*A_n(0, 0) + Ar(1, 0)*A_n(1, 0))*(pm[0] == 0 ? u0.at(0) : 1);
    //……
    M(F(pm[0] * 2 - 1), F(pm[1] * 2 - 1)) += (A_n(0, 0)*A_n(0, 1) + A_n(1, 0)*A_n(1, 1))*(pm[0] == 0 ? u0.at(0) : 1)*(pm[1] == 0 ? u0.at(0) : 1);
    //……
    M(0, 0) += Ar(0, 1)*Ar(0, 1) + Ar(1, 1)*Ar(1, 1);
    M(F(pm[0] * 2), F(pm[0] * 2)) += (A_n(0, 0)*A_n(0, 0) + A_n(1, 0)*A_n(1, 0))*(pm[0] == 0 ? u0.at(1) : 1);
    //……
    M(F(pm[0] * 2), 0) += -(Ar(0, 1)*A_n(0, 0) + Ar(1, 1)*A_n(1, 0))*(pm[0] == 0 ? u0.at(1) : 1);
    M(0, F(pm[0] * 2)) += -(Ar(0, 1)*A_n(0, 0) + Ar(1, 1)*A_n(1, 0))*(pm[0] == 0 ? u0.at(1) : 1);
    //……
    M(F(pm[0] * 2), F(pm[1] * 2)) += (A_n(0, 0)*A_n(0, 1) + A_n(1, 0)*A_n(1, 1))*(pm[0] == 0 ? u0.at(1) : 1)*(pm[1] == 0 ? u0.at(1) : 1);
    //……
    }
    

    2.2.2 能量函数求导

    对于公式$\ref{10}$,其求导等于零的结果为: $$ HV=H \begin{pmatrix} v_{2x}(t) \\ v_{2y}(t) \\ \vdots \end{pmatrix} =-G\label{12}\tag{12} $$

  • 下面来简单推导一下这个公式的由来:

令$F=WMW^T$,$H_{i}$表示$H$矩阵的第$i$列

$$ \begin{align} F & = WMW^T=(1,v_{2x},v_{2y},\cdots,v_{2n-1},v_{2n}) \begin{pmatrix} c & g_1 & g_2 & \cdots & g_{2n-3} & g_{2n-2} \\ g_1 \\ g_2 & \\ \cdots & & & H \\ g_{2n-3} \\ g_{2n-2} \end{pmatrix}W^T\\ & =(c+g_{1}v_{2x}+g_{2}v_{2y}+\cdots+g_{2n-3}v_{2n-1}+g_{2n-2}v_{2n},\\ & \qquad g_{1}+VH_{1},g_{2}+VH_{2},\cdots,g_{2n-3}VH_{2n-3}+g_{2n-2}+VH_{2n-2})W^T\\ & = c+g_{1}v_{2x}+g_{2}v_{2y}+\cdots+g_{2n-3}v_{2n-1}+g_{2n-2}v_{2n}+\\ & \qquad (g_{1}+VH_{1})v_{2x}+(g_{2}+VH_{2})v_{2y}+\cdots+(g_{2n-3}VH_{2n-3})v_{2n-1}+(g_{2n-2}VH_{2n-2})v_{2n}\\ \Rightarrow F' & = g_{1}+g_{2}+\cdots+g_{2n-3}+g_{2n-2}+\\ & \qquad g_{1}+2v_{2x}H_{1}+g_{2}+2v_{2y}H_{2}+\cdots+g_{2n-3}+2v_{2n-3}H_{2n-3}+g_{2n-2}+2v_{2n-2}H_{2n-2}\\ & = 2G+2VH = 0\\ \Rightarrow V & = -GH' \end{align} $$

$H$和$G$矩阵在公式$\ref{10}$中填写完毕,求一个线性方程组即可得到一个时间$t$的中间结果。要想得到精细的路径,就缩小每次跨越的步长,即把$[0,1]$分成更多份。

  • 矩阵$H$对所有时刻$t$都是一样的,所以只需要求一次该矩阵的逆,然后再计算一次矩阵乘法$\ref{12}$就可以大大加快计算速度。

3 Expriment

3.1 运行结果

3.1 代码

代码可以采用Cmake生成Visual Studio等的工程文件,但是需要自己配置Qt5

为了配置简单,我把boost库的内容去掉了, 如果需要多线程,可以增加boost相关的东西。 运行界面如下:

  1. 文件中含有同构剖分的文件夹,点击FILE->Open,进入duck文件夹,依次打开两个文件。
  2. 点击RUN->ARAP运行

3 thoughts on “As-Rigid-As-Possible Shape Interpolation(ASAP)”

Leave a Reply

Your email address will not be published. Required fields are marked *

This site uses Akismet to reduce spam. Learn how your comment data is processed.