Browse Source

laplacian example, viewer set_vertices better names better centroid

Former-commit-id: 2ec6b22960444b8581f6a064802655b590365ff1
Alec Jacobson 11 years ago
parent
commit
a2ced30016

+ 37 - 24
include/igl/viewer/Viewer.cpp

@@ -162,17 +162,18 @@ Eigen::Matrix4f translate(
   #include "igl/xml/XMLSerializer.h"
 #endif
 
-#include "igl/readOBJ.h"
-#include "igl/readOFF.h"
-#include "igl/per_face_normals.h"
-#include "igl/per_vertex_normals.h"
-#include "igl/per_corner_normals.h"
-#include "igl/vf.h"
-#include "igl/adjacency_list.h"
-#include "igl/writeOBJ.h"
-#include "igl/writeOFF.h"
-#include "igl/file_dialog_open.h"
-#include "igl/file_dialog_save.h"
+#include <igl/readOBJ.h>
+#include <igl/readOFF.h>
+#include <igl/per_face_normals.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/per_corner_normals.h>
+#include <igl/vf.h>
+#include <igl/adjacency_list.h>
+#include <igl/writeOBJ.h>
+#include <igl/writeOFF.h>
+#include <igl/massmatrix.h>
+#include <igl/file_dialog_open.h>
+#include <igl/file_dialog_save.h>
 #include <igl/quat_to_mat.h>
 #include <igl/quat_mult.h>
 #include <igl/axis_angle_to_quat.h>
@@ -777,7 +778,7 @@ namespace igl
     if (data.V_uv.rows() == 0)
       grid_texture();
 
-    alignCameraCenter();
+    align_camera_center();
 
     if (plugin_manager)
       for (unsigned int i = 0; i<plugin_manager->plugin_list.size(); ++i)
@@ -1991,9 +1992,9 @@ namespace igl
     return true;
   }
 
-  void Viewer::alignCameraCenter()
+  void Viewer::align_camera_center()
   {
-    get_scale_and_shift_to_fit_mesh(data.V,options.model_zoom,options.model_translation);
+    get_scale_and_shift_to_fit_mesh(data.V,data.F,options.model_zoom,options.model_translation);
     data.object_scale = (data.V.colwise().maxCoeff() - data.V.colwise().minCoeff()).norm();
   }
 
@@ -2004,16 +2005,21 @@ namespace igl
     viewport = Eigen::Vector4f(0,0,width,height);
   }
 
-  void Viewer::get_scale_and_shift_to_fit_mesh(const Eigen::MatrixXd& vertices,
-                                               float& zoom,
-                                               Eigen::Vector3f& shift)
+  void Viewer::get_scale_and_shift_to_fit_mesh(
+    const Eigen::MatrixXd& V,
+    const Eigen::MatrixXi& F,
+    float& zoom,
+    Eigen::Vector3f& shift)
   {
-    if (vertices.rows() == 0)
+    if (V.rows() == 0)
       return;
     //Compute mesh centroid
-    Eigen::RowVector3d centroid  = vertices.colwise().sum()/vertices.rows();
-    Eigen::RowVector3d min_point = vertices.colwise().minCoeff();
-    Eigen::RowVector3d max_point = vertices.colwise().maxCoeff();
+    Eigen::SparseMatrix<double> M;
+    igl::massmatrix(V,F,igl::MASSMATRIX_VORONOI,M);
+    const auto & MV = M*V;
+    Eigen::RowVector3d centroid  = MV.colwise().sum()/M.diagonal().sum();
+    Eigen::RowVector3d min_point = V.colwise().minCoeff();
+    Eigen::RowVector3d max_point = V.colwise().maxCoeff();
 
     shift = -centroid.cast<float>();
     double x_scale = fabs(max_point[0] - min_point[0]);
@@ -2029,7 +2035,7 @@ namespace igl
   }
   void TW_CALL Viewer::align_camera_center_cb(void *clientData)
   {
-    static_cast<Viewer *>(clientData)->alignCameraCenter();
+    static_cast<Viewer *>(clientData)->align_camera_center();
   }
 
   void TW_CALL Viewer::save_scene_cb(void *clientData)
@@ -2176,7 +2182,7 @@ namespace igl
                      Eigen::Vector3d(255.0/255.0,235.0/255.0,80.0/255.0));
 
       grid_texture();
-      alignCameraCenter();
+      align_camera_center();
     }
     else
     {
@@ -2184,7 +2190,7 @@ namespace igl
       {
         data.V = V_temp;
         data.F = F;
-        alignCameraCenter();
+        align_camera_center();
       }
       else
         cerr << "ERROR (set_mesh): The new mesh has a different number of vertices/faces. Please clear the mesh before plotting.";
@@ -2192,6 +2198,13 @@ namespace igl
     data.dirty |= DIRTY_FACE | DIRTY_POSITION;
   }
 
+  void Viewer::set_vertices(const Eigen::MatrixXd& V)
+  {
+    data.V = V;
+    assert(F.empty() || F.maxCoeff() < V.rows());
+    data.dirty |= DIRTY_POSITION;
+  }
+
   void Viewer::set_normals(const Eigen::MatrixXd& N)
   {
     using namespace std;

+ 7 - 4
include/igl/viewer/Viewer.h

@@ -344,13 +344,14 @@ namespace igl
     void grid_texture(); // Generate a default grid texture
 
     void clear_mesh();      // Clear the mesh data
-    void alignCameraCenter(); // Adjust the view to see the entire model
+    void align_camera_center(); // Adjust the view to see the entire model
 
     // Change the visualization mode, invalidating the cache if necessary
     void set_face_based(bool newvalue);
 
     // Helpers that can draw the most common meshes
     void set_mesh(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F);
+    void set_vertices(const Eigen::MatrixXd& V);
     void set_normals(const Eigen::MatrixXd& N);
     // Set the color of the mesh
     //
@@ -403,9 +404,11 @@ namespace igl
 
     // Determines how much to zoom and shift such that the mesh fills the unit
     // box (centered at the origin)
-    static void get_scale_and_shift_to_fit_mesh(const Eigen::MatrixXd& vertices,
-                                                float & zoom,
-                                                Eigen::Vector3f& shift);
+    static void get_scale_and_shift_to_fit_mesh(
+      const Eigen::MatrixXd& V,
+      const Eigen::MatrixXi& F,
+      float & zoom,
+      Eigen::Vector3f& shift);
 
 
     // Init opengl shaders and VBOs

+ 11 - 0
tutorial/205_Laplacian/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(205_Laplacian)
+
+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})

+ 74 - 0
tutorial/205_Laplacian/main.cpp

@@ -0,0 +1,74 @@
+#include <igl/readOFF.h>
+#include <igl/readDMAT.h>
+#include <igl/cotmatrix.h>
+#include <igl/massmatrix.h>
+#include <igl/jet.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/viewer/Viewer.h>
+
+#include <iostream>
+
+Eigen::MatrixXd V,U;
+Eigen::MatrixXi F;
+Eigen::SparseMatrix<double> L;
+igl::Viewer viewer;
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+
+  // Load a mesh in OFF format
+  igl::readOFF("../shared/cow.off", V, F);
+
+  // Compute Laplace-Beltrami operator: #V by #V
+  igl::cotmatrix(V,F,L);
+
+  const auto &key_down = [](igl::Viewer &viewer,unsigned char key,int mod)->bool
+  {
+    switch(key)
+    {
+      case 'r':
+      case 'R':
+        U = V;
+        break;
+      case ' ':
+      {
+        // Recompute just mass matrix on each step
+        SparseMatrix<double> M;
+        igl::massmatrix(U,F,igl::MASSMATRIX_VORONOI,M);
+        // Solve (M-delta*L) U = M*U
+        const auto & S = (M - 0.001*L);
+        Eigen::SimplicialLLT<Eigen::SparseMatrix<double > > solver(S);
+        assert(solver.info() == Eigen::Success);
+        U = solver.solve(M*U).eval();
+        // Normalize to unit surface area (important for numerics)
+        U.array() /= sqrt(M.diagonal().array().sum());
+        break;
+      }
+      default:
+        return false;
+    }
+    // Send new positions, update normals, recenter
+    viewer.set_vertices(U);
+    viewer.compute_normals();
+    viewer.align_camera_center();
+    return true;
+  };
+
+
+  // Use original normals as pseudo-colors
+  MatrixXd N;
+  igl::per_vertex_normals(V,F,N);
+  MatrixXd C = N.rowwise().normalized().array()*0.5+0.5;
+
+  // Initialize smoothing with base mesh
+  U = V;
+  viewer.set_mesh(U, F);
+  viewer.set_colors(C);
+  viewer.callback_key_down = key_down;
+
+  cout<<"Press [space] to smooth."<<endl;;
+  cout<<"Press [r] to reset."<<endl;;
+  viewer.launch();
+}

+ 1 - 0
tutorial/images/cow-curvature-flow.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+e95f80c272b8beed26cff70d5df644e8441150ce

+ 98 - 0
tutorial/readme.md

@@ -36,6 +36,7 @@ applications for each topic.
 * [202 Gaussian Curvature](#gaus)
 * [203 Curvature Directions](#curv)
 * [204 Gradient](#grad)
+* [204 Laplacian](#lapl)
 
 # Compilation Instructions
 
@@ -192,6 +193,101 @@ 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)
 
+## <a id=lapl></a> Laplacian
+
+The discrete Laplacian is an essential geometry processing tool. Many
+interpretations and flavors of the Laplace and Laplace-Beltrami operator exist. 
+
+In open Euclidean space, the _Laplace_ operator is the usual divergence of gradient
+(or equivalently the Laplacian of a function is the trace of its Hessian):
+
+ $\Delta f = 
+ \frac{\partial^2 f}{\partial x^2} +
+ \frac{\partial^2 f}{\partial y^2} +
+ \frac{\partial^2 f}{\partial z^2}.$
+
+The _Laplace-Beltrami_ operator generalizes this to surfaces.
+
+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'' $\mathbf{L}$, 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 $\mathbf{L}$ is a $n\times n$
+matrix with elements:
+
+$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}$
+
+where $N(i)$ are the vertices adjacent to (neighboring) vertex $i$, and
+$\alpha_{ij},\beta_{ij}$ are the angles opposite edge ${ij}$.
+This oft
+produced formula leads to a typical half-edge style implementation for
+constructing $\mathbf{L}$:
+
+```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));
+    }
+  }
+}
+```
+
+Without a half-edge data-structure it may seem at first glance that looping
+over one-rings, and thus constructing the Laplacian would be inefficient.
+However, 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):
+
+```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));
+  }
+}
+```
+
+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.
+**Alec: cite Ariel reconstruction paper**.
+
+The operator applied to mesh vertex positions amounts to smoothing by _flowing_
+the surface along the mean curvature normal direction. This is equivalent to
+minimizing surface area.
+
+![The `Laplacian` example computes conformalized mean curvature flow using the
+cotangent Laplacian [#kazhdan_2012][].](images/cow-curvature-flow.jpg)
+
+### Mass matrix
+The mass matrix $\mathbf{M}$ is another $n \times n$ 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,
+$\mathbf{M}$ is often a diagonal matrix, such that $M_{ii}$ is the barycentric
+or voronoi area around vertex $i$ in the mesh [#meyer_2003][]. The inverse of
+this matrix is also very useful as it transforms integrated quantities into
+point-wise quantities, e.g.:
+
+ $\nabla f \approx \mathbf{M}^{-1} \mathbf{L} \mathbf{f}.$
+
+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:
+
+ $\int_S x\, y\ dA \approx \mathbf{x}^T\mathbf{M}\,\mathbf{y}.$
+
 [#meyer_2003]: Mark Meyer and Mathieu Desbrun and Peter Schröder and Alan H.  Barr,
 "Discrete Differential-Geometry Operators for Triangulated
 2-Manifolds," 2003.
@@ -200,3 +296,5 @@ visualizes the vector field.](images/cheburashka-gradient.jpg)
 [#jacobson_thesis_2013]: Alec Jacobson,
 _Algorithms and Interfaces for Real-Time Deformation of 2D and 3D Shapes_,
 2013.
+[#kazhdan_2012]: Michael Kazhdan, Jake Solomon, Mirela Ben-Chen,
+"Can Mean-Curvature Flow Be Made Non-Singular," 2012.

+ 1 - 0
tutorial/shared/cow.off.REMOVED.git-id

@@ -0,0 +1 @@
+3751246241296df73a713893cb9b873d2c6bc444