floor/ceil, slice and sort exampels

Former-commit-id: 6b245d62eab3dbc0f3e7652fe1ef75a2b89a36b5

include/igl/ceil.cpp

@@ -0,0 +1,33 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// 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 "ceil.h"
+#include <cmath>
+
+template < typename DerivedX, typename DerivedY>
+IGL_INLINE void igl::ceil(
+  const Eigen::PlainObjectBase<DerivedX>& X,
+  Eigen::PlainObjectBase<DerivedY>& Y)
+{
+  using namespace std;
+  //Y = Eigen::PlainObjectBase<DerivedY>::Zero(m,n);
+//#pragma omp parallel for
+  //for(int i = 0;i<m;i++)
+  //{
+  //  for(int j = 0;j<n;j++)
+  //  {
+  //    Y(i,j) = std::ceil(X(i,j));
+  //  }
+  //}
+  typedef typename DerivedX::Scalar Scalar;
+  Y = X.unaryExpr([](const Scalar &x)->Scalar{return std::ceil(x);}).template cast<typename DerivedY::Scalar >();
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit instanciation
+template void igl::ceil<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/ceil.h

@@ -0,0 +1,30 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// 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/.
+#ifndef IGL_CEIL_H
+#define IGL_CEIL_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+namespace igl
+{
+  // Ceil a given matrix to nearest integers 
+  //
+  // Inputs:
+  //   X  m by n matrix of scalars
+  // Outputs:
+  //   Y  m by n matrix of ceiled integers
+  template < typename DerivedX, typename DerivedY>
+  IGL_INLINE void ceil(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    Eigen::PlainObjectBase<DerivedY>& Y);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "ceil.cpp"
+#endif
+
+#endif

include/igl/diag.h

@@ -16,8 +16,9 @@ namespace igl
   // http://forum.kde.org/viewtopic.php?f=74&t=117476&p=292388#p292388
   //
   // This is superceded by 
-  //   VectorXd d = X.diagonal() and 
-  //   SparseVector<double> d = X.diagonal().sparseView()
+  //   VectorXd V = X.diagonal() and 
+  //   SparseVector<double> V = X.diagonal().sparseView()
+  //   SparseMatrix<double> X = V.asDiagonal().sparseView()
   //
   //
   // Either extracts the main diagonal of a matrix as a vector OR converts a

include/igl/floor.cpp

@@ -0,0 +1,35 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// 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 "floor.h"
+#include <cmath>
+#include <iostream>
+
+template < typename DerivedX, typename DerivedY>
+IGL_INLINE void igl::floor(
+  const Eigen::PlainObjectBase<DerivedX>& X,
+  Eigen::PlainObjectBase<DerivedY>& Y)
+{
+  using namespace std;
+  //Y = Eigen::PlainObjectBase<DerivedY>::Zero(m,n);
+//#pragma omp parallel for
+  //for(int i = 0;i<m;i++)
+  //{
+  //  for(int j = 0;j<n;j++)
+  //  {
+  //    Y(i,j) = std::floor(X(i,j));
+  //  }
+  //}
+  typedef typename DerivedX::Scalar Scalar;
+  Y = X.unaryExpr([](const Scalar &x)->Scalar{return std::floor(x);}).template cast<typename DerivedY::Scalar >();
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit instanciation
+template void igl::floor<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::floor<Eigen::Array<double, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Array<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/floor.h

@@ -0,0 +1,30 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// 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/.
+#ifndef IGL_FLOOR_H
+#define IGL_FLOOR_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+namespace igl
+{
+  // Floor a given matrix to nearest integers 
+  //
+  // Inputs:
+  //   X  m by n matrix of scalars
+  // Outputs:
+  //   Y  m by n matrix of floored integers
+  template < typename DerivedX, typename DerivedY>
+  IGL_INLINE void floor(
+    const Eigen::PlainObjectBase<DerivedX>& X,
+    Eigen::PlainObjectBase<DerivedY>& Y);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "floor.cpp"
+#endif
+
+#endif

include/igl/jet.cpp

@@ -111,7 +111,7 @@ IGL_INLINE void igl::jet(
   const double max_z = (normalize?Z.maxCoeff():-1);
   for(int r = 0;r<Z.rows();r++)
   {
-    jet((Z(r,0)-min_z)/(max_z-min_z),C(r,0),C(r,1),C(r,2));
+    jet((-min_z+Z(r,0))/(max_z-min_z),C(r,0),C(r,1),C(r,2));
   }
 }
 
@@ -123,4 +123,5 @@ 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> >&);
+template void igl::jet<Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/reorder.h

@@ -15,7 +15,7 @@
 
 namespace igl
 {
-  // Act like matlab's Y = X[I] for std vectors
+  // Act like matlab's Y = X(I) for std vectors
   // where I contains a vector of indices so that after,
   // Y[j] = X[I[j]] for index j
   // this implies that Y.size() == I.size()

include/igl/slice_into.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 "slice_into.h"
+#include "colon.h"
 
 // Bug in unsupported/Eigen/SparseExtra needs iostream first
 #include <iostream>
@@ -79,6 +80,40 @@ IGL_INLINE void igl::slice_into(
   }
 }
 
+template <typename Mat>
+IGL_INLINE void igl::slice_into(
+  const Mat& X,
+  const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+  const int dim,
+  Mat& Y)
+{
+  Eigen::VectorXi C;
+  switch(dim)
+  {
+    case 1:
+      assert(R.size() == X.rows());
+      // boring base case
+      if(X.cols() == 0)
+      {
+        return;
+      }
+      igl::colon(0,X.cols()-1,C);
+      return slice_into(X,R,C,Y);
+    case 2:
+      assert(R.size() == X.cols());
+      // boring base case
+      if(X.rows() == 0)
+      {
+        return;
+      }
+      igl::colon(0,X.rows()-1,C);
+      return slice_into(X,C,R,Y);
+    default:
+      assert(false && "Unsupported dimension");
+      return;
+  }
+}
+
 template <typename DerivedX>
 IGL_INLINE void igl::slice_into(
   const Eigen::PlainObjectBase<DerivedX> & X,
@@ -97,4 +132,6 @@ IGL_INLINE void igl::slice_into(
 // generated by autoexplicit.sh
 template void igl::slice_into<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::slice_into<Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::slice_into<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, int, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::slice_into<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/slice_into.h

@@ -36,6 +36,18 @@ namespace igl
     const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
     const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
     Eigen::PlainObjectBase<DerivedX> & Y);
+  // Wrapper to only slice in one direction
+  //
+  // Inputs:
+  //   dim  dimension to slice in 1 or 2, dim=1 --> X(R,:), dim=2 --> X(:,R)
+  //
+  // Note: For now this is just a cheap wrapper.
+  template <typename Mat>
+  IGL_INLINE void slice_into(
+    const Mat& X,
+    const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+    const int dim,
+    Mat& Y);
 
   template <typename DerivedX>
   IGL_INLINE void slice_into(

include/igl/viewer/Viewer.cpp.REMOVED.git-id

@@ -1 +1 @@
-5381ba4d384ee8a7bda552c5bf95a26335b93b48
+be827c7827ae18c27874d6e620ad13f84f9b37de

include/igl/viewer/Viewer.h

@@ -31,7 +31,7 @@ namespace igl
   {
   public:
 
-    void launch(std::string filename = "");
+    int launch(std::string filename = "");
     void init(Plugin_manager* pm);
 
     class Options

matlab-to-eigen.html

@@ -11,7 +11,7 @@
       </tr>
 
       <tr class=d0>
-        <td><pre><code>[Y,IX] = sort(Y,dim,mode)</code></pre></td>
+        <td><pre><code>[Y,IX] = sort(X,dim,mode)</code></pre></td>
         <td><pre><code>igl::sort(X,dim,mode,Y,IX)</code></pre></td>
         <td>MATLAB version allows Y to be a multidimensional matrix, but the
         Eigen version is only for 1D or 2D matrices.</td>
@@ -45,7 +45,7 @@
 
       <tr class=d1>
         <td><pre><code>A(i:(i+w),j:(j+h)) = eye(w,h)</code></pre></td>
-        <td><pre><code>A.setIdentity()</code></pre></td>
+        <td><pre><code>A.block(i,j,w,h).setIdentity()</code></pre></td>
         <td></td>
       </tr>
 

tutorial/205_Laplacian/main.cpp

@@ -2,6 +2,9 @@
 #include <igl/readDMAT.h>
 #include <igl/cotmatrix.h>
 #include <igl/massmatrix.h>
+#include <igl/gradMat.h>
+#include <igl/doublearea.h>
+#include <igl/repdiag.h>
 #include <igl/jet.h>
 #include <igl/per_vertex_normals.h>
 #include <igl/viewer/Viewer.h>
@@ -24,6 +27,20 @@ int main(int argc, char *argv[])
   // Compute Laplace-Beltrami operator: #V by #V
   igl::cotmatrix(V,F,L);
 
+  // Alternative construction of same Laplacian
+  SparseMatrix<double> G,K;
+  // Gradient/Divergence
+  igl::gradMat(V,F,G);
+  // Diagonal per-triangle "mass matrix"
+  VectorXd dblA;
+  igl::doublearea(V,F,dblA);
+  // Place areas along diagonal #dim times
+  const auto & T = 1.*(dblA.replicate(3,1)*0.5).asDiagonal();
+  // Laplacian K built as discrete divergence of gradient or equivalently
+  // discrete Dirichelet energy Hessian
+  K = -G.transpose() * T * G;
+  cout<<"|K-L|: "<<(K-L).norm()<<endl;
+
   const auto &key_down = [](igl::Viewer &viewer,unsigned char key,int mod)->bool
   {
     switch(key)
@@ -70,5 +87,5 @@ int main(int argc, char *argv[])
 
   cout<<"Press [space] to smooth."<<endl;;
   cout<<"Press [r] to reset."<<endl;;
-  viewer.launch();
+  return viewer.launch();
 }

tutorial/301_Slice/CMakeLists.txt

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

tutorial/301_Slice/main.cpp

@@ -0,0 +1,40 @@
+#include <igl/readOFF.h>
+#include <igl/slice.h>
+#include <igl/floor.h>
+#include <igl/slice_into.h>
+#include <igl/viewer/Viewer.h>
+#include <iostream>
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd V;
+  MatrixXi F;
+  igl::readOFF("../shared/decimated-knight.off",V,F);
+
+  // 100 random indicies into rows of F
+  VectorXi I;
+  igl::floor((0.5*(VectorXd::Random(100,1).array()+1.)*F.rows()).eval(),I);
+  
+  // 50 random indicies into rows of I
+  VectorXi J;
+  igl::floor((0.5*(VectorXd::Random(50,1).array()+1.)*I.rows()).eval(),J);
+  
+  // K = I(J);
+  VectorXi K;
+  igl::slice(I,J,K);
+
+  // default green for all faces
+  MatrixXd C = RowVector3d(0.4,0.8,0.3).replicate(F.rows(),1);
+  // Red for each in K
+  MatrixXd R = RowVector3d(1.0,0.3,0.3).replicate(K.rows(),1);
+  // C(K,:) = R
+  igl::slice_into(R,K,1,C);
+
+  // Plot the mesh with pseudocolors
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+  viewer.set_colors(C);
+  viewer.launch();
+}

tutorial/302_Sort/CMakeLists.txt

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

tutorial/302_Sort/main.cpp

@@ -0,0 +1,37 @@
+#include <igl/readOFF.h>
+#include <igl/slice_into.h>
+#include <igl/colon.h>
+#include <igl/sortrows.h>
+#include <igl/barycenter.h>
+#include <igl/jet.h>
+#include <igl/viewer/Viewer.h>
+#include <iostream>
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd V;
+  MatrixXi F;
+  igl::readOFF("../shared/decimated-knight.off",V,F);
+
+  // Sort barycenters lexicographically
+  MatrixXd BC,sorted_BC;
+  igl::barycenter(V,F,BC);
+  VectorXi I,J;
+  igl::sortrows(BC,true,sorted_BC,I);
+  // Get sorted "place" from sorted indices
+  J.resize(I.rows());
+  // J(I) = 1:numel(I)
+  igl::slice_into(igl::colon<int>(0,I.size()-1),I,J);
+
+  // Pseudo-color based on sorted place
+  MatrixXd C;
+  igl::jet(J,true,C);
+
+  // Plot the mesh with pseudocolors
+  igl::Viewer viewer;
+  viewer.set_mesh(V, F);
+  viewer.set_colors(C);
+  viewer.launch();
+}

tutorial/images/decimated-knight-slice-color.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+480151754ff0e57c7111cc288fd6c0b6b163f404

tutorial/images/decimated-knight-sort-color.jpg.REMOVED.git-id

@@ -0,0 +1 @@
+7172cceb268b1fb836116ae527caa912200b416c

tutorial/readme.md

@@ -9,35 +9,39 @@ html header:   <script type="text/javascript" src="http://cdn.mathjax.org/mathja
 
 # Introduction
 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
-
-* **100_FileIO**: Example of reading/writing mesh files
-* **101_Serialization**: Example of using the XML serialization framework
-* **102_DrawMesh**: Example of plotting a mesh
-* [201 Normals](#norm)
-* [202 Gaussian Curvature](#gaus)
-* [203 Curvature Directions](#curv)
-* [204 Gradient](#grad)
-* [204 Laplacian](#lapl)
+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 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.  Each section
+of these lecture notes links to a cross-platform example application.
+
+# Table of Contents
+
+* Basic Usage
+    * **100_FileIO**: Example of reading/writing mesh files
+    * **101_Serialization**: Example of using the XML serialization framework
+    * **102_DrawMesh**: Example of plotting a mesh
+* [Chapter 2: Discrete Geometric Quantities and
+  Operators](#chapter2:discretegeometricquantitiesandoperators)
+    * [201 Normals](#normals)
+        * [Per-face](#per-face)
+        * [Per-vertex](#per-vertex)
+        * [Per-corner](#per-corner)
+    * [202 Gaussian Curvature](#gaussiancurvature)
+    * [203 Curvature Directions](#curvaturedirections)
+    * [204 Gradient](#gradient)
+    * [204 Laplacian](#laplacian)
+        * [Mass matrix](#massmatrix)
+        * [Alternative construction of
+          Laplacian](#alternativeconstuctionoflaplacian)
+* [Chapter 3: Matrices and Linear Algebra](#chapter3:matricesandlinearalgebra)
+    * [301 Slice](#slice)
+    * [301 Sort](#sort)
+        * [Other Matlab-style functions](#othermatlab-stylefunctions) 
 
 # Compilation Instructions
 
@@ -72,18 +76,23 @@ mesh. This also provides an introduction to basic drawing and coloring routines
 in our example viewer. Finally, we construct popular discrete differential
 geometry operators.
 
-## <a id=norm></a> Normals
+## Normals
+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.
 
-### Per-face 
-Well defined normal orthogonal to triangle's plane. Produces piecewise-flat
-rending: not smooth.
+### Per-face
+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.
 
 ### Per-vertex
-Using per-vertex normals, Phong or Gouraud shading will produce smooth(er)
-renderings. Most techniques for computing per-vertex normals take an average of
-incident face normals. The techniques vary with respect to their different
-weighting schemes. Uniform weighting is heavily biased by the discretization
-choice, where as area-based or angle-based weighting is more forgiving.
+Storing normals at vertices, Phong or Gouraud shading will interpolate shading
+inside mesh triangles to produce smooth(er) renderings. Most techniques for
+computing per-vertex normals take an average of incident face normals. The
+techniques vary with respect to their different weighting schemes. Uniform
+weighting is heavily biased by the discretization choice, where as area-based
+or angle-based weighting is more forgiving.
 
 The typical half-edge style computation of area-based weights might look
 something like this:
@@ -130,7 +139,7 @@ specified dihedral angle (e.g. 20°).
 ![The `Normals` example computes per-face (left), per-vertex (middle) and
 per-corner (right) normals](images/fandisk-normals.jpg)
 
-## <a id=gaus></a> Gaussian Curvature
+## Gaussian Curvature
 Gaussian curvature on a continuous surface is defined as the product of the
 principal curvatures:
 
@@ -158,7 +167,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=curv></a> Curvature Directions
+## 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
@@ -216,7 +225,7 @@ int main(int argc, char * argv[])
 }
 ```
 
-## <a id=grad></a> Gradient
+## Gradient
 Scalar functions on a surface can be discretized as a piecewise linear function
 with values defined at each mesh vertex:
 
@@ -243,7 +252,7 @@ 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
+where $\mathbf{f}$ is $n\times 1$ and $\mathbf{G}$ is an $md\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
@@ -252,7 +261,7 @@ 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
+## Laplacian
 
 The discrete Laplacian is an essential geometry processing tool. Many
 interpretations and flavors of the Laplace and Laplace-Beltrami operator exist. 
@@ -347,6 +356,146 @@ sampling function values at vertices:
 
  $\int_S x\, y\ dA \approx \mathbf{x}^T\mathbf{M}\,\mathbf{y}.$
 
+An alternative mass matrix $\mathbf{T}$ is a $md \times md$ 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).
+
+### Alternative construction of Laplacian
+
+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):
+
+  $\int_S \nabla f \nabla f dA = \int_S f \Delta f dA$
+
+Or in matrix form which is immediately translatable to code:
+
+  $\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}.$
+
+So we have that $\mathbf{L} = \mathbf{G}^T \mathbf{T} \mathbf{G}$. This also
+hints that we may consider $\mathbf{G}^T$ as a discrete _divergence_ operator,
+since the Laplacian is the divergence of gradient. Naturally, $\mathbf{G}^T$ is
+$n \times md$ sparse matrix which takes vector values stored at triangle faces
+to scalar divergence values at vertices.
+
+# Chapter 3: Matrices and Linear Algebra
+Libigl relies heavily on the Eigen library for dense and sparse linear algebra
+routines. Besides geometry processing routines, libigl has a few linear algebra
+routines which bootstrap Eigen and make Eigen feel even more like a high-level
+algebra library like Matlab.
+
+## Slice
+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:
+
+```matlab
+B = A(R,C);
+```
+
+If `A` is a $m \times n$ matrix and `R` is a $j$-long list of row-indices
+(between 1 and $m$) and `C` is a $k$-long list of column-indices, then as a
+result `B` will be a $j \times k$ matrix drawing elements from `A` according to
+`R` and `C`. In libigl, the same functionality is provided by the `slice`
+function:
+
+```cpp
+VectorXi R,C;
+MatrixXd A,B;
+...
+igl::slice(A,R,C,B);
+```
+
+`A` and `B` could also be sparse matrices.
+
+Similarly, consider the matlab code:
+
+```matlab
+A(R,C) = B;
+```
+
+Now, the selection is on the left-hand side so the $j \times k$ matrix  `B` is
+being _written into_ the submatrix of `A` determined by `R` and `C`. This
+functionality is provided in libigl using `slice_into`:
+
+```cpp
+igl::slice_into(B,R,C,A);
+```
+
+![The example `Slice` shows how to use `igl::slice` to change the colors for triangles
+on a mesh.](images/decimated-knight-slice-color.jpg)
+
+## Sort
+
+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:
+
+```matlab
+[Y,I] = sort(X,1,'ascend');
+```
+
+so if `X` is a $m \times n$ matrix then `Y` will also be an $m \times n$ matrix
+with entries sorted along dimension `1` in `'ascend'`ing order. The second
+output `I` is a $m \times n$ matrix of indices such that `Y(i,j) =
+X(I(i,j),j);`. That is, `I` reveals how `X` is sorted into `Y`.
+
+This same functionality is supported in libigl:
+
+```cpp
+igl::sort(X,1,true,Y,I);
+```
+
+Similarly, sorting entire rows can be accomplished in matlab using:
+
+```matlab
+[Y,I] = sortrows(X,'ascend');
+```
+
+where now `I` is a $m$ vector of indices such that `Y = X(I,:)`.
+
+In libigl, this is supported with
+
+```cpp
+igl::sortrows(X,true,Y,I);
+```
+where again `I` reveals the index of sort so that it can be reproduced with
+`igl::slice(X,I,1,Y)`.
+
+Analogous functions are available in libigl for: `max`, `min`, and `unique`.
+
+![The example `Sort` shows how to use `igl::sortrows` to 
+pseudocolor triangles according to their barycenters' sorted
+order.](images/decimated-knight-sort-color.jpg)
+
+
+### Other Matlab-style functions
+Libigl implements a variety of other routines with the same api and
+functionality as common matlab functions.
+
+- `igl::any_of` Whether any elements are non-zero (true)
+- `igl::cat` Concatenate two matrices (especially useful for dealing with Eigen
+  sparse matrices)
+- `igl::ceil` Round entries up to nearest integer
+- `igl::cumsum` Cumulative sum of matrix elements
+- `igl::colon` Act like Matlab's `:`, similar to Eigen's `LinSpaced`
+- `igl::cross` Cross product per-row
+- `igl::dot` dot product per-row
+- `igl::find` Find subscripts of non-zero entries
+- `igl::floot` Round entries down to nearest integer
+- `igl::histc` Counting occurrences for building a histogram
+- `igl::hsv_to_rgb` Convert HSV colors to RGB (cf. Matlab's `hsv2rgb`)
+- `igl::intersect` Set intersection of matrix elements.
+- `igl::jet` Quantized colors along the rainbow.
+- `igl::kronecker_product` Compare to Matlab's `kronprod`
+- `igl::median` Compute the median per column
+- `igl::mode` Compute the mode per column
+- `igl::orth` Orthogonalization of a basis
+- `igl::speye` Identity as sparse matrix
+
+
 [#meyer_2003]: Mark Meyer and Mathieu Desbrun and Peter Schröder and Alan H.  Barr,
 "Discrete Differential-Geometry Operators for Triangulated
 2-Manifolds," 2003.