3910 lines
201 KiB
HTML
3910 lines
201 KiB
HTML
<!DOCTYPE html>
|
||
<html>
|
||
<head>
|
||
<meta charset="utf-8"/>
|
||
<title>libigl Tutorial</title>
|
||
<meta name="author" content="Daniele Panozzo and Alec Jacobson"/>
|
||
<meta name="date" content="07 November 2015"/>
|
||
<link type="text/css" rel="stylesheet" href="style.css"/>
|
||
<script type="text/javascript" src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML"></script>
|
||
<link rel="stylesheet" href="http://yandex.st/highlightjs/7.3/styles/default.min.css">
|
||
<script src="http://yandex.st/highlightjs/7.3/highlight.min.js"></script>
|
||
<script>hljs.initHighlightingOnLoad();</script>
|
||
</head>
|
||
<body>
|
||
|
||
<h1 id="libigltutorialnotes">libigl tutorial notes</h1>
|
||
|
||
<h4 id="originallypresentedbydanielepanozzoandalecjacobsonatsgpgraduateschool2014">originally presented by Daniele Panozzo and Alec Jacobson at SGP Graduate School 2014</h4>
|
||
|
||
<figure>
|
||
<img src="images/libigl-logo.jpg" alt="" />
|
||
</figure>
|
||
|
||
<p>Libigl is an open source C++ library for geometry processing research and
|
||
development. Dropping the heavy data structures of tradition geometry
|
||
libraries, libigl is a simple header-only library of encapsulated functions.
|
||
This combines the rapid prototyping familiar to Matlab or Python programmers
|
||
with the performance and versatility of C++. The tutorial is a self-contained,
|
||
hands-on introduction to libigl. Via interactive, step-by-step examples, we
|
||
demonstrate how to accomplish common geometry processing tasks such as
|
||
computation of differential quantities and operators, real-time deformation,
|
||
parametrization, numerical optimization and remeshing. Each section of the
|
||
lecture notes links to a cross-platform example application.</p>
|
||
|
||
<h1 id="tableofcontents">Table of contents</h1>
|
||
|
||
<ul>
|
||
<li><a href="#chapter1:introductiontolibigl">Chapter 1: Introduction to libigl</a>
|
||
|
||
<ul>
|
||
<li><a href="#libigldesignprinciples">Libigl design principles</a></li>
|
||
<li><a href="#meshrepresentation">101 Mesh representation</a></li>
|
||
<li><a href="#visualizingsurfaces">102 Visualizing surfaces</a></li>
|
||
<li><a href="#interactionwithkeyboardandmouse">103 Interaction with keyboard and mouse</a></li>
|
||
<li><a href="#scalarfieldvisualization">104 Scalar field visualization</a></li>
|
||
<li><a href="#overlays">105 Overlays</a></li>
|
||
<li><a href="#viewermenu">106 Viewer Menu</a></li>
|
||
</ul></li>
|
||
<li><a href="#chapter2:discretegeometricquantitiesandoperators">Chapter 2: Discrete Geometric Quantities and
|
||
Operators</a>
|
||
|
||
<ul>
|
||
<li><a href="#normals">201 Normals</a>
|
||
|
||
<ul>
|
||
<li><a href="#per-face">Per-face</a></li>
|
||
<li><a href="#per-vertex">Per-vertex</a></li>
|
||
<li><a href="#per-corner">Per-corner</a></li>
|
||
</ul></li>
|
||
<li><a href="#gaussiancurvature">202 Gaussian Curvature</a></li>
|
||
<li><a href="#curvaturedirections">203 Curvature Directions</a></li>
|
||
<li><a href="#gradient">204 Gradient</a></li>
|
||
<li><a href="#laplacian">204 Laplacian</a>
|
||
|
||
<ul>
|
||
<li><a href="#massmatrix">Mass matrix</a></li>
|
||
<li><a href="#alternativeconstructionoflaplacian">Alternative construction of
|
||
Laplacian</a></li>
|
||
</ul></li>
|
||
</ul></li>
|
||
<li><a href="#chapter3:matricesandlinearalgebra">Chapter 3: Matrices and Linear Algebra</a>
|
||
|
||
<ul>
|
||
<li><a href="#slice">301 Slice</a></li>
|
||
<li><a href="#sort">302 Sort</a>
|
||
|
||
<ul>
|
||
<li><a href="#othermatlab-stylefunctions">Other Matlab-style functions</a></li>
|
||
</ul></li>
|
||
<li><a href="#laplaceequation">303 Laplace Equation</a>
|
||
|
||
<ul>
|
||
<li><a href="#quadraticenergyminimization">Quadratic energy minimization</a></li>
|
||
</ul></li>
|
||
<li><a href="#linearequalityconstraints">304 Linear Equality Constraints</a></li>
|
||
<li><a href="#quadraticprogramming">305 Quadratic Programming</a></li>
|
||
<li><a href="#eigendecomposition">306 Eigen Decomposition</a></li>
|
||
</ul></li>
|
||
<li><a href="#chapter4:shapedeformation">Chapter 4: Shape Deformation</a>
|
||
|
||
<ul>
|
||
<li><a href="#biharmonicdeformation">401 Biharmonic Deformation</a></li>
|
||
<li><a href="#polyharmonicdeformation">402 Polyharmonic Deformation</a></li>
|
||
<li><a href="#boundedbiharmonicweights">403 Bounded Biharmonic Weights</a></li>
|
||
<li><a href="#dualquaternionskinning">404 Dual Quaternion Skinning</a></li>
|
||
<li><a href="#as-rigid-as-possible">405 As-rigid-as-possible</a></li>
|
||
<li><a href="#fastautomaticskinningtransformations">406 Fast automatic skinning
|
||
transformations</a>
|
||
|
||
<ul>
|
||
<li><a href="#arapwithgroupededge-sets">ARAP with grouped edge-sets</a></li>
|
||
</ul></li>
|
||
<li><a href="#biharmoniccoordinates">407 Biharmonic Coordinates</a></li>
|
||
</ul></li>
|
||
<li><a href="#chapter5:parametrization">Chapter 5: Parametrization</a>
|
||
|
||
<ul>
|
||
<li><a href="#harmonicparametrization">501 Harmonic parametrization</a></li>
|
||
<li><a href="#leastsquareconformalmaps">502 Least-Square Conformal Maps</a></li>
|
||
<li><a href="#asrigidaspossible">503 As-Rigid-As-Possible</a></li>
|
||
<li><a href="#nrotationallysymmetrictangetfields">504 N-Rotationally symmetric tangent fields</a></li>
|
||
<li><a href="#globalseamlessintegergridparametrization">505 Global, seamless integer-grid parametrization</a></li>
|
||
<li><a href="#anisotropicremeshingusingframefields">506 Anisotropic remeshing using frame fields</a></li>
|
||
<li><a href="#npolyvectorfields">507 N-PolyVector fields</a></li>
|
||
<li><a href="#conjugatevectorfields">508 Conjugate vector fields</a></li>
|
||
<li><a href="#planarization">509 Planarization</a></li>
|
||
<li><a href="#integrable">510 Integrable PolyVector Fields</a></li>
|
||
<li><a href="#npolyvectorfields_general">511 General N-PolyVector Fields</a></li>
|
||
</ul></li>
|
||
<li><a href="#chapter6:externallibraries">Chapter 6: External libraries</a>
|
||
|
||
<ul>
|
||
<li><a href="#stateserialization">601 State serialization</a></li>
|
||
<li><a href="#mixingmatlabcode">602 Mixing Matlab code</a>
|
||
|
||
<ul>
|
||
<li><a href="#savingamatlabworkspace">Saving a Matlab workspace</a></li>
|
||
<li><a href="#dumpingeigenmatricestocopyandpasteintomatlab">Dumping Eigen matrices to copy and paste into
|
||
Matlab</a></li>
|
||
</ul></li>
|
||
<li><a href="#callinglibiglfunctionsfrommatlab">603 Calling libigl functions from Matlab</a></li>
|
||
<li><a href="#triangulationofclosedpolygons">604 Triangulation of closed polygons</a></li>
|
||
<li><a href="#tetrahedralizationofclosedsurfaces">605 Tetrahedralization of closed surfaces</a></li>
|
||
<li><a href="#bakingambientocclusion">606 Baking ambient occlusion</a></li>
|
||
<li><a href="#screencapture">607 Screen Capture</a></li>
|
||
<li><a href="#locallyinjectivemaps">608 Locally Injective Maps</a></li>
|
||
<li><a href="#booleanoperationsonmeshes">609 Boolean Operations on Meshes</a></li>
|
||
<li><a href="#csgtree">610 CSG Tree</a></li>
|
||
</ul></li>
|
||
<li><a href="#chapter7:miscellaneous">Chapter 7: Miscellaneous</a>
|
||
|
||
<ul>
|
||
<li><a href="#meshstatistics">701 Mesh Statistics</a></li>
|
||
<li><a href="#generalizedwindingnumber">702 Generalized Winding Number</a></li>
|
||
<li><a href="#meshdecimation">703 Mesh Decimation</a></li>
|
||
<li><a href="#signeddistances">704 Signed Distances</a></li>
|
||
<li><a href="#marchingcubes">705 Marching Cubes</a></li>
|
||
<li><a href="#facetorientation">706 Facet Orientation</a></li>
|
||
<li><a href="#sweptvolume">707 Swept Volume</a></li>
|
||
<li><a href="#pickingverticesandfaces">708 Picking Vertices and Faces</a></li>
|
||
<li><a href="#vectorfieldvisualizer">709 Vector Field Visualization</a></li>
|
||
<li><a href="#slim">710 Scalable Locally Injective Maps</a></li>
|
||
<li><a href="#subdivision">711 Subdivision surfaces</a></li>
|
||
</ul></li>
|
||
<li><a href="#future">Chapter 8: Outlook for continuing development</a></li>
|
||
</ul>
|
||
|
||
<h1 id="chapter1:introductiontolibigl">Chapter 1</h1>
|
||
|
||
<p>We introduce libigl with a series of self-contained examples. The purpose of
|
||
each example is to showcase a feature of libigl while applying to a practical
|
||
problem in geometry processing. In this chapter, we will present the basic
|
||
concepts of libigl and introduce a simple mesh viewer that allows to
|
||
visualize a surface mesh and its attributes. All the tutorial examples are
|
||
cross-platform and can be compiled on MacOSX, Linux and Windows.</p>
|
||
|
||
<h2 id="libigldesignprinciples"><a href="#libigldesignprinciples">libigl design principles</a></h2>
|
||
|
||
<p>Before getting into the examples, we summarize the main design principles in
|
||
libigl:</p>
|
||
|
||
<ol>
|
||
<li><p><strong>No complex data types.</strong> We mostly use matrices and vectors. This greatly
|
||
favors code reusability and forces the function authors to expose all the
|
||
parameters used by the algorithm. </p></li>
|
||
<li><p><strong>Minimal dependencies.</strong> We use external libraries only when necessary and
|
||
we wrap them in a small set of functions.</p></li>
|
||
<li><p><strong>Header-only.</strong> It is straight forward to use our library since it is only
|
||
one additional include directory in your project. (if you are worried about
|
||
compilation speed, it is also possible to build the library as a <a href="../optional/">static
|
||
library</a>)</p></li>
|
||
<li><p><strong>Function encapsulation.</strong> Every function (including its full
|
||
implementation) is contained in a pair of .h/.cpp files with the same name of
|
||
the function.</p></li>
|
||
</ol>
|
||
|
||
<h3 id="downloadinglibigl">Downloading libigl</h3>
|
||
|
||
<p>libigl can be downloaded from our <a href="https://github.com/libigl/libigl">github
|
||
repository</a> or cloned with git:</p>
|
||
|
||
<pre><code class="bash">git clone --recursive https://github.com/libigl/libigl.git
|
||
</code></pre>
|
||
|
||
<p>The core libigl functionality only depends on the C++ Standard Library and
|
||
Eigen.</p>
|
||
|
||
<p>To build all the examples in the tutorial, you can use the CMakeLists.txt in
|
||
the tutorial folder:</p>
|
||
|
||
<pre><code class="bash">cd tutorial
|
||
mkdir build
|
||
cd build
|
||
cmake -DCMAKE_BUILD_TYPE=Release ../
|
||
make
|
||
</code></pre>
|
||
|
||
<p>The examples can also be built independently using the CMakeLists.txt
|
||
inside each example folder.</p>
|
||
|
||
<p><em>Note for linux users</em>: Many linux distributions do not include gcc and the basic development tools
|
||
in their default installation. On Ubuntu, you need to install the following packages:</p>
|
||
|
||
<pre><code class="bash">sudo apt-get install git
|
||
sudo apt-get install build-essential
|
||
sudo apt-get install cmake
|
||
sudo apt-get install libx11-dev
|
||
sudo apt-get install mesa-common-dev libgl1-mesa-dev libglu1-mesa-dev
|
||
sudo apt-get install libxrandr-dev
|
||
sudo apt-get install libxi-dev
|
||
sudo apt-get install libxmu-dev
|
||
sudo apt-get install libblas-dev
|
||
sudo apt-get install libxinerama-dev
|
||
sudo apt-get install libxcursor-dev
|
||
</code></pre>
|
||
|
||
<p><em>Note for windows users</em>: libigl only supports the Microsoft Visual Studio 2015 compiler in 64bit mode. It will not work with a 32bit build and it will not work
|
||
with older versions of visual studio.</p>
|
||
|
||
<p>A few examples in Chapter 5 requires the <a href="http://www.graphics.rwth-aachen.de/software/comiso">CoMiSo
|
||
solver</a>. We provide a
|
||
mirror of CoMISo that works out of the box with libigl. To install it:</p>
|
||
|
||
<pre><code class="bash">cd libigl/external
|
||
git clone --recursive https://github.com/libigl/CoMISo.git
|
||
</code></pre>
|
||
|
||
<p>You can then build the tutorials again and it libigl will automatically find and
|
||
compile CoMISo.</p>
|
||
|
||
<p><em>Note 1</em>: CoMISo is distributed under the GPL3 license, it does impose restrictions on commercial usage.</p>
|
||
|
||
<p><em>Note 2</em>: CoMISo requires a blas implementation. We use the built-in blas in macosx and linux, and we bundle a precompiled binary for VS2015 64 bit. Do NOT compile the tutorials
|
||
in 32 bit on windows.</p>
|
||
|
||
<h3 id="libiglexampleproject">libigl example project</h3>
|
||
|
||
<p>We provide a <a href="https://github.com/libigl/libigl-example-project">blank project example</a> showing how to use libigl and cmake. Feel free and encouraged to copy or fork this project as a way of starting a new personal project using libigl.</p>
|
||
|
||
<h2 id="meshrepresentation"><a href="#meshrepresentation">Mesh representation</a></h2>
|
||
|
||
<p>libigl uses the <a href="http://eigen.tuxfamily.org/">Eigen</a> library to encode vector
|
||
and matrices. We suggest that you keep the
|
||
<a href="http://eigen.tuxfamily.org/dox/group__QuickRefPage.html">dense</a> and
|
||
<a href="http://eigen.tuxfamily.org/dox/group__SparseQuickRefPage.html">sparse</a> quick
|
||
reference guides at hand while you read the examples in this tutorial.</p>
|
||
|
||
<p>A triangular mesh is encoded as a pair of matrices:</p>
|
||
|
||
<pre><code class="cpp">Eigen::MatrixXd V;
|
||
Eigen::MatrixXi F;
|
||
</code></pre>
|
||
|
||
<p><code>V</code> is a #N by 3 matrix which stores the coordinates of the vertices. Each
|
||
row stores the coordinate of a vertex, with its x,y and z coordinates in the first,
|
||
second and third column, respectively. The matrix <code>F</code> stores the triangle
|
||
connectivity: each line of <code>F</code> denotes a triangle whose 3 vertices are
|
||
represented as indices pointing to rows of <code>V</code>.</p>
|
||
|
||
<figure>
|
||
<img src="images/VF.png" alt="A simple mesh made of 2 triangles and 4 vertices." />
|
||
<figcaption>A simple mesh made of 2 triangles and 4 vertices.</figcaption>
|
||
</figure>
|
||
|
||
<p>Note that the order of the vertex indices in <code>F</code> determines the orientation of
|
||
the triangles and it should thus be consistent for the entire surface.
|
||
This simple representation has many advantages:</p>
|
||
|
||
<ol>
|
||
<li>it is memory efficient and cache friendly</li>
|
||
<li>the use of indices instead of pointers greatly simplifies debugging</li>
|
||
<li>the data can be trivially copied and serialized</li>
|
||
</ol>
|
||
|
||
<p>libigl provides input [output] functions to read [write] many common mesh formats.
|
||
The IO functions are contained in the files read*.h and write*.h. As a general
|
||
rule each libigl function is contained in a pair of .h/.cpp files with the same name.
|
||
By default, the .h files include the corresponding cpp files, making the library header-only.</p>
|
||
|
||
<p>Reading a mesh from a file requires a single libigl function call:</p>
|
||
|
||
<pre><code class="cpp">igl::readOFF(TUTORIAL_SHARED_PATH "/cube.off", V, F);
|
||
</code></pre>
|
||
|
||
<p>The function reads the mesh cube.off and it fills the provided <code>V</code> and <code>F</code> matrices.
|
||
Similarly, a mesh can be written in an OBJ file using:</p>
|
||
|
||
<pre><code class="cpp">igl::writeOBJ("cube.obj",V,F);
|
||
</code></pre>
|
||
|
||
<p><a href="101_FileIO/main.cpp">Example 101</a> contains a simple mesh
|
||
converter from OFF to OBJ format.</p>
|
||
|
||
<h2 id="visualizingsurfaces"><a href="#visualizingsurfaces">Visualizing surfaces</a></h2>
|
||
|
||
<p>Libigl provides an glfw-based OpenGL 3.2 viewer to visualize surfaces, their
|
||
properties and additional debugging informations.</p>
|
||
|
||
<p>The following code (<a href="102_DrawMesh/main.cpp">Example 102</a>) is a basic skeleton
|
||
for all the examples that will be used in the tutorial.
|
||
It is a standalone application that loads a mesh and uses the viewer to
|
||
render it.</p>
|
||
|
||
<pre><code class="cpp">#include <igl/readOFF.h>
|
||
#include <igl/viewer/Viewer.h>
|
||
|
||
Eigen::MatrixXd V;
|
||
Eigen::MatrixXi F;
|
||
|
||
int main(int argc, char *argv[])
|
||
{
|
||
// Load a mesh in OFF format
|
||
igl::readOFF(TUTORIAL_SHARED_PATH "/bunny.off", V, F);
|
||
|
||
// Plot the mesh
|
||
igl::viewer::Viewer viewer;
|
||
viewer.data.set_mesh(V, F);
|
||
viewer.launch();
|
||
}
|
||
</code></pre>
|
||
|
||
<p>The function <code>set_mesh</code> copies the mesh into the viewer.
|
||
<code>Viewer.launch()</code> creates a window, an OpenGL context and it starts the draw loop.
|
||
Additional properties can be plotted on the mesh (as we will see later),
|
||
and it is possible to extend the viewer with standard OpenGL code.
|
||
Please see the documentation in
|
||
<a href="../include/igl/Viewer/Viewer.h">Viewer.h</a> for more details.</p>
|
||
|
||
<figure>
|
||
<img src="images/102_DrawMesh.png" alt="(Example 102) loads and draws a
|
||
mesh." />
|
||
<figcaption>(<a href="102_DrawMesh/main.cpp">Example 102</a>) loads and draws a
|
||
mesh.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="interactionwithkeyboardandmouse"><a href="#interactionwithkeyboardandmouse">Interaction with keyboard and mouse</a></h2>
|
||
|
||
<p>Keyboard and mouse events triggers callbacks that can be registered in the
|
||
viewer. The viewer supports the following callbacks:</p>
|
||
|
||
<pre><code class="cpp">bool (*callback_pre_draw)(Viewer& viewer);
|
||
bool (*callback_post_draw)(Viewer& viewer);
|
||
bool (*callback_mouse_down)(Viewer& viewer, int button, int modifier);
|
||
bool (*callback_mouse_up)(Viewer& viewer, int button, int modifier);
|
||
bool (*callback_mouse_move)(Viewer& viewer, int mouse_x, int mouse_y);
|
||
bool (*callback_mouse_scroll)(Viewer& viewer, float delta_y);
|
||
bool (*callback_key_down)(Viewer& viewer, unsigned char key, int modifiers);
|
||
bool (*callback_key_up)(Viewer& viewer, unsigned char key, int modifiers);
|
||
</code></pre>
|
||
|
||
<p>A keyboard callback can be used to visualize multiple meshes or different
|
||
stages of an algorithm, as demonstrated in <a href="103_Events/main.cpp">Example 103</a>, where
|
||
the keyboard callback changes the visualized mesh depending on the key pressed:</p>
|
||
|
||
<pre><code class="cpp">bool key_down(igl::viewer::Viewer& viewer, unsigned char key, int modifier)
|
||
{
|
||
if (key == '1')
|
||
{
|
||
viewer.data.clear();
|
||
viewer.data.set_mesh(V1, F1);
|
||
viewer.core.align_camera_center(V1,F1);
|
||
}
|
||
else if (key == '2')
|
||
{
|
||
viewer.data.clear();
|
||
viewer.data.set_mesh(V2, F2);
|
||
viewer.core.align_camera_center(V2,F2);
|
||
}
|
||
return false;
|
||
}
|
||
</code></pre>
|
||
|
||
<p>The callback is registered in the viewer as follows:</p>
|
||
|
||
<pre><code class="cpp">viewer.callback_key_down = &key_down;
|
||
</code></pre>
|
||
|
||
<p>Note that the mesh is cleared before using set_mesh. This has to be called
|
||
every time the number of vertices or faces of the plotted mesh changes. Every
|
||
callback returns a boolean value that tells the viewer if the event has been
|
||
handled by the plugin, or if the viewer should process it normally. This is
|
||
useful, for example, to disable the default mouse event handling if you want to
|
||
control the camera directly in your code.</p>
|
||
|
||
<p>The viewer can be extended using plugins, which are classes that implements all
|
||
the viewer’s callbacks. See the
|
||
<a href="../include/igl/viewer/ViewerPlugin.h">Viewer_plugin</a> for more details.</p>
|
||
|
||
<h2 id="scalarfieldvisualization"><a href="#scalarfieldvisualization">Scalar field visualization</a></h2>
|
||
|
||
<p>Colors and normals can be associated to faces or vertices using the
|
||
set_colors function:</p>
|
||
|
||
<pre><code class="cpp">viewer.data.set_colors(C);
|
||
</code></pre>
|
||
|
||
<p><code>C</code> is a #C by 3 matrix with one RGB color per row. <code>C</code> must have as many
|
||
rows as the number of faces <strong>or</strong> the number of vertices of the mesh.
|
||
Depending on the size of <code>C</code>, the viewer applies the color to the faces or to
|
||
the vertices.</p>
|
||
|
||
<p>Colors can be used to visualize a scalar function defined on a surface. The
|
||
scalar function is converted to colors using a color transfer function, which
|
||
maps a scalar value between 0 and 1 to a color. A simple example of a scalar
|
||
field defined on a surface is the z coordinate of each point, which can be
|
||
extract from our mesh representation by taking the last column of <code>V</code>
|
||
(<a href="104_Colors/main.cpp">Example 104</a>). The function <code>igl::jet</code> can be used to
|
||
convert it to colors:</p>
|
||
|
||
<pre><code class="cpp">Eigen::VectorXd Z = V.col(2);
|
||
igl::jet(Z,true,C);
|
||
</code></pre>
|
||
|
||
<p>The first row extracts the third column from <code>V</code> (the z coordinate of each
|
||
vertex) and the second calls a libigl functions that converts a scalar field to colors. The second parameter of jet normalizes the scalar field to lie between 0 and 1 before applying the transfer function.</p>
|
||
|
||
<figure>
|
||
<img src="images/104_Colors.png" alt="(Example 104) igl::jet converts a scalar field to a
|
||
color field." />
|
||
<figcaption>(<a href="104_Colors/main.cpp">Example 104</a>) igl::jet converts a scalar field to a
|
||
color field.</figcaption>
|
||
</figure>
|
||
|
||
<p><code>igl::jet</code> is an example of a standard function in libigl: it takes simple
|
||
types and can be easily reused for many different tasks. Not committing to
|
||
heavy data structures types favors simplicity, ease of use and reusability.</p>
|
||
|
||
<h2 id="overlays"><a href="#overlays">Overlays</a></h2>
|
||
|
||
<p>In addition to plotting the surface, the viewer supports the visualization of points, lines and text labels: these overlays can be very helpful while developing geometric processing algorithms to plot debug informations.</p>
|
||
|
||
<pre><code class="cpp">viewer.data.add_points(P,Eigen::RowVector3d(r,g,b));
|
||
</code></pre>
|
||
|
||
<p>Draws a point of color r,g,b for each row of P. The point is placed at the coordinates specified in each row of P, which is a #P by 3 matrix.</p>
|
||
|
||
<pre><code class="cpp">viewer.data.add_edges(P1,P2,Eigen::RowVector3d(r,g,b);
|
||
</code></pre>
|
||
|
||
<p>Draws a line of color r,g,b for each row of P1 and P2, which connects the 3D point in to the point in P2. Both P1 and P2 are of size #P by 3.</p>
|
||
|
||
<pre><code class="cpp">viewer.data.add_label(p,str);
|
||
</code></pre>
|
||
|
||
<p>Draws a label containing the string str at the position p, which is a vector of length 3.</p>
|
||
|
||
<p>These functions are demonstrate in <a href="105_Overlays/main.cpp">Example 105</a> where
|
||
the bounding box of a mesh is plotted using lines and points.
|
||
Using matrices to encode the mesh and its attributes allows to write short and
|
||
efficient code for many operations, avoiding to write for loops. For example,
|
||
the bounding box of a mesh can be found by taking the colwise maximum and minimum of <code>V</code>:</p>
|
||
|
||
<pre><code class="cpp">Eigen::Vector3d m = V.colwise().minCoeff();
|
||
Eigen::Vector3d M = V.colwise().maxCoeff();
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/105_Overlays.png" alt="(Example 105) The bounding box of a mesh is shown
|
||
using overlays." />
|
||
<figcaption>(<a href="105_Overlays/main.cpp">Example 105</a>) The bounding box of a mesh is shown
|
||
using overlays.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="viewermenu"><a href="#viewermenu">Viewer Menu</a></h2>
|
||
|
||
<p>As of version 1.2 the viewer uses a new menu and completely replaces
|
||
<a href="http://anttweakbar.sourceforge.net/doc/">AntTweakBar</a>. It is based on the
|
||
open-source projects <a href="https://github.com/memononen/nanovg">nanovg</a> and
|
||
<a href="https://github.com/wjakob/nanogui">nanogui</a>. To extend the default menu of the
|
||
viewer and to expose more user defined variables you have to define a callback
|
||
function:</p>
|
||
|
||
<pre><code class="cpp">igl::viewer::Viewer viewer;
|
||
|
||
bool boolVariable = true;
|
||
float floatVariable = 0.1f;
|
||
enum Orientation { Up=0,Down,Left,Right } dir = Up;
|
||
|
||
// Extend viewer menu
|
||
viewer.callback_init = [&](igl::viewer::Viewer& viewer)
|
||
{
|
||
// Add new group
|
||
viewer.ngui->addGroup("New Group");
|
||
|
||
// Expose a variable directly ...
|
||
viewer.ngui->addVariable("float",floatVariable);
|
||
|
||
// Expose an enumaration type
|
||
viewer.ngui->addVariable<Orientation>("Direction",dir)->setItems({"Up","Down","Left","Right"});
|
||
|
||
// Add a button
|
||
viewer.ngui->addButton("Print Hello",[](){ std::cout << "Hello\n"; });
|
||
|
||
// call to generate menu
|
||
viewer.screen->performLayout();
|
||
return false;
|
||
};
|
||
|
||
// start viewer
|
||
viewer.launch();
|
||
</code></pre>
|
||
|
||
<p>If you need a separate new menu window use:</p>
|
||
|
||
<pre><code class="cpp">viewer.ngui->addWindow(Eigen::Vector2i(220,10),"New Window");
|
||
</code></pre>
|
||
|
||
<p>If you do not want to expose variables directly but rather use the get/set functionality:</p>
|
||
|
||
<pre><code class="cpp">// ... or using a custom callback
|
||
viewer.ngui->addVariable<bool>("bool",[&](bool val) {
|
||
boolVariable = val; // setter
|
||
},[&]() {
|
||
return boolVariable; // getter
|
||
});
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/106_ViewerMenu.png" alt="(Example 106) The UI of the viewer can be easily customized." />
|
||
<figcaption>(<a href="106_ViewerMenu/main.cpp">Example 106</a>) The UI of the viewer can be easily customized.</figcaption>
|
||
</figure>
|
||
|
||
<h1 id="chapter2:discretegeometricquantitiesandoperators">Chapter 2: Discrete Geometric Quantities and Operators</h1>
|
||
|
||
<p>This chapter illustrates a few discrete quantities that libigl can compute on a
|
||
mesh and the libigl functions that construct popular discrete differential
|
||
geometry operators. It also provides an introduction to basic drawing and coloring routines of our viewer.</p>
|
||
|
||
<h2 id="normals">Normals</h2>
|
||
|
||
<p>Surface normals are a basic quantity necessary for rendering a surface. There
|
||
are a variety of ways to compute and store normals on a triangle mesh. <a href="201_Normals/main.cpp">Example 201</a> demonstrates how to compute and visualize normals with libigl.</p>
|
||
|
||
<h3 id="per-face">Per-face</h3>
|
||
|
||
<p>Normals are well defined on each triangle of a mesh as the vector orthogonal to
|
||
triangle’s plane. These piecewise-constant normals produce piecewise-flat
|
||
renderings: the surface appears non-smooth and reveals its underlying
|
||
discretization.</p>
|
||
|
||
<h3 id="per-vertex">Per-vertex</h3>
|
||
|
||
<p>Normals can be computed and stored on vertices, and interpolated in the interior of the triangles to produce smooth renderings (<a href="http://en.wikipedia.org/wiki/Phong_shading">Phong shading</a>).
|
||
Most techniques for computing per-vertex normals take an average of incident face normals. The main difference between these techniques is their weighting scheme: Uniform
|
||
weighting is heavily biased by the discretization choice, whereas area-based
|
||
or angle-based weighting is more forgiving.</p>
|
||
|
||
<p>The typical half-edge style computation of area-based weights has this structure:</p>
|
||
|
||
<pre><code class="cpp">N.setZero(V.rows(),3);
|
||
for(int i : vertices)
|
||
{
|
||
for(face : incident_faces(i))
|
||
{
|
||
N.row(i) += face.area * face.normal;
|
||
}
|
||
}
|
||
N.rowwise().normalize();
|
||
</code></pre>
|
||
|
||
<p>At first glance, it might seem inefficient to loop over incident faces—and thus constructing the per-vertex normals— without using an half-edge data structure. However, per-vertex normals may be <em>throwing</em> each face normal to
|
||
running sums on its corner vertices:</p>
|
||
|
||
<pre><code class="cpp">N.setZero(V.rows(),3);
|
||
for(int f = 0; f < F.rows();f++)
|
||
{
|
||
for(int c = 0; c < 3;c++)
|
||
{
|
||
N.row(F(f,c)) += area(f) * face_normal.row(f);
|
||
}
|
||
}
|
||
N.rowwise().normalize();
|
||
</code></pre>
|
||
|
||
<h3 id="per-corner">Per-corner</h3>
|
||
|
||
<p>Storing normals per-corner is an efficient and convenient way of supporting both
|
||
smooth and sharp (e.g. creases and corners) rendering. This format is common to
|
||
OpenGL and the .obj mesh file format. Often such normals are tuned by the mesh
|
||
designer, but creases and corners can also be computed automatically. Libigl
|
||
implements a simple scheme which computes corner normals as averages of
|
||
normals of faces incident on the corresponding vertex which do not deviate by more than a specified dihedral angle (e.g. 20°).</p>
|
||
|
||
<figure>
|
||
<img src="images/fandisk-normals.jpg" alt="The Normals example computes per-face (left), per-vertex (middle) and
|
||
per-corner (right) normals" />
|
||
<figcaption>The <code>Normals</code> example computes per-face (left), per-vertex (middle) and
|
||
per-corner (right) normals</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="gaussiancurvature">Gaussian curvature</h2>
|
||
|
||
<p>Gaussian curvature on a continuous surface is defined as the product of the
|
||
principal curvatures:</p>
|
||
|
||
<p><span class="math">\(k_G = k_1 k_2.\)</span></p>
|
||
|
||
<p>As an <em>intrinsic</em> measure, it depends on the metric and
|
||
not the surface’s embedding.</p>
|
||
|
||
<p>Intuitively, Gaussian curvature tells how locally spherical or <em>elliptic</em> the
|
||
surface is ( <span class="math">\(k_G>0\)</span> ), how locally saddle-shaped or <em>hyperbolic</em> the surface
|
||
is ( <span class="math">\(k_G<0\)</span> ), or how locally cylindrical or <em>parabolic</em> ( <span class="math">\(k_G=0\)</span> ) the
|
||
surface is.</p>
|
||
|
||
<p>In the discrete setting, one definition for a “discrete Gaussian curvature”
|
||
on a triangle mesh is via a vertex’s <em>angular deficit</em>:</p>
|
||
|
||
<p><span class="math">\(k_G(v_i) = 2π - \sum\limits_{j\in N(i)}θ_{ij},\)</span></p>
|
||
|
||
<p>where <span class="math">\(N(i)\)</span> are the triangles incident on vertex <span class="math">\(i\)</span> and <span class="math">\(θ_{ij}\)</span> is the angle
|
||
at vertex <span class="math">\(i\)</span> in triangle <span class="math">\(j\)</span> <a class="citation" href="#fn:1" title="Jump to citation">[1]<span class="citekey" style="display:none">meyer_2003</span></a>.</p>
|
||
|
||
<p>Just like the continuous analog, our discrete Gaussian curvature reveals
|
||
elliptic, hyperbolic and parabolic vertices on the domain, as demonstrated in <a href="202GaussianCurvature/main.cpp">Example 202</a>.</p>
|
||
|
||
<figure>
|
||
<img src="images/bumpy-gaussian-curvature.jpg" alt="The GaussianCurvature example computes discrete Gaussian curvature and
|
||
visualizes it in pseudocolor." />
|
||
<figcaption>The <code>GaussianCurvature</code> example computes discrete Gaussian curvature and
|
||
visualizes it in pseudocolor.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="curvaturedirections">Curvature directions</h2>
|
||
|
||
<p>The two principal curvatures <span class="math">\((k_1,k_2)\)</span> at a point on a surface measure how
|
||
much the surface bends in different directions. The directions of maximum and
|
||
minimum (signed) bending are called principal directions and are always
|
||
orthogonal.</p>
|
||
|
||
<p>Mean curvature is defined as the average of principal curvatures:</p>
|
||
|
||
<p><span class="math">\(H = \frac{1}{2}(k_1 + k_2).\)</span></p>
|
||
|
||
<p>One way to extract mean curvature is by examining the Laplace-Beltrami operator
|
||
applied to the surface positions. The result is a so-called mean-curvature
|
||
normal:</p>
|
||
|
||
<p><span class="math">\(-\Delta \mathbf{x} = H \mathbf{n}.\)</span></p>
|
||
|
||
<p>It is easy to compute this on a discrete triangle mesh in libigl using the
|
||
cotangent Laplace-Beltrami operator <a class="citation" href="#fn:1" title="Jump to citation">[1]<span class="citekey" style="display:none">meyer_2003</span></a>.</p>
|
||
|
||
<pre><code class="cpp">#include <igl/cotmatrix.h>
|
||
#include <igl/massmatrix.h>
|
||
#include <igl/invert_diag.h>
|
||
...
|
||
MatrixXd HN;
|
||
SparseMatrix<double> L,M,Minv;
|
||
igl::cotmatrix(V,F,L);
|
||
igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_VORONOI,M);
|
||
igl::invert_diag(M,Minv);
|
||
HN = -Minv*(L*V);
|
||
H = HN.rowwise().norm(); //up to sign
|
||
</code></pre>
|
||
|
||
<p>Combined with the angle defect definition of discrete Gaussian curvature, one
|
||
can define principal curvatures and use least squares fitting to find
|
||
directions <a class="citation" href="#fn:1" title="Jump to citation">[1]<span class="citekey" style="display:none">meyer_2003</span></a>.</p>
|
||
|
||
<p>Alternatively, a robust method for determining principal curvatures is via
|
||
quadric fitting <a class="citation" href="#fn:2" title="Jump to citation">[2]<span class="citekey" style="display:none">panozzo_2010</span></a>. In the neighborhood around every vertex, a
|
||
best-fit quadric is found and principal curvature values and directions are
|
||
analytically computed on this quadric (<a href="203_curvatureDirections/main.cpp">Example
|
||
203</a>).</p>
|
||
|
||
<figure>
|
||
<img src="images/fertility-principal-curvature.jpg" alt="The CurvatureDirections example computes principal curvatures via quadric
|
||
fitting and visualizes mean curvature in pseudocolor and principal directions
|
||
with a cross field." />
|
||
<figcaption>The <code>CurvatureDirections</code> example computes principal curvatures via quadric
|
||
fitting and visualizes mean curvature in pseudocolor and principal directions
|
||
with a cross field.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="gradient">Gradient</h2>
|
||
|
||
<p>Scalar functions on a surface can be discretized as a piecewise linear function
|
||
with values defined at each mesh vertex:</p>
|
||
|
||
<p><span class="math">\(f(\mathbf{x}) \approx \sum\limits_{i=1}^n \phi_i(\mathbf{x})\, f_i,\)</span></p>
|
||
|
||
<p>where <span class="math">\(\phi_i\)</span> is a piecewise linear hat function defined by the mesh so that
|
||
for each triangle <span class="math">\(\phi_i\)</span> is <em>the</em> linear function which is one only at
|
||
vertex <span class="math">\(i\)</span> and zero at the other corners.</p>
|
||
|
||
<figure>
|
||
<img src="images/hat-function.jpg" alt="Hat function $\phi_i$ is one at vertex $i$, zero at all other vertices, and
|
||
linear on incident triangles." />
|
||
<figcaption>Hat function <span class="math">\(\phi_i\)</span> is one at vertex <span class="math">\(i\)</span>, zero at all other vertices, and
|
||
linear on incident triangles.</figcaption>
|
||
</figure>
|
||
|
||
<p>Thus gradients of such piecewise linear functions are simply sums of gradients
|
||
of the hat functions:</p>
|
||
|
||
<p><span class="math">\(\nabla f(\mathbf{x}) \approx
|
||
\nabla \sum\limits_{i=1}^n \phi_i(\mathbf{x})\, f_i =
|
||
\sum\limits_{i=1}^n \nabla \phi_i(\mathbf{x})\, f_i.\)</span></p>
|
||
|
||
<p>This reveals that the gradient is a linear function of the vector of <span class="math">\(f_i\)</span>
|
||
values. Because the <span class="math">\(\phi_i\)</span> are linear in each triangle, their gradients are
|
||
<em>constant</em> in each triangle. Thus our discrete gradient operator can be written
|
||
as a matrix multiplication taking vertex values to triangle values:</p>
|
||
|
||
<p><span class="math">\(\nabla f \approx \mathbf{G}\,\mathbf{f},\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{f}\)</span> is <span class="math">\(n\times 1\)</span> and <span class="math">\(\mathbf{G}\)</span> is an <span class="math">\(md\times n\)</span> sparse
|
||
matrix. This matrix <span class="math">\(\mathbf{G}\)</span> can be derived geometrically, e.g.
|
||
<a class="citation" href="#fn:3" title="Jump to citation">[<span class="locator">ch. 2</span>, 3]<span class="citekey" style="display:none">jacobson_thesis_2013</span></a>.
|
||
Libigl’s <code>grad</code> function computes <span class="math">\(\mathbf{G}\)</span> for
|
||
triangle and tetrahedral meshes (<a href="204_Gradient/main.cpp">Example 204</a>):</p>
|
||
|
||
<figure>
|
||
<img src="images/cheburashka-gradient.jpg" alt="The Gradient example computes gradients of an input function on a mesh and
|
||
visualizes the vector field." />
|
||
<figcaption>The <code>Gradient</code> example computes gradients of an input function on a mesh and
|
||
visualizes the vector field.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="laplacian">Laplacian</h2>
|
||
|
||
<p>The discrete Laplacian is an essential geometry processing tool. Many
|
||
interpretations and flavors of the Laplace and Laplace-Beltrami operator exist.</p>
|
||
|
||
<p>In open Euclidean space, the <em>Laplace</em> operator is the usual divergence of
|
||
gradient (or equivalently the Laplacian of a function is the trace of its
|
||
Hessian):</p>
|
||
|
||
<p><span class="math">\(\Delta f =
|
||
\frac{\partial^2 f}{\partial x^2} +
|
||
\frac{\partial^2 f}{\partial y^2} +
|
||
\frac{\partial^2 f}{\partial z^2}.\)</span></p>
|
||
|
||
<p>The <em>Laplace-Beltrami</em> operator generalizes this to surfaces.</p>
|
||
|
||
<p>When considering piecewise-linear functions on a triangle mesh, a discrete
|
||
Laplacian may be derived in a variety of ways. The most popular in geometry
|
||
processing is the so-called “cotangent Laplacian” <span class="math">\(\mathbf{L}\)</span>, arising
|
||
simultaneously from FEM, DEC and applying divergence theorem to vertex
|
||
one-rings. As a linear operator taking vertex values to vertex values, the
|
||
Laplacian <span class="math">\(\mathbf{L}\)</span> is a <span class="math">\(n\times n\)</span> matrix with elements:</p>
|
||
|
||
<p><span class="math">\(L_{ij} = \begin{cases}j \in N(i) &\cot \alpha_{ij} + \cot \beta_{ij},\\
|
||
j \notin N(i) & 0,\\
|
||
i = j & -\sum\limits_{k\neq i} L_{ik},
|
||
\end{cases}\)</span></p>
|
||
|
||
<p>where <span class="math">\(N(i)\)</span> are the vertices adjacent to (neighboring) vertex <span class="math">\(i\)</span>, and
|
||
<span class="math">\(\alpha_{ij},\beta_{ij}\)</span> are the angles opposite to edge <span class="math">\({ij}\)</span>.
|
||
This formula leads to a typical half-edge style implementation for
|
||
constructing <span class="math">\(\mathbf{L}\)</span>:</p>
|
||
|
||
<pre><code class="cpp">for(int i : vertices)
|
||
{
|
||
for(int j : one_ring(i))
|
||
{
|
||
for(int k : triangle_on_edge(i,j))
|
||
{
|
||
L(i,j) = cot(angle(i,j,k));
|
||
L(i,i) -= cot(angle(i,j,k));
|
||
}
|
||
}
|
||
}
|
||
</code></pre>
|
||
|
||
<p>Similarly as before, it may seem to loop over one-rings without having an half-edge data structure. However, this is not the case, since the Laplacian may be built by summing together contributions for each triangle, much in spirit with its FEM discretization
|
||
of the Dirichlet energy (sum of squared gradients):</p>
|
||
|
||
<pre><code class="cpp">for(triangle t : triangles)
|
||
{
|
||
for(edge i,j : t)
|
||
{
|
||
L(i,j) += cot(angle(i,j,k));
|
||
L(j,i) += cot(angle(i,j,k));
|
||
L(i,i) -= cot(angle(i,j,k));
|
||
L(j,j) -= cot(angle(i,j,k));
|
||
}
|
||
}
|
||
</code></pre>
|
||
|
||
<p>Libigl implements discrete “cotangent” Laplacians for triangles meshes and
|
||
tetrahedral meshes, building both with fast geometric rules rather than “by the
|
||
book” FEM construction which involves many (small) matrix inversions, cf.
|
||
<a class="citation" href="#fn:4" title="Jump to citation">[4]<span class="citekey" style="display:none">sharf_2007</span></a>.</p>
|
||
|
||
<p>The operator applied to mesh vertex positions amounts to smoothing by <em>flowing</em>
|
||
the surface along the mean curvature normal direction (<a href="205_Laplacian/main.cpp">Example 205</a>). Note that this is equivalent to minimizing surface area.</p>
|
||
|
||
<figure>
|
||
<img src="images/cow-curvature-flow.jpg" alt="The Laplacian example computes conformalized mean curvature flow using the
|
||
cotangent Laplacian ." />
|
||
<figcaption>The <code>Laplacian</code> example computes conformalized mean curvature flow using the
|
||
cotangent Laplacian <a class="citation" href="#fn:5" title="Jump to citation">[5]<span class="citekey" style="display:none">kazhdan_2012</span></a>.</figcaption>
|
||
</figure>
|
||
|
||
<h3 id="massmatrix">Mass matrix</h3>
|
||
|
||
<p>The mass matrix <span class="math">\(\mathbf{M}\)</span> is another <span class="math">\(n \times n\)</span> matrix which takes vertex
|
||
values to vertex values. From an FEM point of view, it is a discretization of
|
||
the inner-product: it accounts for the area around each vertex. Consequently,
|
||
<span class="math">\(\mathbf{M}\)</span> is often a diagonal matrix, such that <span class="math">\(M_{ii}\)</span> is the barycentric
|
||
or voronoi area around vertex <span class="math">\(i\)</span> in the mesh <a class="citation" href="#fn:1" title="Jump to citation">[1]<span class="citekey" style="display:none">meyer_2003</span></a>. The inverse of
|
||
this matrix is also very useful as it transforms integrated quantities into
|
||
point-wise quantities, e.g.:</p>
|
||
|
||
<p><span class="math">\(\Delta f \approx \mathbf{M}^{-1} \mathbf{L} \mathbf{f}.\)</span></p>
|
||
|
||
<p>In general, when encountering squared quantities integrated over the surface,
|
||
the mass matrix will be used as the discretization of the inner product when
|
||
sampling function values at vertices:</p>
|
||
|
||
<p><span class="math">\(\int_S x\, y\ dA \approx \mathbf{x}^T\mathbf{M}\,\mathbf{y}.\)</span></p>
|
||
|
||
<p>An alternative mass matrix <span class="math">\(\mathbf{T}\)</span> is a <span class="math">\(md \times md\)</span> matrix which takes
|
||
triangle vector values to triangle vector values. This matrix represents an
|
||
inner-product accounting for the area associated with each triangle (i.e. the
|
||
triangles true area).</p>
|
||
|
||
<h3 id="alternativeconstructionoflaplacian">Alternative construction of Laplacian</h3>
|
||
|
||
<p>An alternative construction of the discrete cotangent Laplacian is by
|
||
“squaring” the discrete gradient operator. This may be derived by applying
|
||
Green’s identity (ignoring boundary conditions for the moment):</p>
|
||
|
||
<p><span class="math">\(\int_S \|\nabla f\|^2 dA = \int_S f \Delta f dA\)</span></p>
|
||
|
||
<p>Or in matrix form which is immediately translatable to code:</p>
|
||
|
||
<p><span class="math">\(\mathbf{f}^T \mathbf{G}^T \mathbf{T} \mathbf{G} \mathbf{f} =
|
||
\mathbf{f}^T \mathbf{M} \mathbf{M}^{-1} \mathbf{L} \mathbf{f} =
|
||
\mathbf{f}^T \mathbf{L} \mathbf{f}.\)</span></p>
|
||
|
||
<p>So we have that <span class="math">\(\mathbf{L} = \mathbf{G}^T \mathbf{T} \mathbf{G}\)</span>. This also
|
||
hints that we may consider <span class="math">\(\mathbf{G}^T\)</span> as a discrete <em>divergence</em> operator,
|
||
since the Laplacian is the divergence of the gradient. Naturally, <span class="math">\(\mathbf{G}^T\)</span> is
|
||
a <span class="math">\(n \times md\)</span> sparse matrix which takes vector values stored at triangle faces
|
||
to scalar divergence values at vertices.</p>
|
||
|
||
<h1 id="chapter3:matricesandlinearalgebra">Chapter 3: Matrices and linear algebra</h1>
|
||
|
||
<p>Libigl relies heavily on the Eigen library for dense and sparse linear algebra
|
||
routines. Besides geometry processing routines, libigl has linear algebra
|
||
routines which bootstrap Eigen and make it feel even more similar to a high-level
|
||
algebra library such as Matlab.</p>
|
||
|
||
<h2 id="slice">Slice</h2>
|
||
|
||
<p>A very familiar and powerful routine in Matlab is array slicing. This allows
|
||
reading from or writing to a possibly non-contiguous sub-matrix. Let’s consider
|
||
the Matlab code:</p>
|
||
|
||
<pre><code class="matlab">B = A(R,C);
|
||
</code></pre>
|
||
|
||
<p>If <code>A</code> is a <span class="math">\(m \times n\)</span> matrix and <code>R</code> is a <span class="math">\(j\)</span>-long list of row-indices
|
||
(between 1 and <span class="math">\(m\)</span>) and <code>C</code> is a <span class="math">\(k\)</span>-long list of column-indices, then as a
|
||
result <code>B</code> will be a <span class="math">\(j \times k\)</span> matrix drawing elements from <code>A</code> according to
|
||
<code>R</code> and <code>C</code>. In libigl, the same functionality is provided by the <code>slice</code>
|
||
function (<a href="301_Slice/main.cpp">Example 301</a>):</p>
|
||
|
||
<pre><code class="cpp">VectorXi R,C;
|
||
MatrixXd A,B;
|
||
...
|
||
igl::slice(A,R,C,B);
|
||
</code></pre>
|
||
|
||
<p>Note that <code>A</code> and <code>B</code> could also be sparse matrices.</p>
|
||
|
||
<p>Similarly, consider the Matlab code:</p>
|
||
|
||
<pre><code class="matlab">A(R,C) = B;
|
||
</code></pre>
|
||
|
||
<p>Now, the selection is on the left-hand side so the <span class="math">\(j \times k\)</span> matrix <code>B</code> is
|
||
being <em>written into</em> the submatrix of <code>A</code> determined by <code>R</code> and <code>C</code>. This
|
||
functionality is provided in libigl using <code>slice_into</code>:</p>
|
||
|
||
<pre><code class="cpp">igl::slice_into(B,R,C,A);
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/decimated-knight-slice-color.jpg" alt="The example Slice shows how to use igl::slice to change the colors for
|
||
triangles on a mesh." />
|
||
<figcaption>The example <code>Slice</code> shows how to use <code>igl::slice</code> to change the colors for
|
||
triangles on a mesh.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="sort">Sort</h2>
|
||
|
||
<p>Matlab and other higher-level languages make it very easy to extract indices of
|
||
sorting and comparison routines. For example in Matlab, one can write:</p>
|
||
|
||
<pre><code class="matlab">[Y,I] = sort(X,1,'ascend');
|
||
</code></pre>
|
||
|
||
<p>so if <code>X</code> is a <span class="math">\(m \times n\)</span> matrix then <code>Y</code> will also be an <span class="math">\(m \times n\)</span> matrix
|
||
with entries sorted along dimension <code>1</code> in <code>'ascend'</code>ing order. The second
|
||
output <code>I</code> is a <span class="math">\(m \times n\)</span> matrix of indices such that <code>Y(i,j) =
|
||
X(I(i,j),j);</code>. That is, <code>I</code> reveals how <code>X</code> is sorted into <code>Y</code>.</p>
|
||
|
||
<p>This same functionality is supported in libigl:</p>
|
||
|
||
<pre><code class="cpp">igl::sort(X,1,true,Y,I);
|
||
</code></pre>
|
||
|
||
<p>Similarly, sorting entire rows can be accomplished in Matlab using:</p>
|
||
|
||
<pre><code class="matlab">[Y,I] = sortrows(X,'ascend');
|
||
</code></pre>
|
||
|
||
<p>where now <code>I</code> is a <span class="math">\(m\)</span> vector of indices such that <code>Y = X(I,:)</code>.</p>
|
||
|
||
<p>In libigl, this is supported with</p>
|
||
|
||
<pre><code class="cpp">igl::sortrows(X,true,Y,I);
|
||
</code></pre>
|
||
|
||
<p>where again <code>I</code> reveals the index of sort so that it can be reproduced with
|
||
<code>igl::slice(X,I,1,Y)</code>.</p>
|
||
|
||
<p>Analogous functions are available in libigl for: <code>max</code>, <code>min</code>, and <code>unique</code>.</p>
|
||
|
||
<figure>
|
||
<img src="images/decimated-knight-sort-color.jpg" alt="The example Sort shows how to use igl::sortrows to
|
||
pseudocolor triangles according to their barycenters' sorted
|
||
order (Example 302)." />
|
||
<figcaption>The example <code>Sort</code> shows how to use <code>igl::sortrows</code> to
|
||
pseudocolor triangles according to their barycenters’ sorted
|
||
order (<a href="302_Sort/main.cpp">Example 302</a>).</figcaption>
|
||
</figure>
|
||
|
||
<h3 id="othermatlab-stylefunctions">Other Matlab-style functions</h3>
|
||
|
||
<p>Libigl implements a variety of other routines with the same api and
|
||
functionality as common Matlab functions.</p>
|
||
|
||
<table>
|
||
<colgroup>
|
||
<col style="text-align:left;"/>
|
||
<col style="text-align:left;"/>
|
||
</colgroup>
|
||
|
||
<thead>
|
||
<tr>
|
||
<th style="text-align:left;">Name</th>
|
||
<th style="text-align:left;">Description</th>
|
||
</tr>
|
||
</thead>
|
||
|
||
<tbody>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::all</code></td>
|
||
<td style="text-align:left;">Whether all elements are non-zero (true)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::any</code></td>
|
||
<td style="text-align:left;">Whether any elements are non-zero (true)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::cat</code></td>
|
||
<td style="text-align:left;">Concatenate two matrices (especially useful for dealing with Eigen sparse matrices)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::ceil</code></td>
|
||
<td style="text-align:left;">Round entries up to nearest integer</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::cumsum</code></td>
|
||
<td style="text-align:left;">Cumulative sum of matrix elements</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::colon</code></td>
|
||
<td style="text-align:left;">Act like Matlab’s <code>:</code>, similar to Eigen’s <code>LinSpaced</code></td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::components</code></td>
|
||
<td style="text-align:left;">Connected components of graph (cf. Matlab’s <code>graphconncomp</code>)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::count</code></td>
|
||
<td style="text-align:left;">Count non-zeros in rows or columns</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::cross</code></td>
|
||
<td style="text-align:left;">Cross product per-row</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::cumsum</code></td>
|
||
<td style="text-align:left;">Cummulative summation</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::dot</code></td>
|
||
<td style="text-align:left;">dot product per-row</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::eigs</code></td>
|
||
<td style="text-align:left;">Solve sparse eigen value problem</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::find</code></td>
|
||
<td style="text-align:left;">Find subscripts of non-zero entries</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::floor</code></td>
|
||
<td style="text-align:left;">Round entries down to nearest integer</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::histc</code></td>
|
||
<td style="text-align:left;">Counting occurrences for building a histogram</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::hsv_to_rgb</code></td>
|
||
<td style="text-align:left;">Convert HSV colors to RGB (cf. Matlab’s <code>hsv2rgb</code>)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::intersect</code></td>
|
||
<td style="text-align:left;">Set intersection of matrix elements.</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::isdiag</code></td>
|
||
<td style="text-align:left;">Determine whether matrix is diagonal</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::ismember</code></td>
|
||
<td style="text-align:left;">Determine whether elements in A occur in B</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::jet</code></td>
|
||
<td style="text-align:left;">Quantized colors along the rainbow.</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::max</code></td>
|
||
<td style="text-align:left;">Compute maximum entry per row or column</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::median</code></td>
|
||
<td style="text-align:left;">Compute the median per column</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::min</code></td>
|
||
<td style="text-align:left;">Compute minimum entry per row or column</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::mod</code></td>
|
||
<td style="text-align:left;">Compute per element modulo</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::mode</code></td>
|
||
<td style="text-align:left;">Compute the mode per column</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::null</code></td>
|
||
<td style="text-align:left;">Compute the null space basis of a matrix</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::nchoosek</code></td>
|
||
<td style="text-align:left;">Compute all k-size combinations of n-long vector</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::orth</code></td>
|
||
<td style="text-align:left;">Orthogonalization of a basis</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::parula</code></td>
|
||
<td style="text-align:left;">Generate a quantized colormap from blue to yellow</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::pinv</code></td>
|
||
<td style="text-align:left;">Compute Moore-Penrose pseudoinverse</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::randperm</code></td>
|
||
<td style="text-align:left;">Generate a random permutation of [0,…,n–1]</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::rgb_to_hsv</code></td>
|
||
<td style="text-align:left;">Convert RGB colors to HSV (cf. Matlab’s <code>rgb2hsv</code>)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::repmat</code></td>
|
||
<td style="text-align:left;">Repeat a matrix along columns and rows</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::round</code></td>
|
||
<td style="text-align:left;">Per-element round to whole number</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::setdiff</code></td>
|
||
<td style="text-align:left;">Set difference of matrix elements</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::setunion</code></td>
|
||
<td style="text-align:left;">Set union of matrix elements</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::setxor</code></td>
|
||
<td style="text-align:left;">Set exclusive “or” of matrix elements</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::slice</code></td>
|
||
<td style="text-align:left;">Slice parts of matrix using index lists: (cf. Matlab’s <code>B = A(I,J)</code>)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::slice_mask</code></td>
|
||
<td style="text-align:left;">Slice parts of matrix using boolean masks: (cf. Matlab’s <code>B = A(M,N)</code>)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::slice_into</code></td>
|
||
<td style="text-align:left;">Slice left-hand side of matrix assignment using index lists (cf. Matlab’s <code>B(I,J) = A</code>)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::sort</code></td>
|
||
<td style="text-align:left;">Sort elements or rows of matrix</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::speye</code></td>
|
||
<td style="text-align:left;">Identity as sparse matrix</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::sum</code></td>
|
||
<td style="text-align:left;">Sum along columns or rows (of sparse matrix)</td>
|
||
</tr>
|
||
<tr>
|
||
<td style="text-align:left;"><code>igl::unique</code></td>
|
||
<td style="text-align:left;">Extract unique elements or rows of matrix</td>
|
||
</tr>
|
||
</tbody>
|
||
</table>
|
||
|
||
<h2 id="laplaceequation">Laplace equation</h2>
|
||
|
||
<p>A common linear system in geometry processing is the Laplace equation:</p>
|
||
|
||
<p><span class="math">\(∆z = 0\)</span></p>
|
||
|
||
<p>subject to some boundary conditions, for example Dirichlet boundary conditions
|
||
(fixed value):</p>
|
||
|
||
<p><span class="math">\(\left.z\right|_{\partial{S}} = z_{bc}\)</span></p>
|
||
|
||
<p>In the discrete setting, the linear system can be written as:</p>
|
||
|
||
<p><span class="math">\(\mathbf{L} \mathbf{z} = \mathbf{0}\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{L}\)</span> is the <span class="math">\(n \times n\)</span> discrete Laplacian and <span class="math">\(\mathbf{z}\)</span> is a
|
||
vector of per-vertex values. Most of <span class="math">\(\mathbf{z}\)</span> correspond to interior
|
||
vertices and are unknown, but some of <span class="math">\(\mathbf{z}\)</span> represent values at boundary
|
||
vertices. Their values are known so we may move their corresponding terms to
|
||
the right-hand side.</p>
|
||
|
||
<p>Conceptually, this is very easy if we have sorted <span class="math">\(\mathbf{z}\)</span> so that interior
|
||
vertices come first and then boundary vertices:</p>
|
||
|
||
<p><span class="math">\[\left(\begin{array}{cc}
|
||
\mathbf{L}_{in,in} & \mathbf{L}_{in,b}\\
|
||
\mathbf{L}_{b,in} & \mathbf{L}_{b,b}\end{array}\right)
|
||
\left(\begin{array}{c}
|
||
\mathbf{z}_{in}\\
|
||
\mathbf{z}_{b}\end{array}\right) =
|
||
\left(\begin{array}{c}
|
||
\mathbf{0}_{in}\\
|
||
\mathbf{z}_{bc}\end{array}\right)\]</span></p>
|
||
|
||
<p>The bottom block of equations is no longer meaningful so we’ll only consider
|
||
the top block:</p>
|
||
|
||
<p><span class="math">\[\left(\begin{array}{cc}
|
||
\mathbf{L}_{in,in} & \mathbf{L}_{in,b}\end{array}\right)
|
||
\left(\begin{array}{c}
|
||
\mathbf{z}_{in}\\
|
||
\mathbf{z}_{b}\end{array}\right) =
|
||
\mathbf{0}_{in}\]</span></p>
|
||
|
||
<p>We can move the known values to the right-hand side:</p>
|
||
|
||
<p><span class="math">\[\mathbf{L}_{in,in}
|
||
\mathbf{z}_{in} = -
|
||
\mathbf{L}_{in,b}
|
||
\mathbf{z}_{b}\]</span></p>
|
||
|
||
<p>Finally we can solve this equation for the unknown values at interior vertices
|
||
<span class="math">\(\mathbf{z}_{in}\)</span>.</p>
|
||
|
||
<p>However, our vertices will often not be sorted in this way. One option would be to sort <code>V</code>,
|
||
then proceed as above and then <em>unsort</em> the solution <code>Z</code> to match <code>V</code>. However,
|
||
this solution is not very general.</p>
|
||
|
||
<p>With array slicing no explicit sort is needed. Instead we can <em>slice-out</em>
|
||
submatrix blocks (<span class="math">\(\mathbf{L}_{in,in}\)</span>, <span class="math">\(\mathbf{L}_{in,b}\)</span>, etc.) and follow
|
||
the linear algebra above directly. Then we can slice the solution <em>into</em> the
|
||
rows of <code>Z</code> corresponding to the interior vertices (<a href="303_LaplaceEquation/main.cpp">Example 303</a>).</p>
|
||
|
||
<figure>
|
||
<img src="images/camelhead-laplace-equation.jpg" alt="The LaplaceEquation example solves a Laplace equation with Dirichlet
|
||
boundary conditions." />
|
||
<figcaption>The <code>LaplaceEquation</code> example solves a Laplace equation with Dirichlet
|
||
boundary conditions.</figcaption>
|
||
</figure>
|
||
|
||
<h3 id="quadraticenergyminimization">Quadratic energy minimization</h3>
|
||
|
||
<p>The same Laplace equation may be equivalently derived by minimizing Dirichlet
|
||
energy subject to the same boundary conditions:</p>
|
||
|
||
<p><span class="math">\(\mathop{\text{minimize }}_z \frac{1}{2}\int\limits_S \|\nabla z\|^2 dA\)</span></p>
|
||
|
||
<p>On our discrete mesh, recall that this becomes</p>
|
||
|
||
<p><span class="math">\(\mathop{\text{minimize }}_\mathbf{z} \frac{1}{2}\mathbf{z}^T \mathbf{G}^T \mathbf{D}
|
||
\mathbf{G} \mathbf{z} \rightarrow \mathop{\text{minimize }}_\mathbf{z} \mathbf{z}^T \mathbf{L} \mathbf{z}\)</span></p>
|
||
|
||
<p>The general problem of minimizing some energy over a mesh subject to fixed
|
||
value boundary conditions is so wide spread that libigl has a dedicated api for
|
||
solving such systems.</p>
|
||
|
||
<p>Let us consider a general quadratic minimization problem subject to different
|
||
common constraints:</p>
|
||
|
||
<p><span class="math">\[\mathop{\text{minimize }}_\mathbf{z} \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} +
|
||
\mathbf{z}^T \mathbf{B} + \text{constant},\]</span></p>
|
||
|
||
<p>subject to</p>
|
||
|
||
<p><span class="math">\[\mathbf{z}_b = \mathbf{z}_{bc} \text{ and } \mathbf{A}_{eq} \mathbf{z} =
|
||
\mathbf{B}_{eq},\]</span></p>
|
||
|
||
<p>where</p>
|
||
|
||
<ul>
|
||
<li><span class="math">\(\mathbf{Q}\)</span> is a (usually sparse) <span class="math">\(n \times n\)</span> positive semi-definite
|
||
matrix of quadratic coefficients (Hessian),</li>
|
||
<li><span class="math">\(\mathbf{B}\)</span> is a <span class="math">\(n \times 1\)</span> vector of linear coefficients,</li>
|
||
<li><span class="math">\(\mathbf{z}_b\)</span> is a <span class="math">\(|b| \times 1\)</span> portion of
|
||
<span class="math">\(\mathbf{z}\)</span> corresponding to boundary or <em>fixed</em> vertices,</li>
|
||
<li><span class="math">\(\mathbf{z}_{bc}\)</span> is a <span class="math">\(|b| \times 1\)</span> vector of known values corresponding to
|
||
<span class="math">\(\mathbf{z}_b\)</span>,</li>
|
||
<li><span class="math">\(\mathbf{A}_{eq}\)</span> is a (usually sparse) <span class="math">\(m \times n\)</span> matrix of linear
|
||
equality constraint coefficients (one row per constraint), and</li>
|
||
<li><span class="math">\(\mathbf{B}_{eq}\)</span> is a <span class="math">\(m \times 1\)</span> vector of linear equality constraint
|
||
right-hand side values.</li>
|
||
</ul>
|
||
|
||
<p>This specification is overly general as we could write <span class="math">\(\mathbf{z}_b =
|
||
\mathbf{z}_{bc}\)</span> as rows of <span class="math">\(\mathbf{A}_{eq} \mathbf{z} =
|
||
\mathbf{B}_{eq}\)</span>, but these fixed value constraints appear so often that they
|
||
merit a dedicated place in the API.</p>
|
||
|
||
<p>In libigl, solving such quadratic optimization problems is split into two
|
||
routines: precomputation and solve. Precomputation only depends on the
|
||
quadratic coefficients, known value indices and linear constraint coefficients:</p>
|
||
|
||
<pre><code class="cpp">igl::min_quad_with_fixed_data mqwf;
|
||
igl::min_quad_with_fixed_precompute(Q,b,Aeq,true,mqwf);
|
||
</code></pre>
|
||
|
||
<p>The output is a struct <code>mqwf</code> which contains the system matrix factorization
|
||
and is used during solving with arbitrary linear terms, known values, and
|
||
constraint in the right-hand sides:</p>
|
||
|
||
<pre><code class="cpp">igl::min_quad_with_fixed_solve(mqwf,B,bc,Beq,Z);
|
||
</code></pre>
|
||
|
||
<p>The output <code>Z</code> is a <span class="math">\(n \times 1\)</span> vector of solutions with fixed values
|
||
correctly placed to match the mesh vertices <code>V</code>.</p>
|
||
|
||
<h2 id="linearequalityconstraints">Linear equality constraints</h2>
|
||
|
||
<p>We saw above that <code>min_quad_with_fixed_*</code> in libigl provides a compact way to
|
||
solve general quadratic programs. Let’s consider another example, this time
|
||
with active linear equality constraints. Specifically let’s solve the
|
||
<code>bi-Laplace equation</code> or equivalently minimize the Laplace energy:</p>
|
||
|
||
<p><span class="math">\[\Delta^2 z = 0 \leftrightarrow \mathop{\text{minimize }}\limits_z \frac{1}{2}
|
||
\int\limits_S (\Delta z)^2 dA\]</span></p>
|
||
|
||
<p>subject to fixed value constraints and a linear equality constraint:</p>
|
||
|
||
<p><span class="math">\(z_{a} = 1, z_{b} = -1\)</span> and <span class="math">\(z_{c} = z_{d}\)</span>.</p>
|
||
|
||
<p>Notice that we can rewrite the last constraint in the familiar form from above:</p>
|
||
|
||
<p><span class="math">\(z_{c} - z_{d} = 0.\)</span></p>
|
||
|
||
<p>Now we can assembly <code>Aeq</code> as a <span class="math">\(1 \times n\)</span> sparse matrix with a coefficient
|
||
<span class="math">\(1\)</span> in the column corresponding to vertex <span class="math">\(c\)</span> and a <span class="math">\(-1\)</span> at <span class="math">\(d\)</span>. The right-hand
|
||
side <code>Beq</code> is simply zero.</p>
|
||
|
||
<p>Internally, <code>min_quad_with_fixed_*</code> solves using the Lagrange Multiplier
|
||
method. This method adds additional variables for each linear constraint (in
|
||
general a <span class="math">\(m \times 1\)</span> vector of variables <span class="math">\(\lambda\)</span>) and then solves the
|
||
saddle problem:</p>
|
||
|
||
<p><span class="math">\[\mathop{\text{find saddle }}_{\mathbf{z},\lambda}\, \frac{1}{2}\mathbf{z}^T \mathbf{Q} \mathbf{z} +
|
||
\mathbf{z}^T \mathbf{B} + \text{constant} + \lambda^T\left(\mathbf{A}_{eq}
|
||
\mathbf{z} - \mathbf{B}_{eq}\right)\]</span></p>
|
||
|
||
<p>This can be rewritten in a more familiar form by stacking <span class="math">\(\mathbf{z}\)</span> and
|
||
<span class="math">\(\lambda\)</span> into one <span class="math">\((m+n) \times 1\)</span> vector of unknowns:</p>
|
||
|
||
<p><span class="math">\[\mathop{\text{find saddle }}_{\mathbf{z},\lambda}\,
|
||
\frac{1}{2}
|
||
\left(
|
||
\mathbf{z}^T
|
||
\lambda^T
|
||
\right)
|
||
\left(
|
||
\begin{array}{cc}
|
||
\mathbf{Q} & \mathbf{A}_{eq}^T\\
|
||
\mathbf{A}_{eq} & 0
|
||
\end{array}
|
||
\right)
|
||
\left(
|
||
\begin{array}{c}
|
||
\mathbf{z}\\
|
||
\lambda
|
||
\end{array}
|
||
\right) +
|
||
\left(
|
||
\mathbf{z}^T
|
||
\lambda^T
|
||
\right)
|
||
\left(
|
||
\begin{array}{c}
|
||
\mathbf{B}\\
|
||
-\mathbf{B}_{eq}
|
||
\end{array}
|
||
\right)
|
||
+ \text{constant}\]</span></p>
|
||
|
||
<p>Differentiating with respect to <span class="math">\(\left( \mathbf{z}^T \lambda^T \right)\)</span> reveals
|
||
a linear system and we can solve for <span class="math">\(\mathbf{z}\)</span> and <span class="math">\(\lambda\)</span>. The only
|
||
difference from the straight quadratic <em>minimization</em> system, is that this
|
||
saddle problem system will not be positive definite. Thus, we must use a
|
||
different factorization technique (LDLT rather than LLT): libigl’s
|
||
<code>min_quad_with_fixed_precompute</code> automatically chooses the correct solver in
|
||
the presence of linear equality constraints (<a href="304_LinearEqualityConstraints/main.cpp">Example 304</a>).</p>
|
||
|
||
<figure>
|
||
<img src="images/cheburashka-biharmonic-leq.jpg" alt="The example LinearEqualityConstraints first solves with just fixed value
|
||
constraints (left: 1 and -1 on the left hand and foot respectively), then
|
||
solves with an additional linear equality constraint (right: points on right
|
||
hand and foot constrained to be equal)." />
|
||
<figcaption>The example <code>LinearEqualityConstraints</code> first solves with just fixed value
|
||
constraints (left: 1 and –1 on the left hand and foot respectively), then
|
||
solves with an additional linear equality constraint (right: points on right
|
||
hand and foot constrained to be equal).</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="quadraticprogramming">Quadratic programming</h2>
|
||
|
||
<p>We can generalize the quadratic optimization in the previous section even more
|
||
by allowing inequality constraints. Specifically box constraints (lower and
|
||
upper bounds):</p>
|
||
|
||
<p><span class="math">\(\mathbf{l} \le \mathbf{z} \le \mathbf{u},\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{l},\mathbf{u}\)</span> are <span class="math">\(n \times 1\)</span> vectors of lower and upper
|
||
bounds
|
||
and general linear inequality constraints:</p>
|
||
|
||
<p><span class="math">\(\mathbf{A}_{ieq} \mathbf{z} \le \mathbf{B}_{ieq},\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{A}_{ieq}\)</span> is a <span class="math">\(k \times n\)</span> matrix of linear coefficients and
|
||
<span class="math">\(\mathbf{B}_{ieq}\)</span> is a <span class="math">\(k \times 1\)</span> matrix of constraint right-hand sides.</p>
|
||
|
||
<p>Again, we are overly general as the box constraints could be written as
|
||
rows of the linear inequality constraints, but bounds appear frequently enough
|
||
to merit a dedicated api.</p>
|
||
|
||
<p>Libigl implements its own active set routine for solving <em>quadratric programs</em>
|
||
(QPs). This algorithm works by iteratively “activating” violated inequality
|
||
constraints by enforcing them as equalities and “deactivating” constraints
|
||
which are no longer needed.</p>
|
||
|
||
<p>After deciding which constraints are active at each iteration, the problem
|
||
reduces to a quadratic minimization subject to linear <em>equality</em> constraints,
|
||
and the method from the previous section is invoked. This is repeated until convergence.</p>
|
||
|
||
<p>Currently the implementation is efficient for box constraints and sparse
|
||
non-overlapping linear inequality constraints.</p>
|
||
|
||
<p>Unlike alternative interior-point methods, the active set method benefits from
|
||
a warm-start (initial guess for the solution vector <span class="math">\(\mathbf{z}\)</span>).</p>
|
||
|
||
<pre><code class="cpp">igl::active_set_params as;
|
||
// Z is optional initial guess and output
|
||
igl::active_set(Q,B,b,bc,Aeq,Beq,Aieq,Bieq,lx,ux,as,Z);
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/cheburashka-multiscale-biharmonic-kernels.jpg" alt=" Example 305 uses an active set solver to optimize
|
||
discrete biharmonic kernels at multiple scales
|
||
." />
|
||
<figcaption> <a href="305_QuadraticProgramming/main.cpp">Example 305</a> uses an active set solver to optimize
|
||
discrete biharmonic kernels <a class="citation" href="#fn:6" title="Jump to citation">[6]<span class="citekey" style="display:none">rustamov_2011</span></a> at multiple scales
|
||
.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="eigendecomposition">Eigen Decomposition</h2>
|
||
|
||
<p>Libigl has rudimentary support for extracting eigen pairs of a generalized
|
||
eigen value problem:</p>
|
||
|
||
<p><span class="math">\(Ax = \lambda B x\)</span></p>
|
||
|
||
<p>where <span class="math">\(A\)</span> is a sparse symmetric matrix and <span class="math">\(B\)</span> is a sparse positive definite
|
||
matrix. Most commonly in geometry processing, we let <span class="math">\(A=L\)</span> the cotangent
|
||
Laplacian and <span class="math">\(B=M\)</span> the per-vertex mass matrix (e.g. <a class="citation" href="#fn:7" title="Jump to citation">[7]<span class="citekey" style="display:none">vallet_2008</span></a>).
|
||
Typically applications will make use of the <em>low frequency</em> eigen modes.
|
||
Analagous to the Fourier decomposition, a function <span class="math">\(f\)</span> on a surface can be
|
||
represented via its spectral decomposition of the eigen modes of the
|
||
Laplace-Beltrami:</p>
|
||
|
||
<p><span class="math">\(f = \sum\limits_{i=1}^\infty a_i \phi_i\)</span></p>
|
||
|
||
<p>where each <span class="math">\(\phi_i\)</span> is an eigen function satisfying: <span class="math">\(\Delta \phi_i = \lambda_i
|
||
\phi_i\)</span> and <span class="math">\(a_i\)</span> are scalar coefficients. For a discrete triangle mesh, a
|
||
completely analogous decomposition exists, albeit with finite sum:</p>
|
||
|
||
<p><span class="math">\(\mathbf{f} = \sum\limits_{i=1}^n a_i \phi_i\)</span></p>
|
||
|
||
<p>where now a column vector of values at vertices <span class="math">\(\mathbf{f} \in \mathcal{R}^n\)</span>
|
||
specifies a piecewise linear function and <span class="math">\(\phi_i \in \mathcal{R}^n\)</span> is an
|
||
eigen vector satisfying:</p>
|
||
|
||
<p><span class="math">\(\mathbf{L} \phi_i = \lambda_i \mathbf{M} \phi_i\)</span>.</p>
|
||
|
||
<p>Note that Vallet & Levy <a class="citation" href="#fn:7" title="Jump to citation">[7]<span class="citekey" style="display:none">vallet_2008</span></a> propose solving a symmetrized
|
||
<em>standard</em> eigen problem <span class="math">\(\mathbf{M}^{-1/2}\mathbf{L}\mathbf{M}^{-1/2} \phi_i
|
||
= \lambda_i \phi_i\)</span>. Libigl implements a generalized eigen problem solver so
|
||
this unnecessary symmetrization can be avoided.</p>
|
||
|
||
<p>Often the sum above is <em>truncated</em> to the first <span class="math">\(k\)</span> eigen vectors. If the low
|
||
frequency modes are chosen, i.e. those corresponding to small <span class="math">\(\lambda_i\)</span>
|
||
values, then this truncation effectively <em>regularizes</em> <span class="math">\(\mathbf{f}\)</span> to smooth,
|
||
slowly changing functions over the mesh (e.g. <a class="citation" href="#fn:8" title="Jump to citation">[8]<span class="citekey" style="display:none">hildebrandt_2011</span></a>). Modal
|
||
analysis and model subspaces have been used frequently in real-time deformation
|
||
(e.g. <a class="citation" href="#fn:9" title="Jump to citation">[9]<span class="citekey" style="display:none">barbic_2005</span></a>).</p>
|
||
|
||
<p>In <a href="306_EigenDecomposition/main.cpp">Example 306</a>), the first 5 eigen vectors
|
||
of the discrete Laplace-Beltrami operator are computed and displayed in
|
||
pseudo-color atop the beetle. Eigen vectors are computed using <code>igl::eigs</code>
|
||
(mirroring MATLAB’s <code>eigs</code>). The 5 eigen vectors are placed into the columns
|
||
of <code>U</code> and the eigen values are placed into the entries of <code>S</code>:</p>
|
||
|
||
<pre><code class="cpp">SparseMatrix<double> L,M;
|
||
igl::cotmatrix(V,F,L);
|
||
igl::massmatrix(V,F,igl::MASSMATRIX_TYPE_DEFAULT,M);
|
||
Eigen::MatrixXd U;
|
||
Eigen::VectorXd S;
|
||
igl::eigs(L,M,5,igl::EIGS_TYPE_SM,U,S);
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/beetle-eigen-decomposition.gif" alt="(Example 306) Low frequency eigen vectors
|
||
of the discrete Laplace-Beltrami operator vary smoothly and slowly over the
|
||
Beetle." />
|
||
<figcaption>(<a href="306_EigenDecomposition/main.cpp">Example 306</a>) Low frequency eigen vectors
|
||
of the discrete Laplace-Beltrami operator vary smoothly and slowly over the
|
||
<em>Beetle</em>.</figcaption>
|
||
</figure>
|
||
|
||
<h1 id="chapter4:shapedeformation">Chapter 4: Shape deformation</h1>
|
||
|
||
<p>Modern mesh-based shape deformation methods satisfy user deformation
|
||
constraints at handles (selected vertices or regions on the mesh) and propagate
|
||
these handle deformations to the rest of shape <em>smoothly</em> and <em>without removing
|
||
or distorting details</em>. Libigl provides implementations of a variety of
|
||
state-of-the-art deformation techniques, ranging from quadratic mesh-based
|
||
energy minimizers, to skinning methods, to non-linear elasticity-inspired
|
||
techniques.</p>
|
||
|
||
<h2 id="biharmonicdeformation">Biharmonic deformation</h2>
|
||
|
||
<p>The period of research between 2000 and 2010 produced a collection of
|
||
techniques that cast the problem of handle-based shape deformation as a
|
||
quadratic energy minimization problem or equivalently the solution to a linear
|
||
partial differential equation.</p>
|
||
|
||
<p>There are many flavors of these techniques, but a prototypical subset are those
|
||
that consider solutions to the bi-Laplace equation, that is a biharmonic
|
||
function <a class="citation" href="#fn:10" title="Jump to citation">[10]<span class="citekey" style="display:none">botsch_2004</span></a>. This fourth-order PDE provides sufficient
|
||
flexibility in boundary conditions to ensure <span class="math">\(C^1\)</span> continuity at handle
|
||
constraints (in the limit under refinement) <a class="citation" href="#fn:11" title="Jump to citation">[11]<span class="citekey" style="display:none">jacobson_mixed_2010</span></a>.</p>
|
||
|
||
<h3 id="biharmonicsurfaces">Biharmonic surfaces</h3>
|
||
|
||
<p>Let us first begin our discussion of biharmonic <em>deformation</em>, by considering
|
||
biharmonic <em>surfaces</em>. We will casually define biharmonic surfaces as surface
|
||
whose <em>position functions</em> are biharmonic with respect to some initial
|
||
parameterization:</p>
|
||
|
||
<p><span class="math">\(\Delta^2 \mathbf{x}' = 0\)</span></p>
|
||
|
||
<p>and subject to some handle constraints, conceptualized as “boundary
|
||
conditions”:</p>
|
||
|
||
<p><span class="math">\(\mathbf{x}'_{b} = \mathbf{x}_{bc}.\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{x}'\)</span> is the unknown 3D position of a point on the surface. So we
|
||
are asking that the bi-Laplacian of each of spatial coordinate function to be
|
||
zero.</p>
|
||
|
||
<p>In libigl, one can solve a biharmonic problem with <code>igl::harmonic</code>
|
||
and setting <span class="math">\(k=2\)</span> (<em>bi</em>-harmonic):</p>
|
||
|
||
<pre><code class="cpp">// U_bc contains deformation of boundary vertices b
|
||
igl::harmonic(V,F,b,U_bc,2,U);
|
||
</code></pre>
|
||
|
||
<p>This produces a smooth surface that interpolates the handle constraints, but all
|
||
original details on the surface will be <em>smoothed away</em>. Most obviously, if the
|
||
original surface is not already biharmonic, then giving all handles the
|
||
identity deformation (keeping them at their rest positions) will <strong>not</strong>
|
||
reproduce the original surface. Rather, the result will be the biharmonic
|
||
surface that does interpolate those handle positions.</p>
|
||
|
||
<p>Thus, we may conclude that this is not an intuitive technique for shape
|
||
deformation.</p>
|
||
|
||
<h3 id="biharmonicdeformationfields">Biharmonic deformation fields</h3>
|
||
|
||
<p>Now we know that one useful property for a deformation technique is “rest pose
|
||
reproduction”: applying no deformation to the handles should apply no
|
||
deformation to the shape.</p>
|
||
|
||
<p>To guarantee this by construction we can work with <em>deformation fields</em> (ie.
|
||
displacements)
|
||
<span class="math">\(\mathbf{d}\)</span> rather
|
||
than directly with positions <span class="math">\(\mathbf{x}\)</span>. Then the deformed positions can be
|
||
recovered as</p>
|
||
|
||
<p><span class="math">\(\mathbf{x}' = \mathbf{x}+\mathbf{d}.\)</span></p>
|
||
|
||
<p>A smooth deformation field <span class="math">\(\mathbf{d}\)</span> which interpolates the deformation
|
||
fields of the handle constraints will impose a smooth deformed shape
|
||
<span class="math">\(\mathbf{x}'\)</span>. Naturally, we consider <em>biharmonic deformation fields</em>:</p>
|
||
|
||
<p><span class="math">\(\Delta^2 \mathbf{d} = 0\)</span></p>
|
||
|
||
<p>subject to the same handle constraints, but rewritten in terms of their implied
|
||
deformation field at the boundary (handles):</p>
|
||
|
||
<p><span class="math">\(\mathbf{d}_b = \mathbf{x}_{bc} - \mathbf{x}_b.\)</span></p>
|
||
|
||
<p>Again we can use <code>igl::harmonic</code> with <span class="math">\(k=2\)</span>, but this time solve for the
|
||
deformation field and then recover the deformed positions:</p>
|
||
|
||
<pre><code class="cpp">// U_bc contains deformation of boundary vertices b
|
||
D_bc = U_bc - igl::slice(V,b,1);
|
||
igl::harmonic(V,F,b,D_bc,2,D);
|
||
U = V+D;
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/max-biharmonic.jpg" alt="The BiharmonicDeformation example deforms a statues head as a biharmonic
|
||
surface (top) and using a biharmonic displacements
|
||
(bottom)." />
|
||
<figcaption>The <a href="401_BiharmonicDeformation/main.cpp">BiharmonicDeformation</a> example deforms a statue’s head as a <em>biharmonic
|
||
surface</em> (top) and using a <em>biharmonic displacements</em>
|
||
(bottom).</figcaption>
|
||
</figure>
|
||
|
||
<h4 id="relationshiptodifferentialcoordinatesandlaplaciansurfaceediting">Relationship to “differential coordinates” and Laplacian surface editing</h4>
|
||
|
||
<p>Biharmonic functions (whether positions or displacements) are solutions to the
|
||
bi-Laplace equation, but also minimizers of the “Laplacian energy”. For
|
||
example, for displacements <span class="math">\(\mathbf{d}\)</span>, the energy reads</p>
|
||
|
||
<p><span class="math">\(\int\limits_S \|\Delta \mathbf{d}\|^2 dA,\)</span></p>
|
||
|
||
<p>where we define <span class="math">\(\Delta \mathbf{d}\)</span> to simply apply the Laplacian
|
||
coordinate-wise.</p>
|
||
|
||
<p>By linearity of the Laplace(-Beltrami) operator we can reexpress this energy in
|
||
terms of the original positions <span class="math">\(\mathbf{x}\)</span> and the unknown positions
|
||
<span class="math">\(\mathbf{x}' = \mathbf{x} - \mathbf{d}\)</span>:</p>
|
||
|
||
<p><span class="math">\(\int\limits_S \|\Delta (\mathbf{x}' - \mathbf{x})\|^2 dA = \int\limits_S
|
||
\|\Delta \mathbf{x}' - \Delta \mathbf{x})\|^2 dA.\)</span></p>
|
||
|
||
<p>In the early work of Sorkine et al., the quantities <span class="math">\(\Delta \mathbf{x}'\)</span> and
|
||
<span class="math">\(\Delta \mathbf{x}\)</span> were dubbed “differential coordinates” <a class="citation" href="#fn:12" title="Jump to citation">[12]<span class="citekey" style="display:none">sorkine_2004</span></a>.
|
||
Their deformations (without linearized rotations) is thus equivalent to
|
||
biharmonic deformation fields.</p>
|
||
|
||
<h2 id="polyharmonicdeformation">Polyharmonic deformation</h2>
|
||
|
||
<p>We can generalize biharmonic deformation by considering different powers of
|
||
the Laplacian, resulting in a series of PDEs of the form:</p>
|
||
|
||
<p><span class="math">\(\Delta^k \mathbf{d} = 0.\)</span></p>
|
||
|
||
<p>with <span class="math">\(k\in{1,2,3,\dots}\)</span>. The choice of <span class="math">\(k\)</span> determines the level of continuity
|
||
at the handles. In particular, <span class="math">\(k=1\)</span> implies <span class="math">\(C^0\)</span> at the boundary, <span class="math">\(k=2\)</span>
|
||
implies <span class="math">\(C^1\)</span>, <span class="math">\(k=3\)</span> implies <span class="math">\(C^2\)</span> and in general <span class="math">\(k\)</span> implies <span class="math">\(C^{k-1}\)</span>.</p>
|
||
|
||
<pre><code class="cpp">int k = 2;// or 1,3,4,...
|
||
igl::harmonic(V,F,b,bc,k,Z);
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/bump-k-harmonic.jpg" alt="The PolyharmonicDeformation example deforms a flat domain (left) into a bump as a
|
||
solution to various $k$-harmonic PDEs." />
|
||
<figcaption>The <a href="402_PolyharmonicDeformation/main.cpp">PolyharmonicDeformation</a> example deforms a flat domain (left) into a bump as a
|
||
solution to various <span class="math">\(k\)</span>-harmonic PDEs.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="boundedbiharmonicweights">Bounded biharmonic weights</h2>
|
||
|
||
<p>In computer animation, shape deformation is often referred to as “skinning”.
|
||
Constraints are posed as relative rotations of internal rigid “bones” inside a
|
||
character. The deformation method, or skinning method, determines how the
|
||
surface of the character (i.e. its skin) should move as a function of the bone
|
||
rotations.</p>
|
||
|
||
<p>The most popular technique is linear blend skinning. Each point on the shape
|
||
computes its new location as a linear combination of bone transformations:</p>
|
||
|
||
<p><span class="math">\(\mathbf{x}' = \sum\limits_{i = 1}^m w_i(\mathbf{x}) \mathbf{T}_i
|
||
\left(\begin{array}{c}\mathbf{x}_i\\1\end{array}\right),\)</span></p>
|
||
|
||
<p>where <span class="math">\(w_i(\mathbf{x})\)</span> is the scalar <em>weight function</em> of the ith bone evaluated at
|
||
<span class="math">\(\mathbf{x}\)</span> and <span class="math">\(\mathbf{T}_i\)</span> is the bone transformation as a <span class="math">\(4 \times 3\)</span>
|
||
matrix.</p>
|
||
|
||
<p>This formula is embarassingly parallel (computation at one point does not
|
||
depend on shared data need by computation at another point). It is often
|
||
implemented as a vertex shader. The weights and rest positions for each vertex
|
||
are sent as vertex shader <em>attributes</em> and bone transformations are sent as
|
||
<em>uniforms</em>. Then vertices are transformed within the vertex shader, just in
|
||
time for rendering.</p>
|
||
|
||
<p>As the skinning formula is linear (hence its name), we can write it as matrix
|
||
multiplication:</p>
|
||
|
||
<p><span class="math">\(\mathbf{X}' = \mathbf{M} \mathbf{T},\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{X}'\)</span> is <span class="math">\(n \times 3\)</span> stack of deformed positions as row
|
||
vectors, <span class="math">\(\mathbf{M}\)</span> is a <span class="math">\(n \times m\cdot dim\)</span> matrix containing weights and
|
||
rest positions and <span class="math">\(\mathbf{T}\)</span> is a <span class="math">\(m\cdot (dim+1) \times dim\)</span> stack of
|
||
transposed bone transformations.</p>
|
||
|
||
<p>Traditionally, the weight functions <span class="math">\(w_j\)</span> are painted manually by skilled
|
||
rigging professionals. Modern techniques now exist to compute weight functions
|
||
automatically given the shape and a description of the skeleton (or in general
|
||
any handle structure such as a cage, collection of points, selected regions,
|
||
etc.).</p>
|
||
|
||
<p>Bounded biharmonic weights are one such technique that casts weight computation
|
||
as a constrained optimization problem <a class="citation" href="#fn:13" title="Jump to citation">[13]<span class="citekey" style="display:none">jacobson_2011</span></a>. The weights enforce
|
||
smoothness by minimizing the familiar Laplacian energy:</p>
|
||
|
||
<p><span class="math">\(\sum\limits_{i = 1}^m \int_S (\Delta w_i)^2 dA\)</span></p>
|
||
|
||
<p>subject to constraints which enforce interpolation of handle constraints:</p>
|
||
|
||
<p><span class="math">\(w_i(\mathbf{x}) = \begin{cases} 1 & \text{ if } \mathbf{x} \in H_i\\ 0 &
|
||
\text{ otherwise } \end{cases},\)</span></p>
|
||
|
||
<p>where <span class="math">\(H_i\)</span> is the ith handle, and constraints which enforce non-negativity,
|
||
parition of unity and encourage sparsity:</p>
|
||
|
||
<p><span class="math">\(0\le w_i \le 1\)</span> and <span class="math">\(\sum\limits_{i=1}^m w_i = 1.\)</span></p>
|
||
|
||
<p>This is a quadratic programming problem and libigl solves it using its active
|
||
set solver or by calling out to <a href="http://www.mosek.com">Mosek</a>.</p>
|
||
|
||
<figure>
|
||
<img src="images/hand-bbw.jpg" alt="The example BoundedBiharmonicWeights computes weights for a tetrahedral
|
||
mesh given a skeleton (top) and then animates a linear blend skinning
|
||
deformation (bottom)." />
|
||
<figcaption>The example <a href="403_BoundedBiharmonicWeights/main.cpp">BoundedBiharmonicWeights</a> computes weights for a tetrahedral
|
||
mesh given a skeleton (top) and then animates a linear blend skinning
|
||
deformation (bottom).</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="dualquaternionskinning">Dual quaternion skinning</h2>
|
||
|
||
<p>Even with high quality weights, linear blend skinning is limited. In
|
||
particular, it suffers from known artifacts stemming from blending rotations as
|
||
as matrices: a weight combination of rotation matrices is not necessarily a
|
||
rotation. Consider an equal blend between rotating by <span class="math">\(-\pi/2\)</span> and by <span class="math">\(\pi/2\)</span>
|
||
about the <span class="math">\(z\)</span>-axis. Intuitively one might expect to get the identity matrix,
|
||
but instead the blend is a degenerate matrix scaling the <span class="math">\(x\)</span> and <span class="math">\(y\)</span>
|
||
coordinates by zero:</p>
|
||
|
||
<p><span class="math">\(0.5\left(\begin{array}{ccc}0&-1&0\\1&0&0\\0&0&1\end{array}\right)+
|
||
0.5\left(\begin{array}{ccc}0&1&0\\-1&0&0\\0&0&1\end{array}\right)=
|
||
\left(\begin{array}{ccc}0&0&0\\0&0&0\\0&0&1\end{array}\right)\)</span></p>
|
||
|
||
<p>In practice, this means the shape shrinks and collapses in regions where bone
|
||
weights overlap: near joints.</p>
|
||
|
||
<p>Dual quaternion skinning presents a solution <a class="citation" href="#fn:14" title="Jump to citation">[14]<span class="citekey" style="display:none">kavan_2008</span></a>. This method
|
||
represents rigid transformations as a pair of unit quaternions,
|
||
<span class="math">\(\hat{\mathbf{q}}\)</span>. The linear blend skinning formula is replaced with a
|
||
linear blend of dual quaternions:</p>
|
||
|
||
<p><span class="math">\(\mathbf{x}' =
|
||
\cfrac{\sum\limits_{i=1}^m w_i(\mathbf{x})\hat{\mathbf{q}_i}}
|
||
{\left\|\sum\limits_{i=1}^m w_i(\mathbf{x})\hat{\mathbf{q}_i}\right\|}
|
||
\mathbf{x},\)</span></p>
|
||
|
||
<p>where <span class="math">\(\hat{\mathbf{q}_i}\)</span> is the dual quaternion representation of the rigid
|
||
transformation of bone <span class="math">\(i\)</span>. The normalization forces the result of the linear
|
||
blending to again be a unit dual quaternion and thus also a rigid
|
||
transformation.</p>
|
||
|
||
<p>Like linear blend skinning, dual quaternion skinning is best performed in the
|
||
vertex shader. The only difference being that bone transformations are sent as
|
||
dual quaternions rather than affine transformation matrices. Libigl supports
|
||
CPU-side dual quaternion skinning with the <code>igl::dqs</code> function, which takes a
|
||
more traditional representation of rigid transformations as input and
|
||
internally converts to the dual quaternion representation before blending:</p>
|
||
|
||
<pre><code class="cpp">// vQ is a list of rotations as quaternions
|
||
// vT is a list of translations
|
||
igl::dqs(V,W,vQ,vT,U);
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/arm-dqs.jpg" alt="The example DualQuaternionSkinning compares linear blend skinning (top) to dual
|
||
quaternion skinning (bottom), highlighting LBSs candy wrapper effect (middle)
|
||
and joint collapse (right)." />
|
||
<figcaption>The example <a href="404_DualQuaternionSkinning/main.cpp">DualQuaternionSkinning</a> compares linear blend skinning (top) to dual
|
||
quaternion skinning (bottom), highlighting LBS’s candy wrapper effect (middle)
|
||
and joint collapse (right).</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="as-rigid-as-possible">As-rigid-as-possible</h2>
|
||
|
||
<p>Skinning and other linear methods for deformation are inherently limited.
|
||
Difficult arises especially when large rotations are imposed by the handle
|
||
constraints.</p>
|
||
|
||
<p>In the context of energy-minimization approaches, the problem stems from
|
||
comparing positions (our displacements) in the coordinate frame of the
|
||
undeformed shape. These quadratic energies are at best invariant to global
|
||
rotations of the entire shape, but not smoothly varying local rotations. Thus
|
||
linear techniques will not produce non-trivial bending and twisting.</p>
|
||
|
||
<p>Furthermore, when considering solid shapes (e.g. discretized with tetrahedral
|
||
meshes) linear methods struggle to maintain local volume, and they often suffer from
|
||
shrinking and bulging artifacts.</p>
|
||
|
||
<p>Non-linear deformation techniques present a solution to these problems.
|
||
They work by comparing the deformation of a mesh
|
||
vertex to its rest position <em>rotated</em> to a new coordinate frame which best
|
||
matches the deformation. The non-linearity stems from the mutual dependence of
|
||
the deformation and the best-fit rotation. These techniques are often labeled
|
||
“as-rigid-as-possible” as they penalize the sum of all local deformations’
|
||
deviations from rotations.</p>
|
||
|
||
<p>To arrive at such an energy, let’s consider a simple per-triangle energy:</p>
|
||
|
||
<p><span class="math">\(E_\text{linear}(\mathbf{X}') = \sum\limits_{t \in T} a_t \sum\limits_{\{i,j\}
|
||
\in t} w_{ij} \left\|
|
||
\left(\mathbf{x}'_i - \mathbf{x}'_j\right) -
|
||
\left(\mathbf{x}_i - \mathbf{x}_j\right)\right\|^2\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{X}'\)</span> are the mesh’s unknown deformed vertex positions, <span class="math">\(t\)</span> is a
|
||
triangle in a list of triangles <span class="math">\(T\)</span>, <span class="math">\(a_t\)</span> is the area of triangle <span class="math">\(t\)</span> and
|
||
<span class="math">\(\{i,j\}\)</span> is an edge in triangle <span class="math">\(t\)</span>. Thus, this energy measures the norm of
|
||
change between an edge vector in the original mesh <span class="math">\(\left(\mathbf{x}_i -
|
||
\mathbf{x}_j\right)\)</span> and the unknown mesh <span class="math">\(\left(\mathbf{x}'_i -
|
||
\mathbf{x}'_j\right)\)</span>.</p>
|
||
|
||
<p>This energy is <strong>not</strong> rotation invariant. If we rotate the mesh by 90 degrees
|
||
the change in edge vectors not aligned with the axis of rotation will be large,
|
||
despite the overall deformation being perfectly rigid.</p>
|
||
|
||
<p>So, the “as-rigid-as-possible” solution is to append auxiliary variables
|
||
<span class="math">\(\mathbf{R}_t\)</span>
|
||
for each triangle <span class="math">\(t\)</span> which are constrained to be rotations. Then the energy is
|
||
rewritten, this time comparing deformed edge vectors to their rotated rest
|
||
counterparts:</p>
|
||
|
||
<p><span class="math">\(E_\text{arap}(\mathbf{X}',\{\mathbf{R}_1,\dots,\mathbf{R}_{|T|}\}) = \sum\limits_{t \in T} a_t \sum\limits_{\{i,j\}
|
||
\in t} w_{ij} \left\|
|
||
\left(\mathbf{x}'_i - \mathbf{x}'_j\right)-
|
||
\mathbf{R}_t\left(\mathbf{x}_i - \mathbf{x}_j\right)\right\|^2.\)</span></p>
|
||
|
||
<p>The separation into the primary vertex position variables <span class="math">\(\mathbf{X}'\)</span> and the
|
||
rotations <span class="math">\(\{\mathbf{R}_1,\dots,\mathbf{R}_{|T|}\}\)</span> lead to strategy for
|
||
optimization, too. If the rotations <span class="math">\(\{\mathbf{R}_1,\dots,\mathbf{R}_{|T|}\}\)</span>
|
||
are held fixed then the energy is quadratic in the remaining variables
|
||
<span class="math">\(\mathbf{X}'\)</span> and can be optimized by solving a (sparse) global linear system.
|
||
Alternatively, if <span class="math">\(\mathbf{X}'\)</span> are held fixed then each rotation is the
|
||
solution to a localized <em>Procrustes</em> problem (found via <span class="math">\(3 \times 3\)</span> SVD or
|
||
polar decompostion). These two steps—local and global—each weakly decrease
|
||
the energy, thus we may safely iterate them until convergence.</p>
|
||
|
||
<p>The different flavors of “as-rigid-as-possible” depend on the dimension and
|
||
codimension of the domain and the edge-sets <span class="math">\(T\)</span>. The proposed surface
|
||
manipulation technique by Sorkine and Alexa <a class="citation" href="#fn:15" title="Jump to citation">[15]<span class="citekey" style="display:none">sorkine_2007</span></a>, considers <span class="math">\(T\)</span> to
|
||
be the set of sets of edges emanating from each vertex (spokes). Later, Chao et
|
||
al. derived the relationship between “as-rigid-as-possible” mesh energies and
|
||
co-rotational elasticity considering 0-codimension elements as edge-sets:
|
||
triangles in 2D and tetrahedra in 3D <a class="citation" href="#fn:16" title="Jump to citation">[16]<span class="citekey" style="display:none">chao_2010</span></a>. They also showed how
|
||
Sorkine and Alexa’s edge-sets are not a discretization of a continuous energy,
|
||
proposing instead edge-sets for surfaces containing all edges of elements
|
||
incident on a vertex (spokes and rims). They show that this amounts to
|
||
measuring bending, albeit in a discretization-dependent way.</p>
|
||
|
||
<p>Libigl, supports these common flavors. Selecting one is a matter of setting the
|
||
energy type before the precompuation phase:</p>
|
||
|
||
<pre><code class="cpp">igl::ARAPData arap_data;
|
||
arap_data.energy = igl::ARAP_ENERGY_TYPE_SPOKES;
|
||
//arap_data.energy = igl::ARAP_ENERGY_TYPE_SPOKES_AND_RIMS;
|
||
//arap_data.energy = igl::ARAP_ENERGY_TYPE_ELEMENTS; //triangles or tets
|
||
igl::arap_precomputation(V,F,dim,b,arap_data);
|
||
</code></pre>
|
||
|
||
<p>Just like <code>igl::min_quad_with_fixed_*</code>, this precomputation phase only depends
|
||
on the mesh, fixed vertex indices <code>b</code> and the energy parameters. To solve with
|
||
certain constraints on the positions of vertices in <code>b</code>, we may call:</p>
|
||
|
||
<pre><code class="cpp">igl::arap_solve(bc,arap_data,U);
|
||
</code></pre>
|
||
|
||
<p>which uses <code>U</code> as an initial guess and then computes the solution into it.</p>
|
||
|
||
<p>Libigl’s implementation of as-rigid-as-possible deformation takes advantage of
|
||
the highly optimized singular value decomposition code from McAdams et al.
|
||
<a class="citation" href="#fn:17" title="Jump to citation">[17]<span class="citekey" style="display:none">mcadams_2011</span></a> which leverages SSE intrinsics.</p>
|
||
|
||
<figure>
|
||
<img src="images/decimated-knight-arap.jpg" alt="The example AsRigidAsPossible deforms a surface as if it were made of an
|
||
elastic material" />
|
||
<figcaption>The example <a href="405_AsRigidAsPossible/main.cpp">AsRigidAsPossible</a> deforms a surface as if it were made of an
|
||
elastic material</figcaption>
|
||
</figure>
|
||
|
||
<p>The concept of local rigidity will be revisited shortly in the context of
|
||
surface parameterization.</p>
|
||
|
||
<h2 id="fastautomaticskinningtransformations">Fast automatic skinning transformations</h2>
|
||
|
||
<p>Non-linear optimization is, unsurprisingly, slower than its linear cousins. In
|
||
the case of the as-rigid-as-possible optimization, the bottleneck is typically
|
||
the large number of polar decompositions necessary to recover best fit
|
||
rotations for each edge-set (i.e. for each triangle, tetrahedron, or vertex
|
||
cell). Even if this code is optimized, the number of primary degrees of freedom
|
||
is tied to the discretization level, despite the deformations’ low frequency
|
||
behavior.</p>
|
||
|
||
<p>This invites two routes toward fast non-linear optimization. First, is it
|
||
necessary (or even advantageous) to find so many best-fit rotations? Second,
|
||
can we reduce the degrees of freedom to better reflect the frequency of the
|
||
desired deformations.</p>
|
||
|
||
<p>Taken in turn, these optimizations culminate in a method which optimizes over
|
||
the space of linear blend skinning deformations spanned by high-quality weights
|
||
(i.e. manually painted ones or bounded biharmonic weights). This space is a
|
||
low-dimensional subspace of all possible mesh deformations, captured by writing
|
||
linear blend skinning in matrix form:</p>
|
||
|
||
<p><span class="math">\(\mathbf{X}' = \mathbf{M}\mathbf{T}\)</span></p>
|
||
|
||
<p>where the mesh vertex positions in the <span class="math">\(n \times 3\)</span> matrix <span class="math">\(\mathbf{X}'\)</span> are
|
||
replaced by a linear combination of a small number of degrees of freedom in the
|
||
<span class="math">\((3+1)m \times 3\)</span> stack of transposed “handle” transformations. Swapping in
|
||
<span class="math">\(\mathbf{M}\mathbf{T}\)</span> for <span class="math">\(\mathbf{X}'\)</span> in the ARAP energies above immediately
|
||
sees performance gains during the global solve step as <span class="math">\(m << n\)</span>.</p>
|
||
|
||
<p>The complexity of the local step—fitting rotations—is still bound
|
||
to the original mesh discretization. However, if the skinning is well behaved,
|
||
we can make the assumption that places on the shape with similar skinning
|
||
weights will deform similarly and thus imply similar best-fit rotations.
|
||
Therefore, we cluster edge-sets according to their representation in
|
||
<em>weight-space</em>: where a vertex <span class="math">\(\mathbf{x}\)</span> takes the coordinates
|
||
<span class="math">\([w_1(\mathbf{x}),w_2(\mathbf{x}),\dots,w_m(\mathbf{x})]\)</span>. The number of
|
||
clustered edge-sets show diminishing returns on the deformation quality so we
|
||
may choose a small number of clusters, proportional to the number of skinning
|
||
weight functions (rather than the number of discrete mesh vertices).</p>
|
||
|
||
<p>This proposed deformation model <a class="citation" href="#fn:18" title="Jump to citation">[18]<span class="citekey" style="display:none">jacobson_2012</span></a>, can simultaneously be seen as a
|
||
fast, subspace optimization for ARAP and as an automatic method for finding
|
||
<em>the best</em> skinning transformation degrees of freedom.</p>
|
||
|
||
<p>A variety of user interfaces are supported via linear equality constraints on
|
||
the skinning transformations associated with handles. To fix a transformation
|
||
entirely we simply add the constraint:</p>
|
||
|
||
<p><span class="math">\(\left(\begin{array}{cccc}
|
||
1 & 0 & 0 & 0\\
|
||
0 & 1 & 0 & 0\\
|
||
0 & 0 & 1 & 0\\
|
||
0 & 0 & 0 & 1\end{array}\right)
|
||
\mathbf{T}_i^T = \hat{\mathbf{T}}_i^T,\)</span></p>
|
||
|
||
<p>where <span class="math">\(\hat{\mathbf{T}}_i^T\)</span> is the <span class="math">\((3+1) \times 3\)</span> transposed fixed
|
||
transformation for handle <span class="math">\(i\)</span>.</p>
|
||
|
||
<p>To fix only the origin of a handle, we add a constraint requiring the
|
||
transformation to interpolate a point in space (typically the centroid of all
|
||
points with <span class="math">\(w_i = 1\)</span>:</p>
|
||
|
||
<p><span class="math">\(\mathbf{c}'^T\mathbf{T}_i^T = \mathbf{c}^T,\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{c}^T\)</span> is the <span class="math">\(1 \times (3+1)\)</span> position of the point at rest in
|
||
transposed homogeneous coordinates, and <span class="math">\(\mathbf{c}'^T\)</span> the point given by the
|
||
user.</p>
|
||
|
||
<p>We can similarly fix just the linear part of the transformation at a handle,
|
||
freeing the translation component (producing a “chickenhead” effect):</p>
|
||
|
||
<p><span class="math">\(\left(\begin{array}{cccc}
|
||
1&0&0&0\\
|
||
0&1&0&0\\
|
||
0&0&1&0\end{array}\right)
|
||
\mathbf{T}_i^T = \hat{\mathbf{L}}_i^T,\)</span></p>
|
||
|
||
<p>where <span class="math">\(\hat{\mathbf{L}}_i^T\)</span> is the fixed <span class="math">\(3 \times 3\)</span> linear part of the
|
||
transformation at handle <span class="math">\(i\)</span>.</p>
|
||
|
||
<p>And lastly we can allow the user to entirely <em>free</em> the transformation’s
|
||
degrees of freedom, delegating the optimization to find the best possible
|
||
values for all elements. To do this, we simply abstain from adding a
|
||
corresponding constraint.</p>
|
||
|
||
<h3 id="arapwithgroupededge-sets">ARAP with grouped edge-sets</h3>
|
||
|
||
<p>Being a subspace method, an immediate disadvantage is the reduced degrees of
|
||
freedom. This brings performance, but in some situations limits behavior too
|
||
much. In such cases one can use the skinning subspace to build an effective
|
||
clustering of rotation edge-sets for a traditional ARAP optimization: forgoing
|
||
the subspace substitution. This has an two-fold effect. The cost of the
|
||
rotation fitting, local step drastically reduces, and the deformations are
|
||
“regularized” according the clusters. From a high level point of view, if the
|
||
clusters are derived from skinning weights, then they will discourage bending,
|
||
especially along isolines of the weight functions. If handles are not known in
|
||
advance, one could also cluster according to a “geodesic embedding” like the
|
||
biharmonic distance embedding.</p>
|
||
|
||
<p>In this light, we can think of the “spokes+rims” style surface ARAP as a (slight and
|
||
redundant) clustering of the per-triangle edge-sets.</p>
|
||
|
||
<figure>
|
||
<img src="images/armadillo-fast.jpg" alt="The example FastAutomaticSkinningTransformations compares a full (slow)
|
||
ARAP deformation on a detailed shape (left of middle), to ARAP with grouped
|
||
rotation edge sets (right of middle), to the very fast subpsace method
|
||
(right)." />
|
||
<figcaption>The example <a href="406_FastAutomaticSkinningTransformations/main.cpp">FastAutomaticSkinningTransformations</a> compares a full (slow)
|
||
ARAP deformation on a detailed shape (left of middle), to ARAP with grouped
|
||
rotation edge sets (right of middle), to the very fast subpsace method
|
||
(right).</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="biharmoniccoordinates">Biharmonic Coordinates</h2>
|
||
|
||
<p>Linear blend skinning (as <a href="#boundedbiharmonicweights">above</a>) deforms a mesh by
|
||
propogating <em>full affine transformations</em> at handles (bones, points, regions,
|
||
etc.) to the rest of the shape via weights. Another deformation framework,
|
||
called “generalized barycentric coordinates”, is a special case of linear blend
|
||
skinning <a class="citation" href="#fn:19" title="Jump to citation">[19]<span class="citekey" style="display:none">jacobson_skinning_course_2014</span></a>: transformations are restricted to
|
||
<em>pure translations</em> and weights are required to retain <em>affine precision</em>. This
|
||
latter requirement means that we can write the rest-position of any vertex in
|
||
the mesh as the weighted combination of the control handle locations:</p>
|
||
|
||
<p><span class="math">\(\mathbf{x} = \sum\limits_{i=1}^m w_i(\mathbf{x}) * \mathbf{c}_i,\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{c}_i\)</span> is the rest position of the $i$th control point. This
|
||
simplifies the deformation formula at run-time. We can simply take the new
|
||
position of each point of the shape to be the weighted combination of the
|
||
<em>translated</em> control point positions:</p>
|
||
|
||
<p><span class="math">\(\mathbf{x}' = \sum\limits_{i=1}^m w_i(\mathbf{x}) * \mathbf{c}_i'.\)</span></p>
|
||
|
||
<p>There are <em>many</em> different flavors of “generalized barycentric coordinates”
|
||
(see table in “Automatic Methods” section,
|
||
<a class="citation" href="#fn:19" title="Jump to citation">[19]<span class="citekey" style="display:none">jacobson_skinning_course_2014</span></a>). The vague goal of “generalized barycentric
|
||
coordinates” is to capture as many properties of simplicial barycentric
|
||
coordinates (e.g. for triangles in 2D and tetrahedral in 3D) for larger sets of
|
||
points or polyhedra. Some generalized barycentric coordinates can be computed
|
||
in closed form; others require optimization-based precomputation. Nearly all
|
||
flavors require connectivity information describing how the control points form
|
||
a external polyhedron around the input shape: a cage. However, a recent
|
||
techinique does not require a cage <a class="citation" href="#fn:20" title="Jump to citation">[20]<span class="citekey" style="display:none">wang_bc_2015</span></a>. This method ensures
|
||
affine precision during optimization over weights of a smoothness energy with
|
||
affine functions in its kernel:</p>
|
||
|
||
<p><span class="math">\(\mathop{\text{min}}_\mathbf{W}\,\, \text{trace}(\frac{1}{2}\mathbf{W}^T \mathbf{A}
|
||
\mathbf{W}), \text{subject to: } \mathbf{C} = \mathbf{W}\mathbf{C}\)</span></p>
|
||
|
||
<p>subject to interpolation constraints at selected vertices. If <span class="math">\(\mathbf{A}\)</span> has
|
||
affine functions in its kernel—that is, if <span class="math">\(\mathbf{A}\mathbf{V} = 0\)</span>—then
|
||
the weights <span class="math">\(\mathbf{W}\)</span> will retain affine precision and we’ll have that:</p>
|
||
|
||
<p><span class="math">\(\mathbf{V} = \mathbf{W}\mathbf{C}\)</span></p>
|
||
|
||
<p>the matrix form of the equality above. The proposed way to define <span class="math">\(\mathbf{A}\)</span>
|
||
is to construct a matrix <span class="math">\(\mathbf{K}\)</span> that measures the Laplacian at all
|
||
interior vertices <em>and at all boundary vertices</em>. The <em>usual</em> definition of the
|
||
discrete Laplacian (e.g. what libigl returns from <code>igl::cotmatrix</code>), measures
|
||
the Laplacian of a function for interior vertices, but measures the Laplacian
|
||
of a function <em>minus</em> the normal derivative of a function for boundary
|
||
vertices. Thus, we can let:</p>
|
||
|
||
<p><span class="math">\(\mathbf{K} = \mathbf{L} + \mathbf{N}\)</span></p>
|
||
|
||
<p>where <span class="math">\(\mathbf{L}\)</span> is the <em>usual</em> Laplacian and <span class="math">\(\mathbf{N}\)</span> is matrix that
|
||
computes normal derivatives of a piecewise-linear function at boundary vertices
|
||
of a mesh. Then <span class="math">\(\mathbf{A}\)</span> is taken as quadratic form computing the square of
|
||
the integral-average of <span class="math">\(\mathbf{K}\)</span> applied to a function and integrated over
|
||
the mesh:</p>
|
||
|
||
<p><span class="math">\(\mathbf{A} = (\mathbf{M}^{-1}\mathbf{K})^2_\mathbf{M} = \mathbf{K}^T \mathbf{M}^{-1}
|
||
\mathbf{K}.\)</span></p>
|
||
|
||
<p>Since the Laplacian <span class="math">\(\mathbf{K}\)</span> is a second-order derivative it measures zero on affine
|
||
functions, thus <span class="math">\(\mathbf{A}\)</span> has affine functions in its null space. A short
|
||
derivation proves that this implies <span class="math">\(\mathbf{W}\)</span> will be affine precise (see
|
||
<a class="citation" href="#fn:20" title="Jump to citation">[20]<span class="citekey" style="display:none">wang_bc_2015</span></a>).</p>
|
||
|
||
<p>Minimizers of this “squared Laplacian” energy are in some sense <em>discrete
|
||
biharmonic functions</em>. Thus they’re dubbed “biharmonic coordinates” (not the
|
||
same as <em>bounded biharmonic weights</em>, which are <em>not</em> generalized barycentric
|
||
coordinates).</p>
|
||
|
||
<p>In libigl, one can compute biharmonic coordinates given a mesh <code>(V,F)</code> and a
|
||
list <code>S</code> of selected control points or control regions (which act like skinning
|
||
handles):</p>
|
||
|
||
<pre><code class="cpp">igl::biharmonic_coordinates(V,F,S,W);
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/octopus-biharmonic-coordinates-physics.gif" alt="(Example 407) shows a physics
|
||
simulation on a coarse orange mesh. The vertices of this mesh become control
|
||
points for a biharmonic coordinates deformation of the blue high-resolution
|
||
mesh." />
|
||
<figcaption>(<a href="407_BiharmonicCoordinates/main.cpp">Example 407</a>) shows a physics
|
||
simulation on a coarse orange mesh. The vertices of this mesh become control
|
||
points for a biharmonic coordinates deformation of the blue high-resolution
|
||
mesh.</figcaption>
|
||
</figure>
|
||
|
||
<h1 id="chapter5:parametrization">Chapter 5: Parametrization</h1>
|
||
|
||
<p>In computer graphics, we denote as surface parametrization a map from the
|
||
surface to <span class="math">\(\mathbf{R}^2\)</span>. It is usually encoded by a new set of 2D
|
||
coordinates for each vertex of the mesh (and possibly also by a new set of
|
||
faces in one to one correspondence with the faces of the original surface).
|
||
Note that
|
||
this definition is the <em>inverse</em> of the classical differential geometry
|
||
definition.</p>
|
||
|
||
<p>A parametrization has many applications, ranging from texture mapping to
|
||
surface remeshing. Many algorithms have been proposed, and they can be broadly
|
||
divided in four families:</p>
|
||
|
||
<ol>
|
||
<li><p><strong>Single patch, fixed boundary</strong>: these algorithm can parametrize a
|
||
disk-like part of the surface given fixed 2D positions for its boundary. These
|
||
algorithms are efficient and simple, but they usually produce high-distortion maps due to the fixed boundary.</p></li>
|
||
<li><p><strong>Single patch, free boundary:</strong> these algorithms let the boundary
|
||
deform freely, greatly reducing the map distortion. Care should be taken to
|
||
prevent the border to self-intersect.</p></li>
|
||
<li><p><strong>Global parametrization</strong>: these algorithms work on meshes with arbitrary
|
||
genus. They initially cut the mesh in multiple patches that can be separately parametrized. The generated maps are discontinuous on the cuts (often referred as <em>seams</em>).</p></li>
|
||
<li><p><strong>Global seamless parametrization</strong>: these are global parametrization algorithm that hides the seams, making the parametrization “continuous”, under specific assumptions that we will discuss later.</p></li>
|
||
</ol>
|
||
|
||
<h2 id="harmonicparametrization"><a href="#harmonicparametrization">Harmonic parametrization</a></h2>
|
||
|
||
<p>Harmonic parametrization <a class="citation" href="#fn:21" title="Jump to citation">[21]<span class="citekey" style="display:none">eck_2005</span></a> is a single patch, fixed boundary parametrization
|
||
algorithm that computes the 2D coordinates of the flattened mesh as two
|
||
harmonic functions.</p>
|
||
|
||
<p>The algorithm is divided in 3 steps:</p>
|
||
|
||
<ol>
|
||
<li>Detect of the boundary vertices</li>
|
||
<li>Map the boundary vertices to a circle</li>
|
||
<li>Compute two harmonic functions (one for u and one for the v coordinate). The harmonic functions use the fixed vertices on the circle as boundary constraints.</li>
|
||
</ol>
|
||
|
||
<p>The algorithm can be coded using libigl as follows:</p>
|
||
|
||
<pre><code class="cpp">Eigen::VectorXi bnd;
|
||
igl::boundary_loop(V,F,bnd);
|
||
|
||
Eigen::MatrixXd bnd_uv;
|
||
igl::map_vertices_to_circle(V,bnd,bnd_uv);
|
||
|
||
igl::harmonic(V,F,bnd,bnd_uv,1,V_uv);
|
||
</code></pre>
|
||
|
||
<p>where <code>bnd</code> contains the indices of the boundary vertices, bnd_uv their position on the UV plane, and “1” denotes that we want to compute an harmonic function (2 will be for biharmonic, 3 for triharmonic, etc.). Note that each of the three
|
||
functions is designed to be reusable in other parametrization algorithms.</p>
|
||
|
||
<p>A UV parametrization can be visualized in the viewer with:</p>
|
||
|
||
<pre><code class="cpp">viewer.data.set_uv(V_uv);
|
||
</code></pre>
|
||
|
||
<p>The UV coordinates are then used to apply a procedural checkerboard texture to the
|
||
mesh (<a href="501_HarmonicParam/main.cpp">Example 501</a>).</p>
|
||
|
||
<figure>
|
||
<img src="images/501_HarmonicParam.png" alt="(Example 501) Harmonic parametrization. (left)
|
||
mesh with texture, (right) UV parametrization with
|
||
texture" />
|
||
<figcaption>(<a href="501_HarmonicParam/main.cpp">Example 501</a>) Harmonic parametrization. (left)
|
||
mesh with texture, (right) UV parametrization with
|
||
texture</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="leastsquareconformalmaps"><a href="#leastsquareconformalmaps">Least squares conformal maps</a></h2>
|
||
|
||
<p>Least squares conformal maps parametrization <a class="citation" href="#fn:22" title="Jump to citation">[22]<span class="citekey" style="display:none">levy_2002</span></a> minimizes the
|
||
conformal (angular) distortion of the parametrization. Differently from
|
||
harmonic parametrization, it does not need to have a fixed boundary.</p>
|
||
|
||
<p>LSCM minimizes the following energy:</p>
|
||
|
||
<p><span class="math">\[ E_{LSCM}(\mathbf{u},\mathbf{v}) = \int_X \frac{1}{2}| \nabla \mathbf{u}^{\perp} - \nabla \mathbf{v} |^2 dA \]</span></p>
|
||
|
||
<p>which can be rewritten in matrix form as <a class="citation" href="#fn:23" title="Jump to citation">[23]<span class="citekey" style="display:none">mullen_2008</span></a>:</p>
|
||
|
||
<p><span class="math">\[ E_{LSCM}(\mathbf{u},\mathbf{v}) = \frac{1}{2} [\mathbf{u},\mathbf{v}]^t (L_c - 2A) [\mathbf{u},\mathbf{v}] \]</span></p>
|
||
|
||
<p>where <span class="math">\(L_c\)</span> is the cotangent Laplacian matrix and <span class="math">\(A\)</span> is a matrix such that
|
||
<span class="math">\([\mathbf{u},\mathbf{v}]^t A [\mathbf{u},\mathbf{v}]\)</span> is equal to the <a href="http://en.wikipedia.org/wiki/Vector_area">vector
|
||
area</a> of the mesh.</p>
|
||
|
||
<p>Using libigl, this matrix energy can be written in a few lines of code. The
|
||
cotangent matrix can be computed using <code>igl::cotmatrix</code>:</p>
|
||
|
||
<pre><code class="cpp">SparseMatrix<double> L;
|
||
igl::cotmatrix(V,F,L);
|
||
</code></pre>
|
||
|
||
<p>Note that we want to apply the Laplacian matrix to the u and v coordinates at
|
||
the same time, thus we need to extend it taking the left
|
||
Kronecker product with a 2x2 identity matrix:</p>
|
||
|
||
<pre><code class="cpp">SparseMatrix<double> L_flat;
|
||
igl::repdiag(L,2,L_flat);
|
||
</code></pre>
|
||
|
||
<p>The area matrix is computed with <code>igl::vector_area_matrix</code>:</p>
|
||
|
||
<pre><code class="cpp">SparseMatrix<double> A;
|
||
igl::vector_area_matrix(F,A);
|
||
</code></pre>
|
||
|
||
<p>The final energy matrix is <span class="math">\(L_{flat} - 2A\)</span>. Note that in this
|
||
case we do not need to fix the boundary. To remove the null space of the energy and make the minimum unique, it is sufficient to fix two arbitrary
|
||
vertices to two arbitrary positions. The full source code is provided in <a href="502_LSCMParam/main.cpp">Example 502</a>.</p>
|
||
|
||
<figure>
|
||
<img src="images/502_LSCMParam.png" alt="(Example 502) LSCM parametrization. (left) mesh
|
||
with texture, (right) UV parametrization" />
|
||
<figcaption>(<a href="502_LSCMParam/main.cpp">Example 502</a>) LSCM parametrization. (left) mesh
|
||
with texture, (right) UV parametrization</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="asrigidaspossible"><a href="#asrigidaspossible">As-rigid-as-possible parametrization</a></h2>
|
||
|
||
<p>As-rigid-as-possible parametrization <a class="citation" href="#fn:24" title="Jump to citation">[24]<span class="citekey" style="display:none">liu_2008</span></a> is a powerful single-patch,
|
||
non-linear algorithm to compute a parametrization that strives to preserve
|
||
distances (and thus angles). The idea is very similar to ARAP surface
|
||
deformation: each triangle is mapped to the plane trying to preserve its
|
||
original shape, up to a rigid rotation.</p>
|
||
|
||
<p>The algorithm can be implemented reusing the functions discussed in the
|
||
deformation chapter: <code>igl::arap_precomputation</code> and <code>igl::arap_solve</code>. The only
|
||
difference is that the optimization has to be done in 2D instead of 3D and that
|
||
we need to compute a starting point. While for 3D deformation the optimization
|
||
is bootstrapped with the original mesh, this is not the case for ARAP
|
||
parametrization since the starting point must be a 2D mesh. In <a href="503_ARAPParam/main.cpp">Example
|
||
503</a>, we initialize the optimization with harmonic
|
||
parametrization. Similarly to LSCM, the boundary is free to deform to minimize
|
||
the distortion.</p>
|
||
|
||
<figure>
|
||
<img src="images/503_ARAPParam.png" alt="(Example 503) As-Rigid-As-Possible parametrization.
|
||
(left) mesh with texture, (right) UV parametrization with
|
||
texture" />
|
||
<figcaption>(<a href="502_ARAPParam/main.cpp">Example 503</a>) As-Rigid-As-Possible parametrization.
|
||
(left) mesh with texture, (right) UV parametrization with
|
||
texture</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="nrotationallysymmetrictangetfields"><a href="#nrotationallysymmetrictangetfields">N-rotationally symmetric tangent fields</a></h2>
|
||
|
||
<p>The design of tangent fields is a basic tool used to design guidance fields for
|
||
uniform quadrilateral and hexahedral remeshing. Libigl contains an
|
||
implementation of all the state-of-the-art algorithms to design N-RoSy fields
|
||
and their generalizations.</p>
|
||
|
||
<p>In libigl, tangent unit-length vector fields are piece-wise constant on the
|
||
faces of a triangle mesh, and they are described by one or more vectors per-face. The function</p>
|
||
|
||
<pre><code class="cpp">igl::nrosy(V,F,b,bc,b_soft,b_soft_weight,bc_soft,N,0.5,
|
||
output_field,output_singularities);
|
||
</code></pre>
|
||
|
||
<p>creates a smooth unit-length vector field (N=1) starting from a sparse set of
|
||
constrained faces, whose indices are listed in b and their constrained value is
|
||
specified in bc. The functions supports soft_constraints (b_soft,
|
||
b_soft_weight, bc_soft), and returns the interpolated field for each face of
|
||
the triangle mesh (output_field), plus the singularities of the field
|
||
(output_singularities).</p>
|
||
|
||
<figure>
|
||
<img src="images/504_vector_field.png" alt="Design of a unit-length vector field" />
|
||
<figcaption>Design of a unit-length vector field</figcaption>
|
||
</figure>
|
||
|
||
<p>The singularities are vertices where the field vanishes (highlighted in red in
|
||
the figure above). <code>igl::nrosy</code> can also generate N-RoSy fields <a class="citation" href="#fn:25" title="Jump to citation">[25]<span class="citekey" style="display:none">levy_2008</span></a>,
|
||
which are a generalization of vector fields where in every face the vector is
|
||
defined up to a constant rotation of <span class="math">\(2\pi / N\)</span>. As can be observed in
|
||
the following figure, the singularities of the fields generated with different
|
||
N are of different types and they appear in different positions.</p>
|
||
|
||
<figure>
|
||
<img src="images/504_nrosy_field.png" alt="Design of a 2-,4- and 9-RoSy field" />
|
||
<figcaption>Design of a 2-,4- and 9-RoSy field</figcaption>
|
||
</figure>
|
||
|
||
<p>We demonstrate how to call and plot N-RoSy fields in <a href="504_NRosyDesign/main.cpp">Example
|
||
504</a>, where the degree of the field can be change
|
||
pressing the number keys. <code>igl::nrosy</code> implements the algorithm proposed in
|
||
<a class="citation" href="#fn:26" title="Jump to citation">[26]<span class="citekey" style="display:none">bommes_2009</span></a>. N-RoSy fields can also be interpolated with the algorithm
|
||
proposed in <a class="citation" href="#fn:27" title="Jump to citation">[27]<span class="citekey" style="display:none">knoppel_2013</span></a>, see Section <a href="#npolyvectorfields">npolyvectorfields</a> for more details
|
||
(<a href="../include/igl/n_polyvector.h">igl::n_polyvector</a>).</p>
|
||
|
||
<h3 id="globalseamlessintegergridparametrization"><a href="#globalseamlessintegergridparametrization">Global, seamless integer-grid parametrization</a></h3>
|
||
|
||
<p>The previous parametrization methods were focusing on creating parametrizations
|
||
of surface patches aimed at texture mapping or baking of other surface
|
||
properties such as normals and high-frequency details. Global, seamless
|
||
parametrization aims at parametrizing complex shapes with a parametrization
|
||
that is aligned with a given set of directions for the purpose of surface
|
||
remeshing. In libigl, we provide a reference implementation of the pipeline
|
||
proposed in the mixed integer quadrangulation paper <a class="citation" href="#fn:26" title="Jump to citation">[26]<span class="citekey" style="display:none">bommes_2009</span></a>.</p>
|
||
|
||
<p>The first step involves the design of a 4-RoSy field (sometimes called <em>cross</em>
|
||
field) that describes the alignment of the edges of the desired quadrilateral
|
||
remeshing. The field constraints are usually manually specified or extracted
|
||
from the principal curvature directions. In [<a href="506_FrameField/main.cpp">Example
|
||
506</a>], we simply fix one face in a random direction.</p>
|
||
|
||
<figure>
|
||
<img src="images/505_MIQ_1.png" alt="Initial cross field prescribing the edge alignment." />
|
||
<figcaption>Initial cross field prescribing the edge alignment.</figcaption>
|
||
</figure>
|
||
|
||
<h3 id="combingandcutting">Combing and cutting</h3>
|
||
|
||
<p>Given the cross field, we now want to cut the surface so that it becomes
|
||
homeomorphic to a disk. While this could be done directly on the cross-field, we
|
||
opt to perform this operation on its bisector field (a copy of the field
|
||
rotated by 45 degrees) since it is more stable and generic. Working on the
|
||
bisectors allow us to take as input generalized, non-orthogonal and non-unit
|
||
length cross fields.</p>
|
||
|
||
<p>We thus rotate the field,</p>
|
||
|
||
<figure>
|
||
<img src="images/505_MIQ_2.png" alt="Bisector field." />
|
||
<figcaption>Bisector field.</figcaption>
|
||
</figure>
|
||
|
||
<p>and we remove the rotation ambiguity by assigning to each face a u and a v
|
||
direction. The assignment is done with a breadth-first search starting from a
|
||
random face.</p>
|
||
|
||
<figure>
|
||
<img src="images/505_MIQ_3.png" alt="Combed bisector field." />
|
||
<figcaption>Combed bisector field.</figcaption>
|
||
</figure>
|
||
|
||
<p>You can imagine this process as combing an hairy surface: you will be able to
|
||
comb part of it, but at some point you will not be able to consistently comb
|
||
the entire surface (<a href="http://en.wikipedia.org/wiki/Hairy_ball_theorem">Hairy ball
|
||
theorem</a>). The discontinuities
|
||
in the combing define the cut graph:</p>
|
||
|
||
<figure>
|
||
<img src="images/505_MIQ_4.png" alt="Cut graph." />
|
||
<figcaption>Cut graph.</figcaption>
|
||
</figure>
|
||
|
||
<p>Finally, we rotate the combed field by 45 degrees to undo the initial degrees
|
||
rotation:</p>
|
||
|
||
<figure>
|
||
<img src="images/505_MIQ_5.png" alt="Combed cross field." />
|
||
<figcaption>Combed cross field.</figcaption>
|
||
</figure>
|
||
|
||
<p>The combed cross field can be seen as the ideal Jacobian of the parametrization
|
||
that will be computed in the next section.</p>
|
||
|
||
<h3 id="poissonparametrization">Poisson parametrization</h3>
|
||
|
||
<p>The mesh is cut along the seams and a parametrization is computed trying to
|
||
find two scalar functions whose gradient matches the combed cross field
|
||
directions. This is a classical Poisson problem, that is solved minimizing the
|
||
following quadratic energy:</p>
|
||
|
||
<p><span class="math">\[ E(\mathbf{u},\mathbf{v}) = |\nabla \mathbf{u} - X_u|^2 + |\nabla \mathbf{v} - X_v|^2 \]</span></p>
|
||
|
||
<p>where <span class="math">\(X_u\)</span> and <span class="math">\(X_u\)</span> denotes the combed cross field. Solving this
|
||
problem generates a parametrization whose u and v isolines are aligned with the
|
||
input cross field.</p>
|
||
|
||
<figure>
|
||
<img src="images/505_MIQ_8.png" alt="Poisson parametrization." />
|
||
<figcaption>Poisson parametrization.</figcaption>
|
||
</figure>
|
||
|
||
<p>We hide the seams by adding integer constraints to the Poisson problem
|
||
that align the isolines on both sides of each seam <a class="citation" href="#fn:26" title="Jump to citation">[26]<span class="citekey" style="display:none">bommes_2009</span></a>.</p>
|
||
|
||
<figure>
|
||
<img src="images/505_MIQ_7.png" alt="Seamless Poisson parametrization." />
|
||
<figcaption>Seamless Poisson parametrization.</figcaption>
|
||
</figure>
|
||
|
||
<p>Note that this parametrization can only be used for remeshing purposes, since
|
||
it contains many overlaps.</p>
|
||
|
||
<figure>
|
||
<img src="images/505_MIQ_6.png" alt="Seamless Poisson parametrization (in 2D)." />
|
||
<figcaption>Seamless Poisson parametrization (in 2D).</figcaption>
|
||
</figure>
|
||
|
||
<p>A quad mesh can be extracted from this parametrization using
|
||
<a href="https://github.com/hcebke/libQEx">libQEx</a> (not included in libigl).
|
||
The full pipeline is implemented in <a href="505_MIQ/main.cpp">Example 505</a>.</p>
|
||
|
||
<h2 id="anisotropicremeshingusingframefields"><a href="#anisotropicremeshingusingframefields">Anisotropic remeshing</a></h2>
|
||
|
||
<p>Anisotropic and non-uniform quad remeshing is important to concentrate the
|
||
elements in the regions with more details. It is possible to extend the MIQ
|
||
quad meshing framework to generate anisotropic quad meshes using a mesh
|
||
deformation approach <a class="citation" href="#fn:28" title="Jump to citation">[28]<span class="citekey" style="display:none">panozzo_2014</span></a>.</p>
|
||
|
||
<p>The input of the anisotropic remeshing algorithm is a sparse set of constraints
|
||
that define the shape and scale of the desired quads. This can be encoded as a
|
||
frame field, which is a pair of non-orthogonal and non-unit length vectors. The
|
||
frame field can be interpolated by decomposing it in a 4-RoSy field and a
|
||
unique affine transformation. The two parts can then be interpolated
|
||
separately, using <code>igl::nrosy</code> for the cross field, and an harmonic interpolant
|
||
for the affine part.</p>
|
||
|
||
<figure>
|
||
<img src="images/506_FrameField_1.png" alt="Interpolation of a frame field. Colors on the vectors denote the desired
|
||
scale. The red faces contains the frame field
|
||
constraints." />
|
||
<figcaption>Interpolation of a frame field. Colors on the vectors denote the desired
|
||
scale. The red faces contains the frame field
|
||
constraints.</figcaption>
|
||
</figure>
|
||
|
||
<p>After the interpolation, the surface is warped to transform each frame into an
|
||
orthogonal and unit length cross (i.e. removing the scaling and skewness from
|
||
the frame). This deformation defines a new embedding (and a new metric) for the
|
||
surface.</p>
|
||
|
||
<figure>
|
||
<img src="images/506_FrameField_2.png" alt="The surface is deformed to transform the frame field in a cross
|
||
field." />
|
||
<figcaption>The surface is deformed to transform the frame field in a cross
|
||
field.</figcaption>
|
||
</figure>
|
||
|
||
<p>The deformed surface can the be isotropically remeshed using the MIQ algorithm
|
||
that has been presented in the previous section.</p>
|
||
|
||
<figure>
|
||
<img src="images/506_FrameField_3.png" alt="The deformed surface is isotropically remeshed." />
|
||
<figcaption>The deformed surface is isotropically remeshed.</figcaption>
|
||
</figure>
|
||
|
||
<p>The UV coordinates of the deformed surface can then be used to transport the
|
||
parametrization to the original surface, where the isolines will trace a quad
|
||
mesh whose elements are similar to the shape prescribed in the input frame
|
||
field.</p>
|
||
|
||
<figure>
|
||
<img src="images/506_FrameField_4.png" alt="The global parametrization is lifted to the original surface to create the
|
||
anisotropic quad meshing." />
|
||
<figcaption>The global parametrization is lifted to the original surface to create the
|
||
anisotropic quad meshing.</figcaption>
|
||
</figure>
|
||
|
||
<p>Our implementation (<a href="506_FrameField/main.cpp">Example 506</a>) uses MIQ to
|
||
generate the UV parametrization, but other algorithms could be applied: the
|
||
only desiderata is that the generated quad mesh should be as isotropic as
|
||
possible.</p>
|
||
|
||
<h2 id="npolyvectorfields"><a href="#npolyvectorfields">N-PolyVector fields</a></h2>
|
||
|
||
<p>N-RoSy vector fields can be further generalized to represent arbitrary
|
||
vector-sets, with arbitrary angles between them and with arbitrary lengths
|
||
<a class="citation" href="#fn:29" title="Jump to citation">[29]<span class="citekey" style="display:none">diamanti_2014</span></a>. This generalization is called N-PolyVector field, and
|
||
libigl provides the function <code>igl::n_polyvector</code> to design them starting from a
|
||
sparse set of constraints (<a href="507_PolyVectorField/main.cpp">Example 507</a>).</p>
|
||
|
||
<figure>
|
||
<img src="images/507_PolyVectorField.png" alt="Interpolation of a 6-PolyVector field (right) and a 12-PolyVector field from a sparse set of random constraints." />
|
||
<figcaption>Interpolation of a 6-PolyVector field (right) and a 12-PolyVector field from a sparse set of random constraints.</figcaption>
|
||
</figure>
|
||
|
||
<p>The core idea is to represent the vector set as the roots of a complex
|
||
polynomial: The polynomial coefficients are then harmonically interpolated
|
||
leading to polynomials whose roots smoothly vary over the surface.</p>
|
||
|
||
<p>Globally optimal direction fields <a class="citation" href="#fn:27" title="Jump to citation">[27]<span class="citekey" style="display:none">knoppel_2013</span></a> are a special case of
|
||
PolyVector fields. If the constraints are taken from an N-RoSy field,
|
||
<code>igl::n_polyvector</code> generates a field that is equivalent, after normalization,
|
||
to a globally optimal direction field.</p>
|
||
|
||
<h2 id="conjugatevectorfields"><a href="#conjugatevectorfields">Conjugate vector fields</a></h2>
|
||
|
||
<p>Two tangent vectors lying on a face of a triangle mesh are conjugate if</p>
|
||
|
||
<p><span class="math">\[ k_1 (u^T d_1)(v^T d_1) + k_2(u^T d_2)(v^T d_2) = 0. \]</span></p>
|
||
|
||
<p>This condition is very important in architectural geometry: The faces of an
|
||
infinitely dense quad mesh whose edges are aligned with a conjugate field are
|
||
planar. Thus, a quad mesh whose edges follow a conjugate field are easier to
|
||
planarize <a class="citation" href="#fn:30" title="Jump to citation">[30]<span class="citekey" style="display:none">liu_2011</span></a>.</p>
|
||
|
||
<p>Finding a conjugate vector field that satisfies given directional constraints
|
||
is a standard problem in architectural geometry, which can be tackled by
|
||
deforming a Poly-Vector field to the closest conjugate field.</p>
|
||
|
||
<p>This algorithm <a class="citation" href="#fn:29" title="Jump to citation">[29]<span class="citekey" style="display:none">diamanti_2014</span></a> alternates a global step, which enforces
|
||
smoothness, with a local step, that projects the field on every face to the
|
||
closest conjugate field (<a href="508_ConjugateField/main.cpp">Example 508</a>).</p>
|
||
|
||
<figure>
|
||
<img src="images/508_ConjugateField.png" alt="A smooth 4-PolyVector field (left) is deformed to become a conjugate field
|
||
(right)." />
|
||
<figcaption>A smooth 4-PolyVector field (left) is deformed to become a conjugate field
|
||
(right).</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="planarization"><a href="#planarization">Planarization</a></h2>
|
||
|
||
<p>A quad mesh can be transformed in a planar quad mesh with Shape-Up
|
||
<a class="citation" href="#fn:31" title="Jump to citation">[31]<span class="citekey" style="display:none">bouaziz_2012</span></a>, a local/global approach that uses the global step to enforce
|
||
surface continuity and the local step to enforce planarity.</p>
|
||
|
||
<p><a href="509_Planarization/main.cpp">Example 509</a> planarizes a quad mesh until it
|
||
satisfies a user-given planarity threshold.</p>
|
||
|
||
<figure>
|
||
<img src="images/509_Planarization.png" alt="A non-planar quad mesh (left) is planarized using the libigl function
|
||
igl::palanarize (right). The colors represent the planarity of the
|
||
quads." />
|
||
<figcaption>A non-planar quad mesh (left) is planarized using the libigl function
|
||
igl::palanarize (right). The colors represent the planarity of the
|
||
quads.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="integrable"><a href="#integrable">Integrable PolyVector Fields</a></h2>
|
||
|
||
<p>Vector-field guided surface parameterization is based on the idea of designing the gradients
|
||
of the parameterization functions (which are tangent vector fields on the surface) instead of the functions themselves. Thus, vector-set fields (N-Rosy, frame fields, and polyvector fields) that are to be used for parameterization (and subsequent remeshing) need to be integrable: it must be possible to break them down into individual vector fields that are gradients of scalar functions. Fields obtained by most smoothness-based design methods (eg. <a class="citation" href="#fn:25" title="Jump to citation">[25]<span class="citekey" style="display:none">levy_2008</span></a>, <a class="citation" href="#fn:27" title="Jump to citation">[27]<span class="citekey" style="display:none">knoppel_2013</span></a>, <a class="citation" href="#fn:29" title="Jump to citation">[29]<span class="citekey" style="display:none">diamanti_2014</span></a>, <a class="citation" href="#fn:26" title="Jump to citation">[26]<span class="citekey" style="display:none">bommes_2009</span></a>, <a class="citation" href="#fn:28" title="Jump to citation">[28]<span class="citekey" style="display:none">panozzo_2014</span></a>) do not have this property. In <a class="citation" href="#fn:32" title="Jump to citation">[32]<span class="citekey" style="display:none">diamanti_2015</span></a>, a method for creating integrable polyvector fields was introduced. This method takes as input a given field and improves its integrability by removing the vector field curl, thus turning it into a gradient of a function (<a href="510_Integrable/main.cpp">Example 510</a>).</p>
|
||
|
||
<figure>
|
||
<img src="images/510_Integrable.png" alt="Integration error is removed from a frame field to produce a field aligned parameterization free of triangle flips." />
|
||
<figcaption>Integration error is removed from a frame field to produce a field aligned parameterization free of triangle flips.</figcaption>
|
||
</figure>
|
||
|
||
<p>This method retains much of the core principles of the polyvector framework - it expresses the condition for zero discrete curl condition (which typically requires integers for the vector matchings) into a condition involving continuous variables only. This is done using coefficients of appropriately defined polynomials. The parameterizations generated by the resulting fields are exactly aligned to the field directions and contain no inverted triangles.</p>
|
||
|
||
<h2 id="npolyvectorfields_general"><a href="#npolyvectorfields_general">General N-PolyVector fields</a></h2>
|
||
|
||
<p>While mostly applicable for the design of symmetric fields (i.e. fields that comprise of vector sets with symmetries between them at each point, e.g. N-RoSy or frame-fields), the framework presented in <a class="citation" href="#fn:29" title="Jump to citation">[29]<span class="citekey" style="display:none">diamanti_2014</span></a> can be used to design completely general fields, with possibly no such symmetries. For example, one can design fields that at each point comprise of an arbitrary number of vectors, not required to be collinear - as opposed e.g. to the case of the 4 pairwise-collinear vectors designed in the example (<a href="507_PolyVectorField/main.cpp">Example 507</a>). This capability is implemented in the function igl::n_polyvector_general, and is illustrated in the example (<a href="511_PolyVectorFieldGeneral/main.cpp">Example 511</a>).</p>
|
||
|
||
<figure>
|
||
<img src="images/511_PolyVectorFieldGeneral.png" alt="Interpolation of a general field with 3 (left) and 9 vectors per point field from a sparse set of random constraints (in red). The field is defined on all mesh faces, but is only shown on a subset for clarity. " />
|
||
<figcaption>Interpolation of a general field with 3 (left) and 9 vectors per point field from a sparse set of random constraints (in red). The field is defined on all mesh faces, but is only shown on a subset for clarity. </figcaption>
|
||
</figure>
|
||
|
||
<p>The design of these general directional fields (also called vector-set fields) is based on the same polynomial framework and includes the symmetric fields as a special case. Note that in the case that some symmetries do exist in the constraints, the final field is not guaranteed to have these symmetries everywhere else on the mesh. For example, designing a field with 3 vectors per point where, at the constrained faces, two of the vectors are on a line opposite to each other, we are not guaranteed to always have two pairwise-collinear vectors everywhere in the result, as can be seen in the picture. In some cases however (as is the case of the frame field in the previous example <a href="507_PolyVectorField/main.cpp">Example 507</a>) these symmetries are in fact guaranteed due to the particular nature of the polynomial that applies in that case (two coefficients are 0).</p>
|
||
|
||
<p>For a complete categorization of fields used in various applications (including these general ones) see Vaxman et al. 2016 <a class="citation" href="#fn:33" title="Jump to citation">[33]<span class="citekey" style="display:none">vaxman_2016</span></a>.</p>
|
||
|
||
<h1 id="chapter6:externallibraries">Chapter 6: External libraries</h1>
|
||
|
||
<p>An additional positive side effect of using matrices as basic types is that it
|
||
is easy to exchange data between libigl and other softwares and libraries.</p>
|
||
|
||
<h2 id="stateserialization"><a href="#stateserialization">State serialization</a></h2>
|
||
|
||
<p>Geometry processing applications often require a considerable amount of
|
||
computational time and/or manual input. Serializing the state of the application is a simple strategy to greatly increase the development efficiency. It allows to quickly start debugging just
|
||
before the crash happens, avoiding to wait for the precomputation to take place
|
||
every time and it also makes your experiments reproducible, allowing to quickly test algorithms variants on the same input data.</p>
|
||
|
||
<p>Serialization is often not considered in geometry processing due
|
||
to the extreme difficulty in serializing pointer-based data structured, such as
|
||
an half-edge data structure (<a href="http://openmesh.org">OpenMesh</a>, <a href="http://www.cgal.org">CGAL</a>), or a pointer based indexed structure (<a href="http://vcg.isti.cnr.it/~cignoni/newvcglib/html/">VCG</a>).</p>
|
||
|
||
<p>In libigl, serialization is much simpler, since the majority of the functions use basic types, and pointers are used in very rare cases (usually to interface
|
||
with external libraries). Libigl bundles a simple and self-contained binary and XML serialization framework, that drastically reduces the overhead required to add
|
||
serialization to your applications.</p>
|
||
|
||
<p>To de-/serialize a set of variables use the following method:</p>
|
||
|
||
<pre><code class="cpp">#include "igl/serialize.h"
|
||
|
||
bool b = true;
|
||
unsigned int num = 10;
|
||
std::vector<float> vec = {0.1,0.002,5.3};
|
||
|
||
// use overwrite = true for the first serialization to create or overwrite an existing file
|
||
igl::serialize(b,"B","filename",true);
|
||
// append following serialization to existing file
|
||
igl::serialize(num,"Number","filename");
|
||
igl::serialize(vec,"VectorName","filename");
|
||
|
||
// deserialize back to variables
|
||
igl::deserialize(b,"B","filename");
|
||
igl::deserialize(num,"Number","filename");
|
||
igl::deserialize(vec,"VectorName","filename");
|
||
</code></pre>
|
||
|
||
<p>Currently all fundamental data types (bool, int, float, double, …) are supported, as well as std::string, basic <code>STL</code> containers, dense and sparse Eigen matrices and nestings of those.
|
||
Some limitations apply to pointers. Currently, loops or many to one type of link structures are not handled correctly. Each pointer is assumed to point to a different independent object.
|
||
Uninitialized pointers must be set to <code>nullptr</code> before de-/serialization to avoid memory leaks. Cross-platform issues like little-, big-endianess is currently not supported.
|
||
To make user defined types serializable, just derive from <code>igl::Serializable</code> and trivially implementing the <code>InitSerialization</code> method.</p>
|
||
|
||
<p>Assume that the state of your application is a mesh and a set of integer ids:</p>
|
||
|
||
<pre><code class="cpp">#include "igl/serialize.h"
|
||
|
||
struct State : public igl::Serializable
|
||
{
|
||
Eigen::MatrixXd V;
|
||
Eigen::MatrixXi F;
|
||
std::vector<int> ids;
|
||
|
||
void InitSerialization()
|
||
{
|
||
this->Add(V , "V");
|
||
this->Add(F , "F");
|
||
this->Add(ids, "ids");
|
||
}
|
||
};
|
||
</code></pre>
|
||
|
||
<p>If you need more control over the serialization of your types, you can override the following functions or directly inherit from the interface <code>igl::SerializableBase</code>.</p>
|
||
|
||
<pre><code class="cpp">bool Serializable::PreSerialization() const;
|
||
void Serializable::PostSerialization() const;
|
||
bool Serializable::PreDeserialization();
|
||
void Serializable::PostDeserialization();
|
||
</code></pre>
|
||
|
||
<p>Alternatively, if you want a non-intrusive way of serializing your state you can overload the following functions:</p>
|
||
|
||
<pre><code class="cpp">namespace igl
|
||
{
|
||
namespace serialization
|
||
{
|
||
template <> inline void serialize(const State& obj,std::vector<char>& buffer){
|
||
::igl::serialize(obj.V,std::string("V"),buffer);
|
||
::igl::serialize(obj.F,std::string("F"),buffer);
|
||
::igl::serialize(obj.ids,std::string("ids"),buffer);
|
||
}
|
||
template <> inline void deserialize(State& obj,const std::vector<char>& buffer){
|
||
::igl::deserialize(obj.V,std::string("V"),buffer);
|
||
::igl::deserialize(obj.F,std::string("F"),buffer);
|
||
::igl::deserialize(obj.ids,std::string("ids"),buffer);
|
||
}
|
||
}
|
||
}
|
||
</code></pre>
|
||
|
||
<p>Equivalently, you can use the following macros:</p>
|
||
|
||
<pre><code class="cpp">SERIALIZE_TYPE(State,
|
||
SERIALIZE_MEMBER(V)
|
||
SERIALIZE_MEMBER(F)
|
||
SERIALIZE_MEMBER_NAME(ids,"ids")
|
||
)
|
||
</code></pre>
|
||
|
||
<p>All the former code is for binary serialization which is especially useful if you have to handle larger data where the loading and saving times become more important.
|
||
For cases where you want to read and edit the serialized data by hand we provide a serialization to XML files which is based on the library <a href="https://github.com/leethomason/tinyxml2">tinyxml2</a>.
|
||
There you also have the option to create a partial binary serialization of your data by using the binary parameter, exposed in the function <code>serialize_xml()</code>:</p>
|
||
|
||
<pre><code class="cpp">#include "igl/xml/serialize_xml.h"
|
||
|
||
int number;
|
||
|
||
// binary = false, overwrite = true
|
||
igl::serialize_xml(vec,"VectorXML",xmlFile,false,true);
|
||
// binary = true, overwrite = true
|
||
igl::serialize_xml(vec,"VectorBin",xmlFile,true,true);
|
||
igl::deserialize_xml(vec,"VectorXML",xmlFile);
|
||
igl::deserialize_xml(vec,"VectorBin",xmlFile);
|
||
</code></pre>
|
||
|
||
<p>For user defined types derive from <code>XMLSerializable</code>.</p>
|
||
|
||
<p>The code snippets above are extracted from <a href="601_Serialization/main.cpp">Example
|
||
601</a>. We strongly suggest that you make the entire
|
||
state of your application always serializable since it will save you a lot of
|
||
troubles when you will be preparing figures for a scientific report. It is very
|
||
common to have to do small changes to figures, and being able to serialize the
|
||
entire state just before you take screenshots will save you many painful hours
|
||
before a submission deadline.</p>
|
||
|
||
<h2 id="mixingmatlabcode"><a href="#mixingmatlabcode">Mixing Matlab code</a></h2>
|
||
|
||
<p>Libigl can be interfaced with Matlab to offload numerically heavy computation
|
||
to a Matlab script. The major advantage of this approach is that you will be
|
||
able to develop efficient and complex user-interfaces in C++, while exploting
|
||
the syntax and fast protototyping features of matlab. In particular, the use of
|
||
an external Matlab script in a libigl application allows to change the Matlab
|
||
code while the C++ application is running, greatly increasing coding
|
||
efficiency.</p>
|
||
|
||
<p>We demonstrate how to integrate Matlab in a libigl application in <a href="602_Matlab/main.cpp">Example
|
||
602</a>. The example uses Matlab to compute the
|
||
Eigenfunctions of the discrete Laplacian operator, relying on libigl for mesh
|
||
IO, visualization and for computing the Laplacian operator.</p>
|
||
|
||
<p>Libigl can connect to an existing instance of Matlab (or launching a new one on
|
||
Linux/MacOSX) using:</p>
|
||
|
||
<pre><code class="cpp">igl::mlinit(&engine);
|
||
</code></pre>
|
||
|
||
<p>The cotangent Laplacian is computed using igl::cotmatrix and uploaded to the
|
||
Matlab workspace:</p>
|
||
|
||
<pre><code class="cpp">igl::cotmatrix(V,F,L);
|
||
igl::mlsetmatrix(&engine,"L",L);
|
||
</code></pre>
|
||
|
||
<p>It is now possible to use any Matlab function on the data. For example, we can
|
||
see the sparsity pattern of L using spy:</p>
|
||
|
||
<pre><code class="cpp">igl::mleval(&engine,"spy(L)");
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/602_Matlab_1.png" alt="The Matlab spy function is called from a libigl-based
|
||
application." />
|
||
<figcaption>The Matlab spy function is called from a libigl-based
|
||
application.</figcaption>
|
||
</figure>
|
||
|
||
<p>The results of Matlab computations can be returned back to the C++ application</p>
|
||
|
||
<pre><code class="cpp">igl::mleval(&engine,"[EV,~] = eigs(-L,10,'sm')");
|
||
igl::mlgetmatrix(&engine,"EV",EV);
|
||
</code></pre>
|
||
|
||
<p>and plotted using the libigl viewer.</p>
|
||
|
||
<figure>
|
||
<img src="images/602_Matlab_2.png" alt="4 Eigenfunctions of the Laplacian plotted in the libigl
|
||
viewer." />
|
||
<figcaption>4 Eigenfunctions of the Laplacian plotted in the libigl
|
||
viewer.</figcaption>
|
||
</figure>
|
||
|
||
<h3 id="savingamatlabworkspace">Saving a Matlab workspace</h3>
|
||
|
||
<p>To aid debugging, libigl also supplies functions to write Matlab <code>.mat</code>
|
||
“Workspaces”. This C++ snippet saves a mesh and it’s sparse Laplacian matrix to
|
||
a file:</p>
|
||
|
||
<pre><code class="cpp">igl::readOFF(TUTORIAL_SHARED_PATH "/fertility.off", V, F);
|
||
igl::cotmatrix(V,F,L);
|
||
igl::MatlabWorkspace mw;
|
||
mw.save(V,"V");
|
||
mw.save_index(F,"F");
|
||
mw.save(L,"L");
|
||
mw.write("fertility.mat");
|
||
</code></pre>
|
||
|
||
<p>Then this workspace can be loaded into a Matlab IDE:</p>
|
||
|
||
<pre><code class="matlab">load fertility.mat
|
||
</code></pre>
|
||
|
||
<p>The <code>igl::MatlabWorkspace</code> depends on Matlab libraries to compile and run,
|
||
but—in contrast to the engine routines above—will avoid launching a Matlab
|
||
instance upon execution.</p>
|
||
|
||
<h3 id="dumpingeigenmatricestocopyandpasteintomatlab">Dumping Eigen matrices to copy and paste into Matlab</h3>
|
||
|
||
<p>Eigen supplies a sophisticated API for printing its matrix types to the screen.
|
||
Libigl has wrapped up a particularly useful formatting which makes it simple to
|
||
copy standard output from a C++ program into a Matlab IDE. The code:</p>
|
||
|
||
<pre><code class="cpp">igl::readOFF(TUTORIAL_SHARED_PATH "/2triangles.off", V, F);
|
||
igl::cotmatrix(V,F,L);
|
||
std::cout<<igl::matlab_format(V,"V")<<std::endl;
|
||
std::cout<<igl::matlab_format((F.array()+1).eval(),"F")<<std::endl;
|
||
std::cout<<igl::matlab_format(L,"L")<<std::endl;
|
||
</code></pre>
|
||
|
||
<p>produces the output:</p>
|
||
|
||
<pre><code class="matlab">V = [
|
||
0 0 0
|
||
1 0 0
|
||
1 1 1
|
||
2 1 0
|
||
];
|
||
F = [
|
||
1 2 3
|
||
2 4 3
|
||
];
|
||
LIJV = [
|
||
1 1 -0.7071067811865476
|
||
2 1 0.7071067811865475
|
||
3 1 1.570092458683775e-16
|
||
1 2 0.7071067811865475
|
||
2 2 -1.638010440969447
|
||
3 2 0.6422285251880865
|
||
4 2 0.2886751345948129
|
||
1 3 1.570092458683775e-16
|
||
2 3 0.6422285251880865
|
||
3 3 -0.9309036597828995
|
||
4 3 0.2886751345948129
|
||
2 4 0.2886751345948129
|
||
3 4 0.2886751345948129
|
||
4 4 -0.5773502691896258
|
||
];
|
||
L = sparse(LIJV(:,1),LIJV(:,2),LIJV(:,3));
|
||
</code></pre>
|
||
|
||
<p>which is easily copied and pasted into Matlab for debugging, etc.</p>
|
||
|
||
<h2 id="callinglibiglfunctionsfrommatlab"><a href="#callinglibiglfunctionsfrommatlab">Calling libigl functions from Matlab</a></h2>
|
||
|
||
<p>It is also possible to call libigl functions from matlab, compiling them as MEX
|
||
functions. This can be used to offload to C++ code the computationally
|
||
intensive parts of a Matlab application.</p>
|
||
|
||
<p>We provide a wrapper for <code>igl::readOBJ</code> in <a href="603_MEX/compileMEX.m">Example 603</a>.
|
||
We plan to provide wrappers for all our functions in the future, if you are
|
||
interested in this feature (or if you want to help implementing it) please let
|
||
us know.</p>
|
||
|
||
<h2 id="triangulationofclosedpolygons"><a href="#triangulationofclosedpolygons">Triangulation of closed polygons</a></h2>
|
||
|
||
<p>The generation of high-quality triangle and tetrahedral meshes is a very common
|
||
task in geometry processing. We provide wrappers in libigl to
|
||
<a href="http://www.cs.cmu.edu/~quake/triangle.html">triangle</a> and
|
||
<a href="http://wias-berlin.de/software/tetgen/">Tetgen</a>.</p>
|
||
|
||
<p>A triangle mesh with a given boundary can be created with:</p>
|
||
|
||
<pre><code class="cpp">igl::triangulate(V,E,H,V2,F2,"a0.005q");
|
||
</code></pre>
|
||
|
||
<p>where <code>E</code> is a set of boundary edges (#E by 2), <code>H</code> is a set of 2D positions of
|
||
points contained in holes of the triangulation (#H by 2) and (<code>V2</code>,<code>F2</code>) is the
|
||
generated triangulation. Additional parameters can be passed to <code>triangle</code>, to
|
||
control the quality: <code>"a0.005q"</code> enforces a bound on the maximal area of the
|
||
triangles and a minimal angle of 20 degrees. In <a href="604_Triangle/main.cpp">Example
|
||
604</a>, the interior of a square (excluded a smaller square
|
||
in its interior) is triangulated.</p>
|
||
|
||
<figure>
|
||
<img src="images/604_Triangle.png" alt="Triangulation of the interior of a polygon." />
|
||
<figcaption>Triangulation of the interior of a polygon.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="tetrahedralizationofclosedsurfaces"><a href="#tetrahedralizationofclosedsurfaces">Tetrahedralization of closed surfaces</a></h2>
|
||
|
||
<p>Similarly, the interior of a closed manifold surface can be tetrahedralized
|
||
using the function <code>igl::tetrahedralize</code> which wraps the Tetgen library (<a href="605_Tetgen/main.cpp">Example
|
||
605</a>):</p>
|
||
|
||
<pre><code class="cpp">igl::tetrahedralize(V,F,"pq1.414", TV,TT,TF);
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/605_Tetgen.png" alt="Tetrahedralization of the interior of a surface mesh." />
|
||
<figcaption>Tetrahedralization of the interior of a surface mesh.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="bakingambientocclusion"><a href="#bakingambientocclusion">Baking ambient occlusion</a></h2>
|
||
|
||
<p><a href="http://en.wikipedia.org/wiki/Ambient_occlusion">Ambient occlusion</a> is a
|
||
rendering technique used to calculate the exposure of each point in a surface
|
||
to ambient lighting. It is usually encoded as a scalar (normalized between 0
|
||
and 1) associated with the vertice of a mesh.</p>
|
||
|
||
<p>Formally, ambient occlusion is defined as:</p>
|
||
|
||
<p><span class="math">\[ A_p = \frac{1}{\pi} \int_\omega V_{p,\omega}(n \cdot \omega) d\omega \]</span></p>
|
||
|
||
<p>where <span class="math">\(V_{p,\omega}\)</span> is the visibility function at p, defined to be zero if p
|
||
is occluded in the direction <span class="math">\(\omega\)</span> and one otherwise, and <span class="math">\(d\omega\)</span> is the
|
||
infinitesimal solid angle step of the integration variable <span class="math">\(\omega\)</span>.</p>
|
||
|
||
<p>The integral is usually approximated by casting rays in random directions
|
||
around each vertex. This approximation can be computed using the function:</p>
|
||
|
||
<pre><code class="cpp">igl::ambient_occlusion(V,F,V_samples,N_samples,500,AO);
|
||
</code></pre>
|
||
|
||
<p>that given a scene described in <code>V</code> and <code>F</code>, computes the ambient occlusion of
|
||
the points in <code>V_samples</code> whose associated normals are <code>N_samples</code>. The
|
||
number of casted rays can be controlled (usually at least 300–500 rays are
|
||
required to get a smooth result) and the result is returned in <code>AO</code>, as a
|
||
single scalar for each sample.</p>
|
||
|
||
<p>Ambient occlusion can be used to darken the surface colors, as shown in
|
||
<a href="606_AmbientOcclusion/main.c">Example 606</a></p>
|
||
|
||
<figure>
|
||
<img src="images/606_AmbientOcclusion.png" alt="A mesh rendered without (left) and with (right) ambient
|
||
occlusion." />
|
||
<figcaption>A mesh rendered without (left) and with (right) ambient
|
||
occlusion.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="screencapture"><a href="#screencapture">Screen Capture</a></h2>
|
||
|
||
<p>Libigl supports read and writing to .png files via the
|
||
<a href="http://nothings.org/stb_image.h">stb image</a> code.</p>
|
||
|
||
<p>With the viewer used in this tutorial, it is possible to render the scene in a
|
||
memory buffer using the function, <code>igl::viewer::ViewerCore::draw_buffer</code>:</p>
|
||
|
||
<pre><code class="cpp">// Allocate temporary buffers for 1280x800 image
|
||
Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> R(1280,800);
|
||
Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> G(1280,800);
|
||
Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> B(1280,800);
|
||
Eigen::Matrix<unsigned char,Eigen::Dynamic,Eigen::Dynamic> A(1280,800);
|
||
|
||
// Draw the scene in the buffers
|
||
viewer.core.draw_buffer(viewer.data,viewer.opengl,false,R,G,B,A);
|
||
|
||
// Save it to a PNG
|
||
igl::png::writePNG(R,G,B,A,"out.png");
|
||
</code></pre>
|
||
|
||
<p>In <a href="607_ScreenCapture/main.cpp">Example 607</a> a scene is rendered in a temporary
|
||
png and used to texture a quadrilateral.</p>
|
||
|
||
<h2 id="locallyinjectivemaps"><a href="#locallyinjectivemaps">Locally Injective Maps</a></h2>
|
||
|
||
<p>Extreme deformations or parametrizations with high-distortion might flip
|
||
elements. This is undesirable in many applications, and it is possible to
|
||
avoid it by introducing a non-linear constraints that guarantees that the area
|
||
of every element remain positive.</p>
|
||
|
||
<p>Libigl can be used to compute Locally Injective Maps <a class="citation" href="#fn:34" title="Jump to citation">[34]<span class="citekey" style="display:none">schuller_2013</span></a> using a variety of
|
||
deformation energies. A simple deformation of a 2D grid is computed in <a href="608_LIM/main.cpp">Example
|
||
608</a>.</p>
|
||
|
||
<figure>
|
||
<img src="images/608_LIM.png" alt="A mesh (left) deformed using Laplacian editing (middle) and with Laplacian
|
||
editing plus the anti-flipping constraints (right)." />
|
||
<figcaption>A mesh (left) deformed using Laplacian editing (middle) and with Laplacian
|
||
editing plus the anti-flipping constraints (right).</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="booleanoperationsonmeshes"><a href="#booleanoperationsonmeshes">Boolean operations on meshes</a></h2>
|
||
|
||
<p>Constructive solid geometry (CSG) is a technique to define a complex surface as
|
||
the result of a number of set operations on solid regions of space: union,
|
||
intersection, set difference, symmetric difference, complement. Typically, CSG
|
||
libraries represent the inputs and outputs to these operations <em>implicitly</em>:
|
||
the solid <span class="math">\(A\)</span> is defined as the open set of points <span class="math">\(\mathbf{x}\)</span> for which some
|
||
function <span class="math">\(a(\mathbf{x})\)</span> “returns true”. The surface of this shape is the
|
||
<em>closure</em> of all points <span class="math">\(x\)</span> in <span class="math">\(A\)</span>.</p>
|
||
|
||
<p>With this sort of representation, boolean
|
||
operations are straightforward. For example, the union of solids <span class="math">\(A\)</span> and <span class="math">\(B\)</span>
|
||
is simply</p>
|
||
|
||
<p><span class="math">\(A \cup B = \{\mathbf{x} \left.\right|
|
||
a(\mathbf{x}) \text{ or } b(\mathbf{x})\},\)</span></p>
|
||
|
||
<p>the intersection is</p>
|
||
|
||
<p><span class="math">\(A \cap B = \{\mathbf{x} \left.\right|
|
||
a(\mathbf{x}) \text{ and } b(\mathbf{x})\},\)</span></p>
|
||
|
||
<p>the difference <span class="math">\(A\)</span> <em>minus</em> <span class="math">\(B\)</span> is</p>
|
||
|
||
<p><span class="math">\(A \setminus B = \{\mathbf{x} \left.\right|
|
||
a(\mathbf{x}) \text{ and _not_ } b(\mathbf{x})\},\)</span></p>
|
||
|
||
<p>and the symmetric difference (XOR) is</p>
|
||
|
||
<p><span class="math">\(A \triangle B = \{\mathbf{x} \left.\right|
|
||
\text{either } a(\mathbf{x}) \text{ or } b(\mathbf{x}) \text{ but not both }\}.\)</span></p>
|
||
|
||
<p>Stringing together many of these operations, one can design quite complex
|
||
shapes. A typical CSG library might only keep explicit <em>base-case</em>
|
||
representations of canonical shapes: half-spaces, quadrics, etc.</p>
|
||
|
||
<p>In libigl, we do currently <em>not</em> have an implicit surface representation.
|
||
Instead we expect our users to be working with <em>explicit</em> triangle mesh
|
||
<em>boundary representations</em> of solid shapes. CSG operations are much hard to
|
||
compute robustly with boundary representations, but are nonetheless useful.</p>
|
||
|
||
<p>To compute a boolean operation on a triangle mesh with vertices <code>VA</code> and
|
||
triangles <code>FA</code> and another mesh <code>VB</code> and <code>FB</code>, libigl first computes a unified
|
||
“mesh arrangement” (see <a class="citation" href="#fn:35" title="Jump to citation">[35]<span class="citekey" style="display:none">zhou_2016</span></a>) with vertices <code>V</code> and triangles <code>F</code> where all triangle-triangle
|
||
intersections have been “resolved”. That is, edges and vertices are added
|
||
exactly at the intersection lines, so the resulting <em>non-manifold</em> mesh <code>(V,F)</code>
|
||
has no self-intersections.</p>
|
||
|
||
<p>Then libigl labels each “cell” bounded by surfaces of the arrangement according
|
||
to its <em>winding number vector</em>: winding number with respect to each input mesh
|
||
<span class="math">\((w_A,w_B)\)</span>. Finally, according to the desired operation (e.g. union,
|
||
intersection) the boundary of the corresponding cells are extracted.</p>
|
||
|
||
<p>Calling libigl’s boolean operations is simple. To compute the union of
|
||
<code>(VA,FA)</code> and <code>(VB,FB)</code> into a new mesh <code>(VC,FC)</code>, use:</p>
|
||
|
||
<pre><code class="cpp">igl::copyleft::cgal::mesh_boolean(VA,FA,VB,FB,MESH_BOOLEAN_TYPE_UNION,VC,FC);
|
||
</code></pre>
|
||
|
||
<p>The following figure shows each boolean operation on two meshes.</p>
|
||
|
||
<figure>
|
||
<img src="images/cheburashka-knight-boolean.jpg" alt="The example Boolean conducts
|
||
boolean operations on the Cheburashka (red) and Knight (green). From left
|
||
to right: union, intersection, set minus, symmetric difference (XOR),
|
||
resolve. Bottom row reveals inner surfaces, darker color indicates
|
||
back-facing triangles." />
|
||
<figcaption>The example <a href="609_Boolean/main.cpp">Boolean</a> conducts
|
||
boolean operations on the <em>Cheburashka</em> (red) and <em>Knight</em> (green). From left
|
||
to right: union, intersection, set minus, symmetric difference (XOR),
|
||
“resolve”. Bottom row reveals inner surfaces, darker color indicates
|
||
back-facing triangles.</figcaption>
|
||
</figure>
|
||
|
||
<p>The union, symmetric difference and “resolve” have the same outward
|
||
appearance, but differ in their treatment of internal structures. The union has
|
||
no internal surfaces: the triangles are not included in the output. The
|
||
symmetric difference is the same set of triangles as the “resolve”, but
|
||
internal surfaces have been reversed in orientation, indicating that the solid
|
||
result of the operation. The “resolve” operation is not really a boolean
|
||
operation, it is simply the result of resolving all intersections and gluing
|
||
together coincident vertices, maintaining original triangle orientations.</p>
|
||
|
||
<p>Libigl also provides a wrapper <code>igl::copyleft::cork::mesh_boolean</code> to the
|
||
<a href="https://github.com/gilbo/cork">cork</a>, which is typically faster, but is not
|
||
always robust.</p>
|
||
|
||
<h2 id="csgtree"><a href="#csgtree">CSG Tree</a></h2>
|
||
|
||
<p>The <a href="#booleanoperationsonmeshes">previous section</a> discusses using
|
||
<code>igl::copyleft::cgal::mesh_boolean</code> to compute the result of a <em>single</em> boolean
|
||
operation on two input triangle meshes. When employing constructive solid
|
||
geometry (CSG) as a modeling paradigm, shapes are represented as the result of
|
||
many such binary operations. The sequence is stored in a binary tree.</p>
|
||
|
||
<p>Libigl uses exact arithmetic internally to construct the intermediary boolean
|
||
results robustly. “Rounding” this result to floating point (even double
|
||
precision) would cause problems if re-injected into a further boolean
|
||
operation. To facilitate CSG tree operations and encourage callers <em>not</em> to
|
||
call <code>igl::copyleft::cgal::mesh_boolean</code> multiple times explicitly, libigl implements
|
||
a class <code>igl::copyleft::cgal::CSGTree</code>. Leaf nodes of this class are simply “solid”
|
||
meshes (otherwise good input to <code>igl::copyleft::cgal::mesh_boolean</code>). Interior nodes
|
||
of the tree combine two children with a boolean operation. Using the intializer
|
||
list constructor it is easy to hard-code specific tree constructions. Here’s an
|
||
example taking the <em>intersection</em> of a cube A and sphere B <em>minus</em> the <em>union</em>
|
||
of three cylinders:</p>
|
||
|
||
<pre><code class="cpp">// Compute result of (A ∩ B) \ ((C ∪ D) ∪ E)
|
||
igl::copyleft::cgal::CSGTree<MatrixXi> CSGTree =
|
||
{{{VA,FA},{VB,FB},"i"},{{{VC,FC},{VD,FD},"u"},{VE,FE},"u"},"m"};
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/cube-sphere-cylinders-csg-tree.jpg" alt="A CSG Tree represents a shape as a combination of binary boolean
|
||
operations" />
|
||
<figcaption>A CSG Tree represents a shape as a combination of binary boolean
|
||
operations</figcaption>
|
||
</figure>
|
||
|
||
<p>Example <a href="610_CSGTree/main.cpp">610</a> computes each intermediary CSG result and
|
||
then the final composite.</p>
|
||
|
||
<figure>
|
||
<img src="images/cube-sphere-cylinders-csg.gif" alt="Example 610 computes complex CSG Tree operation on 5
|
||
input meshes." />
|
||
<figcaption>Example <a href="610_CSGTree/main.cpp">610</a> computes complex CSG Tree operation on 5
|
||
input meshes.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="meshstatistics"><a href="#meshstatistics">Mesh Statistics</a></h2>
|
||
|
||
<p>Libigl contains various mesh statistics, including face angles, face areas and
|
||
the detection of singular vertices, which are vertices with more or less than 6
|
||
neighbours in triangulations or 4 in quadrangulations.</p>
|
||
|
||
<p>The example <a href="701_Statistics/main.cpp">Statistics</a> computes these quantities and
|
||
does a basic statistic analysis that allows to estimate the isometry and
|
||
regularity of a mesh:</p>
|
||
|
||
<pre><code class="bash">Irregular vertices:
|
||
136/2400 (5.67%)
|
||
Areas (Min/Max)/Avg_Area Sigma:
|
||
0.01/5.33 (0.87)
|
||
Angles in degrees (Min/Max) Sigma:
|
||
17.21/171.79 (15.36)
|
||
</code></pre>
|
||
|
||
<p>The first row contains the number and percentage of irregular vertices, which
|
||
is particularly important for quadrilateral meshes when they are used to define
|
||
subdivision surfaces: every singular point will result in a point of the
|
||
surface that is only C<sup>1</sup>.</p>
|
||
|
||
<p>The second row reports the area of the minimal element, maximal element and the
|
||
standard deviation. These numbers are normalized by the mean area, so in the
|
||
example above 5.33 max area means that the biggest face is 5 times larger than
|
||
the average face. An ideal isotropic mesh would have both min and max area
|
||
close to 1.</p>
|
||
|
||
<p>The third row measures the face angles, which should be close to 60 degrees (90
|
||
for quads) in a perfectly regular triangulation. For FEM purposes, the closer
|
||
the angles are to 60 degrees the more stable will the optimization be. In this
|
||
case, it is clear that the mesh is of bad quality and it will probably result
|
||
in artifacts if used for solving PDEs.</p>
|
||
|
||
<h2 id="generalizedwindingnumber"><a href="#generalizedwindingnumber">Generalized Winding Number</a></h2>
|
||
|
||
<p>The problem of tetrahedralizing the interior of closed watertight surface mesh
|
||
is a difficult, but well-posed problem (see our <a href="#tetrahedralizationofclosedsurfaces">Tetgen wrappers</a>). But
|
||
black-box tet-meshers like TetGen will <em>refuse</em> input triangle meshes with
|
||
self-intersections, open boundaries, non-manifold edges from multiple connected
|
||
components.
|
||
The problem is two-fold: self-intersections present contradictory facet
|
||
constraints and self-intersections/open-boundaries/non-manifold edges make the
|
||
problem of determining inside from outside ill-posed without further
|
||
assumptions.</p>
|
||
|
||
<p>The first problem is <em>easily</em> solved by “resolving” all self-intersections.
|
||
That is, meshing intersecting triangles so that intersects occur exactly at
|
||
edges and vertices. This is accomplished using <code>igl::selfintersect</code>.</p>
|
||
|
||
<p>TetGen can usually tetrahedralize the convex hull of this “resolved” mesh, and
|
||
then the problem becomes determining which of these tets are <em>inside</em> the input
|
||
mesh and which are outside. That is, which should be kept and which should be
|
||
removed.</p>
|
||
|
||
<p>The “Generalized Winding Number” is a robust method for determined
|
||
inside and outside for troublesome meshes <a class="citation" href="#fn:36" title="Jump to citation">[36]<span class="citekey" style="display:none">jacobson_2013</span></a>. The generalized
|
||
winding number with respect to <code>(V,F)</code> at some point <span class="math">\(\mathbf{p} \in
|
||
\mathcal{R}^3\)</span> is defined as scalar function:</p>
|
||
|
||
<p><span class="math">\[w(\mathbf{p}) = \sum\limits_{f_i\in F} \frac{1}{4\pi}\Omega_{f_i}(\mathbf{p})\]</span></p>
|
||
|
||
<p>where <span class="math">\(\Omega_{f_i}\)</span> is the <em>solid angle</em> subtended by <span class="math">\(f_i\)</span> (the ith face in
|
||
<code>F</code>) at the point <span class="math">\(\mathbf{p}\)</span>. This solid angle contribution is a simple,
|
||
closed-form expression involving <code>atan2</code> and some dot-products.</p>
|
||
|
||
<p>If <code>(V,F)</code> <em>does</em> form a closed watertight surface, then <span class="math">\(w(\mathbf{p})=1\)</span> if
|
||
<span class="math">\(\mathbf{p}\)</span> lies inside <code>(V,F)</code> and <span class="math">\(w(\mathbf{p})=0\)</span> if outside <code>(V,F)</code>. If
|
||
<code>(V,F)</code> is closed but overlaps itself then <span class="math">\(w(\mathbf{p})\)</span> is an integer value
|
||
counting how many (signed) times <code>(V,F)</code> <em>wraps</em> around <span class="math">\(\mathbf{p}\)</span>. Finally,
|
||
if <code>(V,F)</code> is not closed or not even manifold (but at least consistently
|
||
oriented), then <span class="math">\(w(\mathbf{p})\)</span> tends smoothly toward 1 as <span class="math">\(\mathbf{p}\)</span> is
|
||
<em>more</em> inside <code>(V,F)</code>, and toward 0 as <span class="math">\(\mathbf{p}\)</span> is more outside.</p>
|
||
|
||
<figure>
|
||
<img src="images/big-sigcat-winding-number.gif" alt="Example 702 computes the
|
||
generalized winding number function for a tetrahedral mesh inside a cat with
|
||
holes and self intersections (gold). The silver mesh is surface of the
|
||
extracted interior tets, and slices show the winding number function on all
|
||
tets in the convex hull: blue (~0), green (~1), yellow
|
||
(~2)." />
|
||
<figcaption>Example <a href="702_WindingNumber/main.cpp">702</a> computes the
|
||
generalized winding number function for a tetrahedral mesh inside a cat with
|
||
holes and self intersections (gold). The silver mesh is surface of the
|
||
extracted interior tets, and slices show the winding number function on all
|
||
tets in the convex hull: blue (~0), green (~1), yellow
|
||
(~2).</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="meshdecimation"><a href="#meshdecimation">Mesh Decimation</a></h2>
|
||
|
||
<p>The study of mesh simplification or <em>decimation</em> is nearly as old as meshes
|
||
themselves. Given a high resolution mesh with too many triangles, find a “well
|
||
approximating” low resolution mesh with far fewer triangles. By now there are a
|
||
variety of different paradigms for solving this problem and state-of-the-art
|
||
methods are fairly advanced.</p>
|
||
|
||
<p>One family of mesh decimation methods operates by successively remove elements
|
||
from the mesh. In particular, Hoppe advocates for successively remove or rather
|
||
collapsing edges <a class="citation" href="#fn:37" title="Jump to citation">[37]<span class="citekey" style="display:none">hoppe_1996</span></a>. The generic form of this technique is to
|
||
construct a sequence of n meshes from the initial high-resolution mesh <span class="math">\(M_0\)</span> to
|
||
the lowest resolution mesh <span class="math">\(M_n\)</span> by collapsing a single edge:</p>
|
||
|
||
<p><span class="math">\(M_0 \mathop{\longrightarrow}_\text{edge collapse}
|
||
M_1 \mathop{\longrightarrow}_\text{edge collapse}
|
||
\dots \mathop{\longrightarrow}_\text{edge collapse}
|
||
M_{n-1} \mathop{\longrightarrow}_\text{edge collapse} M_n.\)</span></p>
|
||
|
||
<p>Hoppe’s original method and subsequent follow-up works propose various ways to
|
||
choose the next edge to collapse in this sequence. Using a cost-based paradigm,
|
||
one can maintain a priority queue of edges based on their “cost” (how much
|
||
“worse” will my approximation be if I remove this edge?). The cheapest edge is
|
||
collapsed and costs of neighboring edges are updated.</p>
|
||
|
||
<p>In order to maintain the topology (e.g. if the mesh is combinatorially as
|
||
sphere or a torus etc.), one should assign infinite cost to edges whose
|
||
collapse would alter the mesh topology. Indeed this happens if and only if the
|
||
number of mutual neighbors of the endpoints of the collapsing edge is not
|
||
exactly two!</p>
|
||
|
||
<p>If there exists a third shared vertex, then another face will be removed, but 2
|
||
edges will be removed. This can result in unwanted holes or non-manifold
|
||
“flaps”.</p>
|
||
|
||
<figure>
|
||
<img src="images/edge-collapse.jpg" alt="A valid edge collapse and an invalid edge collapse." />
|
||
<figcaption>A valid edge collapse and an invalid edge collapse.</figcaption>
|
||
</figure>
|
||
|
||
<blockquote>
|
||
<p>There is also a one-off condition that no edges of a tetrahedron should be
|
||
collapsed.</p>
|
||
</blockquote>
|
||
|
||
<p>Because libigl (purposefully) does not center its implementations around a
|
||
dynamic mesh data structure (e.g. half-edge datastructure), support for
|
||
topology changes are limited. Nonetheless, libigl has support for isolated edge
|
||
collapses, sequences of edge-collapses (each in O(log) time) and priority queue
|
||
based decimation.</p>
|
||
|
||
<p>The simplest is <code>igl::decimation</code>. By calling</p>
|
||
|
||
<pre><code class="cpp">igl::decimate(V,F,1000,U,G);
|
||
</code></pre>
|
||
|
||
<p>the mesh <code>(V,F)</code> will be decimated to a new mesh <code>(U,G)</code> so that <code>G</code> has at
|
||
most <code>1000</code> faces. This uses default (naive) criteria for determining the cost
|
||
of an edge collapse and the placement of the merged vertex. Shortest edges are
|
||
collapsed first, and merged vertices are placed at edge midpoints.</p>
|
||
|
||
<p>One can also provide function handles (<code>c++</code> lambda functions are convenient
|
||
here) <code>cost_and_placement</code> and <code>stopping_condition</code> for determining the
|
||
cost/placement of an edge collapse and the stopping condition respectively. For
|
||
example, the default version above is implemented as:</p>
|
||
|
||
<pre><code class="cpp">igl::decimate(V,F,shortest_edge_and_midpoint,max_m,U,G);
|
||
</code></pre>
|
||
|
||
<p>where <code>shortest_edge_and_midpoint</code> assign the edge’s length as cost and its
|
||
midpoint as the merged vertex placement and <code>max_m</code> counts the current number
|
||
of faces (valid collapses decrease count by 2) and returns <code>true</code> if the count
|
||
drops below <code>m=1000</code>.</p>
|
||
|
||
<p>One can also scratch deeper inside the decimation loop and call
|
||
<code>igl::collapse_edge</code> directly. In order to operate efficiently, this routine
|
||
needs more than the usual <code>(V,F)</code> mesh representation. We need <code>E</code> a list of
|
||
edge indices, where <code>E.row(i) --> [s,d]</code>; we need <code>EMAP</code> which maps the
|
||
“half”-edges of each triangle in <code>F</code> to its corresponding edge in <code>E</code> so that
|
||
<code>E.row(EMAP(f+i*F.rows)) --> [s,d]</code> if the edge across from the ith corner of the
|
||
fth face is <code>[s,d]</code> (up to orientation); we need <code>EF</code> and <code>EI</code> which keep track
|
||
of the faces incident on each edge and across from which corner of those faces
|
||
the edges appears, so that <code>EF(e,o) = f</code> and <code>EI(e,o) = i</code> means that the edge
|
||
<code>E.row(e) --> [s,d]</code> appears in the fth face across from its ith corner (for
|
||
<code>o=0</code> the edge orientations should match, for <code>o=1</code> the orientations are
|
||
opposite).</p>
|
||
|
||
<p>When a collapse occurs, the sizes of the <code>F</code>,<code>E</code>, etc. matrices do not change.
|
||
Rather rows corresponding to “removed” faces and edges are set to a special
|
||
constant value <code>IGL_COLLAPSE_EDGE_NULL</code>. Doing this ensures that we’re able to
|
||
remove edges in truly constant time O(1).</p>
|
||
|
||
<blockquote>
|
||
<p>Conveniently <code>IGL_COLLAPSE_EDGE_NULL==0</code>. This means most OPENGL style renderings of <code>F</code>
|
||
will simply draw a bunch of 0-area triangles at the first vertex.</p>
|
||
</blockquote>
|
||
|
||
<p>The following will collapse the first
|
||
edge and place its merged vertex at the origin:</p>
|
||
|
||
<pre><code class="cpp">igl::collapse_edge(0,RowVector3d(0,0,0),V,F,E,EMAP,EF,EI);
|
||
</code></pre>
|
||
|
||
<p>If valid, then <code>V</code>,<code>F</code>,<code>E</code>,<code>EF</code>,<code>EI</code> are adjusted accordingly.</p>
|
||
|
||
<p>This is powerful, but low level. To build a decimator around this you’d need to
|
||
keep track which edges are left to collapse and which to collapse next.
|
||
Fortunately, libigl also exposes a priority queue based edge collapse with
|
||
function handles to adjust costs and placements.</p>
|
||
|
||
<p>The priority queue is implemented as a (ordered) set <code>Q</code> or (cost,edge index)
|
||
pairs and a list of iterators <code>Qit</code> so that <code>Qit[e]</code> reveals the iterator in
|
||
<code>Q</code> corresponding to the eth edge. Placements are stored in a #E list of
|
||
positions <code>C</code>. When the following is called:</p>
|
||
|
||
<pre><code class="cpp">igl::collapse_edge(cost_and_placement,V,F,E,EMAP,EF,EI,Q,Qit,C);
|
||
</code></pre>
|
||
|
||
<p>the lowest cost edge collapse according to <code>Q</code> is attempted. If valid, then
|
||
<code>V</code>,<code>F</code>,etc. are adjusted accordingly and that edge is “popped” from <code>Q</code>. Using
|
||
<code>Qit</code> its neighboring edges are also popped from <code>Q</code> and re-inserted after
|
||
updating their costs according to <code>cost_and_placement</code>, new placements are
|
||
remembered in <code>C</code>. If not valid, then the edge is “popped” from <code>Q</code> and
|
||
reinserted with infinite cost.</p>
|
||
|
||
<figure>
|
||
<img src="images/fertility-edge-collapse.gif" alt="Example 703 conducts edge collapses on the fertility
|
||
model." />
|
||
<figcaption>Example 703 conducts edge collapses on the fertility
|
||
model.</figcaption>
|
||
</figure>
|
||
|
||
<p>The <a href="./703_Decimation/main.cpp">Example 703</a> demonstrates using this priority
|
||
queue based approach with the simple shortest-edge-midpoint cost/placement
|
||
strategy discussed above.</p>
|
||
|
||
<h2 id="signeddistances"><a href="#signeddistances">Signed Distances</a></h2>
|
||
|
||
<p>In the <a href="#generalizedwindingnumber">Generalized Winding Number section</a>, we
|
||
examined a robust method for determining whether points lie inside or outside
|
||
of a given triangle soup mesh. Libigl complements this algorithm with
|
||
accelerated signed and unsigned distance queries and “in element” queries for
|
||
planar triangle meshes and 3D tetrahedral meshes. These routines make use of
|
||
libigl’s general purpose axis-aligned bounding box hierarchy (<code>igl/AABB.h</code>).
|
||
This class is lightweight and—by design—does not store a copy of the mesh
|
||
(taking it as inputs to its member functions instead).</p>
|
||
|
||
<h3 id="pointlocation">Point location</h3>
|
||
|
||
<p>For tetrahedral meshes, this is useful for “in element” or “point location”
|
||
queries: given a point <span class="math">\(\mathbf{q}\in\mathcal{R}^3\)</span> and a tetrahedral mesh
|
||
<span class="math">\((V,T)\)</span> determine in which tetrahedron <span class="math">\(\mathbf{q}\)</span> lies. This is accomplished
|
||
in libigl for a tet mesh <code>V,T</code> and a list of query points in the rows of <code>Q</code>
|
||
via the <code>igl::in_element()</code>:</p>
|
||
|
||
<pre><code class="cpp">// Initialize AABB tree
|
||
igl::AABB<MatrixXd,3> tree;
|
||
tree.init(V,T);
|
||
VectorXi I;
|
||
igl::in_element(V,T,Q,tree,I);
|
||
</code></pre>
|
||
|
||
<p>the resulting vector <code>I</code> is a list of indices into <code>T</code> revealing the <em>first</em>
|
||
tetrahedron found to contain the corresponding point in <code>Q</code>.</p>
|
||
|
||
<p>For overlapping meshes, a point <span class="math">\(\mathbf{q}\)</span> may belong to more than one
|
||
tetrahedron. In those cases, one can find them all (not just the first) by
|
||
using the <code>igl::in_element</code> overload with a <code>SparseMatrix</code> as the output:</p>
|
||
|
||
<pre><code class="cpp">SparseMatrix<int> I;
|
||
igl::in_element(V,T,Q,tree,I);
|
||
</code></pre>
|
||
|
||
<p>now each row of <code>I</code> reveals whether each tet contains the corresponding row in
|
||
<code>Q</code>: <code>I(q,e)!=0</code> means that point <code>q</code> is in element <code>e</code>.</p>
|
||
|
||
<h3 id="closestpoints">Closest points</h3>
|
||
|
||
<p>For Triangle meshes, we use the AABB tree to accelerate point-mesh closest
|
||
point queries: given a mesh <span class="math">\((V,F)\)</span> and a query point
|
||
<span class="math">\(\mathbf{q}\in\mathcal{R}^3\)</span> find the closest point <span class="math">\(\mathbf{c} \in (V,F)\)</span>
|
||
(where <span class="math">\(\mathbf{c}\)</span> is not necessarily a vertex of <span class="math">\((V,F)\)</span>). This is
|
||
accomplished for a triangle mesh <code>V,F</code> and a list of points in the rows of <code>P</code>
|
||
via <code>igl::point_mesh_squared_distance</code>:</p>
|
||
|
||
<pre><code class="cpp">VectorXd sqrD;
|
||
VectorXi I;
|
||
MatrixXd C;
|
||
igl::point_mesh_squared_distance(P,V,F,sqrD,I,C);
|
||
</code></pre>
|
||
|
||
<p>the output <code>sqrD</code> contains the (unsigned) squared distance from each point in
|
||
<code>P</code> to its closest point given in <code>C</code> which lies on the element in <code>F</code> given by
|
||
<code>I</code> (e.g. from which one could recover barycentric coordinates, using
|
||
<code>igl::barycentric_coordinates</code>).</p>
|
||
|
||
<p>If the mesh <code>V,F</code> is static, but the point set <code>P</code> is changing dynamically then
|
||
it’s best to reuse the AABB hierarchy that’s being built during
|
||
<code>igl::point_mesh_squared_distance</code>:</p>
|
||
|
||
<pre><code class="cpp">igl::AABB tree;
|
||
tree.init(V,F);
|
||
tree.squared_distance(V,F,P,sqrD,I,C);
|
||
... // P changes, but (V,F) does not
|
||
tree.squared_distance(V,F,P,sqrD,I,C);
|
||
</code></pre>
|
||
|
||
<h3 id="signeddistance">Signed distance</h3>
|
||
|
||
<p>Finally, from the closest point or the winding number it’s possible to <em>sign</em>
|
||
this distance. In <code>igl::signed_distance</code> we provide two methods for signing:
|
||
the so-called “pseudo-normal test” <a class="citation" href="#fn:38" title="Jump to citation">[38]<span class="citekey" style="display:none">baerentzen_2005</span></a> and the generalized
|
||
winding number <a class="citation" href="#fn:36" title="Jump to citation">[36]<span class="citekey" style="display:none">jacobson_2013</span></a>.</p>
|
||
|
||
<p>The pseudo-normal test (see also <code>igl::pseudonormal_test</code>) assumes the input
|
||
mesh is a watertight (closed, non-self-intersecting, manifold) mesh. Then given
|
||
a query point <span class="math">\(\mathbf{q}\)</span> and its closest point <span class="math">\(\mathbf{c} \in (V,F)\)</span>, it
|
||
carefully chooses an outward normal <span class="math">\(\mathbf{n}\)</span> at <span class="math">\(\mathbf{c}\)</span> so that
|
||
<span class="math">\(\text{sign}(\mathbf{q}-\mathbf{c})\cdot \mathbf{n}\)</span> reveals whether
|
||
<span class="math">\(\mathbf{q}\)</span> is inside <span class="math">\((V,F)\)</span>: –1, or outside: +1. This is a fast <span class="math">\(O(1)\)</span> test
|
||
once <span class="math">\(\mathbf{c}\)</span> is located, but may fail if <code>V,F</code> is not watertight.</p>
|
||
|
||
<p>An alternative is to use the <a href="#generalizedwindingnumber">generalized winding
|
||
number</a> to determine the sign. This is very robust to
|
||
unclean meshes <code>V,F</code> but slower: something like <span class="math">\(O(\sqrt{n})\)</span> once <span class="math">\(\mathbf{c}\)</span>
|
||
is located.</p>
|
||
|
||
<p>In either case, the interface via <code>igl::signed_distance</code> is:</p>
|
||
|
||
<pre><code class="cpp">// Choose type of signing to use
|
||
igl::SignedDistanceType type = SIGNED_DISTANCE_TYPE_PSEUDONORMAL;
|
||
igl::signed_distance(P,V,F,sign_type,S,I,C,N);
|
||
</code></pre>
|
||
|
||
<p>the outputs are as above for <code>igl::point_mesh_squared_distance</code> but now <code>S</code>
|
||
contains signed (unsquared) distances and the extra output <code>N</code> (only set when
|
||
<code>type == SIGNED_DISTANCE_TYPE_PSEUDON</code>) contains the normals used for signing
|
||
with the pseudo-normal test.</p>
|
||
|
||
<figure>
|
||
<img src="images/bunny-signed-distance.gif" alt="Example 704 computes signed distance on
|
||
slices through the bunny." />
|
||
<figcaption>Example <a href="704_SignedDistance/main.cpp">704</a> computes signed distance on
|
||
slices through the bunny.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="marchingcubes"><a href="#marchingcubes">Marching Cubes</a></h2>
|
||
|
||
<p>Often 3D data is captured as scalar field defined over space <span class="math">\(f(\mathbf{x}) :
|
||
\mathcal{R}^3 \rightarrow \mathcal{R}\)</span>. Lurking within this field,
|
||
<em>iso-surfaces</em> of the scalar field are often salient geometric objects. The
|
||
iso-surface at value <span class="math">\(v\)</span> is composed of all points <span class="math">\(\mathbf{x}\)</span> in
|
||
<span class="math">\(\mathcal{R}^3\)</span> such that <span class="math">\(f(\mathbf{x}) = v\)</span>. A core problem in geometry
|
||
processing is to extract an iso-surface as a triangle mesh for further
|
||
mesh-based processing or visualization. This is referred to as iso-contouring.</p>
|
||
|
||
<p>“Marching Cubes” <a class="citation" href="#fn:39" title="Jump to citation">[39]<span class="citekey" style="display:none">lorensen_1987</span></a> is a <a href="https://en.wikipedia.org/wiki/Marching_cubes">famous
|
||
method</a> for iso-contouring
|
||
tri-linear functions <span class="math">\(f\)</span> on a regular lattice (aka grid). The core idea of this
|
||
method is to contour the iso-surface passing through each cell (if it does at
|
||
all) with a predefined topology (aka connectivity) chosen from a look up table
|
||
depending on the function values at each vertex of the cell. The method
|
||
iterates (“marches”) over all cells (“cubes”) in the grid and stitches together
|
||
the final, watertight mesh.</p>
|
||
|
||
<p>In libigl, <code>igl::marching_cubes</code> constructs a triangle mesh <code>(V,F)</code> from an
|
||
input scalar field <code>S</code> sampled at vertex locations <code>GV</code> of a <code>nx</code> by <code>ny</code> by
|
||
<code>nz</code> regular grid:</p>
|
||
|
||
<pre><code class="cpp">igl::marching_cubes(S,GV,nx,ny,nz,V,F);
|
||
</code></pre>
|
||
|
||
<figure>
|
||
<img src="images/armadillo-marching-cubes.jpg" alt="(Example 705) samples signed distance to the
|
||
input mesh (left) and then reconstructs the surface using
|
||
marching cubes to contour the 0-level set (center). For comparison, clamping
|
||
this signed distance field to an indicator function and contouring reveals
|
||
serious aliasing artifacts." />
|
||
<figcaption>(<a href="705_MarchingCubes/main.cpp">Example 705</a>) samples signed distance to the
|
||
input mesh (left) and then reconstructs the surface using
|
||
marching cubes to contour the 0-level set (center). For comparison, clamping
|
||
this signed distance field to an indicator function and contouring reveals
|
||
serious aliasing artifacts.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="facetorientation"><a href="#facetorientation">Facet Orientation</a></h2>
|
||
|
||
<p>Models from the web occasionally arrive <em>unorientated</em> in the sense that
|
||
the orderings of each triangles vertices do not consistently agree. Determining
|
||
a consistent facet orientation for a mesh is essential for two-sided lighting
|
||
(e.g., a cloth with red velvet on one side and gold silk on the other side) and
|
||
for inside-outside determination(e.g., using <a href="#generalizedwindingnumber">generalized winding
|
||
numbers</a>).</p>
|
||
|
||
<p>For (open) surfaces representing two-sided sheets, libigl provides a routine to
|
||
force consistent orientations within each orientable patch
|
||
(<code>igl::orientable_patches</code>) of a mesh:</p>
|
||
|
||
<pre><code class="cpp">igl::bfs_orient(F,FF,C);
|
||
</code></pre>
|
||
|
||
<p>This simple routine will use breadth-first search on each patch of the mesh to
|
||
enforce a consistent facet orientation in the output faces <code>FF</code>.</p>
|
||
|
||
<p>For (closed or nearly closed) surfaces representing the boundary of a solid
|
||
object, libigl provides a routine to reorient faces so that the vertex ordering
|
||
corresponds to a counter-clockwise ordering of the vertices with a
|
||
right-hand-rule normal pointing outward. This method <a class="citation" href="#fn:40" title="Jump to citation">[40]<span class="citekey" style="display:none">takayama14</span></a> assumes
|
||
that <a href="https://www.reddit.com/r/askscience/comments/32otgx/which_as_a_is_more_empty_an_atom_or_the_universe/">most of the universe is
|
||
empty</a>.
|
||
That is, most points in space are outside of the solid object than inside.
|
||
Points are sampled over surface patches. For each sample point, rays are shot
|
||
into both hemispheres to compute average of the (distance weighted) ambient
|
||
occlusion on each side. A patch is oriented so that the outward side is <em>less
|
||
occluded</em> (lighter, i.e., facing more void space).</p>
|
||
|
||
<pre><code class="cpp">igl::embree::reorient_facets_raycast(V,F,FF,I);
|
||
</code></pre>
|
||
|
||
<p>The boolean vector <code>I</code> reveals which rows of <code>F</code> have been flipped in <code>FF</code>.</p>
|
||
|
||
<figure>
|
||
<img src="images/truck-facet-orientation.jpg" alt="(Example 706) loads a truck model with
|
||
inconsistent orientations (back facing triangles shown darker). Orientable
|
||
patches are uniquely colored and then oriented to face outward (middle left).
|
||
Alternatively, each individual triangle is considered a patch (middle right)
|
||
and oriented outward independently." />
|
||
<figcaption>(<a href="706_FacetOrientation/main.cpp">Example 706</a>) loads a truck model with
|
||
inconsistent orientations (back facing triangles shown darker). Orientable
|
||
patches are uniquely colored and then oriented to face outward (middle left).
|
||
Alternatively, each individual triangle is considered a “patch” (middle right)
|
||
and oriented outward independently.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="sweptvolume"><a href="#sweptvolume">Swept Volume</a></h2>
|
||
|
||
<p>The swept volume <span class="math">\(S\)</span> of a moving solid object <span class="math">\(A\)</span> can be defined as any point in
|
||
space such that at one moment in time the point lies inside the solid. In other
|
||
words, it is the union of the solid object transformed by the rigid motion
|
||
<span class="math">\(f(t)\)</span> over time:</p>
|
||
|
||
<p><span class="math">\(S = \bigcup \limits_{t\in [0,1]} f(t) A.\)</span></p>
|
||
|
||
<p>The surface of the swept volume of a solid bounded by a triangle mesh
|
||
undergoing a rigid motion with non-trivial rotation is <em><strong>not</strong></em> a surface
|
||
exactly representably by triangle mesh: it will be a piecewise-ruled surface.</p>
|
||
|
||
<p>To see this, consider the surface swept by a single edge’s line segment as it
|
||
performs a screw motion.</p>
|
||
|
||
<p>This means that if we’d like to the surface of the swept volume of a triangle
|
||
mesh undergoing a rigid motion and we’d like the output to be another triangle
|
||
mesh, then we’re going to have to be happy with some amount of approximation
|
||
error.</p>
|
||
|
||
<p>With this in mind, the simplest method for computing an approximate swept
|
||
volume is by exploiting an alternative definition of the swept volume based on
|
||
signed distances:</p>
|
||
|
||
<p><span class="math">\(S = \left\{ \mathbf{p}\ \middle| \ d(\mathbf{p},\partial S) < 0 \right\} = \left\{ \mathbf{p}\
|
||
\middle|\
|
||
\min\limits_{t \in [0,1]} d(\mathbf{p},f(t)\ \partial A) < 0 \right\}\)</span></p>
|
||
|
||
<p>If <span class="math">\(\partial A\)</span> is a triangle mesh, then we can approximate this by 1)
|
||
discretizing time at a finite step of steps <span class="math">\([0,\Delta t,2\Delta t, \dots, 1]\)</span>
|
||
and by 2) discretizing space with a regular grid and representing the distance
|
||
field using trilinear interpolation of grid values. Finally the output mesh,
|
||
<span class="math">\(\partial S\)</span> is approximated by contouring using Marching Cubes
|
||
<a class="citation" href="#fn:39" title="Jump to citation">[39]<span class="citekey" style="display:none">lorensen_1987</span></a>.</p>
|
||
|
||
<p>This method is similar to one described by Schroeder et al. in 1994
|
||
<a class="citation" href="#fn:41" title="Jump to citation">[41]<span class="citekey" style="display:none">schroeder_1994</span></a>, and the one used in conjunction with boolean operations by
|
||
Garg et al. 2016 <a class="citation" href="#fn:42" title="Jump to citation">[42]<span class="citekey" style="display:none">garg_2016</span></a>.</p>
|
||
|
||
<p>In libigl, if your input solid’s surface is represented by <code>(V,F)</code> then the
|
||
output surface mesh will be <code>(SV,SF)</code> after calling:</p>
|
||
|
||
<pre><code class="cpp">igl::copyleft::swept_volume(V,F,num_time_steps,grid_size,isolevel,SV,SF);
|
||
</code></pre>
|
||
|
||
<p>The <code>isolevel</code> parameter can be set to zero to approximate the exact swept
|
||
volume, greater than zero to approximate a positive offset of the swept volume
|
||
or less than zero to approximate a negative offset.</p>
|
||
|
||
<figure>
|
||
<img src="images/bunny-swept-volume.gif" alt="(Example 707) computes
|
||
the surface of the swept volume (silver) of the bunny model undergoing a rigid
|
||
motion (gold)." />
|
||
<figcaption>(<a href="707_SweptVolume/main.cpp">Example 707</a>) computes
|
||
the surface of the swept volume (silver) of the bunny model undergoing a rigid
|
||
motion (gold).</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="pickingverticesandfaces"><a href="#pickingverticesandfaces">Picking</a></h2>
|
||
|
||
<p>Picking vertices and faces using the mouse is very common in geometry
|
||
processing applications. While this might seem a simple operation, its
|
||
implementation is not straighforward. Libigl contains a function that solves this problem using the
|
||
<a href="https://software.intel.com/en-us/articles/embree-photo-realistic-ray-tracing-kernels">Embree</a>
|
||
raycaster. Its usage is demonstrated in <a href="708_Picking/main.cpp">Example 708</a>:</p>
|
||
|
||
<pre><code class="cpp">bool hit = igl::unproject_onto_mesh(
|
||
Vector2f(x,y),
|
||
F,
|
||
viewer.core.view * viewer.core.model,
|
||
viewer.core.proj,
|
||
viewer.core.viewport,
|
||
*ei,
|
||
fid,
|
||
vid);
|
||
</code></pre>
|
||
|
||
<p>This function casts a ray from the view plane in the view direction. Variables
|
||
<code>x</code> and <code>y</code> are
|
||
the mouse screen coordinates; <code>view</code>, <code>model</code>, <code>proj</code> are the view, model and
|
||
projection matrix respectively; <code>viewport</code> is the viewport in OpenGL format;
|
||
<code>ei</code>
|
||
contains a <a href="http://en.wikipedia.org/wiki/Bounding_volume_hierarchy">Bounding Volume
|
||
Hierarchy</a> constructed
|
||
by Embree, and <code>fid</code> and <code>vid</code> are the picked face and vertex, respectively.</p>
|
||
|
||
<figure>
|
||
<img src="images/607_Picking.png" alt="(Example 708) Picking via ray casting. The selected
|
||
vertices are colored in red." />
|
||
<figcaption>(<a href="708_Picking/main.cpp">Example 708</a>) Picking via ray casting. The selected
|
||
vertices are colored in red.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="vectorfieldvisualizer"><a href="#vectorfieldvisualizer">Vector Field Visualization</a></h2>
|
||
|
||
<p>Vector fields on surfaces are commonly visualized by tracing <a href="https://en.wikipedia.org/wiki/Streamlines,_streaklines,_and_pathlines">streamlines</a>. Libigl
|
||
supports the seeding and tracing of streamlines, for both simple vector fields
|
||
and for N-rosy fields. The seeds for the streamlines are initialized using <code>streamlines_init</code>,
|
||
and the lines are traced using <code>streamlines_next</code>. Each call to <code>streamlines_next</code> extends
|
||
each line by one triangle, allowing interactive rendering of the traced lines, as demonstrated
|
||
in <a href="709_VectorFieldVisualizer/main.cpp">Example 709</a>.</p>
|
||
|
||
<figure>
|
||
<img src="images/streamlines.jpg" alt="(Example 709) Interactive streamlines tracing." />
|
||
<figcaption>(<a href="709_VectorFieldVisualizer/main.cpp">Example 709</a>) Interactive streamlines tracing.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="slim"><a href="#slim">Scalable Locally Injective Maps</a></h2>
|
||
|
||
<p>The Scalable Locally Injective Maps <a class="citation" href="#fn:43" title="Jump to citation">[43]<span class="citekey" style="display:none">rabinovich_2016</span></a> algorithm allows to
|
||
compute locally injective maps on massive datasets. The algorithm shares many
|
||
similarities with ARAP, but uses a reweighting scheme to minimize arbitrary
|
||
distortion energies, including those that prevent the introduction of flips.</p>
|
||
|
||
<p><a href="710_SLIM/main.cpp">Example 710</a> contains three demos: (1) an example of large
|
||
scale 2D parametrization, (2) an example of 2D deformation with soft
|
||
constraints, and (3) an example of 3D deformation with soft constraints. The
|
||
implementation in libigl is self-contained and relies on Eigen for the solution
|
||
of the linear system used in the global step. An optimized version that relies
|
||
on Pardiso is available
|
||
<a href="https://github.com/MichaelRabinovich/Scalable-Locally-Injective-Mappings">here</a>.</p>
|
||
|
||
<figure>
|
||
<img src="images/slim.png" alt="A locally injective parametrization of a mesh with 50k faces is computed
|
||
using the SLIM algorithm in 10 iterations." />
|
||
<figcaption>A locally injective parametrization of a mesh with 50k faces is computed
|
||
using the SLIM algorithm in 10 iterations.</figcaption>
|
||
</figure>
|
||
|
||
<h2 id="subdivision"><a href="#subdivision">Subdivision surfaces</a></h2>
|
||
|
||
<p>Given a coarse mesh (aka cage) with vertices <code>V</code> and faces <code>F</code>, one can createa
|
||
higher-resolution mesh with more vertices and faces by <em>subdividing</em> every
|
||
face. That is, each coarse triangle in the input is replaced by many smaller
|
||
triangles. Libigl has three different methods for subdividing a triangle mesh.</p>
|
||
|
||
<p>An “in plane” subdivision method will not change the point set or carrier
|
||
surface of the mesh. New vertices are added on the planes of existing triangles
|
||
and vertices surviving from the original mesh are not moved.</p>
|
||
|
||
<p>By adding new faces, a subdivision algorithm changes the <em>combinatorics</em> of the
|
||
mesh. The change in combinatorics and the formula for positioning the
|
||
high-resolution vertices is called the “subdivision rule”.</p>
|
||
|
||
<p>For example, in the <em>in plane</em> subdivision method of <code>igl::upsample</code>, vertices
|
||
are added at the midpoint of every edge: <span class="math">\(v_{ab} = \frac{1}{2}(v_a + v_b)\)</span> and
|
||
each triangle <span class="math">\((i_a,i_b,i_c)\)</span> is replaced with four triangles:
|
||
<span class="math">\((i_a,i_{ab},i_{ca})\)</span>, <span class="math">\((i_b,i_{bc},i_{ab})\)</span>, <span class="math">\((i_{ab},i_{bc},i_{ca})\)</span>, and
|
||
<span class="math">\((i_{bc},i_{c},i_{ca})\)</span>. This process may be applied recursively, resulting in
|
||
a finer and finer mesh.</p>
|
||
|
||
<p>The subdivision method of <code>igl::loop</code> is not in plane. The vertices of the
|
||
refined mesh are moved to weight combinations of their neighbors: the mesh is
|
||
smoothed as it is refined <a class="citation" href="#fn:44" title="Jump to citation">[44]<span class="citekey" style="display:none">loop_1987</span></a>. This and other <em>smooth subdivision</em>
|
||
methods can be understood as generalizations of spline curves to surfaces. In
|
||
particular the Loop subdivision method will converge to a <span class="math">\(C^1\)</span> surface as we
|
||
consider the limit of recursive applications of subdivision. Away from
|
||
“irregular” or “extraordinary” vertices (vertices of the original cage with
|
||
valence not equal to 6), the surface is <span class="math">\(C^2\)</span>. The combinatorics (connectivity
|
||
and number of faces) of <code>igl::loop</code> and <code>igl::upsample</code> are identical: the only
|
||
difference is that the vertices have been smoothed in <code>igl::loop</code>.</p>
|
||
|
||
<p>Finally, libigl also implements a form of <em>in plane</em> “false barycentric
|
||
subdivision” in <code>igl::false_barycentric_subdivision</code>. This method simply adds
|
||
the barycenter of every triangle as a new vertex <span class="math">\(v_{abc}\)</span> and replaces each
|
||
triangle with three triangles <span class="math">\((i_a,i_b,i_{abc})\)</span>, <span class="math">\((i_b,i_c,i_{abc})\)</span>, and
|
||
<span class="math">\((i_c,i_a,i_{abc})\)</span>. In contrast to <code>igl::upsample</code>, this method will create
|
||
triangles with smaller and smaller internal angles and new vertices will sample
|
||
the carrier surfaces with extreme bias.</p>
|
||
|
||
<figure>
|
||
<img src="images/decimated-knight-subdivision.gif" alt="The original coarse mesh and three different subdivision methods:
|
||
igl::upsample, igl::loop and
|
||
igl::false_barycentric_subdivision." />
|
||
<figcaption>The original coarse mesh and three different subdivision methods:
|
||
<code>igl::upsample</code>, <code>igl::loop</code> and
|
||
<code>igl::false_barycentric_subdivision</code>.</figcaption>
|
||
</figure>
|
||
|
||
<h1 id="chapter7:miscellaneous">Miscellaneous</h1>
|
||
|
||
<p>Libigl contains a <em>wide</em> variety of geometry processing tools and functions for
|
||
dealing with meshes and the linear algebra related to them: far too many to
|
||
discuss in this introductory tutorial. We’ve pulled out a couple of the
|
||
interesting functions in this chapter to highlight.</p>
|
||
|
||
<h1 id="future">Outlook for continuing development</h1>
|
||
|
||
<p>Libigl is in active development, and we plan to focus on the following features
|
||
in the next months:</p>
|
||
|
||
<ul>
|
||
<li><p>A better and more consistent <strong>documentation</strong>, plus extending this tutorial
|
||
to cover more libigl features.</p></li>
|
||
<li><p>Implement a <strong>mixed-integer solver</strong> which only uses Eigen to remove the
|
||
dependency on CoMiSo.</p></li>
|
||
<li><p>Improve the robustness and performance of the active set QP solver. In
|
||
particular, handle linearly dependent constraints.</p></li>
|
||
<li><p>Implement more mesh analysis functions, including structural analysis for
|
||
masonry and <em>3D-printability</em> analysis.</p></li>
|
||
<li><p>Increase support for point clouds and general polygonal meshes.</p></li>
|
||
<li><p>What would you like to see in libigl? <a href="mailto:alecjacobson@gmail.com">Contact
|
||
us!</a> or post a <a href="https://github.com/libigl/libigl/issues/new">feature
|
||
request</a>.</p></li>
|
||
</ul>
|
||
|
||
<p>We encourage you to contribute to the library and to report problems and bugs.
|
||
The best way to contribute new feature or bug fixes is to fork the libigl
|
||
repository and to open a <a href="https://help.github.com/articles/using-pull-requests">pull
|
||
request</a> on <a href="https://github.com/libigl/libigl">our github
|
||
repository</a>.</p>
|
||
|
||
<div class="footnotes">
|
||
<hr />
|
||
<ol>
|
||
|
||
<li id="fn:1" class="citation"><span class="citekey" style="display:none">meyer_2003</span><p>Mark Meyer, Mathieu Desbrun, Peter Schröder and Alan H. Barr,
|
||
<a href="https://www.google.com/search?q=Discrete+Differential-Geometry+Operators+for+Triangulated+2-Manifolds">Discrete Differential-Geometry Operators for Triangulated
|
||
2-Manifolds</a>,
|
||
2003.</p>
|
||
</li>
|
||
|
||
<li id="fn:2" class="citation"><span class="citekey" style="display:none">panozzo_2010</span><p>Daniele Panozzo, Enrico Puppo, Luigi Rocca, <a href="https://www.google.com/search?q=Efficient+Multi-scale+Curvature+and+Crease+Estimation">Efficient
|
||
Multi-scale Curvature and Crease
|
||
Estimation</a>,
|
||
2010.</p>
|
||
</li>
|
||
|
||
<li id="fn:3" class="citation"><span class="citekey" style="display:none">jacobson_thesis_2013</span><p>Alec Jacobson,
|
||
<a href="https://www.google.com/search?q=Algorithms+and+Interfaces+for+Real-Time+Deformation+of+2D+and+3D+Shapes"><em>Algorithms and Interfaces for Real-Time Deformation of 2D and 3D
|
||
Shapes</em></a>,
|
||
2013.</p>
|
||
</li>
|
||
|
||
<li id="fn:4" class="citation"><span class="citekey" style="display:none">sharf_2007</span><p>Andrei Sharf, Thomas Lewiner, Gil Shklarski, Sivan Toledo, and
|
||
Daniel Cohen-Or. <a href="https://www.google.com/search?q=Interactive+topology-aware+surface+reconstruction">Interactive topology-aware surface
|
||
reconstruction</a>,
|
||
2007.</p>
|
||
</li>
|
||
|
||
<li id="fn:5" class="citation"><span class="citekey" style="display:none">kazhdan_2012</span><p>Michael Kazhdan, Jake Solomon, Mirela Ben-Chen,
|
||
<a href="https://www.google.com/search?q=Can+Mean-Curvature+Flow+Be+Made+Non-Singular">Can Mean-Curvature Flow Be Made
|
||
Non-Singular</a>,
|
||
2012.</p>
|
||
</li>
|
||
|
||
<li id="fn:6" class="citation"><span class="citekey" style="display:none">rustamov_2011</span><p>Raid M. Rustamov, <a href="https://www.google.com/search?q=Multiscale+Biharmonic+Kernels">Multiscale Biharmonic
|
||
Kernels</a>, 2011.</p>
|
||
</li>
|
||
|
||
<li id="fn:7" class="citation"><span class="citekey" style="display:none">vallet_2008</span><p>Bruno Vallet and Bruno Lévy. <a href="https://www.google.com/search?q=Spectral+Geometry+Processing+with+Manifold+Harmonics">Spectral Geometry Processing with
|
||
Manifold
|
||
Harmonics</a>,
|
||
2008.</p>
|
||
</li>
|
||
|
||
<li id="fn:8" class="citation"><span class="citekey" style="display:none">hildebrandt_2011</span><p>Klaus Hildebrandt, Christian Schulz, Christoph von
|
||
Tycowicz, and Konrad Polthier. <a href="https://www.google.com/search?q=Interactive+Surface+Modeling+using+Modal+Analysis">Interactive Surface Modeling using Modal
|
||
Analysis</a>,
|
||
2011.</p>
|
||
</li>
|
||
|
||
<li id="fn:9" class="citation"><span class="citekey" style="display:none">barbic_2005</span><p>Jernej Barbic and Doug James. <a href="https://www.google.com/search?q=Real-Time+Subspace+Integration+for+St.Venant-Kirchhoff+Deformable+Models">Real-Time Subspace Integration
|
||
for St.Venant-Kirchhoff Deformable
|
||
Models</a>,
|
||
2005.</p>
|
||
</li>
|
||
|
||
<li id="fn:10" class="citation"><span class="citekey" style="display:none">botsch_2004</span><p>Matrio Botsch and Leif Kobbelt.
|
||
<a href="https://www.google.com/search?q=An+Intuitive+Framework+for+Real-Time+Freeform+Modeling">An Intuitive Framework for Real-Time Freeform
|
||
Modeling</a>,
|
||
2004.</p>
|
||
</li>
|
||
|
||
<li id="fn:11" class="citation"><span class="citekey" style="display:none">jacobson_mixed_2010</span><p>Alec Jacobson, Elif Tosun, Olga Sorkine, and Denis
|
||
Zorin. <a href="https://www.google.com/search?q=Mixed+Finite+Elements+for+Variational+Surface+Modeling">Mixed Finite Elements for Variational Surface
|
||
Modeling</a>,
|
||
2010.</p>
|
||
</li>
|
||
|
||
<li id="fn:12" class="citation"><span class="citekey" style="display:none">sorkine_2004</span><p>Olga Sorkine, Yaron Lipman, Daniel Cohen-Or, Marc Alexa,
|
||
Christian Rössl and Hans-Peter Seidel. <a href="https://www.google.com/search?q=Laplacian+Surface+Editing">Laplacian Surface
|
||
Editing</a>, 2004.</p>
|
||
</li>
|
||
|
||
<li id="fn:13" class="citation"><span class="citekey" style="display:none">jacobson_2011</span><p>Alec Jacobson, Ilya Baran, Jovan Popović, and Olga Sorkine.
|
||
<a href="https://www.google.com/search?q=Bounded+biharmonic+weights+for+real-time+deformation">Bounded Biharmonic Weights for Real-Time
|
||
Deformation</a>,
|
||
2011.</p>
|
||
</li>
|
||
|
||
<li id="fn:14" class="citation"><span class="citekey" style="display:none">kavan_2008</span><p>Ladislav Kavan, Steven Collins, Jiri Zara, and Carol O’Sullivan.
|
||
<a href="https://www.google.com/search?q=Geometric+Skinning+with+Approximate+Dual+Quaternion+Blending">Geometric Skinning with Approximate Dual Quaternion
|
||
Blending</a>,
|
||
2008.</p>
|
||
</li>
|
||
|
||
<li id="fn:15" class="citation"><span class="citekey" style="display:none">sorkine_2007</span><p>Olga Sorkine and Marc Alexa. <a href="https://www.google.com/search?q=As-rigid-as-possible+Surface+Modeling">As-rigid-as-possible Surface
|
||
Modeling</a>, 2007.</p>
|
||
</li>
|
||
|
||
<li id="fn:16" class="citation"><span class="citekey" style="display:none">chao_2010</span><p>Isaac Chao, Ulrich Pinkall, Patrick Sanan, Peter Schröder.
|
||
<a href="https://www.google.com/search?q=A+Simple+Geometric+Model+for+Elastic+Deformations">A Simple Geometric Model for Elastic
|
||
Deformations</a>,
|
||
2010.</p>
|
||
</li>
|
||
|
||
<li id="fn:17" class="citation"><span class="citekey" style="display:none">mcadams_2011</span><p>Alexa McAdams, Andrew Selle, Rasmus Tamstorf, Joseph Teran,
|
||
Eftychios Sifakis. <a href="https://www.google.com/search?q=Computing+the+Singular+Value+Decomposition+of+3x3+matrices+with+minimal+branching+and+elementary+floating+point+operations">Computing the Singular Value Decomposition of 3x3
|
||
matrices with minimal branching and elementary floating point
|
||
operations</a>,
|
||
2011.</p>
|
||
</li>
|
||
|
||
<li id="fn:18" class="citation"><span class="citekey" style="display:none">jacobson_2012</span><p>Alec Jacobson, Ilya Baran, Ladislav Kavan, Jovan Popović, and
|
||
Olga Sorkine. <a href="https://www.google.com/search?q=Fast+Automatic+Skinning+Transformations">Fast Automatic Skinning
|
||
Transformations</a>,
|
||
2012.</p>
|
||
</li>
|
||
|
||
<li id="fn:19" class="citation"><span class="citekey" style="display:none">jacobson_skinning_course_2014</span><p>Alec Jacobson, Zhigang Deng, Ladislav Kavan,
|
||
J.P. Lewis. <a href="https://www.google.com/search?q=Skinning+Real-Time+Shape+Deformation"><em>Skinning: Real-Time Shape
|
||
Deformation</em></a>,
|
||
2014.</p>
|
||
</li>
|
||
|
||
<li id="fn:20" class="citation"><span class="citekey" style="display:none">wang_bc_2015</span><p>Yu Wang, Alec Jacobson, Jernej Barbic, Ladislav Kavan. <a href="https://www.google.com/search?q=Linear+Subspace+Design+for+Real-Time+Shape+Deformation">Linear
|
||
Subspace Design for Real-Time Shape
|
||
Deformation</a>,
|
||
2015</p>
|
||
</li>
|
||
|
||
<li id="fn:21" class="citation"><span class="citekey" style="display:none">eck_2005</span><p>Matthias Eck, Tony DeRose, Tom Duchamp, Hugues Hoppe, Michael Lounsbery, Werner
|
||
Stuetzle. <a href="http://research.microsoft.com/en-us/um/people/hoppe/mra.pdf">Multiresolution Analysis of Arbitrary
|
||
Meshes</a>, 2005.</p>
|
||
</li>
|
||
|
||
<li id="fn:22" class="citation"><span class="citekey" style="display:none">levy_2002</span><p>Bruno Lévy, Sylvain Petitjean, Nicolas Ray, Jérome Maillot.
|
||
<a href="http://www.cs.jhu.edu/~misha/Fall09/Levy02.pdf">Least Squares Conformal Maps, for Automatic Texture Atlas
|
||
Generation,</a>, 2002.</p>
|
||
</li>
|
||
|
||
<li id="fn:23" class="citation"><span class="citekey" style="display:none">mullen_2008</span><p>Patrick Mullen, Yiying Tong, Pierre Alliez, Mathieu Desbrun.
|
||
<a href="http://www.geometry.caltech.edu/pubs/MTAD08.pdf">Spectral Conformal
|
||
Parameterization</a>, 2008.</p>
|
||
</li>
|
||
|
||
<li id="fn:24" class="citation"><span class="citekey" style="display:none">liu_2008</span><p>Ligang Liu, Lei Zhang, Yin Xu, Craig Gotsman, Steven J. Gortler.
|
||
<a href="http://cs.harvard.edu/~sjg/papers/arap.pdf">A Local/Global Approach to Mesh
|
||
Parameterization</a>, 2008.</p>
|
||
</li>
|
||
|
||
<li id="fn:25" class="citation"><span class="citekey" style="display:none">levy_2008</span><p>Nicolas Ray, Bruno Vallet, Wan Chiu Li, Bruno Lévy.
|
||
<a href="http://alice.loria.fr/publications/papers/2008/DGF/NSDFD-TOG.pdf">N-Symmetry Direction Field
|
||
Design</a>,
|
||
2008.</p>
|
||
</li>
|
||
|
||
<li id="fn:26" class="citation"><span class="citekey" style="display:none">bommes_2009</span><p>David Bommes, Henrik Zimmer, Leif Kobbelt.
|
||
<a href="http://www-sop.inria.fr/members/David.Bommes/publications/miq.pdf">Mixed-integer
|
||
quadrangulation</a>,
|
||
2009.</p>
|
||
</li>
|
||
|
||
<li id="fn:27" class="citation"><span class="citekey" style="display:none">knoppel_2013</span><p>Felix Knöppel, Keenan Crane, Ulrich Pinkall, and Peter
|
||
Schröder. <a href="http://www.cs.columbia.edu/~keenan/Projects/GloballyOptimalDirectionFields/paper.pdf">Globally Optimal Direction
|
||
Fields</a>,
|
||
2013.</p>
|
||
</li>
|
||
|
||
<li id="fn:28" class="citation"><span class="citekey" style="display:none">panozzo_2014</span><p>Daniele Panozzo, Enrico Puppo, Marco Tarini, Olga
|
||
Sorkine-Hornung. <a href="http://cs.nyu.edu/~panozzo/papers/frame-fields-2014.pdf">Frame Fields: Anisotropic and Non-Orthogonal Cross
|
||
Fields</a>,
|
||
2014.</p>
|
||
</li>
|
||
|
||
<li id="fn:29" class="citation"><span class="citekey" style="display:none">diamanti_2014</span><p>Olga Diamanti, Amir Vaxman, Daniele Panozzo, Olga
|
||
Sorkine-Hornung. <a href="http://igl.ethz.ch/projects/complex-roots/">Designing N-PolyVector Fields with Complex
|
||
Polynomials</a>, 2014</p>
|
||
</li>
|
||
|
||
<li id="fn:30" class="citation"><span class="citekey" style="display:none">liu_2011</span><p>Yang Liu, Weiwei Xu, Jun Wang, Lifeng Zhu, Baining Guo, Falai Chen, Guoping
|
||
Wang. <a href="http://research.microsoft.com/en-us/um/people/yangliu/publication/cdf.pdf">General Planar Quadrilateral Mesh Design Using Conjugate Direction
|
||
Field</a>,
|
||
2008.</p>
|
||
</li>
|
||
|
||
<li id="fn:31" class="citation"><span class="citekey" style="display:none">bouaziz_2012</span><p>Sofien Bouaziz, Mario Deuss, Yuliy Schwartzburg, Thibaut Weise, Mark Pauly
|
||
<a href="http://lgg.epfl.ch/publications/2012/shapeup.pdf">Shape-Up: Shaping Discrete Geometry with
|
||
Projections</a>, 2012</p>
|
||
</li>
|
||
|
||
<li id="fn:32" class="citation"><span class="citekey" style="display:none">diamanti_2015</span><p>Olga Diamanti, Amir Vaxman, Daniele Panozzo, Olga
|
||
Sorkine-Hornung. <a href="http://igl.ethz.ch/projects/integrable/">Integrable PolyVector Fields</a>, 2015</p>
|
||
</li>
|
||
|
||
<li id="fn:33" class="citation"><span class="citekey" style="display:none">vaxman_2016</span><p>Amir Vaxman, Marcel Campen, Olga Diamanti, Daniele Panozzo,
|
||
David Bommes, Klaus Hildebrandt, Mirela Ben-Chen. <a href="https://www.google.com/search?q=Directional+Field+Synthesis+Design+and+Processing">Directional Field
|
||
Synthesis, Design, and
|
||
Processing</a>,
|
||
2016</p>
|
||
</li>
|
||
|
||
<li id="fn:34" class="citation"><span class="citekey" style="display:none">schuller_2013</span><p>Christian Schüller, Ladislav Kavan, Daniele Panozzo, Olga
|
||
Sorkine-Hornung. <a href="http://igl.ethz.ch/projects/LIM/">Locally Injective
|
||
Mappings</a>, 2013.</p>
|
||
</li>
|
||
|
||
<li id="fn:35" class="citation"><span class="citekey" style="display:none">zhou_2016</span><p>Qingnan Zhou, Eitan Grinspun, Denis Zorin. <a href="https://www.google.com/search?q=Mesh+Arrangements+for+Solid+Geometry">Mesh Arrangements for
|
||
Solid
|
||
Geometry</a>,
|
||
2016</p>
|
||
</li>
|
||
|
||
<li id="fn:36" class="citation"><span class="citekey" style="display:none">jacobson_2013</span><p>Alec Jacobson, Ladislav Kavan, and Olga Sorkine.
|
||
<a href="https://www.google.com/search?q=Robust+Inside-Outside+Segmentation+using+Generalized+Winding+Numbers">Robust Inside-Outside Segmentation using Generalized Winding
|
||
Numbers</a>,
|
||
2013.</p>
|
||
</li>
|
||
|
||
<li id="fn:37" class="citation"><span class="citekey" style="display:none">hoppe_1996</span><p>Hugues Hoppe. <a href="https://www.google.com/search?q=Progressive+meshes">Progressive
|
||
Meshes</a>, 1996</p>
|
||
</li>
|
||
|
||
<li id="fn:38" class="citation"><span class="citekey" style="display:none">baerentzen_2005</span><p>J Andreas Baerentzen and Henrik Aanaes.
|
||
<a href="https://www.google.com/search?q=Signed+distance+computation+using+the+angle+weighted+pseudonormal">Signed distance computation using the angle weighted
|
||
pseudonormal</a>,
|
||
2005.</p>
|
||
</li>
|
||
|
||
<li id="fn:39" class="citation"><span class="citekey" style="display:none">lorensen_1987</span><p>W.E. Lorensen and Harvey E. Cline. <a href="https://www.google.com/search?q=Marching+cubes:+A+high+resolution+3d+surface+construction+algorithm">Marching cubes: A high
|
||
resolution 3d surface construction
|
||
algorithm</a>,
|
||
1987.</p>
|
||
</li>
|
||
|
||
<li id="fn:40" class="citation"><span class="citekey" style="display:none">takayama14</span><p>Kenshi Takayama, Alec Jacobson, Ladislav Kavan, Olga
|
||
Sorkine-Hornung. <a href="https://www.google.com/search?q=A+Simple+Method+for+Correcting+Facet+Orientations+in+Polygon+Meshes+Based+on+Ray+Casting">A Simple Method for Correcting Facet Orientations in
|
||
Polygon Meshes Based on Ray
|
||
Casting</a>,
|
||
2014.</p>
|
||
</li>
|
||
|
||
<li id="fn:41" class="citation"><span class="citekey" style="display:none">schroeder_1994</span><p>William J. Schroeder, William E. Lorensen, and Steve
|
||
Linthicum. <a href="https://www.google.com/search?q=implicit+modeling+of+swept+surfaces+and+volumes">Implicit Modeling of Swept Surfaces and
|
||
Volumes</a>,
|
||
1994.</p>
|
||
</li>
|
||
|
||
<li id="fn:42" class="citation"><span class="citekey" style="display:none">garg_2016</span><p>Akash Garg, Alec Jacobson, Eitan Grinspun. <a href="https://www.google.com/search?q=Computational+Design+of+Reconfigurables">Computational Design
|
||
of
|
||
Reconfigurables</a>,
|
||
2016</p>
|
||
</li>
|
||
|
||
<li id="fn:43" class="citation"><span class="citekey" style="display:none">rabinovich_2016</span><p>Michael Rabinovich, Roi Poranne, Daniele Panozzo, Olga
|
||
Sorkine-Hornung. <a href="http://cs.nyu.edu/~panozzo/papers/SLIM-2016.pdf">Scalable Locally Injective
|
||
Mappings</a>, 2016.</p>
|
||
</li>
|
||
|
||
<li id="fn:44" class="citation"><span class="citekey" style="display:none">loop_1987</span><p>Charles Loop. <a href="https://www.google.com/search?q=smooth+subdivision+surfaces+based+on+triangles">Smooth Subdivision Surfaces Based on
|
||
Triangles</a>,
|
||
1987.</p>
|
||
</li>
|
||
|
||
</ol>
|
||
</div>
|
||
|
||
|
||
</body>
|
||
</html>
|