Merge branch 'master' of ssh://github.com/libigl/libigl

Former-commit-id: b368b2dfeec66a8eab991101a96edfeac1cca4c4

Makefile

@@ -17,19 +17,26 @@ CFLAGS += $(OPTFLAGS)
 # We use well-supported features of c++11
 
 EXTRA_DIRS=
-ifeq ($(IGL_WITH_TETGEN),1)
-	# append tetgen extra dir
-	EXTRA_DIRS+=include/igl/tetgen
-	EXTRAS += tetgen
+ifeq ($(IGL_WITH_BBW),1)
+	EXTRA_DIRS+=include/igl/bbw
+	EXTRAS += bbw
+endif
+ifeq ($(IGL_WITH_BOOST),1)
+	EXTRA_DIRS+=include/igl/boost
+	EXTRAS += boost
+endif
+ifeq ($(IGL_WITH_CGAL),1)
+	EXTRA_DIRS+=include/igl/cgal
+	EXTRAS += cgal
+endif
+ifeq ($(IGL_WITH_EMBREE),1)
+	EXTRA_DIRS+=include/igl/embree
+	EXTRAS += embree
 endif
 ifeq ($(IGL_WITH_MATLAB),1)
 	EXTRA_DIRS+=include/igl/matlab
 	EXTRAS += matlab
 endif
-ifeq ($(IGL_WITH_BBW),1)
-	EXTRA_DIRS+=include/igl/bbw
-	EXTRAS += bbw
-endif
 ifeq ($(IGL_WITH_MOSEK),1)
 	EXTRA_DIRS+=include/igl/mosek
 	EXTRAS += mosek
@@ -38,22 +45,19 @@ ifeq ($(IGL_WITH_PNG),1)
 	EXTRA_DIRS+=include/igl/png
 	EXTRAS += png
 endif
-ifeq ($(IGL_WITH_XML),1)
-	EXTRA_DIRS+=include/igl/xml
-	EXTRAS += xml
-endif
-ifeq ($(IGL_WITH_EMBREE),1)
-	EXTRA_DIRS+=include/igl/embree
-	EXTRAS += embree
-endif
-ifeq ($(IGL_WITH_BOOST),1)
-	EXTRA_DIRS+=include/igl/boost
-	EXTRAS += boost
-endif
 ifeq ($(IGL_WITH_SVD3X3),1)
 	EXTRA_DIRS+=include/igl/svd3x3
 	EXTRAS += svd3x3
 endif
+ifeq ($(IGL_WITH_TETGEN),1)
+	# append tetgen extra dir
+	EXTRA_DIRS+=include/igl/tetgen
+	EXTRAS += tetgen
+endif
+ifeq ($(IGL_WITH_XML),1)
+	EXTRA_DIRS+=include/igl/xml
+	EXTRAS += xml
+endif
 
 .PHONY: examples
 .PHONY: extras

Makefile.conf

@@ -51,6 +51,7 @@ ifeq ($(IGL_USERNAME),ajx)
 	IGL_WITH_EMBREE=1
 	IGL_WITH_MATLAB=1
 	IGL_WITH_MOSEK=1
+	IGL_WITH_CGAL=1
 	IGL_WITH_BBW=1
 	IGL_WITH_SVD3X3=1
 	IGL_WITH_PNG=1

README.md

@@ -24,7 +24,11 @@ listed below.
 - libpng  libiglpng extra only
 - Mosek  libiglmosek extra only
 - Matlab  libiglmatlab extra only
-- boost  libboost extra only
+- boost  libiglboost, libiglcgal extra only
+- CGAL  libiglcgal extra only
+    * boost 
+    * gmp
+    * mpfr
 
 ### Optional (included in external/) ###
 - TetGen  libigltetgen extra only

RELEASE_HISTORY.txt

@@ -1,3 +1,4 @@
+0.4.5  CGAL extra: mesh selfintersection
 0.4.4  STL file format support
 0.4.3  ARAP implementation
 0.4.1  Migrated much of the FAST code including extra for Sifakis' 3x3 svd

VERSION.txt

@@ -3,4 +3,4 @@
 # Anyone may increment Minor to indicate a small change.
 # Major indicates a large change or large number of changes (upload to website)
 # World indicates a substantial change or release
-0.4.4
+0.4.5

documentation/active-set.tex

@@ -25,7 +25,7 @@
 \url{http://www.math.uh.edu/~rohop/fall_06/Chapter2.pdf}, and
 \url{http://www.cs.cornell.edu/courses/cs322/2007sp/notes/qr.pdf}. Mostly this
 is to put everything in the same place and use notation compatible with the
-libigl implementation. The ideas and descriptions, however, are not novel.}
+\textsc{libigl} implementation. The ideas and descriptions, however, are not novel.}
 \end{pullout}
 
 Quadratic programming problems (QPs) can be written in general as:

documentation/implemented-papers.tex

@@ -0,0 +1,47 @@
+\documentclass[12pt]{diary}
+\immediate\write18{bibtex \jobname}
+
+\title{Papers implemented in \textsc{libigl}}
+\author{Alec Jacobson}
+\date{last revised 22 April 2014}
+
+\begin{document}
+This document serves as a companion reference to better list the references to
+scientific articles implemented within \textsc{libigl}. It will no doubt be
+incomplete.
+
+\paragraph{\texttt{cotmatrix}, \texttt{massmatrix}}
+build discrete operators on triangle and tetrahedral meshes. 
+\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
+\cite{HarmonicCoodinates07}. There's also an option to solve
+``higher order Laplace problems'' (bi-Laplace, tri-Laplace, etc.)
+\cite{Botsch:2004:AIF,sorkine04lsm,Jacobson:MixedFEM:2010}.
+
+\paragraph{\texttt{bbw/}} implements ``bounded biharmonic
+weights'' \cite{Jacobson:BBW:2011}.
+
+\paragraph{\texttt{svd3x3/arap}} is a generalized implementation
+for solving ``as-rigid-as-possible'' (ARAP) mesh deformation or parameterization
+problems \cite{ARAP_modeling:2007,Liu:2008:ALA,Chao:2010:ASG}.
+
+\paragraph{\texttt{svd3x3/arap\_dof}} implements ``FAST'',
+which is simultaneously a reduced form of ARAP and a method for automatically
+choosing skinning transformations \cite{Jacobson:FAST:2012}.
+
+\paragraph{\texttt{dqs}} implements ``Dual quaternion skinning''
+\cite{Kavan:2008:GSW}.
+
+\paragraph{\texttt{lbs}} implements ``linear blend skinning'', also known as
+``skeletal subspace deformation'', or ``enveloping''. This technique is often
+attributed to \cite{Magnenat-Thalmann:1988:JLD}.
+
+\bibliographystyle{acmsiggraph}
+\bibliography{references} 
+
+\end{document}
+__END__

documentation/references.bib.REMOVED.git-id

@@ -0,0 +1 @@
+92b6b20756fd71dd62d7b52066a7bcf07e4acd28

examples/bbw/main.cpp

@@ -257,7 +257,7 @@ int main(int argc, char * argv[])
   // Default bbw data and flags
   BBWData bbw_data;
   bbw_data.qp_solver = QP_SOLVER_IGL_ACTIVE_SET;
-  bbw_data.qp_solver = QP_SOLVER_MOSEK;
+  //bbw_data.qp_solver = QP_SOLVER_MOSEK;
   // Weights matrix
   MatrixXd W;
   if(!bbw(VV,TT,b,bc,bbw_data,W))

examples/skeleton/example.cpp

@@ -252,7 +252,10 @@ void display()
   glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, SILVER_SPECULAR);
   glMaterialf (GL_FRONT_AND_BACK, GL_SHININESS, 128);
 
-  vector<Quaterniond> dQ(BE.rows(),Quaterniond::Identity()),vQ;
+  typedef std::vector<
+    Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> >
+    RotationList;
+  RotationList dQ(BE.rows(),Quaterniond::Identity()),vQ;
   vector<Vector3d> vT;
   Matrix3d A = Matrix3d::Identity();
   for(int e = 0;e<BE.rows();e++)

include/igl/active_set.h

@@ -80,7 +80,7 @@ struct igl::active_set_params
 {
   // Input parameters for active_set:
   //   Auu_pd  whether Auu is positive definite {false}
-  //   max_iter  Maximum number of iterations ({0} = Infinity)
+  //   max_iter  Maximum number of iterations (0 = Infinity, {100})
   //   inactive_threshold  Threshold on Lagrange multiplier values to determine
   //     whether to keep constraints active {EPS}
   //   constraint_threshold  Threshold on whether constraints are violated (0
@@ -94,7 +94,7 @@ struct igl::active_set_params
   double solution_diff_threshold;
   active_set_params():
     Auu_pd(false),
-    max_iter(-1),
+    max_iter(100),
     inactive_threshold(igl::DOUBLE_EPS),
     constraint_threshold(igl::DOUBLE_EPS),
     solution_diff_threshold(igl::DOUBLE_EPS)

include/igl/cgal/CGAL_includes.hpp

@@ -0,0 +1,37 @@
+#ifndef IGL_CGAL_INCLUDES_H
+#define IGL_CGAL_INCLUDES_H
+
+// Triangle triangle intersection
+#include <CGAL/intersections.h>
+// THIS CANNOT BE INCLUDED IN THE SAME FILE AS <CGAL/intersections.h>
+// #include <CGAL/Boolean_set_operations_2.h>
+
+// Constrained Delaunay Triangulation types
+#include <CGAL/Constrained_Delaunay_triangulation_2.h>
+#include <CGAL/Constrained_triangulation_plus_2.h>
+
+// Axis-align boxes for all-pairs self-intersection detection
+#include <CGAL/point_generators_3.h>
+#include <CGAL/Bbox_3.h>
+#include <CGAL/box_intersection_d.h>
+#include <CGAL/function_objects.h>
+#include <CGAL/Join_input_iterator.h>
+#include <CGAL/algorithm.h>
+#include <vector>
+
+// Axis-aligned bounding box tree for tet tri intersection
+#include <CGAL/AABB_tree.h>
+#include <CGAL/AABB_traits.h>
+#include <CGAL/AABB_triangle_primitive.h>
+
+// Boolean operations
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/Nef_polyhedron_3.h>
+
+// Delaunay Triangulation in 3D
+#include <CGAL/Delaunay_triangulation_3.h>
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Exact_predicates_exact_constructions_kernel.h>
+
+#endif

include/igl/cgal/Makefile

@@ -0,0 +1,42 @@
+
+.PHONY: all
+all: libiglcgal
+debug: libiglcgal
+
+include ../../../Makefile.conf
+all: CFLAGS += -O3 -DNDEBUG 
+debug: CFLAGS += -g -Wall -Werror
+
+.PHONY: libcgal
+libiglcgal: obj ../../../lib/libiglcgal.a
+
+CPP_FILES=$(wildcard *.cpp)
+OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o)))
+
+# include igl headers
+INC+=-I../../../include/
+
+# EXPECTS THAT CFLAGS IS ALREADY SET APPROPRIATELY 
+
+# Eigen dependency
+EIGEN3_INC=-I$(DEFAULT_PREFIX)/include/eigen3 -I$(DEFAULT_PREFIX)/include/eigen3/unsupported
+INC+=$(EIGEN3_INC)
+
+# Tetgen dependency
+CGAL=$(DEFAULT_PREFIX)
+CGAL_INC=-I$(CGAL)/include
+INC+=$(CGAL_INC)
+
+obj: 
+	mkdir -p obj
+
+../../../lib/libiglcgal.a: $(OBJ_FILES)
+	rm -f $@
+	ar cqs $@ $(OBJ_FILES)
+
+obj/%.o: %.cpp %.h
+	g++ $(AFLAGS) $(CFLAGS) -c -o $@ $< $(INC)
+
+clean:
+	rm -f obj/*.o
+	rm -f ../../../lib/libiglcgal.a

include/igl/cgal/SelfIntersectMesh.h

@@ -0,0 +1,931 @@
+#ifndef IGL_SELFINTERSECTMESH_H
+#define IGL_SELFINTERSECTMESH_H
+
+#include "CGAL_includes.hpp"
+#include "selfintersect.h"
+
+#include <Eigen/Dense>
+#include <list>
+#include <map>
+#include <vector>
+
+// The easiest way to keep track of everything is to use a class
+
+namespace igl
+{
+  // Kernel is a CGAL kernel like:
+  //     CGAL::Exact_predicates_inexact_constructions_kernel
+  // or 
+  //     CGAL::Exact_predicates_exact_constructions_kernel
+
+  template <typename Kernel>
+  class SelfIntersectMesh
+  {
+    public:
+      // 3D Primitives
+      typedef CGAL::Point_3<Kernel>    Point_3;
+      typedef CGAL::Segment_3<Kernel>  Segment_3; 
+      typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+      typedef CGAL::Plane_3<Kernel>    Plane_3;
+      typedef CGAL::Tetrahedron_3<Kernel> Tetrahedron_3; 
+      typedef CGAL::Polyhedron_3<Kernel> Polyhedron_3; 
+      typedef CGAL::Nef_polyhedron_3<Kernel> Nef_polyhedron_3; 
+      // 2D Primitives
+      typedef CGAL::Point_2<Kernel>    Point_2;
+      typedef CGAL::Segment_2<Kernel>  Segment_2; 
+      typedef CGAL::Triangle_2<Kernel> Triangle_2; 
+      // 2D Constrained Delaunay Triangulation types
+      typedef CGAL::Triangulation_vertex_base_2<Kernel>  TVB_2;
+      typedef CGAL::Constrained_triangulation_face_base_2<Kernel> CTFB_2;
+      typedef CGAL::Triangulation_data_structure_2<TVB_2,CTFB_2> TDS_2;
+      typedef CGAL::Exact_intersections_tag Itag;
+      typedef CGAL::Constrained_Delaunay_triangulation_2<Kernel,TDS_2,Itag> 
+        CDT_2;
+      typedef CGAL::Constrained_triangulation_plus_2<CDT_2> CDT_plus_2;
+      // Axis-align boxes for all-pairs self-intersection detection
+      typedef std::vector<Triangle_3> Triangles;
+      typedef typename Triangles::iterator TrianglesIterator;
+      typedef typename Triangles::const_iterator TrianglesConstIterator;
+      typedef 
+        CGAL::Box_intersection_d::Box_with_handle_d<double,3,TrianglesIterator> 
+        Box;
+      // Axis-aligned bounding box tree
+      typedef CGAL::AABB_triangle_primitive<Kernel,TrianglesIterator> Primitive;
+      typedef CGAL::AABB_traits<Kernel, Primitive> Traits;
+      typedef CGAL::AABB_tree<Traits> Tree;
+      typedef typename Tree::Object_and_primitive_id Object_and_primitive_id;
+      typedef typename Tree::Primitive_id Primitive_id;
+      // 3D Delaunay Triangulation
+      typedef CGAL::Delaunay_triangulation_3<Kernel> Delaunay_triangulation_3;
+      typedef typename Delaunay_triangulation_3::Finite_cells_iterator 
+        Delaunay3CellIterator;
+  
+      // Input mesh
+      const Eigen::MatrixXd & V;
+      const Eigen::MatrixXi & F;
+      // Number of self-intersecting triangle pairs
+      int count;
+      std::vector<std::list<CGAL::Object> > F_objects;
+      Triangles T;
+      std::list<int> lIF;
+      std::vector<bool> offensive;
+      std::vector<int> offending_index;
+      std::vector<int> offending;
+      // Make a short name for the edge map's key
+      typedef std::pair<int,int> EMK;
+      // Make a short name for the type stored at each edge, the edge map's
+      // value
+      typedef std::list<int> EMV;
+      // Make a short name for the edge map
+      typedef std::map<EMK,EMV> EdgeMap;
+      EdgeMap edge2faces;
+    public:
+      SelfintersectParam params;
+    public:
+      // Constructs (VV,FF) a new mesh with self-intersections of (V,F)
+      // subdivided
+      inline SelfIntersectMesh(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F,
+        const SelfintersectParam & params,
+        Eigen::MatrixXd & VV,
+        Eigen::MatrixXi & FF,
+        Eigen::MatrixXi & IF);
+    private:
+      // Helper function to mark a face as offensive
+      //
+      // Inputs:
+      //   f  index of face in F
+      inline void mark_offensive(const int f);
+      // Helper function to count intersections between faces
+      //
+      // Input:
+      //   fa  index of face A in F
+      //   fb  index of face B in F
+      inline void count_intersection(const int fa,const int fb);
+      // Helper function for box_intersect. Intersect two triangles A and B,
+      // append the intersection object (point,segment,triangle) to a running
+      // list for A and B
+      //
+      // Inputs:
+      //   A  triangle in 3D
+      //   B  triangle in 3D
+      //   fa  index of A in F (and F_objects)
+      //   fb  index of A in F (and F_objects)
+      // Returns true only if A intersects B
+      //
+      inline bool intersect(
+          const Triangle_3 & A, 
+          const Triangle_3 & B, 
+          const int fa,
+          const int fb);
+      // Helper function for box_intersect. In the case where A and B have
+      // already been identified to share a vertex, then we only want to add
+      // possible segment intersections. Assumes truly duplicate triangles are
+      // not given as input
+      //
+      // Inputs:
+      //   A  triangle in 3D
+      //   B  triangle in 3D
+      //   fa  index of A in F (and F_objects)
+      //   fb  index of B in F (and F_objects)
+      //   va  index of shared vertex in A (and F_objects)
+      //   vb  index of shared vertex in B (and F_objects)
+      //// Returns object of intersection (should be Segment or point)
+      //   Returns true if intersection (besides shared point)
+      //
+      inline bool single_shared_vertex(
+          const Triangle_3 & A,
+          const Triangle_3 & B,
+          const int fa,
+          const int fb,
+          const int va,
+          const int vb);
+      // Helper handling one direction
+      inline bool single_shared_vertex(
+          const Triangle_3 & A,
+          const Triangle_3 & B,
+          const int fa,
+          const int fb,
+          const int va);
+      // Helper function for box_intersect. In the case where A and B have
+      // already been identified to share two vertices, then we only want to add
+      // a possible coplanar (Triangle) intersection. Assumes truly degenerate
+      // facets are not givine as input.
+      inline bool double_shared_vertex(
+          const Triangle_3 & A,
+          const Triangle_3 & B,
+          const int fa,
+          const int fb);
+  
+    public:
+      // Callback function called during box self intersections test. Means
+      // boxes a and b intersect. This method then checks if the triangles in
+      // each box intersect and if so, then processes the intersections
+      //
+      // Inputs:
+      //   a  box containing a triangle
+      //   b  box containing a triangle
+      inline void box_intersect(const Box& a, const Box& b);
+    private:
+      // Compute 2D delaunay triangulation of a given 3d triangle and a list of
+      // intersection objects (points,segments,triangles). CGAL uses an affine
+      // projection rather than an isometric projection, so we're not
+      // guaranteed that the 2D delaunay triangulation here will be a delaunay
+      // triangulation in 3D.
+      //
+      // Inputs:
+      //   A  triangle in 3D
+      //   A_objects_3  updated list of intersection objects for A
+      // Outputs:
+      //   cdt  Contrained delaunay triangulation in projected 2D plane
+      inline void projected_delaunay(
+          const Triangle_3 & A,
+          const std::list<CGAL::Object> & A_objects_3,
+          CDT_plus_2 & cdt);
+      // Getters:
+    public:
+      //const std::list<int>& get_lIF() const{ return lIF;}
+      static inline void box_intersect(
+        SelfIntersectMesh * SIM, 
+        const SelfIntersectMesh::Box &a, 
+        const SelfIntersectMesh::Box &b);
+  };
+}
+
+// Implementation
+
+#include "mesh_to_cgal_triangle_list.h"
+
+#include <igl/REDRUM.h>
+#include <igl/C_STR.h>
+
+#include <boost/function.hpp>
+#include <boost/bind.hpp>
+
+#include <algorithm>
+#include <exception>
+#include <cassert>
+#include <iostream>
+
+// References:
+// http://minregret.googlecode.com/svn/trunk/skyline/src/extern/CGAL-3.3.1/examples/Polyhedron/polyhedron_self_intersection.cpp
+// http://www.cgal.org/Manual/3.9/examples/Boolean_set_operations_2/do_intersect.cpp
+
+// Q: Should we be using CGAL::Polyhedron_3?
+// A: No! Input is just a list of unoriented triangles. Polyhedron_3 requires
+// 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...
+
+//static void box_intersect(SelfIntersectMesh * SIM,const Box & A, const Box & B)
+//{
+//  return SIM->box_intersect(A,B);
+//}
+
+
+// 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
+// be bound as if the prototype was:
+//   static void box_intersect(const Box &a, const Box &b)
+// using boost:
+//  boost::function<void(const Box &a,const Box &b)> cb
+//    = boost::bind(&::box_intersect, this, _1,_2);
+//   
+template <typename Kernel>
+inline void igl::SelfIntersectMesh<Kernel>::box_intersect(
+  igl::SelfIntersectMesh<Kernel> * SIM, 
+  const typename igl::SelfIntersectMesh<Kernel>::Box &a, 
+  const typename igl::SelfIntersectMesh<Kernel>::Box &b)
+{
+  SIM->box_intersect(a,b);
+}
+
+#define FIRST_HIT_EXCEPTION 10
+
+template <typename Kernel>
+inline igl::SelfIntersectMesh<Kernel>::SelfIntersectMesh(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const SelfintersectParam & params,
+  Eigen::MatrixXd & VV,
+  Eigen::MatrixXi & FF,
+  Eigen::MatrixXi & IF):
+  V(V),
+  F(F),
+  count(0),
+  F_objects(F.rows()),
+  T(),
+  lIF(),
+  offensive(F.rows(),false),
+  offending_index(F.rows(),-1),
+  offending(),
+  edge2faces(),
+  params(params)
+{
+  using namespace std;
+  using namespace Eigen;
+  // Compute and process self intersections
+  mesh_to_cgal_triangle_list(V,F,T);
+  // http://www.cgal.org/Manual/latest/doc_html/cgal_manual/Box_intersection_d/Chapter_main.html#Section_63.5 
+  // Create the corresponding vector of bounding boxes
+  std::vector<Box> boxes;
+  boxes.reserve(T.size());
+  for ( 
+    TrianglesIterator tit = T.begin(); 
+    tit != T.end(); 
+    ++tit)
+  {
+    boxes.push_back(Box(tit->bbox(), tit));
+  }
+  // Leapfrog callback
+  boost::function<void(const Box &a,const Box &b)> cb
+    = boost::bind(&box_intersect, this, _1,_2);
+  // Run the self intersection algorithm with all defaults
+  try{
+    CGAL::box_self_intersection_d(boxes.begin(), boxes.end(),cb);
+  }catch(int e)
+  {
+    // Rethrow if not FIRST_HIT_EXCEPTION
+    if(e != FIRST_HIT_EXCEPTION)
+    {
+      throw e;
+    }
+    // Otherwise just fall through
+  }
+
+  // Convert lIF to Eigen matrix
+  assert(lIF.size()%2 == 0);
+  IF.resize(lIF.size()/2,2);
+  {
+    int i=0;
+    for(
+      typename list<int>::const_iterator ifit = lIF.begin();
+      ifit!=lIF.end();
+      )
+    {
+      IF(i,0) = (*ifit);
+      ifit++; 
+      IF(i,1) = (*ifit);
+      ifit++;
+      i++;
+    }
+  }
+
+  if(params.detect_only)
+  {
+    return;
+  }
+
+  int NF_count = 0;
+  // list of new faces, we'll fix F later
+  vector<MatrixXi> NF(offending.size());
+  // list of new vertices
+  list<Point_3> NV;
+  int NV_count = 0;
+  vector<CDT_plus_2> cdt(offending.size());
+  vector<Plane_3> P(offending.size());
+  // Use map for *all* faces
+  map<typename CDT_plus_2::Vertex_handle,int> v2i;
+  // Loop over offending triangles
+  for(int o = 0;o<(int)offending.size();o++)
+  {
+    // index in F
+    const int f = offending[o];
+    projected_delaunay(T[f],F_objects[f],cdt[o]);
+    // Q: Is this also delaunay in 3D?
+    // A: No, because the projection is affine and delaunay is not affine
+    // invariant.
+    // Q: Then, can't we first get the 2D delaunay triangulation, then lift it
+    // to 3D and flip any offending edges?
+    // Plane of projection (also used by projected_delaunay)
+    P[o] = Plane_3(T[f].vertex(0),T[f].vertex(1),T[f].vertex(2));
+    // Build index map
+    {
+      int i=0;
+      for(
+        typename CDT_plus_2::Finite_vertices_iterator vit = cdt[o].finite_vertices_begin();
+        vit != cdt[o].finite_vertices_end();
+        ++vit)
+      {
+        if(i<3)
+        {
+          //cout<<T[f].vertex(i)<<
+          //  (T[f].vertex(i) == P[o].to_3d(vit->point())?" == ":" != ")<<
+          //  P[o].to_3d(vit->point())<<endl;
+#ifndef NDEBUG
+          // I want to be sure that the original corners really show up as the
+          // original corners of the CDT. I.e. I don't trust CGAL to maintain the
+          // order
+          assert(T[f].vertex(i) == P[o].to_3d(vit->point()));
+#endif
+          // For first three, use original index in F
+          v2i[vit] = F(f,i);
+        }else
+        {
+          const Point_3 vit_point_3 = P[o].to_3d(vit->point());
+          // First look up each edge's neighbors to see if exact point has
+          // already been added (This makes everything a bit quadratic)
+          bool found = false;
+          for(int e = 0; e<3 && !found;e++)
+          {
+            // Index of F's eth edge in V
+            int i = F(f,(e+1)%3);
+            int j = F(f,(e+2)%3);
+            // Be sure that i<j
+            if(i>j)
+            {
+              swap(i,j);
+            }
+            assert(edge2faces.count(EMK(i,j))==1);
+            // loop over neighbors
+            for(
+              list<int>::const_iterator nit = edge2faces[EMK(i,j)].begin();
+              nit != edge2faces[EMK(i,j)].end() && !found;
+              nit++)
+            {
+              // don't consider self
+              if(*nit == f)
+              {
+                continue;
+              }
+              // index of neighbor in offending (to find its cdt)
+              int no = offending_index[*nit];
+              // Loop over vertices of that neighbor's cdt (might not have been
+              // processed yet, but then it's OK because it'll just be empty)
+              for(
+                  typename CDT_plus_2::Finite_vertices_iterator uit = cdt[no].finite_vertices_begin();
+                  uit != cdt[no].finite_vertices_end() && !found;
+                  ++uit)
+              {
+                if(vit_point_3 == P[no].to_3d(uit->point()))
+                {
+                  assert(v2i.count(uit) == 1);
+                  v2i[vit] = v2i[uit];
+                  found = true;
+                }
+              }
+            }
+          }
+          if(!found)
+          {
+            v2i[vit] = V.rows()+NV_count;
+            NV.push_back(vit_point_3);
+            NV_count++;
+          }
+        }
+        i++;
+      }
+    }
+    {
+      int i = 0;
+      // Resize to fit new number of triangles
+      NF[o].resize(cdt[o].number_of_faces(),3);
+      NF_count+=NF[o].rows();
+      // Append new faces to NF
+      for(
+        typename CDT_plus_2::Finite_faces_iterator fit = cdt[o].finite_faces_begin();
+        fit != cdt[o].finite_faces_end();
+        ++fit)
+      {
+        NF[o](i,0) = v2i[fit->vertex(0)];
+        NF[o](i,1) = v2i[fit->vertex(1)];
+        NF[o](i,2) = v2i[fit->vertex(2)];
+        i++;
+      }
+    }
+  }
+  assert(NV_count == (int)NV.size());
+  // Build output
+#ifndef NDEBUG
+  {
+    int off_count = 0;
+    for(int f = 0;f<F.rows();f++)
+    {
+      off_count+= (offensive[f]?1:0);
+    }
+    assert(off_count==(int)offending.size());
+  }
+#endif
+  // Append faces
+  FF.resize(F.rows()-offending.size()+NF_count,3);
+  // First append non-offending original faces
+  // There's an Eigen way to do this in one line but I forget
+  int off = 0;
+  for(int f = 0;f<F.rows();f++)
+  {
+    if(!offensive[f])
+    {
+      FF.row(off++) = F.row(f);
+    }
+  }
+  assert(off == (int)(F.rows()-offending.size()));
+  // Now append replacement faces for offending faces
+  for(int o = 0;o<(int)offending.size();o++)
+  {
+    FF.block(off,0,NF[o].rows(),3) = NF[o];
+    off += NF[o].rows();
+  }
+  // Append vertices
+  VV.resize(V.rows()+NV_count,3);
+  VV.block(0,0,V.rows(),3) = V;
+  {
+    int i = 0;
+    for(
+      typename list<Point_3>::const_iterator nvit = NV.begin();
+      nvit != NV.end();
+      nvit++)
+    {
+      for(int d = 0;d<3;d++)
+      {
+        const Point_3 & p = *nvit;
+        VV(V.rows()+i,d) = CGAL::to_double(p[d]);
+        // This distinction does not seem necessary:
+//#ifdef INEXACT_CONSTRUCTION
+//        VV(V.rows()+i,d) = CGAL::to_double(p[d]);
+//#else
+//        VV(V.rows()+i,d) = CGAL::to_double(p[d].exact());
+//#endif
+      }
+      i++;
+    }
+  }
+
+  // Q: Does this give the same result as TETGEN?
+  // A: For the cow and beast, yes.
+
+  // Q: Is tetgen faster than this CGAL implementation?
+  // A: Well, yes. But Tetgen is only solving the detection (predicates)
+  // problem. This is also appending the intersection objects (construction).
+  // But maybe tetgen is still faster for the detection part. For example, this
+  // CGAL implementation on the beast takes 98 seconds but tetgen detection
+  // takes 14 seconds
+
+}
+
+
+template <typename Kernel>
+inline void igl::SelfIntersectMesh<Kernel>::mark_offensive(const int f)
+{
+  using namespace std;
+  lIF.push_back(f);
+  if(!offensive[f])
+  {
+    offensive[f]=true;
+    offending_index[f]=offending.size();
+    offending.push_back(f);
+    // Add to edge map
+    for(int e = 0; e<3;e++)
+    {
+      // Index of F's eth edge in V
+      int i = F(f,(e+1)%3);
+      int j = F(f,(e+2)%3);
+      // Be sure that i<j
+      if(i>j)
+      {
+        swap(i,j);
+      }
+      // Create new list if there is no entry
+      if(edge2faces.count(EMK(i,j))==0)
+      {
+        edge2faces[EMK(i,j)] = list<int>();
+      }
+      // append to list
+      edge2faces[EMK(i,j)].push_back(f);
+    }
+  }
+}
+
+template <typename Kernel>
+inline void igl::SelfIntersectMesh<Kernel>::count_intersection(
+  const int fa,
+  const int fb)
+{
+  mark_offensive(fa);
+  mark_offensive(fb);
+  this->count++;
+  // We found the first intersection
+  if(params.first_only && this->count >= 1)
+  {
+    throw FIRST_HIT_EXCEPTION;
+  }
+}
+
+template <typename Kernel>
+inline bool igl::SelfIntersectMesh<Kernel>::intersect(
+  const Triangle_3 & A, 
+  const Triangle_3 & B, 
+  const int fa,
+  const int fb)
+{
+  // Determine whether there is an intersection
+  if(!CGAL::do_intersect(A,B))
+  {
+    return false;
+  }
+  if(!params.detect_only)
+  {
+    // Construct intersection
+    CGAL::Object result = CGAL::intersection(A,B);
+    F_objects[fa].push_back(result);
+    F_objects[fb].push_back(result);
+  }
+  count_intersection(fa,fb);
+  return true;
+}
+
+template <typename Kernel>
+inline bool igl::SelfIntersectMesh<Kernel>::single_shared_vertex(
+  const Triangle_3 & A,
+  const Triangle_3 & B,
+  const int fa,
+  const int fb,
+  const int va,
+  const int vb)
+{
+  ////using namespace std;
+  //CGAL::Object result = CGAL::intersection(A,B);
+  //if(CGAL::object_cast<Segment_3 >(&result))
+  //{
+  //  // Append to each triangle's running list
+  //  F_objects[fa].push_back(result);
+  //  F_objects[fb].push_back(result);
+  //  count_intersection(fa,fb);
+  //}else
+  //{
+  //  // Then intersection must be at point
+  //  // And point must be at shared vertex
+  //  assert(CGAL::object_cast<Point_3>(&result));
+  //}
+  if(single_shared_vertex(A,B,fa,fb,va))
+  {
+    return true;
+  }
+  return single_shared_vertex(B,A,fb,fa,vb);
+}
+
+template <typename Kernel>
+inline bool igl::SelfIntersectMesh<Kernel>::single_shared_vertex(
+  const Triangle_3 & A,
+  const Triangle_3 & B,
+  const int fa,
+  const int fb,
+  const int va)
+{
+  // This was not a good idea. It will not handle coplanar triangles well.
+  using namespace std;
+  Segment_3 sa(
+    A.vertex((va+1)%3),
+    A.vertex((va+2)%3));
+
+  if(CGAL::do_intersect(sa,B))
+  {
+    CGAL::Object result = CGAL::intersection(sa,B);
+    if(const Point_3 * p = CGAL::object_cast<Point_3 >(&result))
+    {
+      if(!params.detect_only)
+      {
+        // Single intersection --> segment from shared point to intersection
+        CGAL::Object seg = CGAL::make_object(Segment_3(
+          A.vertex(va),
+          *p));
+        F_objects[fa].push_back(seg);
+        F_objects[fb].push_back(seg);
+      }
+      count_intersection(fa,fb);
+      return true;
+    }else if(CGAL::object_cast<Segment_3 >(&result))
+    {
+      //cerr<<REDRUM("Coplanar at: "<<fa<<" & "<<fb<<" (single shared).")<<endl;
+      // Must be coplanar
+      if(params.detect_only)
+      {
+        count_intersection(fa,fb);
+      }else
+      {
+        // WRONG:
+        //// Segment intersection --> triangle from shared point to intersection
+        //CGAL::Object tri = CGAL::make_object(Triangle_3(
+        //  A.vertex(va),
+        //  s->vertex(0),
+        //  s->vertex(1)));
+        //F_objects[fa].push_back(tri);
+        //F_objects[fb].push_back(tri);
+        //count_intersection(fa,fb);
+        // Need to do full test. Intersection could be a general poly.
+        bool test = intersect(A,B,fa,fb);
+        ((void)test);
+        assert(test);
+      }
+      return true;
+    }else
+    {
+      cerr<<REDRUM("Segment ∩ triangle neither point nor segment?")<<endl;
+      assert(false);
+    }
+  }
+
+  return false;
+}
+
+
+template <typename Kernel>
+inline bool igl::SelfIntersectMesh<Kernel>::double_shared_vertex(
+  const Triangle_3 & A,
+  const Triangle_3 & B,
+  const int fa,
+  const int fb)
+{
+  using namespace std;
+  // Cheaper way to do this than calling do_intersect?
+  if(
+    // Can only have an intersection if co-planar
+    (A.supporting_plane() == B.supporting_plane() ||
+    A.supporting_plane() == B.supporting_plane().opposite()) &&
+    CGAL::do_intersect(A,B))
+  {
+    // Construct intersection
+    try
+    {
+      CGAL::Object result = CGAL::intersection(A,B);
+      if(result)
+      {
+        if(CGAL::object_cast<Segment_3 >(&result))
+        {
+          // not coplanar
+          return false;
+        } else if(CGAL::object_cast<Point_3 >(&result))
+        {
+          // this "shouldn't" happen but does for inexact
+          return false;
+        } else
+        {
+          if(!params.detect_only)
+          {
+            // Triangle object
+            F_objects[fa].push_back(result);
+            F_objects[fb].push_back(result);
+          }
+          count_intersection(fa,fb);
+          //cerr<<REDRUM("Coplanar at: "<<fa<<" & "<<fb<<" (double shared).")<<endl;
+          return true;
+        }
+      }else
+      {
+        // CGAL::intersection is disagreeing with do_intersect
+        return false;
+      }
+    }catch(...)
+    {
+      // This catches some cgal assertion:
+      //     CGAL error: assertion violation!
+      //     Expression : is_finite(d)
+      //     File       : /opt/local/include/CGAL/GMP/Gmpq_type.h
+      //     Line       : 132
+      //     Explanation: 
+      // But only if NDEBUG is not defined, otherwise there's an uncaught
+      // "Floating point exception: 8" SIGFPE
+      return false;
+    }
+  }
+  // Shouldn't get here either
+  assert(false);
+  return false;
+}
+
+template <typename Kernel>
+inline void igl::SelfIntersectMesh<Kernel>::box_intersect(
+  const Box& a, 
+  const Box& b)
+{
+  using namespace std;
+  // index in F and T
+  int fa = a.handle()-T.begin();
+  int fb = b.handle()-T.begin();
+  const Triangle_3 & A = *a.handle();
+  const Triangle_3 & B = *b.handle();
+  // I'm not going to deal with degenerate triangles, though at some point we
+  // should
+  assert(!a.handle()->is_degenerate());
+  assert(!b.handle()->is_degenerate());
+  // Number of combinatorially shared vertices
+  int comb_shared_vertices = 0;
+  // Number of geometrically shared vertices (*not* including combinatorially
+  // shared)
+  int geo_shared_vertices = 0;
+  // Keep track of shared vertex indices (we only handles single shared
+  // vertices as a special case, so just need last/first/only ones)
+  int va=-1,vb=-1;
+  int ea,eb;
+  for(ea=0;ea<3;ea++)
+  {
+    for(eb=0;eb<3;eb++)
+    {
+      if(F(fa,ea) == F(fb,eb))
+      {
+        comb_shared_vertices++;
+        va = ea;
+        vb = eb;
+      }else if(A.vertex(ea) == B.vertex(eb))
+      {
+        geo_shared_vertices++;
+        va = ea;
+        vb = eb;
+      }
+    }
+  }
+  const int total_shared_vertices = comb_shared_vertices + geo_shared_vertices;
+  if(comb_shared_vertices== 3)
+  {
+    // Combinatorially duplicate face, these should be removed by preprocessing
+    cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are combinatorial duplicates")<<endl;
+    goto done;
+  }
+  if(total_shared_vertices== 3)
+  {
+    // Geometrically duplicate face, these should be removed by preprocessing
+    cerr<<REDRUM("Facets "<<fa<<" and "<<fb<<" are geometrical duplicates")<<endl;
+    goto done;
+  }
+  //// SPECIAL CASES ARE BROKEN FOR COPLANAR TRIANGLES
+  //if(total_shared_vertices > 0)
+  //{
+  //  bool coplanar = 
+  //    CGAL::coplanar(A.vertex(0),A.vertex(1),A.vertex(2),B.vertex(0)) &&
+  //    CGAL::coplanar(A.vertex(0),A.vertex(1),A.vertex(2),B.vertex(1)) &&
+  //    CGAL::coplanar(A.vertex(0),A.vertex(1),A.vertex(2),B.vertex(2));
+  //  if(coplanar)
+  //  {
+  //    cerr<<MAGENTAGIN("Facets "<<fa<<" and "<<fb<<
+  //      " are coplanar and share vertices")<<endl;
+  //    goto full;
+  //  }
+  //}
+
+  if(total_shared_vertices == 2)
+  {
+    // Q: What about coplanar?
+    //
+    // o    o
+    // |\  /|
+    // | \/ |
+    // | /\ |
+    // |/  \|
+    // o----o
+    double_shared_vertex(A,B,fa,fb);
+
+    goto done;
+  }
+  assert(total_shared_vertices<=1);
+  if(total_shared_vertices==1)
+  {
+    assert(va>=0 && va<3);
+    assert(vb>=0 && vb<3);
+//#ifndef NDEBUG
+//    CGAL::Object result =
+//#endif
+    single_shared_vertex(A,B,fa,fb,va,vb);
+//#ifndef NDEBUG
+//    if(!CGAL::object_cast<Segment_3 >(&result))
+//    {
+//      const Point_3 * p = CGAL::object_cast<Point_3 >(&result);
+//      assert(p);
+//      for(int ea=0;ea<3;ea++)
+//      {
+//        for(int eb=0;eb<3;eb++)
+//        {
+//          if(F(fa,ea) == F(fb,eb))
+//          {
+//            assert(*p==A.vertex(ea));
+//            assert(*p==B.vertex(eb));
+//          }
+//        }
+//      }
+//    }
+//#endif
+  }else
+  {
+//full:
+    // No geometrically shared vertices, do general intersect
+    intersect(*a.handle(),*b.handle(),fa,fb);
+  }
+done:
+  return;
+}
+
+// Compute 2D delaunay triangulation of a given 3d triangle and a list of
+// intersection objects (points,segments,triangles). CGAL uses an affine
+// projection rather than an isometric projection, so we're not guaranteed
+// that the 2D delaunay triangulation here will be a delaunay triangulation
+// in 3D.
+//
+// Inputs:
+//   A  triangle in 3D
+//   A_objects_3  updated list of intersection objects for A
+// Outputs:
+//   cdt  Contrained delaunay triangulation in projected 2D plane
+template <typename Kernel>
+inline void igl::SelfIntersectMesh<Kernel>::projected_delaunay(
+  const Triangle_3 & A,
+  const std::list<CGAL::Object> & A_objects_3,
+  CDT_plus_2 & cdt)
+{
+  using namespace std;
+  // http://www.cgal.org/Manual/3.2/doc_html/cgal_manual/Triangulation_2/Chapter_main.html#Section_2D_Triangulations_Constrained_Plus
+  // Plane of triangle A
+  Plane_3 P(A.vertex(0),A.vertex(1),A.vertex(2));
+  // Insert triangle into vertices
+  typename CDT_plus_2::Vertex_handle corners[3];
+  for(int i = 0;i<3;i++)
+  {
+    corners[i] = cdt.insert(P.to_2d(A.vertex(i)));
+  }
+  // Insert triangle edges as constraints
+  for(int i = 0;i<3;i++)
+  {
+    cdt.insert_constraint( corners[(i+1)%3], corners[(i+2)%3]);
+  }
+  // Insert constraints for intersection objects
+  for(
+    typename list<CGAL::Object>::const_iterator lit = A_objects_3.begin();
+    lit != A_objects_3.end();
+    lit++)
+  {
+    CGAL::Object obj = *lit;
+    if(const Point_3 *ipoint = CGAL::object_cast<Point_3 >(&obj))
+    {
+      // Add point
+      cdt.insert(P.to_2d(*ipoint));
+    } else if(const Segment_3 *iseg = CGAL::object_cast<Segment_3 >(&obj))
+    {
+      // Add segment constraint
+      cdt.insert_constraint(P.to_2d(iseg->vertex(0)),P.to_2d(iseg->vertex(1)));
+    } else if(const Triangle_3 *itri = CGAL::object_cast<Triangle_3 >(&obj))
+    {
+      // Add 3 segment constraints
+      cdt.insert_constraint(P.to_2d(itri->vertex(0)),P.to_2d(itri->vertex(1)));
+      cdt.insert_constraint(P.to_2d(itri->vertex(1)),P.to_2d(itri->vertex(2)));
+      cdt.insert_constraint(P.to_2d(itri->vertex(2)),P.to_2d(itri->vertex(0)));
+    } else if(const std::vector<Point_3 > *polyp = 
+        CGAL::object_cast< std::vector<Point_3 > >(&obj))
+    {
+      //cerr<<REDRUM("Poly...")<<endl;
+      const std::vector<Point_3 > & poly = *polyp;
+      const int m = poly.size();
+      assert(m>=2);
+      for(int p = 0;p<m;p++)
+      {
+        const int np = (p+1)%m;
+        cdt.insert_constraint(P.to_2d(poly[p]),P.to_2d(poly[np]));
+      }
+    }else
+    {
+      cerr<<REDRUM("What is this object?!")<<endl;
+      assert(false);
+    }
+  }
+}
+
+#endif
+

include/igl/cgal/mesh_to_cgal_triangle_list.cpp

@@ -0,0 +1,33 @@
+#include "mesh_to_cgal_triangle_list.h"
+
+#include <cassert>
+
+template <typename Kernel>
+IGL_INLINE void igl::mesh_to_cgal_triangle_list(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  std::vector<CGAL::Triangle_3<Kernel> > & T)
+{
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  // Must be 3D
+  assert(V.cols() == 3);
+  // Must be triangles
+  assert(F.cols() == 3);
+  T.reserve(F.rows());
+  // Loop over faces
+  for(int f = 0;f<(int)F.rows();f++)
+  {
+    T.push_back(
+      Triangle_3(
+        Point_3( V(F(f,0),0), V(F(f,0),1), V(F(f,0),2)),
+        Point_3( V(F(f,1),0), V(F(f,1),1), V(F(f,1),2)),
+        Point_3( V(F(f,2),0), V(F(f,2),1), V(F(f,2),2))));
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::mesh_to_cgal_triangle_list<CGAL::Epeck>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, std::vector<CGAL::Triangle_3<CGAL::Epeck>, std::allocator<CGAL::Triangle_3<CGAL::Epeck> > >&);
+template void igl::mesh_to_cgal_triangle_list<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
+#endif

include/igl/cgal/mesh_to_cgal_triangle_list.h

@@ -0,0 +1,25 @@
+#ifndef IGL_MESH_TO_CGAL_TRIANGLE_LIST_H
+#define IGL_MESH_TO_CGAL_TRIANGLE_LIST_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+#include "CGAL_includes.hpp"
+namespace igl
+{
+  // Convert a mesh (V,F) to a list of CGAL triangles
+  //
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  // Outputs:
+  //   T  #F list of CGAL triangles
+  template <typename Kernel>
+  IGL_INLINE void mesh_to_cgal_triangle_list(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    std::vector<CGAL::Triangle_3<Kernel> > & T);
+}
+#ifdef IGL_HEADER_ONLY
+#  include "mesh_to_cgal_triangle_list.cpp"
+#endif
+
+#endif

include/igl/cgal/selfintersect.cpp

@@ -0,0 +1,41 @@
+#include "selfintersect.h"
+#include "SelfIntersectMesh.h"
+#include <igl/C_STR.h>
+#include <list>
+
+IGL_INLINE void igl::selfintersect(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const SelfintersectParam & params,
+  Eigen::MatrixXd & VV,
+  Eigen::MatrixXi & FF,
+  Eigen::MatrixXi & IF)
+{
+  using namespace std;
+  if(params.detect_only)
+  {
+    // This is probably a terrible idea, but CGAL is throwing floating point
+    // exceptions.
+
+//#ifdef __APPLE__
+//#define IGL_THROW_FPE 11
+//    const auto & throw_fpe = [](int e)
+//    {
+//      throw "IGL_THROW_FPE";
+//    };
+//    signal(SIGFPE,throw_fpe);
+//#endif
+
+    typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
+    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF);
+
+//#ifdef __APPLE__
+//    signal(SIGFPE,SIG_DFL);
+//#endif
+
+  }else
+  {
+    typedef CGAL::Exact_predicates_exact_constructions_kernel Kernel;
+    SelfIntersectMesh<Kernel> SIM = SelfIntersectMesh<Kernel>(V,F,params,VV,FF,IF);
+  }
+}

include/igl/cgal/selfintersect.h

@@ -0,0 +1,58 @@
+#ifndef IGL_SELFINTERSECT_H
+#define IGL_SELFINTERSECT_H
+#include <igl/igl_inline.h>
+
+#include <Eigen/Dense>
+
+// IGL LIB includes
+// Use header only so that asserts get guarded.
+#ifdef MEX
+//#  define IGL_HEADER_ONLY
+#  include <mex.h>
+#  include <cassert>
+#  undef assert
+#  define assert( isOK ) ( (isOK) ? (void)0 : (void) mexErrMsgTxt(C_STR(__FILE__<<":"<<__LINE__<<": failed assertion `"<<#isOK<<"'"<<std::endl) ) )
+#endif
+
+namespace igl
+{
+  // Optional Parameters
+  //   DetectOnly  Only compute IF, leave VV and FF alone
+  struct SelfintersectParam
+  {
+    bool detect_only;
+    bool first_only;
+    SelfintersectParam():detect_only(false),first_only(false){};
+  };
+  
+  // Given a triangle mesh (V,F) compute a new mesh (VV,FF) which is the same as
+  // (V,F) except that any self-intersecting triangles in (V,F) have been
+  // subdivided (new vertices and face created) so that the self-intersection
+  // contour lies exactly on edges in (VV,FF). New vertices will appear in
+  // original faces or on original edges. New vertices on edges are "merged" only
+  // across original faces sharing that edge. This means that if the input
+  // triangle mesh is a closed manifold the output will be too.
+  //
+  // Inputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices into V
+  //   params  struct of optional parameters
+  // Outputs:
+  //   VV  #VV by 3 list of vertex positions
+  //   FF  #FF by 3 list of triangle indices into V
+  //   IF  #intersecting face pairs by 2  list of intersecting face pairs,
+  //     indexing F
+  IGL_INLINE void selfintersect(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const SelfintersectParam & params,
+    Eigen::MatrixXd & VV,
+    Eigen::MatrixXi & FF,
+    Eigen::MatrixXi & IF);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "selfintersect.cpp"
+#endif
+  
+#endif

include/igl/cotangent.cpp

@@ -6,13 +6,20 @@
 // 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 "cotangent.h"
+#include "doublearea.h"
 #include "edge_lengths.h"
+#include "face_areas.h"
+#include "volume.h"
+#include "dihedral_angles.h"
 
 #include "verbose.h"
 
 
-template <class MatV, class MatF, class MatC>
-IGL_INLINE void igl::cotangent(const MatV & V, const MatF & F, MatC & C)
+template <typename DerivedV, typename DerivedF, typename DerivedC>
+IGL_INLINE void igl::cotangent(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedC>& C)
 {
   using namespace igl;
   using namespace std;
@@ -23,127 +30,72 @@ IGL_INLINE void igl::cotangent(const MatV & V, const MatF & F, MatC & C)
   int m = F.rows();
 
   // Law of cosines + law of sines
-  if(simplex_size == 3)
+  switch(simplex_size)
   {
-    // Triangles
-    //Matrix<typename MatC::Scalar,Dynamic,3> l;
-    //edge_lengths(V,F,l);
-    // 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++)
+    case 3:
     {
-      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());
+      // Triangles
+      //Matrix<typename DerivedC::Scalar,Dynamic,3> l;
+      //edge_lengths(V,F,l);
+      // edge lengths numbered same as opposite vertices
+      Matrix<typename DerivedC::Scalar,Dynamic,3> l;
+      igl::edge_lengths(V,F,l);
+      // double area
+      Matrix<typename DerivedC::Scalar,Dynamic,1> dblA;
+      doublearea(l,dblA);
+      // 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;
+      }
+      break;
     }
-    // semiperimeters
-    Matrix<typename MatC::Scalar,Dynamic,1> s = l.rowwise().sum()*0.5;
-    assert(s.rows() == m);
-    // Heron's formula for area
-    Matrix<typename MatC::Scalar,Dynamic,1> dblA(m);
-    for(int i = 0;i<m;i++)
+    case 4:
     {
-      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 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;
+      // edge lengths numbered same as opposite vertices
+      Matrix<typename DerivedC::Scalar,Dynamic,6> l;
+      edge_lengths(V,F,l);
+      Matrix<typename DerivedC::Scalar,Dynamic,4> s;
+      face_areas(l,s);
+      Matrix<typename DerivedC::Scalar,Dynamic,6> cos_theta,theta;
+      dihedral_angles_intrinsic(l,s,theta,cos_theta);
 
-    bool diag_all_pos = true;
-    C.resize(m,6);
+      // volume
+      Matrix<typename DerivedC::Scalar,Dynamic,1> vol;
+      volume(l,vol);
 
-    // 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;
+      // Law of sines
+      // http://mathworld.wolfram.com/Tetrahedron.html
+      Matrix<typename DerivedC::Scalar,Dynamic,6> sin_theta(m,6);
+      sin_theta.col(0) = vol.array() / ((2./(3.*l.col(0).array())).array() * s.col(1).array() * s.col(2).array());
+      sin_theta.col(1) = vol.array() / ((2./(3.*l.col(1).array())).array() * s.col(2).array() * s.col(0).array());
+      sin_theta.col(2) = vol.array() / ((2./(3.*l.col(2).array())).array() * s.col(0).array() * s.col(1).array());
+      sin_theta.col(3) = vol.array() / ((2./(3.*l.col(3).array())).array() * s.col(3).array() * s.col(0).array());
+      sin_theta.col(4) = vol.array() / ((2./(3.*l.col(4).array())).array() * s.col(3).array() * s.col(1).array());
+      sin_theta.col(5) = vol.array() / ((2./(3.*l.col(5).array())).array() * s.col(3).array() * s.col(2).array());
 
-      // 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;
+      // http://arxiv.org/pdf/1208.0354.pdf Page 18
+      C = (1./6.) * l.array() * cos_theta.array() / sin_theta.array();
 
-      // 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);
+      break;
     }
-    if(diag_all_pos)
+    default:
     {
-#ifdef VERBOSE 
-      verbose("cotangent.h: Flipping sign of cotangent, so that cots are positive\n");
-#endif
-      C *= -1.0;
+      fprintf(stderr,
+          "cotangent.h: Error: Simplex size (%d) not supported\n", simplex_size);
+      assert(false);
     }
-  }else
-  {
-    fprintf(stderr,
-      "cotangent.h: Error: Simplex size (%d) not supported\n", simplex_size);
-    assert(false);
   }
 }
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
-// generated by autoexplicit.sh
-template void igl::cotangent<Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-// generated by autoexplicit.sh
-template void igl::cotangent<Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-// generated by autoexplicit.sh
-template void igl::cotangent<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
-template void igl::cotangent<Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>&);
+template void igl::cotangent<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/cotangent.h

@@ -8,6 +8,7 @@
 #ifndef IGL_COTANGENT_H
 #define IGL_COTANGENT_H
 #include "igl_inline.h"
+#include <Eigen/Core>
 namespace igl
 {
   // COTANGENT compute the cotangents of each angle in mesh (V,F)
@@ -22,12 +23,17 @@ namespace igl
   // Outputs:
   //   C  #F by {3|6} list of cotangents corresponding angles
   //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
-  //     for tets, columns correspond to edges [1,2],[2,0],[0,1],[3,0],[3,1],[3,2]
+  //     for tets, columns correspond to edges
+  //     [1,2],[2,0],[0,1],[3,0],[3,1],[3,2] **times corresponding edge
+  //     lengths**
   //
   // Known bug:
   //   This computes 0.5*cotangent
-  template <class MatV, class MatF, class MatC>
-  IGL_INLINE void cotangent(const MatV & V, const MatF & F, MatC & C);
+  template <typename DerivedV, typename DerivedF, typename DerivedC>
+  IGL_INLINE void cotangent(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<DerivedC>& C);
 }
 
 #ifdef IGL_HEADER_ONLY

include/igl/cotmatrix.cpp

@@ -17,8 +17,8 @@
 
 template <typename DerivedV, typename DerivedF, typename Scalar>
 IGL_INLINE void igl::cotmatrix(
-  const Eigen::MatrixBase<DerivedV> & V, 
-  const Eigen::MatrixBase<DerivedF> & F, 
+  const Eigen::PlainObjectBase<DerivedV> & V, 
+  const Eigen::PlainObjectBase<DerivedF> & F, 
   Eigen::SparseMatrix<Scalar>& L)
 {
   using namespace igl;
@@ -81,5 +81,5 @@ IGL_INLINE void igl::cotmatrix(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
-template void igl::cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
+template void igl::cotmatrix<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&);
 #endif

include/igl/cotmatrix.h

@@ -43,8 +43,8 @@ namespace igl
   // arithmetic order in cotangent.h: C(i,e) = (arithmetic)/dblA/4
   template <typename DerivedV, typename DerivedF, typename Scalar>
   IGL_INLINE void cotmatrix(
-    const Eigen::MatrixBase<DerivedV> & V, 
-    const Eigen::MatrixBase<DerivedF> & F, 
+    const Eigen::PlainObjectBase<DerivedV> & V, 
+    const Eigen::PlainObjectBase<DerivedF> & F, 
     Eigen::SparseMatrix<Scalar>& L);
 }
 

include/igl/dihedral_angles.cpp

@@ -0,0 +1,94 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "dihedral_angles.h"
+#include <cassert>
+
+template <
+  typename DerivedV, 
+  typename DerivedT, 
+  typename Derivedtheta,
+  typename Derivedcos_theta>
+IGL_INLINE void igl::dihedral_angles(
+  Eigen::PlainObjectBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedT>& T,
+  Eigen::PlainObjectBase<Derivedtheta>& theta,
+  Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta)
+{
+  using namespace Eigen;
+  assert(T.cols() == 4);
+  Matrix<typename Derivedtheta::Scalar,Dynamic,6> l;
+  edge_lengths(V,T,l);
+  Matrix<typename Derivedtheta::Scalar,Dynamic,4> s;
+  face_areas(l,s);
+  return dihedral_angles_intrinsic(l,s,theta,cos_theta);
+}
+
+template <
+  typename DerivedL, 
+  typename DerivedA, 
+  typename Derivedtheta,
+  typename Derivedcos_theta>
+IGL_INLINE void igl::dihedral_angles_intrinsic(
+  Eigen::PlainObjectBase<DerivedL>& L,
+  Eigen::PlainObjectBase<DerivedA>& A,
+  Eigen::PlainObjectBase<Derivedtheta>& theta,
+  Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta)
+{
+  using namespace Eigen;
+  const int m = L.rows();
+  assert(m == A.rows());
+  // Law of cosines
+  // http://math.stackexchange.com/a/49340/35376
+  Matrix<typename Derivedtheta::Scalar,Dynamic,6> H_sqr(m,6);
+  H_sqr.col(0) = (1./16.) * (4. * L.col(3).array().square() * L.col(0).array().square() - 
+                                ((L.col(1).array().square() + L.col(4).array().square()) -
+                                 (L.col(2).array().square() + L.col(5).array().square())).square());
+  H_sqr.col(1) = (1./16.) * (4. * L.col(4).array().square() * L.col(1).array().square() - 
+                                ((L.col(2).array().square() + L.col(5).array().square()) -
+                                 (L.col(3).array().square() + L.col(0).array().square())).square());
+  H_sqr.col(2) = (1./16.) * (4. * L.col(5).array().square() * L.col(2).array().square() - 
+                                ((L.col(3).array().square() + L.col(0).array().square()) -
+                                 (L.col(4).array().square() + L.col(1).array().square())).square());
+  H_sqr.col(3) = (1./16.) * (4. * L.col(0).array().square() * L.col(3).array().square() - 
+                                ((L.col(4).array().square() + L.col(1).array().square()) -
+                                 (L.col(5).array().square() + L.col(2).array().square())).square());
+  H_sqr.col(4) = (1./16.) * (4. * L.col(1).array().square() * L.col(4).array().square() - 
+                                ((L.col(5).array().square() + L.col(2).array().square()) -
+                                 (L.col(0).array().square() + L.col(3).array().square())).square());
+  H_sqr.col(5) = (1./16.) * (4. * L.col(2).array().square() * L.col(5).array().square() - 
+                                ((L.col(0).array().square() + L.col(3).array().square()) -
+                                 (L.col(1).array().square() + L.col(4).array().square())).square());
+  cos_theta.resize(m,6);
+  cos_theta.col(0) = (H_sqr.col(0).array() - 
+      A.col(1).array().square() - A.col(2).array().square()).array() / 
+      (-2.*A.col(1).array() * A.col(2).array());
+  cos_theta.col(1) = (H_sqr.col(1).array() - 
+      A.col(2).array().square() - A.col(0).array().square()).array() / 
+      (-2.*A.col(2).array() * A.col(0).array());
+  cos_theta.col(2) = (H_sqr.col(2).array() - 
+      A.col(0).array().square() - A.col(1).array().square()).array() / 
+      (-2.*A.col(0).array() * A.col(1).array());
+  cos_theta.col(3) = (H_sqr.col(3).array() - 
+      A.col(3).array().square() - A.col(0).array().square()).array() / 
+      (-2.*A.col(3).array() * A.col(0).array());
+  cos_theta.col(4) = (H_sqr.col(4).array() - 
+      A.col(3).array().square() - A.col(1).array().square()).array() / 
+      (-2.*A.col(3).array() * A.col(1).array());
+  cos_theta.col(5) = (H_sqr.col(5).array() - 
+      A.col(3).array().square() - A.col(2).array().square()).array() / 
+      (-2.*A.col(3).array() * A.col(2).array());
+
+  theta = cos_theta.array().acos();
+
+  cos_theta.resize(m,6);
+
+}
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::dihedral_angles_intrinsic<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+#endif

include/igl/dihedral_angles.h

@@ -0,0 +1,55 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_DIHEDRAL_ANGLES_H
+#define IGL_DIHEDRAL_ANGLES_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // DIHEDRAL_ANGLES Compute dihedral angles for all tets of a given tet mesh
+  // (V,T)
+  //
+  // theta = dihedral_angles(V,T)
+  // theta = dihedral_angles(V,T,'ParameterName',parameter_value,...)
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   T  #V by 4 list of tet indices
+  // Outputs:
+  //   theta  #T by 6 list of dihedral angles (in radians)
+  //   cos_theta  #T by 6 list of cosine of dihedral angles (in radians)
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedT, 
+    typename Derivedtheta,
+    typename Derivedcos_theta>
+  IGL_INLINE void dihedral_angles(
+    Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedT>& T,
+    Eigen::PlainObjectBase<Derivedtheta>& theta,
+    Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta);
+  template <
+    typename DerivedL, 
+    typename DerivedA, 
+    typename Derivedtheta,
+    typename Derivedcos_theta>
+  IGL_INLINE void dihedral_angles_intrinsic(
+    Eigen::PlainObjectBase<DerivedL>& L,
+    Eigen::PlainObjectBase<DerivedA>& A,
+    Eigen::PlainObjectBase<Derivedtheta>& theta,
+    Eigen::PlainObjectBase<Derivedcos_theta>& cos_theta);
+
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "dihedral_angles.cpp"
+#endif
+
+#endif
+

include/igl/doublearea.cpp

@@ -101,6 +101,8 @@ IGL_INLINE void igl::doublearea(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::doublearea<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+// generated by autoexplicit.sh
 template void igl::doublearea<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);

include/igl/edge_lengths.cpp

@@ -6,6 +6,7 @@
 // v. 2.0. If a copy of the MPL was not distributed with this file, You can 
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "edge_lengths.h"
+#include <iostream>
 
 template <typename DerivedV, typename DerivedF, typename DerivedL>
 IGL_INLINE void igl::edge_lengths(
@@ -13,20 +14,61 @@ IGL_INLINE void igl::edge_lengths(
   const Eigen::PlainObjectBase<DerivedF>& F,
   Eigen::PlainObjectBase<DerivedL>& L)
 {
-  assert(F.cols() == 3);
-  L.resize(F.rows(),3);
-  // loop over faces
-  for(int i = 0;i<F.rows();i++)
+  using namespace std;
+  switch(F.cols())
   {
-    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());
+    case 3:
+    {
+      L.resize(F.rows(),3);
+      // loop over faces
+      for(int i = 0;i<F.rows();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());
+      }
+      break;
+    }
+    case 4:
+    {
+      const int m = F.rows();
+      L.resize(m,6);
+      // loop over faces
+      for(int i = 0;i<m;i++)
+      {
+        L(i,0) = sqrt((V.row(F(i,3))-V.row(F(i,0))).array().pow(2).sum());
+        L(i,1) = sqrt((V.row(F(i,3))-V.row(F(i,1))).array().pow(2).sum());
+        L(i,2) = sqrt((V.row(F(i,3))-V.row(F(i,2))).array().pow(2).sum());
+        L(i,3) = sqrt((V.row(F(i,1))-V.row(F(i,2))).array().pow(2).sum());
+        L(i,4) = sqrt((V.row(F(i,2))-V.row(F(i,0))).array().pow(2).sum());
+        L(i,5) = sqrt((V.row(F(i,0))-V.row(F(i,1))).array().pow(2).sum());
+      }
+      break;
+    }
+    default:
+    {
+      cerr<< "edge_lengths.h: Error: Simplex size ("<<F.cols()<<
+        ") not supported"<<endl;
+      assert(false);
+    }
   }
 }
 
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<float, -1, 2, 0, -1, 2>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+// generated by autoexplicit.sh
+template void igl::edge_lengths<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::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 template void igl::edge_lengths<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1>, Eigen::Matrix<float, -1, 3, 1, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&);
 template void igl::edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/edge_lengths.h

@@ -14,7 +14,8 @@
 namespace igl
 {
   // Constructs a list of lengths of edges opposite each index in a face
-  // (triangle) list
+  // (triangle/tet) list
+  //
   // Templates:
   //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
   //   DerivedF derived from face indices matrix type: i.e. MatrixXi
@@ -22,10 +23,14 @@ namespace igl
   // Inputs:
   //   V  eigen matrix #V by 3
   //   F  #F by 3 list of mesh faces (must be triangles)
+  //    or
+  //   T  #T by 4 list of mesh elements (must be tets)
   // Outputs:
-  //   E #E by 2 list of edges in no particular order
+  //   L  #F by {3|6} list of edge lengths 
+  //     for triangles, columns correspond to edges [1,2],[2,0],[0,1]
+  //     for tets, columns correspond to edges
+  //     [1,2],[2,0],[0,1],[3,0],[3,1],[3,2]
   //
-  // See also: adjacency_matrix
   template <typename DerivedV, typename DerivedF, typename DerivedL>
   IGL_INLINE void edge_lengths(
     const Eigen::PlainObjectBase<DerivedV>& V,

include/igl/face_areas.cpp

@@ -0,0 +1,48 @@
+#include "face_areas.h"
+#include "edge_lengths.h"
+#include "doublearea.h"
+
+template <typename DerivedV, typename DerivedT, typename DerivedA>
+IGL_INLINE void igl::face_areas(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedT>& T,
+  Eigen::PlainObjectBase<DerivedA>& A)
+{
+  assert(T.cols() == 4);
+  typename Eigen::PlainObjectBase<DerivedA> L;
+  edge_lengths(V,T,L);
+  return face_areas(L,A);
+}
+
+template <typename DerivedL, typename DerivedA>
+IGL_INLINE void igl::face_areas(
+  const Eigen::PlainObjectBase<DerivedL>& L,
+  Eigen::PlainObjectBase<DerivedA>& A)
+{
+  using namespace Eigen;
+  assert(L.cols() == 6);
+  const int m = L.rows();
+  // (unsigned) face Areas (opposite vertices: 1 2 3 4)
+  Matrix<typename DerivedA::Scalar,Dynamic,1> 
+    A0(m,1), A1(m,1), A2(m,1), A3(m,1);
+  Matrix<typename DerivedA::Scalar,Dynamic,3> 
+    L0(m,3), L1(m,3), L2(m,3), L3(m,3);
+  L0<<L.col(1),L.col(2),L.col(3);
+  L1<<L.col(0),L.col(2),L.col(4);
+  L2<<L.col(0),L.col(1),L.col(5);
+  L3<<L.col(3),L.col(4),L.col(5);
+  doublearea(L0,A0);
+  doublearea(L1,A1);
+  doublearea(L2,A2);
+  doublearea(L3,A3);
+  A.resize(m,4);
+  A.col(0) = 0.5*A0;
+  A.col(1) = 0.5*A1;
+  A.col(2) = 0.5*A2;
+  A.col(3) = 0.5*A3;
+}
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::face_areas<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 4, 0, -1, 4> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 4, 0, -1, 4> >&);
+#endif

include/igl/face_areas.h

@@ -0,0 +1,46 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_FACE_AREAS_H
+#define IGL_FACE_AREAS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+
+namespace igl
+{
+  // Constructs a list of face areas of faces opposite each index in a tet list
+  //
+  // Templates:
+  //   DerivedV derived from vertex positions matrix type: i.e. MatrixXd
+  //   DerivedT derived from face indices matrix type: i.e. MatrixXi
+  //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
+  // Inputs:
+  //   V  eigen matrix #V by 3
+  //   T  #T by 3 list of mesh faces (must be triangles)
+  // Outputs:
+  //   E #E by 2 list of edges in no particular order
+  //
+  // See also: adjacency_matrix
+  template <typename DerivedV, typename DerivedT, typename DerivedA>
+  IGL_INLINE void face_areas(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedT>& T,
+    Eigen::PlainObjectBase<DerivedA>& A);
+  template <typename DerivedL, typename DerivedA>
+  IGL_INLINE void face_areas(
+    const Eigen::PlainObjectBase<DerivedL>& L,
+    Eigen::PlainObjectBase<DerivedA>& A);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "face_areas.cpp"
+#endif
+
+#endif
+
+

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
+

include/igl/pos.h

@@ -5,6 +5,7 @@
 // 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_POS_H
 #define IGL_POS_H
 
@@ -29,30 +30,13 @@ template <typename S>
         )
     : F(F), FF(FF), FFi(FFi), fi(fi), ei(ei), reverse(reverse)
     {}
-
-//    // Init the pos by specifying Face,Vertex Index and Orientation
-//    Pos(Eigen::MatrixXi& F, 
-//        Eigen::MatrixXi& FF, 
-//        Eigen::MatrixXi& FFi, 
-//        int fi,
-//        int vi,
-//        bool reverse = false
-//        )
-//    : F(F), FF(FF), FFi(FFi), fi(fi), reverse(reverse)
-//    {
-//      ei = -1;
-//      for (int i=0;i<3;++i)
-//        if (F(fi,i) == vi)
-//          ei = i;
-//      assert(ei != -1);
-//      
-//      if (reverse)
-//        ei = (ei-1)%3;
-//    }
     
     // Change Face
     void flipF()
     {
+      if (isBorder())
+        return;
+      
       int fin = (*FF)(fi,ei);
       int ein = (*FFi)(fi,ei);
       int reversen = !reverse;
@@ -79,6 +63,42 @@ template <typename S>
       reverse = !reverse;
     }
     
+    bool isBorder()
+    {
+      return (*FF)(fi,ei) == -1;
+    }
+    
+    /*!
+     * Returns the next edge skipping the border
+     *      _________
+     *     /\ c | b /\
+     *    /  \  |  /  \
+     *   / d  \ | / a  \
+     *  /______\|/______\
+     *          v
+     * In this example, if a and d are of-border and the pos is iterating counterclockwise, this method iterate through the faces incident on vertex v,
+     * producing the sequence a, b, c, d, a, b, c, ...
+     */
+    bool NextFE()
+    {
+      if ( isBorder() ) // we are on a border
+      {
+        do
+        {
+					flipF();
+					flipE();
+        } while (!isBorder());
+        flipE();
+        return false;
+      }
+      else
+      {
+        flipF();
+        flipE();
+        return true;
+      }
+    }
+    
     // Get vertex index
     int Vi()
     {
@@ -98,6 +118,13 @@ template <typename S>
     {
       return fi;
     }
+
+    // Get edge index
+    int Ei()
+    {
+      return ei;
+    }
+
     
     bool operator==(Pos& p2)
     {
@@ -112,13 +139,14 @@ template <typename S>
        );
     }
     
+  private:
     int fi;
     int ei;
     bool reverse;
     
     const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic>*     F;
-    Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic>*     FF;
-    Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic>*     FFi;
+    const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic>*     FF;
+    const Eigen::Matrix<S, Eigen::Dynamic, Eigen::Dynamic>*     FFi;
   };
   
 }

include/igl/sort.cpp

@@ -213,6 +213,12 @@ IGL_INLINE void igl::sort(
 #ifndef IGL_HEADER_ONLY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<float, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<float, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<float, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
+template void igl::sort<Eigen::Matrix<double, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::sort<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::sort<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > >(std::vector<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > > const&, bool, std::vector<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> >, std::allocator<SortableRow<Eigen::Matrix<int, -1, 1, 0, -1, 1> > > >&, std::vector<unsigned long, std::allocator<unsigned long> >&);
@@ -225,4 +231,5 @@ template void igl::sort<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<in
 template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::sort<Eigen::Matrix<int, -1, 2, 0, -1, 2>, Eigen::Matrix<int, -1, 2, 0, -1, 2> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 2, 0, -1, 2> >&);
+template void igl::sort<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/vf.h

@@ -14,14 +14,21 @@
 
 namespace igl 
 {
-  // Constructs the vertex-face topology of a given mesh (V,F)
+  // TODO: This function is badly named.
+  //
+  // VF Constructs the vertex-face topology of a given mesh (V,F)
+  //
   // Inputs:
   //   V  #V by 3 list of vertex coordinates
   //   F  #F by dim list of mesh faces (must be triangles)
   // Outputs: 
-  // 
+  //   VF  #V list of lists of incident faces (adjacency list)
+  //   VI  #V list of lists of index of incidence within incident faces listed
+  //     in VF
   //
   // See also: edges, cotmatrix, diag, vv
+  //
+  // Known bugs: this should not take V as an input parameter.
     
   template <typename DerivedV, typename DerivedF, typename IndexType>
   IGL_INLINE void vf(

include/igl/volume.cpp

@@ -0,0 +1,72 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "volume.h"
+template <
+  typename DerivedV, 
+  typename DerivedT, 
+  typename Derivedvol>
+IGL_INLINE void igl::volume(
+  Eigen::PlainObjectBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedT>& T,
+  Eigen::PlainObjectBase<Derivedvol>& vol)
+{
+  using namespace Eigen;
+  const int m = T.rows();
+  vol.resize(m,1);
+  for(int t = 0;t<m;t++)
+  {
+    const auto & a = V.row(T(t,0));
+    const auto & b = V.row(T(t,1));
+    const auto & c = V.row(T(t,2));
+    const auto & d = V.row(T(t,3));
+    vol(t) = -(a-d).dot((b-d).cross(c-d))/6.;
+  }
+}
+
+template <
+  typename DerivedL, 
+  typename Derivedvol>
+IGL_INLINE void igl::volume(
+  Eigen::PlainObjectBase<DerivedL>& L,
+  Eigen::PlainObjectBase<Derivedvol>& vol)
+{
+  using namespace Eigen;
+  const int m = L.rows();
+  typedef typename Derivedvol::Scalar ScalarS;
+  vol.resize(m,1);
+  for(int t = 0;t<m;t++)
+  {
+    const ScalarS u = L(t,0);
+    const ScalarS v = L(t,1);
+    const ScalarS w = L(t,2);
+    const ScalarS U = L(t,3);
+    const ScalarS V = L(t,4);
+    const ScalarS W = L(t,5);
+    const ScalarS X = (w - U + v)*(U + v + w);
+    const ScalarS x = (U - v + w)*(v - w + U);
+    const ScalarS Y = (u - V + w)*(V + w + u);
+    const ScalarS y = (V - w + u)*(w - u + V);
+    const ScalarS Z = (v - W + u)*(W + u + v);
+    const ScalarS z = (W - u + v)*(u - v + W);
+    const ScalarS a = sqrt(x*Y*Z); 
+    const ScalarS b = sqrt(y*Z*X); 
+    const ScalarS c = sqrt(z*X*Y); 
+    const ScalarS d = sqrt(x*y*z); 
+    vol(t) = sqrt(
+       (-a + b + c + d)*
+       ( a - b + c + d)*
+       ( a + b - c + d)*
+       ( a + b + c - d))/
+       (192.*u*v*w);
+  }
+}
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::volume<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/volume.h

@@ -0,0 +1,51 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2014 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_VOLUME_H
+#define IGL_VOLUME_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // VOLUME Compute volume for all tets of a given tet mesh
+  // (V,T)
+  //
+  // vol = volume(V,T)
+  //
+  // Inputs:
+  //   V  #V by dim list of vertex positions
+  //   T  #V by 4 list of tet indices
+  // Outputs:
+  //   vol  #T list of dihedral angles (in radians)
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedT, 
+    typename Derivedvol>
+  IGL_INLINE void volume(
+    Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedT>& T,
+    Eigen::PlainObjectBase<Derivedvol>& vol);
+  // Intrinsic version:
+  //
+  // Inputs:
+  //   L  #V by 6 list of edge lengths (see edge_lengths)
+  template <
+    typename DerivedL, 
+    typename Derivedvol>
+  IGL_INLINE void volume(
+    Eigen::PlainObjectBase<DerivedL>& L,
+    Eigen::PlainObjectBase<Derivedvol>& vol);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "volume.cpp"
+#endif
+
+#endif
+
+

include/igl/xml/XMLSerializable.h

@@ -16,7 +16,7 @@
 #define XML_SERIALIZABLE_H
 
 #include <iostream>
-#include "tinyxml2.h"
+#include <tinyxml2.h>
 
 namespace igl
 {

include/igl/xml/XMLSerialization.h

@@ -152,7 +152,7 @@ namespace igl
         Name = obj.Name;
         if(xmlSerializer != NULL)
         {
-          delete obj.xmlSerializer;
+          delete xmlSerializer;
           xmlSerializer = NULL;
         }
       }

include/igl/xml/XMLSerializer.h

@@ -1,2187 +0,0 @@
-// This file is part of libigl, a simple c++ geometry processing library.
-// 
-// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
-// 
-// This Source Code Form is subject to the terms of the Mozilla Public License 
-// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
-// obtain one at http://mozilla.org/MPL/2.0/.
-/* ---------------------------------------------------------------------------
- // XMLSerializer.h
- // Author: Christian Schüller on 08/05/13.
- ------------------------------------------------------------------------------
- 
- Small neader library which allows to save and load a serialization of basic c++ data
- types like char, char*, std::string, bool, uint, int, float, double to and from a xml file.
- Containers like std::vector, std::std::pair, Eigen dense and sparse matrices are supported
- as well as combinations of them (like vector<pair<string,bool>> or vector<vector<int>>).
- 
- To serialize an arbitary object use the XMLSerializable interface or even simpler the
- XMLSerialization class.
- 
- The serialized objects are organised in groups in the xml file and have
- their own names which must be unique within one group.
- 
- You can find examples how to use it in the test case class XMLSerializerTest.
- 
- ----------------------------------------------------------------------------*/
-#ifndef XML_SERIALIZER_H
-#define XML_SERIALIZER_H
-
-#include <igl/igl_inline.h>
-#include <igl/xml/XMLSerializable.h>
-
-#include <iostream>
-//#include <array>
-#include <vector>
-#include <set>
-#include <map>
-
-#include <Eigen/Dense>
-#include <Eigen/Sparse>
-
-#include "tinyxml2.h"
-
-namespace igl
-{
-  namespace
-  {
-    // Utility functions
-    void EncodeXMLElementName(std::string& name);
-    void DecodeXMLElementName(std::string& name);
-    void ReplaceSubString(std::string& str, const std::string& search, const std::string& replace);
-    
-    // Forward declaration
-    class XMLSerializer;
-    
-    /**
-     * class XMLSerializableObject
-     * internal usage
-     */
-    class XMLSerializableObject : public ::igl::XMLSerializable
-    {
-    public:
-      
-      XMLSerializableObject(const std::string& name, const std::string& group);
-      virtual ~XMLSerializableObject();
-      
-      // set attribute conversion functions
-      void SetAttribute(tinyxml2::XMLElement* element, const char* name, char& dest);
-      void SetAttribute(tinyxml2::XMLElement* element, const char* name, char*& dest);
-      void SetAttribute(tinyxml2::XMLElement* element, const char* name, std::string& dest);
-      void SetAttribute(tinyxml2::XMLElement* element, const char* name, bool& dest);
-      void SetAttribute(tinyxml2::XMLElement* element, const char* name, unsigned int& dest);
-      void SetAttribute(tinyxml2::XMLElement* element, const char* name, int& dest);
-      void SetAttribute(tinyxml2::XMLElement* element, const char* name, float& dest);
-      void SetAttribute(tinyxml2::XMLElement* element, const char* name, double& dest);
-      
-      // get attribute conversion functions
-      void GetAttribute(const char* src, char& dest);
-      void GetAttribute(const char* src, char*& dest);
-      void GetAttribute(const char* src, std::string& dest);
-      void GetAttribute(const char* src, bool& dest);
-      void GetAttribute(const char* src, unsigned int& dest);
-      void GetAttribute(const char* src, int& dest);
-      void GetAttribute(const char* src, float& dest);
-      void GetAttribute(const char* src, double& dest);
-      
-      // Initialization
-      
-      // Basic data types
-      using XMLSerializable::Init;
-      void Init(char& val);
-      void Init(char*& val);
-      void Init(std::string& val);
-      void Init(bool& val);
-      void Init(unsigned int& val);
-      void Init(int& val);
-      void Init(float& val);
-      void Init(double& val);
-      
-      // XMLSerializable*
-      template<typename T>
-      void Init(T& obj);
-      template<typename T>
-      void Init(T*& obj);
-      
-      // STL containers
-      /*template<typename T, int S>
-       void Init(std::array<T,S>& obj);*/
-      template<typename T0, typename T1>
-      void Init(std::pair<T0,T1>& obj);
-      template<typename T>
-      void Init(std::vector<T>& obj);
-      template<typename T>
-      void Init(std::set<T>& obj);
-      template<typename T0, typename T1>
-      void Init(std::map<T0,T1>& obj);
-      
-      // Eigen types
-      template<typename T, int R, int C>
-      void Init(Eigen::Matrix<T,R,C>& obj);
-      template<typename T>
-      void Init(Eigen::SparseMatrix<T>& obj);
-      
-      // Serialization
-      
-      // Basic data types
-      using XMLSerializable::Serialize;
-      bool Serialize(char& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(char*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(std::string& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(std::string*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(bool obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(bool*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(unsigned int& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(unsigned int*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(int& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(int*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(float& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(float*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(double& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(double*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      
-      // XMLSerializable*
-      template<typename T>
-      bool Serialize(T& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool Serialize(T*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      
-      // STL containers
-      
-      /*template<typename T, size_t S>
-       bool Serialize(std::array<T,S>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-       template<typename T, size_t S>
-       bool Serialize(std::array<T,S>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);*/
-      
-      template<typename T0, typename T1>
-      bool Serialize(std::pair<T0,T1>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T0, typename T1>
-      bool Serialize(std::pair<T0,T1>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      
-      template<typename T>
-      bool Serialize(std::vector<T>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool Serialize(std::vector<T>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(std::vector<bool>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool Serialize(std::vector<bool>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-
-      template<typename T>
-      bool Serialize(std::set<T>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool Serialize(std::set<T>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-
-      template<typename T0, typename T1>
-      bool Serialize(std::map<T0,T1>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T0, typename T1>
-      bool Serialize(std::map<T0,T1>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      
-      // Eigen types
-      template<typename T, int R, int C>
-      bool Serialize(Eigen::Matrix<T,R,C>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T, int R, int C>
-      bool Serialize(Eigen::Matrix<T,R,C>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      
-      template<typename T>
-      bool Serialize(Eigen::SparseMatrix<T>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool Serialize(Eigen::SparseMatrix<T>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      
-      // Deserialization
-      
-      // Basic data types
-      using XMLSerializable::Deserialize;
-      bool Deserialize(char& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(char*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(std::string& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(std::string*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(bool& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(bool*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(unsigned int& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(unsigned int*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(int& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(int*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(float& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(float*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(double& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(double*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      
-      // XMLSerializable*
-      template<typename T>
-      bool Deserialize(T& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool Deserialize(T*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      
-      // STL containers
-      
-      /*template<typename T, size_t S>
-       bool Deserialize(std::array<T,S>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-       template<typename T, size_t S>
-       bool Deserialize(std::array<T,S>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);*/
-      
-      template<typename T0, typename T1>
-      bool Deserialize(std::pair<T0,T1>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T0, typename T1>
-      bool Deserialize(std::pair<T0,T1>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      
-      template<typename T>
-      bool Deserialize(std::vector<T>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool Deserialize(std::vector<T>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(std::vector<bool>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      bool Deserialize(std::vector<bool>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-
-      template<typename T>
-      bool Deserialize(std::set<T>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool Deserialize(std::set<T>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-
-      template<typename T0, typename T1>
-      bool Deserialize(std::map<T0,T1>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T0, typename T1>
-      bool Deserialize(std::map<T0,T1>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      
-      // Eigen types
-      template<typename T, int R, int C>
-      bool Deserialize(Eigen::Matrix<T,R,C>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T, int R, int C>
-      bool Deserialize(Eigen::Matrix<T,R,C>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      
-      template<typename T>
-      bool Deserialize(Eigen::SparseMatrix<T>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool Deserialize(Eigen::SparseMatrix<T>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-      
-    private:
-      
-      template<typename T>
-      bool setElementAttribute(T& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      template<typename T>
-      bool getElementAttribute(T& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name);
-    };
-    
-    /**
-     * class XMLSerializableInstance
-     * internal usage
-     */
-    template<typename T>
-    class XMLSerializableInstance : public XMLSerializableObject
-    {
-    public:
-      
-      T& Object;
-      T DefaultValue;
-      
-      XMLSerializableInstance(T& obj, const std::string& name, const std::string group);
-      XMLSerializableInstance(T& obj, const std::string& name, const std::string group, T defaultValue);
-      ~XMLSerializableInstance();
-      
-      // XMLSerializable interface implementation
-      void Init();
-      bool Serialize(tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element);
-      bool Deserialize(tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element);
-    };
-    
-    /**
-     * struct XMLSerializerGroup
-     * internal usage
-     */
-    struct XMLSerializerGroup
-    {
-      std::string Name;
-      std::vector<XMLSerializable*>* Objects;
-    };
-    
-    /**
-     * class XMLSerializer
-     * This is the core class which takes care of saving and loading of serialization of object structures.
-     */
-    class XMLSerializer
-    {
-    public:
-      
-      /**
-       * Serializes an object to a file
-       */
-      template<typename T>
-      static bool SaveObject(T& object, const char* filename);
-      template<typename T>
-      static bool SaveObject(T& object, const std::string& objectName, const std::string& groupName, const char* filename, bool overwrite);
-      
-      /**
-       * Loads the serialization of an object from a file.
-       */
-      template<typename T>
-      static bool LoadObject(T& object, const char* filename);
-      template<typename T>
-      static bool LoadObject(T& object, const std::string& objectName, const std::string& groupName, const char* filename);
-      
-      /**
-       * Constructor which sets the default group
-       */
-      XMLSerializer(const std::string& defaultGroup);
-      ~XMLSerializer();
-      
-      /**
-       * Save the serialization of all groups to file.
-       * Parameter overwrite specifies if file gets overwritten or updated
-       */
-      bool Save(const char* filename, bool overwrite);
-      bool Save(const std::string& groupName, const char* filename, bool overwrite);
-      bool Save(const std::string& objectName, const std::string& groupName, const char* filename, bool overwrite);
-      
-      /**
-       * Save the serialization of all groups to a XMLDocument instance.
-       */
-      bool SaveToXMLDoc(tinyxml2::XMLDocument* doc);
-      bool SaveToXMLDoc(const std::string& groupName, tinyxml2::XMLDocument* doc);
-      bool SaveToXMLDoc(const std::string& objectName, const std::string& groupName, tinyxml2::XMLDocument* doc);
-      
-      /**
-       * Save the serialization of a group with a new provided name to given XMLElement instance.
-       */
-      bool SaveGroupToXMLElement(tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      bool SaveGroupToXMLElement(const std::string& groupName, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name);
-      
-      /**
-       * Load the serialization from a file.
-       */
-      bool Load(const char* filename);
-      bool Load(const std::string& groupName, const char* filename);
-      bool Load(const std::string& objectName, const std::string& groupName, const char* filename);
-      
-      /**
-       * Load the serialization from an XMLDocument instance.
-       */
-      bool LoadFromXMLDoc(tinyxml2::XMLDocument* doc);
-      bool LoadFromXMLDoc(const std::string& groupName, tinyxml2::XMLDocument* doc);
-      bool LoadFromXMLDoc(const std::string& objectName, const std::string& groupName, tinyxml2::XMLDocument* doc);
-      
-      /**
-       * Load the serialization from a XMLElement instance to given group.
-       */
-      bool LoadGroupFromXMLElement(const std::string& groupName, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element);
-      
-      /**
-       * Set/Get current group. Every object which is added afterwards will be in this group, except it specifies another group.
-       */
-      void SetCurrentGroup(const std::string& group);
-      std::string GetCurrentGroup();
-      
-      /**
-       * Add an object to the serializer. Can be simple types like char, char*, string, unsigned int, int, float, double or containers like std::array, std::pair, std::vector.
-       * Also Eigen dense or sparse matrices are supported and all objects of type Serializable* and combinations of thoses types like vector<vector>, vector<pair> or even vector<pair<vector,Serializable*>>>.
-       * Also pointers to those objects can be used (for instance like vector<vector<pair<int,float>*>*>).
-       * char* is also possible as base type and represents a array of chars, but be carefull that the pointer is not just a copy but a valid instance in the current programm scope.
-       */
-      
-      // Basic types
-      bool Add(char& obj, const std::string& name);
-      bool Add(char*& obj, const std::string& name);
-      bool Add(std::string& obj, const std::string& name);
-      bool Add(bool& obj, const std::string& name);
-      bool Add(unsigned int& obj, const std::string& name);
-      bool Add(int& obj, const std::string& name);
-      bool Add(float& obj, const std::string& name);
-      bool Add(double& obj, const std::string& name);
-      
-      bool Add(char& obj, const std::string& name, char defaultValue);
-      bool Add(char*& obj, const std::string& name, char* defaultValue);
-      bool Add(std::string& obj, const std::string& name, std::string defaultValue);
-      bool Add(bool& obj, const std::string& name, bool defaultValue);
-      bool Add(unsigned int& obj, const std::string& name, unsigned int defaultValue);
-      bool Add(int& obj, const std::string& name, int defaultValue);
-      bool Add(float& obj, const std::string& name, float defaultValue);
-      bool Add(double& obj, const std::string& name, double defaultValue);
-      
-      // XMLSerializable*
-      template<typename T>
-      bool Add(T& object, const std::string& name);
-      template<typename T>
-      bool Add(T& object, const std::string& name, T defaultValue);
-      
-      // STL containers
-      /*template<typename T, size_t S>
-       bool Add(std::array<T,S>& obj, const std::string& name);*/
-      template<typename T0, typename T1>
-      bool Add(std::pair<T0,T1>& obj, const std::string& name);
-      template<typename T>
-      bool Add(std::vector<T>& obj, const std::string& name);
-      template<typename T>
-      bool Add(std::set<T>& obj, const std::string& name);
-      template<typename T0, typename T1>
-      bool Add(std::map<T0,T1>& obj, const std::string& name);
-      
-      // Eigen types
-      template<typename T, int R, int C>
-      bool Add(Eigen::Matrix<T,R,C>& obj, const std::string& name);
-      template<typename T>
-      bool Add(Eigen::SparseMatrix<T>& obj, const std::string& name);
-      
-    private:
-      
-      std::map<std::string,XMLSerializerGroup*>::iterator currentGroup;
-      std::map<std::string,XMLSerializerGroup*> groups;
-      
-      template<typename T>
-      bool add(T& object, const std::string& name);
-      template<typename T>
-      bool add(T& object, const std::string& name, T defaultValue);
-      bool addObjectToGroup(XMLSerializable* object, const std::string& group);
-      bool addObjectToGroup(XMLSerializable* object, std::map<std::string,XMLSerializerGroup*>::iterator it);
-      std::map<std::string,XMLSerializerGroup*>::iterator setGetGroup(const std::string& group);
-      tinyxml2::XMLDocument* openDoc(const char* filename);
-      tinyxml2::XMLElement* findAddGroup(tinyxml2::XMLDocument* doc, const char* groupName);
-    };
-    
-    int numForbiddenChars = 8;
-    char forbiddenChars[] = {' ','/','~','#','&','>','<','='};
-    
-    void ReplaceSubString(std::string& str, const std::string& search, const std::string& replace)
-    {
-      size_t pos = 0;
-      while ((pos = str.find(search, pos)) != std::string::npos)
-      {
-        str.replace(pos, search.length(), replace);
-        pos += replace.length();
-      }
-    }
-    
-    void EncodeXMLElementName(std::string& name)
-    {
-      // must not start with a digit
-      if(isdigit(*name.begin()))
-      {
-        name = ":::" + name;
-      }
-      
-      std::stringstream stream;
-      for(int i=0;i<numForbiddenChars;i++)
-      {
-        std::string search;
-        search = forbiddenChars[i];
-        std::stringstream replaces;
-        replaces << ":" << (int)forbiddenChars[i];
-        std::string replace = replaces.str();
-        
-        ReplaceSubString(name,search,replace);
-      }
-    }
-    
-    void DecodeXMLElementName(std::string& name)
-    {
-      if(name.find("::", 0) == 0)
-        name.replace(0,3,"");
-      
-      std::stringstream stream;
-      for(unsigned int i=0;i<numForbiddenChars;i++)
-      {
-        std::stringstream searchs;
-        searchs << ":" << (int)forbiddenChars[i];
-        std::string search = searchs.str();
-        std::string replace;
-        replace = forbiddenChars[i];
-        
-        ReplaceSubString(name,search,replace);
-      }
-    }
-    
-    XMLSerializableObject::XMLSerializableObject(const std::string& name, const std::string& group)
-    {
-      std::string groupName = group;
-      std::string objectName = name;
-      
-      EncodeXMLElementName(groupName);
-      EncodeXMLElementName(objectName);
-      
-      Name = objectName;
-    }
-    
-    XMLSerializableObject::~XMLSerializableObject()
-    {
-    }
-    
-    // set attribute conversion functions
-    void XMLSerializableObject::SetAttribute(tinyxml2::XMLElement* element, const char* name, char& dest)
-    {
-      element->SetAttribute(name,dest);
-    }
-    
-    void XMLSerializableObject::SetAttribute(tinyxml2::XMLElement* element, const char* name, char*& dest)
-    {
-      element->SetAttribute(name,const_cast<const char*>(dest));
-    }
-    
-    void XMLSerializableObject::SetAttribute(tinyxml2::XMLElement* element, const char* name, std::string& dest)
-    {
-      element->SetAttribute(name,dest.c_str());
-    }
-    
-    void XMLSerializableObject::SetAttribute(tinyxml2::XMLElement* element, const char* name, bool& dest)
-    {
-      element->SetAttribute(name,dest);
-    }
-    
-    void XMLSerializableObject::SetAttribute(tinyxml2::XMLElement* element, const char* name, unsigned int& dest)
-    {
-      element->SetAttribute(name,dest);
-    }
-    
-    void XMLSerializableObject::SetAttribute(tinyxml2::XMLElement* element, const char* name, int& dest)
-    {
-      element->SetAttribute(name,dest);
-    }
-    
-    void XMLSerializableObject::SetAttribute(tinyxml2::XMLElement* element, const char* name, float& dest)
-    {
-      element->SetAttribute(name,dest);
-    }
-    
-    void XMLSerializableObject::SetAttribute(tinyxml2::XMLElement* element, const char* name, double& dest)
-    {
-      element->SetAttribute(name,dest);
-    }
-    
-    // get attribute conversion functions
-    void XMLSerializableObject::GetAttribute(const char* src, char& dest)
-    {
-      dest = (char)atoi(src);
-    }
-    
-    void XMLSerializableObject::GetAttribute(const char* src, char*& dest)
-    {
-      unsigned int length = strlen(src)+1;
-      dest = new char[length];
-      strcpy(dest, src);
-    }
-    
-    void XMLSerializableObject::GetAttribute(const char* src, std::string& dest)
-    {
-      dest = src;
-    }
-    
-    void XMLSerializableObject::GetAttribute(const char* src, bool& dest)
-    {
-      tinyxml2::XMLUtil::ToBool(src,&dest);
-    }
-    
-    void XMLSerializableObject::GetAttribute(const char* src, unsigned int& dest)
-    {
-      tinyxml2::XMLUtil::ToUnsigned(src,&dest);
-    }
-    
-    void XMLSerializableObject::GetAttribute(const char* src, int& dest)
-    {
-      tinyxml2::XMLUtil::ToInt(src,&dest);
-    }
-    
-    void XMLSerializableObject::GetAttribute(const char* src, float& dest)
-    {
-      tinyxml2::XMLUtil::ToFloat(src,&dest);
-    }
-    
-    void XMLSerializableObject::GetAttribute(const char* src, double& dest)
-    {
-      tinyxml2::XMLUtil::ToDouble(src,&dest);
-    }
-    
-    // specify default value of types
-    void XMLSerializableObject::Init(char& val)
-    {
-      val = '0';
-    }
-    
-    void XMLSerializableObject::Init(char*& val)
-    {
-      val = NULL;
-    }
-    
-    void XMLSerializableObject::Init(std::string& val)
-    {
-      val = "";
-    }
-    
-    void XMLSerializableObject::Init(bool& val)
-    {
-      val = false;
-    }
-    
-    void XMLSerializableObject::Init(unsigned int& val)
-    {
-      val = 0;
-    }
-    
-    void XMLSerializableObject::Init(int& val)
-    {
-      val = 0;
-    }
-    
-    void XMLSerializableObject::Init(float& val)
-    {
-      val = 0.0f;
-    }
-    
-    void XMLSerializableObject::Init(double& val)
-    {
-      val = 0.000000000000000;
-    }
-
-    template<typename T>
-    void XMLSerializableObject::Init(T& obj)
-    {
-      XMLSerializable* object = static_cast<XMLSerializable*>(&obj);
-      object->Init();
-    }
-    
-    template<typename T>
-    void XMLSerializableObject::Init(T*& obj)
-    {
-      XMLSerializable* object = static_cast<XMLSerializable*>(obj);
-      object->Init();
-    }
-    
-    /*template<typename T, int S>
-     void XMLSerializableObject::Init(std::array<T,S>& obj)
-     {
-     for(unsigned int i=0;i<obj.size();i++)
-     Init(obj[i]);
-     }*/
-    
-    template<typename T0, typename T1>
-    void XMLSerializableObject::Init(std::pair<T0,T1>& obj)
-    {
-      Init(obj.first);
-      Init(obj.second);
-    }
-    
-    template<typename T>
-    void XMLSerializableObject::Init(std::vector<T>& obj)
-    {
-      obj.clear();
-    }
-
-    template<typename T>
-    void XMLSerializableObject::Init(std::set<T>& obj)
-    {
-      obj.clear();
-    }
-
-    template<typename T0, typename T1>
-    void XMLSerializableObject::Init(std::map<T0,T1>& obj)
-    {
-      obj.clear();
-    }
-    
-    template<typename T, int R, int C>
-    void XMLSerializableObject::Init(Eigen::Matrix<T,R,C>& obj)
-    {
-      obj.setZero(obj.rows(),obj.cols());
-    }
-    
-    template<typename T>
-    void  XMLSerializableObject::Init(Eigen::SparseMatrix<T>& obj)
-    {
-      obj.setZero();
-    }
-    
-    bool XMLSerializableObject::Serialize(char& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return setElementAttribute(obj,doc,element,name);
-    }
-    
-    // overload function for char*, it interpreted as char array and can be used to handle strings
-    bool XMLSerializableObject::Serialize(char*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return setElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(std::string& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return setElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(std::string*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(bool obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return setElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(bool*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(unsigned int& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return setElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(unsigned int*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(int& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return setElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(int*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(float& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return setElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(float*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(double& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return setElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Serialize(double*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Serialize(T& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      // Serialize object implementing XMLSerializable interface
-      XMLSerializable* object = static_cast<XMLSerializable*>(&obj);
-      
-      tinyxml2::XMLElement* child = doc->NewElement(name.c_str());
-      element->InsertEndChild(child);
-      
-      return object->Serialize(doc,child);
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Serialize(T*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      // Serialize object implementing XMLSerializable interface
-      XMLSerializable* object = static_cast<XMLSerializable*>(obj);
-      
-      tinyxml2::XMLElement* child = doc->NewElement(name.c_str());
-      element->InsertEndChild(child);
-      
-      return object->Serialize(doc,child);
-    }
-    
-    bool XMLSerializableObject::Deserialize(char& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return getElementAttribute(obj,doc,element,name);
-    }
-    
-    // template specialisation for char*, it interpreted as char array and can be used to handle strings
-    bool XMLSerializableObject::Deserialize(char*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return getElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(std::string& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return getElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(std::string*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(bool& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return getElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(bool*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(unsigned int& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return getElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(unsigned int*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(int& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return getElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(int*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(float& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return getElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(float*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(double& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return getElementAttribute(obj,doc,element,name);
-    }
-    
-    bool XMLSerializableObject::Deserialize(double*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Deserialize(T& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = T();
-      XMLSerializable* object = static_cast<XMLSerializable*>(&obj);
-      
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      
-      object->Name = child->FirstChild()->Value();
-      
-      if(child != NULL)
-      {
-        obj.Deserialize(doc,child);
-      }
-      else
-      {
-        obj.Init();
-        return false;
-      }
-      
-      return true;
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Deserialize(T*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = new T();
-      XMLSerializable* object = static_cast<XMLSerializable*>(obj);
-      
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      
-      object->Name = child->FirstChild()->Value();
-      
-      if(child != NULL)
-      {
-        obj->Deserialize(doc,child);
-      }
-      else
-      {
-        obj->Init();
-        return false;
-      }
-      
-      return true;
-    }
-    
-    /*
-     template<typename T, size_t S>
-     bool XMLSerializableObject::Serialize(std::array<T,S>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-     {
-     tinyxml2::XMLElement* ar = doc->NewElement(name.c_str());
-     element->InsertEndChild(ar);
-     
-     ar->SetAttribute("size",(unsigned int)obj.size());
-     
-     std::stringstream num;
-     for(unsigned int i=0;i<obj.size();i++)
-     {
-     num.str("");
-     num << "value" << i;
-     Serialize(obj[i],doc,ar,num.str());
-     }
-     
-     return true;
-     }
-     
-     template<typename T, size_t S>
-     bool XMLSerializableObject::Serialize(std::array<T,S>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-     {
-     return Serialize(*obj,doc,element,name);
-     }
-     
-     template<typename T, size_t S>
-     bool XMLSerializableObject::Deserialize(std::array<T,S>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-     {
-     bool res = true;
-     
-     const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-     if(child != NULL)
-     {
-     int size = child->UnsignedAttribute("size");
-     size = S < size ? S : size;
-     
-     std::stringstream num;
-     const tinyxml2::XMLAttribute* attribute = NULL;
-     for(unsigned int i=0;i<size;i++)
-     {
-     num.str("");
-     num << "value" << i;
-     
-     res &= Deserialize(obj[i],doc,child,num.str());
-     }
-     }
-     else
-     return false;
-     
-     return res;
-     }
-     
-     template<typename T, size_t S>
-     bool XMLSerializableObject::Deserialize(std::array<T,S>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-     {
-     obj = new std::array<T,S>();
-     return Deserialize(*obj,doc,element,name);
-     }
-     */
-    template<typename T0, typename T1>
-    bool XMLSerializableObject::Serialize(std::pair<T0,T1>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      bool res = true;
-      
-      tinyxml2::XMLElement* pair = doc->NewElement(name.c_str());
-      element->InsertEndChild(pair);
-      
-      res &= Serialize(obj.first,doc,pair,"first");
-      res &= Serialize(obj.second,doc,pair,"second");
-      
-      return res;
-    }
-    
-    template<typename T0, typename T1>
-    bool XMLSerializableObject::Serialize(std::pair<T0,T1>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    template<typename T0, typename T1>
-    bool XMLSerializableObject::Deserialize(std::pair<T0,T1>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      bool res = true;
-      
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      if(child != NULL)
-      {
-        res &= Deserialize(obj.first,doc,child,"first");
-        res &= Deserialize(obj.second,doc,child,"second");
-      }
-      else
-      {
-        Init(obj.first);
-        Init(obj.second);
-        return false;
-      }
-      
-      return res;
-    }
-    
-    template<typename T0, typename T1>
-    bool XMLSerializableObject::Deserialize(std::pair<T0,T1>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = new std::pair<T0,T1>();
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Serialize(std::vector<T>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      tinyxml2::XMLElement* vector = doc->NewElement(name.c_str());
-      element->InsertEndChild(vector);
-      
-      vector->SetAttribute("size",(unsigned int)obj.size());
-      
-      std::stringstream num;
-      for(unsigned int i=0;i<obj.size();i++)
-      {
-        num.str("");
-        num << "value" << i;
-        Serialize(obj[i],doc,vector,num.str());
-      }
-      
-      return true;
-    }
-
-    template<typename T>
-    bool XMLSerializableObject::Serialize(std::vector<T>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-
-    bool XMLSerializableObject::Serialize(std::vector<bool>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      tinyxml2::XMLElement* matrix = doc->NewElement(name.c_str());
-      element->InsertEndChild(matrix);
-      
-      const unsigned int size = obj.size();
-      
-      matrix->SetAttribute("size",size);
-      
-      std::stringstream ms;
-      ms << "\n";
-      for(unsigned int i=0;i<size;i++)
-        ms << obj[i] << ",";
-      ms << "\n";
-      
-      std::string mString = ms.str();
-      if(mString.size() > 1)
-        mString[mString.size()-2] = '\0';
-      
-      matrix->SetAttribute("vector_bool",mString.c_str());
-      
-      return true;
-    }
-
-    bool XMLSerializableObject::Serialize(std::vector<bool>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Deserialize(std::vector<T>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      bool res = true;
-      obj.clear();
-      
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      if(child != NULL)
-      {
-        unsigned int size = child->UnsignedAttribute("size");
-        obj.resize(size);
-        
-        std::stringstream num;
-        for(unsigned int i=0;i<size;i++)
-        {
-          num.str("");
-          num << "value" << i;
-          
-          res &= Deserialize(obj[i],doc,child,num.str());
-        }
-      }
-      else
-        return false;
-      
-      return res;
-    }
-
-    template<typename T>
-    bool XMLSerializableObject::Deserialize(std::vector<T>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = new std::vector<T>();
-      return Deserialize(*obj,doc,element,name);
-    }
-
-    bool XMLSerializableObject::Deserialize(std::vector<bool>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      if(child != NULL)
-      {
-        const unsigned int size = child->UnsignedAttribute("size");
-        
-        obj.resize(size);
-        
-        const tinyxml2::XMLAttribute* attribute = child->FindAttribute("vector_bool");
-        if(attribute == NULL)
-        {
-          Init(obj);
-          return false;
-        }
-        
-        char* matTemp;
-        GetAttribute(attribute->Value(),matTemp);
-        
-        std::string line, ssize;
-        std::stringstream mats;
-        mats.str(matTemp);
-        
-        int r=0;
-        std::string val;
-        // for each line
-        getline(mats,line); // matrix starts with an empty line
-        while(getline(mats,line))
-        {
-          // get current line
-          std::stringstream liness(line);
-          
-          for(unsigned int c=0;c<size-1;c++)
-          {
-            // split line
-            getline(liness, val, ',');
-            
-            // push pack the data if any
-            if(!val.empty())
-            {
-              bool t;
-              GetAttribute(val.c_str(),t);
-              obj[c] = t;
-            }
-          }
-          
-          getline(liness, val);
-          bool t;
-          GetAttribute(val.c_str(),t);
-          obj[size-1] = t;
-          
-          r++;
-        }
-      }
-      else
-      {
-        Init(obj);
-        return false;
-      }
-      
-      return true;
-    }
-
-    bool XMLSerializableObject::Deserialize(std::vector<bool>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = new std::vector<bool>();
-      return Deserialize(*obj,doc,element,name);
-    }
-
-    template<typename T>
-    bool XMLSerializableObject::Serialize(std::set<T>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      tinyxml2::XMLElement* vector = doc->NewElement(name.c_str());
-      element->InsertEndChild(vector);
-      
-      vector->SetAttribute("size",(unsigned int)obj.size());
-      
-      std::stringstream num;
-      typename std::set<T>::iterator iter = obj.begin();
-      for(int i=0;iter!=obj.end();iter++,i++)
-      {
-        num.str("");
-        num << "value" << i;
-        Serialize((T)*iter,doc,vector,num.str());
-      }
-      
-      return true;
-    }
-
-    template<typename T>
-    bool XMLSerializableObject::Serialize(std::set<T>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-
-    template<typename T>
-    bool XMLSerializableObject::Deserialize(std::set<T>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      bool res = true;
-      obj.clear();
-      
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      if(child != NULL)
-      {
-        unsigned int size = child->UnsignedAttribute("size");
-        
-        std::stringstream num;
-        typename std::set<T>::iterator iter = obj.begin();
-        for(int i=0;i<size;i++)
-        {
-          num.str("");
-          num << "value" << i;
-          
-          T val;
-          res &= Deserialize(val,doc,child,num.str());
-          obj.insert(val);
-        }
-      }
-      else
-        return false;
-      
-      return res;
-    }
-
-    template<typename T>
-    bool XMLSerializableObject::Deserialize(std::set<T>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = new std::set<T>();
-      return Deserialize(*obj,doc,element,name);
-    }
-
-    template<typename T0, typename T1>
-    bool XMLSerializableObject::Serialize(std::map<T0,T1>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      tinyxml2::XMLElement* vector = doc->NewElement(name.c_str());
-      element->InsertEndChild(vector);
-      
-      vector->SetAttribute("size",(unsigned int)obj.size());
-      
-      std::stringstream num;
-      typename std::map<T0,T1>::iterator iter = obj.begin();
-      for(int i=0;iter!=obj.end();iter++,i++)
-      {
-        num.str("");
-        num << "value" << i;
-        Serialize((std::pair<T0,T1>)*iter,doc,vector,num.str());
-      }
-      
-      return true;
-    }
-
-    template<typename T0, typename T1>
-    bool XMLSerializableObject::Serialize(std::map<T0,T1>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-
-    template<typename T0, typename T1>
-    bool XMLSerializableObject::Deserialize(std::map<T0,T1>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      bool res = true;
-      obj.clear();
-      
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      if(child != NULL)
-      {
-        unsigned int size = child->UnsignedAttribute("size");
-        
-        std::stringstream num;
-        typename std::map<T0,T1>::iterator iter = obj.begin();
-        for(int i=0;i<size;i++)
-        {
-          num.str("");
-          num << "value" << i;
-          
-          std::pair<T0,T1> pair;
-          res &= Deserialize(pair,doc,child,num.str());
-          obj.insert(pair);
-        }
-      }
-      else
-        return false;
-      
-      return res;
-    }
-
-    template<typename T0, typename T1>
-    bool XMLSerializableObject::Deserialize(std::map<T0,T1>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = new std::map<T0,T1>();
-      return Deserialize(*obj,doc,element,name);
-    }
-
-    template<typename T, int R, int C>
-    bool XMLSerializableObject::Serialize(Eigen::Matrix<T,R,C>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      tinyxml2::XMLElement* matrix = doc->NewElement(name.c_str());
-      element->InsertEndChild(matrix);
-      
-      const unsigned int rows = obj.rows();
-      const unsigned int cols = obj.cols();
-      
-      matrix->SetAttribute("rows",rows);
-      matrix->SetAttribute("cols",cols);
-      
-      std::stringstream ms;
-      ms << "\n";
-      for(unsigned int r=0;r<rows;r++)
-      {
-        for(unsigned int c=0;c<cols;c++)
-        {
-          ms << obj(r,c) << ",";
-        }
-        ms << "\n";
-      }
-      
-      std::string mString = ms.str();
-      if(mString.size() > 1)
-        mString[mString.size()-2] = '\0';
-      
-      matrix->SetAttribute("matrix",mString.c_str());
-      
-      return true;
-    }
-    
-    template<typename T, int R, int C>
-    bool XMLSerializableObject::Serialize(Eigen::Matrix<T,R,C>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    template<typename T, int R, int C>
-    bool XMLSerializableObject::Deserialize(Eigen::Matrix<T,R,C>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      if(child != NULL)
-      {
-        const unsigned int rows = child->UnsignedAttribute("rows");
-        const unsigned int cols = child->UnsignedAttribute("cols");
-        
-        obj.resize(rows,cols);
-        
-        const tinyxml2::XMLAttribute* attribute = child->FindAttribute("matrix");
-        if(attribute == NULL)
-        {
-          Init(obj);
-          return false;
-        }
-        
-        char* matTemp;
-        GetAttribute(attribute->Value(),matTemp);
-        
-        std::string line, srows, scols;
-        std::stringstream mats;
-        mats.str(matTemp);
-        
-        int r=0;
-        std::string val;
-        // for each line
-        getline(mats,line);
-        while(getline(mats,line))
-        {
-          // get current line
-          std::stringstream liness(line);
-          
-          for(unsigned int c=0;c<cols-1;c++)
-          {
-            // split line
-            getline(liness, val, ',');
-            
-            // push pack the data if any
-            if(!val.empty())
-              GetAttribute(val.c_str(),obj.coeffRef(r,c));
-          }
-          
-          getline(liness, val);
-          GetAttribute(val.c_str(),obj.coeffRef(r,cols-1));
-          
-          r++;
-        }
-      }
-      else
-      {
-        Init(obj);
-        return false;
-      }
-      
-      return true;
-    }
-    
-    template<typename T, int R, int C>
-    bool XMLSerializableObject::Deserialize(Eigen::Matrix<T,R,C>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = new Eigen::PlainObjectBase<T>();
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Serialize(Eigen::SparseMatrix<T>& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      tinyxml2::XMLElement* matrix = doc->NewElement(name.c_str());
-      element->InsertEndChild(matrix);
-      
-      const unsigned int rows = obj.rows();
-      const unsigned int cols = obj.cols();
-      
-      matrix->SetAttribute("rows",rows);
-      matrix->SetAttribute("cols",cols);
-      
-      std::stringstream ms;
-      ms << "\n";
-      for (int k=0;k<obj.outerSize();++k)
-      {
-        for (typename Eigen::SparseMatrix<T>::InnerIterator it(obj,k);it;++it)
-        {
-          ms << it.row() << "," << it.col() << "," << it.value() << "\n";
-        }
-      }
-      
-      std::string mString = ms.str();
-      if(mString.size() > 0)
-        mString[mString.size()-1] = '\0';
-      
-      matrix->SetAttribute("matrix",mString.c_str());
-      
-      return true;
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Serialize(Eigen::SparseMatrix<T>*& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return Serialize(*obj,doc,element,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Deserialize(Eigen::SparseMatrix<T>& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      if(child != NULL)
-      {
-        const unsigned int rows = child->UnsignedAttribute("rows");
-        const unsigned int cols = child->UnsignedAttribute("cols");
-        
-        obj.resize(rows,cols);
-        obj.setZero();
-        
-        const tinyxml2::XMLAttribute* attribute = child->FindAttribute("matrix");
-        if(attribute == NULL)
-        {
-          Init(obj);
-          return false;
-        }
-        
-        char* matTemp;
-        GetAttribute(attribute->Value(),matTemp);
-        
-        std::string line, srows, scols;
-        std::stringstream mats;
-        mats.str(matTemp);
-        
-        std::vector<Eigen::Triplet<T> > triplets;
-        int r=0;
-        std::string val;
-        
-        // for each line
-        getline(mats,line);
-        while(getline(mats,line))
-        {
-          // get current line
-          std::stringstream liness(line);
-          
-          // row
-          getline(liness, val, ',');
-          int row = atoi(val.c_str());
-          // col
-          getline(liness, val, ',');
-          int col = atoi(val.c_str());
-          // val
-          getline(liness, val);
-          T value;
-          GetAttribute(val.c_str(),value);
-          
-          triplets.push_back(Eigen::Triplet<T>(row,col,value));
-          
-          r++;
-        }
-        
-        obj.setFromTriplets(triplets.begin(),triplets.end());
-      }
-      else
-      {
-        Init(obj);
-        return false;
-      }
-      
-      return true;
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::Deserialize(Eigen::SparseMatrix<T>*& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      obj = new Eigen::SparseMatrix<T>();
-      return Deserialize(*obj,doc,element,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::setElementAttribute(T& obj, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      tinyxml2::XMLElement* child = doc->NewElement(name.c_str());
-      element->InsertEndChild(child);
-      SetAttribute(child,"val",obj);
-      return true;
-    }
-    
-    template<typename T>
-    bool XMLSerializableObject::getElementAttribute(T& obj, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element, const std::string& name)
-    {
-      // basic data type
-      const tinyxml2::XMLElement* child = element->FirstChildElement(name.c_str());
-      if(child != NULL)
-      {
-        XMLSerializableObject::GetAttribute(child->Attribute("val"),obj);
-        return true;
-      }
-      else
-      {
-        Init(obj);
-        return false;
-      }
-    }
-    
-    template<typename T>
-    XMLSerializableInstance<T>::XMLSerializableInstance(T& obj, const std::string& name, const std::string group)
-    : XMLSerializableObject(name, group), Object(obj)
-    {
-      XMLSerializableObject::Init(DefaultValue);
-    }
-    
-    template<typename T>
-    XMLSerializableInstance<T>::XMLSerializableInstance(T& obj, const std::string& name, const std::string group, T defaultValue)
-    : XMLSerializableObject(name, group), Object(obj), DefaultValue(defaultValue)
-    {
-    }
-    
-    template<typename T>
-    XMLSerializableInstance<T>::~XMLSerializableInstance()
-    {
-    }
-    
-    template<typename T>
-    void XMLSerializableInstance<T>::Init()
-    {
-      XMLSerializableObject::Init(DefaultValue);
-    }
-    
-    template<typename T>
-    bool XMLSerializableInstance<T>::Serialize(tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element)
-    {
-      return XMLSerializableObject::Serialize(Object,doc,element,Name);
-    }
-    
-    template<typename T>
-    bool XMLSerializableInstance<T>::Deserialize(tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element)
-    {
-      return XMLSerializableObject::Deserialize(Object,doc,element,Name);
-    }
-    
-    template<typename T>
-    bool XMLSerializer::SaveObject(T& object, const char* filename)
-    {
-      return SaveObject(object,"Object","Serialization",filename,true);
-    }
-    
-    template<typename T>
-    bool XMLSerializer::SaveObject(T& object, const std::string& objectName, const std::string& groupName, const char* filename, bool overwrite)
-    {
-      bool result = true;
-      XMLSerializer* serializer = new XMLSerializer(groupName);
-      result &= serializer->Add(object,objectName);
-      result &= serializer->Save(objectName,groupName,filename,overwrite);
-      delete serializer;
-      return result;
-    }
-    
-    template<typename T>
-    bool XMLSerializer::LoadObject(T& object, const char* filename)
-    {
-      return LoadObject(object,"Object","Serialization",filename);
-    }
-    
-    template<typename T>
-    bool XMLSerializer::LoadObject(T& object, const std::string& objectName, const std::string& groupName, const char* filename)
-    {
-      bool result = true;
-      XMLSerializer* serializer = new XMLSerializer(groupName);
-      result &= serializer->Add(object,objectName);
-      result &= serializer->Load(objectName,groupName,filename);
-      delete serializer;
-      return result;
-    }
-    
-    XMLSerializer::XMLSerializer(const std::string& defaultGroup)
-    {
-      SetCurrentGroup(defaultGroup);
-    }
-    
-    XMLSerializer::~XMLSerializer()
-    {
-      std::map<std::string,XMLSerializerGroup*>::iterator it;
-      for (it=groups.begin();it!=groups.end();it++)
-      {
-        delete it->second->Objects;
-        delete it->second;
-      }
-    }
-    
-    bool XMLSerializer::Save(const char* filename, bool overwrite)
-    {
-      tinyxml2::XMLDocument* doc = new tinyxml2::XMLDocument();
-      
-      if(overwrite == false)
-      {
-        // Check if file exists
-        tinyxml2::XMLError error = doc->LoadFile(filename);
-        if(error != tinyxml2::XML_NO_ERROR)
-          doc->Clear();
-      }
-      
-      if(SaveToXMLDoc(doc) == false)
-        return false;
-      
-      // Save
-      tinyxml2::XMLError error = doc->SaveFile(filename);
-      if(error != tinyxml2::XML_NO_ERROR)
-      {
-        doc->PrintError();
-        return false;
-      }
-      
-      delete doc;
-      
-      return true;
-    }
-    
-    bool XMLSerializer::SaveToXMLDoc(tinyxml2::XMLDocument* doc)
-    {
-      std::map<std::string,XMLSerializerGroup*>::iterator it;
-      for (it=groups.begin();it!=groups.end();it++)
-      {
-        // Update group
-        tinyxml2::XMLElement* element = doc->FirstChildElement(it->first.c_str());
-        if(element != NULL)
-        {
-          element->DeleteChildren();
-        }
-        else
-        {
-          element = doc->NewElement(it->first.c_str());
-          doc->InsertEndChild(element);
-        }
-        
-        std::vector<XMLSerializable*>* group = it->second->Objects;
-        for(unsigned  int i=0;i<group->size();i++)
-        {
-          if((*group)[i]->Serialize(doc,element) == false)
-            return false;
-        }
-      }
-      
-      return true;
-    }
-    
-    bool XMLSerializer::Save(const std::string& groupName, const char* filename, bool overwrite)
-    {
-      tinyxml2::XMLDocument* doc = new tinyxml2::XMLDocument();
-      
-      if(overwrite == false)
-      {
-        // Check if file exists
-        tinyxml2::XMLError error = doc->LoadFile(filename);
-        if(error != tinyxml2::XML_NO_ERROR)
-          doc->Clear();
-      }
-      
-      if(SaveToXMLDoc(groupName, doc) == false)
-        return false;
-      
-      // Save
-      tinyxml2::XMLError error = doc->SaveFile(filename);
-      if(error != tinyxml2::XML_NO_ERROR)
-      {
-        doc->PrintError();
-        return false;
-      }
-      
-      delete doc;
-      
-      return true;
-    }
-    
-    bool XMLSerializer::SaveToXMLDoc(const std::string& groupName, tinyxml2::XMLDocument* doc)
-    {
-      std::string gn = groupName;
-      EncodeXMLElementName(gn);
-      
-      std::map<std::string,XMLSerializerGroup*>::iterator it = groups.find(gn);
-      if(it == groups.end())
-        return false;
-      
-      // Update group
-      tinyxml2::XMLElement* element = doc->FirstChildElement(it->first.c_str());
-      if(element != NULL)
-      {
-        element->DeleteChildren();
-      }
-      else
-      {
-        element = doc->NewElement(it->first.c_str());
-        doc->InsertEndChild(element);
-      }
-      
-      std::vector<XMLSerializable*>* groups = it->second->Objects;
-      for(unsigned int i=0;i<groups->size();i++)
-      {
-        if((*groups)[i]->Serialize(doc,element) == false)
-          return false;
-      }
-      
-      return true;
-    }
-    
-    bool XMLSerializer::Save(const std::string& objectName, const std::string& groupName, const char* filename, bool overwrite)
-    {
-      tinyxml2::XMLDocument* doc = new tinyxml2::XMLDocument();
-      
-      if(overwrite == false)
-      {
-        // Check if file exists
-        tinyxml2::XMLError error = doc->LoadFile(filename);
-        if(error != tinyxml2::XML_NO_ERROR)
-          doc->Clear();
-      }
-      
-      if(SaveToXMLDoc(objectName, groupName, doc) == false)
-        return false;
-      
-      // Save
-      tinyxml2::XMLError error = doc->SaveFile(filename);
-      if(error != tinyxml2::XML_NO_ERROR)
-      {
-        doc->PrintError();
-        return false;
-      }
-      
-      delete doc;
-      
-      return true;
-    }
-    
-    bool XMLSerializer::SaveToXMLDoc(const std::string& objectName, const std::string& groupName, tinyxml2::XMLDocument* doc)
-    {
-      std::string gn = groupName;
-      EncodeXMLElementName(gn);
-      
-      std::string on = objectName;
-      EncodeXMLElementName(on);
-      
-      std::map<std::string,XMLSerializerGroup*>::iterator it = groups.find(gn);
-      if(it == groups.end())
-        return false;
-      
-      // Get/Add group
-      tinyxml2::XMLElement* element = findAddGroup(doc, it->first.c_str());
-      
-      // Serialize
-      std::vector<XMLSerializable*>* groups = it->second->Objects;
-      bool found = false;
-      for(unsigned int i=0;i<groups->size();i++)
-      {
-        if((*groups)[i]->Name == on)
-        {
-          found = true;
-          
-          tinyxml2::XMLElement* child = element->FirstChildElement(on.c_str());
-          if(child != NULL)
-          {
-            element->DeleteChild(child);
-          }
-          
-          if((*groups)[i]->Serialize(doc,element) == false)
-            return false;
-        }
-      }
-      
-      return found;
-    }
-    
-    bool XMLSerializer::SaveGroupToXMLElement(tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      return SaveGroupToXMLElement(currentGroup->first,doc,element,name);
-    }
-    
-    bool XMLSerializer::SaveGroupToXMLElement(const std::string& groupName, tinyxml2::XMLDocument* doc, tinyxml2::XMLElement* element, const std::string& name)
-    {
-      std::string gn = groupName;
-      EncodeXMLElementName(gn);
-      
-      std::map<std::string,XMLSerializerGroup*>::iterator it = groups.find(gn);
-      if(it == groups.end())
-        return false;
-      
-      // Add new group
-      tinyxml2::XMLElement* group = doc->NewElement(name.c_str());
-      element->InsertEndChild(group);
-      
-      std::vector<XMLSerializable*>* groups = it->second->Objects;
-      for(unsigned int i=0;i<groups->size();i++)
-      {
-        if((*groups)[i]->Serialize(doc,group) == false)
-          return false;
-      }
-      
-      return true;
-    }
-    
-    bool XMLSerializer::Load(const char* filename)
-    {
-      tinyxml2::XMLDocument* doc = openDoc(filename);
-      if(doc == NULL)
-        return false;
-      
-      if(LoadFromXMLDoc(doc) == false)
-        return false;
-      
-      delete doc;
-      
-      return true;
-    }
-    
-    bool XMLSerializer::LoadFromXMLDoc(tinyxml2::XMLDocument* doc)
-    {
-      std::map<std::string,XMLSerializerGroup*>::iterator it;
-      for (it=groups.begin();it!=groups.end();it++)
-      {
-        tinyxml2::XMLElement* element = doc->FirstChildElement(it->first.c_str());
-        if(element == NULL)
-          return false;
-        
-        // Deserialize
-        std::vector<XMLSerializable*>* group = it->second->Objects;
-        for(unsigned int i=0;i<group->size();i++)
-        {
-          if(element == NULL || (*group)[i]->Deserialize(doc,element) == false)
-            (*group)[i]->Init(); // Load default value;
-        }
-      }
-      
-      return true;
-    }
-    
-    bool XMLSerializer::Load(const std::string& groupName, const char* filename)
-    {
-      tinyxml2::XMLDocument* doc = openDoc(filename);
-      if(doc == NULL)
-        return false;
-      
-      if(LoadFromXMLDoc(groupName, doc) == false)
-        return false;
-      
-      delete doc;
-      
-      return true;
-    }
-    
-    bool XMLSerializer::LoadFromXMLDoc(const std::string& groupName, tinyxml2::XMLDocument* doc)
-    {
-      std::string gn = groupName;
-      EncodeXMLElementName(gn);
-      
-      std::map<std::string,XMLSerializerGroup*>::iterator it = groups.find(gn);
-      if(it == groups.end())
-        return false;
-      
-      tinyxml2::XMLElement* element = doc->FirstChildElement(it->first.c_str());
-      if(element == NULL)
-        return false;
-      
-      // Deserialize
-      std::vector<XMLSerializable*>* groups = it->second->Objects;
-      for(unsigned int i=0;i<groups->size();i++)
-      {
-        if(element == NULL || (*groups)[i]->Deserialize(doc,element) == false)
-          (*groups)[i]->Init(); // Load default value;
-      }
-      
-      return true;
-    }
-    
-    bool XMLSerializer::Load(const std::string& objectName, const std::string& groupName, const char* filename)
-    {
-      tinyxml2::XMLDocument* doc = openDoc(filename);
-      if(doc == NULL)
-        return false;
-      
-      if(LoadFromXMLDoc(objectName,groupName,doc) == false)
-        return false;
-      
-      delete doc;
-      
-      return true;
-    }
-    
-    bool XMLSerializer::LoadFromXMLDoc(const std::string& objectName, const std::string& groupName, tinyxml2::XMLDocument* doc)
-    {
-      std::string gn = groupName;
-      EncodeXMLElementName(gn);
-      
-      std::string on = objectName;
-      EncodeXMLElementName(on);
-      
-      std::map<std::string,XMLSerializerGroup*>::iterator it = groups.find(gn);
-      if(it == groups.end())
-        return false;
-      
-      tinyxml2::XMLElement* element = doc->FirstChildElement(it->first.c_str());
-      if(element == NULL)
-        return false;
-      
-      // Deserialize
-      std::vector<XMLSerializable*>* groups = it->second->Objects;
-      bool found = false;
-      for(unsigned int i=0;i<groups->size();i++)
-      {
-        if((*groups)[i]->Name == on)
-        {
-          found = true;
-          if(element == NULL || (*groups)[i]->Deserialize(doc,element) == false)
-            (*groups)[i]->Init(); // Load default value;
-        }
-      }
-      
-      return found;
-    }
-    
-    bool XMLSerializer::LoadGroupFromXMLElement(const std::string& groupName, tinyxml2::XMLDocument* doc, const tinyxml2::XMLElement* element)
-    {
-      std::string gn = groupName;
-      EncodeXMLElementName(gn);
-      
-      std::map<std::string,XMLSerializerGroup*>::iterator it = groups.find(gn);
-      
-      if(it == groups.end())
-        return false;
-      
-      const tinyxml2::XMLElement* group = element->FirstChildElement(groupName.c_str());
-      if(group == NULL)
-        return false;
-      
-      // Deserialize
-      std::vector<XMLSerializable*>* groups = it->second->Objects;
-      for(unsigned int i=0;i<groups->size();i++)
-      {
-        if(element == NULL || (*groups)[i]->Deserialize(doc,group) == false)
-          (*groups)[i]->Init(); // Load default value;
-      }
-      
-      return true;
-    }
-    
-    void XMLSerializer::SetCurrentGroup(const std::string& group)
-    {
-      currentGroup = setGetGroup(group);
-    }
-    
-    std::string XMLSerializer::GetCurrentGroup()
-    {
-      return currentGroup->first;
-    }
-    
-    template<typename T>
-    bool XMLSerializer::Add(T& obj, const std::string& name)
-    {
-      XMLSerializable* object = static_cast<XMLSerializable*>(&obj);
-      
-      object->Name = name;
-      return addObjectToGroup(object,currentGroup);
-    }
-    
-    bool XMLSerializer::Add(char& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    bool XMLSerializer::Add(char*& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    bool XMLSerializer::Add(std::string& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    bool XMLSerializer::Add(bool& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    bool XMLSerializer::Add(unsigned int& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    bool XMLSerializer::Add(int& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    bool XMLSerializer::Add(float& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    bool XMLSerializer::Add(double& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    /*template<typename T, size_t S>
-     bool XMLSerializer::Add(std::array<T,S>& obj, const std::string& name)
-     {
-     return add(obj,name);
-     }*/
-    
-    template<typename T0, typename T1>
-    bool XMLSerializer::Add(std::pair<T0,T1>& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializer::Add(std::vector<T>& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-
-    template<typename T>
-    bool XMLSerializer::Add(std::set<T>& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-
-    template<typename T0, typename T1>
-    bool XMLSerializer::Add(std::map<T0,T1>& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    template<typename T, int R, int C>
-    bool XMLSerializer::Add(Eigen::Matrix<T,R,C>& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializer::Add(Eigen::SparseMatrix<T>& obj, const std::string& name)
-    {
-      return add(obj,name);
-    }
-    
-    template<typename T>
-    bool XMLSerializer::Add(T& object, const std::string& name, T defaultValue)
-    {
-      return false;
-    }
-    
-    bool XMLSerializer::Add(char& obj, const std::string& name, char defaultValue)
-    {
-      return add(obj,name,defaultValue);
-    }
-    
-    bool XMLSerializer::Add(char*& obj, const std::string& name, char* defaultValue)
-    {
-      return add(obj,name,defaultValue);
-    }
-    
-    bool XMLSerializer::Add(std::string& obj, const std::string& name, std::string defaultValue)
-    {
-      return add(obj,name,defaultValue);
-    }
-    
-    bool XMLSerializer::Add(bool& obj, const std::string& name, bool defaultValue)
-    {
-      return add(obj,name,defaultValue);
-    }
-    
-    bool XMLSerializer::Add(unsigned int& obj, const std::string& name, unsigned int defaultValue)
-    {
-      return add(obj,name,defaultValue);
-    }
-    
-    bool XMLSerializer::Add(int& obj, const std::string& name, int defaultValue)
-    {
-      return add(obj,name,defaultValue);
-    }
-    
-    bool XMLSerializer::Add(float& obj, const std::string& name, float defaultValue)
-    {
-      return add(obj,name,defaultValue);
-    }
-    
-    bool XMLSerializer::Add(double& obj, const std::string& name, double defaultValue)
-    {
-      return add(obj,name,defaultValue);
-    }
-    
-    template<typename T>
-    bool XMLSerializer::add(T& obj, const std::string& name)
-    {
-      XMLSerializable* object = new XMLSerializableInstance<T>(obj,name,currentGroup->first);
-      return addObjectToGroup(object,currentGroup);
-    }
-    
-    template<typename T>
-    bool XMLSerializer::add(T& obj, const std::string& name, T defaultValue)
-    {
-      XMLSerializable* object = new XMLSerializableInstance<T>(obj,name,currentGroup->first,defaultValue);
-      return addObjectToGroup(object,currentGroup);
-    }
-    
-    bool XMLSerializer::addObjectToGroup(XMLSerializable* obj, const std::string& group)
-    {
-      std::map<std::string,XMLSerializerGroup*>::iterator it = setGetGroup(group);
-      return addObjectToGroup(obj, it);
-    }
-    
-    bool XMLSerializer::addObjectToGroup(XMLSerializable* object, std::map<std::string,XMLSerializerGroup*>::iterator it)
-    {
-      std::vector<XMLSerializable*>* objects = it->second->Objects;
-      for(unsigned int i=0;i<objects->size();i++)
-      {
-        if((*objects)[i]->Name == object->Name)
-          return false;
-      }
-      
-      objects->push_back(object);
-      
-      return true;
-    }
-    
-    std::map<std::string,XMLSerializerGroup*>::iterator XMLSerializer::setGetGroup(const std::string& group)
-    {
-      std::string groupName = group;
-      EncodeXMLElementName(groupName);
-      
-      std::map<std::string,XMLSerializerGroup*>::iterator it = groups.find(groupName);
-      if(it == groups.end())
-      {
-        XMLSerializerGroup* newGroup = new XMLSerializerGroup();
-        newGroup->Objects = new std::vector<XMLSerializable*>();
-        groups[groupName] = newGroup;
-        it = groups.find(groupName);
-      }
-      
-      return it;
-    }
-    
-    tinyxml2::XMLDocument* XMLSerializer::openDoc(const char* filename)
-    {
-      tinyxml2::XMLDocument* doc = new tinyxml2::XMLDocument();
-      
-      tinyxml2::XMLError error = doc->LoadFile(filename);
-      if(error != tinyxml2::XML_NO_ERROR)
-      {
-        doc->PrintError();
-        doc = NULL;
-      }
-      
-      return doc;
-    }
-    
-    tinyxml2::XMLElement* XMLSerializer::findAddGroup(tinyxml2::XMLDocument* doc, const char* groupName)
-    {
-      tinyxml2::XMLElement* group = doc->FirstChildElement(groupName);
-      if(group == NULL)
-      {
-        group = doc->NewElement(groupName);
-        doc->InsertEndChild(group);
-      }
-      return group;
-    }
-  }
-}
-#endif

include/igl/xml/XMLSerializer.h.REMOVED.git-id

@@ -0,0 +1 @@
+ef40e013b092c2c55d5032d93ad436ea05230a88

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.