The joy of computer graphics programming

Bruno Levy

More Related Content

1.Mathématiques - InformatiqueEurographics 2017Bruno LévyALICE Géométrie & LumièreCENTRE INRIA Nancy Grand-EstThe Joy ofComputer Graphics ProgrammingBruno Lévy, Inria researcher
2.Introducing myself…First name: BrunoFamilly name: LévyAge: approaching 45Occupation: Inria ALICE team lead (since 2004)ALICE = Geometry processing + Fabrication (S. Lefebvre)Research: navigating between graphics, math. and physicsCode: open-source (geogram),commercial (GOCAD, Vorpaline)
3.Introducing myself…First name: BrunoFamilly name: LévyAge: approaching 45 (born in 1972 = 4004+1)Occupation: Inria ALICE team lead (since 2004)ALICE = Geometry processing + Fabrication (S. Lefebvre)Research: navigating between graphics, math. and physicsCode: open-source (geogram),commercial (GOCAD, Vorpaline)
4.OVERVIEW1. Introduction, let s talk about programming2. Graphics: eye candy with GLUP3. Numerics: cranking the number cruncher4. Geometry: predicates without the agonizing painWhat s next ? Physics+Mathematics+Computing=?
5.Introduction1
6.Part. 1 On Software DesignWar stories on the design and implementationof geogram, graphite and vorpaline
7.Part. 1 On Software DesignWar stories on the design and implementationof geogram, graphite and vorpalineDocumented open-source* implementations of reference algos.+ the most important results from my group (2000 to 2017)Algorithms from >20 research articlesand two ERC projects (GOODSHAPE and VORPALINE)*Geogram (BSD) and Graphite (GLP), Vorpaline is proprietary
8.Part. 1 On Software DesignWar stories on the design and implementationof geogram, graphite and vorpalineDocumented open-source* implementations of reference algos.,+ the most important results from my group (2000 to 2017)Algorithms from >20 research articlesand two ERC projects (GOODSHAPE and VORPALINE)•Mesh parameterization (LSCM, ABF++, PGP)•Spectral mesh processing (Manifold Harmonics)•Newton Centroidal Voronoi Tesselation (surfaces and volumes)•Remeshing, reconstruction, Optimal Transport•Low-level algorithms (Delaunay, Voronoi, predicates) …*Geogram (BSD) and Graphite (GLP), Vorpaline is proprietary
9.Part. 1 On Software DesignFrom my attic: my first computer !
10.Part. 1 On Software Design1979 Apple ][6502 processor, 1MHz64Kb RAMApprox. 10 FLOPs
11.Part. 1 On Software Design1979 Apple ][6502 processor, 1MHz64Kb RAMApprox. 10 FLOPs2017 PCCore i7 gen3, 3 GHz16Gb RAMApprox. 100 GFLOPs
12.Part. 1 On Software Design1979 Apple ][6502 processor, 1MHz64Kb RAMApprox. 10 FLOPs2017 PCCore i7 gen3, 3 GHz16Gb RAMApprox. 100 GFLOPsX 1 million !!!!38 years
13.Part. 1 On Software DesignBoots in 20 seconds Boots in 3 minutes
14.Part. 1 On Software DesignBoots in 20 seconds Boots in 3 minutesWhere did the 1 million acceleration factor go ?
15.Part. 1 On Software DesignWhat can you do in 20 seconds ?3 GHz, 4 cores = 240 billions instructions !!
16.Part. 1 On Software DesignIf you cannot do the job in less than 20 secondson a modern PC, then there is probably aproblem somewhereWhat can you do in 20 seconds ?3 GHz, 4 cores = 240 billions instructions !!
17.Part. 1 On Software DesignBoots in 20 seconds Boots in 3 minutesWhere did the 1 million acceleration factor go ?- Lost in abstraction -
18.Part. 1 On Software DesignAbstraction in Programming:+ Separates concepts+ Separates specification from Implementation
19.Part. 1 On Software DesignAbstraction in Programming:+ Separates concepts+ Separates specification from Implementation- Sometimes separates thingsthat should be considered together !!
20.Part. 1 On Software DesignBackground on Futuristic ProgrammingPaul Haeberli - 1994www.graficaobscura.com/future/index.html
21.Part. 1 On Software DesignBackground on Futuristic ProgrammingPaul Haeberli - 1994
22.Part. 1 On Software DesignFuturistic Programming PrioritiesPaul Haeberli - 19941. It is something that has NEVER BEEN DONE BEFORE.2. The USER LIKES to use the program.3. The program is as FAST as it can be.4. The program is as SMALL as it can be.5. The program is BUG-FREE.6. The program needs NO USER MAINTENANCE.7. The program requires NO USER DOCUMENTATION.8. The program requires NO SYSTEM ADMINISTRATOR
23.Part. 1 On Software DesignGeogram/Graphite Programming Priorities1. Make it as simple as possible (but not simpler)
24.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possibleGeogram/Graphite Programming Priorities
25.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possible3. Make it as easy to compile as possibleGeogram/Graphite Programming Priorities
26.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possible3. Make it as easy to compile as possible4. Maximize speedGeogram/Graphite Programming Priorities
27.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possible3. Make it as easy to compile as possible4. Maximize speed5. Minimize memory consumptionGeogram/Graphite Programming Priorities
28.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possible3. Make it as easy to compile as possible4. Maximize speed5. Minimize memory consumption6. Minimize number of lines of codeGeogram/Graphite Programming Priorities
29.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possible3. Make it as easy to compile as possible4. Maximize speed5. Minimize memory consumption6. Minimize number of lines of code7. Minimize number of C++ classesGeogram/Graphite Programming Priorities
30.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possible3. Make it as easy to compile as possible4. Maximize speed5. Minimize memory consumption6. Minimize number of lines of code7. Minimize number of C++ classesGeogram/Graphite Programming PrioritiesSimplicity is theultimate sophistication
31.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possible3. Make it as easy to compile as possible4. Maximize speed5. Minimize memory consumption6. Minimize number of lines of code7. Minimize number of C++ classesGeogram/Graphite Programming PrioritiesSimplicity is theultimate sophistication
32.Part. 1 On Software Design1. Make it as simple as possible (but not simpler)2. Make it as easy to use as possible3. Make it as easy to compile as possible4. Maximize speed5. Minimize memory consumption6. Minimize number of lines of code7. Minimize number of C++ classes8. Systematically document all classes, all functions, [+Biblio.]9. Systematically document the implementation of all algorithms10. Assertion checks everywhere11. Zero warnings with all compilers / platforms12. Perform systematic non-regression testing and mem. check.Geogram/Graphite Programming Priorities
33.Part. 1 On Software DesignFuturistic programmingUsefullness (and coolness) are primary !
34.Part. 1 On Software Design• Jonathan Shewchuk s Triangle and exact predicates• Tetgen, MGTetra, MMG3d• Omar Cornut s ImGUI library• David Mount s ANN library• The LUA prog. languageFuturistic programmingExamples of Futuristic codesUsefullness (and coolness) are primary !
35.Part. 1 On Software DesignCase study: mesh data structuresFrom several Computer Graphics / Mesh Processing 101 courses- including (earlier versions of) mine -
36.Halfedges Design Principles1. Individual combinatorial elementscan be created/destroyed at any timein constant time2. Basic operations (collapse, split, )3. Higher-level operations on top of themPart. 1 On Software DesignCase study: mesh data structuresHalfedges and edgeuses, harmful or useful ?
37.+ Benefit of abstraction: layered designPart. 1 On Software DesignCase study: mesh data structuresHalfedges Design Principles1. Individual combinatorial elementscan be created/destroyed at any timein constant time2. Basic operations (collapse, split, )3. Higher-level operations on top of themHalfedges and edgeuses, harmful or useful ?
38.+ Benefit of abstraction: layered design- Separates things that should have been considered togetherPart. 1 On Software DesignCase study: mesh data structuresHalfedges and edgeuses, harmful or useful ?Halfedges Design Principles1. Individual combinatorial elementscan be created/destroyed at any timein constant time2. Basic operations (collapse, split, )3. Higher-level operations on top of them
39.+ Benefit of abstraction: layered design- Separates things that should have been considered togetherNot the optimal level of granularityPart. 1 On Software DesignCase study: mesh data structuresHalfedges and edgeuses, harmful or useful ?Halfedges Design Principles1. Individual combinatorial elementscan be created/destroyed at any timein constant time2. Basic operations (collapse, split, )3. Higher-level operations on top of them
40.Part. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
41.Part. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
42.Part. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
43.Part. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
44.Part. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
45.Part. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
46.Objection ! (?) Deletion of one individual element is O(n) instead of constantIn fact: deleting 1 element = wrong granularity level !Part. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
47.Deleting a bunch of elementsPart. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
48.Deleting a bunch of elementsPart. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
49.Deleting a bunch of elementsPart. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
50.Deleting a bunch of elements – operates also in O(n) – right granularityPart. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)
51.Deleting a bunch of elements – operates also in O(n) – right granularityPart. 1 On Software DesignCase study: mesh data structuresNeeds array permutation (not in the STL unfortunately,but easy to implement, see GEOGRAM::permutation).Indexed Mesh data structure (no data structure)
52.Indexed Mesh data structure (no data structure)Summary:•An array of vertices coordinates•An array of triangle vertices indicesPart. 1 On Software DesignCase study: mesh data structures
53.•An array of vertices coordinates•An array of triangle vertices indices•An array of triangle adjacencies (optional)•An array of facet first indices (optional)Part. 1 On Software DesignCase study: mesh data structuresIndexed Mesh data structure (no data structure)Summary:
54.Benefits of such anon-datastructure, non-object, non-oriented, (non-)programming1. Simpler codePart. 1 On Software DesignCase study: mesh data structures
55.Benefits of such anon-datastructure, non-object, non-oriented, (non-)programming1. Simpler code2. Less memory;Part. 1 On Software DesignCase study: mesh data structures
56.Benefits of such anon-datastructure, non-object, non-oriented, (non-)programming1. Simpler code2. Less memory;3. Parallelization #pragma omp parallel forPart. 1 On Software DesignCase study: mesh data structures
57.Benefits of such anon-datastructure, non-object, non-oriented, (non-)programming1. Simpler code2. Less memory;3. Parallelization #pragma omp parallel for4. Easy copy: memcpy()Part. 1 On Software DesignCase study: mesh data structures
58.Benefits of such anon-datastructure, non-object, non-oriented, (non-)programming1. Simpler code2. Less memory;3. Parallelization #pragma omp parallel for4. Easy copy: memcpy()5. Easy I/O: fread()/fwrite()Part. 1 On Software DesignCase study: mesh data structures
59.Benefits of such anon-datastructure, non-object, non-oriented, (non-)programming1. Simpler code2. Less memory;3. Parallelization #pragma omp parallel for4. Easy copy: memcpy()5. Easy I/O: fread()/fwrite()6. Properties/attributes are simply additional arraysPart. 1 On Software DesignCase study: mesh data structures
60.Benefits of such anon-datastructure, non-object, non-oriented, (non-)programming1. Simpler code2. Less memory;3. Parallelization #pragma omp parallel for4. Easy copy: memcpy()5. Easy I/O: fread()/fwrite()6. Properties/attributes are simply additional arrays7. Directly understood by OpenGL: VertexBufferObjectPart. 1 On Software DesignCase study: mesh data structures
61.A mesh = a bunch of arrays (std::vectors)Part. 1 On Software DesignCase study: mesh data structures
62.A mesh = a bunch of arrays (std::vectors)How do you iterate on a vector ?Part. 1 On Software DesignCase study: mesh data structures
63.Pre-2011:for(std::vector<Thing>::iterator it = V.begin(); it!=V.end(); ++it) {do something with *it}Part. 1 On Software DesignCase study: mesh data structuresHow do you iterate on a vector ?
64.Pre-2011:for(std::vector<Thing>::iterator it = V.begin(); it!=V.end(); ++it) {do something with *it}Part. 1 On Software DesignCase study: mesh data structuresHow do you iterate on a vector ?It s a pain to typeClutters the source-code (less legible)
65.Pre-2011:for(std::vector<Thing>::iterator it = V.begin(); it!=V.end(); ++it) {do something with *it}2011:for(auto it = V.begin(); it!=V.end(); ++it) {do something with *it}Part. 1 On Software DesignCase study: mesh data structuresHow do you iterate on a vector ?
66.Pre-2011:for(std::vector<Thing>::iterator it = V.begin(); it!=V.end(); ++it) {do something with *it}2011:for(auto it = V.begin(); it!=V.end(); ++it) {do something with *it}now:for(auto&& i : V) {do something with i}Part. 1 On Software DesignCase study: mesh data structuresHow do you iterate on a vector ?
67.A mesh = a bunch of arrays (std::vectors)How do you iterate on a vector ?Part. 1 On Software DesignCase study: mesh data structuresWarning: flying tomatoes alert ahead,(modern C++ lovers might throw tomatoes at me)!
68.Pre-2011:for(std::vector<Thing>::iterator it = V.begin(); it!=V.end(); ++it) {do something with *it}2011:for(auto it = V.begin(); it!=V.end(); ++it) {do something with *it}now:for(auto&& i : V) {do something with I}Part. 1 On Software DesignCase study: mesh data structuresHow do you iterate on a vector ?
69.How I iterate on a vector:for(uint i=0; i<V.size(); ++i) {do something with V[i];}Part. 1 On Software DesignCase study: mesh data structuresHow do you iterate on a vector ?
70.How I iterate on a vector:for(uint i=0; i<V.size(); ++i) {do something with V[i];}Part. 1 On Software DesignCase study: mesh data structures+ Easy to understand, even by C-only and Fortran programmers+ Compatible with all compilers+ #pragma omp parallel-for friendly* (and also better vectorization)- 15 additional keystrokes as compared to modern C++ range-forHow do you iterate on a vector ?*But use ints instead of uints, omp does not like uints
71.How I iterate on a vector:for(uint i=0; i<V.size(); ++i) {do something with V[i];}Part. 1 On Software DesignCase study: mesh data structures+ Easy to understand, even by C-only programmers+ Compatible with all compilers+ #pragma omp parallel-for friendly- 15 additional keystrokes as compared to modern C++ range-forThe following code has been approved forAPPROPRIATE AUDIENCESPG-13
72.How I iterate on a vector:for(uint i=0; i<V.size(); ++i) {do something with V[i];}Part. 1 On Software DesignCase study: mesh data structures+ Easy to understand, even by C-only programmers+ Compatible with all compilers+ #pragma omp parallel-for friendly- 15 additional keystrokes as compared to modern C++ range-for#define FOR(i,N) for(uint i=0; i<(N); ++i)FOR(i,V.size()) {do something with V[i];}+ Easy to understand, legible+ Trivial iterations easy to spot- Macros are evil[Nicolas Ray]
73.Part. 1 On Software DesignCase study: mesh data structuresObjection:It is bad because it is not flexible,what if you want to adapt your algorithm to another container ?
74.Part. 1 On Software DesignCase study: mesh data structuresObjection:It is bad because it is not flexible,what if you want to adapt your algorithm to another container ?Answer:I m not going to use something else than a vector, because it is the bestchoice for the mesh data structure. Why keeping a tuning knob on thedash board if it is always on the same position ?modern/generic != futuristic
75.Part. 1 On Software DesignCase study: mesh data structuresObjection:It is bad because it is not flexible,what if you want to adapt your algorithm to another container ?Answer:I m not going to use something else than a vector, because it is the bestchoice for the mesh data structure. Why keeping a tuning knob on thedash board if it is always on the same position ? (think of the I-Phone)The Nokia N95, amodern/generic phone.The I-phone,a futuristic phone.modern/generic != futuristic
76.Part. 1 On Software DesignCase study: mesh data structuresmodern/generic != futuristicIt is good because it haseverything that you needThe Nokia N95, amodern/generic phone.The I-phone,a futuristic phone.
77.Part. 1 On Software DesignCase study: mesh data structuresmodern/generic != futuristicIt is good because it haseverything that you needIt is even better because it hasnothing else than what you needThe Nokia N95, amodern/generic phone.The I-phone,a futuristic phone.
78.Part. 1 On Software DesignCase study: mesh data structuresmodern/generic != futuristicIt is good because it haseverything that you needIt is even better because it hasnothing else than what you needWork is finished when you havenothing to add and nothing to remove !The Nokia N95, amodern/generic phone.The I-phone,a futuristic phone.
79.Part. 1 On Software DesignCase study: mesh data structuresWhy keeping a tuning knob on the dash board if it is alwayson the same position ? (think of the I-Phone)A parameter that always has the same value is not a parameter and should beremoved from the API.A template that is instanced only once should not be a template.
80.Part. 1 On Software DesignCase study: mesh data structuresA parameter that always has the same value is not a parameter and should beremoved from the API.A template that is instanced only once should not be a template.Objection: but we loose genericity if we do that ???Why keeping a tuning knob on the dash board if it is alwayson the same position ? (think of the I-Phone)
81.Part. 1 On Software DesignCase study: mesh data structuresA parameter that always has the same value is not a parameter and should beremoved from the API.A template that is instanced only once should not be a template.Objection: but we loose genericity if we do that ???Answer to objection: but we gain a lot, it declutters thecode, makes it more legible, reduces compilation times,makes C++ compilation error messages more legible.Why keeping a tuning knob on the dash board if it is alwayson the same position ? (think of the I-Phone)
82.Part. 1 On Software DesignCase study: mesh data structuresRule of thumb:Make it a parameter not before you need it with at least two different values.Make it a template not before you need at least two different instanciations.*regarding compilation time, legibiity of error messages and run-time flex.
83.Part. 1 On Software DesignCase study: mesh data structuresRule of thumb:Make it a parameter not before you need it with at least two different values.Make it a template not before you need at least two different instanciations.About templates, consider less annoying* alternatives, such as(1) object oriented programming / virtual functions(2) or simply a parameter and if() statements*regarding compilation time, legibility of error messages and run-time flex.
84.Graphics2Fast and Easy eye-candy with GLUP
85.Part. 2 Graphics – endangered speciesThe Windows start menu The I-Phone jackglBegin(GL_TRIANGLES)glVertex(… )glEnd()OpenGL immediate modeglPush/Pop/MultMatrix()glLight()OpenGL fixed functionality pipeline
86.Part. 2 Graphics – endangered speciesglBegin(GL_TRIANGLES)glVertex(… )glEnd()OpenGL immediate modeglPush/Pop/MultMatrix()glLight()OpenGL fixed functionality pipelineDifficulties for an undergraduate who starts:(1) Assemble vertex buffer objects(2) Design vertex and fragment shaders(3) (+ the Professor, I see nothing syndrom, nothing new here)R.I.P. R.I.P.
87.Part. 2 Graphics – GLUPImmediate mode + fixed functionality pipeline implementedin modern OpenGL (4.x) (nearly source-compatible)glupBegin(), glupEnd(), glupVertex(), glupPushMatrix()…
88.Part. 2 Graphics – GLUPImmediate mode + fixed functionality pipeline implementedin modern OpenGL (4.x) (nearly source-compatible)glupBegin(), glupEnd(), glupVertex(), glupPushMatrix()…normals not needed (per-fragment shading)
89.Part. 2 Graphics – GLUPImmediate mode + fixed functionality pipeline implementedin modern OpenGL (4.x) (nearly source-compatible)glupBegin(), glupEnd(), glupVertex(), glupPushMatrix()…normals not needed (per-fragment shading)new volumetric primitives:GLUP_TETRAHEDRA, GLUP_HEXAHEDRA,GLUP_PRISMS, GLUP_PYRAMIDS, GLUP_SPHERES
90.Part. 2 Graphics – GLUPImmediate mode + fixed functionality pipeline implementedin modern OpenGL (4.x) (nearly source-compatible)glupBegin(), glupEnd(), glupVertex(), glupPushMatrix()…normals not needed (per-fragment shading)new volumetric primitives:GLUP_TETRAHEDRA, GLUP_HEXAHEDRA,GLUP_PRISMS, GLUP_PYRAMIDS, GLUP_SPHERESnew clipping modes (on-GPU marching cells algorithm)
91.Part. 2 Graphics – GLUPImmediate mode + fixed functionality pipeline implementedin modern OpenGL (4.x) (nearly source-compatible)glupBegin(), glupEnd(), glupVertex(), glupPushMatrix()…normals not needed (per-fragment shading)new volumetric primitives:GLUP_TETRAHEDRA, GLUP_HEXAHEDRA,GLUP_PRISMS, GLUP_PYRAMIDS, GLUP_SPHERESnew clipping modes (on-GPU marching cells algorithm)glupEnable(GLUP_DRAW_MESH)
92.Part. 2 Graphics – GLUPglupDrawArrays(), glupDrawElements(), VBOsVanillaGL (1.x)OpenGL ES,WebGLGLSL 1.5 (GL3.3)GLSL 4.4 (GL4.4)POINTS LINES TRGLS QUADS SPHERES T ETS PRISMS HEXES PYRAMIDSglupDrawArrays()/glupDrawElements()/VBOs and immediate modeImmediate mode only
93.Part. 2 Graphics – GLUP
94.Part. 2 Graphics – GLUP
95.Part. 2 Graphics – GLUPglupDrawElements(GLUP_TETRAHEDRA)1 single draw call with VBOs, everything occurs on GPUglupSetClippingMode(GLUP_CLIP_SLICE)
96.Part. 2 Graphics – GLUP
97.Part. 2 Graphics – GLUPglupDrawElements(GLUP_SPHERES)1 single draw call with VBOs, everything occurs on GPU(GLUP + screen-space ambient occlusion in Graphite)
98.Part. 2 Graphics – GLUP
99.Part. 2 Graphics – GLUPglupDrawElements(GLUP_HEXAHEDRA)glupDrawElements(GLUP_TETRAHEDRA)2 draw calls with VBOs, everything occurs on GPU(GLUP + screen-space ambient occlusion in Graphite)Hexahedral-dominant mesh [Ray, Sokolov, Unterreiner, L 2016]
100.Part. 2 Graphics – GLUPglupDrawElements(GLUP_HEXAHEDRA)glupDrawElements(GLUP_TETRAHEDRA)2 draw calls with VBOs, everything occurs on GPU(GLUP + screen-space ambient occlusion in Graphite)Cross-section computed with on-GPU marching cells with interpolated attrib.
101.Numerics3Cranking the number cruncher
102.Part. 3. NumericsNumerical problems: neighborhoodsF = sum of terms, attached to neighborhoodsDiscrete fairing[Kobbelt98, Mallet95]Parameterization[Desbrun02]Deformations[CohenOr], [Sorkine]Curv. Estimation[Cohen-Steiner 03]Texture mapping[L01]Discrete fairing[Desbrun], ...Parameterization[Haker00][L02][Eck]
103.Part. 3. NumericsNumerical problems: mesh parameterizationij1j2j…Ui = ai,jUjj  Ni i,j ai,j > 0ai,i = -  ai,jThe border is mapped toa convex polygon[Tutte], [Floater]
104.Part. 3. NumericsNumerical problems: mesh fairingij1j2j…F(p)=  pi - ai,jpj2j  Nii[Mallet], [Kobbelt], [Sorkine]
105.Part. 3. NumericsNumerical problems: mesh fairing
106.Part. 3. NumericsNumerical problems: mesh fairingF(x) =2A x - b
107.Part. 3. NumericsNumerical problems: mesh fairingF(xf) =xfxl[ Af Al] - b2F(x) =2A x - b
108.Part. 3. NumericsNumerical problems: least squaresF(xf) = A.x - d = Al.xl + Af.xf - d2 2
109.Part. 3. NumericsNumerical problems: least squaresF(xf) = A.x - d = Al.xl + Af.xf - d2 2F(xf) minimum Aft.Af.xf = Aft.d - AftAl.xlM.x = b}}
110.Part. 3. NumericsNumerical problems: least squaresF(xf) = A.x - d = Al.xl + Af.xf - d2 2F(xf) minimum Aft.Af.xf = Aft.d - AftAl.xlM.x = b}}The problem:(1) construct the linear system (assembly)(2) solve a linear system
111.Part. 3. Numerics – the OpenNL library.Numerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear system
112.Part. 3. Numerics – the OpenNL library.Numerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear systemSimilarity between a mesh and a sparse matrix:
113.Part. 3. Numerics – the OpenNL library.Numerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear systemnlBegin(NL_ROW)nlAddCoefficient(I,j,val)nlRightHandSide(val)nlEnd(NL_ROW)Similarity between a mesh and a sparse matrix:OpenNL: API inspired by OpenGL immediate mode.(+ automatic construction of AtA for least-squares)
114.Part. 3. Numerics – the OpenNL library.Numerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear systemNLSparseMatrix: dynamic data structure (can grow)OpenNL internals……………
115.Part. 3. Numerics – the OpenNL library.Numerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear systemNLSparseMatrix NLCRSMatrixOpenNL internalsnlCompress()
116.Part. 3. Numerics – the OpenNL library.Numerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear systemNLSparseMatrix NLCRSMatrix NLCusparseMatrixOpenNL internals GPUnlCompress() (optionnal)
117.Part. 3. Numerics – the OpenNL libraryNumerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear systemnlSolve() Iterative solversCGGMResBiCGStabCPU/GPUAbstraction layerOpenMPmulticoreCUDACuBLASCuSparse
118.Part. 3. Numerics – the OpenNL libraryNumerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear systemnlSolve() Iterative solversDirect solversCGGMResBiCGStabSuperLUCHOLDMODCPU/GPUAbstraction layerOpenMPmulticoreCUDACuBLASCuSparse
119.Part. 3. Numerics – the OpenNL libraryNumerical problems: least squaresThe problem:(1) construct the linear system (assembly)(2) solve a linear systemnlSolve()Iterative solversDirect solversCGGMResBiCGStabSuperLUCHOLDMODCPU/GPUAbstraction layerOpenMPmulticoreCUDACuBLASCuSparseEigen solver ARPACK
120.Part. 3. Numerics – the OpenNL libraryOpenNL as a pluggable software moduleFollowing futuristic programming principles:* OpenNL also available as a single .c,.h pair, portableto all architectures, easy to compile.* Object-Oriented abstract matrix interface in C (achievesrun-time CPU/GPU flexibility)* Dynamically loads CUDA CuBLAS and CuSparseon demand (as well as CHOLMOD, SUPERLU,ARPACK)* Double precision (and no single precision).
121.Part. 3. Numerics – the OpenNL libraryOpenNL at work Everything available in geogramGUI in graphite
122.Part. 3. Numerics – the OpenNL libraryOpenNL at workLSCM,Spectral parameterization.(600 lines of code for both)ABF++ (800 lines of code)PGP, QuadCover, MIP.(700 lines of code)Manifold Harmonics, HKS, ADF. (400 lines of code)Everything available in geogramGUI in graphite
123.Part. 3. Numerics – the OpenNL libraryOpenNL at work Everything available in geogramGUI in graphiteBaby groot © Marvel (this version by Byambaa on Thingyverse)
124.Part. 3. Numerics – the OpenNL libraryOpenNL at work Everything available in geogramGUI in graphiteAutomatic decimation and normalmap generation.Baby groot © Marvel (this version by Byambaa on Thingyverse)
125.Geometry4Predicates without the agonizing pain**Section title is a tribute to Jonathan Shewchuk.
127.4. Geometry - Motivations[Martinez, Dumas, Lefebvre]
128.4. Geometry - Motivations
129.1. Motivations
130.4.Geometry MotivationsRené DescartesWonderShapesStructuresSymmetry
131.Part. 4. GeometryComputing (restricted) Voronoi diagrams
132.Part. 4. GeometryComputing (restricted) Voronoi diagrams
133.Part. 4 GeometryPointset X Simplicial complex Seither triangulated surfaceor tetrahedralized solidThe input
134.Part. 4 GeometryVor(X)|S (Intersection between Vor(X) and S)The output
135.Part 4. Geometry – a difficult datasetLots of degeneracies:Voronoi diagram with degree 4 vertices.Voronoi cell faces match exactly facets of the initial surface.
136.Part. 4 GeometryVoronoi cells as iterative convex clippingHalf-space clippingxix1x2x3x4x5x6x7x8x9x10x11
137.Part. 4 Difficulties – predicatesxixjElementary operation: cut a polygon(or polyhedron) with a bisector
138.Part. 4 Difficulties – predicatesxixjElementary operation: cut a polygon(or polyhedron) with a bisectorClassify the vertices of the polygon
139.Part. 4 Difficulties – predicatesxixjElementary operation: cut a polygon(or polyhedron) with a bisectorClassify the vertices of the polygonSign( d(p,xj) – d(p,xi)) > 0
140.Part. 4 Difficulties – predicatesxixjElementary operation: cut a polygon(or polyhedron) with a bisectorClassify the vertices of the polygonSign( d(p,xj) – d(p,xi)) > 0Sign( d(p,xj) – d(p,xi)) < 0
141.Part. 4 Difficulties – predicatesxixjElementary operation: cut a polygon(or polyhedron) with a bisectorClassify the vertices of the polygonCompute the intersectionsSign( d(p,xj) – d(p,xi)) > 0Sign( d(p,xj) – d(p,xi)) < 0
142.Part. 4 Difficulties – predicatesxixjSign( d(p,xj) – d(p,xi)) > 0Sign( d(p,xj) – d(p,xi)) < 0Elementary operation: cut a polygon(or polyhedron) with a bisectorClassify the vertices of the polygonCompute the intersections - discard
143.Part. 4 Difficulties – predicatesxixjxk Now clipping with the bisector of (xi, xk)
144.Part. 4 Difficulties – predicatesxixjxk Now clipping with the bisector of (xi, xk)We need to classify all the points, includingthese ones !
145.Part. 4 Difficulties – predicatesxixjxk Now clipping with the bisector of (xi, xk)We need to classify all the points, includingthese ones !(they are intersections between asegment and a bisector)
146.Part. 4 Difficulties – predicatesxixjxk Now clipping with the bisector of (xi, xk)We need to classify all the points, includingthese ones !(they are intersections between asegment and a bisector)This generates a new intersection(between a facet and two bisector)
147.Part. 4 Difficulties – predicatesThree configurations1) Side(xi,xj,q) where q is a vertex of S
148.Part. 4 Difficulties – predicatesThree configurations1) Side1(xi,xj,q)
149.Part. 4 Difficulties – predicatesThree configurations1) Side1(xi,xj,q)2) Side(xi,xj,q) where q = π(i,k) ∩ [p1,p2]
150.Part. 4 Difficulties – predicatesThree configurations1) Side1(xi,xj,q)2) Side2(xi,xj,xk,p1,p2)
151.Part. 4 Difficulties – predicatesThree configurations1) Side1(xi,xj,q)2) Side2(xi,xj,xk,p1,p2)3) Side(xi,xj,q) whereq = π(i,k) ∩ π(i,l) ∩ [p1,p2,p3]
152.Part. 4 Difficulties – predicatesThree configurations1) Side1(xi,xj,q)2) Side2(xi,xj,xk,p1,p2)3) Side3(xi,xj,xk, xl,p1,p2,p3)
153.Part. 4 Difficulties – predicatesThree configurations1) Side1(xi,xj,q)2) Side2(xi,xj,xk,p1,p2)3) Side3(xi,xj,xk, xl,p1,p2,p3)Implementations of exact predicates:- J. Shewchuk s code- CGAL (Pion, Meyer)
154.Part. 4 Difficulties – predicatesThree configurations1) Side1(xi,xj,q)2) Side2(xi,xj,xk,p1,p2)3) Side3(xi,xj,xk, xl,p1,p2,p3)Implementations of exact predicates:- J. Shewchuk s code- CGAL (Pion, Meyer)They do not have Side1(), Side2(), Side3() ( exotic predicates )
155.Part. 4 Difficulties – predicatesThree configurations1) Side1(xi,xj,q)2) Side2(xi,xj,xk,p1,p2)3) Side3(xi,xj,xk, xl,p1,p2,p3)How to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()
156.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Wish list:• Easy to use(no Guru needed for each new predicate)• Reasonably efficient• Easy to compile/integrate ( Futuristic programming )(multi_precision.h, multi_precision.cpp and that s all)
157.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #1: (dense) multi-precision (GMP)a0a1a2a3x 20x 21*32x 22*32x 23*32…Each number is an array of (32 bits) integers:Implement +,-,* (reasonably easy)Sign: look at the leading non-zero component
158.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #1: (dense) multi-precision (GMP)a limitation :c = a10 * 210*32 + b0
159.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #1: (dense) multi-precision (GMP)a limitation :b000a10x 20x 210*32c = a10 * 210*32 + b0…
160.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #2: (sparse) multi-precisionb0 | 0a10 | 10c = a10 * 210*32 + b0Store the exponents of the components
161.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #2: (sparse) multi-precisionb0 | 0a10 | 10c = a10 * 210*32 + b0Exp. Exp.
162.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #2: (sparse) multi-precisionb0 | 0a10 | 10c = a10 * 210*32 + b0Exp. Exp.mantissa mantissa
163.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #2: (sparse) multi-precisionb0 | 0a10 | 10c = a10 * 210*32 + b0Exp. Exp.mantissa mantissaThese are floating point numbers !!!
164.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)x3 x2 x1…Each number is represented by the sum of an array of componentsThese are floating point numbers !!!
165.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)x3 x2 x1…Each number is represented by the sum of an array of componentsThey are sorted in decreasing exponentsThese are floating point numbers !!!
166.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)x3 x2 x1…Each number is represented by the sum of an array of componentsThey are sorted in decreasing exponentsThey are non-overlappingThese are floating point numbers !!!
167.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)x3 x2 x1…Each number is represented by the sum of an array of componentsThey are sorted in decreasing exponentsThey are non-overlappingThe sign is determined by the first component (highest exponent)These are floating point numbers !!!
168.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)x2 x1Two_sum(double a, double b)
169.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)x2 x1Two_sum(double a, double b)+Length:l+mLength l Length m
170.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)x2 x1Two_sum(double a, double b)+*Length:l+mLength:2*lA doubleLength l
171.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)* …Length: 2*l*mExpansion * Expansion product implemented by a recursive function( distillation )Length l Length m
172.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)* …Expansion * Expansion product implemented by a recursive function( distillation )Performance ? 10 to 40 times slower than standard doublesLength: 2*l*mLength l Length m
173.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)* …Expansion * Expansion product implemented by a recursive function( distillation )Performance ? 10 to 40 times slower than standard doublesUse arithmetic filtersAdaptive precision [Shewchuk] ? Too complicated to get rightLength: 2*l*mLength l Length m
174.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()Idea #3: expansions (Shewchuk)* …Expansion * Expansion product implemented by a recursive function( distillation )Performance ? 10 to 40 times slower than standard doublesUse arithmetic filtersAdaptive precision [Shewchuk] ? Too complicated to get rightQuasi-static filters [Meyer and Pion] – FPG generator (easy to use)Length: 2*l*mLength l Length m
175.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()PCK (Predicate Construction Kit)multi_precision.h / multi_precision.cppA low-level expansion class (allocates expansions on stack)A high-level C++ number type (+,-,*,Sign overloads)a compiler that generates the filter with FPG [Meyer and Pion] and theexact precision version with expansions
176.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()PCK (Predicate Construction Kit)multi_precision.h / multi_precision.cppA low-level expansion class (allocates expansions on stack)A high-level C++ number type (+,-,*,Sign overloads)a compiler that generates the filter with FPG [Meyer and Pion] and theexact precision version with expansions(why not templates / metaprogramming ?would be possible, but specialized language is much better here).
177.Part. 4 Exact ArithmeticsHow to implement Side1(), Side2(), Side3() ?We need an exact number type with +,-,*,Sign()PCK (Predicate Construction Kit)multi_precision.h / multi_precision.cppA low-level expansion class (allocates expansions on stack)A high-level C++ number type (+,-,*,Sign overloads)a compiler that generates the filter with FPG [Meyer and Pion] and theexact precision version with expansionsSo we are done ?
178.Part. 4 Symbolic PerturbationHow to implement Side1(), Side2(), Side3() ?xixjNot yet !!
179.Part. 4 Symbolic PerturbationHow to implement Side1(), Side2(), Side3() ?xixjNot yet !!
180.Part. 4 Symbolic Perturbation
181.Part. 4 Symbolic Perturbation
182.Part. 4 Symbolic Perturbation
183.Part. 4 Symbolic Perturbation
184.Part. 4 Symbolic Perturbation
185.Part. 4 Symbolic Perturbation
186.Part. 4 Symbolic Perturbation
187.Part. 4 Symbolic Perturbation
188.Part. 4 Symbolic Perturbation
189.Part. 4 Symbolic Perturbation
190.Part. 4 Symbolic Perturbation
191.Part. 4 Symbolic Perturbation
192.Part. 4 Symbolic Perturbation
193.Part. 4 Symbolic Perturbation
194.Part. 4 Symbolic Perturbation
195.Part. 4 Symbolic Perturbation
196.[Voronoi][Edelsbrunner et.al][Devillers et.al]Part. 4 Symbolic Perturbation
197.[Voronoi][Edelsbrunner et.al][Devillers et.al]Part. 4 Symbolic PerturbationEach point pi is replacedby a trajectory pi(ε)xi(ε) = xi + ε3iyi(ε) = yi + ε3i+1zi(ε) = zi + ε3i+2For instance:
198.[Voronoi][Edelsbrunner et.al][Devillers et.al]Part. 4 Symbolic PerturbationEach point pi is replacedby a trajectory pi(ε)Take the limit whenε tends to 0
199.xixjπ(i,j) = {p | d2(p,xi) = d2(p,xj)}[Voronoi][Edelsbrunner et.al][Devillers et.al]Part. 4 Symbolic PerturbationIn our case, perturb theweights of a power diagram.
200.xixjπw(i,j) = {p | d2(p,xi) - wi = d2(p,xj) - wj}[Voronoi][Edelsbrunner et.al][Devillers et.al]Part. 4 Symbolic PerturbationIn our case, perturb theweights of a power diagram.
201.xixjπw(i,j) = {p | d2(p,xi) - wi = d2(p,xj) - wj}[Voronoi][Edelsbrunner et.al][Devillers et.al]The Voronoi diagram is replaced with a power diagramPart. 4 Symbolic Perturbation
202.xixjπw(i,j) = {p | d2(p,xi) - wi = d2(p,xj) - wj}[Voronoi][Edelsbrunner et.al][Devillers et.al]The Voronoi diagram is replaced with a power diagramSymbolic perturbation – Simulation of Simplicity:Define the weight as a function of ε: wi = εiPart. 4 Symbolic Perturbation
203.xixjπw(i,j) = {p | d2(p,xi) - wi = d2(p,xj) - wj}[Voronoi][Edelsbrunner et.al][Devillers et.al]The Voronoi diagram is replaced with a power diagramSymbolic perturbation – Simulation of Simplicity:Define the weight as a function of ε: wi = εiThe combinatorics is determined by the limit ε →0Part. 4 Symbolic Perturbation
204.How to write side1(), side2(), side3(), side4() ?d2(q,pj)-wj – d2(q,pi) + wi where:q = πw(i,k1) ∩ …πw(i,kd) ∩ [p1,p2,p3…pd]Part. 4 Symbolic Perturbation
205.How to write side1(), side2(), side3(), side4() ?d2(q,pj)-wj – d2(q,pi) + wi where:q = πw(i,k1) ∩ …πw(i,kd) ∩ [p1,p2,p3…pd]Solve for q in:q Є πw(i,k1)…q Є πw(i,kd)q Є [p1,p2,p3…pd]{Part. 4 Symbolic Perturbation
206.How to write side1(), side2(), side3(), side4() ?d2(q,pj)-wj – d2(q,pi) + wi where:q = πw(i,k1) ∩ …πw(i,kd) ∩ [p1,p2,p3…pd]q = (1/d) QKeep numerator and denom.separate (remember, we arenot allowed to divide)Solve for q in:q Є πw(i,k1)…q Є πw(i,kd)q Є [p1,p2,p3…pd]{Part. 4 Symbolic Perturbation
207.How to write side1(), side2(), side3(), side4() ?d2(q,pj)-wj – d2(q,pi) + wi where:q = πw(i,k1) ∩ …πw(i,kd) ∩ [p1,p2,p3…pd]q = (1/d) QKeep numerator and denom.separate (remember, we arenot allowed to divide)Inject q=(1/d) Q in side1() and mutliply by d to remove the divisionSolve for q in:q Є πw(i,k1)…q Є πw(i,kd)q Є [p1,p2,p3…pd]{Part. 4 Symbolic Perturbation
208.How to write side1(), side2(), side3(), side4() ?d2(q,pj)-wj – d2(q,pi) + wi where:q = πw(i,k1) ∩ …πw(i,kd) ∩ [p1,p2,p3…pd]q = (1/d) QKeep numerator and denom.separate (remember, we arenot allowed to divide)Inject q=(1/d) Q in side1() and mutliply by d to remove the divisionOrder the terms in wi = εiSolve for q in:q Є πw(i,k1)…q Є πw(i,kd)q Є [p1,p2,p3…pd]{Part. 4 Symbolic Perturbation
209.How to write side1(), side2(), side3(), side4() ?d2(q,pj)-wj – d2(q,pi) + wi where:q = πw(i,k1) ∩ …πw(i,kd) ∩ [p1,p2,p3…pd]q = (1/d) QKeep numerator and denom.separate (remember, we arenot allowed to divide)Inject q=(1/d) Q in side1() and mutliply by d to remove the divisionOrder the terms in wi = εi the constant one is non-perturbed predicateif zero, the first non-zero one determines the signSolve for q in:q Є πw(i,k1)…q Є πw(i,kd)q Є [p1,p2,p3…pd]{Part. 4 Symbolic Perturbation
210.Part. 4 Symbolic Perturbation – side1#include "kernel.pckh"Sign predicate(side1)(point(p0), point(p1), point(q0) DIM) {scalar r = sq_dist(p0,p1) ;r -= 2*dot_at(p1,q0,p0) ;generic_predicate_result(sign(r)) ;begin_sos2(p0,p1)sos(p0,POSITIVE)sos(p1,NEGATIVE)end_sos}Source PCK
211.Part. 4 Symbolic Perturbation – side1#include "kernel.pckh"Sign predicate(side1)(point(p0), point(p1), point(q0) DIM) {scalar r = sq_dist(p0,p1) ;r -= 2*dot_at(p1,q0,p0) ;generic_predicate_result(sign(r)) ;begin_sos2(p0,p1)sos(p0,POSITIVE)sos(p1,NEGATIVE)end_sos}Source PCK(p1-p0).(q0-p0)
212.Part. 4 Symbolic Perturbation – side2#include "kernel.pckh“Sign predicate(side2)( point(p0), point(p1), point(p2), point(q0), point(q1) DIM) {scalar l1 = 1*sq_dist(p1,p0) ;scalar l2 = 1*sq_dist(p2,p0) ;scalar a10 = 2*dot_at(p1,q0,p0);scalar a11 = 2*dot_at(p1,q1,p0);scalar a20 = 2*dot_at(p2,q0,p0);scalar a21 = 2*dot_at(p2,q1,p0);scalar Delta = a11 - a10 ;scalar DeltaLambda0 = a11 - l1 ;scalar DeltaLambda1 = l1 - a10 ;scalar r = Delta*l2 - a20*DeltaLambda0 - a21*DeltaLambda1 ;Sign Delta_sign = sign(Delta) ;Sign r_sign = sign(r) ;generic_predicate_result(Delta_sign*r_sign) ;begin_sos3(p0,p1,p2)sos(p0, Sign(Delta_sign*sign(Delta-a21+a20)))sos(p1, Sign(Delta_sign*sign(a21-a20)))sos(p2, NEGATIVE)end_sos} Source PCK
213.Part. 4 Symbolic Perturbation – side2#include "kernel.pckh“Sign predicate(side2)( point(p0), point(p1), point(p2), point(q0), point(q1) DIM) {scalar l1 = 1*sq_dist(p1,p0) ;scalar l2 = 1*sq_dist(p2,p0) ;scalar a10 = 2*dot_at(p1,q0,p0);scalar a11 = 2*dot_at(p1,q1,p0);scalar a20 = 2*dot_at(p2,q0,p0);scalar a21 = 2*dot_at(p2,q1,p0);scalar Delta = a11 - a10 ;scalar DeltaLambda0 = a11 - l1 ;scalar DeltaLambda1 = l1 - a10 ;scalar r = Delta*l2 - a20*DeltaLambda0 - a21*DeltaLambda1 ;Sign Delta_sign = sign(Delta) ;Sign r_sign = sign(r) ;generic_predicate_result(Delta_sign*r_sign) ;begin_sos3(p0,p1,p2)sos(p0, Sign(Delta_sign*sign(Delta-a21+a20)))sos(p1, Sign(Delta_sign*sign(a21-a20)))sos(p2, NEGATIVE)end_sos} Source PCKDenominator
214.Part. 4 Symbolic Perturbation – side2#include "kernel.pckh“Sign predicate(side2)( point(p0), point(p1), point(p2), point(q0), point(q1) DIM) {scalar l1 = 1*sq_dist(p1,p0) ;scalar l2 = 1*sq_dist(p2,p0) ;scalar a10 = 2*dot_at(p1,q0,p0);scalar a11 = 2*dot_at(p1,q1,p0);scalar a20 = 2*dot_at(p2,q0,p0);scalar a21 = 2*dot_at(p2,q1,p0);scalar Delta = a11 - a10 ;scalar DeltaLambda0 = a11 - l1 ;scalar DeltaLambda1 = l1 - a10 ;scalar r = Delta*l2 - a20*DeltaLambda0 - a21*DeltaLambda1 ;Sign Delta_sign = sign(Delta) ;Sign r_sign = sign(r) ;generic_predicate_result(Delta_sign*r_sign) ;begin_sos3(p0,p1,p2)sos(p0, Sign(Delta_sign*sign(Delta-a21+a20)))sos(p1, Sign(Delta_sign*sign(a21-a20)))sos(p2, NEGATIVE)end_sos} Source PCKBarycentriccoords. of q
215.Part. 4 Symbolic Perturbation – side2#include "kernel.pckh“Sign predicate(side2)( point(p0), point(p1), point(p2), point(q0), point(q1) DIM) {scalar l1 = 1*sq_dist(p1,p0) ;scalar l2 = 1*sq_dist(p2,p0) ;scalar a10 = 2*dot_at(p1,q0,p0);scalar a11 = 2*dot_at(p1,q1,p0);scalar a20 = 2*dot_at(p2,q0,p0);scalar a21 = 2*dot_at(p2,q1,p0);scalar Delta = a11 - a10 ;scalar DeltaLambda0 = a11 - l1 ;scalar DeltaLambda1 = l1 - a10 ;scalar r = Delta*l2 - a20*DeltaLambda0 - a21*DeltaLambda1 ;Sign Delta_sign = sign(Delta) ;Sign r_sign = sign(r) ;generic_predicate_result(Delta_sign*r_sign) ;begin_sos3(p0,p1,p2)sos(p0, Sign(Delta_sign*sign(Delta-a21+a20)))sos(p1, Sign(Delta_sign*sign(a21-a20)))sos(p2, NEGATIVE)end_sos} Source PCKBarycentriccoords. of qSolely dependon dot productsbetween inputpoints
216.Part. 4 Symbolic Perturbation – side2#include "kernel.pckh“Sign predicate(side2)( point(p0), point(p1), point(p2), point(q0), point(q1) DIM) {scalar l1 = 1*sq_dist(p1,p0) ;scalar l2 = 1*sq_dist(p2,p0) ;scalar a10 = 2*dot_at(p1,q0,p0);scalar a11 = 2*dot_at(p1,q1,p0);scalar a20 = 2*dot_at(p2,q0,p0);scalar a21 = 2*dot_at(p2,q1,p0);scalar Delta = a11 - a10 ;scalar DeltaLambda0 = a11 - l1 ;scalar DeltaLambda1 = l1 - a10 ;scalar r = Delta*l2 - a20*DeltaLambda0 - a21*DeltaLambda1 ;Sign Delta_sign = sign(Delta) ;Sign r_sign = sign(r) ;generic_predicate_result(Delta_sign*r_sign) ;begin_sos3(p0,p1,p2)sos(p0, Sign(Delta_sign*sign(Delta-a21+a20)))sos(p1, Sign(Delta_sign*sign(a21-a20)))sos(p2, NEGATIVE)end_sos} Source PCKResult when ingeneric position
217.Part. 4 Symbolic Perturbation – side2#include "kernel.pckh“Sign predicate(side2)( point(p0), point(p1), point(p2), point(q0), point(q1) DIM) {scalar l1 = 1*sq_dist(p1,p0) ;scalar l2 = 1*sq_dist(p2,p0) ;scalar a10 = 2*dot_at(p1,q0,p0);scalar a11 = 2*dot_at(p1,q1,p0);scalar a20 = 2*dot_at(p2,q0,p0);scalar a21 = 2*dot_at(p2,q1,p0);scalar Delta = a11 - a10 ;scalar DeltaLambda0 = a11 - l1 ;scalar DeltaLambda1 = l1 - a10 ;scalar r = Delta*l2 - a20*DeltaLambda0 - a21*DeltaLambda1 ;Sign Delta_sign = sign(Delta) ;Sign r_sign = sign(r) ;generic_predicate_result(Delta_sign*r_sign) ;begin_sos3(p0,p1,p2)sos(p0, Sign(Delta_sign*sign(Delta-a21+a20)))sos(p1, Sign(Delta_sign*sign(a21-a20)))sos(p2, NEGATIVE)end_sos} Source PCKSymbolicperturbation
218.Part. 4 Symbolic Perturbation – side2Sign side2_exact_SOS(const double* p0, const double* p1, const double* p2,const double* q0, const double* q1,coord_index_t dim) {const expansion& l1 = expansion_sq_dist(p1, p0, dim);const expansion& l2 = expansion_sq_dist(p2, p0, dim);const expansion& a10 = expansion_dot_at(p1, q0, p0, dim).scale_fast(2.0);const expansion& a11 = expansion_dot_at(p1, q1, p0, dim).scale_fast(2.0);const expansion& a20 = expansion_dot_at(p2, q0, p0, dim).scale_fast(2.0);const expansion& a21 = expansion_dot_at(p2, q1, p0, dim).scale_fast(2.0);const expansion& Delta = expansion_diff(a11, a10);Sign Delta_sign = Delta.sign();vor_assert(Delta_sign != ZERO);const expansion& DeltaLambda0 = expansion_diff(a11, l1);const expansion& DeltaLambda1 = expansion_diff(l1, a10);const expansion& r0 = expansion_product(Delta, l2);const expansion& r1 = expansion_product(a20, DeltaLambda0).negate();const expansion& r2 = expansion_product(a21, DeltaLambda1).negate();const expansion& r = expansion_sum3(r0, r1, r2);Sign r_sign = r.sign();………..Exact version with expansions
219.Part. 4 Symbolic Perturbation – side2if(r_sign == ZERO) {const double* p_sort[3];p_sort[0] = p0;p_sort[1] = p1;p_sort[2] = p2;std::sort(p_sort, p_sort + 3);for(index_t i = 0; i < 3; ++i) {if(p_sort[i] == p0) {const expansion& z1 = expansion_diff(Delta, a21);const expansion& z = expansion_sum(z1, a20);Sign z_sign = z.sign();len_side2_SOS = vor_max(len_side2_SOS, z.length());if(z_sign != ZERO) {return Sign(Delta_sign * z_sign);}}if(p_sort[i] == p1) {const expansion& z = expansion_diff(a21, a20);Sign z_sign = z.sign();len_side2_SOS = vor_max(len_side2_SOS, z.length());if(z_sign != ZERO) {return Sign(Delta_sign * z_sign);}}if(p_sort[i] == p2) {return NEGATIVE;}}vor_assert_not_reached;}return Sign(Delta_sign * r_sign);} Exact version with expansions
220.Part. 4 Symbolic Perturbation – side3#include "kernel.pckh"Sign predicate(side3)(point(p0), point(p1), point(p2), point(p3),point(q0), point(q1), point(q2) DIM) {scalar l1 = 1*sq_dist(p1,p0);scalar l2 = 1*sq_dist(p2,p0);scalar l3 = 1*sq_dist(p3,p0);scalar a10 = 2*dot_at(p1,q0,p0);scalar a11 = 2*dot_at(p1,q1,p0);scalar a12 = 2*dot_at(p1,q2,p0);scalar a20 = 2*dot_at(p2,q0,p0);scalar a21 = 2*dot_at(p2,q1,p0);scalar a22 = 2*dot_at(p2,q2,p0);scalar a30 = 2*dot_at(p3,q0,p0);scalar a31 = 2*dot_at(p3,q1,p0);scalar a32 = 2*dot_at(p3,q2,p0);scalar b00 = a11*a22-a12*a21;scalar b01 = a21-a22;scalar b02 = a12-a11;scalar b10 = a12*a20-a10*a22;scalar b11 = a22-a20;scalar b12 = a10-a12;scalar b20 = a10*a21-a11*a20;scalar b21 = a20-a21;scalar b22 = a11-a10;scalar Delta = b00+b10+b20;scalar DeltaLambda0 =b01*l1+b02*l2+b00 ;scalar DeltaLambda1 =b11*l1+b12*l2+b10 ;scalar DeltaLambda2 =b21*l1+b22*l2+b20 ;scalar r = Delta*l3 - (a30 * DeltaLambda0 +a31 * DeltaLambda1 +a32 * DeltaLambda2) ;Sign Delta_sign = sign(Delta) ;Sign r_sign = sign(r) ;generic_predicate_result(Delta_sign*r_sign) ;begin_sos4(p0,p1,p2,p3)sos(p0, Sign(Delta_sign*sign(Delta-((b01+b02)*a30+(b11+b12)*a31+(b21+b22)*a32))))sos(p1, Sign(Delta_sign*sign((a30*b01)+(a31*b11)+(a32*b21))))sos(p2, Sign(Delta_sign*sign((a30*b02)+(a31*b12)+(a32*b22))))sos(p3, NEGATIVE)end_sos}Source PCK
221.Part. 4 Symbolic Perturbation – side4……….scalar b00= det3x3(a11,a12,a13,a21,a22,a23,a31,a32,a33);scalar b01=-det_111_2x3(a21,a22,a23,a31,a32,a33);scalar b02= det_111_2x3(a11,a12,a13,a31,a32,a33);scalar b03=-det_111_2x3(a11,a12,a13,a21,a22,a23);scalar b10=-det3x3(a10,a12,a13,a20,a22,a23,a30,a32,a33);scalar b11= det_111_2x3(a20,a22,a23,a30,a32,a33);scalar b12=-det_111_2x3(a10,a12,a13,a30,a32,a33);scalar b13= det_111_2x3(a10,a12,a13,a20,a22,a23);scalar b20= det3x3(a10,a11,a13,a20,a21,a23,a30,a31,a33);scalar b21=-det_111_2x3(a20,a21,a23,a30,a31,a33);scalar b22= det_111_2x3(a10,a11,a13,a30,a31,a33);scalar b23=-det_111_2x3(a10,a11,a13,a20,a21,a23);scalar b30=-det3x3(a10,a11,a12,a20,a21,a22,a30,a31,a32);scalar b31= det_111_2x3(a20,a21,a22,a30,a31,a32);scalar b32=-det_111_2x3(a10,a11,a12,a30,a31,a32);scalar b33= det_111_2x3(a10,a11,a12,a20,a21,a22);scalar Delta=b00+b10+b20+b30;scalar DeltaLambda0 = b01*l1+b02*l2+b03*l3+b00;scalar DeltaLambda1 = b11*l1+b12*l2+b13*l3+b10;scalar DeltaLambda2 = b21*l1+b22*l2+b23*l3+b20;scalar DeltaLambda3 = b31*l1+b32*l2+b33*l3+b30;……….scalar r = Delta*l4 - (a40*DeltaLambda0+a41*DeltaLambda1+a42*DeltaLambda2+a43*DeltaLambda3);Sign Delta_sign = sign(Delta);generic_predicate_result(Delta_sign*sign(r)) ;begin_sos5(p0,p1,p2,p3,p4)sos(p0, Sign( Delta_sign*sign(Delta - ((b01+b02+b03)*a30 +(b11+b12+b13)*a31 +(b21+b22+b23)*a32 +(b31+b32+b33)*a33))))sos(p1, Sign( Delta_sign*sign(a30*b01+a31*b11+a32*b21+a33*b31)))sos(p2, Sign( Delta_sign*sign(a30*b02+a31*b12+a32*b22+a33*b32)))sos(p3, Sign( Delta_sign*sign(a30*b03+a31*b13+a32*b23+a33*b33)))sos(p4, NEGATIVE)end_sos}q is the intersection of three bisectors and a tetrahedron embedded in nDNote: There is a special case in 3d (no need for the tetrahedron) = insphere3d()The code of the generalnD version (excerpt) :Source PCK
222.Part 4. Geometry – Crash test 1/2
223.Part 4. Geometry – Crash test 1/2
224.Part 4. Geometry – Crash test 1/2
225.Part 4. Geometry – Crash test 1/2
226.Part 4. Geometry – Crash test 2/2
227.Part 4. Geometry – Crash test 2/2
228.Part 4. Geometry - RVD
229.Part 4. Geometry - RVD
230.Part 4. Geometry - RVD
231.Part 4. Geometry - RVD
232.Part 4. Geometry - RVDIf we eliminate the zero components during computations, the length of theexpansions remain reasonable (see observation in Shewchuk’s paper)
233.Part 4. Clipped Voronoi diagram
234.Part 4. Clipped Voronoi diagramPolyhedral elements forFinite Elements Analysis>vorpaline/vorpalite profile=poly fusee.meshb
235.Part 4. Geometry – 3D Anisotropic VD
236.Part 4. Geometry – 3D Anisotropic VD
237.The PCK (Predicate Construction Kit)Features:• Expansion number type – low level API (allocation on stack, efficient)• High level API with operators (easy to use, less efficient)• Script to generate FPG filter and exact version with SOS• Standard predicates (orient2d, 3d, insphere3d, orient4d)• More exotic predicates for RVD (side1, side2, side3, side 4 in dim 3,4,6,7)• 3D Delaunay triangulation• 3D weighted Delaunay triangulation• Only 6 source files• multi_precision.(h,cpp)• predicates.(h,cpp)• delaunay3d.(h,cpp)• Fully documented• No dependency (compiles and runs everywhere*, I got a version on my phone :-)BSD license (do what you want with it, including business)*In the IEEE754 world
238.What s next5Physics + Mathematics + Computing = ?
239.?Tρ1 ρ2Minimize C(T) =∫Ω|| x – T(x) ||2 dxSubject to: T is measure-preservingρ1(x)Optimal Transport
240.A map T is a transport map between μ and ν ifμ(T-1(B)) = ν(B) for any Borel subset BB(X;μ) (Y;ν)Optimal Transport
241.A map T is a transport map between μ and ν ifμ(T-1(B)) = ν(B) for any Borel subset BBT-1(B)(X;μ) (Y;ν)Optimal Transport
242.Continuous(X;μ) (Y;ν)Optimal Transport : probability measures
243.ContinuousSemi-discreteDiscrete(X;μ) (Y;ν)Optimal Transport : probability measures
244.ContinuousSemi-discreteDiscrete(X;μ) (Y;ν)Optimal Transport : probability measures
245.Optimal Transport – semi-discrete∫X ψc (x)dμ + ∫Y ψ(y)dνSupψ Є ψc(DMK)∑j ∫Lagψ(yj) || x – yj ||2 - ψ(yj) dμ∫X inf yj ЄY [ || x – yj ||2 - ψ(yj) ] dμ∑j ψ(yj) vj
246.Voronoi diagram: Vor(xi) = { x | d2(x,xi) < d2(x,xj) }Power Diagrams
247.Power diagram: Pow(xi) = { x | d2(x,xi) – ψi < d2(x,xj) – ψj }Voronoi diagram: Vor(xi) = { x | d2(x,xi) < d2(x,xj) }Power Diagrams
248.Shape Interpolation
249.Shape Interpolation
250.Shape Fitting
251.Shape FittingTopOpt2FEA<<<<<<
252.Scan2FEA3D scan from ARTEC
253.Remesh with CVTScan2FEA
254.Periodic Global ParameterizationScan2FEA
255.Subd control meshScan2FEA
256.Subd surface fitted using Optimal TransportScan2FEA
257.Finite Element Analysis (vibration modes)Scan2FEA
258.Think Big
259.Think Big: Simulating the entire universe
260.The millenium simulation project, Max Planck Institute fur Astrophysikpc/h : parsec (= 3.2 années lumières)Large-Scale structure
261.Optimal TransportEarly Universe ReconstructionThe Data
262.Optimal TransportInvert Newton Einstein Eqn to go back 14 milliards years in timeThe millenium simulation project,Max Planck Institute fur Astrophysikpc/h : parsec (= 3.2 light years)
263.Optimal TransportThe millenium simulation project,Max Planck Institute fur Astrophysikpc/h : parsec (= 3.2 années lumières)In 2002, 5 hours of computationwith a supercomputer / 5000 points
264.Optimal TransportThe millenium simulation project,Max Planck Institute fur Astrophysikpc/h : parsec (= 3.2 années lumières)In 2002, 5 hours of computationwith a supercomputer / 5000 pointsCan we do it with 1 000 000 points ?
265.Optimal TransportThe millenium simulation project,Max Planck Institute fur Astrophysikpc/h : parsec (= 3.2 années lumières)In 2002, 5 hours of computationwith a supercomputer / 5000 pointsCan we do it with 1 000 000 points ?Yes, if we wait 4500 years
266.Optimal TransportIn 2002, 5 hours of computationwith a supercomputer / 5000 pointsCan we do it with 1 000 000 points ?Yes, if we wait 4500 years
267.Optimal TransportIn 2002, 5 hours of computationwith a supercomputer / 5000 pointsCan we do it with 1 000 000 points ?Yes, if we wait 4500 yearsWe need a new algorithm.
268.Towards Early Universe ReconstructionNumerical Experiment: Performances2002: several hours of supercomputer time were neededfor computing OT with a few thousand Dirac masses, with a combinatorialalgorithm in O(n2log(n))
269.Towards Early Universe ReconstructionNumerical Experiment: Performances2002: several hours of supercomputer time were neededfor computing OT with a few thousand Dirac masses, with a combinatorialalgorithm in O(n2log(n))2015:3D version of [Mérigot] (multilevel + BFGS) + several tricks [L 2015]
270.Towards Early Universe ReconstructionNumerical Experiment: Performances2002: several hours of supercomputer time were neededfor computing OT with a few thousand Dirac masses, with a combinatorialalgorithm in O(n2log(n))2015:3D version of [Mérigot] (multilevel + BFGS) + several tricks [L 2015]
271.Towards Early Universe ReconstructionNumerical Experiment: Performances2002: several hours of supercomputer time were neededfor computing OT with a few thousand Dirac masses, with a combinatorialalgorithm in O(n2log(n))2015:3D version of [Mérigot] (multilevel + BFGS) + several tricks [L 2015]2016: Damped Newton [Mérigot, Thibert] + several tricks for 3D:1 million Dirac masses in 240 seconds
272.Towards Early Universe ReconstructionNumerical Experiment: Performances2002: several hours of supercomputer time were neededfor computing OT with a few thousand Dirac masses, with a combinatorialalgorithm in O(n2log(n))2015:3D version of [Mérigot] (multilevel + BFGS) + several tricks [L 2015]2016: Damped Newton [Mérigot, Thibert] + several tricks for 3D:1 million Dirac masses in 240 seconds10 million Dirac masses in 90 minutes
273.Towards Early Universe ReconstructionNumerical Experiment: Performances2002: several hours of supercomputer time were neededfor computing OT with a few thousand Dirac masses, with a combinatorialalgorithm in O(n2log(n))2015:3D version of [Mérigot] (multilevel + BFGS) + several tricks [L 2015]2016: Damped Newton [Mérigot, Thibert] + several tricks for 3D:1 million Dirac masses in 240 seconds10 million Dirac masses in 90 minutes2017: 10 million Dirac masses in 2 minutes (for specific configurations)Semi-discrete OT is now scalable ! (new tool in Num. Ana. Toolbox)
274.Some concluding words
275.Take-home messageProgramming is a great source of fun
276.Take-home messageProgramming is a great source of funFuturistic programming principles(make it fun for your users as well)
277.Take-home messageProgramming is a great source of funFuturistic programming principlesProgram speedLow memory consumptionLines of codeClassesTemplates
278.Take-home messageProgramming is a great source of funFuturistic programming principlesProgram speedLow memory consumptionLines of codeClassesTemplatesGains Costs
279.Take-home messageProgramming is a great source of funFuturistic programming principlesProgram speedLow memory consumptionLines of codeClassesTemplatesGains CostsMeasure it !(profiler, continuous integration)
280.Take-home message: Usefulness is primary.h, .cpp.h, .cpp.h, .cpp.h, .cpp.h, .cpp.h, .cpp.h, .cppObject-oriented designGeneric design
281.Take-home message: Usefulness is primary.h, .cpp.h, .cpp.h, .cpp.h, .cpp.h, .cpp.h, .cpp.h, .cppObject-oriented designGeneric design
282.Take-home message: Usefulness is primaryFuturistic design.cpp.hvoid the_functionality(Data& )
283.Take-home message: Usefulness is primaryFuturistic design.cpp.hvoid the_functionality(Data& )* Make it easy to compile* Do not bother the userwith internal details (even ifyou are proud of them)
284.Take-home message: Usefulness is primary.cppGLUP.h – C APIglupBegin(), glupEnd(), glupVertex()GLUP statevariablesVanillaGLOpenGL ES /WebGLGLSL 1.5GLSL 4.4matrix stacksbuffersShadersOpenNL.h – C APInlBegin(), nlEnd(), nlCoeff(), nlSolve().cCGGMResBiCGStabSuperLUCHOLDMODCPU/GPUAbstraction layerOpenMPmulticoreCUDACuBLASCuSparseARPACKNLSparseMatrix NLCRSMatrixnlCompress()
285.Take-home message: Usefulness is primary.cppGLUP.h – C APIglupBegin(), glupEnd(), glupVertex()GLUP statevariablesVanillaGLOpenGL ES /WebGLGLSL 1.5GLSL 4.4matrix stacksbuffersShadersOpenNL.h – C APInlBegin(), nlEnd(), nlCoeff(), nlSolve().cppCGGMResBiCGStabSuperLUCHOLDMODCPU/GPUAbstraction layerOpenMPmulticoreCUDACuBLASCuSparseARPACKNLSparseMatrix NLCRSMatrixnlCompress()Algorithms and Data StructuresMostly arrays (std::vector) and for() loopsObject oriented / virtual functionsRun-time selection of algorithmGeneric programmingFor dimension-independent code (RVD)If() statementsPredicates
286.Futuristic programming – the productGLUP.h – C APIglupBegin(), glupEnd(), glupVertex()OpenNL.h – C APInlBegin(), nlEnd(), nlCoeff(), nlSolve()Geogram – C++ APIMesh class, Delaunay, Voronoi,Remesh, param., repair, reconstruct.Predicate Construction KitCompiler for arbitrary precisionGeometric predicates (C/C++ API)Graphite: Qt GUI + LUA scriptingDownload from gforge.inria.fr, see also links on my webpage
287.AcknowledgementsEurographics AssociationEuropean Research CouncilGOODSHAPE ERC-StG-205693VORPALINE ERC-PoC-334829ANR MORPHO (Vision), ANR BECASIM (Physics)ANR MAGA (P.I.: Q. Merigot),ANR ROOT (P.I.: N. Bonneel) Optimal Transport stuff.Inria EXPLORAGRAMKey contributors to Geogram/Graphite:N. Ray, W.-C. Li, B. Vallet, D. Sokolov, N. Bonneel, R. Zayer, A. ShefferMOKA team (Monge-Ampere-Kantorovich)