数字几何处理

点云#

有符号距离场#

隐式函数#

空间网格与层次结构#

四叉树，八叉树-hierarchy结构

如何构造数据结构呢？

对立方体八叉树的一个优化：

用平面对曲面进行逼近，如果逼近程度已经很好了，就不用继续细分了。得到的是一些平面的块？意思是直接用平面逼近曲面

三角网格表示#

流形#

半边结构#

OFF格式

思考#

法向#

梯度离散化#

计算梯度

可以看到这个梯度的值和x的值是无关的，也就是每个三角形面中的梯度都是一样的constant，f(x)是一个线性函数

拉普拉斯算子#

每个三角形内部的梯度是一个constant，而每个顶点所连接的三角形的梯度是不一样的，因此可以对顶点求梯度的梯度

均匀 Laplace 算子#

余切 Laplace 算子#

利用格林公式，把面积分转化为线积分

知道为什么是除以4了，因为有一个除以2是求得。蓝色区域的面积，而不是某一个点的一阶邻域的总面积。 那么余切是怎么来的：上面（xi-xk）(xj-xk)点乘是cos，而下面的AT是1/2ab sin(),上下一比就是cot

公式里有两部分的cot，是因为一条边可能对应的是两个三角形的cot

高斯曲率#

平均曲率就是（最大曲率加上最小曲率）/2,H是这个点的平均曲率,两个主曲率除以2，δx就是算出来的离散拉普拉斯算子

高斯曲率是主曲率相乘(k1 * k2)，平均曲率是相加除以2(k1 + k2) / 2

基于滤波的方法#

空间和时间的离散：

时间上：由于上面给出的扩散方程，所以下面这个式子成立

曲率和laplace转换公式的条件：需要时cot权重

uniform的laplace不满足线性LIN性质的，因为他是把每个点都移动到了一阶邻域的中心

网格光顺（Fairing）#

高斯滤波#

这个Kp是为了实现一个归一化，Ws权值归一化把p和q是度量两点之间的距离权重

不能保边缘

解决方法：双边滤波，多了一个对于颜色跳变权重的度量，考虑颜色上的相似性

加入对颜色跳变的估计

如何把这个对像素的处理，推广到对网格上的处理？，高度 + 距离

去找一个顶点，找到法向，做一个切平面，对一阶邻域中的点求出高度，作为双边滤波中，二维像素的灰度值

类比于对图像处理的ws，wr

法向滤波与顶点位置恢复#

那直接对曲面上的点的法向量做处理，然后对点的位置进行恢复，这样会很不精准。 因此，想要从别的信号恢复点的位置，瞄准了三角形的法向量是无歧义的，从这个法向量进行恢复

这样可以算出来新的法向量，更加自然

傅里叶变换#

本质是空域和时域之间的转化

如何推广到二维流形曲面？

效果：

相关论文#

效果：

技术实现#

要使输出和输入之间的距离最小，然后输出的梯度为0

beta是自己设置的，需要逐渐变大

为什么1方法会存在解析解？

得到一种分片，某种颜色比较光滑的结果

推广到三角网格#

怎么做？

cotangent laplace效果不好的原因是 他没有把flat rigion给抹平

如何推广定义一条边的laplace算子？

基于面积的边算子#

问题：

基于数据驱动的方法#

data-driven methods

通过拟合函数自适应地去找一个合适的结果,找s，r，k

参数化方法#

Tutte 重心映射#

LSCM#

不能保证双射和不自交

conformal-》保角

在每个输入的三角形网格建立一个局部的坐标系

这个工作的目标：是找一套uv坐标，使得空间三角网格映射到二维坐标上的畸变程度最小

参数化坐标（u，v）是定义在三角形网格上的一个函数，也就是一个分片线性函数

由于中心坐标表达，三角形上每个点的位置都可以有三个顶点的重心坐标公式表示出来，因此是线性的

把uv映射看成函数，ui = u（x），由于三角形中心坐标那里，可以得出u的梯度

要使这个连续函数是一个相似变换： scale 乘上一个旋转矩阵？这是目标，找出共性变换的

由于3D曲面有高斯曲率（比如一个球帽），想把它完美地无畸变地摊平成平面是不可能的（除非把它撕开）。

ABF#

Angle-Based flattening

约束：

最后一个约束，为了保证最后一个三角形放进来之后可以使得整个曲面合得上，所以根据正弦定理，可以得到这个约束式子

这个式子的理论基础是正弦定理：

推导

非线性等式约束，不是很好解

解决方案：

ARAP参数化#

变形类型

Local/Global 求解#

为什么能收敛：

每次求解，Jt到Lt的距离都在减少，能量方程是单调下降的，所以是收敛的

相关论文#

思路：

效果：

优化：

形变区域与约束#

基于传播的方法（Propagation Based）#

黑色红色线是等值线，怎么算扩散？

从handle(1)到fix(0)，中间的变形区域设为s，那么如何计算s。实际上变形区域内部的点就是x = x + s * displacement(位移，绿色那部分的位移) 方法：

线性差值距离。

方法2：直接解出来s，通过laplace矩阵。R上面的那些顶点相当于是变量，H上设置为1？ 接一个laplace方程，上面矩阵保证顶点关系，多出来的约束用做是变形后的顶点的位置？ 接出来的应该是s（距离值吧），那就把变形对应的点右侧b设置为1，其他都是0，然后解出来这个矩阵方程，得到的每个点的解就应该是距离handle那一圈的距离s的值

问题：

这里中间不能保证是大于1的，会往下凹。

多尺度方法（Multi-Scale Based）#

多尺度的方法

分成两（或者多）个频率带，一个是低频的，一个是高频的。直接把两个部分分开，把低频部分做一个deform，然后直接把高频部分加回来。因为高频部分可能会有一些细节，对低频部分做一个操作，细节比较少。高频加回来之后细节扭曲比较少

步骤：

如何表示这个滤波： base就是滤波光滑后网格顶点位置，h（details）是每个顶点上的displacement（位移）

hi不是正常的在三维空间中的表示，而是用局部的标准正交基构成的，local frame，ni是base曲面的法向量，ti1，ti2是垂直于ni的两个向量

alpha beta gamma三个数字是不变的

微分坐标（Differential Coordinates）#

基于微分坐标

先恢复出来微分坐标（梯度或者laplace）

过程是：先计算出来原始曲面的梯度，然后根据这个梯度去恢复出来变形后的曲面的位置

步骤：

如何算这个Mi呢？相当于用户操作的形变的表达

Mi一开始由用户给定，现在的问题变成如何把Mi从一开始扩散到网格外面

这是一个默认的好用的方法：直接对r和s进行差值

对平移的敏感程度不好

下面这个方法是点变换的同时对rotation矩阵Mi持续优化

形变迁移（Deformation Transfer）#

变形参考，把马的变形转移到骆驼上去，对比着transfer，有相似的变形量

第一步要算出来deformation的表示

只有三个点的话是不足以算出来3个点的仿射变换的，所以在空间中加一个v4

ARAP 去噪#

希望：变形后的那个向量和变形前的那个向量只经过一个旋转

这个第一个式子，因为前后两项跟变量Ri没有关系，Ri是local中要求的旋转矩阵，所以直接不考虑，只考虑带有Ri的项。 矩阵的迹是一个数字，它等于主对角线上所有元素的和。所以第二行可以这么写。

想象一下(1, 2, 3) * (1, 2, 3)T 是一个数字，把它转置一下就是(1, 2, 3)T * (1, 2, 3),这是一个矩阵，迹就是原来的数字

课堂小测：推导

空间形变#

基于笼的自由形变#

Cl是网格的控制顶点，φl是第l个顶点的基函数，这个奇函数在l上pi那一点的取值

问题：#

重心坐标#

均值坐标（Mean Value Coordinates）#

坐标推导#

条件：

证明：

效果#

坐标权重效果#

负数角度解决方法：

Harmonic Coordinates#

harmonic 调和坐标

学思路：层次结构/多分辨率结构，从下往上迭代，为了在求解过程中及早得陷入一个局部最小。

Bounded Biharmonic Weights#

什么是实时的？如果要解一个方程那就肯定不是实时的，而基于这种wij，或者linear blend skinning的方式，算是实时的，速度快

三类形变控制方式#

修复问题分类#

经常出现的数据类型#

修复方法#

法向处理

补洞

合缝

拓扑简化

找到环丙-切开网格，补上洞

嵌入#

嵌入到体素会出现的问题：

解决方法：

映射#

度量#

优势#

要取得的是黄色圆盘的外部

投影方法#

PolyCube#

定义就是横平竖直的，边界的那个法向量要和轴对齐，轴有6个方向那种

PolyCube 生成方法#

流程

把表面上的顶点旋转矩阵扩散到内部

基于体素的方法#

人为地搭建出来一个polycube

新的方式：

先把多边形网格反映射回到rotation-dirven模型上，再头回到原始网格上进行mapping

条件：

基于聚类的方法#

复数角

欧拉角

曲面映射问题#

构造polycube和原始输入模型之间的映射

连接关系，顶点个数，三角形个数一样，顶点位置不一样，这不一定是双射的，比如说有个三角形退化了或者自交

算法：#

几何渐变（Morphing）#

形状差值的过程

要求是compatible meshes，就是顶点个数和连接关系是一样的，只不过顶点位置不一样

作业9#

实现morphing

性质：#

ARAP插值#

效果：

输入：compatible input output meshes

数据驱动方法#

找到最短路径，直接进行一个arap

点集配准#

Atlas 生成#

什么是Atlas genertaion？

纹理和obj对应

一个缝对应的两块之间缝合后的颜色不一致性问题

应用：

网格切割方法#

point to path文章：

基于segmentation方法：

Bijection#

scaffold(脚手架方法)

轴向形变#

分解与装箱#

畸变优化#

altas 完结

简化目标#

基本定义#

局部操作#

二次误差度量（QEM）#

二次误差度量

作业9：#

如果得到的quadratic矩阵不是满秩的，那就直接用v1v2这条边的中点来代替，最小化顶点应该也在这条边上

把error放到最小堆里，然后每次把cost最小的边拿出来。每次更新后都要更新涉及到更新的边的cost，新顶点的cost‘就继承到了原来的那个Q矩阵

书上115页求解最优解

VSA#

把网格分割得到region，每个region用平面去逼近

球面参数化问题#

把一个亏格为0的映射到一个球面上

应用：

challenge：

如何实现球面参数化：

层次化方法#

网格细化#

方向场定义#

定义

方向场中，每个方向都相差90°？

作业10#

离散化#

方向场表示#

方向场的表示：

三角函数

作业#

基本性质#

约束#

对齐

对称

可积场#

参数化的梯度可以看做是两个分离的向量场

标架场（Frame Field）#

这个polar decomposition 把frame filed转为cross filed 通过向量的形式

Delaunay 条件#

生成方法#

proof：

全部点都做一遍flip？

最大化最小角

直接比角度，去代替做外接圆看另一个点是不是在外接圆外

最优 Delaunay 三角剖分#

这里会动这个顶点位置

很有用的方法（上面的图）#

求中间变量，再求目标变量

作业3#

对于rock这个obj，如果不加这两行处理，会导致部分曲面无法显示

1. 你没对向量归一化，dot(a,b) 可能远大于 1。sqrt(1 - cos^2) 里出现负数 → NaN 或极小值，cot = cos / sin 爆炸，顶点位移失控，面退化。
2. 加 std::min(…,0.9) 把过大的 dot 强行压到 0.9，避免 1 - cos^2 变成负数或极小，数值不再爆炸，所以“好了”。

作业4#

边界点没有按照边界的拓扑顺序排序，导致参数化之后，连接三角面时开始交叉

如何找边界点的拓扑顺序呢？

1
for (auto& v : mesh.vertices) {
2
  v->boundary_index = -1; // 重置
3
}
4
// 遍历所有边界半边，标记边界顶点的 boundary_index
5
// 按拓扑顺序遍历边界环** 很重要
6
int boundary_size = 0;          // 总边界顶点数
7
int boundary_cycle = 0;         // 边界环个数
8
std::unordered_set<geometry::HalfEdge*> visitedBoundaryHE;// 访问过的边界半边
9
// 遍历所有半边，找到边界半边
10
for (auto& hePtr : mesh.halfEdges) {
11
  geometry::HalfEdge* start = hePtr.get();
12
  if (start->pair) continue;// 不是边界半边
13
  if (visitedBoundaryHE.count(start)) continue;// 如果遍历过了，就跳过
14

15
  geometry::HalfEdge* he = start;
16
  int step = 0;
17
  do {
18
    visitedBoundaryHE.insert(he);
19
    geometry::Vertex* v = he->vertex;
20
    if (v->boundary_index < 0) {
21
      v->boundary_index = boundary_size++; // 从 0 开始递增
22
    }
23

24
    // 跳到下一条边界半边
25
    geometry::HalfEdge* candidate = he->next;
26
    while (candidate && candidate->pair) {
27
      candidate = candidate->pair->next;
28
    }
29
    he = candidate;
30
    step++;
31
    if (step > (int)mesh.halfEdges.size()) { // 安全退出：异常拓扑
32
      std::cerr << "Warning: boundary walk aborted (non-manifold?)" << std::endl;
33
      break;
34
    }
35
  } while (he && he != start);
36
  boundary_cycle++;
37
}
38

39
std::cout << "Boundary cycles: " << boundary_cycle
40
  << " boundary vertices: " << boundary_size << std::endl;

作业5：#

1
1.  你构造的是标准 cotangent 拉普拉斯 L：对角等于该行所有邻边权值之和，邻接是负权值，故每行元素和为 0 → 常量向量 (1,1,...,1) 满足 L·c=0，形成零空间。
2
2.  如果网格有多个连通分量，每个分量的常量向量都在零空间 → 零空间维度≥分量数。
3
3.  跳过退化三角形或边界半边会让某些顶点行近似全零或只剩很小对角，增加额外近似零方向。
4
4.  未加任何 Dirichlet 约束（固定顶点）或对角正则前，平移自由度未被去除，因此 L 不可逆。
5
5.  若再存在孤立顶点（度=0），对应行全 0，直接再增一维零空间。
6
解决：固定至少一个（更常用：一圈边界）顶点行或加极小正则 εI；有多分量需各自锚点。

1
2.  Laplacian A 构造：跳过 boundary 半边 (ithe->isBoundary()) 使某些仅与边界接触、拓扑度很低的顶点行几乎全 0（只被正则 1e-8 微调）。SparseLU 在近奇异情况下把它们解成接近 0 的数 → 出现从原点辐射的线。解决：对所有 boundary 顶点加入硬 Dirichlet 行 (A(ii)=1, b=原参数化坐标)，或只排除 boundary 顶点出线性系统。
2

3
这行 continue 的意思是：“如果这是一条边界边，就完全忽略它，不要为它计算任何权重，也不要在矩阵 A 中建立任何连接。”
4

5
这导致了什么后果？
6

7
边界顶点在矩阵A中是“孤立”的：考虑一个边界顶点。它连接的所有边都是边界边（除了极少数情况）。因为您的代码跳过了所有边界边，所以在构建矩阵A时，这个边界顶点与它的任何邻居之间都没有建立连接！
8

9
软约束无力回天：之后，您对边界点施加了一个软约束 A.coeffRef(i, i) += 1e-8;。此时，一个边界顶点 i 在矩阵 A 中对应的整行，除了对角线上的 1e-8 之外，几乎全是零！
10

11
方程爆炸：这意味着，对于一个边界顶点 i，它在全局步骤中需要求解的方程，实际上简化成了： 1e-8 * u_i = b_i 这里的 b_i 是从局部步骤计算出的一个非零的“理想拉力”。为了解出 u_i，系统只能计算 u_i = b_i / 1e-8。这是一个极其巨大的数值！
12

13
视觉结果：因此，所有边界顶点都被计算到了一个离原点非常遥远的位置。而内部顶点由于受到正确的约束，位置相对正常。最终，您就看到了从内部“飞”向遥远边界的那些奇怪直线。

1
3. “累加”的魔力：从“贡献”到完整的矩阵
2
最关键的一点是，这个过程会对网格中的每一条内部边都执行一次。
3

4
想象一下顶点 i，它有三个邻居 j, k, l。
5

6
当循环处理边 (i, j) 时，A(i, i) 会被加上 w_ij。
7

8
当循环处理边 (i, k) 时，A(i, i) 会被再次加上 w_ik。
9

10
当循环处理边 (i, l) 时，A(i, i) 会被第三次加上 w_il。
11

12
当所有边都处理完毕，A.setFromTriplets(...) 会执行累加。最终，矩阵 A 中 A(i, i) 的值，就是所有与顶点 i 相连的边的 cotangent 权重之和：w_ij + w_ik + w_il。

1
那现在后面那个b矩阵我没有理解，为什么是要这么做：
2

3
就算你说参考了力学的作用，对每个顶点有这样一个力，我还是没理解为什么要这么写得出来的结果是arap的？为什么要用cots，global我知道是旋转矩阵，那为什么能根据这个公式得到b，b是什么？b是拍平后的3维三角形，每个向量的大小，乘以到tutte参数化上的那个旋转矩阵，乘以cots？为什么这么写，cots是什么意义，我这么乘是啥意思，这不是有缩放了吗，我不是要保证尽可能isometric吗，我认为右边就是一个缩放矩阵就行了

1
1. b 是什么？—— “理想的内部拉力”集合
2
首先，也是最重要的一点：b 向量不是坐标，也不是一个变换矩阵。
3

4
您可以把 b 向量想象成一个列表，记录了施加在每一个顶点上的一个**“理想的合力”**。整个全局步骤的目标就是求解一个巨大的平衡方程：A * u = b
5

6
u: 我们想要求解的、新的、更好的顶点位置。
7

8
A: 我们已经知道，这是余切拉普拉斯矩阵。它可以被看作是描述当前网格 u 内部弹力的系统。A*u 计算出的就是当前形状 u 下，每个顶点受到的内部合力。
9

10
b: 这是我们期望的、理想状态下每个顶点应该受到的内部合力。
11

12
所以，求解 A*u = b 的物理意义就是：“请找到这样一个新的形状 u，使得它自身的内部弹力 (A*u)，与我们计算出的‘理想拉力’ (b) 完全相等。”

1
2. b 的构建：一步步拆解 cots * global * e
2
现在我们来拆解这个“理想拉力”是如何计算出来的。它分为三步，对应三个乘法：
3

4
第一步: e —— “完美”的边 (e01, e12, e21)
5
是什么？ e01, e12, e21 是从 local_coords 中重建出来的2D边向量。而 local_coords 来自于对原始3D模型的“拍平”。
6
意义？ e 代表了在没有任何变形、拉伸或压缩的情况下，这条边最完美的、最原始的2D形态。它是我们追求“保长度”的**“黄金标准”**。
7

8
第二步: global[i] * e —— “理想”的边
9
是什么？ global[i] 是我们在局部步骤中为这个三角形计算出的最佳旋转矩阵 R。
10
意义？ R * e 的意思是：“取那条‘完美’的边 e，将它旋转到当前参数化平面上最合适的位置。”
11
结果？ R * e 是一个**“理想的边向量”。它具有我们想要的完美的长度**（来自e）和理想的方向（来自R）。
12

13
第三步: cots * (global[i] * e) —— “有影响力的理想拉力”
14
是什么？ cots 是我们计算出的余切权重。
15
为什么要有它？ 这就是最核心、也最违反直觉的一步。cots 在这里不是为了缩放几何，而是为了赋予这个“理想拉力”一个正确的“影响力”或“权重”。
16
意义？ cots 来自ARAP能量函数的数学推导（求导）。
17
ARAP的能量衡量了当前形状与理想形状的差异。
18
这个能量公式中，余切权重 cots 自然地作为加权因子出现。一个形状良好、面积大的三角形，其cotangent值通常更大，它对总能量的贡献也更大。这很合理：我们应该更关注那些“重要”的三角形。
19
当我们对这个能量函数求导以找到最小值时，cots 这个权重项就同时保留在了等式的左边（构成了矩阵 A）和右边（构成了向量 b）。
20
所以，b 向量必须乘以 cots，不是为了改变几何，而是为了与矩阵 A 的构建方式保持数学上的一致性！ 如果 b 不乘以 cots，那么 A*u=b 这个方程本身就是错误的，它不再是原始能量函数的正确梯度。

1
解答您的疑问
2
现在，我们可以直接回答您的问题了：
3

4
“为什么要用cots？”
5
因为 cots 是ARAP能量函数中的自然权重。为了让我们的线性系统 A*u=b 成为该能量函数的正确梯度，cots 必须同时作为权重出现在矩阵A和向量b中。它衡量的是每条边在整个系统中的“重要性”。
6

7
“我这么乘是啥意思，这不是有缩放了吗？”
8
您看得完全正确，这是一个缩放。但关键在于，我们缩放的不是几何本身，而是这个几何元素（边）在全局平衡方程中的**“话语权”或“影响力”**。一条权重（cots）大的边，它产生的“理想拉力”就更大，求解器在寻找最终平衡点时，就会更“努力”地去满足这条边的长度和方向要求。
9

10
“我不是要保证尽可能isometric吗？”
11
是的！而 A*u=b 这个方程，正是实现这个目标的数学工具。求解这个方程得到的新 u，其每个三角形的形状会趋向于 global[i] * e 所定义的那个只有旋转、没有缩放的“理想形状”。cots 在这个过程中扮演的角色，是告诉求解器在成千上万个三角形中，应该优先满足哪些三角形的“保长度”要求。
12

13
“我认为右边就是一个缩放矩阵就行了”
14
这里的混淆点在于，b 不是一个作用于顶点的变换。它是一个目标值。整个过程更像是一个巨大的“最小二乘”拟合问题：我们有一个描述内部连接的矩阵A，我们想找到一个位置向量u，使得A*u的结果尽可能地接近我们的目标b。而这个b，正是由“旋转过的完美边”和“边的重要性(cots)”共同决定的。
15

16
希望这个从物理意义到数学原理的拆解，能够彻底解开您关于 b 向量的疑惑！这是一个非常深刻的问题，理解了它，就真正理解了ARAP算法的精髓。

1
ARAP能量函数与梯度的推导
2

3
您的问题非常深刻：Au = b 这个线性系统，是如何从“尽可能保持等距(isometric)”这个目标推导出来的？为什么cotangent权重会同时出现在A和b中？
4

5
答案在于，ARAP算法的核心是定义一个**“变形能量” (Deformation Energy)**，然后通过最小化这个能量来找到最佳的2D顶点位置。我们求解的线性系统Au=b，正是这个能量函数梯度为零的数学表达。
6

7
1. 核心思想：定义“变形能量”
8

9
ARAP的能量衡量的是当前2D形状与理想2D形状之间的差异。
10

11
当前2D形状: 由我们要求解的变量——顶点位置 u 决定。
12

13
理想2D形状: 它就是原始3D三角形经过一次纯旋转后得到的形状。它具有完美的边长和角度，只是方向变了。
14

15
因此，ARAP的总能量 E(U) 是所有三角形的“局部变形能量” E_T 的总和：
16

17
$$E(U) = \sum_{T \in \text{Mesh}} E_T$$
18

19
而每个三角形 T 的局部变形能量，是其三条边的“变形程度”的加权和：
20

21
$$E_T(u_i, u_j, u_k) = \sum_{e_{ij} \in T} w_{ij} \left\| (u_i - u_j) - R_T(p_i - p_j) \right\|^2$$
22

23
让我们拆解这个关键公式：
24

25
u_i, u_j: 当前2D参数化平面上的顶点坐标（这是我们的变量）。
26

27
p_i, p_j: 原始3D模型中的顶点坐标（这是常量）。
28

29
(u_i - u_j): 当前2D网格中的边向量。
30

31
(p_i - p_j): 原始3D网格中的边向量。
32

33
R_T: 在局部步骤中为这个三角形T计算出的最佳旋转矩阵。
34

35
R_T(p_i - p_j): 这就是我们说的**“理想的边向量”**——它有完美的长度和角度，只是被旋转到了最佳方向。
36

37
w_{ij}: 这就是余切权重 (cotangent weight)。它衡量了这条边在整个能量计算中的“重要性”或“影响力”。
38

39
2. 最小化能量：梯度为零
40

41
为了找到使总能量 E(U) 最小的顶点位置 U，我们需要使用微积分的基本原理：在最小值点，能量函数对所有变量的梯度（偏导数）都必须为零。
42

43
我们的变量是所有顶点的坐标 u_k (具体来说是 u_k.x 和 u_k.y)。所以，我们需要求解以下方程组：
44

45
$$\frac{\partial E(U)}{\partial u_k} = 0 \quad (\text{对于所有顶点 } k)$$
46

47
3. 开始推导：对单个顶点求偏导
48

49
我们来关注能量函数 E 对某一个特定顶点 u_k 的偏导数。
50

51
首先，注意到 u_k 只会出现在与它相邻的边的能量项中。所以我们只需要对包含 u_k 的那些能量项求导。
52

53
让我们看一条边 (k, j) 的能量项：
54
$E_{kj} = w_{kj} \left\| (u_k - u_j) - R_{kj}(p_k - p_j) \right\|^2$
55
(这里的 R_kj 是包含边 (k,j) 的某个三角形的旋转矩阵)
56

57
对这个能量项求 u_k 的偏导数，根据链式法则 (d/dx) ||v(x)||^2 = 2v(x) · (dv/dx)，我们得到：
58

59
$$\frac{\partial E_{kj}}{\partial u_k} = 2 w_{kj} \left( (u_k - u_j) - R_{kj}(p_k - p_j) \right)$$
60

61
现在，我们将所有与 u_k 相连的边的贡献加起来，并让总和为零：
62

63
$$\frac{\partial E}{\partial u_k} = \sum_{j \in N(k)} 2 w_{kj} \left( (u_k - u_j) - R_{kj}(p_k - p_j) \right) = 0$$
64

65

66
(这里的 N(k) 代表顶点 k 的所有邻居顶点 j)
67

68
4. 整理方程：分离变量 u
69

70
我们可以把上面这个方程中的 u 项和 p, R 项分开。首先，去掉常数 2：
71

72
$$\sum_{j \in N(k)} w_{kj} (u_k - u_j) - \sum_{j \in N(k)} w_{kj} R_{kj}(p_k - p_j) = 0$$
73

74
然后，移项得到我们最终的平衡方程：
75

76
$$\sum_{j \in N(k)} w_{kj} (u_k - u_j) = \sum_{j \in N(k)} w_{kj} R_{kj}(p_k - p_j)$$
77

78
这个方程必须对每一个内部顶点 k 都成立。
79

80
5. 连接数学与代码：Au = b 的诞生
81

82
现在，我们来看这个方程是如何变成您的代码的。
83

84
等式左边 (LHS): A * u
85

86
$$\text{LHS} = \sum_{j \in N(k)} w_{kj} (u_k - u_j)$$
87

88
我们把它展开：
89

90

91
$$\text{LHS} = \left(\sum_{j \in N(k)} w_{kj}\right) u_k - \sum_{j \in N(k)} w_{kj} u_j$$
92

93
这正是您的余切拉普拉斯矩阵 A 的定义！
94

95
对角线元素 A(k, k): 它的值是所有与 k 相连的边的权重之和，即 (∑ w_kj)。
96

97
非对角线元素 A(k, j): 它的值是连接 k 和 j 的那条边的权重的负数，即 -w_kj。
98

99
当矩阵 A 的第 k 行乘以顶点位置向量 u 时，得到的结果 (Au)_k 就精确地等于上面的展开式。
100

101
等式右边 (RHS): b
102

103
$$\text{RHS} = \sum_{j \in N(k)} w_{kj} R_{kj}(p_k - p_j)$$
104

105
这正是您的**“理想拉力”向量 b**！
106

107
向量 b 的第 k 个元素 b_k，就是所有与 k 相连的“理想边” R(p_k-p_j)，经过它们各自的余切权重 w_kj 加权后的总和。
108

109
您的代码通过遍历每个三角形，并将计算出的 cots * global * e 的贡献累加到相应的顶点上，最终构建出了这个 b 向量。
110

111
结论
112

113
我们已经完整地从ARAP的能量函数出发，推导出了您代码中求解的线性系统 Au = b。
114

115
这个推导过程清晰地表明：
116

117
cotangent 是能量函数中固有的权重，它衡量了每条边的重要性。
118

119
对能量函数求导（寻找最小值）时，这个 cotangent 权重 w 自然地、同时地出现在了包含变量 u 的项（最终构成了矩阵A）和不包含变量 u 的项（最终构成了向量b）中。
120

121
因此，为了正确地最小化ARAP能量，cotangent 必须同时作为构建A和b的权重。您看到的 b 向量中的 cots 乘法，不是一个随意的缩放，而是能量梯度推导的必然结果。

1
https://blog.csdn.net/why18767183086/article/details/108034725?ops_request_misc=%257B%2522request%255Fid%2522%253A%25227b52c04b5fa98b59c44a9cf535761a34%2522%252C%2522scm%2522%253A%252220140713.130102334..%2522%257D&request_id=7b52c04b5fa98b59c44a9cf535761a34&biz_id=0&utm_medium=distribute.pc_search_result.none-task-blog-2~all~sobaiduend~default-1-108034725-null-null.142^v102^pc_search_result_base2&utm_term=arap%E5%8F%82%E6%95%B0%E5%8C%96&spm=1018.2226.3001.4187

作业6#

问题：

进行变换之后，所有顶点的Z轴方向改变，全都缩到原点去了，解决：把所有的Z轴不做变换。 我记得一开始还出现了一旦拖动顶点，所有顶点都会往中心去，是因为Z轴的坐标突然迭代变很大。这是什么原因造成的？

后面一个问题：fixed顶点在变形过程中还会改变：这是因为相机的缘故，每次通过meshupdate()的方法都会重置相机位置

现在的问题是，变形过程中这个模型会发生萎缩。原因是我只做了一次迭代，没有做多次迭代，效果不到位

重点错误：#

这里做出来的重点错误就是，每次变换后，把这次变换前的坐标（上一次变换的new_position）作为old_position保存，然后进行下次迭代，这是严重的理解错误。正确做法是把模型一开始的顶点位置作为old_position保存起来，然后每次都是依据这个Old_position进行迭代的。

作业 7#

我的代码好像确实是没问题的，问题在于结果的展示有问题，理论上来讲，参数化后应该是按照原来模型的分布进行展示，而我现在的样子是把他逆时针转了75度左右。是因为我边界半边参数化做的有问题吗

下面这个是普通的tuttes embedding

正交矩阵列向量是单位向量#