6. Deformation(变形)
引言
所谓变形,是指交互式地改变网格。
形式化定义:
S
′
=
{
p
+
d
(
p
)
∣
p
∈
S
}
\mathcal{S}' =\{p+d(p)|p\in \mathcal{S}\}
S′={p+d(p)∣p∈S}
对离散三角网格而言,位移函数
d
d
d是分段线性的,完全由原始网格的顶点对应的位移向量定义。
- 用户预先定义一个位移区域, handle points,同时约束一些区域保持不变, f i x e d p o i n t s fixed points fixedpoints
d ( p i ) = d i ˉ , p i ∈ H d ( p i ) = 0 , p i ∈ F d(p_i)=\bar{d_i}, p_i \in \mathcal{H}\\ d(p_i)=0, p_i \in \mathcal{F} d(pi)=diˉ,pi∈Hd(pi)=0,pi∈F
- 剩余无约束的点所在区域 R \mathcal{R} R,需要计算决定
两类变形方法:
- Surface-based deformation,位移函数依赖原始表面 S \mathcal{S} S,并且在三角网格上计算。这类方法控制程度高,每个顶点都可以独立约束。相应的,鲁棒性和计算效率受网格复杂度和三角网格质量的影响。
- Space deformation,隐式变形,计算整个嵌入空间的位移函数,然后warp。不需要在三角网格上计算,受网格质量和复杂度影响小。相应的,变形的控制程度就低了。
上述方法都是线性变形,需要求解一个线性方程系统。这些方法的好处在于求解高效;相应的,对大尺度变形,效果可能违反直觉。非线性变形可以解决这个问题,相应数值计算复杂。
Transformation Propagation
最简单且常见的变形方法就是transformation propagation,通过一个标量场 s : S → [ 0 , 1 ] s: \mathcal{S}\rightarrow[0,1] s:S→[0,1],1表示完全变形(handle region),0表示无变形(fixed region),中间区域可以sommthly blending
怎样构建这个标量场呢?
-
一种方式就是计算变形点 p p p到fixed region和handle region的距离,然后插值。这里的距离可以是geodesic distance on the surface,也可以是Eucidean distance on the space。前者效果好但计算复杂。
-
另一种方法是使用调和场计算,手工区域和固定区域有Dirichlet约束 s ( p i ) = 1 , p i ∈ H , s ( p i ) = 1 , p i ∈ F s(p_i)=1, p_i \in \mathcal{H}, s(p_i)=1, p_i \in \mathcal{F} s(pi)=1,pi∈H,s(pi)=1,pi∈F,变形区域满足 Δ s ( p i ) = 0 \Delta s(p_i)=0 Δs(pi)=0, 然后求解 s ( p i ) s(p_i) s(pi)(求解线性Laplacian系统)。这种方法计算更复杂,但效果平滑。这种方法还可以通过一个transfer function t : [ 0 , 1 ] → [ 0 , 1 ] t:[0,1]\rightarrow [0,1] t:[0,1]→[0,1],进一步控制blending过程。
计算出标量场后,对每个顶点插值即可。
p
i
′
=
s
(
o
i
)
T
(
p
i
)
+
(
1
−
s
(
p
i
)
)
p
i
p_i' = s(o_i)T(p_i) + (1-s(p_i))p_i
pi′=s(oi)T(pi)+(1−s(pi))pi
这里
T
T
T是用户定义的变形操作。
这类方法的问题是基于距离传播,可能出现违反几何直觉的方案,例如中间塌陷
Shell-based Deformation
受物理启发的,最小化能量的方法。表面看做物理的,拉伸力和张力作用在上面。
设 S \mathcal{S} S和 S ′ \mathcal{S}' S′都是光滑参数化表面。位移函数也是 d : Ω → R 3 d:\Omega\rightarrow \mathbb{R}^3 d:Ω→R3
First fundamental form 和 second fundamental form
I
(
u
,
v
)
,
I
I
(
u
,
v
)
I(u,v), II(u,v)
I(u,v),II(u,v)可以用来度量几何本证性质(与参数化无关)。当表面形变时,基本型的变化可以看组总是一个elastic thin shell energy, 表征了stretching和bending:
E
(
S
′
)
=
∬
Ω
k
s
∣
∣
I
′
(
u
,
v
)
−
I
(
u
,
v
)
∣
∣
F
2
+
k
b
∣
∣
I
I
′
(
u
,
v
)
−
I
I
(
u
,
v
)
∣
∣
F
2
E(S') = \iint_{\Omega}k_s ||I'(u,v)-I(u,v)||_F^2 + k_b ||II'(u,v)-II(u,v)||_F^2
E(S′)=∬Ωks∣∣I′(u,v)−I(u,v)∣∣F2+kb∣∣II′(u,v)−II(u,v)∣∣F2
需要最小化这个elastic energy,在给定用户defomation constraints时。
上式计算复杂,把基本型的差简化为位移函数
d
d
d的偏导,定义了thin-shell energy:
E
(
d
)
=
∬
Ω
k
s
(
∣
∣
d
u
(
u
,
v
)
∣
∣
2
+
∣
∣
d
v
(
u
,
v
)
∣
∣
2
)
+
k
b
(
∣
∣
d
u
u
(
u
,
v
)
2
∣
∣
+
2
∣
∣
d
u
v
(
u
,
v
)
∣
∣
2
+
∣
∣
d
v
v
(
u
,
v
)
∣
∣
2
)
d
u
d
v
E(d)=\iint_{\Omega}k_s(||d_u(u,v)||^2 + ||d_v(u,v)||^2) + k_b(||d_{uu}(u,v)^2||+2||d_{uv}(u,v)||^2 + ||d_{vv}(u,v)||^2)dudv
E(d)=∬Ωks(∣∣du(u,v)∣∣2+∣∣dv(u,v)∣∣2)+kb(∣∣duu(u,v)2∣∣+2∣∣duv(u,v)∣∣2+∣∣dvv(u,v)∣∣2)dudv
再使用变分法构造Euler-Langrange方程:
−
k
s
Δ
d
+
k
b
Δ
2
d
=
0
-k_s\Delta d + k_b \Delta^2 d = 0
−ksΔd+kbΔ2d=0
上式是一个PDE。然后,我们就可以将其离散化(例如,使用discrete cotangent Laplacian,以及二阶差分bi-Laplacian)。这样对表面的每个顶点都能给出条件方程。然后求解线性方程系统就可以了(未知元就是自由顶点的个数)。
上述优化方法得到的变形结果是 C 1 C^1 C1连续的。
基于shell能量的方法的最大问题还是计算复杂度
多尺度变形法
基于Shell的方法是线性的,是对非线性shell energy的化简。在处理fine-scale surface details时,变形效果会违反几何直觉。因为几何细节的局部旋转是一种非线性行为,不可能通过纯线性方法实现。
一种更好的保留几何细节,同时仍然采用线性变形的方法,叫多尺度方法。
这种方法的核心思想是把物体分解为两部分,依照频率带,分成低频部分(平滑全局形状)和高频部分(精细),对低频部分做变形,同时保留高频细节。
一个典型的多尺度方法包括三部分:decomposition operator,editing operator,和reconstruction operator。这里只是增加了一个defomation operator,把基底表面 B \mathcal{B} B变形为 B ′ \mathcal{B}' B′。
同时,分解和重建可以是层级的。
- Decomposition可以用mesh smoothing或fairing的方法
- Deformation可以用前述deformation的方法
- Geometric detail representation比较麻烦
Displacement Vectors
在原始网格上,计算表面 S \mathcal{S} S上每一个顶点到基底表面 B \mathcal{B} B的位移向量 h i = ( p i − b i ) h_i=(p_i-b_i) hi=(pi−bi)
因此有:
p
i
=
b
i
+
h
i
,
h
i
∈
R
3
p_i=b_i+h_i, h_i \in \mathbb{R}^3
pi=bi+hi,hi∈R3
h
i
h_i
hi使用局部坐标系(local frame)编码,具体地,就是基底表面
B
\mathcal{B}
B上对应的顶点
b
i
b_i
bi的法向量
n
i
n_i
ni和两个切向量
t
i
,
1
,
t
i
,
2
t_{i,1}, t_{i,2}
ti,1,ti,2:
h
i
=
α
i
n
i
+
β
i
t
i
,
1
+
γ
i
t
i
,
2
h_i=\alpha_i n_i + \beta_i t_{i,1} + \gamma_i t_{i,2}
hi=αini+βiti,1+γiti,2
这里
α
i
,
β
i
,
γ
i
\alpha_i, \beta_i, \gamma_i
αi,βi,γi就是局部坐标系的坐标。当基底表面变形后,假设位移向量随着基地表面局部坐标系的旋转而旋转,那么就有:
p
i
′
=
b
i
′
+
α
i
n
i
′
+
β
i
t
i
,
1
′
+
γ
i
t
i
,
2
′
p_i' = b_i' + \alpha_i n_i' + \beta_i t_{i,1}' + \gamma_i t_{i,2}'
pi′=bi′+αini′+βiti,1′+γiti,2′
这里法向好定义,切向量不太好定义。一般是把点
b
i
b_i
bi的一条相邻边投影到切平面上,作为一个切向量,然后叉乘得到另一个切向量。
Normal Displacement
长距离的位移向量可能导致不稳定。因此,位移向量应当尽可能短。
当位移向量连结的不再是点
p
p
p和基底表面对应点
b
b
b,而是它在基底表面上的最近邻点(cloest)时,位移向量的长度确保最短。基于此,诞生了normal displacements(法向偏移):
p
i
=
b
i
+
h
i
⋅
n
i
,
h
i
∈
R
p_i = b_i + h_i \cdot n_i, h_i \in \mathbb{R}
pi=bi+hi⋅ni,hi∈R
因为之前定义的
h
i
h_i
hi并不是和表面法向平行,因此,需要re-sampling原始表面和基底表面,构建点-点关系。
这种re-sampling既既可以从基底表面上的顶点 b i b_i bi出发,寻找在原始表面上的对应顶点 p i p_i pi;也可以反过来。但是原始表面一般比基底表面细节更多,前者可能引入走样。
对后者方法,使用局部牛顿迭代法找到 p i p_i pi在基底表面上对应的点 b i b_i bi。这里, b i b_i bi不一定非要是基底表面的顶点,任意点都可以。使用重心坐标插值即可。插值点的法向也可以使用重心插值计算得到。(下面的重心坐标就是讲的这个,法向可以看成是定义在顶点上的函数)。
这第二种方法避免了 S \mathcal{S} S的resampling,保留了它的细节。
高级技术
基于向量位移的方法没有考虑邻居位移向量,在处理一些角度变化较大的形变时,会出现一些几何失真,甚至出现表面自交。
Displacement volumes
使用位移体取代位移向量。原始表面上的三个点和对应基底表面的三个点构成一个三棱柱。这个三棱柱的体积用来表示细节系数,在变形时保持不变。代价是高昂的计算开销。
Deformation transfer
把 B → B ′ \mathcal{B}\rightarrow \mathcal{B}' B→B′的transformation给迁移到 S → S ′ \mathcal{S} \rightarrow \mathcal{S}' S→S′上。这个做法在其他应用商店很常见。效果和时间复杂度介于displacement vector和displacement volume之间。
微分坐标法
几何模型复杂时,多尺度方法的分解也变得麻烦起来。为了避免显式的多尺度分解,另一类方法通过修改微分表面的性质取代空间坐标,然后重建变形后的表面。
基于梯度的变形
基于梯度的方法分两步走:第一步,首先对原始表面的梯度进行操作(变形);第二步,找到和目标梯度场匹配的表面(使用最小二乘优化)。
第一步
如何对梯度进行操作呢?
对每个面片,梯度是一个3*3的Jacobian矩阵 J T J_T JT.给它左乘变换矩阵 M T M_T MT,就可以得到新的梯度 J T ′ J_T' JT′.
怎么从用户交互得到变换矩阵 M T M_T MT呢?见后文。
通过对每个面片的梯度进行变换,破坏了原始mesh的连结。
第二步
下一步就是要找到新的顶点
p
i
′
p_i'
pi′,使得变形后的网格的梯度
∇
p
′
∣
T
\nabla p'|_T
∇p′∣T尽可能和目标梯度
J
T
′
J_T'
JT′接近。这是一个能量最小化问题:
E
(
f
)
=
∬
Ω
∣
∣
∇
f
(
u
,
v
)
−
g
(
u
,
v
)
∣
∣
2
d
u
d
v
E(f) = \iint_{\Omega} ||\nabla f(u,v) - g(u,v)||^2 dudv
E(f)=∬Ω∣∣∇f(u,v)−g(u,v)∣∣2dudv
通过变分法构造Euler-Lagrange方程
Δ
f
=
d
i
v
g
\Delta f = div \bf{g}
Δf=divg
相较前面基于Shell能量的方法,这种方法只需要计算一次Poisson系统,比计算bi-Laplacian高效一点。不过代价是,这种方法只能保证变形区域边界的
C
0
C^0
C0连续,而Shell-based可以保证
C
1
C^1
C1连续。
基于Laplacian的变形
Laplacian editing是另一种基于微分坐标的方法。和基于梯度的方法类似,只不过操作的不是per-face的梯度,而是per-vertex的Laplacians。
第一步
计算每个顶点的初始Laplace坐标 δ i = Δ ( p i ) \delta_i = \Delta(p_i) δi=Δ(pi),变形为 δ i ′ = M i δ i \delta_i'=M_i\delta_i δi′=Miδi。
第二步
找到新的坐标 p i ′ p_i' pi′来匹配目标Laplacian坐标。
这个问题也是一个能量最小化问题:
E
(
p
′
)
=
∬
Ω
∣
∣
Δ
p
′
(
u
,
v
)
−
δ
′
(
u
,
v
)
∣
∣
2
d
u
d
v
E(p') = \iint_{\Omega}||\Delta p'(u,v) - \delta'(u,v)||^2 dudv
E(p′)=∬Ω∣∣Δp′(u,v)−δ′(u,v)∣∣2dudv
变为Euler-Lagrange方程
Δ
2
p
′
=
Δ
δ
′
\Delta^2 p' = \Delta \delta'
Δ2p′=Δδ′
对离散三角网格,这是个bi-Laplacian系统。这里使用contangent Laplacian效果比uniform Laplacian好一点。
Laplacian方法和Shell-based方法有一定关系。Laplacian方法也是 C 1 C^1 C1连续的。
局部变形(计算 M i M_i Mi或 M T M_T MT)
计算 M T M_T MT
用户在handle region上的操作定义为一个放射变换
T
(
x
)
=
M
x
+
t
T(x) = Mx + t
T(x)=Mx+t
T
(
x
)
T(x)
T(x)的梯度
M
M
M表示了旋转、尺度/错切等变换。
怎么传播(propagate)呢?类似尺度场的做法,从handle region的变换 M M M衰减到固定区域的 I d Id Id(identity)。
但需要注意的是,旋转的插值和尺度的插值不一样,因此要先把这个矩阵给分解开来。把一个矩阵分解为旋转和尺度/错切的方法叫polar decomposition(极分解)。具体地,对
M
M
M做SVD分解
M
=
U
Σ
V
T
M=U\Sigma V^T
M=UΣVT,然后
R
=
U
V
T
S
=
V
Σ
V
T
R = UV^T\\ S = V\Sigma V^T
R=UVTS=VΣVT
现在可以分别对旋转和尺度插值了。
这种方法对旋转表现很好,但是对平移很无力。当平移变化大的时候,变形结果会反直觉。
Implicit optimization (计算 M i M_i Mi)
联合优化问题
E
(
p
′
)
=
∑
i
=
1
n
A
i
∣
∣
M
i
δ
i
−
Δ
p
i
′
∣
∣
2
E(p') = \sum_{i=1}^n A_i ||M_i\delta_i - \Delta p_i'||^2
E(p′)=i=1∑nAi∣∣Miδi−Δpi′∣∣2
这里
A
i
A_i
Ai是local vertex area。
M
i
M_i
Mi依赖于新顶点位置
p
j
′
p_j'
pj′.
为了防止非线性优化,局部变换被限制为线性相似变换,也就是说变换矩阵是个反对称矩阵。
Freeform Deformation (自由变形)
Space deformation, deform the ambient space and thus implicitly deform the embedded objects.
Lattice-based FFD
传统FFD把空间变形表示为三变量的张量积样条函数:
d
(
u
,
v
,
w
)
=
∑
i
∑
j
∑
k
δ
c
i
j
k
N
i
(
u
)
N
j
(
v
)
N
k
(
w
)
d(u,v,w) = \sum_i\sum_j\sum_k \delta c_{ijk}N_i(u)N_j(v)N_k(w)
d(u,v,w)=i∑j∑k∑δcijkNi(u)Nj(v)Nk(w)
N
i
N_i
Ni是B-样条基函数,
δ
c
i
j
k
\delta c_{ijk}
δcijk是控制点
c
i
j
k
c_{ijk}
cijk的位移。
可以简写为
d
(
u
)
=
∑
l
=
1
n
δ
c
l
N
l
(
u
)
d(u) = \sum_{l=1}^n \delta c_l N_l(u)
d(u)=l=1∑nδclNl(u)
对原始表面的每个点
p
i
p_i
pi,只要知道它的参数值
u
i
=
(
u
i
,
v
i
,
w
i
)
u_i=(u_i,v_i,w_i)
ui=(ui,vi,wi),满足
p
i
=
∑
l
c
l
N
l
(
u
i
)
p_i=\sum_l c_lN_l(u_i)
pi=∑lclNl(ui)。顶点就被变换到
p
i
′
=
p
i
+
d
(
u
i
)
p_i'=p_i+d(u_i)
pi′=pi+d(ui)。这个计算相当高效,因为
N
l
(
u
i
)
N_l(u_i)
Nl(ui)是预计算且不变的。
从上式还可以知道,变形是通过控制点的位移 δ c l \delta c_l δcl描述的。这种通过控制点的方式有时候显得很无力;支持的变形也很难预测。
通过用户界面,让用户指定表面点的位移而不是控制点的位移可以简化变形处理。但也有不足。
此外,这类方法依赖规则grid的基函数。如果变形不能很好地和grid axes对齐,效果也不好。
Cage-based FFD
- A generalization of the Lattice-based FFD
- Cage is a coarse, arbitrary triangle mesh enclosing the object to be modified.
径向基函数方法
在surface-based方法里,高质量的结果是通过一个变换函数 d : S → R 3 d:\mathcal{S} \rightarrow \mathbb{R}^3 d:S→R3,插值用户的位移约束实现的。一般是最小化一些能量。类似的,space deformation也可以通过一个光滑插值三变量函数 d : R 3 → R 3 d: \mathbb{R}^3 \rightarrow \mathbb{R}^3 d:R3→R3,最小化类似的fairness能量。
一个三变量径向基变换可以使用一个superposition 径向对称核
ϕ
i
(
x
)
\phi_i(x)
ϕi(x)定义:
d
(
x
)
=
∑
j
=
1
n
w
j
ϕ
(
∣
∣
c
j
−
x
∣
∣
)
+
π
(
x
)
d(x) = \sum_{j=1}^n w_j \phi(||c_j-x||)+\pi(x)
d(x)=j=1∑nwjϕ(∣∣cj−x∣∣)+π(x)
ϕ
j
(
x
)
=
ϕ
(
∣
∣
c
j
−
x
∣
∣
)
\phi_j(x) = \phi(||c_j-x||)
ϕj(x)=ϕ(∣∣cj−x∣∣)是基函数,对应第j个中心
c
j
∈
R
3
c_j \in \mathbb{R}^3
cj∈R3.
怎么找到RBF函数,能够在位移约束 d ( p i ) = d ˉ i d(p_i)=\bar{d}_i d(pi)=dˉi之间插值出连续光滑的结果?
求解线性方程系统。
注意到kernel function的选择对计算复杂度和结果的光顺有很大影响。紧致支撑的径向基函数导致一个稀疏的线性系统,可以用来插值大规模的约束。但是不能提供和具有global support的RBF 一样的光顺程度。Global supported 基函数和surface-based方法提供同样程度的光顺。不过global support导致一个致密的线性系统,数值上求解比较困难。
线性方法的局限性
- Shell-based 方法对平移比较好,细节保持较好;但对大的旋转不好;
- Gradient-based editing对旋转较好,但对平移不行;
- Laplacian-based editing 平移和旋转都挺好。只不过对大的变换可能有失真。
Barycentric coordinates(重心坐标)
- 三角形的重心坐标(面积比)
- 推广到四面体(体积比)
- 推广到任意简单形
Applications
- Function interpolation
- Function composite
- …
Generalized barycentric coordinates (GBC)
-
Convex polygon
-
Conditions
- ϕ i ( x ) ≥ 0 \phi_i(x) \geq 0 ϕi(x)≥0
- ∑ i = 1 n ϕ i ( x ) = 1 , ∑ i = 1 n ϕ i ( x ) v i = x \sum_{i=1}^n \phi_i(x)=1, \sum_{i=1}^n \phi_{i}(x)v_i=x ∑i=1nϕi(x)=1,∑i=1nϕi(x)vi=x
-
Properties
- Lagrange property ϕ i ( v j ) = δ i j \phi_i(v_j)=\delta_{ij} ϕi(vj)=δij
- Piecewise linearity on ∂ P \partial P ∂P
- Interpolation
- Linear preision
Bilinear interpolation
在矩形或正方形上定义的特殊的重心坐标
Convex quadrilaterals
- View P P P as the image of a bilinear map from the unit square [ 0 , 1 ] × [ 0 , 1 ] [0,1]\times[0,1] [0,1]×[0,1]
- For each x ∈ P x\in P x∈P, there exist unique λ , μ ∈ ( 0 , 1 ) \lambda, \mu \in (0,1) λ,μ∈(0,1),满足类似重心坐标的关系。这时, x x x的重心坐标为:
ϕ 1 ( x ) = ( 1 − λ ) ( 1 − μ ) ϕ 2 ( x ) = λ ( 1 − μ ) ϕ 3 ( x ) = λ μ ϕ 4 ( x ) = ( 1 − λ ) μ \phi_1(x)=(1-\lambda)(1-\mu)\\ \phi_2(x)=\lambda(1-\mu)\\ \phi_3(x)=\lambda\mu\\ \phi_4(x)=(1-\lambda)\mu ϕ1(x)=(1−λ)(1−μ)ϕ2(x)=λ(1−μ)ϕ3(x)=λμϕ4(x)=(1−λ)μ
How to solve?
Inverse of the bilinear map
构造性方法证明