320 likes | 467 Views
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. 1970’s. 2000’s. Motivations:. Need for scalability in GP.
E N D
Mesh Parameterization:Theory and PracticeMaking it work in practice - Numerics Bruno Lévy - INRIA
Overview 1. Numerical Problems 2. Linear and Quadratic 3. Non-linear
1970’s 2000’s Motivations: Need for scalability in GP
1. Numerical Problems in GPMesh Parameterization Ui = S ai,jUj j Ni i i,j ai,j > 0 ai,i = - S ai,j j1 j… The border is mapped to a convex polygon j2 [Tutte], [Floater]
1. Numerical Problems in GPDiscrete Fairing 2 F(p)=Spi - S ai,jpj i j Ni i j1 j… j2 [Mallet], [Kobbelt], [Sorkine]…
1. Numerical Problems in GPNeighborhoods F = sum of terms, attached to neighborhoods Parameterization [Haker00] [Levy02] Curv. Estimation [Cohen-Steiner 03] Texture mapping [Levy01] Discrete fairing [Desbrun], ... Discrete fairing [Kobbelt98, Mallet95] Parameterization [Desbrun02] Deformations [CohenOr], [Sorkine] [Eck]
2. Linear and Quadratic GPRemoving degrees of freedom 2 F(x) = A x - b 2 xf xl [ Af Al] - b F(xf) =
2 2 F(xf) = A.x - d = Al.xl + Af.xf - d F(xf) minimum Aft.Af.xf = Aft.d - AftAl.xl } } M.x = b 2. Linear and Quadratic GPRemoving degrees of freedom The problem: (1) construct a linear system (2) solve a linear system
NlSolve() 2. Linear and Quadratic GPThe 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) … For each stencil instance (one-rings): NlBeginRow(); NlAddCoefficient(i, a); … NlEndRow(); • Need for • Dynamic Matrix DS • Updating formula
U M = n=106 10 centuries n=100 0.01 s (1) solve U x = b Y = b L L L U x = Y (2) solve 2. Linear and Quadratic GPDirect Solvers (LU) A Textbook solver: LU factorization (and Cholesky) a ‘small’ problem: O(n3) !!
mi,j xfj … … … mn,1 m1,1 mi,1 … … … m1,j mi,j mn,j m1,n mn,n mi,n j = i ) ( 1 ci- xfi mi,i 2. Linear and Quadratic GPSuccessive Over-Relaxation (Gauss-Seidel) xf1 …. …. …. …. …. = xfi ci …. …. …. …. …. xfnf used in [Taubin95] …
2. Linear and Quadratic DGPSuccessive Over-Relaxation (Gauss-Seidel) 1000 iterations S.O.R.
Complicated ops: Matrix x vector (see course notes) 2. Linear and Quadratic GPWhite 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) [Shewchuck: CG without the agonizing pain]
Precond. Conj. Grad. 5 - log|G.x + c | Conj. Grad. 2.5 S.O.R. 0 0 45 90 Time(s) Iterative Solvers • Successive Over-Relaxation • Sparse C.G. >100x speedup
Demo j, gi,j Iterative Solvers The Sparse Conjugate Gradient Solver Sparse storage of G = AftAf 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 !!
2 1 4 1 2 3 Cascadic multigrid Coarse to fine [Bornemann96] 4 3 White magic: Multigrid [Lee],[Schroeder], [Kobbelt],[Hackbuch]
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
U M = L j, gi,j 2. Linear and Quadratic DGPBlack 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 TAUCS, SuperLU, CHOLDMOD
[ ] Demos Black Magic: Sparse direct Interactivity TAUCS, SuperLU, CHOLDMOD
LS with reduced degrees of freedom • Built-in (CG, GMRES, BICGSTAB) • SuperLU • MUMPS • TAUCS • … j, gi,j 2. Linear and Quadratic GPOpenNL architecture NlLockVariables(i,a) … NlBeginRow() NlAddCoefficient(i,a) … NlEndRow() NlSolve()
2. Linear and Quadratic GPApplications Maya Gocad: Meshing for oil-exploration Blender (OpenSource) VSP-Technology ATARI-Infogrammes
2. Linear and Quadratic GPApplications 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]
Newton’s method: Conquer the non-linear world We want to optimize a function F(x) What can we do when F is non-linear ?
Conquer the non-linear world Non-linear shapes functionals Connectivity shapes Angle Based Flattening ++ Demos
[ ] Direct Iterative S.O.R. Naive Sparse C.G. sparse direct Multi-grid Conclusionsa map to the solvers jungle Numerical Solvers
Sparse C.G. Multi-grid Sparse direct Non-linear Conclusions 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 …
ConclusionTake 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