Laplacian on Triangle Mesh
局部平均区域
基本思想是mesh面上的点x的值设置为其局部领域的平均。下图是几种具体的情况:
其中蓝色区域为局部平均区域。
- Barycentric cell:三角形的重心和三角形边的中点相连
- Voronoi cell:三角形的外心和三角形的中点相连(对于钝角三角形可能会出现外心在外部)
- Mixed Voronoi cell:在Voronoi cell的基础上,将钝角三角形的外心替换成钝角对边的中点
梯度
拉普拉斯算子是梯度的散度,所以我们需要先求出梯度。
属性插值函数:
f
(
x
)
=
α
f
i
+
β
f
j
+
γ
f
k
f(x)=\alpha f_i+\beta f_j+\gamma f_k \\
f(x)=αfi+βfj+γfk
其中:
f
i
,
f
j
,
f
k
f_i,f_j,f_k
fi,fj,fk是三角形三个顶点的属性值,x为三角形内部的点,
α
,
β
,
γ
\alpha , \beta , \gamma
α,β,γ 是点x在三角形的重心坐标,
α
+
β
+
γ
=
1
\alpha + \beta + \gamma = 1
α+β+γ=1。
梯度:
∇
x
f
(
x
)
=
f
i
∇
x
α
+
f
j
∇
x
β
+
f
k
∇
x
γ
\nabla_xf(x)=f_i\nabla_x\alpha +f_j\nabla_x\beta +f_k\nabla_x\gamma
∇xf(x)=fi∇xα+fj∇xβ+fk∇xγ
对于
α
\alpha
α来说:
α
=
A
i
A
T
=
(
(
x
−
x
j
)
⋅
(
x
k
−
x
j
)
⊥
∣
∣
x
k
−
x
j
∣
∣
2
)
∣
∣
x
k
−
x
j
∣
∣
2
2
A
T
=
(
x
−
x
j
)
⋅
(
x
k
−
x
j
)
⊥
2
A
T
\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}
α=ATAi=2AT((x−xj)⋅∣∣xk−xj∣∣2(xk−xj)⊥)∣∣xk−xj∣∣2=2AT(x−xj)⋅(xk−xj)⊥
对
x
x
x求导:
∇
x
α
=
(
x
k
−
x
j
)
⊥
2
A
T
∇
x
β
=
(
x
i
−
x
k
)
⊥
2
A
T
∇
x
γ
=
(
x
j
−
x
i
)
⊥
2
A
T
\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}
∇xα=2AT(xk−xj)⊥∇xβ=2AT(xi−xk)⊥∇xγ=2AT(xj−xi)⊥
得到梯度公式:
∇
x
f
(
x
)
=
f
i
(
x
k
−
x
j
)
⊥
2
A
T
+
f
j
(
x
i
−
x
k
)
⊥
2
A
T
+
f
k
(
x
j
−
x
i
)
⊥
2
A
T
\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}
∇xf(x)=fi2AT(xk−xj)⊥+fj2AT(xi−xk)⊥+fk2AT(xj−xi)⊥
因为
(
x
k
−
x
j
)
⊥
+
(
x
i
−
x
k
)
⊥
+
(
x
j
−
x
i
)
⊥
=
0
(x_k-x_j)^\bot+(x_i-x_k)^\bot+(x_j-x_i)^\bot=0
(xk−xj)⊥+(xi−xk)⊥+(xj−xi)⊥=0
所以
∇
x
f
(
x
)
=
(
f
j
−
f
i
)
(
x
i
−
x
k
)
⊥
2
A
T
+
(
f
k
−
f
i
)
(
x
j
−
x
i
)
⊥
2
A
T
\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}
∇xf(x)=(fj−fi)2AT(xi−xk)⊥+(fk−fi)2AT(xj−xi)⊥
离散拉普拉斯算子
均匀拉普拉斯
Δ f ( v i ) = 1 ∣ N 1 ( v i ) ∣ ∑ v j ∈ N 1 ( v i ) ( f j − f i ) \Delta f(v_i) =\frac {1}{|N_1(v_i)|}\sum_{v_j \in N_1(v_i)}(f_j-f_i) Δf(vi)=∣N1(vi)∣1vj∈N1(vi)∑(fj−fi)
此公式可以用于各向同性的remesh上。
余切形式的拉普拉斯
使用混合的有限元/有限体方法(mixed finite element/finite volume method),拉普拉斯算子可以更加精确地进行离散推导。目标是在局部平均区域A_i=A(v_i)上,对逐片线性函数梯度地散度进行积分。
∫
A
i
Δ
f
d
A
=
∫
A
i
d
i
v
∇
f
d
A
=
∫
∂
A
i
(
∇
f
)
⋅
n
d
s
\int_{A_i}\Delta fdA=\int_{A_i}div\nabla fdA=\int_{\partial{A_i}}(\nabla f)\cdot n ds
∫AiΔfdA=∫Aidiv∇fdA=∫∂Ai(∇f)⋅nds
将积分分解到每个三角形上,并且局部的Voronoi区域通过三角形两条边的中点
a
a
a和
b
b
b,并且每个三角形中
∇
f
(
x
)
\nabla f(x)
∇f(x)是常数,那么在该三角形
T
T
T上的积分为
∫
∂
A
i
∩
T
∇
f
⋅
n
d
s
=
∇
f
⋅
(
a
−
b
)
⊥
=
1
2
∇
f
(
u
)
⋅
(
x
j
−
x
k
)
⊥
\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=∇f⋅(a−b)⊥=21∇f(u)⋅(xj−xk)⊥
将梯度公式代入:
∫
∂
A
i
∩
T
∇
f
⋅
n
d
s
=
(
f
j
−
f
i
)
(
x
i
−
x
k
)
⊥
⋅
(
x
j
−
x
k
)
⊥
2
A
T
+
(
f
k
−
f
i
)
(
x
j
−
x
i
)
⊥
⋅
(
x
j
−
x
k
)
⊥
2
A
T
\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}
∫∂Ai∩T∇f⋅nds=(fj−fi)2AT(xi−xk)⊥⋅(xj−xk)⊥+(fk−fi)2AT(xj−xi)⊥⋅(xj−xk)⊥
因为
A
T
=
1
2
s
i
n
γ
j
∣
∣
x
j
−
x
i
∣
∣
∣
∣
x
j
−
x
k
∣
∣
=
1
2
s
i
n
γ
k
∣
∣
x
i
−
x
k
∣
∣
∣
∣
x
j
−
x
k
∣
∣
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||
AT=21sinγj∣∣xj−xi∣∣∣∣xj−xk∣∣=21sinγk∣∣xi−xk∣∣∣∣xj−xk∣∣
c o s γ j = ( x j − x k ) ⋅ ( x j − x k ) ∣ ∣ x j − x i ∣ ∣ ∣ ∣ x j − x k ∣ ∣ c o s γ k = ( x i − x k ) ⋅ ( x j − x k ) ∣ ∣ x i − x k ∣ ∣ ∣ ∣ x j − x k ∣ ∣ 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||} cosγj=∣∣xj−xi∣∣∣∣xj−xk∣∣(xj−xk)⋅(xj−xk)cosγk=∣∣xi−xk∣∣∣∣xj−xk∣∣(xi−xk)⋅(xj−xk)
于是得到:
∫
∂
A
i
∩
T
∇
f
⋅
n
d
s
=
1
2
(
c
o
t
γ
k
(
f
j
−
f
i
)
+
c
o
t
γ
j
(
f
k
−
f
i
)
)
\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))
∫∂Ai∩T∇f⋅nds=21(cotγk(fj−fi)+cotγj(fk−fi))
在整个局部平均区域
A
i
A_i
Ai上积分得到:
∫
∂
A
i
∩
T
∇
f
⋅
n
d
s
=
1
2
∑
j
∈
Ω
(
i
)
(
c
o
t
α
i
j
+
c
o
t
β
i
j
)
(
f
j
−
f
i
)
\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)
∫∂Ai∩T∇f⋅nds=21j∈Ω(i)∑(cotαij+cotβij)(fj−fi)
所以在顶点
v
i
v_i
vi处,函数
f
f
f的拉普拉斯算子的离散平均为:
Δ
f
(
v
i
)
=
1
2
A
i
∑
j
∈
Ω
(
i
)
(
c
o
t
α
i
j
+
c
o
t
β
i
j
)
(
f
j
−
f
i
)
\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)=2Ai1j∈Ω(i)∑(cotαij+cotβij)(fj−fi)
矩阵化求解拉普拉斯
Δ
f
(
v
i
)
=
w
i
∑
v
j
∈
N
1
(
v
i
)
w
i
j
(
f
(
v
j
)
−
f
(
v
i
)
)
\Delta f(v_i) = w_i \sum_{v_j \in \mathcal{N}_1(v_i)} w_{ij} (f(v_j) - f(v_i))
Δf(vi)=wivj∈N1(vi)∑wij(f(vj)−f(vi))
其中
w
i
=
1
2
A
i
,
w
i
j
=
c
o
t
α
i
j
+
c
o
t
β
i
j
w_i=\frac{1}{2A_i},w_{ij}=cot\alpha_{ij}+cot\beta_{ij}
wi=2Ai1,wij=cotαij+cotβij。将线性和写成矩阵形式:
(
Δ
f
(
v
1
)
⋮
Δ
f
(
v
n
)
)
=
D
M
⏟
L
(
f
(
v
1
)
⋮
f
(
v
n
)
)
\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}
⎝⎜⎛Δf(v1)⋮Δf(vn)⎠⎟⎞=L
DM⎝⎜⎛f(v1)⋮f(vn)⎠⎟⎞
其中
D
=
d
i
a
g
(
w
1
,
.
.
.
,
w
n
)
D=diag(w_1,...,w_n)
D=diag(w1,...,wn),M是对称矩阵:
m
i
,
j
=
{
−
∑
v
k
∈
N
1
(
v
i
)
w
i
k
,
i
=
j
w
i
j
,
v
j
∈
N
i
(
v
i
)
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}
mi,j=⎩⎪⎨⎪⎧−∑vk∈N1(vi)wik,wij,0,i=jvj∈Ni(vi)otherwise
更高阶的拉普拉斯可以递归得到:
Δ
k
f
(
v
i
)
=
w
i
∑
v
j
∈
N
1
(
v
i
)
w
i
j
(
Δ
k
−
1
f
(
v
j
)
−
Δ
k
−
1
f
(
v
i
)
)
\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))
Δkf(vi)=wivj∈N1(vi)∑wij(Δk−1f(vj)−Δk−1f(vi))
那么这个矩阵表示就很简单的对应拉普拉斯矩阵
L
L
L的
k
k
k次幂
L
k
=
(
D
M
)
k
L^k=(DM)^k
Lk=(DM)k。
因此在n个顶点的mesh上更高阶的偏微分方程
Δ
k
f
=
b
Δ^kf=b
Δkf=b离散化之后会得到一个
(
n
×
n
)
(n×n)
(n×n)的线性系统:
L
k
x
=
b
L^kx=b
Lkx=b
其中
x
=
(
f
(
v
1
)
,
.
.
.
,
f
(
v
n
)
)
T
x=(f(v_1),...,f(v_n))^T
x=(f(v1),...,f(vn))T,且
b
=
(
(
b
(
v
1
)
,
.
.
.
,
b
(
v
n
)
)
)
T
b=((b(v_1),...,b(v_n)))^T
b=((b(v1),...,b(vn)))T
矩阵性质
稀疏性 简单说一下吧。拉普拉斯矩阵中只有对角元和边对应的元素是非0的,其余全是0。根据欧拉公式,平均每一行大约7个非0元素。
对称性 由于对角矩阵
D
D
D捣乱,所以矩阵
L
=
D
M
L=DM
L=DM不是对称的。不过对于任意k阶的拉普拉斯系统很容易转换为对称系统,即
M
(
D
M
)
k
−
1
x
=
D
−
1
b
M(DM)^{k-1}x=D^{-1}b
M(DM)k−1x=D−1b
正定性 一般来讲,解偏微分方程都需要边界条件。即对于一个线性系统
A
x
=
b
Ax=b
Ax=b,有部分变量
x
i
x_i
xi保持不变,此时它们不再是未知量了。此时将对应的列
a
i
a_i
ai移到右手边,即
b
←
b
−
x
i
a
i
b←b−x_ia_i
b←b−xiai,同时对应的行从整个系统中移除。注意到此时整个系统还是保持对称的!!!在移除边界条件之后,可以证明矩阵
L
L
L是负定的!负定很简单,直接在每个
L
L
L上乘一个-1就行,因此我们最后得到
(
−
1
)
k
M
(
D
M
)
k
−
1
x
=
(
−
1
)
k
D
−
1
b
(-1)^kM(DM)^{k-1}x=(-1)^kD^{-1}b
(−1)kM(DM)k−1x=(−1)kD−1b
由此我们就得到了一个稀疏、对称且正定的线性系统。
Reference
[1]http://staff.ustc.edu.cn/~fuxm/course/2019_Spring_DGP/index.html
[2]https://optsolution.github.io/archives/19907.html