ONE variables, cotangent for each edge of each face, snotes in zum

Former-commit-id: e4f1aee51772230fa3d31712d88653b0d52a54a5

ONE.h

@@ -0,0 +1,15 @@
+#ifndef ONE_H
+#define ONE_H
+// Often one needs a reference to a dummy variable containing zero as its
+// value, for example when using AntTweakBar's
+// TwSetParam( "3D View", "opened", TW_PARAM_INT32, 1, &INT_ONE);
+namespace igl
+{
+  const char CHAR_ONE = 1;
+  const int INT_ONE = 1;
+  const unsigned int UNSIGNED_INT_ONE = 1;
+  const double DOUBLE_ONE = 1;
+  const float FLOAT_ONE = 1;
+}
+#endif
+

ZERO.h

@@ -1,6 +1,8 @@
+#ifndef ZERO_H
+#define ZERO_H
 // Often one needs a reference to a dummy variable containing zero as its
 // value, for example when using AntTweakBar's
-// TwSetParam( "3D View", "opened", TW_PARAM_INT32, 1, &TW_ZERO);
+// TwSetParam( "3D View", "opened", TW_PARAM_INT32, 1, &INT_ZERO);
 namespace igl
 {
   const char CHAR_ZERO = 0;
@@ -9,3 +11,4 @@ namespace igl
   const double DOUBLE_ZERO = 0;
   const float FLOAT_ZERO = 0;
 }
+#endif

cotangent.h

@@ -0,0 +1,146 @@
+#ifndef IGL_COTANGENT_H
+#define IGL_COTANGENT_H
+namespace igl
+{
+  // COTANGENT compute the cotangents of each angle in mesh (V,F)
+  // 
+  // Templates:
+  //   MatV  vertex position matrix, e.g. Eigen::MatrixXd
+  //   MatF  face index matrix, e.g. Eigen::MatrixXd
+  //   MatC  cotangent weights matrix, e.g. Eigen::MatrixXd
+  // Inputs:
+  //   V  #V by dim list of rest domain positions
+  //   F  #F by {3|4} list of {triangle|tetrahedra} indices into V
+  // Outputs:
+  //   C  #F by {3|6} list of cotangents corresponding angles
+  //     for triangles, columns correspond to edges 23,31,12
+  //     for tets, columns correspond to edges 23,31,12,41,42,43
+  template <class MatV, class MatF, class MatC>
+  void cotangent(const MatV & V, const MatF & F, MatC & C);
+}
+
+// Implementation
+#include <verbose.h>
+#include <Eigen/Dense>
+
+template <class MatV, class MatF, class MatC>
+void igl::cotangent(const MatV & V, const MatF & F, MatC & C)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  // simplex size (3: triangles, 4: tetrahedra)
+  int simplex_size = F.cols();
+  // Number of elements
+  int m = F.rows();
+
+  if(simplex_size == 3)
+  {
+    // Triangles
+    // edge lengths numbered same as opposite vertices
+    Matrix<typename MatC::Scalar,Dynamic,3> l(m,3);
+    // loop over faces
+    for(int i = 0;i<m;i++)
+    {
+      l(i,0) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
+      l(i,1) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
+      l(i,2) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
+    }
+    // semiperimeters
+    Matrix<typename MatC::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
+    assert(s.rows() == m);
+    // Heron's forumal for area
+    Matrix<typename MatC::Scalar,Dynamic,1> dblA(m);
+    for(int i = 0;i<m;i++)
+    {
+      dblA(i) = 2.0*sqrt(s(i)*(s(i)-l(i,0))*(s(i)-l(i,1))*(s(i)-l(i,2)));
+    }
+    // cotangents and diagonal entries for element matrices
+    // correctly divided by 4 (alec 2010)
+    C.resize(m,3);
+    for(int i = 0;i<m;i++)
+    {
+      C(i,0) = (l(i,1)*l(i,1) + l(i,2)*l(i,2) - l(i,0)*l(i,0))/dblA(i)/4.0;
+      C(i,1) = (l(i,2)*l(i,2) + l(i,0)*l(i,0) - l(i,1)*l(i,1))/dblA(i)/4.0;
+      C(i,2) = (l(i,0)*l(i,0) + l(i,1)*l(i,1) - l(i,2)*l(i,2))/dblA(i)/4.0;
+    }
+  }else if(simplex_size == 4)
+  {
+    // Tetrahedra
+    typedef Matrix<typename MatC::Scalar,3,1> Vec3;
+    typedef Matrix<typename MatC::Scalar,3,3> Mat3;
+    typedef Matrix<typename MatC::Scalar,3,4> Mat3x4;
+    typedef Matrix<typename MatC::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)
+    {
+      verbose("cotangent.h: Flipping sign of cotangent, so that cots are positive\n");
+      C *= -1.0;
+    }
+  }else
+  {
+    fprintf(stderr,
+      "cotangent.h: Error: Simplex size (%d) not supported\n", simplex_size);
+    assert(false);
+  }
+}
+#endif

examples/ReAntTweakBar/static/libAntTweakBar.a.REMOVED.git-id

@@ -1 +1 @@
-c99c7d486e0f5e86a03def2cfb181bd2c5ff421a
+fe412ef5ceaf9aa2985e8a41fe97550bf8460f62

min_quad_dense.h

@@ -32,7 +32,10 @@ namespace igl
 		Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& S)
 	{
 		typedef Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Mat;
-		const T treshold = igl::DOUBLE_EPS;
+                // This threshold seems to matter a lot but I'm not sure how to
+                // set it
+		const T treshold = igl::FLOAT_EPS;
+		//const T treshold = igl::DOUBLE_EPS;
 
 		const int n = A.rows();
 		assert(A.cols() == n);

readMESH.h

@@ -0,0 +1,294 @@
+#ifndef IGL_READMESH_H
+#define IGL_READMESH_H
+
+#include <string>
+#include <vector>
+
+namespace igl
+{
+  // load a tetrahedral volume mesh from a .mesh file
+  //
+  // Templates:
+  //   Scalar  type for positions and vectors (will be read as double and cast
+  //     to Scalar)
+  //   Index  type for indices (will be read as int and cast to Index)
+  // Input:
+  //   mesh_file_name  path of .mesh file
+  // Outputs:
+  //   V  double matrix of vertex positions  #V by 3
+  //   T  #T list of tet indices into vertex positions
+  //   F  #F list of face indices into vertex positions
+  template <typename Scalar, typename Index>
+  inline bool readMESH(
+    const std::string mesh_file_name,
+    std::vector<std::vector<Scalar > > & V,
+    std::vector<std::vector<Index > > & T,
+    std::vector<std::vector<Index > > & F);
+
+  // Input:
+  //   mesh_file_name  path of .mesh file
+  // Outputs:
+  //   V  eigen double matrix #V by 3
+  //   T  eigen int matrix #T by 4
+  //   F  eigen int matrix #F by 3
+  inline bool readMESH(
+    const std::string str,
+    Eigen::MatrixXd& V,
+    Eigen::MatrixXi& T,
+    Eigen::MatrixXi& F);
+}
+
+// Implementation
+#include <cstdio>
+#include <verbose.h>;
+
+template <typename Scalar, typename Index>
+inline bool igl::readMESH(
+  const std::string mesh_file_name,
+  std::vector<std::vector<Scalar > > & V,
+  std::vector<std::vector<Index > > & T,
+  std::vector<std::vector<Index > > & F)
+{
+  using namespace std;
+  using namespace igl;
+  FILE * mesh_file = fopen(mesh_file_name.c_str(),"r");
+  if(NULL==mesh_file)
+  {
+    fprintf(stderr,"IOError: %s could not be opened...",mesh_file_name.c_str());
+    return false;
+  }
+#ifndef LINE_MAX
+#  define LINE_MAX 2048
+#endif
+  char line[LINE_MAX];
+  bool still_comments;
+  V.clear();
+  T.clear();
+  F.clear();
+
+  // eat comments at beginning of file
+  still_comments= true;
+  while(still_comments)
+  {
+    fgets(line,LINE_MAX,mesh_file);
+    still_comments = (line[0] == '#' || line[0] == '\n');
+  }
+
+  char str[LINE_MAX];
+  sscanf(line," %s",str);
+  // check that first word is MeshVersionFormatted
+  if(0!=strcmp(str,"MeshVersionFormatted"))
+  {
+    fprintf(stderr,
+      "Error: first word should be MeshVersionFormatted not %s\n",str);
+    fclose(mesh_file);
+    return false;
+  }
+  int one = -1;
+  if(2 != sscanf(line,"%s %d",str,&one))
+  {
+    // 1 appears on next line?
+    fscanf(mesh_file," %d",&one);
+  }
+  if(one != 1)
+  {
+    fprintf(stderr,"Error: second word should be 1 not %d\n",one);
+    fclose(mesh_file);
+    return false;
+  }
+
+  // eat comments
+  still_comments= true;
+  while(still_comments)
+  {
+    fgets(line,LINE_MAX,mesh_file);
+    still_comments = (line[0] == '#' || line[0] == '\n');
+  }
+
+  sscanf(line," %s",str);
+  // check that third word is Dimension
+  if(0!=strcmp(str,"Dimension"))
+  {
+    fprintf(stderr,"Error: third word should be Dimension not %s\n",str);
+    fclose(mesh_file);
+    return false;
+  }
+  int three = -1;
+  if(2 != sscanf(line,"%s %d",str,&three))
+  {
+    // 1 appears on next line?
+    fscanf(mesh_file," %d",&three);
+  }
+  if(three != 3)
+  {
+    fprintf(stderr,"Error: only Dimension 3 supported not %d\n",three);
+    fclose(mesh_file);
+    return false;
+  }
+
+  // eat comments
+  still_comments= true;
+  while(still_comments)
+  {
+    fgets(line,LINE_MAX,mesh_file);
+    still_comments = (line[0] == '#' || line[0] == '\n');
+  }
+
+  sscanf(line," %s",str);
+  // check that fifth word is Vertices
+  if(0!=strcmp(str,"Vertices"))
+  {
+    fprintf(stderr,"Error: fifth word should be Vertices not %s\n",str);
+    fclose(mesh_file);
+    return false;
+  }
+  size_t number_of_vertices;
+  if(1 != fscanf(mesh_file," %ld",&number_of_vertices))
+  {
+    fprintf(stderr,"Error: expecting number of vertices...\n");
+    fclose(mesh_file);
+    return false;
+  }
+  // allocate space for vertices
+  V.resize(number_of_vertices,vector<Scalar>(3,0));
+  size_t extra;
+  for(size_t i = 0;i<number_of_vertices;i++)
+  {
+    double x,y,z;
+    if(4 != fscanf(mesh_file," %lg %lg %lg %ld",&x,&y,&z,&extra))
+    {
+      fprintf(stderr,"Error: expecting vertex position...\n");
+      fclose(mesh_file);
+      return false;
+    }
+    V[i][0] = x;
+    V[i][1] = y;
+    V[i][2] = z;
+  }
+
+  // eat comments
+  still_comments= true;
+  while(still_comments)
+  {
+    fgets(line,LINE_MAX,mesh_file);
+    still_comments = (line[0] == '#' || line[0] == '\n');
+  }
+
+  sscanf(line," %s",str);
+  // check that sixth word is Triangles
+  if(0!=strcmp(str,"Triangles"))
+  {
+    fprintf(stderr,"Error: sixth word should be Triangles not %s\n",str);
+    fclose(mesh_file);
+    return false;
+  }
+  size_t number_of_triangles;
+  if(1 != fscanf(mesh_file," %ld",&number_of_triangles))
+  {
+    fprintf(stderr,"Error: expecting number of triangles...\n");
+    fclose(mesh_file);
+    return false;
+  }
+  // allocate space for triangles
+  F.resize(number_of_triangles,vector<Index>(3));
+  // triangle indices
+  size_t tri[3];
+  for(size_t i = 0;i<number_of_triangles;i++)
+  {
+    if(4 != fscanf(mesh_file," %ld %ld %ld %ld",&tri[0],&tri[1],&tri[2],&extra))
+    {
+      printf("Error: expecting triangle indices...\n");
+      return false;
+    }
+    for(size_t j = 0;j<3;j++)
+    {
+      F[i][j] = tri[j]-1;
+    }
+  }
+
+  // eat comments
+  still_comments= true;
+  while(still_comments)
+  {
+    fgets(line,LINE_MAX,mesh_file);
+    still_comments = (line[0] == '#' || line[0] == '\n');
+  }
+
+  sscanf(line," %s",str);
+  // check that sixth word is Triangles
+  if(0!=strcmp(str,"Tetrahedra"))
+  {
+    fprintf(stderr,"Error: seventh word should be Tetrahedra not %s\n",str);
+    fclose(mesh_file);
+    return false;
+  }
+  size_t number_of_tetrahedra;
+  if(1 != fscanf(mesh_file," %ld",&number_of_tetrahedra))
+  {
+    fprintf(stderr,"Error: expecting number of tetrahedra...\n");
+    fclose(mesh_file);
+    return false;
+  }
+  // allocate space for tetrahedra
+  T.resize(number_of_tetrahedra,vector<Index>(4));
+  // tet indices
+  size_t a,b,c,d;
+  for(size_t i = 0;i<number_of_tetrahedra;i++)
+  {
+    if(5 != fscanf(mesh_file," %ld %ld %ld %ld %ld",&a,&b,&c,&d,&extra))
+    {
+      fprintf(stderr,"Error: expecting tetrahedra indices...\n");
+      fclose(mesh_file);
+      return false;
+    }
+    T[i][0] = a-1;
+    T[i][1] = b-1;
+    T[i][2] = c-1;
+    T[i][3] = d-1;
+  }
+  fclose(mesh_file);
+  return true;
+}
+
+#include <Eigen/Core>
+#include "list_to_matrix.h"
+
+inline bool igl::readMESH(
+  const std::string str,
+  Eigen::MatrixXd& V,
+  Eigen::MatrixXi& T,
+  Eigen::MatrixXi& F)
+{
+  std::vector<std::vector<double> > vV,vT,vF;
+  bool success = igl::readMESH(str,vV,vT,vF);
+  if(!success)
+  {
+    // readOBJ(str,vV,vTC,vN,vF,vFTC,vFN) should have already printed an error
+    // message to stderr
+    return false;
+  }
+  bool V_rect = igl::list_to_matrix(vV,V);
+  if(!V_rect)
+  {
+    // igl::list_to_matrix(vV,V) already printed error message to std err
+    return false;
+  }
+  bool T_rect = igl::list_to_matrix(vT,T);
+  if(!T_rect)
+  {
+    // igl::list_to_matrix(vT,T) already printed error message to std err
+    return false;
+  }
+  bool F_rect = igl::list_to_matrix(vF,F);
+  if(!F_rect)
+  {
+    // igl::list_to_matrix(vF,F) already printed error message to std err
+    return false;
+  }
+  assert(V.cols() == 3);
+  assert(T.cols() == 4);
+  assert(F.cols() == 3);
+  return true;
+}
+
+#endif

readOBJ.h

@@ -17,19 +17,6 @@
 
 namespace igl 
 {
-  //! Read a mesh from an ascii obj file
-  // Inputs:
-  //   str  path to .obj file
-  // Outputs:
-  //   V  eigen double matrix #V by 3
-  //   F  eigen int matrix #F by 3
-  //
-  // KNOWN BUG: This only knows how to read *triangle* meshes. It will probably
-  // crash or give garbage on anything else.
-  //
-  // KNOWN BUG: This only knows how to face lines without normal or texture
-  // indices. It will probably crash or give garbage on anything else.
-  inline bool readOBJ(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F);
   // Read a mesh from an ascii obj file, filling in vertex positions, normals
   // and texture coordinates. Mesh may have faces of any number of degree
   //
@@ -56,6 +43,20 @@ namespace igl
     std::vector<std::vector<Index > > & F,
     std::vector<std::vector<Index > > & FTC,
     std::vector<std::vector<Index > > & FN);
+
+  //! Read a mesh from an ascii obj file
+  // Inputs:
+  //   str  path to .obj file
+  // Outputs:
+  //   V  eigen double matrix #V by 3
+  //   F  eigen int matrix #F by 3
+  //
+  // KNOWN BUG: This only knows how to read *triangle* meshes. It will probably
+  // crash or give garbage on anything else.
+  //
+  // KNOWN BUG: This only knows how to face lines without normal or texture
+  // indices. It will probably crash or give garbage on anything else.
+  inline bool readOBJ(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F);
 }
 
 // Implementation
@@ -64,42 +65,6 @@ namespace igl
 #include <iostream>
 #include <fstream>
 
-inline bool igl::readOBJ(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
-{
-  std::vector<std::vector<double> > vV,vTC,vN;
-  std::vector<std::vector<int> > vF,vFTC,vFN;
-  bool success = igl::readOBJ(str,vV,vTC,vN,vF,vFTC,vFN);
-  if(!success)
-  {
-    // readOBJ(str,vV,vTC,vN,vF,vFTC,vFN) should have already printed an error
-    // message to stderr
-    return false;
-  }
-  bool V_rect = igl::list_to_matrix(vV,V);
-  if(!V_rect)
-  {
-    // igl::list_to_matrix(vV,V) already printed error message to std err
-    return false;
-  }
-  bool F_rect = igl::list_to_matrix(vF,F);
-  if(!F_rect)
-  {
-    // igl::list_to_matrix(vF,F) already printed error message to std err
-    return false;
-  }
-  // Legacy
-  if(F.cols() != 3)
-  {
-    fprintf(stderr,
-      "Error: readOBJ(filename,V,F) is meant for reading triangle-only"
-      " meshes. This mesh has faces all with size %d. See readOBJ.h for other"
-      " options.\n",
-      (int)F.cols());
-    return false;
-  }
-  return true;
-}
-
 template <typename Scalar, typename Index>
 inline bool igl::readOBJ(
   const std::string obj_file_name, 
@@ -292,4 +257,39 @@ inline bool igl::readOBJ(
   return true;
 }
 
+inline bool igl::readOBJ(const std::string str, Eigen::MatrixXd& V, Eigen::MatrixXi& F)
+{
+  std::vector<std::vector<double> > vV,vTC,vN;
+  std::vector<std::vector<int> > vF,vFTC,vFN;
+  bool success = igl::readOBJ(str,vV,vTC,vN,vF,vFTC,vFN);
+  if(!success)
+  {
+    // readOBJ(str,vV,vTC,vN,vF,vFTC,vFN) should have already printed an error
+    // message to stderr
+    return false;
+  }
+  bool V_rect = igl::list_to_matrix(vV,V);
+  if(!V_rect)
+  {
+    // igl::list_to_matrix(vV,V) already printed error message to std err
+    return false;
+  }
+  bool F_rect = igl::list_to_matrix(vF,F);
+  if(!F_rect)
+  {
+    // igl::list_to_matrix(vF,F) already printed error message to std err
+    return false;
+  }
+  // Legacy
+  if(F.cols() != 3)
+  {
+    fprintf(stderr,
+      "Error: readOBJ(filename,V,F) is meant for reading triangle-only"
+      " meshes. This mesh has faces all with size %d. See readOBJ.h for other"
+      " options.\n",
+      (int)F.cols());
+    return false;
+  }
+  return true;
+}
 #endif

sum.h

@@ -4,8 +4,10 @@
 
 namespace igl
 {
-  // Ideally, this becomes a super overloaded function that works with sparse
-  // and dense matrices like the matlab sum function
+  // Note: If your looking for dense matrix matlab like sum for eigen matrics
+  // just use:
+  //   M.colwise().sum() or M.rowwise().sum()
+  // 
 
   // Sum the columns or rows of a sparse matrix
   // Templates: