Merge branch 'master' of github.com:libigl/libigl into alecjacobson

Former-commit-id: 9a805a9764880a1c51d6ccbb5240dc6d0f8399bb

RELEASE_HISTORY.md

@@ -37,7 +37,7 @@ Version | Short description
 0.3.1   | Linearly dependent constraints in min_quad_with_fixed, SparseQR buggy
 0.3.0   | Better active set method support
 0.2.3   | More explicits, active set method, opengl/anttweakbar guards
-0.2.2   | More explicit instanciations, faster sorts and uniques
+0.2.2   | More explicit instantiations, faster sorts and uniques
 0.2.1   | Bug fixes in barycenter and doublearea found by Martin Bisson
 0.2.0   | XML serializer more stable and fixed bug in remove_duplicate_vertices
 0.1.8   | Embree and xml (windows only) extras

file-formats/dmat.html

@@ -12,7 +12,7 @@
     <p>
 A .dmat file contains a dense matrix in column major order. It can contain
 ASCII or binary data. Note that it is uncompressed so binary only reduces the
-file size by 50%. But writing and reading binary is usualy faster. In MATLAB,
+file size by 50%. But writing and reading binary is usually faster. In MATLAB,
 binary is almost 100x faster.
     </p>
     <h2>ASCII</h2>
@@ -33,7 +33,7 @@ Then there should be another header containing the size of the binary part:
     </p>
     <pre><code>[#cols] [#rows]</code></pre>
     <p>
-Then coefficents are written in column-major order in Little-endian 8-byte
+Then coefficients are written in column-major order in Little-endian 8-byte
 double precision IEEE floating point format.
     </p>
     <p><strong>Note:</strong> Line endings must be <code>'\n'</code> aka char(10) aka line feeds.</p>

file-formats/index.html

@@ -20,7 +20,7 @@
     <li><a href=http://tetgen.berlios.de/fformats.face.html
     class=missing>.face</a> TetGen's file format for simplicial facets.</li>
     <li><a href="http://www.ann.jussieu.fr/frey/publications/RT-0253.pdf#page=33">.mesh</a> Medit's triangle surface mesh + tetrahedral volume mesh file format, see page 33, section 7.2.1</li>
-    <li><i>.node</i> List of points (vertices). Described indentically (upto
+    <li><i>.node</i> List of points (vertices). Described identically (upto
         accepted dimensions, use of attributes and boundary markers) by <a
     href="https://www.cs.cmu.edu/~quake/triangle.node.html">Triangle</a>, <a
     href=http://tetgen.berlios.de/fformats.node.html>TetGen</a>, and <a
@@ -37,7 +37,7 @@
     href=https://en.wikipedia.org/wiki/Portable_Network_Graphics>.png</a>
     Portable Network Graphics image file. IGLLIB (in the libiglpng extra)
   supports png image files via the <a href=https://github.com/yig/yimg>yimg</a>
-  library. Alpha channels and compression are suppported.</li>
+  library. Alpha channels and compression are supported.</li>
     <li><i>.poly</i> Piecewise-linear complex. This format comes in many similar but importantly different flavors: 
       <a href="https://www.cs.cmu.edu/~quake/triangle.poly.html">triangle's</a>, <a href="http://tetgen.berlios.de/fformats.poly.html">tetgen's</a>, <a href="http://sparse-meshing.com/svr/0.2.1/format-poly.html">pyramid/SVR's</a></li>
     <li><a href=./rbr.html>.rbr</a> ASCII files for saving state of ReAntTweakBar</li>

file-formats/tgf.html

@@ -11,7 +11,7 @@
     <hr>
     <p>
 A .tgf file contains a graph of describing a set of control handles/structures:
-point controls, bones of a skelton and cages made of "cage edges".
+point controls, bones of a skeleton and cages made of "cage edges".
     </p>
     <p>
 The first part of the file consists of lines regarding each vertex of the graph. Each line reads:

file-formats/xml.html

@@ -31,7 +31,7 @@ STL containers: <b>std::array, std::vector, std::pair</b><br/>
 Eigen types: <b>Eigen::Matrix, Eigen::SparseMatrix</b><br/>
 User defined types: <b>XMLSerializable*.</b><br/>
 <br/>
-There can also be a hierachical structure like vector&lt;int&gt;, this will result in the following serialization:
+There can also be a hierarchical structure like vector&lt;int&gt;, this will result in the following serialization:
     </p>
     <pre><code>&lt;group&gt;
   &lt;vector size="3"&gt;

include/igl/AABB.cpp

@@ -587,7 +587,7 @@ IGL_INLINE typename igl::AABB<DerivedV,DIM>::Scalar
     return sqr_d;
   }
 
-  //// Exact minimum squared distance between arbitary primitives inside this and
+  //// Exact minimum squared distance between arbitrary primitives inside this and
   //// othre's bounding boxes
   //const auto & min_squared_distance = [&](
   //  const AABB<DerivedV,DIM> * A,

include/igl/AtA_cached.cpp

@@ -0,0 +1,130 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@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 "AtA_cached.h"
+
+#include <iostream>
+#include <vector>
+#include <utility>
+
+template <typename Scalar>
+IGL_INLINE void igl::AtA_cached_precompute(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    igl::AtA_cached_data& data)
+{
+  // 1 Compute At (this could be avoided, but performance-wise it will not make a difference)
+  std::vector<std::vector<int> > Col_RowPtr;
+  std::vector<std::vector<int> > Col_IndexPtr;
+
+  Col_RowPtr.resize(A.cols());
+  Col_IndexPtr.resize(A.cols());
+
+  for (unsigned k=0; k<A.outerSize(); ++k)
+  {
+    unsigned outer_index = *(A.outerIndexPtr()+k);
+    unsigned next_outer_index = (k+1 == A.outerSize()) ? A.nonZeros() : *(A.outerIndexPtr()+k+1); 
+    
+    for (unsigned l=outer_index; l<next_outer_index; ++l)
+    {
+      int col = k;
+      int row = *(A.innerIndexPtr()+l);
+      int value_index = l;
+      assert(col < A.cols());
+      assert(col >= 0);
+      assert(row < A.rows());
+      assert(row >= 0);
+      assert(value_index >= 0);
+      assert(value_index < A.nonZeros());
+
+      Col_RowPtr[col].push_back(row);
+      Col_IndexPtr[col].push_back(value_index);
+    }
+  }
+
+  Eigen::SparseMatrix<Scalar> At = A.transpose();
+  At.makeCompressed();
+  AtA = At * A;
+  AtA.makeCompressed();
+
+  assert(AtA.isCompressed());
+
+  // If weights are not provided, use 1
+  if (data.W.size() == 0)
+    data.W = Eigen::VectorXd::Ones(A.rows());
+  assert(data.W.size() == A.rows());
+
+  data.I_outer.reserve(AtA.outerSize());
+  data.I_row.reserve(2*AtA.nonZeros());
+  data.I_col.reserve(2*AtA.nonZeros());
+  data.I_w.reserve(2*AtA.nonZeros());
+
+  // 2 Construct the rules
+  for (unsigned k=0; k<AtA.outerSize(); ++k)
+  {
+    unsigned outer_index = *(AtA.outerIndexPtr()+k);
+    unsigned next_outer_index = (k+1 == AtA.outerSize()) ? AtA.nonZeros() : *(AtA.outerIndexPtr()+k+1); 
+    
+    for (unsigned l=outer_index; l<next_outer_index; ++l)
+    {
+      int col = k;
+      int row = *(AtA.innerIndexPtr()+l);
+      int value_index = l;
+      assert(col < AtA.cols());
+      assert(col >= 0);
+      assert(row < AtA.rows());
+      assert(row >= 0);
+      assert(value_index >= 0);
+      assert(value_index < AtA.nonZeros());
+
+      data.I_outer.push_back(data.I_row.size());
+
+      // Find correspondences
+      unsigned i=0;
+      unsigned j=0;
+      while (i<Col_RowPtr[row].size() && j<Col_RowPtr[col].size())
+      {
+          if (Col_RowPtr[row][i] == Col_RowPtr[col][j])
+          {
+            data.I_row.push_back(Col_IndexPtr[row][i]);
+            data.I_col.push_back(Col_IndexPtr[col][j]);
+            data.I_w.push_back(Col_RowPtr[col][j]);
+            ++i;
+            ++j;
+          } else 
+          if (Col_RowPtr[row][i] > Col_RowPtr[col][j])
+            ++j;
+          else
+            ++i;
+
+      }
+    }
+  }
+  data.I_outer.push_back(data.I_row.size()); // makes it more efficient to iterate later on
+
+  igl::AtA_cached(A,AtA,data);
+}
+
+template <typename Scalar>
+IGL_INLINE void igl::AtA_cached(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    const igl::AtA_cached_data& data)
+{
+  for (unsigned i=0; i<data.I_outer.size()-1; ++i)
+  {
+    *(AtA.valuePtr() + i) = 0;
+    for (unsigned j=data.I_outer[i]; j<data.I_outer[i+1]; ++j)
+      *(AtA.valuePtr() + i) += *(A.valuePtr() + data.I_row[j]) * data.W[data.I_w[j]] * *(A.valuePtr() + data.I_col[j]);
+  }
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::AtA_cached<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, igl::AtA_cached_data const&);
+template void igl::AtA_cached_precompute<double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, igl::AtA_cached_data&);
+#endif

include/igl/AtA_cached.h

@@ -0,0 +1,68 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@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_ATA_CACHED_H
+#define IGL_ATA_CACHED_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{  
+  struct AtA_cached_data
+  {
+    // Weights
+    Eigen::VectorXd W;
+
+    // Flatten composition rules
+    std::vector<int> I_row;
+    std::vector<int> I_col;
+    std::vector<int> I_w;
+
+    // For each entry of AtA, points to the beginning
+    // of the composition rules
+    std::vector<int> I_outer;
+  };
+
+  // Computes At * W * A, where A is sparse and W is diagonal. Divides the 
+  // construction in two phases, one
+  // for fixing the sparsity pattern, and one to populate it with values. Compared to
+  // evaluating it directly, this version is slower for the first time (since it requires a
+  // precomputation), but faster to the subsequent evaluations.
+  //
+  // Input:
+  //   A m x n sparse matrix
+  //   data stores the precomputed sparsity pattern, data.W contains the optional diagonal weights (stored as a dense vector). If W is not provided, it is replaced by the identity.
+  // Outputs:
+  //   AtA  m by m matrix computed as AtA * W * A
+  //
+  // Example:
+  // AtA_data = igl::AtA_cached_data();
+  // AtA_data.W = W;
+  // if (s.AtA.rows() == 0)
+  //   igl::AtA_cached_precompute(s.A,s.AtA,s.AtA_data);
+  // else
+  //   igl::AtA_cached(s.A,s.AtA,s.AtA_data);
+  template <typename Scalar>
+  IGL_INLINE void AtA_cached_precompute(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    AtA_cached_data& data);
+
+  template <typename Scalar>
+  IGL_INLINE void AtA_cached(
+    const Eigen::SparseMatrix<Scalar>& A,
+    Eigen::SparseMatrix<Scalar>& AtA,
+    const AtA_cached_data& data);
+  
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "AtA_cached.cpp"
+#endif
+
+#endif

include/igl/Camera.h

@@ -23,7 +23,7 @@ namespace igl
 
   // A simple camera class. The camera stores projection parameters (field of
   // view angle, aspect ratio, near and far clips) as well as a rigid
-  // tranformation *of the camera as if it were also a scene object*. Thus, the
+  // transformation *of the camera as if it were also a scene object*. Thus, the
   // **inverse** of this rigid transformation is the modelview transformation.
   class Camera
   {

include/igl/REDRUM.h

@@ -13,7 +13,7 @@
 // A: I guess the right way is to not use a macro but a proper function with
 // streams as input and output.
 
-// ANSI color codes for formating iostream style output
+// ANSI color codes for formatting iostream style output
 
 #ifdef IGL_REDRUM_NOOP
 

include/igl/SortableRow.h

@@ -15,7 +15,7 @@
 namespace igl
 {
   // Templates:
-  //   T  should be a matrix that implments .size(), and operator(int i)
+  //   T  should be a matrix that implements .size(), and operator(int i)
   template <typename T>
   class SortableRow
   {

include/igl/WindingNumberTree.h

@@ -52,7 +52,7 @@ namespace igl
       MatrixXS SV;
       // Facets in this bounding volume
       MatrixXF F;
-      // Tesselated boundary curve
+      // Tessellated boundary curve
       MatrixXF cap;
       // Upper Bound on radius of enclosing ball
       typename DerivedV::Scalar radius;
@@ -97,7 +97,7 @@ namespace igl
       // Same as above, but always computes winding number using exact method
       // (sum over every facet)
       inline typename DerivedV::Scalar winding_number_all(const Point & p) const;
-      // Same as above, but always computes using sum over tesslated boundary
+      // Same as above, but always computes using sum over tessllated boundary
       inline typename DerivedV::Scalar winding_number_boundary(const Point & p) const;
       //// Same as winding_number above, but if max_simple_abs_winding_number is
       //// less than some threshold min_max_w just return 0 (colloquially the "fast
@@ -493,7 +493,7 @@ igl::WindingNumberTree<Point,DerivedV,DerivedF>::cached_winding_number(
   return 0;
 }
 
-// Explicit instanciation of static variable
+// Explicit instantiation of static variable
 template <
   typename Point,
   typename DerivedV, 

include/igl/arap_dof.cpp

@@ -34,7 +34,7 @@
 // defined if no early exit is supported, i.e., always take a fixed number of iterations
 #define IGL_ARAP_DOF_FIXED_ITERATIONS_COUNT
 
-// A carefull derivation of this implementation is given in the corresponding
+// A careful derivation of this implementation is given in the corresponding
 // matlab function arap_dof.m
 template <typename LbsMatrixType, typename SSCALAR>
 IGL_INLINE bool igl::arap_dof_precomputation(
@@ -339,7 +339,7 @@ namespace igl
   }
   
   // converts CSM_M_SSCALAR[0], CSM_M_SSCALAR[1], CSM_M_SSCALAR[2] into one
-  // "condensed" matrix CSM while checking we're not loosing any information by
+  // "condensed" matrix CSM while checking we're not losing any information by
   // this process; specifically, returns maximal difference from scaled 3x3
   // identity blocks, which should be pretty small number
   template <typename MatrixXS>
@@ -449,7 +449,7 @@ namespace igl
   }
   
   // converts "Solve1" the "rotations" part of FullSolve matrix (the first part)
-  // into one "condensed" matrix CSolve1 while checking we're not loosing any
+  // into one "condensed" matrix CSolve1 while checking we're not losing any
   // information by this process; specifically, returns maximal difference from
   // scaled 3x3 identity blocks, which should be pretty small number
   template <typename MatrixXS>

include/igl/arap_dof.h

@@ -127,7 +127,7 @@ namespace igl
   //   L0  #handles * dim * dim+1 list of initial guess transformation entries,
   //     also holds fixed transformation entries for fixed handles
   //   max_iters  maximum number of iterations
-  //   tol  stopping critera parameter. If variables (linear transformation
+  //   tol  stopping criteria parameter. If variables (linear transformation
   //     matrix entries) change by less than 'tol' the optimization terminates,
   //       0.75 (weak tolerance)
   //       0.0 (extreme tolerance)

include/igl/bbw.h

@@ -49,7 +49,7 @@ namespace igl
   //   Ele  #Elements by simplex-size list of element indices
   //   b  #b boundary indices into V
   //   bc #b by #W list of boundary values
-  //   data  object containing options, intial guess --> solution and results
+  //   data  object containing options, initial guess --> solution and results
   // Outputs:
   //   W  #V by #W list of *unnormalized* weights to normalize use
   //    igl::normalize_row_sums(W,W);

include/igl/bijective_composite_harmonic_mapping.h

@@ -29,7 +29,7 @@ namespace igl
   //   bc  #b by 2 list of boundary conditions corresponding to b
   // Outputs:
   //   U  #V by 2 list of output mesh vertex locations
-  // Returns true if and only if U contains a successfull bijectie mapping
+  // Returns true if and only if U contains a successful bijectie mapping
   //
   // 
   template <
@@ -46,12 +46,12 @@ namespace igl
     Eigen::PlainObjectBase<DerivedU> & U);
   //
   // Inputs:
-  //   min_steps  mininum number of steps to take from V(b,:) to bc
-  //   max_steps  mininum number of steps to take from V(b,:) to bc (if
+  //   min_steps  minimum number of steps to take from V(b,:) to bc
+  //   max_steps  minimum number of steps to take from V(b,:) to bc (if
   //     max_steps == min_steps then no further number of steps will be tried)
   //   num_inner_iters  number of iterations of harmonic solves to run after
   //     for each morph step (to try to push flips back in)
-  //   test_for_flips  wether to check if flips occured (and trigger more
+  //   test_for_flips  whether to check if flips occurred (and trigger more
   //     steps). if test_for_flips = false then this function always returns
   //     true
   // 

include/igl/boundary_conditions.cpp

@@ -146,7 +146,7 @@ IGL_INLINE bool igl::boundary_conditions(
     bc(bim[bci[i]],bcj[i]) = bcv[i];
   }
 
-  // Normalize accross rows so that conditions sum to one
+  // Normalize across rows so that conditions sum to one
   for(i = 0;i<bc.rows();i++)
   {
     double sum = bc.row(i).sum();

include/igl/cat.h

@@ -25,7 +25,7 @@ namespace igl
   // Template:
   //   Scalar  scalar data type for sparse matrices like double or int
   //   Mat  matrix type for all matrices (e.g. MatrixXd, SparseMatrix)
-  //   MatC  matrix type for ouput matrix (e.g. MatrixXd) needs to support
+  //   MatC  matrix type for output matrix (e.g. MatrixXd) needs to support
   //     resize
   // Inputs:
   //   A  first input matrix

include/igl/connect_boundary_to_infinity.h

@@ -11,7 +11,7 @@
 #include <Eigen/Core>
 namespace igl
 {
-  // Connect all boundary edges to a ficticious point at infinity.
+  // Connect all boundary edges to a fictitious point at infinity.
   //
   // Inputs:
   //   F  #F by 3 list of face indices into some V

include/igl/copyleft/cgal/SelfIntersectMesh.h

@@ -236,7 +236,7 @@ namespace igl
 // a 2-manifold.
 // A: But! It seems we could use CGAL::Triangulation_3. Though it won't be easy
 // to take advantage of functions like insert_in_facet because we want to
-// constrain segments. Hmmm. Actualy Triangulation_3 doesn't look right...
+// constrain segments. Hmmm. Actually Triangulation_3 doesn't look right...
 
 // CGAL's box_self_intersection_d uses C-style function callbacks without
 // userdata. This is a leapfrog method for calling a member function. It should

include/igl/copyleft/cgal/cell_adjacency.h

@@ -27,7 +27,7 @@ namespace igl
       //   num_cells        number of cells.
       //
       // Outputs:
-      //   adjacency_list  #C array of list of adjcent cell informations.  If
+      //   adjacency_list  #C array of list of adjcent cell information.  If
       //                   cell i and cell j are adjacent via patch x, where i
       //                   is on the positive side of x, and j is on the
       //                   negative side.  Then,

include/igl/copyleft/cgal/closest_facet.cpp

@@ -124,7 +124,7 @@ IGL_INLINE void igl::copyleft::cgal::closest_facet(
       case CGAL::COPLANAR:
         // Warning:
         // This can only happen if fid contains a boundary edge.
-        // Catergorized this ambiguous case as negative side.
+        // Categorized this ambiguous case as negative side.
         return false;
       default:
         throw std::runtime_error("Unknown CGAL state.");

include/igl/copyleft/cgal/order_facets_around_edge.h

@@ -31,7 +31,7 @@ namespace igl
       //   V  #V by 3 list of vertices.
       //   F  #F by 3 list of faces
       //   s  Index of source vertex.
-      //   d  Index of desination vertex.
+      //   d  Index of destination vertex.
       //   adj_faces  List of adjacent face signed indices.
       // Output:
       //   order  List of face indices that orders adjacent faces around edge
@@ -50,7 +50,7 @@ namespace igl
           Eigen::PlainObjectBase<DerivedI>& order,
           bool debug=false);
 
-      // This funciton is a wrapper around the one above.  Since the ordering
+      // This function is a wrapper around the one above.  Since the ordering
       // is circular, the pivot point is used to define a starting point.  So
       // order[0] is the index into adj_face that is immediately after the
       // pivot face (s, d, pivot point) in clockwise order.

include/igl/copyleft/cgal/outer_facet.h

@@ -26,7 +26,7 @@ namespace igl
         // See cgal::remesh_self_intersections.h for how to obtain such input.
         //
         // This function differ from igl::outer_facet() in the fact this
-        // funciton is more robust because it does not rely on facet normals.
+        // function is more robust because it does not rely on facet normals.
         //
         // Inputs:
         //   V  #V by 3 list of vertex positions

include/igl/copyleft/cgal/point_mesh_squared_distance.cpp

@@ -81,7 +81,7 @@ IGL_INLINE void igl::copyleft::cgal::point_mesh_squared_distance_precompute(
   // accelerate_distance_queries doesn't seem actually to do _all_ of the
   // precomputation. the tree (despite being const) will still do more
   // precomputation and reorganizing on the first call of `closest_point` or
-  // `closest_point_and_primitive`. Therefor, call it once here.
+  // `closest_point_and_primitive`. Therefore, call it once here.
   tree.closest_point_and_primitive(Point_3(0,0,0));
 }
 

include/igl/copyleft/cgal/points_inside_component.cpp

@@ -321,7 +321,7 @@ IGL_INLINE void igl::copyleft::cgal::points_inside_component(
                 inside(i,0) = determine_point_face_orientation(V, F, I, query, fid);
                 break;
             default:
-                throw "Unknow closest element type!";
+                throw "Unknown closest element type!";
         }
     }
 }

include/igl/copyleft/cgal/remesh_intersections.cpp

@@ -453,7 +453,7 @@ IGL_INLINE void igl::copyleft::cgal::remesh_intersections(
           FF.data(), [&vv_to_unique](const typename DerivedFF::Scalar& a)
           { return vv_to_unique[a]; });
       IM.resize(unique_vv.rows());
-      // Have to use << instead of = becasue Eigen's PlainObjectBase is annoying
+      // Have to use << instead of = because Eigen's PlainObjectBase is annoying
       IM << igl::LinSpaced<
         Eigen::Matrix<typename DerivedIM::Scalar, Eigen::Dynamic,1 >
         >(unique_vv.rows(), 0, unique_vv.rows()-1);

include/igl/copyleft/cgal/remesh_intersections.h

@@ -30,7 +30,7 @@ namespace igl
       //   offending #offending map taking face indices into F to pairs of order
       //     of first finding and list of intersection objects from all
       //     intersections
-      //   stitch_all  if true, merge all vertices with thte same coordiante.
+      //   stitch_all  if true, merge all vertices with the same coordinate.
       // Outputs:
       //   VV  #VV by 3 list of vertex positions, if stitch_all = false then
       //     first #V vertices will always be V

include/igl/copyleft/cgal/trim_with_solid.cpp

@@ -56,7 +56,7 @@ IGL_INLINE void igl::copyleft::cgal::trim_with_solid(
   igl::slice_mask(Eigen::MatrixXi(F),A,1,F);
   igl::slice_mask(Eigen::VectorXi(P),A,1,P);
   igl::slice_mask(Eigen::VectorXi(J),A,1,J);
-  // Agregate representative query points for each patch
+  // Aggregate representative query points for each patch
   std::vector<bool> flag(num_patches);
   std::vector<std::vector<CGAL::Epeck::FT> > vQ;
   Eigen::VectorXi P2Q(num_patches);

include/igl/copyleft/comiso/miq.h

@@ -43,7 +43,7 @@ namespace igl
     //   UV             #UV by 2 list of vertices in 2D
     //   FUV            #FUV by 3 list of face indices in UV
     //
-    // TODO: rename the parameters name in the cpp consistenly
+    // TODO: rename the parameters name in the cpp consistently
     //       improve the handling of hard_features, right now it might fail in difficult cases
 
     template <typename DerivedV, typename DerivedF, typename DerivedU>

include/igl/copyleft/comiso/nrosy.h

@@ -32,7 +32,7 @@ namespace igl
     //   w_soft  #S by 1 weight for the soft constraints (0-1)
     //   bc_soft #S by 3 bc for soft constraints
     //   N       the degree of the N-RoSy vector field
-    //   soft    the strenght of the soft contraints w.r.t. smoothness
+    //   soft    the strength of the soft constraints w.r.t. smoothness
     //           (0 -> smoothness only, 1->constraints only)
     // Outputs:
     //   R       #F by 3 the representative vectors of the interpolated field

include/igl/copyleft/offset_surface.h

@@ -8,7 +8,7 @@ namespace igl
 {
   namespace copyleft
   {
-    // Compute a triangulated offset surface using maching cubes on a gride of
+    // Compute a triangulated offset surface using matching cubes on a grid of
     // signed distance values from the input triangle mesh.
     //
     // Inputs:

include/igl/copyleft/opengl2/render_to_tga.h

@@ -20,7 +20,7 @@ namespace igl
     //   width  width of scene and resulting image
     //   height height of scene and resulting image
     ///  alpha  whether to include alpha channel
-    // Returns true only if no errors occured
+    // Returns true only if no errors occurred
     //
     // See also: png/render_to_png which is slower but writes .png files
     IGL_INLINE bool render_to_tga(

include/igl/copyleft/tetgen/cdt.h

@@ -34,7 +34,7 @@ namespace igl
         // Flags to tetgen. Do not include the "c" flag here! {"Y"}
         std::string flags = "Y";
       };
-      // Create a constrained delaunay tesselation containing convex hull of the
+      // Create a constrained delaunay tessellation containing convex hull of the
       // given **non-selfintersecting** mesh.
       //
       // Inputs:

include/igl/copyleft/tetgen/read_into_tetgenio.h

@@ -30,8 +30,8 @@ namespace igl
       //   .medit
       //   .vtk
       //   etc.
-      // Noteably it does not support .obj which is loaded by hand here (also
-      // demonstrating how to load points/faces programatically)
+      // Notably it does not support .obj which is loaded by hand here (also
+      // demonstrating how to load points/faces programmatically)
       //
       // If the file extension is not recognized the filename is assumed to be
       // the basename of a collection describe a tetmesh, (of which at least

include/igl/cotmatrix.h

@@ -16,7 +16,7 @@
 //  Used const references rather than copying the entire mesh 
 //    Alec 9 October 2011
 //  removed cotan (uniform weights) optional parameter it was building a buggy
-//    half of the uniform laplacian, please see adjacency_matrix istead 
+//    half of the uniform laplacian, please see adjacency_matrix instead 
 //    Alec 9 October 2011
 
 namespace igl 

include/igl/cotmatrix_entries.cpp

@@ -35,10 +35,10 @@ IGL_INLINE void igl::cotmatrix_entries(
     case 3:
     {
       // Triangles
-      //Compute Squared Edge lenghts 
+      //Compute Squared Edge lengths 
       Matrix<typename DerivedC::Scalar,Dynamic,3> l2;
       igl::squared_edge_lengths(V,F,l2);
-      //Compute Edge lenghts 
+      //Compute Edge lengths 
       Matrix<typename DerivedC::Scalar,Dynamic,3> l;
       l = l2.array().sqrt();
       

include/igl/crouzeix_raviart_cotmatrix.cpp

@@ -19,7 +19,7 @@ void igl::crouzeix_raviart_cotmatrix(
   Eigen::PlainObjectBase<DerivedE> & E,
   Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
 {
-  // All occurances of directed "facets"
+  // All occurrences of directed "facets"
   Eigen::MatrixXi allE;
   oriented_facets(F,allE);
   Eigen::VectorXi _1;

include/igl/crouzeix_raviart_massmatrix.cpp

@@ -24,7 +24,7 @@ void igl::crouzeix_raviart_massmatrix(
     Eigen::PlainObjectBase<DerivedE> & E,
     Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
 {
-  // All occurances of directed "facets"
+  // All occurrences of directed "facets"
   Eigen::MatrixXi allE;
   oriented_facets(F,allE);
   Eigen::VectorXi _1;

include/igl/cut_mesh.cpp

@@ -70,11 +70,11 @@ namespace igl {
     IGL_INLINE void FindInitialPos(const int vert, int &edge, int &face);
 
 
-    // intialize the mapping given an initial pos
+    // initialize the mapping given an initial pos
     // whih must be initialized with FindInitialPos
     IGL_INLINE void MapIndexes(const int  vert, const int edge_init, const int f_init);
 
-    // intialize the mapping for a given vertex
+    // initialize the mapping for a given vertex
     IGL_INLINE void InitMappingSeam(const int vert);
 
   };
@@ -187,7 +187,7 @@ FindInitialPos(const int vert,
 
 
 
-///intialize the mapping given an initial pos
+///initialize the mapping given an initial pos
 ///whih must be initialized with FindInitialPos
 template <typename DerivedV, typename DerivedF, typename VFType, typename DerivedTT, typename DerivedC>
 IGL_INLINE void igl::MeshCutterMini<DerivedV, DerivedF, VFType, DerivedTT, DerivedC>::

include/igl/directed_edge_parents.h

@@ -13,7 +13,7 @@
 
 namespace igl
 {
-  // Recover "parents" (preceeding edges) in a tree given just directed edges.
+  // Recover "parents" (preceding edges) in a tree given just directed edges.
   //
   // Inputs:
   //   E  #E by 2 list of directed edges

include/igl/embree/EmbreeIntersector.h

@@ -209,7 +209,7 @@ inline void igl::embree::EmbreeIntersector::global_init()
   {
     rtcInit();
     if(rtcGetError() != RTC_NO_ERROR)
-      std::cerr << "Embree: An error occured while initialiting embree core!" << std::endl;
+      std::cerr << "Embree: An error occurred while initializing embree core!" << std::endl;
 #ifdef IGL_VERBOSE
     else
       std::cerr << "Embree: core initialized." << std::endl;
@@ -324,7 +324,7 @@ inline void igl::embree::EmbreeIntersector::init(
   rtcCommit(scene);
 
   if(rtcGetError() != RTC_NO_ERROR)
-      std::cerr << "Embree: An error occured while initializing the provided geometry!" << endl;
+      std::cerr << "Embree: An error occurred while initializing the provided geometry!" << endl;
 #ifdef IGL_VERBOSE
   else
     std::cerr << "Embree: geometry added." << endl;
@@ -348,7 +348,7 @@ void igl::embree::EmbreeIntersector::deinit()
 
     if(rtcGetError() != RTC_NO_ERROR)
     {
-        std::cerr << "Embree: An error occured while resetting!" << std::endl;
+        std::cerr << "Embree: An error occurred while resetting!" << std::endl;
     }
 #ifdef IGL_VERBOSE
     else
@@ -374,7 +374,7 @@ inline bool igl::embree::EmbreeIntersector::intersectRay(
   rtcIntersect(scene,ray);
 #ifdef IGL_VERBOSE
   if(rtcGetError() != RTC_NO_ERROR)
-      std::cerr << "Embree: An error occured while resetting!" << std::endl;
+      std::cerr << "Embree: An error occurred while resetting!" << std::endl;
 #endif
 
   if((unsigned)ray.geomID != RTC_INVALID_GEOMETRY_ID)

include/igl/exterior_edges.cpp

@@ -67,7 +67,7 @@ IGL_INLINE void igl::exterior_edges(
   {
     int e = 0;
     const size_t nue = uE.rows();
-    // Append each unique edge with a non-zero amount of signed occurances
+    // Append each unique edge with a non-zero amount of signed occurrences
     for(size_t ue = 0; ue<nue; ue++)
     {
       const int count = counts(ue);

include/igl/exterior_edges.h

@@ -12,7 +12,7 @@
 namespace igl
 {
   // EXTERIOR_EDGES Determines boundary "edges" and also edges with an
-  // odd number of occurances where seeing edge (i,j) counts as +1 and seeing
+  // odd number of occurrences where seeing edge (i,j) counts as +1 and seeing
   // the opposite edge (j,i) counts as -1
   //
   // Inputs:

include/igl/facet_components.h

@@ -16,7 +16,7 @@ namespace igl
   //
   // Inputs:
   //   F  #F by 3 list of triangle indices
-  // Ouputs:
+  // Outputs:
   //   C  #F list of connected component ids
   template <typename DerivedF, typename DerivedC>
   IGL_INLINE void facet_components(
@@ -25,7 +25,7 @@ namespace igl
   // Inputs:
   //   TT  #TT by 3 list of list of adjacency triangles (see
   //   triangle_triangle_adjacency.h)
-  // Ouputs:
+  // Outputs:
   //   C  #F list of connected component ids
   //   counts #C list of number of facets in each components
   template <

include/igl/hausdorff.h

@@ -59,7 +59,7 @@ namespace igl
   //   V   3 by 3 list of corner positions so that V.row(i) is the position of the
   //     ith corner
   //   dist_to_B  function taking the x,y,z coordinate of a query position and
-  //     outputing the closest-point distance to some point-set B
+  //     outputting the closest-point distance to some point-set B
   // Outputs:
   //   l  lower bound on Hausdorff distance 
   //   u  upper bound on Hausdorff distance

include/igl/histc.h

@@ -13,7 +13,7 @@
 
 namespace igl
 {
-  // Like matlab's histc. Count occurances of values in X between consecutive
+  // Like matlab's histc. Count occurrences of values in X between consecutive
   // entries in E
   //
   // Inputs:

include/igl/integrable_polyvector_fields.h

@@ -15,7 +15,7 @@
 
 namespace igl {
   // Compute a curl-free frame field from user constraints, optionally starting
-  // from a gived frame field (assumed to be interpolating the constraints).
+  // from a given frame field (assumed to be interpolating the constraints).
   // Implementation of the paper "Integrable PolyVector Fields", SIGGRAPH 2015.
 
   // Set of parameters used during solve
@@ -51,7 +51,7 @@ namespace igl {
                                                     igl::IntegrableFieldSolverData<DerivedV, DerivedF, DerivedFF, DerivedC> &data);
 
 
-  // Given the current estimate of the field, performes one round of optimization
+  // Given the current estimate of the field, performs one round of optimization
   // iterations and updates the current estimate. The intermediate data is saved
   // and returned for the next iteration.
   // Inputs:

include/igl/is_boundary_edge.cpp

@@ -47,14 +47,14 @@ void igl::is_boundary_edge(
   MatrixXi uE;
   VectorXi EMAP;
   unique_rows(sorted_EallE,uE,_,EMAP);
-  // Counts of occurances
+  // Counts of occurrences
   VectorXi N = VectorXi::Zero(uE.rows());
   for(int e = 0;e<EMAP.rows();e++)
   {
     N(EMAP(e))++;
   }
   B.resize(E.rows());
-  // Look of occurances of 2: one for original and another for boundary
+  // Look of occurrences of 2: one for original and another for boundary
   for(int e = 0;e<E.rows();e++)
   {
     B(e) = (N(EMAP(e)) == 2);
@@ -97,14 +97,14 @@ void igl::is_boundary_edge(
   sort(allE,2,true,sorted_allE,_);
   // Determine unique undirected edges E and map to directed edges EMAP
   unique_rows(sorted_allE,E,_,EMAP);
-  // Counts of occurances
+  // Counts of occurrences
   VectorXi N = VectorXi::Zero(E.rows());
   for(int e = 0;e<EMAP.rows();e++)
   {
     N(EMAP(e))++;
   }
   B.resize(E.rows());
-  // Look of occurances of 1
+  // Look of occurrences of 1
   for(int e = 0;e<E.rows();e++)
   {
     B(e) = N(e) == 1;

include/igl/is_planar.h

@@ -13,7 +13,7 @@
 
 namespace igl
 {
-  // Determin if a set of points lies on the XY plane
+  // Determine if a set of points lies on the XY plane
   //
   // Inputs:
   //   V  #V by dim list of vertex positions

include/igl/lim/lim.h

@@ -59,7 +59,7 @@ namespace igl
     //--------------------------------------------------------------------------
     // Return values:
     //  Succeeded : Successful optimization with fulfilled tolerance
-    //  LocalMinima : Convergenged to a local minima / tolerance not fullfilled
+    //  LocalMinima : Convergenged to a local minima / tolerance not fulfilled
     //  IterationLimit : Max iteration reached before tolerance was fulfilled
     //  Infeasible : not feasible -> has inverted elements (decrease eps?)
   

include/igl/matlab/validate_arg.h

@@ -16,9 +16,9 @@ namespace igl
     // Throw an error if arg i+1 is not a scalar
     //
     // Inputs:
-    //   i  index of current arguement
+    //   i  index of current argument
     //   nrhs  total number of arguments
-    //   prhs  pointer to arguements array
+    //   prhs  pointer to arguments array
     //   name   name of current argument
     IGL_INLINE void validate_arg_scalar(
       const int i, const int nrhs, const mxArray * prhs[], const char * name);

include/igl/matlab_format.h

@@ -22,9 +22,9 @@ namespace igl
   // Templates:
   //   DerivedM  e.g. derived from MatrixXd
   // Input:
-  //   input  some matrix to be formated
+  //   input  some matrix to be formatted
   //   name  name of matrix
-  // Returns  Formated matrix
+  // Returns  Formatted matrix
   //
   // Example:
   // // M := [1 2 3;4 5 6];
@@ -44,7 +44,7 @@ namespace igl
   IGL_INLINE const Eigen::WithFormat< DerivedM > matlab_format(
     const Eigen::DenseBase<DerivedM> & M,
     const std::string name = "");
-  // Same but for sparse matrices. Print IJV format into an auxillary variable
+  // Same but for sparse matrices. Print IJV format into an auxiliary variable
   // and then print a call to sparse which will construct the sparse matrix
   // Example:
   // // S := [0 2 3;4 5 0];

include/igl/min_quad_with_fixed.cpp

@@ -204,7 +204,7 @@ IGL_INLINE bool igl::min_quad_with_fixed_precompute(
       data.preY.resize(data.unknown_lagrange.size(),0);
     }
 
-    // Positive definite and no equality constraints (Postive definiteness
+    // Positive definite and no equality constraints (Positive definiteness
     // implies symmetric)
 #ifdef MIN_QUAD_WITH_FIXED_CPP_DEBUG
     cout<<"    factorize"<<endl;

include/igl/mosek/bbw.h

@@ -30,7 +30,7 @@ namespace igl
     //   Ele  #Elements by simplex-size list of element indices
     //   b  #b boundary indices into V
     //   bc #b by #W list of boundary values
-    //   data  object containing options, intial guess --> solution and results
+    //   data  object containing options, initial guess --> solution and results
     //   mosek_data  object containing mosek options
     // Outputs:
     //   W  #V by #W list of *unnormalized* weights to normalize use

include/igl/mosek/mosek_quadprog.cpp

@@ -31,7 +31,7 @@ igl::mosek::MosekData::MosekData()
 #if MSK_VERSION_MAJOR >= 8
   douparam[MSK_DPAR_INTPNT_QO_TOL_REL_GAP]=1e-12;
 #endif
-  // Force using multiple threads, not sure if MOSEK is properly destorying
+  // Force using multiple threads, not sure if MOSEK is properly destroying
   //extra threads...
 #if MSK_VERSION_MAJOR >= 7
   intparam[MSK_IPAR_NUM_THREADS] = 6;
@@ -53,7 +53,7 @@ igl::mosek::MosekData::MosekData()
   //   choose the right ordering method when really any of them are
   //   instantaneous
   intparam[MSK_IPAR_INTPNT_ORDER_METHOD] = MSK_ORDER_METHOD_NONE;
-  // 1.0 means optimizer is very leniant about declaring model infeasible
+  // 1.0 means optimizer is very lenient about declaring model infeasible
   douparam[MSK_DPAR_INTPNT_TOL_INFEAS] = 1e-8;
   // Hard to say if this is doing anything, probably nothing dramatic
   douparam[MSK_DPAR_INTPNT_TOL_PSAFE]= 1e2;

include/igl/opengl/glfw/background_window.cpp

@@ -25,7 +25,7 @@ IGL_INLINE bool igl::opengl::glfw::background_window(GLFWwindow* & window)
       /* Problem: glewInit failed, something is seriously wrong. */
      fprintf(stderr, "Error: %s\n", glewGetErrorString(err));
     }
-    glGetError(); // pull and savely ignonre unhandled errors like GL_INVALID_ENUM
+    glGetError(); // pull and safely ignore unhandled errors like GL_INVALID_ENUM
   #endif
   return true;
 }

include/igl/opengl2/MouseController.h

@@ -24,7 +24,7 @@ namespace igl
     {
       public:
         typedef Eigen::VectorXi VectorXb;
-        // Propogate selection to descendants so that selected bones and their
+        // Propagate selection to descendants so that selected bones and their
         // subtrees are all selected.
         //
         // Input:
@@ -495,7 +495,7 @@ inline void igl::opengl2::MouseController::set_selection(
   // Combine upward, group rigid parts, repeat
   while(true)
   {
-    // Spread selection accross rigid pieces
+    // Spread selection across rigid pieces
     VectorXb SRP(VectorXb::Zero(RP.maxCoeff()+1));
     for(int e = 0;e<BE.rows();e++)
     {

include/igl/opengl2/lens_flare.cpp

@@ -115,7 +115,7 @@ IGL_INLINE void igl::opengl2::lens_flare_draw(
   using namespace Eigen;
   using namespace std;
 
-  //// view_dir  direction from eye to position is is looking at
+  //// view_dir  direction from eye to position it is looking at
   //const Vector3f view_dir =  (at - from).normalized();
 
   //// near_clip  distance from eye to near clipping plane along view_dir

include/igl/opengl2/lens_flare.h

@@ -21,7 +21,7 @@ namespace igl
     struct Flare{
       int type;             /* flare texture index, 0..5 */
       float scale;
-      float loc;            /* postion on axis */
+      float loc;            /* position on axis */
       float color[3];
       Flare():
         type(-1),

include/igl/opengl2/unproject_to_zero_plane.h

@@ -14,7 +14,7 @@ namespace igl
   namespace opengl2
   {
   // Wrapper for gluUnproject that uses the current GL_MODELVIEW_MATRIX,
-    // GL_PROJECTION_MATRIX, and GL_VIEWPORT to unproject a screen postion
+    // GL_PROJECTION_MATRIX, and GL_VIEWPORT to unproject a screen position
     // (winX,winY) to a 3d location at same depth as the current origin.
     // Inputs:
     //   win*  screen space x, y, and z coordinates respectively

include/igl/parallel_for.h

@@ -31,7 +31,7 @@ namespace igl
   //
   // Inputs:
   //   loop_size  number of iterations. I.e. for(int i = 0;i<loop_size;i++) ...
-  //   func  function handle taking iteration index as only arguement to compute
+  //   func  function handle taking iteration index as only argument to compute
   //     inner block of for loop I.e. for(int i ...){ func(i); }
   //   min_parallel  min size of loop_size such that parallel (non-serial)
   //     thread pooling should be attempted {0}
@@ -65,7 +65,7 @@ namespace igl
   //   prep_func function handle taking n >= number of threads as only
   //     argument 
   //   func  function handle taking iteration index i and thread id t as only
-  //     arguements to compute inner block of for loop I.e. 
+  //     arguments to compute inner block of for loop I.e. 
   //     for(int i ...){ func(i,t); }
   //   accum_func  function handle taking thread index as only argument, to be
   //     called after all calls of func, e.g., for serial accumulation across

include/igl/pathinfo.cpp

@@ -9,7 +9,7 @@
 
 #include "dirname.h"
 #include "basename.h"
-// Verbose should be removed once everythings working correctly
+// Verbose should be removed once everything working correctly
 #include "verbose.h"
 #include <algorithm>
 

include/igl/per_edge_normals.cpp

@@ -35,7 +35,7 @@ IGL_INLINE void igl::per_edge_normals(
   assert(F.cols() == 3 && "Faces must be triangles");
   // number of faces
   const int m = F.rows();
-  // All occurances of directed edges
+  // All occurrences of directed edges
   MatrixXi allE;
   oriented_facets(F,allE);
   // Find unique undirected edges and mapping

include/igl/ply.h.REMOVED.git-id

@@ -1 +1 @@
-e4826d12221d5fca54d3cc95e6a30cbb50fc7fe9
+75a3bab471417c278eec124122f42d7a6f8cfee2

include/igl/png/render_to_png.h

@@ -22,7 +22,7 @@ namespace igl
     //   height height of scene and resulting image
     //   alpha  whether to include alpha channel
     //   fast  sacrifice compression ratio for speed
-    // Returns true only if no errors occured
+    // Returns true only if no errors occurred
     //
     // See also: igl/render_to_tga which is faster but writes .tga files
     IGL_INLINE bool render_to_png(

include/igl/png/render_to_png_async.h

@@ -26,7 +26,7 @@ namespace igl
     //   height height of scene and resulting image
     //   alpha  whether to include alpha channel
     //   fast  sacrifice compression ratio for speed
-    // Returns true only if no errors occured
+    // Returns true only if no errors occurred
     //
     // See also: igl/render_to_tga which is faster but writes .tga files
     IGL_INLINE std::thread render_to_png_async(

include/igl/png/writePNG.cpp

@@ -21,7 +21,7 @@ IGL_INLINE bool igl::png::writePNG(
   assert((R.cols() == G.cols()) && (G.cols() == B.cols()) && (B.cols() == A.cols()));
 
   const int comp = 4;                                  // 4 Channels Red, Green, Blue, Alpha
-  const int stride_in_bytes = R.rows()*comp;           // Lenght of one row in bytes
+  const int stride_in_bytes = R.rows()*comp;           // Length of one row in bytes
   std::vector<unsigned char> data(R.size()*comp,0);     // The image itself;
 
   for (unsigned i = 0; i<R.rows();++i)

include/igl/polyvector_field_matchings.cpp

@@ -26,7 +26,7 @@ IGL_INLINE void igl::polyvector_field_matching(
   // make sure the matching preserve ccw order of the vectors across the edge
   // 1) order vectors in a, ccw  e.g. (0,1,2,3)_a not ccw --> (0,3,2,1)_a ccw
   // 2) order vectors in b, ccw  e.g. (0,1,2,3)_b not ccw --> (0,2,1,3)_b ccw
-  // 3) the vectors in b that match the ordered vectors in a (in this case  (0,3,2,1)_a ) must be a circular shift of the ccw ordered vectors in b  - so we need to explicitely check only these matchings to find the best ccw one, there are N of them
+  // 3) the vectors in b that match the ordered vectors in a (in this case  (0,3,2,1)_a ) must be a circular shift of the ccw ordered vectors in b  - so we need to explicitly check only these matchings to find the best ccw one, there are N of them
   int hN = _ua.cols()/3;
   int N;
   if (is_symmetric)

include/igl/principal_curvature.h

@@ -27,7 +27,7 @@ namespace igl
   // Inputs:
   //   V       eigen matrix #V by 3
   //   F       #F by 3 list of mesh faces (must be triangles)
-  //   radius  controls the size of the neighbourhood used, 1 = average edge lenght
+  //   radius  controls the size of the neighbourhood used, 1 = average edge length
   //
   // Outputs:
   //   PD1 #V by 3 maximal curvature direction for each vertex.

include/igl/project_to_line.h

@@ -12,7 +12,7 @@
 
 namespace igl
 {
-  // PROJECT_TO_LINE  project points onto vectors, that is find the paramter
+  // PROJECT_TO_LINE  project points onto vectors, that is find the parameter
   // t for a point p such that proj_p = (y-x).*t, additionally compute the
   // squared distance from p to the line of the vector, such that 
   // |p - proj_p|² = sqr_d

include/igl/project_to_line_segment.h

@@ -12,7 +12,7 @@
 
 namespace igl
 {
-  // PROJECT_TO_LINE_SEGMENT project points onto vectors, that is find the paramter
+  // PROJECT_TO_LINE_SEGMENT project points onto vectors, that is find the parameter
   // t for a point p such that proj_p = (y-x).*t, additionally compute the
   // squared distance from p to the line of the vector, such that 
   // |p - proj_p|² = sqr_d

include/igl/raytri.c

@@ -129,7 +129,7 @@ inline int intersect_triangle1(double orig[3], double dir[3],
       if (*v > 0.0 || *u + *v < det)
 	 return 0;
    }
-   else return 0;  /* ray is parallell to the plane of the triangle */
+   else return 0;  /* ray is parallel to the plane of the triangle */
 
 
    inv_det = 1.0 / det;
@@ -196,7 +196,7 @@ inline int intersect_triangle2(double orig[3], double dir[3],
       if (*v > 0.0 || *u + *v < det)
 	 return 0;
    }
-   else return 0;  /* ray is parallell to the plane of the triangle */
+   else return 0;  /* ray is parallel to the plane of the triangle */
 
    /* calculate t, ray intersects triangle */
    *t = IGL_RAY_TRI_DOT(edge2, qvec) * inv_det;
@@ -256,7 +256,7 @@ inline int intersect_triangle3(double orig[3], double dir[3],
       if (*v > 0.0 || *u + *v < det)
 	 return 0;
    }
-   else return 0;  /* ray is parallell to the plane of the triangle */
+   else return 0;  /* ray is parallel to the plane of the triangle */
 
    *t = IGL_RAY_TRI_DOT(edge2, qvec) * inv_det;
    (*u) *= inv_det;

include/igl/readOBJ.cpp

@@ -48,7 +48,7 @@ IGL_INLINE bool igl::readOBJ(
   std::vector<std::vector<Index > > & FTC,
   std::vector<std::vector<Index > > & FN)
 {
-  // File open was succesfull so clear outputs
+  // File open was successful so clear outputs
   V.clear();
   TC.clear();
   N.clear();
@@ -56,7 +56,7 @@ IGL_INLINE bool igl::readOBJ(
   FTC.clear();
   FN.clear();
 
-  // variables an constants to assist parsing the .obj file
+  // variables and constants to assist parsing the .obj file
   // Constant strings to compare against
   std::string v("v");
   std::string vn("vn");

include/igl/readTGF.h

@@ -25,7 +25,7 @@ namespace igl
   //
   // Input:
   //  filename  .tgf file name
-  // Ouput:
+  // Output:
   //  V  # vertices by 3 list of vertex positions
   //  E  # edges by 2 list of edge indices
   //  P  # point-handles list of point handle indices

include/igl/setdiff.cpp

@@ -42,7 +42,7 @@ IGL_INLINE void igl::setdiff(
       }
     }
     assert(k == C.size());
-    // Have to use << instead of = becasue Eigen's PlainObjectBase is annoying
+    // Have to use << instead of = because Eigen's PlainObjectBase is annoying
     IA << igl::LinSpaced<Eigen::Matrix<typename DerivedIA::Scalar,Eigen::Dynamic,1> >(
       C.size(),0,C.size()-1);
   }

include/igl/slice.cpp

@@ -354,6 +354,7 @@ template void igl::slice<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<in
 template Eigen::Matrix<double, -1, -1, 0, -1, -1> igl::slice<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::DenseBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
 template void igl::slice<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, int, Eigen::Matrix<double, -1, 3, 0, -1, 3>&);
 template void igl::slice<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::DenseBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::DenseBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+template void igl::slice<unsigned int, unsigned int>(Eigen::SparseMatrix<unsigned int, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<unsigned int, 0, int>&);
 #ifdef WIN32
 template void igl::slice<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>>>(class Eigen::Matrix<__int64, -1, 1, 0, -1, 1> const &, class Eigen::DenseBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> &);
 template void igl::slice<class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>>, class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>>>(class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> const &, class Eigen::DenseBase<class Eigen::Matrix<__int64, -1, 1, 0, -1, 1>> const &, int, class Eigen::PlainObjectBase<class Eigen::Matrix<int, -1, -1, 0, -1, -1>> &);

include/igl/slice_cached.cpp

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@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 "slice_cached.h"
+
+#include <iostream>
+#include <vector>
+#include <utility>
+#include <igl/slice.h>
+
+template <typename TX, typename TY>
+IGL_INLINE void igl::slice_cached_precompute(
+  const Eigen::SparseMatrix<TX>& X,
+  const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+  const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
+  Eigen::SparseMatrix<TY>& Y,
+  Eigen::VectorXi& data)
+{
+  // Create a sparse matrix whose entries are the ids
+  Eigen::SparseMatrix<unsigned> TS = X.template cast<unsigned>();
+
+  TS.makeCompressed();
+  for (unsigned i=0;i<TS.nonZeros();++i)
+    *(TS.valuePtr() + i) = i;
+
+  Eigen::SparseMatrix<unsigned> TS_sliced;
+  igl::slice(TS,R,C,TS_sliced);
+  Y = TS_sliced.cast<TY>();
+
+  data.resize(TS_sliced.nonZeros());
+  for (unsigned i=0;i<data.size();++i)
+  {
+    data[i] = *(TS_sliced.valuePtr() + i);
+    *(Y.valuePtr() + i) = *(X.valuePtr() + data[i]);
+  }
+}
+
+template <typename TX, typename TY>
+IGL_INLINE void igl::slice_cached(
+  const Eigen::SparseMatrix<TX>& X,
+  Eigen::SparseMatrix<TY>& Y,
+  const Eigen::VectorXi& data)
+{
+  for (unsigned i=0; i<data.size(); ++i)
+    *(Y.valuePtr() + i) = *(X.valuePtr() + data[i]);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::slice_cached_precompute<double, double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template void igl::slice_cached<double, double>(Eigen::SparseMatrix<double, 0, int> const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
+#endif

include/igl/slice_cached.h

@@ -0,0 +1,66 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@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_SLICE_CACHED_H
+#define IGL_SLICE_CACHED_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{  
+
+  // Act like the matlab X(row_indices,col_indices) operator, where
+  // row_indices, col_indices are non-negative integer indices. This is a fast version
+  // of igl::slice that can analyze and store the sparsity structure. It is slower at the // first evaluation (slice_cached_precompute), but faster on the subsequent ones.
+  // 
+  // Inputs:
+  //   X  m by n matrix
+  //   R  list of row indices
+  //   C  list of column indices
+  //   
+  // Output:
+  //   Y  #R by #C matrix
+  //   data Temporary data used by slice_cached to repeat this operation
+  //
+  // Usage:
+  //
+  // // Construct and slice up Laplacian
+  // SparseMatrix<double> L,L_sliced;
+  // igl::cotmatrix(V,F,L);
+
+  // // Normal igl::slice call
+  // igl::slice(L,in,in,L_in_in);
+
+  // // Fast version
+  // static VectorXi data; // static or saved in a global state
+  // if (data.size() == 0)
+  //   igl::slice_cached_precompute(L,in,in,L_sliced,data);
+  // else
+  //   igl::slice_cached(L,L_sliced,temp);
+
+  template <typename TX, typename TY>
+  IGL_INLINE void slice_cached_precompute(
+    const Eigen::SparseMatrix<TX>& X,
+    const Eigen::Matrix<int,Eigen::Dynamic,1> & R,
+    const Eigen::Matrix<int,Eigen::Dynamic,1> & C,
+    Eigen::SparseMatrix<TY>& Y,
+    Eigen::VectorXi& data);
+
+  template <typename TX, typename TY>
+  IGL_INLINE void slice_cached(
+    const Eigen::SparseMatrix<TX>& X,
+    Eigen::SparseMatrix<TY>& Y,
+    const Eigen::VectorXi& data);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "slice_cached.cpp"
+#endif
+
+#endif

include/igl/slim.cpp

@@ -34,6 +34,14 @@
 #include <Eigen/SparseCholesky>
 #include <Eigen/IterativeLinearSolvers>
 
+#include <igl/Timer.h>
+#include <igl/sparse_cached.h>
+#include <igl/AtA_cached.h>
+
+#ifdef CHOLMOD
+#include <Eigen/CholmodSupport>
+#endif
+
 namespace igl
 {
   namespace slim
@@ -42,8 +50,8 @@ namespace igl
     IGL_INLINE void compute_surface_gradient_matrix(const Eigen::MatrixXd &V, const Eigen::MatrixXi &F,
                                                     const Eigen::MatrixXd &F1, const Eigen::MatrixXd &F2,
                                                     Eigen::SparseMatrix<double> &D1, Eigen::SparseMatrix<double> &D2);
-    IGL_INLINE void buildA(igl::SLIMData& s, Eigen::SparseMatrix<double> &A);
-    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At);
+    IGL_INLINE void buildA(igl::SLIMData& s, std::vector<Eigen::Triplet<double> > & IJV);
+    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &A);
     IGL_INLINE void add_soft_constraints(igl::SLIMData& s, Eigen::SparseMatrix<double> &L);
     IGL_INLINE double compute_energy(igl::SLIMData& s, Eigen::MatrixXd &V_new);
     IGL_INLINE double compute_soft_const_energy(igl::SLIMData& s,
@@ -395,8 +403,12 @@ namespace igl
       Eigen::SparseMatrix<double> L;
       build_linear_system(s,L);
 
+      igl::Timer t;
+      
+      //t.start();
       // solve
       Eigen::VectorXd Uc;
+#ifndef CHOLMOD
       if (s.dim == 2)
       {
         SimplicialLDLT<Eigen::SparseMatrix<double> > solver;
@@ -406,13 +418,21 @@ namespace igl
       { // seems like CG performs much worse for 2D and way better for 3D
         Eigen::VectorXd guess(uv.rows() * s.dim);
         for (int i = 0; i < s.v_num; i++) for (int j = 0; j < s.dim; j++) guess(uv.rows() * j + i) = uv(i, j); // flatten vector
-        ConjugateGradient<Eigen::SparseMatrix<double>, Eigen::Lower | Upper> solver;
-        solver.setTolerance(1e-8);
-        Uc = solver.compute(L).solveWithGuess(s.rhs, guess);
+        ConjugateGradient<Eigen::SparseMatrix<double>, Lower | Upper> cg;
+        cg.setTolerance(1e-8);
+        cg.compute(L);
+        Uc = cg.solveWithGuess(s.rhs, guess);
       }
-
+#else
+        CholmodSimplicialLDLT<Eigen::SparseMatrix<double> > solver;
+        Uc = solver.compute(L).solve(s.rhs);
+#endif
       for (int i = 0; i < s.dim; i++)
         uv.col(i) = Uc.block(i * s.v_n, 0, s.v_n, 1);
+
+      // t.stop();
+      // std::cerr << "solve: " << t.getElapsedTime() << std::endl;
+
     }
 
 
@@ -478,20 +498,60 @@ namespace igl
     IGL_INLINE void build_linear_system(igl::SLIMData& s, Eigen::SparseMatrix<double> &L)
     {
       // formula (35) in paper
+      std::vector<Eigen::Triplet<double> > IJV;
+      
+      #ifdef SLIM_CACHED
+      buildA(s,IJV);
+      if (s.A.rows() == 0)
+      {
+        s.A = Eigen::SparseMatrix<double>(s.dim * s.dim * s.f_n, s.dim * s.v_n);
+        igl::sparse_cached_precompute(IJV,s.A,s.A_data);
+      }
+      else
+        igl::sparse_cached(IJV,s.A,s.A_data);
+      #else
       Eigen::SparseMatrix<double> A(s.dim * s.dim * s.f_n, s.dim * s.v_n);
-      buildA(s,A);
+      buildA(s,IJV);
+      A.setFromTriplets(IJV.begin(),IJV.end());
+      A.makeCompressed();
+      #endif
 
+      #ifdef SLIM_CACHED
+      #else
       Eigen::SparseMatrix<double> At = A.transpose();
       At.makeCompressed();
+      #endif
+
+      #ifdef SLIM_CACHED
+      Eigen::SparseMatrix<double> id_m(s.A.cols(), s.A.cols());
+      #else
+      Eigen::SparseMatrix<double> id_m(A.cols(), A.cols());
+      #endif
 
-      Eigen::SparseMatrix<double> id_m(At.rows(), At.rows());
       id_m.setIdentity();
 
       // add proximal penalty
+      #ifdef SLIM_CACHED
+      s.AtA_data.W = s.WGL_M;
+      if (s.AtA.rows() == 0)
+        igl::AtA_cached_precompute(s.A,s.AtA,s.AtA_data);
+      else
+        igl::AtA_cached(s.A,s.AtA,s.AtA_data);
+
+      L = s.AtA + s.proximal_p * id_m; //add also a proximal 
+      L.makeCompressed();
+
+      #else
       L = At * s.WGL_M.asDiagonal() * A + s.proximal_p * id_m; //add also a proximal term
       L.makeCompressed();
+      #endif
+
+      #ifdef SLIM_CACHED
+      buildRhs(s, s.A);
+      #else
+      buildRhs(s, A);
+      #endif
 
-      buildRhs(s, At);
       Eigen::SparseMatrix<double> OldL = L;
       add_soft_constraints(s,L);
       L.makeCompressed();
@@ -657,10 +717,9 @@ namespace igl
       return energy;
     }
 
-    IGL_INLINE void buildA(igl::SLIMData& s, Eigen::SparseMatrix<double> &A)
+    IGL_INLINE void buildA(igl::SLIMData& s, std::vector<Eigen::Triplet<double> > & IJV)
     {
       // formula (35) in paper
-      std::vector<Eigen::Triplet<double> > IJV;
       if (s.dim == 2)
       {
         IJV.reserve(4 * (s.Dx.outerSize() + s.Dy.outerSize()));
@@ -780,10 +839,9 @@ namespace igl
           }
         }
       }
-      A.setFromTriplets(IJV.begin(), IJV.end());
     }
 
-    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &At)
+    IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &A)
     {
       Eigen::VectorXd f_rhs(s.dim * s.dim * s.f_n);
       f_rhs.setZero();
@@ -830,7 +888,7 @@ namespace igl
         for (int j = 0; j < s.v_n; j++)
           uv_flat(s.v_n * i + j) = s.V_o(j, i);
 
-      s.rhs = (At * s.WGL_M.asDiagonal() * f_rhs + s.proximal_p * uv_flat);
+      s.rhs = (f_rhs.transpose() * s.WGL_M.asDiagonal() * A).transpose() + s.proximal_p * uv_flat;
     }
 
   }

include/igl/slim.h

@@ -12,6 +12,14 @@
 #include <Eigen/Dense>
 #include <Eigen/Sparse>
 
+// This option makes the iterations faster (all except the first) by caching the 
+// sparsity pattern of the matrix involved in the assembly. It should be on if you plan to do many iterations, off if you have to change the matrix structure at every iteration.
+#define SLIM_CACHED 
+
+#ifdef SLIM_CACHED
+#include <igl/AtA_cached.h>
+#endif
+
 namespace igl
 {
 
@@ -64,6 +72,13 @@ struct SLIMData
   bool first_solve;
   bool has_pre_calc = false;
   int dim;
+
+  #ifdef SLIM_CACHED
+  Eigen::SparseMatrix<double> A;
+  Eigen::VectorXi A_data;
+  Eigen::SparseMatrix<double> AtA;
+  igl::AtA_cached_data AtA_data;
+  #endif
 };
 
 // Compute necessary information to start using SLIM

include/igl/sort_angles.cpp

@@ -18,7 +18,7 @@ IGL_INLINE void igl::sort_angles(
     assert(num_cols >= 2);
 
     R.resize(num_rows);
-    // Have to use << instead of = becasue Eigen's PlainObjectBase is annoying
+    // Have to use << instead of = because Eigen's PlainObjectBase is annoying
     R << igl::LinSpaced<
       Eigen::Matrix<typename DerivedR::Scalar, Eigen::Dynamic, 1> >
       (num_rows, 0, num_rows-1);

include/igl/sort_triangles.cpp

@@ -390,7 +390,7 @@ IGL_INLINE void igl::sort_triangles(
 //  assert(false && 
 //    "THIS WILL NEVER WORK because depth sorting is not a numerical sort where"
 //    "pairwise comparisons of triangles are transitive.  Rather it is a"
-//    "topological sort on a dependecy graph. Dependency encodes 'This triangle"
+//    "topological sort on a dependency graph. Dependency encodes 'This triangle"
 //    "must be drawn before that one'");
 //  using namespace std;
 //  using namespace Eigen;

include/igl/sparse_cached.cpp

@@ -0,0 +1,122 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@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 "sparse_cached.h"
+
+#include <iostream>
+#include <vector>
+#include <array>
+#include <unordered_map>
+#include <map>
+#include <utility>
+
+template <typename DerivedI, typename Scalar>
+IGL_INLINE void igl::sparse_cached_precompute(
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedI> & J,
+  Eigen::SparseMatrix<Scalar>& X,
+  Eigen::VectorXi& data)
+{
+  // Generates the triplets
+  std::vector<Eigen::Triplet<double> > t(I.size());
+  for (unsigned i = 0; i<I.size(); ++i)
+    t[i] = Eigen::Triplet<double>(I[i],J[i],1);
+
+  // Call the triplets version
+  sparse_cached_precompute(t,X,data);
+}
+
+template <typename Scalar>
+IGL_INLINE void igl::sparse_cached_precompute(
+  const std::vector<Eigen::Triplet<Scalar> >& triplets,
+  Eigen::SparseMatrix<Scalar>& X,
+  Eigen::VectorXi& data)
+{
+    // Construct an empty sparse matrix
+    X.setFromTriplets(triplets.begin(),triplets.end());
+    X.makeCompressed();
+
+    std::vector<std::array<int,3> > T(triplets.size());
+    for (unsigned i=0; i<triplets.size(); ++i)
+    {
+      T[i][0] = triplets[i].col();
+      T[i][1] = triplets[i].row();
+      T[i][2] = i;
+    }
+
+    std::sort(T.begin(), T.end());
+
+    data.resize(triplets.size());
+
+    int t = 0;
+
+    for (unsigned k=0; k<X.outerSize(); ++k)
+    {
+      unsigned outer_index = *(X.outerIndexPtr()+k);
+      unsigned next_outer_index = (k+1 == X.outerSize()) ? X.nonZeros() : *(X.outerIndexPtr()+k+1); 
+      
+      for (unsigned l=outer_index; l<next_outer_index; ++l)
+      {
+        int col = k;
+        int row = *(X.innerIndexPtr()+l);
+        int value_index = l;
+        assert(col < X.cols());
+        assert(col >= 0);
+        assert(row < X.rows());
+        assert(row >= 0);
+        assert(value_index >= 0);
+        assert(value_index < X.nonZeros());
+
+        std::pair<int,int> p_m = std::make_pair(row,col);
+
+        while (t<T.size() && (p_m == std::make_pair(T[t][1],T[t][0])))
+          data[T[t++][2]] = value_index;
+      }
+    }
+    assert(t==T.size());
+
+}
+  
+template <typename Scalar>
+IGL_INLINE void igl::sparse_cached(
+  const std::vector<Eigen::Triplet<Scalar> >& triplets,
+  Eigen::SparseMatrix<Scalar>& X,
+  const Eigen::VectorXi& data)
+{
+  assert(triplets.size() == data.size());
+
+  // Clear it first
+  for (unsigned i = 0; i<data.size(); ++i)
+    *(X.valuePtr() + data[i]) = 0;
+ 
+  // Then sum them up
+  for (unsigned i = 0; i<triplets.size(); ++i)
+    *(X.valuePtr() + data[i]) += triplets[i].value();
+}
+
+template <typename DerivedV, typename Scalar>
+IGL_INLINE void igl::sparse_cached(
+  const Eigen::MatrixBase<DerivedV>& V,
+  Eigen::SparseMatrix<Scalar>& X,
+  const Eigen::VectorXi& data)
+{
+  assert(V.size() == data.size());
+
+  // Clear it first
+  for (unsigned i = 0; i<data.size(); ++i)
+    *(X.valuePtr() + data[i]) = 0;
+ 
+  // Then sum them up
+  for (unsigned i = 0; i<V.size(); ++i)
+    *(X.valuePtr() + data[i]) += V[i];
+}
+
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::sparse_cached_precompute<double>(std::vector<Eigen::Triplet<double, Eigen::SparseMatrix<double, 0, int>::Index>, std::allocator<Eigen::Triplet<double, Eigen::SparseMatrix<double, 0, int>::Index> > > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1>&);
+template void igl::sparse_cached<double>(std::vector<Eigen::Triplet<double, Eigen::SparseMatrix<double, 0, int>::Index>, std::allocator<Eigen::Triplet<double, Eigen::SparseMatrix<double, 0, int>::Index> > > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&);
+#endif

include/igl/sparse_cached.h

@@ -0,0 +1,82 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 Daniele Panozzo <daniele.panozzo@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_SPARSE_CACHED_H
+#define IGL_SPARSE_CACHED_H
+#include "igl_inline.h"
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+namespace igl
+{
+  // Build a sparse matrix from list of indices and values (I,J,V), similarly to 
+  // the sparse function in matlab. Divides the construction in two phases, one
+  // for fixing the sparsity pattern, and one to populate it with values. Compared to
+  // igl::sparse, this version is slower for the first time (since it requires a
+  // precomputation), but faster to the subsequent evaluations.
+  //
+  // Templates:
+  //   IndexVector  list of indices, value should be non-negative and should
+  //     expect to be cast to an index. Must implement operator(i) to retrieve
+  //     ith element
+  //   ValueVector  list of values, value should be expect to be cast to type
+  //     T. Must implement operator(i) to retrieve ith element
+  //   T  should be a eigen sparse matrix primitive type like int or double
+  // Input:
+  //   I  nnz vector of row indices of non zeros entries in X
+  //   J  nnz vector of column indices of non zeros entries in X
+  //   V  nnz vector of non-zeros entries in X
+  //   Optional:
+  //     m  number of rows
+  //     n  number of cols
+  // Outputs:
+  //   X  m by n matrix of type T whose entries are to be found 
+  //
+  // Example:
+  //   Eigen::SparseMatrix<double> A;
+  //   std::vector<Eigen::Triplet<double> > IJV;
+  //   buildA(IJV);
+  //   if (A.rows() == 0)
+  //   {
+  //     A = Eigen::SparseMatrix<double>(rows,cols);
+  //     igl::sparse_cached_precompute(IJV,A,A_data);
+  //   }
+  //   else
+  //     igl::sparse_cached(IJV,s.A,s.A_data);
+
+  template <typename DerivedI, typename Scalar>
+  IGL_INLINE void sparse_cached_precompute(
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedI> & J,
+    Eigen::SparseMatrix<Scalar>& X,
+    Eigen::VectorXi& data);
+  
+  template <typename Scalar>
+  IGL_INLINE void sparse_cached_precompute(
+    const std::vector<Eigen::Triplet<Scalar> >& triplets,
+    Eigen::SparseMatrix<Scalar>& X,
+    Eigen::VectorXi& data);
+
+  template <typename Scalar>
+  IGL_INLINE void sparse_cached(
+    const std::vector<Eigen::Triplet<Scalar> >& triplets,
+    Eigen::SparseMatrix<Scalar>& X,
+    const Eigen::VectorXi& data);
+
+  template <typename DerivedV, typename Scalar>
+  IGL_INLINE void sparse_cached(
+    const Eigen::MatrixBase<DerivedV>& V,
+    Eigen::SparseMatrix<Scalar>& X,
+    const Eigen::VectorXi& data);
+  
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "sparse_cached.cpp"
+#endif
+
+#endif

include/igl/triangle_fan.h

@@ -11,13 +11,13 @@
 #include <Eigen/Dense>
 namespace igl
 {
-  // Given a list of faces tesselate all of the "exterior" edges forming another
+  // Given a list of faces tessellate all of the "exterior" edges forming another
   // list of 
   //
   // Inputs:
   //   E  #E by simplex_size-1  list of exterior edges (see exterior_edges.h)
   // Outputs:
-  //   cap  #cap by simplex_size  list of "faces" tessleating the boundary edges
+  //   cap  #cap by simplex_size  list of "faces" tessellating the boundary edges
   IGL_INLINE void triangle_fan(
     const Eigen::MatrixXi & E,
     Eigen::MatrixXi & cap);

include/igl/unique_edge_map.cpp

@@ -25,7 +25,7 @@ IGL_INLINE void igl::unique_edge_map(
 {
   using namespace Eigen;
   using namespace std;
-  // All occurances of directed edges
+  // All occurrences of directed edges
   oriented_facets(F,E);
   const size_t ne = E.rows();
   // This is 2x faster to create than a map from pairs to lists of edges and 5x

include/igl/viewer/OpenGL_state.cpp

@@ -72,7 +72,7 @@ IGL_INLINE void igl::viewer::OpenGL_state::set_data(const igl::viewer::ViewerDat
 
   if (!data.face_based)
   {
-    if (!per_corner_uv)
+    if (!(per_corner_uv || per_corner_normals))
     {
       // Vertex positions
       if (dirty & ViewerData::DIRTY_POSITION)
@@ -115,7 +115,7 @@ IGL_INLINE void igl::viewer::OpenGL_state::set_data(const igl::viewer::ViewerDat
 
       if (dirty & ViewerData::DIRTY_AMBIENT)
       {
-        V_ambient_vbo.resize(3,data.F.rows()*3);
+        V_ambient_vbo.resize(4,data.F.rows()*3);
         for (unsigned i=0; i<data.F.rows();++i)
           for (unsigned j=0;j<3;++j)
             V_ambient_vbo.col (i*3+j) = data.V_material_ambient.row(data.F(i,j)).transpose().cast<float>();
@@ -123,7 +123,7 @@ IGL_INLINE void igl::viewer::OpenGL_state::set_data(const igl::viewer::ViewerDat
 
       if (dirty & ViewerData::DIRTY_DIFFUSE)
       {
-        V_diffuse_vbo.resize(3,data.F.rows()*3);
+        V_diffuse_vbo.resize(4,data.F.rows()*3);
         for (unsigned i=0; i<data.F.rows();++i)
           for (unsigned j=0;j<3;++j)
             V_diffuse_vbo.col (i*3+j) = data.V_material_diffuse.row(data.F(i,j)).transpose().cast<float>();
@@ -131,7 +131,7 @@ IGL_INLINE void igl::viewer::OpenGL_state::set_data(const igl::viewer::ViewerDat
 
       if (dirty & ViewerData::DIRTY_SPECULAR)
       {
-        V_specular_vbo.resize(3,data.F.rows()*3);
+        V_specular_vbo.resize(4,data.F.rows()*3);
         for (unsigned i=0; i<data.F.rows();++i)
           for (unsigned j=0;j<3;++j)
             V_specular_vbo.col(i*3+j) = data.V_material_specular.row(data.F(i,j)).transpose().cast<float>();
@@ -142,7 +142,12 @@ IGL_INLINE void igl::viewer::OpenGL_state::set_data(const igl::viewer::ViewerDat
         V_normals_vbo.resize(3,data.F.rows()*3);
         for (unsigned i=0; i<data.F.rows();++i)
           for (unsigned j=0;j<3;++j)
-            V_normals_vbo.col (i*3+j) = data.V_normals.row(data.F(i,j)).transpose().cast<float>();
+          
+            V_normals_vbo.col (i*3+j) = 
+                         per_corner_normals ?
+               data.F_normals.row(i*3+j).transpose().cast<float>() :
+               data.V_normals.row(data.F(i,j)).transpose().cast<float>();
+
 
         if (invert_normals)
           V_normals_vbo = -V_normals_vbo;

include/igl/viewer/Viewer.cpp

@@ -937,7 +937,7 @@ namespace viewer
         /* Problem: glewInit failed, something is seriously wrong. */
        fprintf(stderr, "Error: %s\n", glewGetErrorString(err));
       }
-      glGetError(); // pull and savely ignonre unhandled errors like GL_INVALID_ENUM
+      glGetError(); // pull and safely ignore unhandled errors like GL_INVALID_ENUM
       fprintf(stdout, "Status: Using GLEW %s\n", glewGetString(GLEW_VERSION));
     #endif
 
@@ -946,7 +946,7 @@ namespace viewer
       major = glfwGetWindowAttrib(window, GLFW_CONTEXT_VERSION_MAJOR);
       minor = glfwGetWindowAttrib(window, GLFW_CONTEXT_VERSION_MINOR);
       rev = glfwGetWindowAttrib(window, GLFW_CONTEXT_REVISION);
-      printf("OpenGL version recieved: %d.%d.%d\n", major, minor, rev);
+      printf("OpenGL version received: %d.%d.%d\n", major, minor, rev);
       printf("Supported OpenGL is %s\n", (const char*)glGetString(GL_VERSION));
       printf("Supported GLSL is %s\n", (const char*)glGetString(GL_SHADING_LANGUAGE_VERSION));
     #endif

include/igl/viewer/ViewerCore.h

@@ -119,7 +119,7 @@ public:
   float model_zoom;
   Eigen::Vector3f model_translation;
 
-  // Model viewing paramters (uv coordinates)
+  // Model viewing parameters (uv coordinates)
   float model_zoom_uv;
   Eigen::Vector3f model_translation_uv;
 

include/igl/viewer/ViewerData.h

@@ -44,7 +44,7 @@ public:
     DIRTY_ALL            = 0x03FF
   };
 
-  // Empy all fields
+  // Empty all fields
   IGL_INLINE void clear();
 
   // Change the visualization mode, invalidating the cache if necessary

include/igl/xml/writeDAE.cpp

@@ -25,7 +25,7 @@ IGL_INLINE bool igl::xml::writeDAE(
 
   const auto & ele = [&doc](
     const std::string tag,
-    // Can't just use `{}` as the default arguement because of a clang-bug
+    // Can't just use `{}` as the default argument because of a clang-bug
     // http://stackoverflow.com/questions/17264067/
     const std::map<std::string,std::string> attribs =
       std::map<std::string,std::string>(),

optional/README.md

@@ -232,7 +232,7 @@ Alternatively, you can also compress everything into a single header file:
 
 * **Hard to debug/edit**: The compressed files are
   automatically generated. They're huge and should not be edited. Thus
-  debugging and editting are near impossible.
+  debugging and editing are near impossible.
 
 * **Compounded dependencies**:
   An immediate disadvantage of this
@@ -241,7 +241,7 @@ Alternatively, you can also compress everything into a single header file:
   `igl.cpp` will require linking to all of `libigl`'s
   dependencies (`OpenGL`, `GLUT`,
   `AntTweakBar`, `BLAS`). However, because all
-  depencies other than Eigen should be encapsulated between
+  dependencies other than Eigen should be encapsulated between
   `#ifndef` guards (e.g. `#ifndef IGL_NO_OPENGL`, it
   is possible to ignore certain functions that have such dependencies.
 

optional/index.html

@@ -283,7 +283,7 @@ exposing templated functions.</li>
 <ul>
 <li><p><strong>Hard to debug/edit</strong>: The compressed files are
  automatically generated. They&#8217;re huge and should not be edited. Thus
- debugging and editting are near impossible.</p></li>
+ debugging and editing are near impossible.</p></li>
 <li><p><strong>Compounded dependencies</strong>:
  An immediate disadvantage of this
  seems to be that even to use a single function (e.g.
@@ -291,7 +291,7 @@ exposing templated functions.</li>
  <code>igl.cpp</code> will require linking to all of <code>libigl</code>&#8217;s
  dependencies (<code>OpenGL</code>, <code>GLUT</code>,
  <code>AntTweakBar</code>, <code>BLAS</code>). However, because all
- depencies other than Eigen should be encapsulated between
+ dependencies other than Eigen should be encapsulated between
  <code>#ifndef</code> guards (e.g. <code>#ifndef IGL_NO_OPENGL</code>, it
  is possible to ignore certain functions that have such dependencies.</p></li>
 <li><p><strong>Long compile</strong>:

python/README.md

@@ -122,11 +122,11 @@ viewer.launch()
 
 ### Remote viewer
 
-Whe using the viewer from an interactive python shell (iPython), it is
+When using the viewer from an interactive python shell (iPython), it is
 inconvenient to let the viewer take control of the main thread for rendering
 purposes. We provide a simple wrapper for the viewer that allows to launch
 a remote process and send meshes to it via a TCP/IP socket. For more
-informations on how to use it see the documentation in [tcpviewer.py](tcpviewer.py)
+information on how to use it see the documentation in [tcpviewer.py](tcpviewer.py)
 
 ## Matlab
 

python/modules/py_vector.cpp

@@ -100,7 +100,7 @@ py::class_<Type> bind_eigen_2(py::module &m, const char *name) {
         .def("rows", [](const Type &m) { return m.rows(); })
         .def("shape", [](const Type &m) { return std::tuple<int,int>(m.rows(), m.cols()); })
 
-        /* Extract rows and colums */
+        /* Extract rows and columns */
         .def("col", [](const Type &m, int i) {
             if (i<0 || i>=m.cols())
               throw std::runtime_error("Column index out of bound.");

python/py_doc.cpp

@@ -126,7 +126,7 @@ const char *__doc_igl_bbw = R"igl_Qu8mg5v7(// Compute Bounded Biharmonic Weights
   //   Ele  #Elements by simplex-size list of element indices
   //   b  #b boundary indices into V
   //   bc #b by #W list of boundary values
-  //   data  object containing options, intial guess --> solution and results
+  //   data  object containing options, initial guess --> solution and results
   // Outputs:
   //   W  #V by #W list of *unnormalized* weights to normalize use
   //    igl::normalize_row_sums(W,W);
@@ -182,7 +182,7 @@ const char *__doc_igl_cat = R"igl_Qu8mg5v7(// Perform concatenation of a two mat
   // Template:
   //   Scalar  scalar data type for sparse matrices like double or int
   //   Mat  matrix type for all matrices (e.g. MatrixXd, SparseMatrix)
-  //   MatC  matrix type for ouput matrix (e.g. MatrixXd) needs to support
+  //   MatC  matrix type for output matrix (e.g. MatrixXd) needs to support
   //     resize
   // Inputs:
   //   A  first input matrix
@@ -322,7 +322,7 @@ const char *__doc_igl_copyleft_comiso_miq = R"igl_Qu8mg5v7(// Inputs:
     //   UV             #UV by 2 list of vertices in 2D
     //   FUV            #FUV by 3 list of face indices in UV
     //
-    // TODO: rename the parameters name in the cpp consistenly
+    // TODO: rename the parameters name in the cpp consistently
     //       improve the handling of hard_features, right now it might fail in difficult cases)igl_Qu8mg5v7";
 const char *__doc_igl_copyleft_comiso_nrosy = R"igl_Qu8mg5v7(// Generate a N-RoSy field from a sparse set of constraints
     //
@@ -336,7 +336,7 @@ const char *__doc_igl_copyleft_comiso_nrosy = R"igl_Qu8mg5v7(// Generate a N-RoS
     //   w_soft  #S by 1 weight for the soft constraints (0-1)
     //   bc_soft #S by 3 bc for soft constraints
     //   N       the degree of the N-RoSy vector field
-    //   soft    the strenght of the soft contraints w.r.t. smoothness
+    //   soft    the strength of the soft constraints w.r.t. smoothness
     //           (0 -> smoothness only, 1->constraints only)
     // Outputs:
     //   R       #F by 3 the representative vectors of the interpolated field
@@ -473,7 +473,7 @@ const char *__doc_igl_directed_edge_orientations = R"igl_Qu8mg5v7(// Determine r
   // Outputs:
   //   Q  #E list of quaternions 
   //)igl_Qu8mg5v7";
-const char *__doc_igl_directed_edge_parents = R"igl_Qu8mg5v7(// Recover "parents" (preceeding edges) in a tree given just directed edges.
+const char *__doc_igl_directed_edge_parents = R"igl_Qu8mg5v7(// Recover "parents" (preceding edges) in a tree given just directed edges.
   //
   // Inputs:
   //   E  #E by 2 list of directed edges
@@ -1040,7 +1040,7 @@ const char *__doc_igl_principal_curvature = R"igl_Qu8mg5v7(// Compute the princi
   // Inputs:
   //   V       eigen matrix #V by 3
   //   F       #F by 3 list of mesh faces (must be triangles)
-  //   radius  controls the size of the neighbourhood used, 1 = average edge lenght
+  //   radius  controls the size of the neighbourhood used, 1 = average edge length
   //
   // Outputs:
   //   PD1 #V by 3 maximal curvature direction for each vertex.
@@ -1123,7 +1123,7 @@ const char *__doc_igl_readTGF = R"igl_Qu8mg5v7(// READTGF
   //
   // Input:
   //  filename  .tgf file name
-  // Ouput:
+  // Output:
   //  V  # vertices by 3 list of vertex positions
   //  E  # edges by 2 list of edge indices
   //  P  # point-handles list of point handle indices

python/tutorial/301_Slice.py

@@ -22,11 +22,11 @@ F = igl.eigen.MatrixXi()
 
 igl.readOFF(TUTORIAL_SHARED_PATH + "decimated-knight.off", V, F)
 
-# 100 random indicies into rows of F
+# 100 random indices into rows of F
 I = igl.eigen.MatrixXi()
 igl.floor((0.5 * (igl.eigen.MatrixXd.Random(100, 1) + 1.) * F.rows()), I)
 
-# 50 random indicies into rows of I
+# 50 random indices into rows of I
 J = igl.eigen.MatrixXi()
 igl.floor((0.5 * (igl.eigen.MatrixXd.Random(50, 1) + 1.) * I.rows()), J)
 

python/tutorial/403_BoundedBiharmonicWeights.py

@@ -27,7 +27,7 @@ def pre_draw(viewer):
         for e in range(len(pose)):
             anim_pose[e] = pose[e].slerp(anim_t, igl.eigen.Quaterniond.Identity())
 
-        # Propogate relative rotations via FK to retrieve absolute transformations
+        # Propagate relative rotations via FK to retrieve absolute transformations
         vQ = igl.RotationList()
         vT = []
         igl.forward_kinematics(C, BE, P, anim_pose, vQ, vT)

python/tutorial/404_DualQuaternionSkinning.py

@@ -34,7 +34,7 @@ def pre_draw(viewer):
         for e in range(len(poses[begin])):
             anim_pose.append(poses[begin][e].slerp(t, poses[end][e]))
 
-        # Propogate relative rotations via FK to retrieve absolute transformations
+        # Propagate relative rotations via FK to retrieve absolute transformations
         vQ = igl.RotationList()
         vT = []
         igl.forward_kinematics(C, BE, P, anim_pose, vQ, vT)