ECCV Tutorial Mesh Processing Numerics Bruno Lvy INRIA

ECCV Tutorial Mesh Processing Numerics Bruno Lévy INRIA - ALICE

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 = a U i, j j j Ni i j 1 j… j 2 [Tutte], [Floater] i, j ai, j > 0 ai, i = - ai, j The border is mapped to a convex polygon

1. Numerical Problems in GP Discrete Fairing F(p)= i pi - ai, jpj j N i j 1 j… i j 2 [Mallet], [Kobbelt], [Sorkine]… 2

1. Numerical Problems in GP Neighborhoods F = sum of terms, attached to neighborhoods Discrete fairing [Kobbelt 98, Mallet 95] Parameterization [Desbrun 02] Deformations [Cohen. Or], [Sorkine] Curv. Estimation [Cohen-Steiner 03] Texture mapping [Levy 01] Discrete fairing [Desbrun], . . . Parameterization [Haker 00] [Levy 02] [Eck]

2. Linear and Quadratic GP Removing degrees of freedom F(x) = A x - b 2 2 [ Af Al ] F(xf) = xf xl -b

2. Linear and Quadratic GP Removing degrees of freedom F(xf) = A. x - d 2 = Al. xl + Af. xf - d } Aft. Af. xf = Aft. d - Aft. Al. xl } F(xf) minimum 2 M. x = b The problem: (1) construct a linear system (2) solve a linear system

2. Linear and Quadratic GP The Open. NL approach (http: //alice. loria. fr/software) The problem: (1) construct a linear system (2) solve a linear system Nl. Lock. Variable(i 1, val 1) Nl. Lock. Variable(i 2, val 2) Nl. Solve() … For each stencil instance (one-rings): Nl. Begin. Row(); Need for Nl. Add. Coefficient(i, a); • Dynamic Matrix DS … • Updating formula Nl. End. Row();

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(n 3) !! 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) mn, 1 … mn, j … xfi …. …. mi, 1 … mi, j … mi, n xfi = mn, n xf nf ci- mi, j xfj) ( m j=i 1 i, i …. ci …. m 1, n xf 1 …. …. m 1, 1 … m 1, j … used in [Taubin 95] …

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 Sparse. Matrix &A, const Vector& b, Vector& x, double eps, int max_iter ){ int N = A. n() ; double t, tau, sig, rho, gam; double bnorm 2 = BLAS: : ddot (N, b, 1) ; BLAS: : ddot(N, b, 1) double err=eps*bnorm 2 eps*bnorm 2 ; mult(A, x, g); BLAS: : daxpy (N, -1. , b, 1, g, 1); BLAS: : daxpy(N, -1. , b, 1, g, 1); BLAS: : dscal (N, -1. , g, 1); BLAS: : dscal(N, -1. , g, 1); BLAS: : dcopy (N, g, 1, r, 1); BLAS: : dcopy(N, g, 1, r, 1); while ( BLAS: : ddot (N, g, 1)>err && its < max_iter) BLAS: : ddot(N, g, 1)> max_iter) { mult(A, r, p); rho=BLAS: : ddot(N, p, 1, p, 1); rho=BLAS: : ddot sig=BLAS: : ddot(N, r, 1, p, 1); sig=BLAS: : ddot tau=BLAS: : ddot (N, g, 1, r, 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); BLAS: : daxpy(N, -t, p, 1, g, 1); gam=(t*t* rho-tau)/tau; gam=(t*t*rho-tau BLAS: : dscal (N, gam, r, 1); BLAS: : dscal(N, BLAS: : daxpy (N, 1. , g, 1, r, 1); BLAS: : daxpy(N, 1. , g, 1, r, 1); ++its; ; } return its ; } Only simple vector ops (BLAS) [Shewchuck: CG without the agonizing pain] Complicated ops: Matrix x vector (see course notes)

Iterative Solvers • Successive Over-Relaxation • Sparse C. G. >100 x 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 = Aft. Af Demo j, gi, j Interactive solver !!!

White magic: Multigrid Remember: direct solver is O(n 3) Sparse Conjugate Gradient is O(n 2) !! That’s much better, but … … we want even more efficiency !!

White magic: Multigrid 1 2 4 1 2 3 Cascadic multigrid Coarse to fine [Bornemann 96] 4 [Lee], [Schroeder], [Kobbelt], [Hackbuch] 3

White magic: Multigrid [MIPS, HLSCM, ABF++…] Step 2: Cascadic multigrid

White Magic: Multigrid direct solver : O(n 3) Sparse CG : O(n 2) Multigrid : O(n) !!

[ ] Black Magic: Sparse Direct We started from O(n 3) We achieved O(n) Can we do better ? -- Interactivity --- Ease of implementation -- [Sheffer et. al] Super. LU 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, Super. LU, CHOLDMOD j, gi, j

[ ] Black Magic: Sparse direct Demos Interactivity TAUCS, Super. LU, CHOLDMOD

2. Linear and Quadratic GP Open. NL architecture Nl. Lock. Variables(i, a) … Nl. Begin. Row() Nl. Add. Coefficient(i, a) … Nl. End. Row() Nl. Solve() LS with reduced degrees of freedom • Built-in (CG, GMRES, BICGSTAB) • Super. LU • MUMPS • TAUCS • … j, gi, j

2. Linear and Quadratic GP Applications Gocad: Maya Meshing for oil-exploration VSP-Technology ATARI-Infogrammes Blender (Open. Source)

2. Linear and Quadratic GP Applications Open. NL in SILO

3. Non-Linear GP MIPS [Hormann], Stretch [Sander] n ABF [Sheffer], ABF++ [Sheffer & Lévy] n n PGP [Ray, Levy, Li, Sheffer, Alliez] n 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 ++ Demos

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 100 x speedup w. r. t. S. O. R. simple to implement Linear O(n) !!! (1000 x) best for huge objects Ultra-fast (best for interactivity) (10000 x) Big memory overhead Difficult to tune …

Conclusion Take home message n In most cases, TAUCS + OOC will do the job. n For small datasets, Pre. CG has a good simplicity/mem. requirement/efficientcy ratio. n Use Open. NL for Matrix Assembly - Solver Abstraction

Resources n Source code & papers on http: //alice. loria. fr/WIKI – – Graphite CGAL Open. Mesh Open. NL

Conclusion Mesh Processing: a wide topic Data structures Mesh repair Mesh analysis Smoothing Parameterization Mesh simplification Remeshing Freeform modeling. . . SIGGRAPH and EUROGRAPHICS tutorials (with M. Botsch, M. Pauly, L. Kobbelt and P. Alliez) http: //alice. loria. fr/WIKI/