您的位置: 首页  >  文章  >  Laplacian on Triangle Mesh

Laplacian on Triangle Mesh

分类: 文章 • 2025-01-05 10:34:10

Laplacian on Triangle Mesh<\center>

局部平均区域

基本思想是mesh面上的点x的值设置为其局部领域的平均。下图是几种具体的情况：

Laplacian on Triangle Mesh

其中蓝色区域为局部平均区域。

  1. Barycentric cell：三角形的重心和三角形边的中点相连
  2. Voronoi cell：三角形的外心和三角形的中点相连（对于钝角三角形可能会出现外心在外部）
  3. Mixed Voronoi cell：在Voronoi cell的基础上，将钝角三角形的外心替换成钝角对边的中点

梯度

拉普拉斯算子是梯度的散度，所以我们需要先求出梯度。

Laplacian on Triangle Mesh

属性插值函数：
f(x)=αfi+βfj+γfk f(x)=\alpha f_i+\beta f_j+\gamma f_k \\
其中：fi,fj,fkf_i,f_j,f_k是三角形三个顶点的属性值，x为三角形内部的点，α,β,γ\alpha , \beta , \gamma 是点x在三角形的重心坐标，α+β+γ=1\alpha + \beta + \gamma = 1。

梯度：
∇xf(x)=fi∇xα+fj∇xβ+fk∇xγ \nabla_xf(x)=f_i\nabla_x\alpha +f_j\nabla_x\beta +f_k\nabla_x\gamma
对于α\alpha来说：
α=AiAT=((x−xj)⋅(xk−xj)⊥∣∣xk−xj∣∣2)∣∣xk−xj∣∣22AT=(x−xj)⋅(xk−xj)⊥2AT \alpha = \frac{A_i}{A_T}=\frac{((x-x_j)\cdot\frac{(x_k-x_j)^\bot}{||x_k-x_j||_2})||x_k-x_j||_2}{2A_T}=\frac{(x-x_j)\cdot(x_k-x_j)^\bot}{2A_T}
对xx求导：
∇xα=(xk−xj)⊥2AT∇xβ=(xi−xk)⊥2AT∇xγ=(xj−xi)⊥2AT \nabla_x\alpha=\frac{(x_k-x_j)^\bot}{2A_T} \\ \nabla_x\beta=\frac{(x_i-x_k)^\bot}{2A_T} \\ \nabla_x\gamma=\frac{(x_j-x_i)^\bot}{2A_T}
得到梯度公式：
∇xf(x)=fi(xk−xj)⊥2AT+fj(xi−xk)⊥2AT+fk(xj−xi)⊥2AT \nabla_xf(x)=f_i\frac{(x_k-x_j)^\bot}{2A_T} +f_j\frac{(x_i-x_k)^\bot}{2A_T}+f_k\frac{(x_j-x_i)^\bot}{2A_T}
因为
(xk−xj)⊥+(xi−xk)⊥+(xj−xi)⊥=0 (x_k-x_j)^\bot+(x_i-x_k)^\bot+(x_j-x_i)^\bot=0
所以
∇xf(x)=(fj−fi)(xi−xk)⊥2AT+(fk−fi)(xj−xi)⊥2AT \nabla_xf(x)=(f_j-f_i)\frac{(x_i-x_k)^\bot}{2A_T}+(f_k-f_i)\frac{(x_j-x_i)^\bot}{2A_T}

离散拉普拉斯算子

均匀拉普拉斯

Δf(vi)=1∣N1(vi)∣∑vj∈N1(vi)(fj−fi) \Delta f(v_i) =\frac {1}{|N_1(v_i)|}\sum_{v_j \in N_1(v_i)}(f_j-f_i)

此公式可以用于各向同性的remesh上。

余切形式的拉普拉斯

使用混合的有限元/有限体方法（mixed finite element/finite volume method），拉普拉斯算子可以更加精确地进行离散推导。目标是在局部平均区域A_i=A(v_i)上，对逐片线性函数梯度地散度进行积分。

Laplacian on Triangle Mesh

∫AiΔfdA=∫Aidiv∇fdA=∫∂Ai(∇f)⋅nds \int_{A_i}\Delta fdA=\int_{A_i}div\nabla fdA=\int_{\partial{A_i}}(\nabla f)\cdot n ds
将积分分解到每个三角形上，并且局部的Voronoi区域通过三角形两条边的中点aa和bb，并且每个三角形中∇f(x)\nabla f(x)是常数，那么在该三角形TT上的积分为
∫∂Ai∩T∇f⋅nds=∇f⋅(a−b)⊥=12∇f(u)⋅(xj−xk)⊥ \int_{\partial{A_i\cap T}}∇f⋅nds = ∇f⋅(a−b) ^⊥= \frac{1}{2}∇f(u)⋅(x_j-x_k)^\bot

将梯度公式代入：
∫∂Ai∩T∇f⋅nds=(fj−fi)(xi−xk)⊥⋅(xj−xk)⊥2AT+(fk−fi)(xj−xi)⊥⋅(xj−xk)⊥2AT \int_{\partial{A_i\cap T}}∇f⋅nds =(f_j-f_i)\frac{(x_i-x_k)^\bot \cdot (x_j-x_k)^\bot}{2A_T}+(f_k-f_i)\frac{(x_j-x_i)^\bot \cdot (x_j-x_k)^\bot}{2A_T}

因为
AT=12sinγj∣∣xj−xi∣∣∣∣xj−xk∣∣=12sinγk∣∣xi−xk∣∣∣∣xj−xk∣∣ A_T=\frac{1}{2}sin\gamma_j||x_j-x_i|| ||x_j-x_k||=\frac{1}{2}sin\gamma_k||x_i-x_k|| ||x_j-x_k||

cosγj=(xj−xk)⋅(xj−xk)∣∣xj−xi∣∣∣∣xj−xk∣∣cosγk=(xi−xk)⋅(xj−xk)∣∣xi−xk∣∣∣∣xj−xk∣∣ cos\gamma_j=\frac{(x_j-x_k)\cdot (x_j-x_k)}{||x_j-x_i||||x_j-x_k||} \\ cos\gamma_k=\frac{(x_i-x_k)\cdot (x_j-x_k)}{||x_i-x_k||||x_j-x_k||}

于是得到：
∫∂Ai∩T∇f⋅nds=12(cotγk(fj−fi)+cotγj(fk−fi)) \int_{\partial{A_i\cap T}}∇f⋅nds = \frac{1}{2}(cot\gamma_k(f_j-f_i)+cot\gamma_j(f_k-f_i))
在整个局部平均区域AiA_i上积分得到：
∫∂Ai∩T∇f⋅nds=12∑j∈Ω(i)(cotαij+cotβij)(fj−fi) \int_{\partial{A_i\cap T}}∇f⋅nds = \frac{1}{2} \sum_{j \in \Omega(i)}(cot\alpha_{ij}+cot\beta_{ij})(f_j-f_i)
所以在顶点viv_i处，函数ff的拉普拉斯算子的离散平均为：
Δf(vi)=12Ai∑j∈Ω(i)(cotαij+cotβij)(fj−fi) \Delta f(v_i)=\frac{1}{2A_i} \sum_{j \in \Omega(i)}(cot\alpha_{ij}+cot\beta_{ij})(f_j-f_i)

矩阵化求解拉普拉斯

Δf(vi)=wi∑vj∈N1(vi)wij(f(vj)−f(vi)) \Delta f(v_i) = w_i \sum_{v_j \in \mathcal{N}_1(v_i)} w_{ij} (f(v_j) - f(v_i))
其中wi=12Ai,wij=cotαij+cotβijw_i=\frac{1}{2A_i},w_{ij}=cot\alpha_{ij}+cot\beta_{ij}。将线性和写成矩阵形式：
(Δf(v1)⋮Δf(vn))=DM⏟L(f(v1)⋮f(vn)) \begin{pmatrix} \Delta f(v_1) \\ \vdots \\ \Delta f(v_n) \end{pmatrix}=\underbrace{DM}_L \begin{pmatrix} f(v_1) \\ \vdots \\ f(v_n) \end{pmatrix}
其中D=diag(w1,...,wn)D=diag(w_1,...,w_n)，M是对称矩阵：
mi,j={−∑vk∈N1(vi)wik,i=jwij,vj∈Ni(vi)0,otherwise m_{i,j} = \begin{cases} -\sum_{v_k \in \mathcal{N}_1(v_i)} w_{ik}, &i = j \\ w_{ij}, &v_j \in \mathcal{N}_i(v_i) \\ 0, & \text{otherwise} \end{cases}
更高阶的拉普拉斯可以递归得到：
Δkf(vi)=wi∑vj∈N1(vi)wij(Δk−1f(vj)−Δk−1f(vi)) \Delta^k f(v_i) = w_i \sum_{v_j \in \mathcal{N}_1(v_i)} w_{ij} ( \Delta^{k-1} f(v_j) - \Delta^{k-1} f(v_i))
那么这个矩阵表示就很简单的对应拉普拉斯矩阵LL的kk次幂Lk=(DM)kL^k=(DM)^k。

因此在n个顶点的mesh上更高阶的偏微分方程Δkf=bΔ^kf=b离散化之后会得到一个(n×n)(n×n)的线性系统：
Lkx=b L^kx=b
其中x=(f(v1),...,f(vn))Tx=(f(v_1),...,f(v_n))^T，且b=((b(v1),...,b(vn)))Tb=((b(v_1),...,b(v_n)))^T

矩阵性质

稀疏性 简单说一下吧。拉普拉斯矩阵中只有对角元和边对应的元素是非0的，其余全是0。根据欧拉公式，平均每一行大约7个非0元素。

对称性 由于对角矩阵DD捣乱，所以矩阵L=DML=DM不是对称的。不过对于任意k阶的拉普拉斯系统很容易转换为对称系统，即
M(DM)k−1x=D−1b M(DM)^{k-1}x=D^{-1}b
正定性 一般来讲，解偏微分方程都需要边界条件。即对于一个线性系统Ax=bAx=b，有部分变量xix_i保持不变，此时它们不再是未知量了。此时将对应的列aia_i移到右手边，即b←b−xiaib←b−x_ia_i，同时对应的行从整个系统中移除。注意到此时整个系统还是保持对称的！！！在移除边界条件之后，可以证明矩阵LL是负定的！负定很简单，直接在每个LL上乘一个-1就行，因此我们最后得到
(−1)kM(DM)k−1x=(−1)kD−1b (-1)^kM(DM)^{k-1}x=(-1)^kD^{-1}b
由此我们就得到了一个稀疏、对称且正定的线性系统。

Reference

[1]http://staff.ustc.edu.cn/~fuxm/course/2019_Spring_DGP/index.html
[2]https://optsolution.github.io/archives/19907.html

相关推荐