Browse Source

Gradient example

Former-commit-id: ded5fc4186b3a48f01c8a779840f8abf6209f3d1
Alec Jacobson 11 years ago
parent
commit
cac3b04a2a

+ 3 - 1
include/igl/gradMat.cpp

@@ -6,6 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "gradMat.h"
+#include <Eigen/Geometry>
 #include <vector>
 
 template <typename T, typename S>
@@ -93,7 +94,7 @@ IGL_INLINE void igl::gradMat(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynam
   // create sparse gradient operator matrix
   G.resize(3*F.rows(),V.rows()); 
   std::vector<Eigen::Triplet<T> > triplets;
-  for (int i=0;i<vs.size();++i)
+  for (int i=0;i<(int)vs.size();++i)
   {
     triplets.push_back(Eigen::Triplet<T>(rs[i],cs[i],vs[i]));
   }
@@ -102,5 +103,6 @@ IGL_INLINE void igl::gradMat(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynam
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
+template void igl::gradMat<double, int>(Eigen::Matrix<double, -1, -1, 0, -1,-1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&,Eigen::SparseMatrix<double, 0, int>&);
 #endif
 

+ 1 - 0
include/igl/jet.cpp

@@ -122,4 +122,5 @@ template void igl::jet<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<dou
 template void igl::jet<double>(double, double*);
 template void igl::jet<double>(double, double&, double&, double&);
 template void igl::jet<float>(float, float*);
+template void igl::jet<Eigen::Array<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Array<double, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

+ 1 - 1
tutorial/201_Normals/CMakeLists.txt

@@ -1,5 +1,5 @@
 cmake_minimum_required(VERSION 2.6)
-project(202_Normals)
+project(201_Normals)
 
 include("../CMakeLists.shared")
 

+ 2 - 3
tutorial/202_GaussianCurvature/main.cpp

@@ -6,18 +6,17 @@
 int main(int argc, char *argv[])
 {
   using namespace Eigen;
-  using namespace igl;
   using namespace std;
   MatrixXd V;
   MatrixXi F;
   igl::readOFF("../shared/bumpy.off",V,F);
 
   VectorXd K;
-  gaussian_curvature(V,F,K);
+  igl::gaussian_curvature(V,F,K);
 
   // Compute pseudocolor
   MatrixXd C;
-  jet(K,true,C);
+  igl::jet(K,true,C);
 
   // Plot the mesh with pseudocolors
   igl::Viewer viewer;

+ 3 - 2
tutorial/203_CurvatureDirections/main.cpp

@@ -4,7 +4,6 @@
 #include <igl/per_face_normals.h>
 #include <igl/per_corner_normals.h>
 #include <igl/avg_edge_length.h>
-#include <igl/barycenter.h>
 #include <igl/principal_curvature.h>
 #include <igl/jet.h>
 #include <igl/cotmatrix.h>
@@ -26,8 +25,10 @@ int main(int argc, char *argv[])
   igl::cotmatrix(V,F,L);
   igl::massmatrix(V,F,igl::MASSMATRIX_VORONOI,M);
   igl::invert_diag(M,Minv);
+  // Laplace-Beltrami of position
   HN = -Minv*(L*V);
-  H = (HN.rowwise().squaredNorm()).array().sqrt();
+  // Extract magnitude as mean curvature
+  H = HN.rowwise().norm();
 
   // Compute curvature directions via quadric fitting
   MatrixXd PD1,PD2;

+ 11 - 0
tutorial/204_Gradient/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(204_Gradient)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${CMAKE_PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${CMAKE_PROJECT_NAME} ${SHARED_LIBRARIES})

+ 56 - 0
tutorial/204_Gradient/main.cpp

@@ -0,0 +1,56 @@
+#include <igl/readOFF.h>
+#include <igl/readDMAT.h>
+#include <igl/gradMat.h>
+#include <igl/avg_edge_length.h>
+#include <igl/jet.h>
+#include <igl/barycenter.h>
+#include <igl/viewer/Viewer.h>
+
+#include <iostream>
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd V;
+  MatrixXi F;
+
+  // Load a mesh in OFF format
+  igl::readOFF("../shared/cheburashka.off", V, F);
+
+  // Read scalar function values from a file, U: #V by 1
+  VectorXd U;
+  igl::readDMAT("../shared/cheburashka-scalar.dmat",U);
+
+  // Compute gradient operator: #F*3 by #V
+  SparseMatrix<double> G;
+  igl::gradMat(V,F,G);
+
+  // Compute gradient of U
+  MatrixXd GU = Map<const MatrixXd>((G*U).eval().data(),F.rows(),3);
+  // Compute gradient magnitude
+  const VectorXd GU_mag = GU.rowwise().norm();
+
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+
+  // Compute pseudocolor for original function
+  MatrixXd C;
+  igl::jet(U,true,C);
+  // // Or for gradient magnitude
+  //igl::jet(GU_mag,true,C);
+  viewer.set_colors(C);
+
+  // Average edge length divided by average gradient (for scaling)
+  const double max_size = igl::avg_edge_length(V,F) / GU_mag.mean();
+  // Draw a black segment in direction of gradient at face barycenters
+  MatrixXd BC;
+  igl::barycenter(V,F,BC);
+  const RowVector3d black(0,0,0);
+  viewer.add_edges(BC,BC+max_size*GU, black);
+
+  // Hide wireframe
+  viewer.options.show_lines = false;
+
+  viewer.launch();
+}

BIN
tutorial/images/background.gif


+ 1 - 0
tutorial/images/cheburashka-gradient.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+f423446740ec2dd88082ee2a47d836d378e07495

+ 1 - 0
tutorial/images/fertility-principal-curvature.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+54d3ee53cc0b43659f330f645a76a79f3b3ab851

BIN
tutorial/images/hat-function.jpg


+ 74 - 7
tutorial/readme.md

@@ -1,3 +1,6 @@
+title: libigl Tutorial
+author: Alec Jacobson, Daniele Pannozo and others
+date: 20 June 2014
 css: style.css
 html header:   <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">
@@ -5,8 +8,25 @@ html header:   <script type="text/javascript" src="http://cdn.mathjax.org/mathja
 <script>hljs.initHighlightingOnLoad();</script>
 
 # Introduction
-
-TODO
+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 \textsc{Matlab} or
+\textsc{Python} programmers
+with the performance and versatility of C++.
+%
+The tutorial is a self-contained, hands-on introduction to \libigl. 
+%
+Via live coding and interactive examples, we demonstrate how to accomplish
+various common geometry processing tasks such as computation of differential
+quantities and operators, real-time deformation, global parametrization,
+numerical optimization and mesh repair.
+%
+Accompanying lecture notes contain further details and cross-platform example
+applications for each topic.
+\end{abstract}
 
 # Index
 
@@ -14,6 +34,8 @@ TODO
 * **101_Serialization**: Example of using the XML serialization framework
 * **102_DrawMesh**: Example of plotting a mesh
 * [202 Gaussian Curvature](#gaus)
+* [203 Curvature Directions](#curv)
+* [204 Gradient](#grad)
 
 # Compilation Instructions
 
@@ -68,7 +90,7 @@ on a triangle mesh is via a vertex's _angular deficit_:
  $k_G(v_i) = 2π - \sum\limits_{j\in N(i)}θ_{ij},$
 
 where $N(i)$ are the triangles incident on vertex $i$ and $θ_{ij}$ is the angle
-at vertex $i$ in triangle $j$. (**Alec: cite Meyer or something**)
+at vertex $i$ in triangle $j$ [][#meyer_2003].
 
 Just like the continuous analog, our discrete Gaussian curvature reveals
 elliptic, hyperbolic and parabolic vertices on the domain.
@@ -76,7 +98,7 @@ elliptic, hyperbolic and parabolic vertices on the domain.
 ![The `GaussianCurvature` example computes discrete Gaussian curvature and visualizes it in
 pseudocolor.](images/bumpy-gaussian-curvature.jpg)
 
-## <a id=principal></a> Curvature Directions
+## <a id=curv></a> Curvature Directions
 The two principal curvatures $(k_1,k_2)$ at a point on a surface measure how much the
 surface bends in different directions. The directions of maximum and minimum
 (signed) bending are call principal directions and are always
@@ -93,7 +115,7 @@ normal:
   $-\Delta \mathbf{x} = H \mathbf{n}.$ 
 
 It is easy to compute this on a discrete triangle mesh in libigl using the cotangent
-Laplace-Beltrami operator (**Alec: cite Meyer**):
+Laplace-Beltrami operator [][#meyer_2003].
 
 ```cpp
 #include <igl/cotmatrix.h>
@@ -111,10 +133,10 @@ H = (HN.rowwise().squaredNorm()).array().sqrt();
 
 Combined with the angle defect definition of discrete Gaussian curvature, one
 can define principal curvatures and use least squares fitting to find
-directions (**Alec: cite meyer**).
+directions [][#meyer_2003].
 
 Alternatively, a robust method for determining principal curvatures is via
-quadric fitting (**Alec: cite whatever we're using**). In the neighborhood
+quadric fitting [][#pannozo_2010]. In the neighborhood
 around every vertex, a best-fit quadric is found and principal curvature values
 and directions are sampled from this quadric. With these in tow, one can
 compute mean curvature and Gaussian curvature as sums and products
@@ -133,3 +155,48 @@ int main(int argc, char * argv[])
   return 0;
 }
 ```
+
+## <a id=grad></a> Gradient
+Scalar functions on a surface can be discretized as a piecewise linear function
+with values defined at each mesh vertex:
+
+ $f(\mathbf{x}) \approx \sum\limits_{i=0}^n \phi_i(\mathbf{x})\, f_i,$
+
+where $\phi_i$ is a piecewise linear hat function defined by the mesh so that
+for each triangle $\phi_i$ is _the_ linear function which is one only at
+vertex $i$ and zero at the other corners.
+
+![Hat function $\phi_i$ is one at vertex $i$, zero at all other vertices, and
+linear on incident triangles.](images/hat-function.jpg)
+
+Thus gradients of such piecewise linear functions are simply sums of gradients
+of the hat functions:
+
+ $\nabla f(\mathbf{x}) \approx 
+ \nabla \sum\limits_{i=0}^n \nabla \phi_i(\mathbf{x})\, f_i = 
+ \sum\limits_{i=0}^n \nabla \phi_i(\mathbf{x})\, f_i.$
+
+This reveals that the gradient is a linear function of the vector of $f_i$
+values. Because $\phi_i$ are linear in each triangle their gradient are
+_constant_ in each triangle. Thus our discrete gradient operator can be written
+as a matrix multiplication taking vertex values to triangle values:
+
+ $\nabla f \approx \mathbf{G}\,\mathbf{f},$
+
+where $\mathbf{f}$ is $n\times 1$ and $\mathbf{G}$ is an $m\times n$ sparse
+matrix. This matrix $\mathbf{G}$ can be derived geometrically, e.g.
+[ch. 2][#jacobson_thesis_2013].
+Libigl's `gradMat`**Alec: check name** function computes $\mathbf{G}$ for
+triangle and tetrahedral meshes:
+
+![The `Gradient` example computes gradients of an input function on a mesh and
+visualizes the vector field.](images/cheburashka-gradient.jpg)
+
+[#meyer_2003]: Mark Meyer and Mathieu Desbrun and Peter Schröder and Alan H.  Barr,
+"Discrete Differential-Geometry Operators for Triangulated
+2-Manifolds," 2003.
+[#pannozo_2010]: Daniele Pannozo, Enrico Puppo, Luigi Rocca,
+"Efficient Multi-scale Curvature and Crease Estimation," 2010.
+[#jacobson_thesis_2013]: Alec Jacobson,
+_Algorithms and Interfaces for Real-Time Deformation of 2D and 3D Shapes_,
+2013.

+ 1 - 0
tutorial/shared/cheburashka-scalar.dmat.REMOVED.git-id

@@ -0,0 +1 @@
+217ce358634f5d12ec16b1e9965d9b1f9ec0cdfa

+ 10 - 4
tutorial/style.css

@@ -40,7 +40,13 @@ body
 
 img
 {
-  max-width: 650px;
+  max-width: 600px;
+}
+
+figure
+{
+  margin: 0px auto;
+  width: 600px;
 }
 
 figcaption
@@ -159,16 +165,16 @@ h1 {
 }
 
 h2 {
-  font-size: 1.4em;
+  font-size: 1.3em;
   border-bottom: 1px solid #eee;
 }
 
 h3 {
-  font-size: 1.3em;
+  font-size: 1.2em;
 }
 
 h4 {
-  font-size: 1.2em;
+  font-size: 1.1em;
 }
 
 h5 {