Analytic Eigensystems for Isotropic Distortion Energies 论文笔记

GAMES301-曲面参数化(GAMES 301: Surface Parameterization)作业2解析与matlab实现。

Analytic Eigensystems for Isotropic Distortion Energies是来自Pixar动画工作室的一篇文章。这篇文章本身是用来处理变形体仿真的问题，其主要贡献在于使用矩阵的不变量来描述变形能量并且推导了各向同性材料在相应能量下的解析解。除了物理仿真外，这项技术同样可以用来处理曲面参数化的问题。我们可以把初始状态的曲面想象成某种各向同性材料组成的外皮，而参数化的过程相当于把这层皮展开并约束到平面上，此时它会通过释放变形能来降低自身的几何扭曲。这篇笔记会详细推导相关的理论并记录下我的matlab实现。

Background

Deformation Gradient

三角形网格的变形可以使用一个线性函数来进行表示。具体来说，在平面上三角形的每一个顶点的位移都是其坐标的线性组合：

[\begin{matrix} x \\ y \end{matrix}] = [\begin{matrix} F_{11} & F_{12} \\ F_{21} & F_{22} \end{matrix}] [\begin{matrix} \bar{x} \\ \bar{y} \end{matrix}] + [\begin{matrix} t_{x} \\ t_{y} \end{matrix}]

上式的矩阵形式为：

x = F \bar{x} + t

其中系数矩阵 $F$ 即为坐标变换的Jacobian矩阵，也称为变形梯度(deformation gradient)，它描述了三角形的旋转、缩放和反射；而向量 $t$ 则描述了三角形的刚体位移。

变形梯度的计算非常简单，我们只需要考虑三角形两条对应边在变形前后的线性变换即可。我们可以把 ${\bar{x}}_{0}$ 作为三角形的局部原点，此时它的两条邻边向量分别为：

{\bar{o}}_{1} = {\bar{x}}_{1} - {\bar{x}}_{0}

{\bar{o}}_{2} = {\bar{x}}_{2} - {\bar{x}}_{0}

带入变形梯度 $F$ 可以得到：

F [\begin{matrix} {\bar{o}}_{1} | {\bar{o}}_{2} \end{matrix}] = [\begin{matrix} o_{1} | o_{2} \end{matrix}]

F D_{m} = D_{s}

因此只需要求解矩阵方程即可：

F = D_{s} D_{m}^{- 1}

计算三角形从初始构形x1到新构形x2的变形梯度(Jacobian矩阵)可以参考如下findJacobian()函数：

 function J = findJacobian(x1, x2)
%% Find Jacobian matrix from rest pose x1 to current pose x2
%% Args:
%%      x1[3, 2]: rest pose
%%      x2[3, 2]: current pose
%% Return:
%%      J[2, 2]: Jacobian matrix (deformation gradient)

Dm = (x1(2:end, :) - x1(1, :))';
Ds = (x2(2:end, :) - x2(1, :))';

%% J * Dm = Ds
J = Ds * matrixInv2x2(Dm);

end
 

Compute Rest Pose

很多函数都需要用到三角网格的初始构型x1，它是三个顶点在三角形所在平面上定义坐标系下的坐标。要计算这个构型有很多种不同的方法，比如说可以计算空间三角形的三条边长然后重新放置到xy平面上。这里project2Plane()函数实现了利用刚体变换把空间三角形放置到平面上的计算方法，其思路在于构造一个刚体变换使得变换后三角形的中心为原点且法向指向z轴，此时只需提取变换后顶点的xy坐标即可。

 function uv = project2Plane(V)
%% Project a triangle to plane
%% Args:
%%      V[3, 3]: vertex coordinates in 3D, each row is a coordinate vector
%% Returns:
%%      uv[3, 2]: vertex coordinates in 2D

%% center in 3D
c = mean(V, 1);

%% normal
e1 = V(2, :) - V(1, :);
e2 = V(3, :) - V(1, :);
n = cross(e1, e2);
n = n / norm(n);

%% rotation
z = [0, 0, 1];
r = vrrotvec(n, z);

%% project to plane
uv = vrrotvec2mat(r) * (V - c)';
uv = uv(1:2, :)';

end
 

Symmetric Dirichlet Energy

当三角形发生变形时会产生变形能，我们可以定义不同的能量函数来描述变形的大小。显然变形能量与三角形的刚体变换无关，因此它应该是只关于变形梯度 $F$ 的函数，而且不包含其中的旋转成分。这里可以使用极分解(polar decomposition)来分离出变形梯度中的旋转部分：

F = R S

上式中矩阵 $R$ 即为变形梯度的旋转成分；矩阵 $S$ 则表示伸缩，称为拉伸张量(stretch tensor)。极分解的构造非常简单，我们只需要利用矩阵的奇异值分解(singular value decomposition, SVD)就可以实现：

F = U Σ V^{T}

R = U V^{T}

S = V Σ V^{T}

容易验证如此定义的极分解满足所需的性质。不过需要注意的是直接利用SVD计算的旋转梯度 $R$ 可能会出现行列式为负的情况，而旋转矩阵一般要求其行列式为1。因此我们需要对标准的SVD进行一定的修改，保证 $U$ 和 $V$ 行列式同号：

 function [U, Sigma, V] = svd_rv(F)
%% rotation-variant SVD

%% standard SVD
[U, Sigma, V] = svd(F);

%% reflection matrix
%% modified for 2D deformation
L = eye(2,2);
L(2,2) = det(U * V');

%% see where to pull the reflection out of
detU = det(U);
detV = det(V);

if (detU < 0 && detV > 0)
    U = U * L;
elseif (detU > 0 && detV < 0)
    V = V * L;
end

%% push the reflection to the diagonal
Sigma = Sigma * L;

end
 

利用svd_rv()函数可以非常容易地定义极分解：

 function [R, S] = polarDecomposition(F)
%% polar decomposition with rotation-variant SVD

%% rotation variant SVD
[U, Sigma, V] = svd_rv(F);

R = U * V';
S = V * Sigma * V';

end
 

分理出旋转后就可以直接使用 $S$ 矩阵来定义变形能。对称Dirichlet能量(symmetric Dirichlet energy)是几何处理中非常流行的能量函数，它是由变形梯度的两个奇异值来定义的：

\begin{aligned} Ψ_{SD} & = \frac{1}{2} (σ_{1}^{2} + σ_{2}^{2} + \frac{1}{σ_{1}^{2}} + \frac{1}{σ_{2}^{2}}) \\ = \frac{1}{2} (‖ F ‖_{F}^{2} + ‖ F^{- 1} ‖_{F}^{2}) \\ = \frac{1}{2} (‖ S ‖_{F}^{2} + ‖ S^{- 1} ‖_{F}^{2}) \end{aligned}

对称Dirichlet能量的特点在于它包含了奇异值的倒数，因此当三角形接近退化时会产生很大的能量作为惩罚项来避免进一步退化。计算对称Dirichlet能量 $Ψ_{SD}$ 的代码如下：

 function E = computeSDEnergy(x1, x2)
%% Symmetric Dirichlet energy from x1 to x2
%% Args:
%%      x1: rest pose
%%      x2: current pose
%% Return:
%%      E: symmetric Dirichlet energy

F = findJacobian(x1, x2);
Finv = matrixInv2x2(F);

E = norm(F, "fro")^2 + norm(Finv, "fro")^2;
E = 0.5*E;

end
 

而整个三角网格上我们一般会对每个三角形的变形能量按照它原来的面积进行加权：

Ψ = \sum_{q} | q | Ψ_{SD} (F_{q})

其中下标 $q$ 表示每个单独的三角形。此时每个三角形的初始构型x1为单元自身定义的局部坐标系下的坐标，而构型x2则为参数平面上的uv坐标。计算整个网格上对称Dirichlet能量的matlab代码如下：

 function E = computeTotalSDEnergy(F, uv, X1, As)
%% Compute total symmetric Dirichlet energy of the mesh with current parameterization
%% Args:
%%      F[nF, 3]: face connectivity
%%      uv[nV, 2]: vertex coordinates in 2D
%%      X1[3, 2, nF]: rest pose of each triangle
%%      As[nF, 1]: area of each triangle
%% Returns:
%%      E: total energy on the mesh

nF = size(F, 1);
E  = 0;

for i=1:nF
    %% F[i, :] at local frame
    x1 = X1(:,:,i);
    q  = As(i);
    
    %% F[i, :] at uv plane
    x2= uv(F(i, :), :);

    Eq = computeSDEnergy(x1, x2);
    E = E + q*Eq;
end

end
 

Projected Newton Method

变形体会通过最小化自身的变形能来恢复原始形状。从数学上来看这等价于一个最优化过程，因此我们可以通过求解一个最优化问题来计算变形。在几何处理中我们一般会使用二阶优化方法，和基于梯度的一阶方法相比二阶方法收敛更快而且无需手动设置优化步长。二阶方法中最经典的是牛顿法(Newton’s method)，也称为Newton-Raphson方法(Newton-Raphson method)。我们对目标函数 $f (x)$ 进行二阶近似有：

f (x + d) \approx f (x) + J^{T} (x) d + \frac{1}{2} d^{T} H (x) d

对上式右端按照 $d$ 进行求导并令导数为0可以得到：

H (x) d = - J (x)

当函数 $f (x)$ 的Hessian矩阵非奇异时通过求解线性方程即可得到优化方向：

d = - H^{- 1} (x) J (x)

因此对于参数化问题我们只需要计算对称Dirichlet能量 $Ψ_{SD}$ 关于平面uv坐标的梯度与Hessian矩阵，然后利用牛顿法进行迭代即可，具体方法会在后面的小节进行介绍。需要额外说明的是牛顿法的二阶收敛性需要保证目标函数的Hessian矩阵是正定矩阵，这在很多情况下是无法严格满足的。因此我们会把Hessian矩阵投影到一个附近的正定矩阵上来进行处理，此时的优化算法称为投影牛顿法(projected Newton method)。

Vec Operator

在计算导数时可能会遇到向量或是矩阵对矩阵求导的情况，此时得到的导数是一个高阶的张量(tensor)。高阶的张量非常不直观，因此这里介绍一下vec运算来将矩阵向量化。对于2×2的矩阵 $A$ ，vec运算将其转换为一个4维向量：

vec (A) = vec ([\begin{matrix} a & c \\ b & d \end{matrix}]) = [\begin{matrix} a \\ b \\ c \\ d \end{matrix}]

vec运算的一大优势在于它与求导运算是可交换的，对张量求导再进行向量化等价于对向量化的张量进行求导。

vec (\frac{\partial A}{\partial x}) = \frac{\partial vec (A)}{\partial x}

vec (\frac{\partial A}{\partial x}) = \frac{\partial A}{\partial vec (x)}

这样利用vec进行向量化就可以把所有的张量运算转换为矩阵和向量的运算，使得推导和编程都更加清晰直观。

 function [a] = vec(A)
%% vec operator on 2nd-order tensor

[rows, cols] = size(A);
a = reshape(A, rows * cols , 1);

end
 

Algorithm Framework

基于最小化变形能量的参数化算法框架如下：

计算一个无翻转无重叠的参数化作为初始uv坐标，一般可由Tutte参数化来实现；
计算变形能量关于当前uv坐标的梯度；
计算变形能量关于当前uv坐标的Hessian矩阵，并投影到正定矩阵空间中；
求解线性方程 $H d = - J$ ，得到uv坐标的更新方向；
计算保证无翻转的更新步长，更新uv坐标；
返回第2步进行迭代，直到梯度的模长足够小算法收敛。

整个算法框架和每一步的迭代分别封装在projectedNewtonSolver()以及projectedNewtonStep()两个函数中：

 function uv = projectedNewtonSolver(V, F, uv, maxSteps, maxIter, lam, c, tau)
%% One step projected Newton solver
%% Args:
%%      V[nV, 3]: vertex coordinates in 3D
%%      F[nF, 3]: face connectivity
%%      uv[nV, 2]: vertex coordinates in 2D
%%      maxSteps: maximum steps in projected Newton
%%      maxIter: maximum iteration in line search
%%      lam: positive-definite parameter
%%      c: search control parameters
%%      tau: shrinkage factor in line search
%% Returns:
%%      uv_new[nV, 2]: updated uv coordinates

nV = size(V, 1);
nF = size(F, 1);

%% reusable parameters
As    = zeros(nF, 1);       %% triangle area
X1    = zeros(3, 2, nF);    %% triangle at local frame
PfPxs = zeros(4, 6, nF);    %% PfPx at different triangles

for i=1:nF  
    x1 = project2Plane(V(F(i, :), :));

    As(i)        = Area(x1);
    X1(:,:,i)    = x1;
    PfPxs(:,:,i) = computePfPx(x1);
end

%% projected Newton step
for i=1:maxSteps
    fprintf('\nstep %d\n', i);
    
    [uv, stop] = projectedNewtonStep(V, F, uv, X1, As, PfPxs, maxIter, lam, c, tau);
    
    if stop
        break
    end
end

end
 

 function [uv_new, stop] = projectedNewtonStep(V, F, uv, X1, As, PfPxs, maxIter, lam, c, tau)
%% One step projected Newton update
%% Args:
%%      V[nV, 3]: vertex coordinates in 3D
%%      F[nF, 3]: face connectivity
%%      uv[nV, 2]: vertex coordinates in 2D
%%      X1[3, 2, nF]: rest pose of each triangle
%%      As[nF, 1]: triangle areas
%%      PfPxs[4, 6, nF]: PfPx on different triangle
%%      maxIter: maximum iteration in line search
%%      lam: positive-definite parameter
%%      c: search control parameter
%%      tau: shrinkage factor in line search
%% Returns:
%%      uv_new[nV, 2]: updated uv coordinates
%%      stop: stop iteration flag

nV = size(V, 1);
nF = size(F, 1);

stop = false;

%% gradient
b = computeGradient(V, F, uv, X1, As, PfPxs);
g = reshape(b, [nV, 2]);

%% early stop
if norm(b) < 1e-4
    uv_new = uv;
    stop = true;
    return
end

%% Hessian
H = projectHessian(V, F, uv, X1, As, PfPxs);
H = H + lam * speye(nV * 2);

%% solve the update direction
d = -H \ b;
d = reshape(d, [nV, 2]);

%% one step update
[uv_new, E, alpha] = lineSearch(F, uv, X1, As, d, g, maxIter, c, tau);

fprintf("SD Energy: %.2e \tgradient: %.2e \talpha: %.2e\n", E, norm(b), alpha);

end
 

Gradient

直接计算 $Ψ$ 关于网格顶点uv坐标的导数往往是比较困难的，因此我们首先利用链式法则和变形梯度对导数进行分解：

\begin{aligned} \frac{\partial Ψ}{\partial x} & = \sum_{q} | q | \frac{\partial}{\partial x} Ψ_{SD} (F_{q}) \\ = \sum_{q} | q | {\frac{\partial f_{q}}{\partial x}}^{T} \frac{\partial Ψ_{q}}{\partial f_{q}} \end{aligned}

注意这里 $f_{q} = vec (F_{q})$ 表示向量化后的变形梯度。根据上式我们只需要分别计算每个三角形上变形梯度关于uv坐标的导数 $\frac{\partial f}{\partial x}$ 以及对称Dirichlet能量关于变形梯度的导数 $\frac{\partial Ψ}{\partial f}$ 两部分，再把它们乘起来按照初始构型的面积加权就得到了整体的梯度。

Invariants

对称Dirichlet能量 $Ψ_{SD}$ 由于存在变形梯度的逆使得直接计算导数会非常的困难，而pixar这篇文章的主要贡献就在于给出了伸缩张量 $S$ 的不变量(invariants)所定义的变形能等价形式，以及能量关于不变量导数和Hessian矩阵的解析式。这里我们首先来介绍不变量的概念，然后利用不变量来重新表示变形能量。

文章中作者定义伸缩张量 $S$ 的3个不变量如下：

I_{1} = tr (S) = \sum_{i} σ_{i}

I_{2} = ‖ S ‖_{F}^{2} = tr (S^{T} S) = \sum_{i} σ_{i}^{2}

I_{3} = det (S) = Π_{i} σ_{i}

需要说明的是二维情况下这3个不变量并不是独立的，实际上 $I_{2}$ 可以由 $I_{1}$ 和 $I_{3}$ 推导出来：

I_{2} = I_{1}^{2} - 2 I_{3}

同时，二维情况下上面定义的不变量也满足变形梯度 $F$ 的特征多项式(characteristic polynomial)：

σ_{i}^{2} - I_{1} σ_{i} + I_{3} = 0

使用不变量来重新表示对称Dirichlet能量如下：

\begin{aligned} Ψ_{SD} & = \frac{1}{2} (σ_{1}^{2} + σ_{2}^{2} + \frac{1}{σ_{1}^{2}} + \frac{1}{σ_{2}^{2}}) \\ = \frac{1}{2} (I_{2} + \frac{I_{2}}{I_{3}^{2}}) \end{aligned}

Compute $\frac{\partial Ψ}{\partial f}$

计算 $Ψ_{SD}$ 关于变形梯度 $f$ 的导数时我们再次利用链式法则将它拆成关于不变量的两部分：

\frac{\partial Ψ_{SD}}{\partial f} = \sum_{i} \frac{\partial Ψ_{SD}}{\partial I_{i}} \frac{\partial I_{i}}{\partial f}

$Ψ_{SD}$ 关于不变量的导数比较容易，直接求导即可：

\frac{\partial Ψ_{SD}}{\partial I_{1}} = 0, \frac{\partial Ψ_{SD}}{\partial I_{2}} = \frac{1}{2} (1 + \frac{1}{I_{3}^{2}}), \frac{\partial Ψ_{SD}}{\partial I_{3}} = - \frac{I_{2}}{I_{3}^{3}}

而不变量关于变形梯度的导数则比较复杂，这里直接给出结论具体推导可以参考原文：

\frac{\partial I_{1}}{\partial f} = vec (R) = r

\frac{\partial I_{2}}{\partial f} = 2 f

\frac{\partial I_{3}}{\partial f} = vec (G) = g, G = U [\begin{matrix} σ_{2} & 0 \\ 0 & σ_{1} \end{matrix}] V^{T}

在计算 $G$ 矩阵时实际上无需对变形梯度 $F$ 进行分解，我们可以使用一些矩阵变换技巧来直接得到：

G = T F T^{T}, T = [\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}]

整理之后就得到了计算能量函数关于变形梯度导数的函数computePphiPf()：

 function PphiPf = computePphiPf(F)
%% Compute partial derivative PphiPf of symmetric Dirichlet energy
%% Args:
%%      F[2, 2]: deformation graident
%% Returns:
%%      PphiPf[4, 1]: partial derivative

%% polar decomposition
[R, S] = polarDecomposition(F);

%% invariants
I1 = trace(S);
I2 = norm(S, "fro")^2;
I3 = det(S);

%% PPhiPI
PphiPI = [0; 0.5*(1+1/I3^2); -I2/I3^3];

%% PIPf
PIPf = zeros(4, 3);
PIPf(:, 1) = vec(R);
PIPf(:, 2) = 2*vec(F);

twist = [0, -1; 
         1,  0];
G = twist * F * twist';
PIPf(:, 3) = vec(G);

%% PphiPf
PphiPf = PIPf * PphiPI;

end
 

Compute $\frac{\partial f}{\partial x}$

整个梯度项的第2部分需要考虑变形梯度 $f$ 关于uv坐标的导数 $\frac{\partial f}{\partial x}$ 。首先回忆变形梯度的计算过程：

F = D_{s} D_{m}^{- 1}

其中 $D_{m}^{- 1}$ 项只与初始构型有关，而与网格的当前uv坐标无关。因此对等式两边同时求导得到：

\frac{\partial F}{\partial x} = \frac{\partial D_{s}}{\partial x} D_{m}^{- 1}

上式说明我们只需要考虑 $D_{s}$ 项的导数再乘以 $D_{m}^{- 1}$ 即可。回忆 $D_{s}$ 项的定义：

\begin{aligned} D_{s} & = [\begin{array}{c} o_{1} | o_{2} \end{array}] \\ = [\begin{array}{c} x_{1} - x_{0} | x_{2} - x_{0} \end{array}] \\ = [\begin{array}{c} u_{x_{1}} - u_{x_{0}} & u_{x_{2}} - u_{x_{0}} \\ v_{x_{1}} - v_{x_{0}} & v_{x_{2}} - v_{x_{0}} \end{array}] \end{aligned}

对其进行求导展开后得到：

\frac{\partial D_{s}}{\partial u_{x_{0}}} = [\begin{matrix} - 1 & - 1 \\ 0 & 0 \end{matrix}], \frac{\partial D_{s}}{\partial v_{x_{0}}} = [\begin{matrix} 0 & 0 \\ - 1 & - 1 \end{matrix}]

\frac{\partial D_{s}}{\partial u_{x_{1}}} = [\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix}], \frac{\partial D_{s}}{\partial v_{x_{1}}} = [\begin{matrix} 0 & 0 \\ 1 & 0 \end{matrix}]

\frac{\partial D_{s}}{\partial u_{x_{2}}} = [\begin{matrix} 0 & 1 \\ 0 & 0 \end{matrix}], \frac{\partial D_{s}}{\partial v_{x_{2}}} = [\begin{matrix} 0 & 1 \\ 0 & 0 \end{matrix}]

再带入 $D_{m}^{- 1}$ 可以得到完整的导数：

D_{m}^{- 1} = [\begin{matrix} a & c \\ b & d \end{matrix}]

\frac{\partial F}{\partial u_{x_{0}}} = [\begin{matrix} - (a + b) & - (c + d) \\ 0 & 0 \end{matrix}], \frac{\partial F}{\partial v_{x_{0}}} = [\begin{matrix} 0 & 0 \\ - (a + b) & - (c + d) \end{matrix}]

\frac{\partial F}{\partial u_{x_{1}}} = [\begin{matrix} a & c \\ 0 & 0 \end{matrix}], \frac{\partial F}{\partial v_{x_{1}}} = [\begin{matrix} 0 & 0 \\ a & c \end{matrix}]

\frac{\partial F}{\partial u_{x_{2}}} = [\begin{matrix} b & d \\ 0 & 0 \end{matrix}], \frac{\partial F}{\partial v_{x_{2}}} = [\begin{matrix} 0 & 0 \\ b & d \end{matrix}]

最后利用向量化的技术就得到了矩阵形式的 $\frac{\partial f}{\partial x}$ ：

\begin{aligned} \frac{\partial f}{\partial x} & = vec (\frac{\partial F}{\partial x}) \\ = [\begin{array}{c} vec (\frac{\partial F}{\partial u_{x_{0}}}) | vec (\frac{\partial F}{\partial v_{x_{0}}}) | vec (\frac{\partial F}{\partial u_{x_{1}}}) | vec (\frac{\partial F}{\partial v_{x_{1}}}) | vec (\frac{\partial F}{\partial u_{x_{2}}}) | vec (\frac{\partial F}{\partial v_{x_{2}}}) \end{array}] \\ = [\begin{array}{c} - (a + b) & 0 & a & 0 & b & 0 \\ 0 & - (a + b) & 0 & a & 0 & b \\ - (c + d) & 0 & c & 0 & d & 0 \\ 0 & - (c + d) & 0 & c & 0 & d \end{array}] \end{aligned}

 function PfPx = computePfPx(x1)
%% Compute partial derivative of deformation gradient
%% Args:
%%      x1: rest pose
%% Returns:
%%      PfPx: partial derivative of deformation gradient

Dm = (x1(2:end, :) - x1(1, :))';
DmInv = matrixInv2x2(Dm);
a = DmInv(1, 1); c = DmInv(1, 2);
b = DmInv(2, 1); d = DmInv(2, 2);

PfPx = zeros(4, 6);

%% hard coded PfPx
PfPx(1, 1) = -(a+b); PfPx(1, 2) =      0; PfPx(1, 3) = a; PfPx(1, 4) = 0; PfPx(1, 5) = b; PfPx(1, 6) = 0;
PfPx(2, 1) =      0; PfPx(2, 2) = -(a+b); PfPx(2, 3) = 0; PfPx(2, 4) = a; PfPx(2, 5) = 0; PfPx(2, 6) = b;
PfPx(3, 1) = -(c+d); PfPx(3, 2) =      0; PfPx(3, 3) = c; PfPx(3, 4) = 0; PfPx(3, 5) = d; PfPx(3, 6) = 0;
PfPx(4, 1) =      0; PfPx(4, 2) = -(c+d); PfPx(4, 3) = 0; PfPx(4, 4) = c; PfPx(4, 5) = 0; PfPx(4, 6) = d;

end
 

Compute Gradient

把上面的代码组合起来就得到了每个小三角形上计算能量函数关于uv坐标梯度项的代码，需要注意的是这里梯度项uv坐标的排布顺序为：

\frac{\partial Ψ_{SD}}{\partial x} |_{q} = {[\begin{matrix} \frac{\partial Ψ_{SD}}{\partial u_{x_{0}}} & \frac{\partial Ψ_{SD}}{\partial v_{x_{0}}} & \frac{\partial Ψ_{SD}}{\partial u_{x_{1}}} & \frac{\partial Ψ_{SD}}{\partial v_{x_{1}}} & \frac{\partial Ψ_{SD}}{\partial u_{x_{2}}} & \frac{\partial Ψ_{SD}}{\partial v_{x_{2}}} \end{matrix}]}^{T}

而在整个三角网格上我们只需要把每个三角形每个顶点上的梯度按照相同编号进行累加即可，具体的代码如下：

 function b = computeGradient(V, F, uv, X1, As, PfPxs)
%% Gradient of symmetric Dirichlet energy
%% Args:
%%      V[nV, 3]: vertex coordinates in 3D
%%      F[nF, 3]: face connectivity
%%      uv[nV, 2]: vertex coordinates in 2D
%%      X1[3, 2, nF]: rest pose of each triangle
%%      As[nF, 1]: triangle areas
%%      PfPxs[4, 6, nF]: PfPx on different triangle
%% Returns:
%%      b[nV*2, 1]: gradient of current parameterization

nV = size(V, 1);
nF = size(F, 1);

b = zeros(nV, 2);

%% loop over all the triangles
for i=1:nF
    %% F[i, :] at local frame
    x1 = X1(:,:,i);
    q  = As(i);

    %% F[i, :] at uv plane
    x2= uv(F(i, :), :);
    
    %% deformation gradient from local frame to uv plane
    f = findJacobian(x1, x2);
    
    %% partial derivative of phi
    PphiPf = computePphiPf(f);
    PfPx   = PfPxs(:,:,i);

    PphiPx = q * PfPx' * PphiPf;

    v1 = F(i, 1); v2 = F(i, 2); v3 = F(i, 3);
    b(v1, 1) = b(v1, 1) + PphiPx(1); b(v1, 2) = b(v1, 2) + PphiPx(2);
    b(v2, 1) = b(v2, 1) + PphiPx(3); b(v2, 2) = b(v2, 2) + PphiPx(4);
    b(v3, 1) = b(v3, 1) + PphiPx(5); b(v3, 2) = b(v3, 2) + PphiPx(6);
end

b = reshape(b, [nV*2, 1]);

end
 

Hessian

Hessian矩阵的计算与梯度项是非常类似的，我们利用链式法则对Hessian矩阵进行展开：

\begin{aligned} \frac{\partial^{2} Ψ}{\partial x^{2}} & = \sum_{q} | q | \frac{\partial^{2}}{\partial x^{2}} Ψ_{SD} (F_{q}) \\ = \sum_{q} | q | {\frac{\partial f_{q}}{\partial x}}^{T} (\frac{\partial^{2} Ψ_{q}}{\partial f_{q}^{2}}) \frac{\partial f_{q}}{\partial x} \end{aligned}

其中 $\frac{\partial f_{q}}{\partial x}$ 已经在Compute $\frac{\partial f}{\partial x}$ 中推导过了，因此我们只需要关注能量函数关于变形梯度的二阶项 $\frac{\partial^{2} Ψ_{q}}{\partial f_{q}^{2}}$ 。

Analytic Eigen System

再次利用链式法则，把二阶项 $\frac{\partial^{2} Ψ}{\partial f^{2}}$ 展开为关于不变量的形式：

\begin{aligned} \frac{\partial^{2} Ψ}{\partial f^{2}} & = \sum_{i} \frac{\partial Ψ}{\partial I_{i}} \frac{\partial^{2} I_{i}}{\partial f^{2}} \\ + \sum_{i} \frac{\partial^{2} Ψ}{\partial I_{i}^{2}} (\frac{\partial I_{i}}{\partial f}) (\frac{\partial I_{i}}{\partial f})^{T} \\ + \sum_{i < j} \frac{\partial^{2} Ψ}{\partial I_{i} \partial I_{j}} [(\frac{\partial I_{i}}{\partial f}) (\frac{\partial I_{j}}{\partial f})^{T} + (\frac{\partial I_{j}}{\partial f}) (\frac{\partial I_{i}}{\partial f})^{T}] \end{aligned}

式中红色的部分为不变量关于变形梯度的二阶导数，而蓝色的部分为不变量关于变形梯度的一阶导数。其中关于变形梯度的一阶导数我们已经在Compute $\frac{\partial Ψ}{\partial f}$ 中推导过，而二阶导数同样具有解析的形式，具体可以参考原文。

实际上我们并不需要真的去计算各个导数来构造Hessian矩阵，利用特征值分解可以把Hessian矩阵表示为：

H = E Λ E^{- 1} = \sum_{i} λ_{i} e_{i} e_{i}^{T}

式中 $λ_{i}$ 为Hessian矩阵的特征值，而 $e_{i}$ 则是对应的特征向量。上式说明Hessian矩阵可以表示为特征值与特征向量的组合，而对于 $Ψ_{SD}$ 其4个特征值与特征向量都具有解析形式：

λ_{1} = 1 + \frac{3}{σ_{1}^{4}}, e_{1} = d_{1}

λ_{2} = 1 + \frac{3}{σ_{2}^{4}}, e_{2} = d_{2}

λ_{3} = 1 + \frac{1}{I_{3}^{2}} + \frac{I_{2}}{I_{3}^{3}}, e_{3} = l

λ_{4} = 1 + \frac{1}{I_{3}^{2}} - \frac{I_{2}}{I_{3}^{3}}, e_{4} = t

d_{1} = vec (D_{1}) = vec (U [\begin{matrix} 1 & 0 \\ 0 & 0 \end{matrix}] V^{T})

d_{2} = vec (D_{2}) = vec (U [\begin{matrix} 0 & 0 \\ 0 & 1 \end{matrix}] V^{T})

l = vec (L) = vec (\frac{1}{\sqrt{2}} U [\begin{matrix} 0 & 1 \\ 1 & 0 \end{matrix}] V^{T})

t = vec (T) = vec (\frac{1}{\sqrt{2}} U [\begin{matrix} 0 & - 1 \\ 1 & 0 \end{matrix}] V^{T})

这些解析形式特征值与特征向量的推导都非常的复杂，具体过程需要参考原文以及补充材料。不过在编程时我们只需要对变形梯度 $F$ 进行SVD和极分解就可以显式地构造出特征系统，见SDEigenSystem()函数：

 function [eigs, lams] = SDEigenSystem(F)
%% Analytic eigensystem of symmetric Dirichlet energy
%% Args:
%%      F[2, 2]: deformation gradient
%% Returns:
%%      eigs[4, 1]: eigen values
%%      lams[4, 4]: each column is an eigen vector

%% SVD
[R, S, U, Sigma, V] = polarSVD(F);

I1 = trace(S);
I2 = norm(S, "fro")^2;
I3 = det(S);

%% eigen values
eigs = zeros(4, 1);
eigs(1) = 1 + 3/Sigma(1,1)^4;
eigs(2) = 1 + 3/Sigma(2,2)^4;
eigs(3) = 1 + 1/I3^2 + I2/I3^3;
eigs(4) = 1 + 1/I3^2 - I2/I3^3;

%% eigen vectors
lams = zeros(4);

D1 = U * [1, 0; 0, 0] * V';
D2 = U * [0, 0; 0, 1] * V';
L  = U * [0, 1; 1, 0] * V' / sqrt(2);
T  = U * [0,-1; 1, 0] * V' / sqrt(2);

lams(:, 1) = vec(D1);
lams(:, 2) = vec(D2);
lams(:, 3) = vec(L);
lams(:, 4) = vec(T);

end
 

Compute Hessian

利用解析形式的特征值和特征向量可以构造每个三角形上6×6的Hessian矩阵，这里需要注意矩阵元素的编号顺序：

\frac{\partial^{2} Ψ_{SD}}{\partial x^{2}} |_{q} = [\begin{matrix} \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}}^{2}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial v_{x_{0}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial u_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial v_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial u_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial v_{x_{2}}} \\ \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial v_{x_{0}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{0}}^{2}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{1}} \partial v_{x_{0}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{1}} \partial v_{x_{0}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{2}} \partial v_{x_{0}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{2}} \partial v_{x_{0}}} \\ \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial u_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{0}} \partial u_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{1}}^{2}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{1}} \partial u_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{2}} \partial u_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{2}} \partial u_{x_{1}}} \\ \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial v_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{0}} \partial v_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{1}} \partial v_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{1}}^{2}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{2}} \partial v_{x_{1}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{2}} \partial v_{x_{1}}} \\ \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial u_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{0}} \partial u_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{1}} \partial u_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{1}} \partial u_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{2}}^{2}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{2}} \partial u_{x_{2}}} \\ \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{0}} \partial v_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{0}} \partial v_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{1}} \partial v_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{1}} \partial v_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial u_{x_{2}} \partial v_{x_{2}}} & \frac{\partial^{2} Ψ_{SD}}{\partial v_{x_{2}}^{2}} \end{matrix}]

最后把每个三角形顶点上的Hessian矩阵按照顶点编号进行累加就得到了整个网格上的Hessian矩阵，当然为了保证Hessian矩阵的正定性我们需要手动令每个三角形上的所有特征值都至少为0。具体代码可参见projectHessian()函数：

 function H = projectHessian(V, F, uv, X1, As, PfPxs)
%% Projected Hessian of symmetric Dirichlet energy
%% Args:
%%      V[nV, 3]: vertex coordinates in 3D
%%      F[nF, 3]: face connectivity
%%      uv[nV, 2]: vertex coordinates in 2D
%%      X1[3, 2, nF]: rest pose of each triangle
%%      As[nF, 1]: triangle areas
%%      PfPxs[4, 6, nF]: PfPx on different triangle
%% Returns:
%%      H[nV*2, nV*2]: Hessian matrix of current parameterization

nV = size(V, 1);
nF = size(F, 1);

%% sparse matrix indices and values
I = zeros(nF*36, 1); J = zeros(nF*36, 1);
Value = zeros(nF*36, 1);

%% update with each triangle
for i=1:nF
    %% F[i, :] at local frame
    x1 = X1(:,:,i);
    q  = As(i);

    %% F[i, :] at uv plane
    x2= uv(F(i, :), :);
    
    %% deformation gradient from local frame to uv plane
    f = findJacobian(x1, x2);

    %% SVD
    [eigs, lams] = SDEigenSystem(f);
    eigs = max(eigs, 0);

    %% local Hessian
    Hq = lams * diag(eigs) * lams';

    %% PfPx
    PfPx = PfPxs(:,:,i);
    Hq   = q * PfPx' * Hq * PfPx;

    %% update indices and values
    v1x = F(i, 1); v2x = F(i, 2); v3x = F(i, 3);
    v1y = v1x+nV;  v2y = v2x+nV;  v3y = v3x+nV;

    indices = [v1x, v1y, v2x, v2y, v3x, v3y];
    idxI = zeros(6) + indices';
    idxJ = zeros(6) + indices;
    
    idx = (i-1)*36+1;   % starting index
    I(idx:idx+35)     = reshape(idxI, 1, []);
    J(idx:idx+35)     = reshape(idxJ, 1, []);
    Value(idx:idx+35) = reshape(Hq,   1, []);
end

H = sparse(I, J, Value, nV*2, nV*2);

end
 

Line Search

有了梯度项和Hessian矩阵后通过求解矩阵方程就得到了网格顶点uv坐标的更新方向 $d$ 。在标准牛顿法中我们可以直接使用 $d$ 来进行更新，但对于参数化问题直接使用 $d$ 进行更新可能会跳过 $Ψ_{SD}$ 的障碍造成三角形的翻转。为了避免这种情况我们需要使用线搜索(line search)算法来计算合适的更新步长 $α$ 。

Flip Avoid

在每个三角形上对于给定的uv坐标和更新方向 $d$ ，出现翻转的条件为顶点相邻边构成的行列式小于0，即：

det ([\begin{matrix} u_{x_{1}} - u_{x_{0}} & u_{x_{2}} - u_{x_{0}} \\ v_{x_{1}} - v_{x_{0}} & v_{x_{2}} - v_{x_{0}} \end{matrix}] + α [\begin{matrix} d_{u_{x_{1}}} - d_{u_{x_{0}}} & d_{u_{x_{2}}} - d_{u_{x_{0}}} \\ d_{v_{x_{1}}} - d_{v_{x_{0}}} & d_{v_{x_{2}}} - d_{v_{x_{0}}} \end{matrix}]) < 0

整理后可以得到关于步长 $α$ 的二次不等式：

A α^{2} + B α + C < 0

其中各个系数满足：

A = det (D), B = det (D_{s}) + det (D) - det (D_{s} - D), C = det (D_{s})

D_{s} = [\begin{matrix} u_{x_{1}} - u_{x_{0}} & u_{x_{2}} - u_{x_{0}} \\ v_{x_{1}} - v_{x_{0}} & v_{x_{2}} - v_{x_{0}} \end{matrix}]

D = [\begin{matrix} d_{u_{x_{1}}} - d_{u_{x_{0}}} & d_{u_{x_{2}}} - d_{u_{x_{0}}} \\ d_{v_{x_{1}}} - d_{v_{x_{0}}} & d_{v_{x_{2}}} - d_{v_{x_{0}}} \end{matrix}]

这3个系数的解析形式非常不直观，但可以使用符号计算来验证其正确性。

得到系数后可以求解二次方程 $A α^{2} + B α + C = 0$ ，其中较小的那个根即为当前三角形上满足无翻转的更新步长上界。通过对网格上所有三角形进行遍历，所有三角形上步长上界中的最小值即为整个网格更新步长的上界。代码可参见initStepSize()和flipCheckStepSize()两个函数。

 function alpha = initStepSize(F, uv, d)
%% A helper function to find initial step size in line search
%% Args:
%%      F[nF, 3]: face connectivity
%%      uv[nV, 2]: vertex coordinates in 2D
%%      d[nV, 2]: update direction
%% Returns:
%%      alpha: maximum allowed step size

alpha = inf;
nF = size(F, 1);

for i=1:nF
    x = uv(F(i, :), :);
    dx= d(F(i, :), :);
    
    alpha = min([alpha, flipCheckStepSize(x, dx)]);
end

end
 

 function alpha = flipCheckStepSize(x, dx)
%% A helper function to find maximum step size to prevent flip
%% Args:
%%      x[3, 2]: triangle coordinates in 2D, each row is a coordinate vector
%%      dx[3, 2]: update direction
%% Return:
%%      alpha: maximum step size

alpha = inf;

Ds= x(2:end, :) - x(1, :);
D = dx(2:end, :) - dx(1, :);

%% quadratic equation
detDs= det(Ds);
detD = det(D);

A = detD;
B = detDs + detD - det(Ds-D);
C = detDs;

delta = B*B - 4*A*C;

if delta < 0
    return
else
    alpha = (-B-sqrt(delta))/(2*A);
    if alpha < 0
        alpha = inf;
    end
end

end
 

Armijo–Goldstein Condition

在计算更新步长时除了需要考虑三角形翻转外还可以结合Armijo–Goldstein准则(Armijo–Goldstein condition)来估计最优步长。Armijo–Goldstein准则有两个参数 $c$ 和 $τ$ ，其中 $c$ 用来衡量目标函数下降的程度而 $τ$ 则用来控制步长 $α$ 的缩减。我们认为当目标函数 $f$ 在 $x$ 位置上沿 $d$ 方向上的步长满足如下条件时就完成了充分地下降：

f (x) - f (x + α d) \geq c \cdot α (- \nabla f (x) \cdot d)

否则需要使用参数 $τ$ 来缩减步长 $α$ ：

α \leftarrow τ \cdot α

这样我们就只需要从满足网格无翻转条件的最大步长 $α_{max}$ 开始不断迭代直到满足Armijo–Goldstein准则就得到了最优更新步长，这一过程在lineSearch()函数中实现。

 function [uv_new, Enew, alpha] = lineSearch(F, uv, X1, As, d, g, maxIter, c, tau)
%% Line search to update parameterization
%% Args:
%%      F[nF, 3]: face connectivity
%%      uv[nV, 2]: vertex coordinates in 2D
%%      X1[3, 2, nF]: rest pose of each triangle
%%      As[nF, 1]: triangle areas
%%      d[nV, 2]: update direction
%%      g[nV, 2]: gradients
%%      maxIter: maximum iteration
%%      c: search control parameter
%%      tau: shrinkage factor in line search
%% Returns:
%%      uv_new[nV, 2]: updated uv coordiantes
%%      Enew: updated energy
%%      alpha: optimal step size

%% initialize step size with flip check
alpha = initStepSize(F, uv, d);
alpha = min(0.99*alpha, 1.0);

%% backtracking line search
Eold = computeTotalSDEnergy(F, uv, X1, As);
Enew = Eold;

uv_new = uv;

for i=1:maxIter
    uv_new = uv + alpha * d;
    Enew = computeTotalSDEnergy(F, uv_new, X1, As);

    %% Armijo–Goldstein condition
    if Enew <= Eold + c * alpha * sum(d .* g, "all")
        return;
    end

    alpha = alpha * tau;
end

end
 