libigl第四章-变形

现代基于网格的形状变形方法满足手柄（网格上选定的顶点或区域）处的用户变形约束，并将这些手柄变形平滑地传播到形状的其余部分，而不会删除或扭曲细节。Libigl 提供了各种最先进的变形技术的实现，从基于二次网格的能量最小化器到蒙皮方法，再到非线性弹性启发技术。

双谐波变形

2000 年至 2010 年期间的研究产生了一系列技术，这些技术将基于手柄的形状变形问题视为二次能量最小化问题或等效于线性偏微分方程的解。这些技术有很多种，但一个典型的子集是那些考虑双拉普拉斯方程的解，即双调和函数。这个四阶 PDE 在边界条件中提供了足够的灵活性，以确保处理手柄约束的${\mathbf{C}}^1$连续性（在细化的限制中）。

双谐波表面

让我们首先通过考虑双谐波表面来开始讨论双谐波变形。我们将大致定义双调和曲面为其位置函数相对于一些初始参数化是双调和的曲面：
${\Delta }^2{x}^{\prime }=0$
并受到一些手柄约束，概念化为“边界条件”：
${\mathbf{x}}_b^{\prime }={\mathbf{x}}_{bc}$
其中${\mathbf{x}}^{\prime }$是表面上一个点的未知 3D 位置。所以我们要求每个空间坐标函数的双拉普拉斯算子为零。
在 libigl 中，可以使用 igl::harmonic 并设置 K=2（双谐波）来解决双谐波问题：

// U_bc contains deformation of boundary vertices b
igl::harmonic(V,F,b,U_bc,2,U);

这会生成一个插入手柄约束的平滑曲面，但曲面上的所有原始细节都将被平滑掉。最明显的是，如果原始曲面还不是双调和的，那么为所有手柄提供恒等变形（将它们保持在静止位置）将不会重现原始曲面。相反，结果将是对那些手柄位置进行插值的双调和曲面。
因此，我们可以得出结论，这不是一种直观的形状变形技术。

双谐波变形场

现在我们知道变形技术的一个有用属性是“静止姿势再现”：对手柄应用不变形应该不会对形状应用变形。
为了通过构造来保证这一点，我们可以使用变形场（即位移）$\mathbf{d}$ 而不是直接使用位置$\mathbf{x}$。
${\mathbf{x}}^{\prime }=\mathbf{x}+\mathbf{d}$.
插入手柄约束的变形场的平滑变形场 $\mathbf{d}$ 将施加平滑变形形状${\mathbf{x}}^{\prime }$。 自然地，我们考虑双谐波变形场：
${\Delta }^2\mathbf{d}=0$.
受到相同的手柄约束，但根据边界处的隐含变形场（手柄）重写：
${\mathbf{d}}_b={\mathbf{x}}_{bc}-{\mathbf{x}}_b$.
我们可以再次使用 igl::harmonic并设置k=2，但这次先求解变形场，然后恢复变形位置：

// U_bc contains deformation of boundary vertices b
D_bc = U_bc - igl::slice(V,b,1);
igl::harmonic(V,F,b,D_bc,2,D);
U = V+D;

双调和函数（无论是位置还是位移）是双拉普拉斯方程的解，也是“拉普拉斯能量”的极小值。对于位移$\mathbf{d}$而言，能量可以表示为
${\int }_S||\Delta \mathbf{d}|{|}^{2}d\mathbf{A}$,
其中，我们将$\Delta \mathbf{d}$定义为简单地应用拉普拉斯坐标。根据拉普拉斯（-Beltrami）算子的线性性，我们可以用原始位置$\mathbf{x}$和未知位置${\mathbf{x}}^{\prime }=\mathbf{x}-\mathbf{d}$重新表达这种能量:
${\int }_S||\Delta ({\mathbf{x}}^{\prime }-\mathbf{x})|{|}^{2}d\mathbf{A}={\int }_S||\Delta {\mathbf{x}}^{\prime }-\Delta x|{|}^{2}d\mathbf{A}$
在 Sorkine 等人的早期工作中，$\Delta {\mathbf{x}}^{\prime }$和$\Delta \mathbf{x}$被称为“微分坐标”。因此，它们的变形（没有线性化旋转）相当于双谐波变形场。

多谐变形

我们可以通过考虑拉普拉斯算子的不同幂来概括双谐波变形，从而得到以下形式的一系列偏微分方程：
${\Delta }^k\mathbf{d}=0$.
其中$k\in 1,2,3...$.k 的选择决定了手柄的连续性水平。尤其是，$k=1$意味着边界处的$C^0$连续，$k=2$意味着边界处的$C^1$连续，$k=3$意味着边界处的$C^2$连续，$k$意味着边界处的$C^{k-1}$连续。

int k = 2;// or 1,3,4,...
igl::harmonic(V,F,b,bc,k,Z);

有界双谐波权重

在计算机动画中，形状变形通常被称为“蒙皮”。约束表现为角色内部刚性“骨骼”的相对旋转。变形方法或蒙皮方法确定角色的表面（即其皮肤）应如何作为骨骼旋转的函数移动。
最流行的技术是线性混合蒙皮。 形状上的每个点都将其新位置计算为骨骼变换的线性组合：
${\mathbf{x}}^{\prime }={\sum }_{i=2}^m{\omega }_i(\mathbf{x}){\mathbf{T}}_i\left(\begin{array}{c}{\mathbf{x}}_i\\ 1\end{array}\right)$.
其中${\omega }_i$为在$\mathbf{x}$处评估第 i 块骨骼的标量权函数，而${\mathbf{T}}_i$是一个表示骨骼变换的$4\times 3$的矩阵。
这个公式是令人尴尬的(why embarassingly?)并行（一点的计算不依赖于另一点计算所需的共享数据）。它通常被实现为顶点着色器。每个顶点的权重和静止位置作为顶点着色器属性发送，骨骼变换作为统一发送。然后顶点在顶点着色器中进行变换，正好赶上渲染.
由于蒙皮公式是线性的（因此得名），我们可以将其写为矩阵乘法：
${\mathbf{X}}^{\prime }=\mathbf{M}\mathbf{T}$,
其中$\mathbf{X}‘$作为行向量变形位置的$n\times 3$的堆栈，$\mathbf{M}$是表示权重以及其他位置信息的$n\times m$的矩阵，$\mathbf{T}$则是变形骨骼的$m\cdot (dim+1)\times dim$的变换堆栈。
传统上，权重函数${\omega }_i$由熟练的索具专业人员手动绘制。现在存在现代技术来自动计算权重函数，给定骨架的形状和描述（或通常任何手柄结构，如笼子、点集合、选定区域等）。
有界双调和权重是一种将权重计算作为约束优化问题的技术。权重通过最小化熟悉的拉普拉斯能量来增强平滑度：
${\sum }_{i=1}^m{\int }_S(\Delta {\omega }_i{)}^{2}dA$
受到强制插值句柄约束的约束：
${\omega }_i(\mathbf{x})=\{\begin{array}{ll}1& if\mathbf{x}\in H_i\\ 0& otherwise\end{array}$,
其中，$H_i$是第i个手柄，以及强制非负性、统一划分和鼓励稀疏性的约束：
$0\le {\omega }_i\le 1and{\sum }_{i=1}^m{\omega }_i=1$.
这是一个二次规划问题，libigl 使用其活动集求解器或调用 Mosek 来解决它。

双四元数蒙皮

即使使用高质量的权重，线性混合蒙皮也是有限的。特别是，它存在由将旋转混合为矩阵而产生的已知伪影：旋转矩阵的权重组合不一定是旋转。考虑一个绕着z轴旋转$-\pi /2$以及$\pi /2$的一个等效的旋转。直觉上，人们可能期望得到单位矩阵，但混合是一个退化矩阵，将 x 和 y 坐标缩放为零：
$0.5\left(\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 1\end{array}\right)+0.5\left(\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 1\end{array}\right)=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right)$
在实践中，这意味着形状会在骨骼权重重叠的区域收缩和塌陷：关节附近。
双四元数蒙皮提供了一种解决方案。该方法将刚性变换表示为一对单位四元数，$\hat{\mathbf{q}}$。线性混合蒙皮公式被双四元数的线性混合取代：
${\mathbf{x}}^{\prime }=\frac{{\sum }_{i=1}^m{\omega }_i(\mathbf{x}){\hat{\mathbf{q}}}_i}{||{\sum }_{i=1}^m{\omega }_i(\mathbf{x}){\hat{\mathbf{q}}}_i||}\mathbf{x}$,
其中${\hat{\mathbf{q}}}_i$是骨骼$i$刚性变换的对偶四元数表示。归一化迫使线性混合的结果再次成为单位对偶四元数，因此也是刚性变换。
与线性混合蒙皮一样，双四元数蒙皮最好在顶点着色器中执行。唯一的区别是骨骼变换是作为对偶四元数而不是仿射变换矩阵发送的。Libigl 使用 igl::dqs 函数支持 CPU 端双四元数蒙皮，该函数采用更传统的刚性变换表示作为输入，并在混合之前在内部转换为双四元数表示：

// vQ is a list of rotations as quaternions
// vT is a list of translations
igl::dqs(V,W,vQ,vT,U);

尽可能刚性

蒙皮和其他线性变形方法本质上是有限的。 尤其是当手柄约束施加大的旋转时会出现困难。
在能量最小化方法的背景下，问题源于比较未变形形状坐标系中的位置（我们的位移）。这些二次能量充其量对于整个形状的全局旋转是不变的，但不是平滑变化的局部旋转。 因此线性技术不会产生非平凡的弯曲和扭曲。
此外，在考虑实体形状（例如用四面体网格离散化）时，线性方法难以保持局部体积，并且它们经常遭受收缩和凸出伪影的困扰。
非线性变形技术为这些问题提供了解决方案。 它们通过将网格顶点的变形与其旋转到最匹配变形的新坐标系的静止位置进行比较来工作。非线性源于变形和最佳拟合旋转的相互依赖。这些技术通常被标记为“尽可能刚性”，因为它们会惩罚所有局部变形与旋转的偏差的总和。
为了得到这样的能量，让我们考虑一个简单的每个三角形能量：
E l i n e a r ( X ′ ) = ∑ t ∈ T α t ∑ { i , j } ∈ t ω i j ∣ ∣ ( x i ′ − x j ′ ) − ( x i − x j ) ∣ ∣ 2 E_{linear}(\mathbf{X'})=\sum_{t\in T}\alpha_t\sum_{\{i,j\}\in t}\omega_{ij}||(\mathbf x_i'-\mathbf x_j')-(\mathbf x_i-\mathbf x_j)||^2 Elinear​(X′)=∑t∈T​αt​∑{i,j}∈t​ωij​∣∣(xi′​−xj′​)−(xi​−xj​)∣∣2
其中${\mathbf{X}}^{\prime }$是网格的未知变形顶点位置，$t$是三角形列表$T$中的一个三角，${\alpha }_t$是三角形$t$的面积， { i , j } \{i,j\} {i,j}是三角形$t$的一条边。因此，该能量测量原始网格$({\mathbf{x}}_i-{\mathbf{x}}_j)$中的边向量与未知网格$({\mathbf{x}}_i^{\prime }-{\mathbf{x}}_j^{\prime })$之间的变化范数。
这种能量不是旋转不变的。如果我们将网格旋转 90 度，则不与旋转轴对齐的边矢量的变化会很大，尽管整体变形是完全刚性的。
因此，“尽可能刚性”的解决方案是为每个限制为旋转的三角形$t$附加辅助变量${\mathbf{R}}_t$。然后重写能量，这次将变形的边缘向量与其旋转的静止对应物进行比较：
E a r a p ( X ′ , ( R 1 , . . . , R ∣ T ∣ ) ) = ∑ t ∈ T α t ∑ { i , j } ∈ t ω i j ∣ ∣ ( x i ′ − x j ′ ) − R t ( x i − x j ) ∣ ∣ 2 E_{arap}(\mathbf{X}',(\mathbf{R}_1,...,\mathbf{R}_{|T|}))=\sum_{t\in T}\alpha_t\sum_{\{i,j\}\in t}\omega_{ij}||(\mathbf x_i'-\mathbf x_j')-\mathbf R_t(\mathbf x_i-\mathbf x_j)||^2 Earap​(X′,(R1​,...,R∣T∣​))=∑t∈T​αt​∑{i,j}∈t​ωij​∣∣(xi′​−xj′​)−Rt​(xi​−xj​)∣∣2.
主要顶点位置变量${\mathbf{X}}^{\prime }$和旋转$({\mathbf{R}}_1,...,{\mathbf{R}}_{|T|})$的分离也导致了优化策略。如果旋转$({\mathbf{R}}_1,...,{\mathbf{R}}_{|T|})$保持固定，则剩余变量${\mathbf{X}}^{\prime }$中的能量是二次的，并且可以通过求解（稀疏）全局线性系统来优化。或者，如果${\mathbf{X}}^{\prime }$ 保持固定，则每次旋转都是局部 procrustes 问题的解决方案（通过$3\times 3$ SVD 或极性分解找到）。这两个步骤——局部的和全局的——每一个都会弱化能量，因此我们可以安全地迭代它们直到收敛。
“尽可能刚性”的不同风格取决于域和边集$T$的维度和余维度。Sorkine 和 Alexa 提出的表面处理技术将$T$视为从每个顶点（辐条）发出的一组边。后来，Chao 等人。 推导出“尽可能刚性”网格能量和共旋转弹性之间的关系，将 0 维元素视为边集：2D 中的三角形和 3D 中的四面体。他们还展示了 Sorkine 和 Alexa 的边集如何不是连续能量的离散化，而是为包含入射在顶点（辐条和轮辋）上的元素的所有边的表面提出了边集。 他们表明，这相当于测量弯曲，尽管是以离散化相关的方式。
Libigl，支持这些常见的风格。 选择一个是在预计算阶段之前设置能量类型的问题：

igl::ARAPData arap_data;
arap_data.energy = igl::ARAP_ENERGY_TYPE_SPOKES;
//arap_data.energy = igl::ARAP_ENERGY_TYPE_SPOKES_AND_RIMS;
//arap_data.energy = igl::ARAP_ENERGY_TYPE_ELEMENTS; //triangles or tets
igl::arap_precomputation(V,F,dim,b,arap_data);

就像 igl::min_quad_with_fixed_* 一样，这个预计算阶段仅取决于网格、固定顶点索引 b 和能量参数。 为了解决b中顶点位置的某些约束，我们可以调用：

igl::arap_solve(bc,arap_data,U);

它使用 U 作为初始猜测，然后计算解决方案。
Libigl 的尽可能刚性变形的实现利用了 McAdams 等人的高度优化的奇异值分解代码。 它利用了 SSE 内在函数。
局部刚性的概念将在表面参数化的背景下很快重新讨论。

快速自动蒙皮转换

毫不奇怪，非线性优化比它的线性优化慢。在尽可能刚性优化的情况下，瓶颈通常是为每个边集（即每个三角形、四面体或顶点单元）恢复最佳拟合旋转所需的大量极坐标分解。即使此代码经过优化，主要自由度的数量也与离散化水平相关，尽管变形的低频行为。
这为快速非线性优化带来了两条路线。 首先，是否有必要（甚至有利）找到这么多最适合的旋转？ 其次，我们能否降低自由度以更好地反映所需变形的频率。
反过来，这些优化最终形成了一种方法，该方法优化了由高质量权重（即手动绘制的权重或有界双谐波权重）跨越的线性混合蒙皮变形空间。该空间是所有可能的网格变形的低维子空间，通过以矩阵形式编写线性混合蒙皮来捕获：
${\mathbf{X}}^{\prime }=\mathbf{M}\mathbf{T}$
其中 $n\times 3$ 矩阵 ${\mathbf{X}}^{\prime }$ 中的网格顶点位置被转置“句柄”的 $(3+1)m\times 3$堆栈中的少量自由度的线性组合替换的转变。
在上面的 ARAP 能量中将 $\mathbf{M}\mathbf{T}$ 换成 ${\mathbf{X}}^{\prime }$ 会立即看到全局求解步骤中的性能提升这是因为 $m<<n$。
局部步骤（拟合旋转）的复杂性仍然与原始网格离散化有关。然而，如果蒙皮表现良好，我们可以假设具有相似蒙皮权重的形状上的位置会发生相似的变形，因此意味着相似的最佳拟合旋转。因此，我们根据它们在权重空间中的表示对边集进行聚类：其中顶点$\mathbf{x}$q取坐标$[{\omega }_1(\mathbf{x}),{\omega }_2(\mathbf{x}),...,{\omega }_n(\mathbf{x})]$.聚类边集的数量显示变形质量的收益递减，因此我们可以选择少量的聚类，与蒙皮权重函数的数量（而不是离散网格顶点的数量）成正比。
这个提出的变形模型可以同时被视为 ARAP 的快速子空间优化和寻找最佳蒙皮变换自由度的自动方法。
通过与句柄相关的蒙皮转换上的线性等式约束支持各种用户界面。 要完全修复转换，我们只需添加约束：
$\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right){\mathbf{T}}_i^T={\hat{\mathbf{T}}}_i^T$,
其中${\hat{\mathbf{T}}}_i^T$是表示手柄的转置固定变换的$(3+1)\times 3$的矩阵
为了仅修复句柄的原点，我们添加了一个约束，要求转换以插入空间中的一个点（通常是所有满足${\omega }_i=0$点的质心):
${\mathbf{c}}^{\prime T}{\mathbf{T}}_i^T={\mathbf{c}}^T$,
其中${\mathbf{c}}^T$是转置齐次坐标中其他点的位置，而${\mathbf{c}}^{\prime T}$是用户提供的顶点信息。
我们可以类似地在句柄处仅修复转换的线性部分，释放平移组件（产生“鸡头”效果）：
$\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\end{array}\right){\mathbf{T}}_i^T={\hat{\mathbf{L}}}_i^T$,
其中${\hat{\mathbf{L}}}_i^T$是手柄$i$处变换的固定线性部分的一个$3\times 3$的矩阵。
最后，我们可以允许用户完全释放转换的自由度，委托优化以找到所有元素的最佳值。
为此，我们只需避免添加相应的约束。

带分组边集的 Arap

作为一种子空间方法，一个直接的缺点是降低了自由度。这会带来性能，但在某些情况下会过多地限制行为。在这种情况下，可以使用蒙皮子空间为传统的 ARAP 优化构建旋转边缘集的有效聚类：放弃子空间替换。这具有双重效果。旋转拟合、局部步长的成本大大降低，并且变形根据聚类“正则化”。从高层次的角度来看，如果集群是从蒙皮权重派生的，那么它们将阻止弯曲，尤其是沿着权重函数的等值线。如果事先不知道句柄，也可以根据“测地线嵌入”进行聚类，例如双谐波距离嵌入。
有鉴于此，我们可以将“轮辐+轮辋”样式的表面 ARAP 视为每个三角形边集的（轻微且冗余的）聚类。

双谐波坐标

线性混合蒙皮（如上）通过将句柄（骨骼、点、区域等）处的完整仿射变换通过权重传播到形状的其余部分来使网格变形。另一个变形框架，称为“广义重心坐标”，是线性混合蒙皮的一个特例：转换仅限于纯平移，并且需要权重来保持仿射精度。后一个要求意味着我们可以将网格中任何顶点的静止位置写为控制手柄位置的加权组合：
$\mathbf{x}={\sum }_{i=1}^m{\omega }_i(\mathbf{x})\ast {\mathbf{c}}_i$,
其中${\mathbf{c}}_i$是第 $i$个控制点的静止位置。这简化了运行时的变形公式。我们可以简单地将形状的每个点的新位置作为转换后的控制点位置的加权组合:
${\mathbf{x}}^{\prime }={\sum }_{i=1}^m{\omega }_i(\mathbf{x})\ast {\mathbf{c}}_i^{\prime }$,
“广义重心坐标”有许多不同的风格（参见“自动方法”部分中的表格）。“广义重心坐标”的模糊目标是为更大的点集或多面体捕获尽可能多的单纯重心坐标的属性（例如，用于 2D 中的三角形和 3D 中的四面体）。一些广义重心坐标可以以封闭形式计算； 其他需要基于优化的预计算。几乎所有风格都需要连接信息来描述控制点如何围绕输入形状形成外部多面体：笼子。然而，最近的一项技术不需要笼子 23。这种方法确保在优化过程中对具有仿射函数的平滑能量的权重进行仿射精度：
$mi{n}_{\mathbf{w}}trace(\frac{1}{2}{\mathbf{W}}^T\mathbf{A}\mathbf{W})$,受制于：$\mathbf{C}=\mathbf{W}\mathbf{C}$
受选定顶点的插值约束。如果 $\mathbf{A}$在其内核中具有仿射函数，也就是说$\mathbf{A}\mathbf{V}=0$，那么权重 $\mathbf{W}$将保持仿射精度，我们将得到：
$\mathbf{V}=\mathbf{W}\mathbf{C}$
上述等式的矩阵形式。定义$\mathbf{A}$的建议方法是构造一个矩阵$\mathbf{K}$，它在所有内部顶点和所有边界顶点处测量拉普拉斯算子。离散拉普拉斯算子的通常定义（例如 libigl 从 igl::cotmatrix 返回的内容）测量内部顶点函数的拉普拉斯算子，但测量函数的拉普拉斯算子减去边界顶点函数的法线导数。因此，我们可以让：
$\mathbf{K}=\mathbf{L}+\mathbf{N}$
其中$\mathbf{L}$是通常的拉普拉斯算子，$\mathbf{N}$是计算网格边界顶点处的分段线性函数的法向导数的矩阵。然后将$\mathbf{A}$视为二次形式，计算$\mathbf{K}$的积分平均值的平方，并将其应用于函数并在网格上积分：
$\mathbf{A}=({\mathbf{M}}^{-1}\mathbf{K}{)}_{\mathbf{M}}^2={\mathbf{K}}^{\mathbf{T}}{\mathbf{M}}^{-1}\mathbf{K}$
由于拉普拉斯 $\mathbf{K}$ 是二阶导数，它在仿射函数上测量为零，因此 $\mathbf{A}$ 在其零空间中具有仿射函数。一个简短的推导证明这意味着 W 将是仿射精确的（参见 23）。
这种“平方拉普拉斯”能量的最小化器在某种意义上是离散的双谐波函数。因此它们被称为“双调和坐标”（与有界双调和权重不同，后者不是广义的重心坐标）。
在 libigl 中，可以计算给定网格 (V,F) 和选定控制点或控制区域列表 S 的双调和坐标（其作用类似于蒙皮手柄）：

igl::biharmonic_coordinates(V,F,S,W);

示例407显示了粗橙色网格上的物理模拟。 该网格的顶点成为蓝色高分辨率网格的双调和坐标变形的控制点。

直接Delta Mush

为了产生平滑的变形，线性混合蒙皮需要平滑蒙皮权重。这些可以手动绘制或自动计算（例如，使用有界双谐波权重 ）。即便如此，线性混合蒙皮由于其固有的线性（见上文）而遭受收缩和塌陷伪影。
“Direct Delta Mush” 18 蒙皮尝试通过提供直接蒙皮方法来解决这两个问题，该方法将具有分段常数权重函数的装备作为输入（权重在任何地方都是 = 0 或 =1）。直接 delta mush 是一种性能较低的方法的改编，简称为“Delta Mush”19。Delta Mush 的计算分为“绑定姿势”预计算和运行时评估。
在绑定时，对绑定姿势进行拉普拉斯平滑，将每个顶点从其静止位置${\mathbf{v}}_i$移动到新位置${\tilde{\mathbf{v}}}_i$。描述撤消此平滑过程的“增量”被计算并存储在与顶点关联的局部坐标系中：
${\delta }_i={\mathbf{T}}_i^{-1}({\mathbf{v}}_i-{\tilde{\mathbf{v}}}_i)$.
在运行时，使用线性混合蒙皮和分段常数权重对网格进行变形。在骨骼附近，变形是完全刚性的，而在骨骼相遇的关节附近，网格会随着突然变化到下一个刚性变换而撕裂。在运行时将相同数量的拉普拉斯平滑应用到这个姿势网格。将每个顶点移动到一个位置 ${\tilde{\mathbf{u}}}_i$ 。在此位置计算本地帧 ${\mathbf{S}}_i$ ，并在此解析帧中添加缓存的增量以恢复形状的原始细节：
${\mathbf{u}}_i={\tilde{\mathbf{u}}}_i+{\mathbf{S}}_i{\delta }_i$.
“Delta Mush”的关键见解是拉普拉斯平滑对其余模型和姿势模型的作用类似。
“Direct Delta Mush”的关键见解是，运行时的拉普拉斯平滑过程几乎是线性的，并且可以使用 SVD（参见 ARAP）以令人惊讶的并行方式计算局部帧。直接 delta mush 将平滑步骤移动到预计算中，导致每个顶点每个骨骼的“向量值”蒙皮权重，存储在矩阵 $\Omega$ 中。在 libigl 中，对于网格 $(V,F)$ 和（例如，分段常数）权重 $W$，此预计算是使用：

igl::direct_delta_mush_precomputation(V, F,Wsparse, p, lambda, kappa, alpha, Omega);

参数 p、lambda、kappa、alpha 控制所产生变形的平滑度和紧凑度。 预计算的输出是矩阵 Omega。
在运行时，$\Omega$ 用于将网格变形到其最终位置。 在 libigl 中，这是使用以下方法计算的：

igl::direct_delta_mush(V, T_list, Omega, U);

其中 T_list 是与每个骨骼相关的输入姿势（仿射）变换，最终位置存储在 $U$ 中。

使用 Kelvinlet 进行网格变形

Kelvinlets24 是一种基于物理的实时虚拟弹性材料体积雕刻技术。 该技术将网格视为由可压缩材料制成的流体，并通过沿位移场的平流点使它们变形。 它依赖于弹性方程的解析解。

线性弹性静力学快速入门

线弹性的平衡状态由使弹性势能最小的位移场决定 $\mathbf{u}:R^3\to R^3$,
$E(\mathbf{u})=\frac{\mu }{2}||\nabla \mathbf{u}|{|}^2+\frac{\mu }{2(1-2\nu )}||\nabla \cdot \mathbf{u}|{|}^2-<\mathbf{b},\mathbf{u}>$
其中 $\mu$ 是弹性剪切模量，$\nu$ 是泊松比，而 $\mathbf{b}$ 代表外力。
第一项控制位移场的平滑度，第二项惩罚无限小的体积变化，最后一项表示要抵消的外部物体力。
可以将最优位移场与上述方程临界点的解相关联，也称为 Navier-Cauchy 方程：
$\mu \Delta \mathbf{u}+\frac{\mu }{1-2\nu }\nabla (\nabla \cdot \mathbf{u})+\mathbf{b}=0$
Kelvinlet 是 在由于点 ${\mathbf{x}}_0$ 处的力矢量 $\mathbf{f}$ 导致体载荷集中的情况下,Navier-Cauchy 方程的解，即 $\mathbf{b}(\mathbf{x})=\mathbf{f}\delta (\mathbf{x}-{\mathbf{x}}_0)$ ，方程可以重写为：
$\mathbf{u}(\mathbf{r})=[\frac{(a-b)}{r}\mathbf{I}+\frac{b}{r^3}{\mathbf{r}\mathbf{r}}^t]\mathbf{f}\equiv \mathbf{K}(\mathbf{r})\mathbf{f}$
其中 $\mathbf{K}$ 是Kelvinlet函数，${\mathbf{r}=\mathbf{x}-\mathbf{x}}_0$ 是从负载位置 ${\mathbf{x}}_0$ 到观察点 $\mathbf{x}$ 的相对位置向量，$r$ 是 $\mathbf{r}$的范数，$a=\frac{1}{4\pi \mu }$，$b=\frac{a}{1-\nu }$
位移场 ${\mathbf{u}(\mathbf{x}-\mathbf{x}}_0)$ 使线弹性材料中的一个点 $\mathbf{x}$ 变形为 ${\mathbf{x}+\mathbf{u}(\mathbf{x}-\mathbf{x}}_0)$。
然后，相关的变形梯度由以下形式的 $3\times 3$ 的矩阵定义：
${\mathbf{G}(\mathbf{x}-\mathbf{x}}_0)=\mathbf{I}+\nabla {\mathbf{u}(\mathbf{x}-\mathbf{x}}_0)$.
这个梯度 $\mathbf{G}(\mathbf{r})$ 决定了位移场 $\mathbf{u}(\mathbf{r})$ 的不同性质。例如，$\nabla \mathbf{u}(\mathbf{r})$ 的斜对称部分表示由 $\mathbf{u}(\mathbf{r})$ 引起的旋转，而其对称部分对应于弹性应变并决定了拉伸。应变张量也可以分解为表示弹性介质体积缩放的有迹项和表示收缩变形的无迹项。
这构成了 Kelvinlet 画笔的基本原理。

正则化Kelvinlet

单个点 ${\mathbf{x}}_0$ 的集中体载荷在 ${\mathbf{x}}_0$ 点为 Kelvinlet 解引入了奇点。因此，将Kelvinlet方程修改为：
${\mathbf{u}}_{ϵ}(\mathbf{r})=[\frac{(a-b)}{r_{ϵ}}\mathbf{I}+\frac{b}{r_{ϵ}^3}{\mathbf{r}}_{ϵ}{\mathbf{r}}_{ϵ}^t+\frac{a}{2}\frac{{ϵ}^2}{r_{ϵ}^3}\mathbf{I}]\equiv {\mathbf{K}}_{ϵ}(\mathbf{r})\mathbf{f}$
其中$r_{ϵ}=\sqrt{r^2+{ϵ}^2}$是正则化的距离，$ϵ>0$ 是雕刻刷的半径。
因此，给定一个力矢量 $\mathbf{f}$，可以计算 $R^3$ 中任意点的位移，从而定义基于物理的空间变形。在实践中，力矢量根据画笔尖端位移 $\bar{\mathbf{u}}$ 参数化。为此，我们可以用约束 ${\mathbf{u}}_{ϵ}(0)=\bar{\mathbf{u}}$ 展开正则化 kelvinlet 方程以得到 ${\mathbf{u}}_{ϵ}=cϵ{\mathbf{K}}_{ϵ}(\mathbf{r})\bar{\mathbf{u}}$，其中 $c=2/(3a-2b)$
不同径向尺度的 Kelvinlets 可以线性组合以构建具有任意快速衰减的画笔：
${\mathbf{u}}_{{ϵ}_1,..{ϵ}_n}(\mathbf{r})=c({\sum }_i\frac{{\omega }_i}{{ϵ}_i}{)}^{-1}[{\sum }_i{\omega }_i{\mathbf{K}}_{ϵi}(\mathbf{r})]\bar{\mathbf{u}}$
其中，${\omega }_i$是权重因子并且有 ${ϵ}_i<{ϵ}_{i+1}$。由于叠加原理，这些复合画笔仍然满足 Navier-Cauchy 方程。
正则化 kelvinlets 可以通过用基于矩阵的分布替换基于向量的负载分布来进一步扩展，以实现如前所述的扭曲、收缩和缩放等非仿射变换。
在 libigl 中，这是使用以下方法计算的：

igl::KelvinletParams<double> brushParams{brushRadius, scale, brushType};
igl::kelvinlets(V, origin, forceVec, forceMatrix, brushParams, result);

其中brushRadius、scale、brushType对应于 $ϵ$、衰减和操作（抓取、捏合、缩放、扭曲）。
(Example 409) pinch, twist, grab and scale in action

.