Transcript Modelisation volumique de 3

Mesh Parameterization:
Theory and Practice
Making it work in practice - Numerics
Bruno Lévy - INRIA
Overview
1. Numerical Problems
2. Linear and Quadratic
3. Non-linear
Motivations:
Need for scalability in GP
1970’s
2000’s
1. Numerical Problems in GP
Mesh Parameterization
Ui =
i,j j
j  Ni
i
j1
j…
Sa U
j2
[Tutte], [Floater]
 i,j ai,j > 0
ai,i = - S ai,j
The border is mapped to
a convex polygon
1. Numerical Problems in GP
Discrete Fairing
i
j1
j…
F(p)= S
pi - S ai,jpj
i
jN
2
i
j2
[Mallet], [Kobbelt], [Sorkine]…
1. Numerical Problems in GP
Neighborhoods
F = sum of terms, attached to neighborhoods
Discrete fairing
[Kobbelt98, Mallet95]
Parameterization
[Desbrun02]
Deformations
[CohenOr], [Sorkine]
Curv. Estimation
[Cohen-Steiner 03]
Texture mapping
[Levy01]
Discrete fairing
[Desbrun], ...
Parameterization
[Haker00]
[Levy02]
[Eck]
2. Linear and Quadratic GP
Removing degrees of freedom
F(x) = A x - b
2
2
F(xf) = [ Af Al]
xf
xl
-b
2. Linear and Quadratic GP
Removing degrees of freedom
2
2
F(xf) = A.x - d
Aft.Af.xf = Aft.d - AftAl.xl
}
}
F(xf) minimum
= Al.xl + Af.xf - d
M.x
=
b
The problem: (1) construct a linear system
(2) solve a linear system
2. Linear and Quadratic GP
The OpenNL approach
(http://alice.loria.fr/software)
The problem: (1) construct a linear system
(2) solve a linear system
NlLockVariable(i1, val1)
NlLockVariable(i2, val2)
NlSolve()
…
For each stencil instance (one-rings):
NlBeginRow();
NlAddCoefficient(i, a); Need for
• Dynamic Matrix DS
…
• Updating formula
NlEndRow();
2. Linear and Quadratic GP
Direct Solvers (LU)
A Textbook solver: LU factorization (and Cholesky)
M
U
= L
(1) solve
(2) solve
L
Y
U
a ‘small’ problem: O(n3) !!
n=100 0.01 s
n=106 10 centuries
= b
x =Y
L
U
x = b
2. Linear and Quadratic GP
Successive Over-Relaxation (Gauss-Seidel)
xfi
i,i
….
….
….
mn,n xf
nf
ci- S mi,j xfj)
(
m
j=i
1
ci
….
mn,1 … mn,j …
mi,n xfi =
….
….
….
mi,1 … mi,j …
m1,n xf1
….
….
….
m1,1 … m1,j …
used in [Taubin95]
…
2. Linear and Quadratic DGP
Successive Over-Relaxation (Gauss-Seidel)
1000 iterations S.O.R.
2. Linear and Quadratic GP
White Magic: The Conjugate Gradient
inline int solve_conjugate_gradient(
const SparseMatrix &A, const Vector& b, Vector& x,
double eps, int max_iter
){
int N = A.n() ;
double t, tau, sig, rho, gam;
double bnorm2 = BLAS::ddot(N,b,1,b,1) ;
double err=eps*eps*bnorm2 ;
mult(A,x,g);
BLAS::daxpy(N,-1.,b,1,g,1);
BLAS::dscal(N,-1.,g,1);
BLAS::dcopy(N,g,1,r,1);
while ( BLAS::ddot(N,g,1,g,1)>err && its < max_iter) {
mult(A,r,p);
rho=BLAS::ddot(N,p,1,p,1);
sig=BLAS::ddot(N,r,1,p,1);
tau=BLAS::ddot(N,g,1,r,1);
t=tau/sig;
BLAS::daxpy(N,t,r,1,x,1);
BLAS::daxpy(N,-t,p,1,g,1);
gam=(t*t*rho-tau)/tau;
BLAS::dscal(N,gam,r,1);
BLAS::daxpy(N,1.,g,1,r,1);
++its;
}
return its ;
}
Only simple
vector ops (BLAS)
Complicated ops:
Matrix x vector
(see course notes)
[Shewchuck: CG without the agonizing pain]
Iterative Solvers
•Successive Over-Relaxation
•Sparse C.G.
>100x speedup
Precond.
Conj. Grad.
5
- log|G.x + c |
Conj. Grad.
2.5
0
0
S.O.R.
Time(s)
45
90
Iterative Solvers
The Sparse Conjugate Gradient Solver
Sparse storage of G = AftAf
j, gi,j
Interactive solver !!!
White magic: Multigrid
Remember: direct solver is O(n3)
Sparse Conjugate Gradient is O(n2) !!
That’s much better, but …
… we want even more efficiency !!
White magic: Multigrid
1
2
4
1
2
3
Cascadic multigrid
Coarse to fine
[Bornemann96]
4
[Lee],[Schroeder],
[Kobbelt],[Hackbuch]
3
White magic: Multigrid
[MIPS, HLSCM, ABF++…]
Step 2: Cascadic multigrid
White Magic: Multigrid
direct solver : O(n3)
Sparse CG : O(n2)
Multigrid
: O(n) !!
[
]
Black Magic: Sparse Direct
We started from O(n3)
We achieved O(n)
Can we do better ?
-- Interactivity --- Ease of implementation --
[Sheffer et.al] SuperLU for ABF
[Botsch et.al] Interactive mesh deformation
2. Linear and Quadratic DGP
Black Magic: Sparse Direct Solvers
Super-nodal: [Demmel et.al 96]
Multi-frontal:[Lexcellent et.al 98]
[Toledo et.al]
Direct method’s revenge:
Super-Nodal data structure
M
= L
U
TAUCS, SuperLU, CHOLDMOD
j, gi,j
[
]
Black Magic: Sparse direct
Interactivity
TAUCS, SuperLU, CHOLDMOD
2. Linear and Quadratic GP
OpenNL architecture
NlLockVariables(i,a)
…
NlBeginRow()
NlAddCoefficient(i,a)
…
NlEndRow()
NlSolve()
LS with
reduced
degrees of
freedom
•Built-in (CG, GMRES, BICGSTAB)
•SuperLU
•MUMPS
•TAUCS
•…
j, gi,j
2. Linear and Quadratic GP
Applications
Gocad:
Maya
Meshing for
oil-exploration
VSP-Technology
ATARI-Infogrammes
Blender
(OpenSource)
2. Linear and Quadratic GP
Applications
OpenNL in SILO
3. Non-Linear GP
MIPS [Hormann], Stretch [Sander]
 ABF [Sheffer], ABF++ [Sheffer & Lévy]


PGP [Ray,Levy,Li,Sheffer,Alliez]

Circle Packings/Patterns [Bobenko],
[Karevych]
Conquer the non-linear world
We want to optimize a function F(x)
What can we do when F is non-linear ?
Newton’s method:
Conquer the non-linear world
Non-linear shapes functionals
Connectivity shapes
Angle Based Flattening ++
Conclusions
a map to the solvers jungle
Numerical Solvers
Direct
Naive
sparse direct
Iterative
[
] S.O.R.
Sparse C.G.
Multi-grid
Conclusions
Sparse C.G.
Multi-grid
Sparse direct
Non-linear
100x speedup w.r.t. S.O.R.
simple to implement
Linear O(n) !!! (1000x)
best for huge objects
Ultra-fast (best for interactivity)
(10000x) Big memory overhead
Difficult to tune …
Conclusion
Take home message

In most cases, TAUCS + OOC will do the job.

For small datasets, PreCG has a good
simplicity/mem. requirement/efficientcy ratio.

Use OpenNL for Matrix Assembly - Solver Abstraction
Resources

Source code & papers
on http://alice.loria.fr
– Graphite
– OpenNL