added cot3 references

Former-commit-id: 6aabba3e84b2c1c27060fd6ad52c497f4a1fd834

documentation/implemented-papers.tex

@@ -3,7 +3,7 @@
 
 \title{Papers implemented in \textsc{libigl}}
 \author{Alec Jacobson}
-\date{4 April 2014}
+\date{last revised 22 April 2014}
 
 \begin{document}
 This document serves as a companion reference to better list the references to
@@ -12,7 +12,9 @@ incomplete.
 
 \paragraph{\texttt{cotmatrix}, \texttt{massmatrix}}
 build discrete operators on triangle and tetrahedral meshes. 
-\cite{Pinkall:1993:CDM,meyer03ddo,Jacobson:THESIS:2013}. 
+\cite{Pinkall:1993:CDM,meyer03ddo,Jacobson:THESIS:2013}. For tet meshes, we no
+longer use the ``by the book'' FEM construction \`a la \cite{Sharf:2007vv},
+rather a purely geometric approach \cite{Barth:1994,Xu:1999}.
 
 \paragraph{\texttt{harmonic}} solves a Laplace problem (equivalently
 minimizes the Dirichlet energy) with some simple boundary conditions
@@ -38,17 +40,6 @@ choosing skinning transformations \cite{Jacobson:FAST:2012}.
 ``skeletal subspace deformation'', or ``enveloping''. This technique is often
 attributed to \cite{Magnenat-Thalmann:1988:JLD}.
 
-  % appears in the appendix of: ``Interactive Topology-aware Surface
-  % Reconstruction,'' by Sharf, A. et al
-  %
-  % versus
-  %
-  % ND derivation given in "A MONOTONE FINITE ELEMENT SCHEME FOR
-  % CONVECTION-DIFFUSION EQUATIONS" [Xu & ZIKATANOV 1999]
-  %
-  % 3D derivation given in "Aspects of unstructured grids and finite-volume
-  % solvers for the Euler and Navier-Stokes equations" [Barth 1992]
-
 \bibliographystyle{acmsiggraph}
 \bibliography{references} 
 

documentation/references.bib.REMOVED.git-id

@@ -1 +1 @@
-d5780af0c6c5e10e0a3b5ffacff51e396a71ff2c
+92b6b20756fd71dd62d7b52066a7bcf07e4acd28

include/igl/cotangent.cpp

@@ -84,79 +84,6 @@ IGL_INLINE void igl::cotangent(
       // http://arxiv.org/pdf/1208.0354.pdf Page 18
       C = (1./6.) * l.array() * cos_theta.array() / sin_theta.array();
 
-// LEGACY
-//      // Tetrahedra
-//      typedef Matrix<typename MatV::Scalar,3,1> Vec3;
-//      typedef Matrix<typename MatV::Scalar,3,3> Mat3;
-//      typedef Matrix<typename MatV::Scalar,3,4> Mat3x4;
-//      typedef Matrix<typename MatV::Scalar,4,4> Mat4x4;
-//
-//      // preassemble right hand side
-//      // COLUMN-MAJOR ORDER FOR LAPACK
-//      Mat3x4 rhs;
-//      rhs <<
-//        1,0,0,-1,
-//        0,1,0,-1,
-//        0,0,1,-1;
-//
-//      bool diag_all_pos = true;
-//      C.resize(m,6);
-//
-//      // loop over tetrahedra
-//      for(int j = 0;j<F.rows();j++)
-//      {
-//        // points a,b,c,d make up the tetrahedra
-//        size_t a = F(j,0);
-//        size_t b = F(j,1);
-//        size_t c = F(j,2);
-//        size_t d = F(j,3);
-//        //const std::vector<double> & pa = vertices[a];
-//        //const std::vector<double> & pb = vertices[b];
-//        //const std::vector<double> & pc = vertices[c];
-//        //const std::vector<double> & pd = vertices[d];
-//        Vec3 pa = V.row(a);
-//        Vec3 pb = V.row(b);
-//        Vec3 pc = V.row(c);
-//        Vec3 pd = V.row(d);
-//
-//        // Following definition that appears in the appendix of: ``Interactive
-//        // Topology-aware Surface Reconstruction,'' by Sharf, A. et al
-//        // http://www.cs.bgu.ac.il/~asharf/Projects/InSuRe/Insure_siggraph_final.pdf
-//
-//        // compute transpose of jacobian Jj
-//        Mat3 JTj;
-//        JTj.row(0) = pa-pd;
-//        JTj.row(1) = pb-pd;
-//        JTj.row(2) = pc-pd;
-//
-//        // compute abs(determinant of JTj)/6 (volume of tet)
-//        // determinant of transpose of A equals determinant of A
-//        double volume = fabs(JTj.determinant())/6.0;
-//        //printf("volume[%d] = %g\n",j+1,volume);
-//
-//        // solve Jj' * Ej = [-I -1], for Ej
-//        // in other words solve JTj * Ej = [-I -1], for Ej
-//        Mat3x4 Ej = JTj.inverse() * rhs;
-//        // compute Ej'*Ej
-//        Mat4x4 EjTEj = Ej.transpose() * Ej;
-//
-//        // Kj =  det(JTj)/6 * Ej'Ej 
-//        Mat4x4 Kj = EjTEj*volume;
-//        diag_all_pos &= ((Kj(0,0)>0) & (Kj(1,1)>0)) & ((Kj(2,2)>0) & (Kj(3,3)>0));
-//        C(j,0) = Kj(1,2);
-//        C(j,1) = Kj(2,0);
-//        C(j,2) = Kj(0,1);
-//        C(j,3) = Kj(3,0);
-//        C(j,4) = Kj(3,1);
-//        C(j,5) = Kj(3,2);
-//      }
-//      if(diag_all_pos)
-//      {
-//#ifdef VERBOSE 
-//        verbose("cotangent.h: Flipping sign of cotangent, so that cots are positive\n");
-//#endif
-//        C *= -1.0;
-//      }
       break;
     }
     default:

include/igl/gradMat.cpp

@@ -5,101 +5,102 @@
 // 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 "gradMat.h"
-#include <vector>
-
-template <typename T, typename S>
-IGL_INLINE void igl::gradMat(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
-  const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> &F,
-  Eigen::SparseMatrix<T> &G )
-{
-  Eigen::PlainObjectBase<Eigen::Matrix<T,Eigen::Dynamic,3> > eperp21, eperp13;
-  eperp21.resize(F.rows(),3);
-  eperp13.resize(F.rows(),3);
-
-  for (int i=0;i<F.rows();++i)
-  {
-    // renaming indices of vertices of triangles for convenience
-    int i1 = F(i,0);
-    int i2 = F(i,1);
-    int i3 = F(i,2);
-    
-    // #F x 3 matrices of triangle edge vectors, named after opposite vertices
-    Eigen::Matrix<T, 1, 3> v32 = V.row(i3) - V.row(i2);
-    Eigen::Matrix<T, 1, 3> v13 = V.row(i1) - V.row(i3);
-    Eigen::Matrix<T, 1, 3> v21 = V.row(i2) - V.row(i1);
-    
-    // area of parallelogram is twice area of triangle
-    // area of parallelogram is || v1 x v2 || 
-    Eigen::Matrix<T, 1, 3> n  = v32.cross(v13); 
-    
-    // This does correct l2 norm of rows, so that it contains #F list of twice
-    // triangle areas
-    double dblA = std::sqrt(n.dot(n));
-    
-    // now normalize normals to get unit normals
-    Eigen::Matrix<T, 1, 3> u = n / dblA;
-    
-    // rotate each vector 90 degrees around normal
-    double norm21 = std::sqrt(v21.dot(v21));
-    double norm13 = std::sqrt(v13.dot(v13));
-    eperp21.row(i) = u.cross(v21);
-    eperp21.row(i) = eperp21.row(i) / std::sqrt(eperp21.row(i).dot(eperp21.row(i)));
-    eperp21.row(i) *= norm21 / dblA;
-    eperp13.row(i) = u.cross(v13);
-    eperp13.row(i) = eperp13.row(i) / std::sqrt(eperp13.row(i).dot(eperp13.row(i)));
-    eperp13.row(i) *= norm13 / dblA;
-  }
-
-  std::vector<int> rs;
-  rs.reserve(F.rows()*4*3);
-  std::vector<int> cs;
-  cs.reserve(F.rows()*4*3);
-  std::vector<double> vs;
-  vs.reserve(F.rows()*4*3);
-
-  // row indices
-  for(int r=0;r<3;r++)
-  {
-    for(int j=0;j<4;j++)
-    {
-      for(int i=r*F.rows();i<(r+1)*F.rows();i++) rs.push_back(i);
-    }
-  }
-
-  // column indices
-  for(int r=0;r<3;r++)
-  {
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,1));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,2));
-    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
-  }
-  
-  // values
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,0));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,1));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,2));
-  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,2));
-
-  // create sparse gradient operator matrix
-  G.resize(3*F.rows(),V.rows()); 
-  std::vector<Eigen::Triplet<T> > triplets;
-  for (int i=0;i<vs.size();++i)
-  {
-    triplets.push_back(Eigen::Triplet<T>(rs[i],cs[i],vs[i]));
-  }
-  G.setFromTriplets(triplets.begin(), triplets.end());
-}
-
-#ifndef IGL_HEADER_ONLY
-// Explicit template specialization
-#endif
+#include "gradMat.h"
+#include <vector>
+
+template <typename T, typename S>
+IGL_INLINE void igl::gradMat(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &V,
+  const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic> &F,
+  Eigen::SparseMatrix<T> &G )
+{
+  Eigen::PlainObjectBase<Eigen::Matrix<T,Eigen::Dynamic,3> > eperp21, eperp13;
+  eperp21.resize(F.rows(),3);
+  eperp13.resize(F.rows(),3);
+
+  for (int i=0;i<F.rows();++i)
+  {
+    // renaming indices of vertices of triangles for convenience
+    int i1 = F(i,0);
+    int i2 = F(i,1);
+    int i3 = F(i,2);
+    
+    // #F x 3 matrices of triangle edge vectors, named after opposite vertices
+    Eigen::Matrix<T, 1, 3> v32 = V.row(i3) - V.row(i2);
+    Eigen::Matrix<T, 1, 3> v13 = V.row(i1) - V.row(i3);
+    Eigen::Matrix<T, 1, 3> v21 = V.row(i2) - V.row(i1);
+    
+    // area of parallelogram is twice area of triangle
+    // area of parallelogram is || v1 x v2 || 
+    Eigen::Matrix<T, 1, 3> n  = v32.cross(v13); 
+    
+    // This does correct l2 norm of rows, so that it contains #F list of twice
+    // triangle areas
+    double dblA = std::sqrt(n.dot(n));
+    
+    // now normalize normals to get unit normals
+    Eigen::Matrix<T, 1, 3> u = n / dblA;
+    
+    // rotate each vector 90 degrees around normal
+    double norm21 = std::sqrt(v21.dot(v21));
+    double norm13 = std::sqrt(v13.dot(v13));
+    eperp21.row(i) = u.cross(v21);
+    eperp21.row(i) = eperp21.row(i) / std::sqrt(eperp21.row(i).dot(eperp21.row(i)));
+    eperp21.row(i) *= norm21 / dblA;
+    eperp13.row(i) = u.cross(v13);
+    eperp13.row(i) = eperp13.row(i) / std::sqrt(eperp13.row(i).dot(eperp13.row(i)));
+    eperp13.row(i) *= norm13 / dblA;
+  }
+
+  std::vector<int> rs;
+  rs.reserve(F.rows()*4*3);
+  std::vector<int> cs;
+  cs.reserve(F.rows()*4*3);
+  std::vector<double> vs;
+  vs.reserve(F.rows()*4*3);
+
+  // row indices
+  for(int r=0;r<3;r++)
+  {
+    for(int j=0;j<4;j++)
+    {
+      for(int i=r*F.rows();i<(r+1)*F.rows();i++) rs.push_back(i);
+    }
+  }
+
+  // column indices
+  for(int r=0;r<3;r++)
+  {
+    for(int i=0;i<F.rows();i++) cs.push_back(F(i,1));
+    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
+    for(int i=0;i<F.rows();i++) cs.push_back(F(i,2));
+    for(int i=0;i<F.rows();i++) cs.push_back(F(i,0));
+  }
+  
+  // values
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,0));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,0));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,0));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,0));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,1));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,1));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,1));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,1));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp13(i,2));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp13(i,2));
+  for(int i=0;i<F.rows();i++) vs.push_back(eperp21(i,2));
+  for(int i=0;i<F.rows();i++) vs.push_back(-eperp21(i,2));
+
+  // create sparse gradient operator matrix
+  G.resize(3*F.rows(),V.rows()); 
+  std::vector<Eigen::Triplet<T> > triplets;
+  for (int i=0;i<vs.size();++i)
+  {
+    triplets.push_back(Eigen::Triplet<T>(rs[i],cs[i],vs[i]));
+  }
+  G.setFromTriplets(triplets.begin(), triplets.end());
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+#endif
+

todos.txt

@@ -1,10 +1,9 @@
+- sort out `grad.*` vs `gradMat.*`
 - clean up externals
   - What's there for convenience?
   - What's there because it's patched/
 - compile on Linux, Mac OS X, Windows
 - unit tests for all functions
-- standardize name of library "libigl", purge IGL LIB, igl_lib, igl lib, IGL
-    Library, etc.
 - clean up examples
 - standardize use of string/char *, add to style conventions
 - standardize Eigen Templates, add to style conventions
@@ -19,13 +18,13 @@
 - fix bugs in examples/Core/example2.cpp
 - replace generic names read.h/write.h with something like read_poly_mesh.h
 - replace DynamicSparseMatrix: coeffRef += ---> setFromTriplets()
-- create libigltga extra
+- create libigltga extra...why? license?
 - rename moveFV.h
 - is_border_vertex.h should not require V
 - consistent checks/asserts for: manifoldness, closedness, triangles-only
-    e.g. check_mesh(V,F,MANIFOLD | CLOSED | TRIANGLES_ONLY)
+    e.g. check_mesh(V,F,IS_MANIFOLD | IS_CLOSED | IS_TRIANGLES_ONLY)
 - clean up Timer.h
-- svd.* depends on lapack, should use eigen...
+x svd.* depends on lapack, should use eigen...
 - use preprocessor macros to automatically include .cpp files at end of .h
 - unify include opengl with convenience includes
 - MatrixBase --> PlainObjectBase
@@ -33,3 +32,5 @@
 - everywhere appropriate change:
   `#pragma omp parallel for` to `#pragma omp parallel for if (n>10000)` where
   `n` and `1000` are replaced accordingly
++ standardize name of library "libigl", purge IGL LIB, igl_lib, igl lib, IGL
+    Library, etc.