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

Former-commit-id: 6da80efe0667ddaf038ca2a1142963fe06495e90

README.md

@@ -25,7 +25,7 @@ small "Hello, World" program:
 ```cpp
 #include <igl/cotmatrix.h>
 #include <Eigen/Dense>
-#include <Eigen/SparseMatrix>
+#include <Eigen/Sparse>
 #include <iostream>
 int main()
 {

RELEASE_HISTORY.txt

@@ -1,3 +1,5 @@
+1.0.2  Bug fix in winding number code
+1.0.1  Bug fixes and more CGAL support
 1.0.0  Major beta release: many renames, tutorial, triangle wrapper, org. build
 0.4.6  Generalized Winding Numbers
 0.4.5  CGAL extra: mesh selfintersection

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
-1.0.0
+1.0.2

build/Makefile.conf

@@ -106,7 +106,8 @@ endif
 ifeq ($(IGL_USERNAME),daniele)
 	IGL_WITH_MATLAB=0
 	AFLAGS=-m64
-	GG=g++-mp-4.7
+	GG=g++-4.7
+	EIGEN3_INC=-I/usr/local/include/eigen3
 endif
 
 ifeq ($(IGL_USERNAME),olkido)

examples/ambient-occlusion-mex/ambient_occlusion.m

@@ -12,3 +12,11 @@
   %    S  #P list of ambient occlusion values between 1 (fully occluded) and 0
   %      (not occluded)
   %
+  % Examples:
+  %   % mesh (V,F), scalar field Z
+  %   AO = ambient_occlusion(V,F,V,per_vertex_normals(V,F),1000);
+  %   tsurf(F,V,'FaceVertexCData', ...
+  %     bsxfun(@times,1-AO, ...
+  %       squeeze(ind2rgb(floor(matrixnormalize(Z)*256),jet(256)))), ...
+  %     fphong,'EdgeColor','none');
+  %

examples/arap/example.cpp

@@ -214,8 +214,8 @@ bool init_arap()
   {
     return false;
   }
-  //arap_data.with_dynamics = true;
-  //arap_data.h = 0.5;
+  arap_data.with_dynamics = true;
+  arap_data.h = 0.05;
   //arap_data.max_iter = 100;
   //partition(W,100,arap_data.G,_S,_D);
   return arap_precomputation(V,*Ele,V.cols(),b,arap_data);

examples/intersections/example.cpp

@@ -7,6 +7,7 @@
 #include <igl/quat_to_mat.h>
 #include <igl/report_gl_error.h>
 #include <igl/readOBJ.h>
+#include <igl/writeOBJ.h>
 #include <igl/readDMAT.h>
 #include <igl/readOFF.h>
 #include <igl/readMESH.h>
@@ -23,6 +24,7 @@
 #include <igl/ReAntTweakBar.h>
 #include <igl/pathinfo.h>
 #include <igl/Camera.h>
+#include <igl/cat.h>
 #include <igl/get_seconds.h>
 #include <igl/cgal/remesh_self_intersections.h>
 #include <igl/cgal/intersect_other.h>
@@ -123,7 +125,7 @@ float light_pos[4] = {0.1,0.1,-0.9,0};
 // C,D  Colors
 // N,W  Normals
 // mid  combined "centroid"
-Eigen::MatrixXd V,N,C,Z,mid,U,W,D;
+Eigen::MatrixXd V,N,C,Z,mid,U,W,D,VU;
 // F,G  faces
 Eigen::MatrixXi F,G;
 bool has_other = false;
@@ -272,6 +274,7 @@ void display()
     const MatrixXd & N,
     const MatrixXd & C)
   {
+    glEnable(GL_COLOR_MATERIAL);
     glEnable(GL_POLYGON_OFFSET_FILL); // Avoid Stitching!
     glPolygonOffset(1.0,1);
     glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
@@ -305,7 +308,7 @@ void display()
   // Draw a nice floor
   glPushMatrix();
   const double floor_offset =
-    -2./bbd*(V.col(1).maxCoeff()-mid(1));
+    -2./bbd*(VU.col(1).maxCoeff()-mid(1));
   glTranslated(0,floor_offset,0);
   const float GREY[4] = {0.5,0.5,0.6,1.0};
   const float DARK_GREY[4] = {0.2,0.2,0.3,1.0};
@@ -475,12 +478,14 @@ void color_selfintersections(
 {
   using namespace igl;
   using namespace Eigen;
+  using namespace std;
   MatrixXd SV;
   MatrixXi SF,IF;
   VectorXi J,IM;
   RemeshSelfIntersectionsParam params;
-  params.detect_only = true;
+  params.detect_only = false;
   remesh_self_intersections(V,F,params,SV,SF,IF,J,IM);
+  writeOBJ("FUCK.obj",SV,SF);
   C.resize(F.rows(),3);
   C.col(0).setConstant(0.4);
   C.col(1).setConstant(0.8);
@@ -633,18 +638,15 @@ int main(int argc, char * argv[])
     {
       return 1;
     }
-    mid = 0.25*(V.colwise().maxCoeff() + V.colwise().minCoeff()) +
-      0.25*(U.colwise().maxCoeff() + U.colwise().minCoeff());
-    bbd = max(
-      (V.colwise().maxCoeff() - V.colwise().minCoeff()).maxCoeff(),
-      (U.colwise().maxCoeff() - U.colwise().minCoeff()).maxCoeff());
+    cat(1,V,U,VU);
     color_intersections(V,F,U,G,C,D);
   }else
   {
-    mid = 0.5*(V.colwise().maxCoeff() + V.colwise().minCoeff());
-    bbd = (V.colwise().maxCoeff() - V.colwise().minCoeff()).maxCoeff();
+    VU = V;
     color_selfintersections(V,F,C);
   }
+  mid = 0.5*(VU.colwise().maxCoeff() + VU.colwise().minCoeff());
+  bbd = (VU.colwise().maxCoeff() - VU.colwise().minCoeff()).maxCoeff();
 
   // Init glut
   glutInit(&argc,argv);

examples/skeleton-builder/Makefile

@@ -0,0 +1,76 @@
+
+.PHONY: all
+
+# Shared flags etc.
+include ../../build/Makefile.conf
+
+all: obj example
+
+.PHONY: example
+
+LIBIGL=../../
+LIBIGL_INC=-I$(LIBIGL)/include
+LIBIGL_LIB=-L$(LIBIGL)/lib -ligl -liglembree
+
+EIGEN3_INC=-I/opt/local/include/eigen3 -I/opt/local/include/eigen3/unsupported
+
+ANTTWEAKBAR_INC=-I$(LIBIGL)/external/AntTweakBar/include
+ANTTWEAKBAR_LIB=-L$(LIBIGL)/external/AntTweakBar/lib -lAntTweakBar -framework AppKit
+
+TETGEN=$(LIBIGL)/external/tetgen
+TETGEN_LIB=-L$(TETGEN) -ltet 
+TETGEN_INC=-I$(TETGEN)
+
+EMBREE=$(LIBIGL)/external/embree
+EMBREE_INC=-I$(EMBREE)/ -I$(EMBREE)/embree
+EMBREE_LIB=-L$(EMBREE)/build -lembree -lsys
+
+CARBON_LIB=-framework Carbon
+
+# Use free glut for mouse scrolling
+#FREE_GLUT=/opt/local/
+#FREE_GLUT_INC=-I$(FREE_GLUT)/include
+#FREE_GLUT_LIB=-L$(FREE_GLUT)/lib -lglut
+GLUT_LIB=-framework GLUT
+GLUT_INC=
+
+ifdef IGL_NO_MOSEK
+CFLAGS+=-DIGL_NO_MOSEK
+else
+# Adjust your mosek paths etc. accordingly
+ifndef MOSEKPLATFORM
+  MOSEKPLATFORM=osx64x86
+endif
+ifndef MOSEKVERSION
+  MOSEKVERSION=6
+endif
+IGLMOSEK=
+IGLMOSEK_INC=
+MOSEK=/usr/local/mosek
+MOSEK_INC=
+MOSEK_LIB=
+endif
+
+INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(MATLAB_INC) $(GLUT_INC) ${TETGEN_INC} $(MOSEK_INC) $(EMBREE_INC)
+LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(MATLAB_LIB) $(CARBON_LIB) $(TETGEN_LIB) $(MOSEK_LIB) $(EMBREE_LIB)
+
+CPP_FILES=$(wildcard ./*.cpp)
+OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o))) 
+
+CFLAGS+=-std=c++11 -g
+
+example: obj $(OBJ_FILES)
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o example $(LIB) $(OBJ_FILES) 
+
+obj:
+	mkdir -p obj
+
+obj/%.o: %.cpp
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
+
+obj/%.o: %.cpp %.h
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
+
+clean:
+	rm -f $(OBJ_FILES)
+	rm -f example

examples/skeleton-builder/example.cpp

@@ -0,0 +1,1161 @@
+#include <igl/draw_skeleton_3d.h>
+#include <igl/draw_skeleton_vector_graphics.h>
+#include <igl/two_axis_valuator_fixed_up.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/readTGF.h>
+#include <igl/writeOBJ.h>
+#include <igl/writeOFF.h>
+#include <igl/report_gl_error.h>
+#include <igl/draw_mesh.h>
+#include <igl/draw_floor.h>
+#include <igl/pathinfo.h>
+#include <igl/list_to_matrix.h>
+#include <igl/quat_to_mat.h>
+#include <igl/per_face_normals.h>
+#include <igl/material_colors.h>
+#include <igl/trackball.h>
+#include <igl/snap_to_canonical_view_quat.h>
+#include <igl/snap_to_fixed_up.h>
+#include <igl/REDRUM.h>
+#include <igl/Camera.h>
+#include <igl/ReAntTweakBar.h>
+#include <igl/get_seconds.h>
+#include <igl/forward_kinematics.h>
+#include <igl/boundary_conditions.h>
+#include <igl/normalize_row_sums.h>
+#include <igl/lbs_matrix.h>
+#include <igl/sort_triangles.h>
+#include <igl/slice.h>
+#include <igl/project.h>
+#include <igl/unproject.h>
+#include <igl/embree/EmbreeIntersector.h>
+#include <igl/embree/unproject_in_mesh.h>
+#include <igl/matlab_format.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/adjacency_list.h>
+#include <igl/adjacency_matrix.h>
+#include <igl/right_axis.h>
+#include <igl/colon.h>
+#include <igl/unique.h>
+#include <igl/REDRUM.h>
+#include <igl/writeTGF.h>
+#include <igl/file_exists.h>
+#include <igl/centroid.h>
+
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+#ifdef __APPLE__
+#include <GLUT/glut.h>
+#include <Carbon/Carbon.h>
+#else
+#include <GL/glut.h>
+#endif
+
+#include <string>
+#include <vector>
+#include <queue>
+#include <stack>
+#include <iostream>
+
+enum SkelStyleType
+{
+  SKEL_STYLE_TYPE_3D = 0,
+  SKEL_STYLE_TYPE_VECTOR_GRAPHICS = 1,
+  NUM_SKEL_STYLE_TYPE = 2
+}skel_style;
+
+Eigen::MatrixXd V,N,sorted_N;
+Eigen::Vector3d Vmid,Vcen;
+double bbd = 1.0;
+Eigen::MatrixXi F,sorted_F;
+Eigen::VectorXi P;
+igl::Camera camera;
+struct State
+{
+  Eigen::MatrixXd C;
+  Eigen::MatrixXi BE;
+  Eigen::VectorXi sel;
+} s;
+
+bool wireframe = false;
+bool skeleton_on_top = false;
+double alpha = 0.5;
+
+// See README for descriptions
+enum RotationType
+{
+  ROTATION_TYPE_IGL_TRACKBALL = 0,
+  ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP = 1,
+  NUM_ROTATION_TYPES = 2,
+} rotation_type = ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP;
+
+std::stack<State> undo_stack;
+std::stack<State> redo_stack;
+
+bool is_rotating = false;
+bool is_dragging = false;
+bool new_leaf_on_drag = false;
+bool new_root_on_drag = false;
+int down_x,down_y;
+Eigen::MatrixXd down_C;
+igl::Camera down_camera;
+std::string output_filename;
+
+bool is_animating = false;
+double animation_start_time = 0;
+double ANIMATION_DURATION = 0.5;
+Eigen::Quaterniond animation_from_quat;
+Eigen::Quaterniond animation_to_quat;
+
+int width,height;
+Eigen::Vector4f light_pos(-0.1,-0.1,0.9,0);
+
+#define REBAR_NAME "temp.rbr"
+igl::ReTwBar rebar;
+igl::EmbreeIntersector ei;
+
+void push_undo()
+{
+  undo_stack.push(s);
+  // Clear
+  redo_stack = std::stack<State>();
+}
+
+// No-op setter, does nothing
+void TW_CALL no_op(const void * /*value*/, void * /*clientData*/)
+{
+}
+
+void TW_CALL set_rotation_type(const void * value, void * clientData)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  const RotationType old_rotation_type = rotation_type;
+  rotation_type = *(const RotationType *)(value);
+  if(rotation_type == ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP &&
+    old_rotation_type != ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP)
+  {
+    push_undo();
+    animation_from_quat = camera.m_rotation_conj;
+    snap_to_fixed_up(animation_from_quat,animation_to_quat);
+    // start animation
+    animation_start_time = get_seconds();
+    is_animating = true;
+  }
+}
+void TW_CALL get_rotation_type(void * value, void *clientData)
+{
+  RotationType * rt = (RotationType *)(value);
+  *rt = rotation_type;
+}
+
+void reshape(int width, int height)
+{
+  ::width = width;
+  ::height = height;
+  glViewport(0,0,width,height);
+  // Send the new window size to AntTweakBar
+  TwWindowSize(width, height);
+  camera.m_aspect = (double)width/(double)height;
+}
+
+void push_scene()
+{
+  using namespace igl;
+  using namespace std;
+  glMatrixMode(GL_PROJECTION);
+  glPushMatrix();
+  glLoadIdentity();
+  gluPerspective(camera.m_angle,camera.m_aspect,camera.m_near,camera.m_far);
+  glMatrixMode(GL_MODELVIEW);
+  glPushMatrix();
+  glLoadIdentity();
+  gluLookAt(
+    camera.eye()(0), camera.eye()(1), camera.eye()(2),
+    camera.at()(0), camera.at()(1), camera.at()(2),
+    camera.up()(0), camera.up()(1), camera.up()(2));
+}
+
+void push_object()
+{
+  using namespace igl;
+  glPushMatrix();
+  glScaled(2./bbd,2./bbd,2./bbd);
+  glTranslated(-Vmid(0),-Vmid(1),-Vmid(2));
+}
+
+void pop_object()
+{
+  glPopMatrix();
+}
+
+void pop_scene()
+{
+  glMatrixMode(GL_PROJECTION);
+  glPopMatrix();
+  glMatrixMode(GL_MODELVIEW);
+  glPopMatrix();
+}
+
+// Set up double-sided lights
+void lights()
+{
+  using namespace std;
+  using namespace Eigen;
+  glEnable(GL_LIGHTING);
+  glLightModelf(GL_LIGHT_MODEL_TWO_SIDE,GL_TRUE);
+  glEnable(GL_LIGHT0);
+  glEnable(GL_LIGHT1);
+  float WHITE[4] =  {0.8,0.8,0.8,1.};
+  float GREY[4] =  {0.4,0.4,0.4,1.};
+  float BLACK[4] =  {0.,0.,0.,1.};
+  Vector4f pos = light_pos;
+  glLightfv(GL_LIGHT0,GL_AMBIENT,GREY);
+  glLightfv(GL_LIGHT0,GL_DIFFUSE,WHITE);
+  glLightfv(GL_LIGHT0,GL_SPECULAR,BLACK);
+  glLightfv(GL_LIGHT0,GL_POSITION,pos.data());
+  pos(0) *= -1;
+  pos(1) *= -1;
+  pos(2) *= -1;
+  glLightfv(GL_LIGHT1,GL_AMBIENT,GREY);
+  glLightfv(GL_LIGHT1,GL_DIFFUSE,WHITE);
+  glLightfv(GL_LIGHT1,GL_SPECULAR,BLACK);
+  glLightfv(GL_LIGHT1,GL_POSITION,pos.data());
+}
+
+void sort()
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  push_scene();
+  push_object();
+  VectorXi I;
+  sort_triangles(V,F,sorted_F,I);
+  slice(N,I,1,sorted_N);
+  pop_object();
+  pop_scene();
+}
+
+void display()
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  const float back[4] = {0.75, 0.75, 0.75,0};
+  glClearColor(back[0],back[1],back[2],0);
+  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+  static bool first = true;
+  if(first)
+  {
+    sort();
+    first = false;
+  }
+
+  if(is_animating)
+  {
+    double t = (get_seconds() - animation_start_time)/ANIMATION_DURATION;
+    if(t > 1)
+    {
+      t = 1;
+      is_animating = false;
+    }
+    Quaterniond q = animation_from_quat.slerp(t,animation_to_quat).normalized();
+    camera.orbit(q.conjugate());
+  }
+
+  glEnable(GL_DEPTH_TEST);
+  glDepthFunc(GL_LEQUAL);
+  glEnable(GL_NORMALIZE);
+  glEnable(GL_BLEND);
+  glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
+  lights();
+  push_scene();
+
+  // Draw a nice floor
+  glEnable(GL_DEPTH_TEST);
+  glPushMatrix();
+  const double floor_offset =
+    -2./bbd*(V.col(1).maxCoeff()-Vmid(1));
+  glTranslated(0,floor_offset,0);
+  const float GREY[4] = {0.5,0.5,0.6,1.0};
+  const float DARK_GREY[4] = {0.2,0.2,0.3,1.0};
+  glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+  glEnable(GL_CULL_FACE);
+  glCullFace(GL_BACK);
+  draw_floor(GREY,DARK_GREY);
+  glDisable(GL_CULL_FACE);
+  glPopMatrix();
+
+  push_object();
+
+  const auto & draw_skeleton = []()
+  {
+    switch(skel_style)
+    {
+      default:
+      case SKEL_STYLE_TYPE_3D:
+      {
+        MatrixXf colors = MAYA_VIOLET.transpose().replicate(s.BE.rows(),1);
+        for(int si=0;si<s.sel.size();si++)
+        {
+          for(int b=0;b<s.BE.rows();b++)
+          {
+            if(s.BE(b,0) == s.sel(si) || s.BE(b,1) == s.sel(si))
+            {
+              colors.row(b) = MAYA_SEA_GREEN;
+            }
+          }
+        }
+        draw_skeleton_3d(s.C,s.BE,MatrixXd(),colors);
+        break;
+      }
+      case SKEL_STYLE_TYPE_VECTOR_GRAPHICS:
+        draw_skeleton_vector_graphics(s.C,s.BE);
+        break;
+    }
+  };
+  
+  if(!skeleton_on_top)
+  {
+    draw_skeleton();
+  }
+
+  // Set material properties
+  glDisable(GL_COLOR_MATERIAL);
+  glMaterialfv(GL_FRONT, GL_AMBIENT,
+    Vector4f(GOLD_AMBIENT[0],GOLD_AMBIENT[1],GOLD_AMBIENT[2],alpha).data());
+  glMaterialfv(GL_FRONT, GL_DIFFUSE,
+    Vector4f(GOLD_DIFFUSE[0],GOLD_DIFFUSE[1],GOLD_DIFFUSE[2],alpha).data());
+  glMaterialfv(GL_FRONT, GL_SPECULAR,
+    Vector4f(GOLD_SPECULAR[0],GOLD_SPECULAR[1],GOLD_SPECULAR[2],alpha).data());
+  glMaterialf (GL_FRONT, GL_SHININESS, 128);
+  glMaterialfv(GL_BACK, GL_AMBIENT,
+    Vector4f(SILVER_AMBIENT[0],SILVER_AMBIENT[1],SILVER_AMBIENT[2],alpha).data());
+  glMaterialfv(GL_BACK, GL_DIFFUSE,
+    Vector4f(FAST_GREEN_DIFFUSE[0],FAST_GREEN_DIFFUSE[1],FAST_GREEN_DIFFUSE[2],alpha).data());
+  glMaterialfv(GL_BACK, GL_SPECULAR,
+    Vector4f(SILVER_SPECULAR[0],SILVER_SPECULAR[1],SILVER_SPECULAR[2],alpha).data());
+  glMaterialf (GL_BACK, GL_SHININESS, 128);
+
+  if(wireframe)
+  {
+    glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
+  }
+  glLineWidth(1.0);
+  draw_mesh(V,sorted_F,sorted_N);
+  glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+
+  if(skeleton_on_top)
+  {
+    glDisable(GL_DEPTH_TEST);
+    draw_skeleton();
+  }
+
+  pop_object();
+  pop_scene();
+
+  report_gl_error();
+
+  TwDraw();
+  glutSwapBuffers();
+  if(is_animating)
+  {
+    glutPostRedisplay();
+  }
+}
+
+void mouse_wheel(int wheel, int direction, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  GLint viewport[4];
+  glGetIntegerv(GL_VIEWPORT,viewport);
+  if(wheel == 0 && TwMouseMotion(mouse_x, viewport[3] - mouse_y))
+  {
+    static double mouse_scroll_y = 0;
+    const double delta_y = 0.125*direction;
+    mouse_scroll_y += delta_y;
+    TwMouseWheel(mouse_scroll_y);
+    return;
+  }
+  push_undo();
+
+  if(wheel==0)
+  {
+    // factor of zoom change
+    double s = (1.-0.01*direction);
+    //// FOV zoom: just widen angle. This is hardly ever appropriate.
+    //camera.m_angle *= s;
+    //camera.m_angle = min(max(camera.m_angle,1),89);
+    camera.push_away(s);
+  }else
+  {
+    // Dolly zoom:
+    camera.dolly_zoom((double)direction*1.0);
+  }
+  glutPostRedisplay();
+}
+
+Eigen::VectorXi selection(const std::vector<bool> & mask)
+{
+  const int count = std::count(mask.begin(),mask.end(),true);
+  Eigen::VectorXi sel(count);
+  int s = 0;
+  for(int c = 0;c<(int)mask.size();c++)
+  {
+    if(mask[c])
+    {
+      sel(s) = c;
+      s++;
+    }
+  }
+  return sel;
+}
+
+std::vector<bool> selection_mask(const Eigen::VectorXi & sel, const int n)
+{
+  std::vector<bool> mask(n,false);
+  for(int si = 0;si<sel.size();si++)
+  {
+    const int i = sel(si);
+    mask[i] = true;
+  }
+  return mask;
+}
+
+bool ss_select(
+  const double mouse_x, 
+  const double mouse_y,
+  const Eigen::MatrixXd & C,
+  const bool accum,
+  Eigen::VectorXi & sel)
+{
+  using namespace igl;
+  using namespace Eigen;
+  using namespace std;
+  //// zap old list
+  //if(!accum)
+  //{
+  //  sel.resize(0,1);
+  //}
+
+  vector<bool> old_mask = selection_mask(s.sel,s.C.rows());
+  vector<bool> mask(old_mask.size(),false);
+
+  double min_dist = 1e25;
+  bool sel_changed = false;
+  bool found = false;
+  for(int c = 0;c<C.rows();c++)
+  {
+    const RowVector3d & Cc = C.row(c);
+    const auto Pc = project(Cc);
+    const double SELECTION_DIST = 18;//pixels
+    const double dist = (Pc.head(2)-RowVector2d(mouse_x,height-mouse_y)).norm();
+    if(dist < SELECTION_DIST && (accum || dist < min_dist))
+    {
+      mask[c] = true;
+      min_dist = dist;
+      found = true;
+      sel_changed |= mask[c] != old_mask[c];
+    }
+  }
+  for(int c = 0;c<C.rows();c++)
+  {
+    if(accum)
+    {
+      mask[c] = mask[c] ^ old_mask[c];
+    }else
+    {
+      if(!sel_changed)
+      {
+        mask[c] = mask[c] || old_mask[c];
+      }
+    }
+  }
+  sel = selection(mask);
+  return found;
+}
+
+void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  bool tw_using = TwEventMouseButtonGLUT(glutButton,glutState,mouse_x,mouse_y);
+  const int mod = (glutButton <=2 ? glutGetModifiers() : 0);
+  const bool option_down = mod & GLUT_ACTIVE_ALT;
+  const bool shift_down = mod & GLUT_ACTIVE_SHIFT;
+  const bool command_down = GLUT_ACTIVE_COMMAND & mod;
+  switch(glutButton)
+  {
+    case GLUT_RIGHT_BUTTON:
+    case GLUT_LEFT_BUTTON:
+    {
+      switch(glutState)
+      {
+        case 1:
+          // up
+          glutSetCursor(GLUT_CURSOR_INHERIT);
+          if(is_rotating)
+          {
+            sort();
+          }
+          is_rotating = false;
+          is_dragging = false;
+          break;
+        case 0:
+          new_leaf_on_drag = false;
+          new_root_on_drag = false;
+          if(!tw_using)
+          {
+            down_x = mouse_x;
+            down_y = mouse_y;
+            if(option_down)
+            {
+              glutSetCursor(GLUT_CURSOR_CYCLE);
+              // collect information for trackball
+              is_rotating = true;
+              down_camera = camera;
+            }else
+            {
+              push_undo();
+              push_scene();
+              push_object();
+              // Zap selection
+              if(shift_down)
+              {
+                s.sel.resize(0,1);
+              }
+              if(ss_select(mouse_x,mouse_y,s.C,
+                command_down && !shift_down,
+                s.sel))
+              {
+                if(shift_down)
+                {
+                  new_leaf_on_drag = true;
+                }
+              }else
+              {
+                new_root_on_drag = true;
+              }
+              is_dragging = !command_down;
+              down_C = s.C;
+              pop_object();
+              pop_scene();
+            }
+          }
+        break;
+      }
+      break;
+    }
+    // Scroll down
+    case 3:
+    {
+      mouse_wheel(0,-1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll up
+    case 4:
+    {
+      mouse_wheel(0,1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll left
+    case 5:
+    {
+      mouse_wheel(1,-1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll right
+    case 6:
+    {
+      mouse_wheel(1,1,mouse_x,mouse_y);
+      break;
+    }
+  }
+  glutPostRedisplay();
+}
+
+void mouse_drag(int mouse_x, int mouse_y)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+
+  if(is_rotating)
+  {
+    glutSetCursor(GLUT_CURSOR_CYCLE);
+    Quaterniond q;
+    switch(rotation_type)
+    {
+      case ROTATION_TYPE_IGL_TRACKBALL:
+      {
+        // Rotate according to trackball
+        igl::trackball<double>(
+          width,
+          height,
+          2.0,
+          down_camera.m_rotation_conj.coeffs().data(),
+          down_x,
+          down_y,
+          mouse_x,
+          mouse_y,
+          q.coeffs().data());
+          break;
+      }
+      case ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP:
+      {
+        // Rotate according to two axis valuator with fixed up vector
+        two_axis_valuator_fixed_up(
+          width, height,
+          2.0,
+          down_camera.m_rotation_conj,
+          down_x, down_y, mouse_x, mouse_y,
+          q);
+        break;
+      }
+      default:
+        break;
+    }
+    camera.orbit(q.conjugate());
+  }
+
+  if(is_dragging)
+  {
+    push_scene();
+    push_object();
+    if(new_leaf_on_drag)
+    {
+      assert(s.C.size() >= 1);
+      // one new node
+      s.C.conservativeResize(s.C.rows()+1,3);
+      const int nc = s.C.rows();
+      assert(s.sel.size() >= 1);
+      s.C.row(nc-1) = s.C.row(s.sel(0));
+      // one new bone
+      s.BE.conservativeResize(s.BE.rows()+1,2);
+      s.BE.row(s.BE.rows()-1) = RowVector2i(s.sel(0),nc-1);
+      // select just last node
+      s.sel.resize(1,1);
+      s.sel(0) = nc-1;
+      // reset down_C
+      down_C = s.C;
+      new_leaf_on_drag = false;
+    }
+    if(new_root_on_drag)
+    {
+      // two new nodes
+      s.C.conservativeResize(s.C.rows()+2,3);
+      const int nc = s.C.rows();
+      Vector3d obj;
+      int nhits = unproject_in_mesh(mouse_x,height-mouse_y,ei,obj);
+      if(nhits == 0)
+      {
+        Vector3d pV_mid = project(Vcen);
+        obj = unproject(Vector3d(mouse_x,height-mouse_y,pV_mid(2)));
+      }
+      s.C.row(nc-2) = obj;
+      s.C.row(nc-1) = obj;
+      // select last node
+      s.sel.resize(1,1);
+      s.sel(0) = nc-1;
+      // one new bone
+      s.BE.conservativeResize(s.BE.rows()+1,2);
+      s.BE.row(s.BE.rows()-1) = RowVector2i(nc-2,nc-1);
+      // reset down_C
+      down_C = s.C;
+      new_root_on_drag = false;
+    }
+    double z = 0;
+    Vector3d obj,win;
+    int nhits = unproject_in_mesh(mouse_x,height-mouse_y,ei,obj);
+    project(obj,win);
+    z = win(2);
+
+    for(int si = 0;si<s.sel.size();si++)
+    {
+      const int c = s.sel(si);
+      Vector3d pc = project((RowVector3d) down_C.row(c));
+      pc(0) += mouse_x-down_x;
+      pc(1) += (height-mouse_y)-(height-down_y);
+      if(nhits > 0)
+      {
+        pc(2) = z;
+      }
+      s.C.row(c) = unproject(pc);
+    }
+    pop_object();
+    pop_scene();
+  }
+
+  glutPostRedisplay();
+}
+
+void init_relative()
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+  per_face_normals(V,F,N);
+  const auto Vmax = V.colwise().maxCoeff();
+  const auto Vmin = V.colwise().minCoeff();
+  Vmid = 0.5*(Vmax + Vmin);
+  centroid(V,F,Vcen);
+  bbd = (Vmax-Vmin).norm();
+  camera.push_away(2);
+}
+
+void undo()
+{
+  using namespace std;
+  if(!undo_stack.empty())
+  {
+    redo_stack.push(s);
+    s = undo_stack.top();
+    undo_stack.pop();
+  }
+}
+
+void redo()
+{
+  using namespace std;
+  if(!redo_stack.empty())
+  {
+    undo_stack.push(s);
+    s = redo_stack.top();
+    redo_stack.pop();
+  }
+}
+
+void symmetrize()
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  if(s.sel.size() == 0)
+  {
+    cout<<YELLOWGIN("Make a selection first.")<<endl;
+    return;
+  }
+  push_undo();
+  push_scene();
+  push_object();
+  Vector3d right;
+  right_axis(right.data(),right.data()+1,right.data()+2);
+  right.normalize();
+  MatrixXd RC(s.C.rows(),s.C.cols());
+  MatrixXd old_C = s.C;
+  for(int c = 0;c<s.C.rows();c++)
+  {
+    const Vector3d Cc = s.C.row(c);
+    const auto A = Cc-Vcen;
+    const auto A1 = A.dot(right) * right;
+    const auto A2 = A-A1;
+    RC.row(c) = Vcen + A2 - A1;
+  }
+  vector<bool> mask = selection_mask(s.sel,s.C.rows());
+  // stupid O(n²) matching
+  for(int c = 0;c<s.C.rows();c++)
+  {
+    // not selected
+    if(!mask[c])
+    {
+      continue;
+    }
+    const Vector3d Cc = s.C.row(c);
+    int min_r = -1;
+    double min_dist = 1e25;
+    double max_dist = 0.1*bbd;
+    for(int r= 0;r<RC.rows();r++)
+    {
+      const Vector3d RCr = RC.row(r);
+      const double dist = (Cc-RCr).norm();
+      if(
+          dist<min_dist &&  // closest
+          dist<max_dist && // not too far away
+          (c==r || (Cc-Vcen).dot(right)*(RCr-Vcen).dot(right) > 0) // on same side
+        )
+      {
+        min_dist = dist;
+        min_r = r;
+      }
+    }
+    if(min_r>=0)
+    {
+      if(mask[min_r])
+      {
+        s.C.row(c) = 0.5*(Cc.transpose()+RC.row(min_r));
+      }else
+      {
+        s.C.row(c) = RC.row(min_r);
+      }
+    }
+  }
+  pop_object();
+  pop_scene();
+}
+
+bool save()
+{
+  using namespace std;
+  using namespace igl;
+  if(writeTGF(output_filename,s.C,s.BE))
+  {
+    cout<<GREENGIN("Current skeleton written to "+output_filename+".")<<endl;
+    return true;
+  }else
+  {
+    cout<<REDRUM("Writing to "+output_filename+" failed.")<<endl;
+    return false;
+  }
+}
+
+void key(unsigned char key, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  int mod = glutGetModifiers();
+  const bool command_down = GLUT_ACTIVE_COMMAND & mod;
+  const bool shift_down = GLUT_ACTIVE_SHIFT & mod;
+  switch(key)
+  {
+    // ESC
+    case char(27):
+      rebar.save(REBAR_NAME);
+    // ^C
+    case char(3):
+      exit(0);
+    case char(127):
+    {
+      push_undo();
+      // delete
+      MatrixXi new_BE(s.BE.rows(),s.BE.cols());
+      int count = 0;
+      for(int b=0;b<s.BE.rows();b++)
+      {
+        bool selected = false;
+        for(int si=0;si<s.sel.size();si++)
+        {
+          if(s.BE(b,0) == s.sel(si) || s.BE(b,1) == s.sel(si))
+          {
+            selected = true;
+            break;
+          }
+        }
+        if(!selected)
+        {
+          new_BE.row(count) = s.BE.row(b);
+          count++;
+        }
+      }
+      new_BE.conservativeResize(count,new_BE.cols());
+      const auto old_C = s.C;
+      VectorXi I;
+      remove_unreferenced(old_C,new_BE,s.C,s.BE,I);
+      s.sel.resize(0,1);
+      break;
+    }
+    case 'A':
+    case 'a':
+    {
+      push_undo();
+      s.sel = colon<int>(0,s.C.rows()-1);
+      break;
+    }
+    case 'C':
+    case 'c':
+    {
+      push_undo();
+      // snap close vertices
+      SparseMatrix<double> A;
+      adjacency_matrix(s.BE,A);
+      VectorXi J = colon<int>(0,s.C.rows()-1);
+      // stupid O(n²) version
+      for(int c = 0;c<s.C.rows();c++)
+      {
+        for(int d = c+1;d<s.C.rows();d++)
+        {
+          if(
+           A.coeff(c,d) == 0 &&  // no edge
+           (s.C.row(c)-s.C.row(d)).norm() < 0.02*bbd //close
+           )
+          {
+            // c < d
+            J(d) = c;
+          }
+        }
+      }
+      for(int e = 0;e<s.BE.rows();e++)
+      {
+        s.BE(e,0) = J(s.BE(e,0));
+        s.BE(e,1) = J(s.BE(e,1));
+      }
+      const auto old_BE = s.BE;
+      const auto old_C = s.C;
+      VectorXi I;
+      remove_unreferenced(old_C,old_BE,s.C,s.BE,I);
+      for(int i = 0;i<s.sel.size();i++)
+      {
+        s.sel(i) = J(s.sel(i));
+      }
+      VectorXi _;
+      igl::unique(s.sel,s.sel,_,_);
+      break;
+    }
+    case 'D':
+    case 'd':
+    {
+      push_undo();
+      s.sel.resize(0,1);
+      break;
+    }
+    case 'P':
+    case 'p':
+    {
+      // add bone to parents (should really only be one)
+      push_undo();
+      vector<int> new_sel;
+      const int old_nbe = s.BE.rows();
+      for(int si=0;si<s.sel.size();si++)
+      {
+        for(int b=0;b<old_nbe;b++)
+        {
+          if(s.BE(b,1) == s.sel(si))
+          {
+            // one new node
+            s.C.conservativeResize(s.C.rows()+1,3);
+            const int nc = s.C.rows();
+            s.C.row(nc-1) = 0.5*(s.C.row(s.BE(b,1)) + s.C.row(s.BE(b,0)));
+            // one new bone
+            s.BE.conservativeResize(s.BE.rows()+1,2);
+            s.BE.row(s.BE.rows()-1) = RowVector2i(nc-1,s.BE(b,1));
+            s.BE(b,1) = nc-1;
+            // select last node
+            new_sel.push_back(nc-1);
+          }
+        }
+      }
+      list_to_matrix(new_sel,s.sel);
+      break;
+    }
+    case 'R':
+    case 'r':
+    {
+      // re-root try at first selected
+      if(s.sel.size() > 0)
+      {
+        push_undo();
+        // only use first
+        s.sel.conservativeResize(1,1);
+        // Ideally this should only effect the connected component of s.sel(0)
+        const auto & C = s.C;
+        auto & BE = s.BE;
+        vector<bool> seen(C.rows(),false);
+        // adjacency list
+        vector<vector< int> > A;
+        adjacency_list(BE,A,false);
+        int e = 0;
+        queue<int> Q;
+        Q.push(s.sel(0));
+        seen[s.sel(0)] = true;
+        while(!Q.empty())
+        {
+          const int c = Q.front();
+          Q.pop();
+          for(const auto & d : A[c])
+          {
+            if(!seen[d])
+            {
+              BE(e,0) = c;
+              BE(e,1) = d;
+              e++;
+              Q.push(d);
+              seen[d] = true;
+            }
+          }
+        }
+        // only keep tree
+        BE.conservativeResize(e,BE.cols());
+      }
+      break;
+    }
+    case 'S':
+    case 's':
+    {
+      save();
+      break;
+    }
+    case 'U':
+    case 'u':
+    {
+      push_scene();
+      push_object();
+      for(int c = 0;c<s.C.rows();c++)
+      {
+        Vector3d P = project((Vector3d)s.C.row(c));
+        Vector3d obj;
+        int nhits = unproject_in_mesh(P(0),P(1),ei,obj);
+        if(nhits > 0)
+        {
+          s.C.row(c) = obj;
+        }
+      }
+      pop_object();
+      pop_scene();
+      break;
+    }
+    case 'Y':
+    case 'y':
+    {
+      symmetrize();
+      break;
+    }
+    case 'z':
+    case 'Z':
+      is_rotating = false;
+      is_dragging = false;
+      if(command_down)
+      {
+        if(shift_down)
+        {
+          redo();
+        }else
+        {
+          undo();
+        }
+        break;
+      }else
+      {
+        push_undo();
+        Quaterniond q;
+        snap_to_canonical_view_quat(camera.m_rotation_conj,1.0,q);
+        camera.orbit(q.conjugate());
+      }
+      break;
+    default:
+      if(!TwEventKeyboardGLUT(key,mouse_x,mouse_y))
+      {
+        cout<<"Unknown key command: "<<key<<" "<<int(key)<<endl;
+      }
+  }
+
+  glutPostRedisplay();
+}
+
+
+int main(int argc, char * argv[])
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  string filename = "../shared/decimated-knight.obj";
+  string skel_filename = "";
+  output_filename = "";
+  switch(argc)
+  {
+    case 4:
+      output_filename = argv[3];
+      //fall through
+    case 3:
+      skel_filename = argv[2];
+      if(output_filename.size() == 0)
+      {
+        output_filename = skel_filename;
+      }
+      //fall through
+    case 2:
+      // Read and prepare mesh
+      filename = argv[1];
+      break;
+    default:
+      cerr<<"Usage:"<<endl<<"    ./example input.obj [input/output.tgf]"<<endl;
+      cout<<endl<<"Opening default mesh..."<<endl;
+  }
+
+  // print key commands
+  cout<<"[Click] and [drag]     Create bone (or select node) and reposition."<<endl;
+  cout<<"⇧ +[Click] and [drag]  Select node (or create one) and _pull out_ new bone."<<endl;
+  cout<<"⌥ +[Click] and [drag]  Rotate secene."<<endl;
+  cout<<"⌫                      Delete selected node(s) and incident bones."<<endl;
+  cout<<"A,a                    Select all."<<endl;
+  cout<<"D,d                    Deselect all."<<endl;
+  cout<<"C,c                    Snap close nodes."<<endl;
+  cout<<"P,p                    Split \"parent\" bone(s) of selection by creating new node(s)."<<endl;
+  cout<<"R,r                    Breadth first search at selection to redirect skeleton into tree."<<endl;
+  cout<<"S,s                    Save current skeleton to output .tgf file."<<endl;
+  cout<<"U,u                    Project then unproject inside mesh (as if dragging each by ε)."<<endl;
+  cout<<"Y,Y                    Symmetrize selection over plane through object centroid and right vector."<<endl;
+  cout<<"Z,z                    Snap to canonical view."<<endl;
+  cout<<"⌘ Z                    Undo."<<endl;
+  cout<<"⇧ ⌘ Z                  Redo."<<endl;
+  cout<<"^C,ESC                 Exit (without saving)."<<endl;
+
+  string dir,_1,_2,name;
+  read_triangle_mesh(filename,V,F,dir,_1,_2,name);
+
+  if(output_filename.size() == 0)
+  {
+    output_filename = dir+"/"+name+".tgf";
+  }
+
+  if(file_exists(output_filename.c_str()))
+  {
+    cout<<YELLOWGIN("Output set to overwrite "<<output_filename)<<endl;
+  }else
+  {
+    cout<<BLUEGIN("Output set to "<<output_filename)<<endl;
+  }
+
+  if(skel_filename.length() > 0)
+  {
+    readTGF(skel_filename,s.C,s.BE);
+  }
+
+  init_relative();
+  ei.init(V.cast<float>(),F.cast<int>());
+
+  // Init glut
+  glutInit(&argc,argv);
+  if( !TwInit(TW_OPENGL, NULL) )
+  {
+    // A fatal error occured
+    fprintf(stderr, "AntTweakBar initialization failed: %s\n", TwGetLastError());
+    return 1;
+  }
+  // Create a tweak bar
+  rebar.TwNewBar("TweakBar");
+  rebar.TwAddVarRW("camera_rotation", TW_TYPE_QUAT4D,
+    camera.m_rotation_conj.coeffs().data(), "open readonly=true");
+  TwType RotationTypeTW = ReTwDefineEnumFromString("RotationType",
+    "igl_trackball,two-a...-fixed-up");
+  rebar.TwAddVarCB( "rotation_type", RotationTypeTW,
+    set_rotation_type,get_rotation_type,NULL,"keyIncr=] keyDecr=[");
+  rebar.TwAddVarRW("skeleton_on_top", TW_TYPE_BOOLCPP,&skeleton_on_top,"key=O");
+  rebar.TwAddVarRW("wireframe", TW_TYPE_BOOLCPP,&wireframe,"key=l");
+  TwType SkelStyleTypeTW = ReTwDefineEnumFromString("SkelStyleType",
+    "3d,vector-graphics");
+  rebar.TwAddVarRW("style",SkelStyleTypeTW,&skel_style,"");
+  rebar.TwAddVarRW("alpha",TW_TYPE_DOUBLE,&alpha,
+    "keyIncr=} keyDecr={ min=0 max=1 step=0.1");
+  rebar.load(REBAR_NAME);
+
+  // Init antweakbar
+  glutInitDisplayString( "rgba depth double samples>=8 ");
+  glutInitWindowSize(glutGet(GLUT_SCREEN_WIDTH)/2.0,glutGet(GLUT_SCREEN_HEIGHT)/2.0);
+  glutCreateWindow("skeleton-builder");
+  glutDisplayFunc(display);
+  glutReshapeFunc(reshape);
+  glutKeyboardFunc(key);
+  glutMouseFunc(mouse);
+  glutMotionFunc(mouse_drag);
+  glutPassiveMotionFunc((GLUTmousemotionfun)TwEventMouseMotionGLUT);
+  glutMainLoop();
+
+  return 0;
+}

examples/skeleton-posing/Makefile

@@ -0,0 +1,83 @@
+
+.PHONY: all
+
+# Shared flags etc.
+include ../../build/Makefile.conf
+
+all: obj example
+
+.PHONY: example
+
+LIBIGL=../../
+LIBIGL_INC=-I$(LIBIGL)/include
+LIBIGL_LIB=-L$(LIBIGL)/lib -ligl -liglembree -liglcgal -ligltetgen -liglbbw -liglmosek
+
+EIGEN3_INC=-I/opt/local/include/eigen3 -I/opt/local/include/eigen3/unsupported
+
+ANTTWEAKBAR_INC=-I$(LIBIGL)/external/AntTweakBar/include
+ANTTWEAKBAR_LIB=-L$(LIBIGL)/external/AntTweakBar/lib -lAntTweakBar -framework AppKit
+
+TETGEN=$(LIBIGL)/external/tetgen
+TETGEN_LIB=-L$(TETGEN) -ltet 
+TETGEN_INC=-I$(TETGEN)
+
+EMBREE=$(LIBIGL)/external/embree
+EMBREE_INC=-I$(EMBREE)/ -I$(EMBREE)/embree
+EMBREE_LIB=-L$(EMBREE)/build -lembree -lsys
+
+CARBON_LIB=-framework Carbon
+
+# Use free glut for mouse scrolling
+#FREE_GLUT=/opt/local/
+#FREE_GLUT_INC=-I$(FREE_GLUT)/include
+#FREE_GLUT_LIB=-L$(FREE_GLUT)/lib -lglut
+GLUT_LIB=-framework GLUT
+GLUT_INC=
+
+ifdef IGL_NO_MOSEK
+CFLAGS+=-DIGL_NO_MOSEK
+else
+# Adjust your mosek paths etc. accordingly
+ifndef MOSEKPLATFORM
+  MOSEKPLATFORM=osx64x86
+endif
+ifndef MOSEKVERSION
+  MOSEKVERSION=7
+endif
+IGLMOSEK=../mosek/
+IGLMOSEK_INC=-I$(IGLMOSEK)/
+MOSEK=/usr/local/mosek
+MOSEK_INC=-I$(MOSEK)/$(MOSEKVERSION)/tools/platform/$(MOSEKPLATFORM)/h
+MOSEK_LIB=-L$(MOSEK)/$(MOSEKVERSION)/tools/platform/$(MOSEKPLATFORM)/bin -lmosek64
+endif
+
+CGAL=/opt/local/
+CGAL_LIB=-L$(CGAL)/lib -lCGAL -lCGAL_Core -lgmp -lmpfr -lboost_thread-mt -lboost_system-mt
+CGAL_INC=-I$(CGAL)/include -I/usr/include/
+# This is absolutely necessary for Exact Construction
+CGAL_FLAGS=-frounding-math -fsignaling-nans 
+CFLAGS+=$(CGAL_FLAGS)
+
+INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(MATLAB_INC) $(GLUT_INC) ${CGAL_INC} ${TETGEN_INC} $(MOSEK_INC) $(EMBREE_INC)
+LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(MATLAB_LIB) ${CGAL_LIB} $(CARBON_LIB) $(TETGEN_LIB) $(MOSEK_LIB) $(EMBREE_LIB)
+
+CPP_FILES=$(wildcard ./*.cpp)
+OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o))) 
+
+CFLAGS+=-std=c++11 -g
+
+example: obj $(OBJ_FILES)
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o example $(LIB) $(OBJ_FILES) 
+
+obj:
+	mkdir -p obj
+
+obj/%.o: %.cpp
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
+
+obj/%.o: %.cpp %.h
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -c $< -o $@ $(INC)
+
+clean:
+	rm -f $(OBJ_FILES)
+	rm -f example

examples/skeleton-posing/example.cpp

@@ -0,0 +1,953 @@
+#include <igl/Camera.h>
+#include <igl/MouseController.h>
+#include <igl/REDRUM.h>
+#include <igl/ReAntTweakBar.h>
+#include <igl/STR.h>
+#include <igl/barycenter.h>
+#include <igl/bbw/bbw.h>
+#include <igl/bone_parents.h>
+#include <igl/boundary_conditions.h>
+#include <igl/boundary_facets.h>
+#include <igl/centroid.h>
+#include <igl/cgal/remesh_self_intersections.h>
+#include <igl/colon.h>
+#include <igl/draw_beach_ball.h>
+#include <igl/draw_floor.h>
+#include <igl/draw_mesh.h>
+#include <igl/draw_skeleton_3d.h>
+#include <igl/draw_skeleton_vector_graphics.h>
+#include <igl/forward_kinematics.h>
+#include <igl/get_seconds.h>
+#include <igl/lbs_matrix.h>
+#include <igl/material_colors.h>
+#include <igl/normalize_row_sums.h>
+#include <igl/pathinfo.h>
+#include <igl/per_face_normals.h>
+#include <igl/quat_to_mat.h>
+#include <igl/readDMAT.h>
+#include <igl/readTGF.h>
+#include <igl/read_triangle_mesh.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/report_gl_error.h>
+#include <igl/snap_to_canonical_view_quat.h>
+#include <igl/snap_to_fixed_up.h>
+#include <igl/tetgen/mesh_with_skeleton.h>
+#include <igl/tetgen/tetrahedralize.h>
+#include <igl/trackball.h>
+#include <igl/two_axis_valuator_fixed_up.h>
+#include <igl/winding_number.h>
+#include <igl/writeDMAT.h>
+#include <igl/writeOBJ.h>
+#include <igl/writeMESH.h>
+#include <igl/writeOFF.h>
+#include <igl/writeTGF.h>
+#include <igl/next_filename.h>
+#include <igl/volume.h>
+
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+#ifdef __APPLE__
+#include <GLUT/glut.h>
+#include <Carbon/Carbon.h>
+#else
+#include <GL/glut.h>
+#endif
+
+#include <string>
+#include <vector>
+#include <queue>
+#include <stack>
+#include <iostream>
+#include <iomanip>
+
+#define VERBOSE
+
+enum SkelStyleType
+{
+  SKEL_STYLE_TYPE_3D = 0,
+  SKEL_STYLE_TYPE_VECTOR_GRAPHICS = 1,
+  NUM_SKEL_STYLE_TYPE = 2
+}skel_style;
+
+Eigen::MatrixXd V,N,W,M;
+Eigen::Vector3d Vmid;
+double bbd = 1.0;
+Eigen::MatrixXi F;
+igl::Camera camera;
+
+Eigen::MatrixXd C;
+Eigen::MatrixXi BE;
+Eigen::VectorXi P,RP;
+
+struct State
+{
+  igl::MouseController mouse;
+  Eigen::MatrixXf colors;
+} s;
+
+bool wireframe = false;
+
+// See README for descriptions
+enum RotationType
+{
+  ROTATION_TYPE_IGL_TRACKBALL = 0,
+  ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP = 1,
+  NUM_ROTATION_TYPES = 2,
+} rotation_type = ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP;
+
+std::stack<State> undo_stack;
+std::stack<State> redo_stack;
+
+bool is_rotating = false;
+bool centroid_is_visible = true;
+int down_x,down_y;
+igl::Camera down_camera;
+std::string output_prefix;
+
+struct CameraAnimation
+{
+  bool is_animating = false;
+  double DURATION = 0.5;
+  double start_time = 0;
+  Eigen::Quaterniond from_quat,to_quat;
+} canim;
+
+typedef std::vector<
+Eigen::Quaterniond,
+  Eigen::aligned_allocator<Eigen::Quaterniond> > RotationList;
+
+struct PoseAnimation
+{
+  bool is_animating = false;
+  double DURATION = 2;
+  double start_time = 0;
+  RotationList pose;
+} panim;
+
+int width,height;
+Eigen::Vector4f light_pos(-0.1,-0.1,0.9,0);
+
+#define REBAR_NAME "temp.rbr"
+igl::ReTwBar rebar;
+
+void push_undo()
+{
+  undo_stack.push(s);
+  // Clear
+  redo_stack = std::stack<State>();
+}
+
+// No-op setter, does nothing
+void TW_CALL no_op(const void * /*value*/, void * /*clientData*/)
+{
+}
+
+void TW_CALL set_rotation_type(const void * value, void * clientData)
+{
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  const RotationType old_rotation_type = rotation_type;
+  rotation_type = *(const RotationType *)(value);
+  if(rotation_type == ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP &&
+    old_rotation_type != ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP)
+  {
+    push_undo();
+    canim.from_quat = camera.m_rotation_conj;
+    snap_to_fixed_up(canim.from_quat,canim.to_quat);
+    // start animation
+    canim.start_time = get_seconds();
+    canim.is_animating = true;
+  }
+}
+void TW_CALL get_rotation_type(void * value, void *clientData)
+{
+  RotationType * rt = (RotationType *)(value);
+  *rt = rotation_type;
+}
+
+void reshape(int width, int height)
+{
+  ::width = width;
+  ::height = height;
+  glViewport(0,0,width,height);
+  // Send the new window size to AntTweakBar
+  TwWindowSize(width, height);
+  camera.m_aspect = (double)width/(double)height;
+  s.mouse.reshape(width,height);
+}
+
+void push_scene()
+{
+  using namespace igl;
+  using namespace std;
+  glMatrixMode(GL_PROJECTION);
+  glPushMatrix();
+  glLoadIdentity();
+  gluPerspective(camera.m_angle,camera.m_aspect,camera.m_near,camera.m_far);
+  glMatrixMode(GL_MODELVIEW);
+  glPushMatrix();
+  glLoadIdentity();
+  gluLookAt(
+    camera.eye()(0), camera.eye()(1), camera.eye()(2),
+    camera.at()(0), camera.at()(1), camera.at()(2),
+    camera.up()(0), camera.up()(1), camera.up()(2));
+}
+
+void push_object()
+{
+  using namespace igl;
+  glPushMatrix();
+  glScaled(2./bbd,2./bbd,2./bbd);
+  glTranslated(-Vmid(0),-Vmid(1),-Vmid(2));
+}
+
+void pop_object()
+{
+  glPopMatrix();
+}
+
+void pop_scene()
+{
+  glMatrixMode(GL_PROJECTION);
+  glPopMatrix();
+  glMatrixMode(GL_MODELVIEW);
+  glPopMatrix();
+}
+
+// Set up double-sided lights
+void lights()
+{
+  using namespace std;
+  using namespace Eigen;
+  glEnable(GL_LIGHTING);
+  glLightModelf(GL_LIGHT_MODEL_TWO_SIDE,GL_TRUE);
+  glEnable(GL_LIGHT0);
+  glEnable(GL_LIGHT1);
+  float WHITE[4] =  {0.8,0.8,0.8,1.};
+  float GREY[4] =  {0.4,0.4,0.4,1.};
+  float BLACK[4] =  {0.,0.,0.,1.};
+  Vector4f pos = light_pos;
+  glLightfv(GL_LIGHT0,GL_AMBIENT,GREY);
+  glLightfv(GL_LIGHT0,GL_DIFFUSE,WHITE);
+  glLightfv(GL_LIGHT0,GL_SPECULAR,BLACK);
+  glLightfv(GL_LIGHT0,GL_POSITION,pos.data());
+  pos(0) *= -1;
+  pos(1) *= -1;
+  pos(2) *= -1;
+  glLightfv(GL_LIGHT1,GL_AMBIENT,GREY);
+  glLightfv(GL_LIGHT1,GL_DIFFUSE,WHITE);
+  glLightfv(GL_LIGHT1,GL_SPECULAR,BLACK);
+  glLightfv(GL_LIGHT1,GL_POSITION,pos.data());
+}
+
+void display()
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+  const float back[4] = {0.75, 0.75, 0.75,0};
+  glClearColor(back[0],back[1],back[2],0);
+  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+  if(canim.is_animating)
+  {
+    double t = (get_seconds() - canim.start_time)/canim.DURATION;
+    if(t > 1)
+    {
+      t = 1;
+      canim.is_animating = false;
+    }
+    Quaterniond q = canim.from_quat.slerp(t,canim.to_quat).normalized();
+    camera.orbit(q.conjugate());
+  }
+
+  glEnable(GL_DEPTH_TEST);
+  glDepthFunc(GL_LEQUAL);
+  glEnable(GL_NORMALIZE);
+  lights();
+  push_scene();
+
+  // Draw a nice floor
+  glEnable(GL_DEPTH_TEST);
+  glPushMatrix();
+  const double floor_offset =
+    -2./bbd*(V.col(1).maxCoeff()-Vmid(1));
+  glTranslated(0,floor_offset,0);
+  const float GREY[4] = {0.5,0.5,0.6,1.0};
+  const float DARK_GREY[4] = {0.2,0.2,0.3,1.0};
+  glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+  glEnable(GL_CULL_FACE);
+  glCullFace(GL_BACK);
+  draw_floor(GREY,DARK_GREY);
+  glDisable(GL_CULL_FACE);
+  glPopMatrix();
+
+  push_object();
+
+  const auto & draw_skeleton = [](const MatrixXd & T)
+  {
+    switch(skel_style)
+    {
+      default:
+      case SKEL_STYLE_TYPE_3D:
+      {
+        draw_skeleton_3d(C,BE,T,s.colors);
+        break;
+      }
+      case SKEL_STYLE_TYPE_VECTOR_GRAPHICS:
+        draw_skeleton_vector_graphics(C,BE,T);
+        break;
+    }
+  };
+  // Set material properties
+  glDisable(GL_COLOR_MATERIAL);
+  glMaterialfv(GL_FRONT, GL_AMBIENT,GOLD_AMBIENT);
+  glMaterialfv(GL_FRONT, GL_DIFFUSE,GOLD_DIFFUSE);
+  glMaterialfv(GL_FRONT, GL_SPECULAR,GOLD_SPECULAR);
+  glMaterialf (GL_FRONT, GL_SHININESS, 128);
+  glMaterialfv(GL_BACK, GL_AMBIENT,SILVER_AMBIENT);
+  glMaterialfv(GL_BACK, GL_DIFFUSE,FAST_GREEN_DIFFUSE);
+  glMaterialfv(GL_BACK, GL_SPECULAR,SILVER_SPECULAR);
+  glMaterialf (GL_BACK, GL_SHININESS, 128);
+  if(wireframe)
+  {
+    glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
+  }
+  glLineWidth(1.0);
+  MatrixXd T;
+  RotationList dQ;
+  if(panim.is_animating)
+  {
+    double t = (get_seconds() - panim.start_time)/panim.DURATION;
+    if(t > 1)
+    {
+      t = 1;
+      panim.is_animating = false;
+    }
+    const auto & ease = [](const double t)
+    {
+      return 3.*t*t-2.*t*t*t;
+    };
+    double f = (t<0.5?ease(2.*t):ease(2.-2.*t));
+    dQ.resize(panim.pose.size());
+    for(int e = 0;e<(int)panim.pose.size();e++)
+    {
+      dQ[e] = panim.pose[e].slerp(f,Quaterniond::Identity()).normalized();
+    }
+  }else
+  {
+    dQ = s.mouse.rotations();
+  }
+  forward_kinematics(C,BE,P,dQ,T);
+  MatrixXd U = M*T;
+  MatrixXd UN;
+  per_face_normals(U,F,UN);
+  draw_mesh(U,F,UN);
+  glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+  glDisable(GL_DEPTH_TEST);
+  draw_skeleton(T);
+
+  if(centroid_is_visible)
+  {
+    Vector3d cen;
+    centroid(U,F,cen);
+    glEnable(GL_DEPTH_TEST);
+    glPushMatrix();
+    glTranslated(cen(0),cen(1),cen(2));
+    glScaled(bbd/2.,bbd/2.,bbd/2.);
+    glScaled(0.1,0.1,0.1);
+    glEnable(GL_POLYGON_OFFSET_FILL);
+    glPolygonOffset(0,-100000);
+    draw_beach_ball();
+    glDisable(GL_POLYGON_OFFSET_FILL);
+    glPopMatrix();
+  }
+
+  // Mouse is always on top
+  glDisable(GL_DEPTH_TEST);
+  if(!panim.is_animating)
+  {
+    s.mouse.draw();
+  }
+  pop_object();
+  pop_scene();
+
+  report_gl_error();
+
+  TwDraw();
+  glutSwapBuffers();
+  if(canim.is_animating || panim.is_animating)
+  {
+    glutPostRedisplay();
+  }
+}
+
+void mouse_wheel(int wheel, int direction, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  GLint viewport[4];
+  glGetIntegerv(GL_VIEWPORT,viewport);
+  if(wheel == 0 && TwMouseMotion(mouse_x, viewport[3] - mouse_y))
+  {
+    static double mouse_scroll_y = 0;
+    const double delta_y = 0.125*direction;
+    mouse_scroll_y += delta_y;
+    TwMouseWheel(mouse_scroll_y);
+    return;
+  }
+  push_undo();
+
+  if(wheel==0)
+  {
+    // factor of zoom change
+    double s = (1.-0.01*direction);
+    //// FOV zoom: just widen angle. This is hardly ever appropriate.
+    //camera.m_angle *= s;
+    //camera.m_angle = min(max(camera.m_angle,1),89);
+    camera.push_away(s);
+  }else
+  {
+    // Dolly zoom:
+    camera.dolly_zoom((double)direction*1.0);
+  }
+  glutPostRedisplay();
+}
+
+void mouse(int glutButton, int glutState, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  bool tw_using = TwEventMouseButtonGLUT(glutButton,glutState,mouse_x,mouse_y);
+  const int mod = (glutButton <=2 ? glutGetModifiers() : 0);
+  const bool option_down = mod & GLUT_ACTIVE_ALT;
+  switch(glutButton)
+  {
+    case GLUT_RIGHT_BUTTON:
+    case GLUT_LEFT_BUTTON:
+    {
+      push_scene();
+      push_object();
+      switch(glutState)
+      {
+        case 1:
+        {
+          // up
+          const bool mouse_was_selecting = s.mouse.is_selecting();
+          is_rotating = false;
+          s.mouse.up(mouse_x,mouse_y);
+          glutSetCursor(GLUT_CURSOR_INHERIT);
+          if(mouse_was_selecting)
+          {
+            s.mouse.set_selection_from_last_drag(C,BE,P,RP);
+            MouseController::VectorXb S;
+            MouseController::propogate_to_descendants_if(
+              s.mouse.selection(),P,S);
+            MouseController::color_if(S,MAYA_SEA_GREEN,MAYA_VIOLET,s.colors);
+          }
+          break;
+        }
+        case 0:
+          if(!tw_using)
+          {
+            down_x = mouse_x;
+            down_y = mouse_y;
+            if(option_down || glutButton==GLUT_RIGHT_BUTTON)
+            {
+              glutSetCursor(GLUT_CURSOR_CYCLE);
+              // collect information for trackball
+              is_rotating = true;
+              down_camera = camera;
+            }else
+            {
+              push_undo();
+              s.mouse.down(mouse_x,mouse_y);
+            }
+          }
+        break;
+      }
+      pop_object();
+      pop_scene();
+      break;
+    }
+    // Scroll down
+    case 3:
+    {
+      mouse_wheel(0,-1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll up
+    case 4:
+    {
+      mouse_wheel(0,1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll left
+    case 5:
+    {
+      mouse_wheel(1,-1,mouse_x,mouse_y);
+      break;
+    }
+    // Scroll right
+    case 6:
+    {
+      mouse_wheel(1,1,mouse_x,mouse_y);
+      break;
+    }
+  }
+  glutPostRedisplay();
+}
+
+void mouse_drag(int mouse_x, int mouse_y)
+{
+  using namespace igl;
+  using namespace std;
+  using namespace Eigen;
+
+  push_scene();
+  push_object();
+  if(is_rotating)
+  {
+    glutSetCursor(GLUT_CURSOR_CYCLE);
+    Quaterniond q;
+    switch(rotation_type)
+    {
+      case ROTATION_TYPE_IGL_TRACKBALL:
+      {
+        // Rotate according to trackball
+        igl::trackball(
+          width,
+          height,
+          2.0,
+          down_camera.m_rotation_conj,
+          down_x,
+          down_y,
+          mouse_x,
+          mouse_y,
+          q);
+          break;
+      }
+      case ROTATION_TYPE_TWO_AXIS_VALUATOR_FIXED_UP:
+      {
+        // Rotate according to two axis valuator with fixed up vector
+        two_axis_valuator_fixed_up(
+          width, height,
+          2.0,
+          down_camera.m_rotation_conj,
+          down_x, down_y, mouse_x, mouse_y,
+          q);
+        break;
+      }
+      default:
+        break;
+    }
+    camera.orbit(q.conjugate());
+  }else if(s.mouse.drag(mouse_x,mouse_y))
+  {
+  }
+  pop_object();
+  pop_scene();
+  glutPostRedisplay();
+}
+
+void init_relative()
+{
+  using namespace Eigen;
+  using namespace igl;
+  using namespace std;
+  per_face_normals(V,F,N);
+  const auto Vmax = V.colwise().maxCoeff();
+  const auto Vmin = V.colwise().minCoeff();
+  Vmid = 0.5*(Vmax + Vmin);
+  bbd = (Vmax-Vmin).norm();
+  camera.push_away(2);
+}
+
+void undo()
+{
+  using namespace std;
+  if(!undo_stack.empty())
+  {
+    redo_stack.push(s);
+    s = undo_stack.top();
+    undo_stack.pop();
+    s.mouse.reshape(width,height);
+  }
+}
+
+void redo()
+{
+  using namespace std;
+  if(!redo_stack.empty())
+  {
+    undo_stack.push(s);
+    s = redo_stack.top();
+    redo_stack.pop();
+    s.mouse.reshape(width,height);
+  }
+}
+
+bool save()
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  string output_filename;
+  next_filename(output_prefix,4,".dmat",output_filename);
+  MatrixXd T;
+  forward_kinematics(C,BE,P,s.mouse.rotations(),T);
+  if(writeDMAT(output_filename,T))
+  {
+    cout<<GREENGIN("Current pose written to "+output_filename+".")<<endl;
+    return true;
+  }else
+  {
+    cout<<REDRUM("Writing to "+output_filename+" failed.")<<endl;
+    return false;
+  }
+}
+
+void key(unsigned char key, int mouse_x, int mouse_y)
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  int mod = glutGetModifiers();
+  const bool command_down = GLUT_ACTIVE_COMMAND & mod;
+  const bool shift_down = GLUT_ACTIVE_SHIFT & mod;
+  switch(key)
+  {
+    // ESC
+    case char(27):
+      rebar.save(REBAR_NAME);
+    // ^C
+    case char(3):
+      exit(0);
+    case 'A':
+    case 'a':
+    {
+      panim.is_animating = !panim.is_animating;
+      panim.pose = s.mouse.rotations();
+      panim.start_time = get_seconds();
+      break;
+    }
+    case 'D':
+    case 'd':
+    {
+      push_undo();
+      s.mouse.clear_selection();
+      break;
+    }
+    case 'R':
+    {
+      push_undo();
+      s.mouse.reset_selected_rotations();
+      break;
+    }
+    case 'r':
+    {
+      push_undo();
+      s.mouse.reset_rotations();
+      break;
+    }
+    case 'S':
+    case 's':
+    {
+      save();
+      break;
+    }
+    case 'z':
+    case 'Z':
+      is_rotating = false;
+      if(command_down)
+      {
+        if(shift_down)
+        {
+          redo();
+        }else
+        {
+          undo();
+        }
+        break;
+      }else
+      {
+        push_undo();
+        Quaterniond q;
+        snap_to_canonical_view_quat(camera.m_rotation_conj,1.0,q);
+        camera.orbit(q.conjugate());
+      }
+      break;
+    default:
+      if(!TwEventKeyboardGLUT(key,mouse_x,mouse_y))
+      {
+        cout<<"Unknown key command: "<<key<<" "<<int(key)<<endl;
+      }
+  }
+
+  glutPostRedisplay();
+}
+
+bool clean(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXd & CV,
+  Eigen::MatrixXi & CF)
+{
+  using namespace igl;
+  using namespace Eigen;
+  using namespace std;
+  {
+    MatrixXi _1;
+    VectorXi _2,IM;
+#ifdef VERBOSE
+    cout<<"remesh_self_intersections"<<endl;
+#endif
+    remesh_self_intersections(V,F,{},CV,CF,_1,_2,IM);
+    for_each(CF.data(),CF.data()+CF.size(),[&IM](int & a){a=IM(a);});
+    MatrixXd oldCV = CV;
+    MatrixXi oldCF = CF;
+    remove_unreferenced(oldCV,oldCF,CV,CF,IM);
+  }
+  MatrixXd TV;
+  MatrixXi TT;
+  {
+    MatrixXi _1;
+    // c  convex hull
+    // Y  no boundary steiners
+    // p  polygon input
+#ifdef VERBOSE
+    cout<<"tetrahedralize"<<endl;
+#endif
+    if(tetrahedralize(CV,CF,"cYpC",TV,TT,_1) != 0)
+    {
+      cout<<REDRUM("CDT failed.")<<endl;
+      return false;
+    }
+  }
+  {
+    MatrixXd BC;
+    barycenter(TV,TT,BC);
+    VectorXd W;
+#ifdef VERBOSE
+    cout<<"winding_number"<<endl;
+#endif
+    winding_number(V,F,BC,W);
+    W = W.array().abs();
+    const double thresh = 0.5;
+    const int count = (W.array()>thresh).cast<int>().sum();
+    MatrixXi CT(count,TT.cols());
+    int c = 0;
+    for(int t = 0;t<TT.rows();t++)
+    {
+      if(W(t)>thresh)
+      {
+        CT.row(c++) = TT.row(t);
+      }
+    }
+    assert(c==count);
+    boundary_facets(CT,CF);
+  }
+  return true;
+}
+
+bool robust_weights(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & C,
+  const Eigen::MatrixXi & BE,
+  Eigen::MatrixXd & W)
+{
+  using namespace igl;
+  using namespace Eigen;
+  using namespace std;
+  // clean mesh
+  MatrixXd CV;
+  MatrixXi CF;
+  if(!clean(V,F,CV,CF))
+  {
+    return false;
+  }
+  MatrixXd TV;
+  MatrixXi TT;
+  // compute tet-mesh
+  {
+    MatrixXi _1;
+#ifdef VERBOSE
+    cout<<"mesh_with_skeleton"<<endl;
+#endif
+    if(!mesh_with_skeleton(CV,CF,C,{},BE,{},10,"pq1.5Y",TV,TT,_1))
+    {
+      cout<<REDRUM("tetgen failed.")<<endl;
+      return false;
+    }
+  }
+  // Finally, tetgen may have still included some insanely small tets.
+  // Just ignore these during weight computation (and hope they don't isolate
+  // any vertices).
+  {
+    const MatrixXi oldTT = TT;
+    VectorXd vol;
+    volume(TV,TT,vol);
+    const double thresh = 1e-17;
+    const int count = (vol.array()>thresh).cast<int>().sum();
+    TT.resize(count,oldTT.cols());
+    int c = 0;
+    for(int t = 0;t<oldTT.rows();t++)
+    {
+      if(vol(t)>thresh)
+      {
+        TT.row(c++) = oldTT.row(t);
+      }
+    }
+  }
+
+  // compute weights
+  VectorXi b;
+  MatrixXd bc;
+  if(!boundary_conditions(TV,TT,C,{},BE,{},b,bc))
+  {
+    cout<<REDRUM("boundary_conditions failed.")<<endl;
+    return false;
+  }
+  // compute BBW
+  // Default bbw data and flags
+  BBWData bbw_data;
+  bbw_data.verbosity = 1;
+#ifdef IGL_NO_MOSEK
+  bbw_data.qp_solver = QP_SOLVER_IGL_ACTIVE_SET;
+  bbw_data.active_set_params.max_iter = 4;
+#else
+  bbw_data.mosek_data.douparam[MSK_DPAR_INTPNT_TOL_REL_GAP]=1e-14;
+  bbw_data.qp_solver = QP_SOLVER_MOSEK;
+#endif
+  // Weights matrix
+  if(!bbw(TV,TT,b,bc,bbw_data,W))
+  {
+    return false;
+  }
+  // Normalize weights to sum to one
+  normalize_row_sums(W,W);
+  // keep first #V weights
+  W.conservativeResize(V.rows(),W.cols());
+  return true;
+}
+
+
+int main(int argc, char * argv[])
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  string filename = "../shared/cheburashka.off";
+  string skel_filename = "../shared/cheburashka.tgf";
+  string weights_filename = "";
+  output_prefix = "";
+  switch(argc)
+  {
+    case 5:
+      output_prefix = argv[4];
+      //fall through
+    case 4:
+      weights_filename = argv[3];
+      //fall through
+    case 3:
+      skel_filename = argv[2];
+      // Read and prepare mesh
+      filename = argv[1];
+      break;
+    default:
+      cerr<<"Usage:"<<endl<<"    ./example model.obj skeleton.tgf [weights.dmat] [pose-prefix]"<<endl;
+      cout<<endl<<"Opening default rig..."<<endl;
+  }
+
+  // print key commands
+  cout<<"[Click] and [drag]     Select a bone/Use onscreen widget to rotate bone."<<endl;
+  cout<<"⌥ +[Click] and [drag]  Rotate secene."<<endl;
+  cout<<"⌫                      Delete selected node(s) and incident bones."<<endl;
+  cout<<"D,d                    Deselect all."<<endl;
+  cout<<"S,s                    Save current pose."<<endl;
+  cout<<"R                      Reset selected rotation."<<endl;
+  cout<<"r                      Reset all rotations."<<endl;
+  cout<<"Z,z                    Snap to canonical view."<<endl;
+  cout<<"⌘ Z                    Undo."<<endl;
+  cout<<"⇧ ⌘ Z                  Redo."<<endl;
+  cout<<"^C,ESC                 Exit (without saving)."<<endl;
+
+  string dir,_1,_2,name;
+  read_triangle_mesh(filename,V,F,dir,_1,_2,name);
+
+  if(output_prefix.size() == 0)
+  {
+    output_prefix = dir+"/"+name+"-pose-";
+  }
+
+  {
+    string output_filename;
+    next_filename(output_prefix,4,".dmat",output_filename);
+    cout<<BLUEGIN("Output set to start with "<<output_filename)<<endl;
+  }
+
+  // Read in skeleton and precompute hierarchy
+  readTGF(skel_filename,C,BE);
+  // initialize mouse interface
+  s.mouse.set_size(BE.rows());
+  // Rigid parts (not used)
+  colon<int>(0,BE.rows()-1,RP);
+  assert(RP.size() == BE.rows());
+  // Bone parents
+  bone_parents(BE,P);
+  if(weights_filename.size() == 0)
+  {
+    robust_weights(V,F,C,BE,W);
+  }else
+  {
+    // Read in weights and precompute LBS matrix
+    readDMAT(weights_filename,W);
+  }
+  lbs_matrix(V,W,M);
+
+  init_relative();
+
+  // Init glut
+  glutInit(&argc,argv);
+  if( !TwInit(TW_OPENGL, NULL) )
+  {
+    // A fatal error occured
+    fprintf(stderr, "AntTweakBar initialization failed: %s\n", TwGetLastError());
+    return 1;
+  }
+  // Create a tweak bar
+  rebar.TwNewBar("TweakBar");
+  rebar.TwAddVarRW("camera_rotation", TW_TYPE_QUAT4D,
+    camera.m_rotation_conj.coeffs().data(), "open readonly=true");
+  TwType RotationTypeTW = ReTwDefineEnumFromString("RotationType",
+    "igl_trackball,two-a...-fixed-up");
+  rebar.TwAddVarCB( "rotation_type", RotationTypeTW,
+    set_rotation_type,get_rotation_type,NULL,"keyIncr=] keyDecr=[");
+  rebar.TwAddVarRW("wireframe", TW_TYPE_BOOLCPP,&wireframe,"key=l");
+  rebar.TwAddVarRW("centroid_is_visible", TW_TYPE_BOOLCPP,&centroid_is_visible,
+    "keyIncr=C keyDecr=c label='centroid visible?'");
+  TwType SkelStyleTypeTW = ReTwDefineEnumFromString("SkelStyleType",
+    "3d,vector-graphics");
+  rebar.TwAddVarRW("style",SkelStyleTypeTW,&skel_style,"");
+  rebar.load(REBAR_NAME);
+
+  // Init antweakbar
+  glutInitDisplayString( "rgba depth double samples>=8 ");
+  glutInitWindowSize(glutGet(GLUT_SCREEN_WIDTH)/2.0,glutGet(GLUT_SCREEN_HEIGHT)/2.0);
+  glutCreateWindow("skeleton-poser");
+  glutDisplayFunc(display);
+  glutReshapeFunc(reshape);
+  glutKeyboardFunc(key);
+  glutMouseFunc(mouse);
+  glutMotionFunc(mouse_drag);
+  glutPassiveMotionFunc((GLUTmousemotionfun)TwEventMouseMotionGLUT);
+  glutMainLoop();
+
+  return 0;
+}
+

examples/skeleton/Makefile

@@ -38,7 +38,7 @@ ifndef MOSEKPLATFORM
   MOSEKPLATFORM=osx64x86
 endif
 ifndef MOSEKVERSION
-  MOSEKVERSION=6
+  MOSEKVERSION=7
 endif
 IGLMOSEK=../mosek/
 IGLMOSEK_INC=-I$(IGLMOSEK)/

include/igl/InElementAABB.h

@@ -280,13 +280,16 @@ inline std::vector<int> igl::InElementAABB::find(
   const int dim = m_bb_max.size();
   assert(q.size() == m_bb_max.size());
   const double epsilon = 1e-14;
+  // Check if outside bounding box
   for(int d = 0;d<q.size()&&inside;d++)
   {
     inside &= (q(d)-m_bb_min(d))>=epsilon;
     inside &= (m_bb_max(d)-q(d))>=epsilon;
   }
+  cout<<"searching..."<<endl;
   if(!inside)
   {
+    cout<<"not in bb"<<endl;
     return std::vector<int>();
   }
   if(m_element != -1)
@@ -314,9 +317,16 @@ inline std::vector<int> igl::InElementAABB::find(
           const Vector2d V1 = V.row(Ele(m_element,0));
           const Vector2d V2 = V.row(Ele(m_element,1));
           const Vector2d V3 = V.row(Ele(m_element,2));
+          double a0 = doublearea_single(V1,V2,V3);
           a1 = doublearea_single(V1,V2,(Vector2d)q);
           a2 = doublearea_single(V2,V3,(Vector2d)q);
           a3 = doublearea_single(V3,V1,(Vector2d)q);
+          cout<<
+            a0<<" "<<
+            a1<<" "<<
+            a2<<" "<<
+            a3<<" "<<
+            endl;
           break;
         }
       default:assert(false);

include/igl/MouseController.h

@@ -80,6 +80,7 @@ namespace igl
       inline const VectorXb & selection() const{return m_selection;};
       //                          〃 m_is_selecting
       inline const bool & is_selecting() const{return m_is_selecting;}
+      inline bool is_widget_down() const{return m_widget.is_down();}
       //                          〃 m_rotations
       inline const RotationList & rotations() const{return m_rotations;}
       // Returns non-const reference to m_root_enabled
@@ -122,6 +123,7 @@ namespace igl
       inline void set_size(const int n);
       // Resets m_rotation elements to identity
       inline void reset_rotations();
+      inline void reset_selected_rotations();
       inline bool set_rotations(const RotationList & vQ);
       // Sets all entries in m_selection to false
       inline void clear_selection();
@@ -500,6 +502,19 @@ inline void igl::MouseController::reset_rotations()
   // cop out. just clear selection
   clear_selection();
 }
+
+inline void igl::MouseController::reset_selected_rotations()
+{
+  using namespace Eigen;
+  for(int e = 0;e<m_selection.size();e++)
+  {
+    if(m_selection(e))
+    {
+      m_rotations[e] = Quaterniond::Identity();
+    }
+  }
+}
+
 inline bool igl::MouseController::set_rotations(const RotationList & vQ)
 {
   if(vQ.size() != m_rotations.size())

include/igl/RotateWidget.h

@@ -238,7 +238,7 @@ inline bool igl::RotateWidget::down(const int x, const int y)
     auto on_meridian = [&](
       const Vector3d & hit, 
       const Quaterniond & rot, 
-      const Quaterniond m,
+      const Quaterniond & m,
       Vector3d & pl_hit) -> bool
     {
       // project onto rotate plane

include/igl/WindingNumberAABB.h

@@ -1,3 +1,15 @@
+// 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/.
+
+// # MUTUAL DEPENDENCY ISSUE FOR HEADER ONLY VERSION
+// MUST INCLUDE winding_number.h first before guard:
+#include "winding_number.h"
+
 #ifndef IGL_WINDINGNUMBERAABB_H
 #define IGL_WINDINGNUMBERAABB_H
 #include "WindingNumberTree.h"
@@ -202,7 +214,10 @@ inline bool igl::WindingNumberAABB<Point>::inside(const Point & p) const
   assert(p.size() == min_corner.size());
   for(int i = 0;i<p.size();i++)
   {
-    if( p(i) < min_corner(i) || p(i) >= max_corner(i))
+    //// Perfect matching is **not** robust
+    //if( p(i) < min_corner(i) || p(i) >= max_corner(i))
+    // **MUST** be conservative!!
+    if( p(i) < min_corner(i) || p(i) > max_corner(i))
     {
       return false;
     }

include/igl/WindingNumberTree.h

@@ -1,5 +1,12 @@
-#ifndef IGL_BOUNDINGTREE_H
-#define IGL_BOUNDINGTREE_H
+// 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_WINDINGNUMBERTREE_H
+#define IGL_WINDINGNUMBERTREE_H
 #include <list>
 #include <map>
 #include <Eigen/Dense>

include/igl/avg_edge_length.cpp

@@ -19,10 +19,10 @@ IGL_INLINE double igl::avg_edge_length(
 
   for (unsigned i=0;i<F.rows();++i)
   {
-    for (unsigned j=0;j<3;++j)
+    for (unsigned j=0;j<F.cols();++j)
     {
       ++count;
-      avg += (V.row(F(i,j)) - V.row(F(i,(j+1)%3))).norm();
+      avg += (V.row(F(i,j)) - V.row(F(i,(j+1)%F.cols()))).norm();
     }
   }
 

include/igl/avg_edge_length.h

@@ -22,7 +22,7 @@ namespace igl
   //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
   // Inputs:
   //   V  eigen matrix #V by 3
-  //   F  #F by 3 list of mesh faces (must be triangles)
+  //   F  #F by simplex-size list of mesh faces (must be simplex)
   // Outputs:
   //   l  average edge length
   //

include/igl/bone_parents.cpp

@@ -0,0 +1,25 @@
+#include "bone_parents.h"
+
+template <typename DerivedBE, typename DerivedP>
+IGL_INLINE void igl::bone_parents(
+  const Eigen::PlainObjectBase<DerivedBE>& BE,
+  Eigen::PlainObjectBase<DerivedP>& P)
+{
+  P.resize(BE.rows(),1);
+  // Stupid O(n²) version
+  for(int e = 0;e<BE.rows();e++)
+  {
+    P(e) = -1;
+    for(int f = 0;f<BE.rows();f++)
+    {
+      if(BE(e,0) == BE(f,1))
+      {
+        P(e) = f;
+      }
+    }
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::bone_parents<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&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/bone_parents.h

@@ -0,0 +1,32 @@
+// 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_BONE_PARENTS_H
+#define IGL_BONE_PARENTS_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // BONE_PARENTS Recover "parent" bones from directed graph representation.
+  // 
+  // Inputs:
+  //   BE  #BE by 2 list of directed bone edges
+  // Outputs:
+  //   P  #BE by 1 list of parent indices into BE, -1 means root.
+  //
+  template <typename DerivedBE, typename DerivedP>
+  IGL_INLINE void bone_parents(
+    const Eigen::PlainObjectBase<DerivedBE>& BE,
+    Eigen::PlainObjectBase<DerivedP>& P);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "bone_parents.cpp"
+#endif
+
+#endif
+

include/igl/boundary_conditions.cpp

@@ -157,6 +157,12 @@ IGL_INLINE bool igl::boundary_conditions(
     bc.row(i).array() /= sum;
   }
 
+  if(bc.size() == 0)
+  {
+    verbose("^%s: Error: boundary conditions are empty.\n",__FUNCTION__);
+    return false;
+  }
+
   // If there's only a single boundary condition, the following tests
   // are overzealous.
   if(bc.rows() == 1)

include/igl/boundary_conditions.h

@@ -13,7 +13,11 @@
 namespace igl
 {
 
-  // Compute boundary conditions for automatic weights computation
+  // Compute boundary conditions for automatic weights computation. This
+  // function expects that the given mesh (V,Ele) has sufficient samples
+  // (vertices) exactly at point handle locations and exactly along bone and
+  // cage edges.
+  //
   // Inputs:
   //   V  #V by dim list of domain vertices
   //   Ele  #Ele by simplex-size list of simplex indices
@@ -24,10 +28,10 @@ namespace igl
   // Outputs:
   //   b  #b list of boundary indices (indices into V of vertices which have
   //     known, fixed values)
-  //   bc #b by #weights list of known/fixed values for boundary vertices (notice
-  //     the #b != #weights in general because #b will include all the
-  //     intermediary samples along each bone, etc.. The ordering of the weights
-  //     corresponds to [P;BE]
+  //   bc #b by #weights list of known/fixed values for boundary vertices
+  //     (notice the #b != #weights in general because #b will include all the
+  //     intermediary samples along each bone, etc.. The ordering of the
+  //     weights corresponds to [P;BE]
   // Returns true if boundary conditions make sense
   IGL_INLINE bool boundary_conditions(
     const Eigen::MatrixXd & V,

include/igl/boundary_facets.cpp

@@ -133,9 +133,10 @@ Ret igl::boundary_facets(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::boundary_facets<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template void igl::boundary_facets<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&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 template void igl::boundary_facets<int, int>(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > >&);
 //template Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > igl::boundary_facets(Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&);
 template Eigen::Matrix<int, -1, -1, 0, -1, -1> igl::boundary_facets<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&);
 #endif
-

include/igl/centroid.cpp

@@ -0,0 +1,63 @@
+// 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 "centroid.h"
+#include <Eigen/Geometry>
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename Derivedc, 
+  typename Derivedvol>
+IGL_INLINE void igl::centroid(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<Derivedc>& cen,
+  Derivedvol & vol)
+{
+  using namespace Eigen;
+  assert(F.cols() == 3 && "F should contain triangles.");
+  assert(V.cols() == 3 && "V should contain 3d points.");
+  const int m = F.rows();
+  cen.setZero();
+  vol = 0;
+  // loop over faces
+  for(int f = 0;f<m;f++)
+  {
+    // "Calculating the volume and centroid of a polyhedron in 3d" [Nuernberg 2013]
+    // http://www2.imperial.ac.uk/~rn/centroid.pdf
+    // rename corners
+    const RowVector3d & a = V.row(F(f,0));
+    const RowVector3d & b = V.row(F(f,1));
+    const RowVector3d & c = V.row(F(f,2));
+    // un-normalized normal
+    const RowVector3d & n = (b-a).cross(c-a);
+    // total volume via divergence theorem: ∫ 1
+    vol += n.dot(a)/6.;
+    // centroid via divergence theorem and midpoint quadrature: ∫ x
+    cen.array() += (1./24.*n.array()*((a+b).array().square() + (b+c).array().square() + 
+        (c+a).array().square()).array());
+  }
+  cen *= 1./(2.*vol);
+}
+
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename Derivedc>
+IGL_INLINE void igl::centroid(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<Derivedc>& c)
+{
+  typename Derivedc::Scalar vol;
+  return centroid(V,F,c,vol);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template void igl::centroid<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, 3, 1, 0, 3, 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, 3, 1, 0, 3, 1> >&);
+#endif

include/igl/centroid.h

@@ -0,0 +1,49 @@
+// 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_CENTROID_H
+#define IGL_CENTROID_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // CENTROID Computes the centroid of a closed mesh using a surface integral.
+  // 
+  // Inputs:
+  //   V  #V by dim list of rest domain positions
+  //   F  #F by 3 list of triangle indices into V
+  // Outputs:
+  //    c  dim vector of centroid coordinates
+  //    vol  total volume of solid.
+  //
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename Derivedc, 
+    typename Derivedvol>
+  IGL_INLINE void centroid(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<Derivedc>& c,
+    Derivedvol & vol);
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename Derivedc>
+  IGL_INLINE void centroid(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    Eigen::PlainObjectBase<Derivedc>& c);
+
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "centroid.cpp"
+#endif
+
+#endif
+

include/igl/cgal/complex_to_mesh.cpp

@@ -0,0 +1,150 @@
+// 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 "complex_to_mesh.h"
+
+#include <igl/centroid.h>
+#include <igl/remove_unreferenced.h>
+
+#include <CGAL/Surface_mesh_default_triangulation_3.h>
+#include <set>
+#include <stack>
+
+template <typename Tr, typename DerivedV, typename DerivedF>
+IGL_INLINE bool igl::complex_to_mesh(
+  const CGAL::Complex_2_in_triangulation_3<Tr> & c2t3,
+  Eigen::PlainObjectBase<DerivedV> & V, 
+  Eigen::PlainObjectBase<DerivedF> & F)
+{
+  using namespace igl;
+  using namespace Eigen;
+  // CGAL/IO/Complex_2_in_triangulation_3_file_writer.h
+  using CGAL::Surface_mesher::number_of_facets_on_surface;
+
+  typedef typename CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
+  typedef typename Tr::Finite_facets_iterator Finite_facets_iterator;
+  typedef typename Tr::Finite_vertices_iterator Finite_vertices_iterator;
+  typedef typename Tr::Facet Facet;
+  typedef typename Tr::Edge Edge;
+  typedef typename Tr::Vertex_handle Vertex_handle;
+
+  // Header.
+  const Tr& tr = c2t3.triangulation();
+
+  bool success = true;
+
+  const int n = tr.number_of_vertices();
+  const int m = c2t3.number_of_facets();
+
+  assert(m == number_of_facets_on_surface(tr));
+  
+  // Finite vertices coordinates.
+  std::map<Vertex_handle, int> v2i;
+  V.resize(n,3);
+  {
+    int v = 0;
+    for(Finite_vertices_iterator vit = tr.finite_vertices_begin();
+        vit != tr.finite_vertices_end();
+        ++vit)
+    {
+      V(v,0) = vit->point().x(); 
+      V(v,1) = vit->point().y(); 
+      V(v,2) = vit->point().z(); 
+      v2i[vit] = v++;
+    }
+  }
+
+  {
+    Finite_facets_iterator fit = tr.finite_facets_begin();
+    std::set<Facet> oriented_set;
+    std::stack<Facet> stack;
+
+    while ((int)oriented_set.size() != m) 
+    {
+      while ( fit->first->is_facet_on_surface(fit->second) == false ||
+        oriented_set.find(*fit) != oriented_set.end() ||
+        oriented_set.find(c2t3.opposite_facet(*fit)) !=
+        oriented_set.end() ) 
+      {
+        ++fit;
+      }
+      oriented_set.insert(*fit);
+      stack.push(*fit);
+      while(! stack.empty() )
+      {
+        Facet f = stack.top();
+        stack.pop();
+        for(int ih = 0 ; ih < 3 ; ++ih)
+        {
+          const int i1  = tr.vertex_triple_index(f.second, tr. cw(ih));
+          const int i2  = tr.vertex_triple_index(f.second, tr.ccw(ih));
+
+          const typename C2t3::Face_status face_status
+            = c2t3.face_status(Edge(f.first, i1, i2));
+          if(face_status == C2t3::REGULAR) 
+          {
+            Facet fn = c2t3.neighbor(f, ih);
+            if (oriented_set.find(fn) == oriented_set.end())
+            {
+              if(oriented_set.find(c2t3.opposite_facet(fn)) == oriented_set.end())
+              {
+                oriented_set.insert(fn);
+                stack.push(fn);
+              }else {
+                success = false; // non-orientable
+              }
+            }
+          }else if(face_status != C2t3::BOUNDARY)
+          {
+            success = false; // non manifold, thus non-orientable
+          }
+        } // end "for each neighbor of f"
+      } // end "stack non empty"
+    } // end "oriented_set not full"
+
+    F.resize(m,3);
+    int f = 0;
+    for(typename std::set<Facet>::const_iterator fit = 
+        oriented_set.begin();
+        fit != oriented_set.end();
+        ++fit)
+    {
+      const typename Tr::Cell_handle cell = fit->first;
+      const int& index = fit->second;
+      const int index1 = v2i[cell->vertex(tr.vertex_triple_index(index, 0))];
+      const int index2 = v2i[cell->vertex(tr.vertex_triple_index(index, 1))];
+      const int index3 = v2i[cell->vertex(tr.vertex_triple_index(index, 2))];
+      // This order is flipped
+      F(f,0) = index1;
+      F(f,1) = index2;
+      F(f,2) = index3;
+      f++;
+    }
+    assert(f == m);
+  } // end if(facets must be oriented)
+
+  // This CGAL code seems to randomly assign the global orientation
+  // Flip based on the signed volume.
+  Eigen::Vector3d cen;
+  double vol;
+  igl::centroid(V,F,cen,vol);
+  if(vol < 0)
+  {
+    // Flip
+    F = F.rowwise().reverse().eval();
+  }
+
+  // CGAL code somehow can end up with unreferenced vertices
+  VectorXi _1;
+  remove_unreferenced( MatrixXd(V), MatrixXi(F), V,F,_1);
+
+  return success;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template bool igl::complex_to_mesh<CGAL::Delaunay_triangulation_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_data_structure_3<CGAL::Surface_mesh_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Triangulation_cell_base_with_circumcenter_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Surface_mesh_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_cell_base_3<void> > > > >, CGAL::Default>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(CGAL::Complex_2_in_triangulation_3<CGAL::Delaunay_triangulation_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_data_structure_3<CGAL::Surface_mesh_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_vertex_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_vertex_base_3<void> > >, CGAL::Triangulation_cell_base_with_circumcenter_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Surface_mesh_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_cell_base_3<CGAL::Robust_circumcenter_traits_3<CGAL::Epick>, CGAL::Triangulation_ds_cell_base_3<void> > > > >, CGAL::Default>, void> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+#endif

include/igl/cgal/complex_to_mesh.h

@@ -0,0 +1,41 @@
+// 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_COMPLEX_TO_MESH_H
+#define IGL_COMPLEX_TO_MESH_H
+#include "../igl_inline.h"
+
+#include <Eigen/Dense>
+#include <CGAL/Complex_2_in_triangulation_3.h>
+
+namespace igl 
+{
+  // Templates:
+  //   Tr  CGAL triangulation type, e.g.
+  //     CGAL::Surface_mesh_default_triangulation_3
+  // Inputs
+  //   c2t3  2-complex (surface) living in a 3d triangulation (e.g. result of
+  //     CGAL::make_surface_mesh)
+  // Outputs:
+  //   V  #V by 3 list of vertex positions
+  //   F  #F by 3 list of triangle indices
+  // Returns true iff conversion was successful, failure can ok if CGAL code
+  // can't figure out ordering.
+  //
+  template <typename Tr, typename DerivedV, typename DerivedF>
+  IGL_INLINE bool complex_to_mesh(
+    const CGAL::Complex_2_in_triangulation_3<Tr> & c2t3,
+    Eigen::PlainObjectBase<DerivedV> & V, 
+    Eigen::PlainObjectBase<DerivedF> & F);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "complex_to_mesh.cpp"
+#endif
+
+#endif
+

include/igl/cgal/point_mesh_squared_distance.cpp

@@ -18,26 +18,72 @@ IGL_INLINE void igl::point_mesh_squared_distance(
   Eigen::MatrixXd & C)
 {
   using namespace std;
-  typedef CGAL::Point_3<Kernel>    Point_3;
   typedef CGAL::Triangle_3<Kernel> Triangle_3; 
   typedef typename std::vector<Triangle_3>::iterator Iterator;
   typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
   typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
   typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
-  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+  Tree tree;
+  vector<Triangle_3> T;
+  point_mesh_squared_distance_precompute(V,F,tree,T);
+  return point_mesh_squared_distance(P,tree,T,sqrD,I,C);
+}
+
+template <typename Kernel>
+IGL_INLINE void igl::point_mesh_squared_distance_precompute(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  std::vector<CGAL::Triangle_3<Kernel> > & T)
+{
+  using namespace std;
+
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
 
   // Must be 3D
   assert(V.cols() == 3);
-  assert(P.cols() == 3);
   // Must be triangles
   assert(F.cols() == 3);
   // Make list of cgal triangles
-  Tree tree;
-  vector<Triangle_3> T;
   mesh_to_cgal_triangle_list(V,F,T);
+  tree.clear();
   tree.insert(T.begin(),T.end());
-
   tree.accelerate_distance_queries();
+}
+
+template <typename Kernel>
+IGL_INLINE void igl::point_mesh_squared_distance(
+  const Eigen::MatrixXd & P,
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const std::vector<CGAL::Triangle_3<Kernel> > & T,
+  Eigen::VectorXd & sqrD,
+  Eigen::VectorXi & I,
+  Eigen::MatrixXd & C)
+{
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+  typedef CGAL::Point_3<Kernel>    Point_3;
+  assert(P.cols() == 3);
   const int n = P.rows();
   sqrD.resize(n,1);
   I.resize(n,1);
@@ -59,5 +105,6 @@ IGL_INLINE void igl::point_mesh_squared_distance(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::point_mesh_squared_distance<CGAL::Epeck>( const Eigen::MatrixXd & P, const Eigen::MatrixXd & V, const Eigen::MatrixXi & F, Eigen::VectorXd & sqrD, Eigen::VectorXi & I, Eigen::MatrixXd & C);
+template void igl::point_mesh_squared_distance_precompute<CGAL::Epick>(Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > >&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >&);
 
 #endif

include/igl/cgal/point_mesh_squared_distance.h

@@ -9,6 +9,7 @@
 #define IGL_POINT_MESH_SQUARED_DISTANCE_H
 #include <igl/igl_inline.h>
 #include <Eigen/Core>
+#include <vector>
 #include "CGAL_includes.hpp"
 namespace igl
 {
@@ -36,7 +37,45 @@ namespace igl
     Eigen::VectorXd & sqrD,
     Eigen::VectorXi & I,
     Eigen::MatrixXd & C);
+  // Probably can do this in a way that we don't pass around `tree` and `T`
+  //
+  // Outputs:
+  //   tree  CGAL's AABB tree
+  //   T  list of CGAL triangles in order of F (for determining which was found
+  //     in computation)
+  template <typename Kernel>
+  IGL_INLINE void point_mesh_squared_distance_precompute(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    std::vector<CGAL::Triangle_3<Kernel> > & T);
+  // Inputs:
+  //  see above
+  // Outputs:
+  //  see above
+  template <typename Kernel>
+  IGL_INLINE void point_mesh_squared_distance(
+    const Eigen::MatrixXd & P,
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const std::vector<CGAL::Triangle_3<Kernel> > & T,
+    Eigen::VectorXd & sqrD,
+    Eigen::VectorXi & I,
+    Eigen::MatrixXd & C);
+
 }
+
 #ifndef IGL_STATIC_LIBRARY
 #  include "point_mesh_squared_distance.cpp"
 #endif

include/igl/cgal/signed_distance.cpp

@@ -0,0 +1,125 @@
+// 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 "signed_distance.h"
+template <typename Kernel>
+IGL_INLINE typename Kernel::FT igl::signed_distance_pseudonormal(
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const std::vector<CGAL::Triangle_3<Kernel> > & T,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & FN,
+  const Eigen::MatrixXd & VN,
+  const Eigen::MatrixXd & EN,
+  const Eigen::VectorXi & EMAP,
+  const typename Kernel::Point_3 & q)
+{
+  using namespace Eigen;
+  using namespace std;
+  typedef typename Kernel::FT FT;
+  typedef typename Kernel::Point_3 Point_3;
+  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+
+  Point_and_primitive_id pp = tree.closest_point_and_primitive(q);
+  Point_3 & p = pp.first;
+  const auto & qp = q-p;
+  const FT sqrd = qp.squared_length();
+  Vector3d v(qp.x(),qp.y(),qp.z());
+  const int f = pp.second - T.begin();
+  const Triangle_3 & t = *pp.second;
+  // barycentric coordinates
+  const auto & area = [&p,&t](const int i, const int j)->FT
+  {
+    return sqrt(Triangle_3(p,t.vertex(i),t.vertex(j)).squared_area());
+  };
+  Vector3d b(area(1,2),area(2,0),area(0,1));
+  b /= b.sum();
+  // Determine which normal to use
+  Vector3d n;
+  const double epsilon = 1e-12;
+  const int type = (b.array()<=epsilon).cast<int>().sum();
+  switch(type)
+  {
+    case 2:
+      // Find vertex
+      for(int c = 0;c<3;c++)
+      {
+        if(b(c)>epsilon)
+        {
+          n = VN.row(F(f,c));
+          break;
+        }
+      }
+      break;
+    case 1:
+      // Find edge
+      for(int c = 0;c<3;c++)
+      {
+        if(b(c)<=epsilon)
+        {
+          n = EN.row(EMAP(F.rows()*c+f));
+          break;
+        }
+      }
+      break;
+    default:
+      assert(false && "all barycentric coords zero.");
+    case 0:
+      n = FN.row(f);
+      break;
+  }
+  const double s = (v.dot(n) >= 0 ? 1. : -1.);
+  return s*sqrt(sqrd);
+}
+
+template <typename Kernel>
+IGL_INLINE typename Kernel::FT igl::signed_distance_winding_number(
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & tree,
+  const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+  const typename Kernel::Point_3 & q)
+{
+  typedef typename Kernel::FT FT;
+  typedef typename Kernel::Point_3 Point_3;
+  typedef typename CGAL::Triangle_3<Kernel> Triangle_3;
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef typename CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef typename CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef typename CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+  using namespace Eigen;
+  using namespace std;
+  using namespace igl;
+  Point_and_primitive_id pp = tree.closest_point_and_primitive(q);
+  Point_3 & p = pp.first;
+  const auto & qp = q-p;
+  const Vector3d eq(q.x(),q.y(),q.z()); 
+  const FT sqrd = qp.squared_length();
+  const double w = hier.winding_number(eq);
+  const FT s = 1.-2.*w;
+  return s*sqrt(sqrd);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+template CGAL::Epick::FT igl::signed_distance_winding_number<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, igl::WindingNumberAABB<Eigen::Matrix<double, 3, 1, 0, 3, 1> > const&, CGAL::Epick::Point_3 const&);
+template CGAL::Epick::FT igl::signed_distance_pseudonormal<CGAL::Epick>(CGAL::AABB_tree<CGAL::AABB_traits<CGAL::Epick, CGAL::AABB_triangle_primitive<CGAL::Epick, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > >::iterator, CGAL::Boolean_tag<false> > > > const&, std::vector<CGAL::Triangle_3<CGAL::Epick>, std::allocator<CGAL::Triangle_3<CGAL::Epick> > > const&, Eigen::Matrix<int, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1> const&, Eigen::Matrix<int, -1, 1, 0, -1, 1> const&, CGAL::Epick::Point_3 const&);
+#endif

include/igl/cgal/signed_distance.h

@@ -0,0 +1,79 @@
+// 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_SIGNED_DISTANCE_H
+#define IGL_SIGNED_DISTANCE_H
+#include <igl/igl_inline.h>
+#include <igl/WindingNumberAABB.h>
+#include <Eigen/Core>
+#include <vector>
+#include "CGAL_includes.hpp"
+namespace igl
+{
+  enum SignedDistanceType
+  {
+    SIGNED_DISTANCE_TYPE_PSEUDONORMAL   = 0,
+    SIGNED_DISTANCE_TYPE_WINDING_NUMBER = 1,
+    SIGNED_DISTANCE_TYPE_DEFAULT        = 2,
+    NUM_SIGNED_DISTANCE_TYPE            = 3
+  };
+  // Computes signed distance to mesh
+  //
+  // Templates:
+  //   Kernal  CGAL computation and construction kernel (e.g.
+  //     CGAL::Simple_cartesian<double>)
+  // Inputs:
+  //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
+  //   T  #F list of CGAL triangles (see point_mesh_squared_distance.h)
+  //   F  #F by 3 list of triangle indices
+  //   FN  #F by 3 list of triangle normals 
+  //   VN  #V by 3 list of vertex normals (ANGLE WEIGHTING)
+  //   EN  #E by 3 list of edge normals (UNIFORM WEIGHTING)
+  //   EMAP  #F*3 mapping edges in F to E
+  //   q  Query point
+  // Returns signed distance to mesh
+  //
+  template <typename Kernel>
+  IGL_INLINE typename Kernel::FT signed_distance_pseudonormal(
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const std::vector<CGAL::Triangle_3<Kernel> > & T,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXd & FN,
+    const Eigen::MatrixXd & VN,
+    const Eigen::MatrixXd & EN,
+    const Eigen::VectorXi & EMAP,
+    const typename Kernel::Point_3 & q);
+  // Inputs:
+  //   tree  AABB acceleration tree (see point_mesh_squared_distance.h)
+  //   hier  Winding number evaluation hierarchy
+  //   q  Query point
+  // Returns signed distance to mesh
+  template <typename Kernel>
+  IGL_INLINE typename Kernel::FT signed_distance_winding_number(
+    const CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    > & tree,
+    const igl::WindingNumberAABB<Eigen::Vector3d> & hier,
+    const typename Kernel::Point_3 & q);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "signed_distance.cpp"
+#endif
+
+#endif
+

include/igl/cgal/signed_distance_isosurface.cpp

@@ -0,0 +1,136 @@
+// 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 "signed_distance_isosurface.h"
+#include "point_mesh_squared_distance.h"
+#include "complex_to_mesh.h"
+#include "signed_distance.h"
+
+#include <igl/per_face_normals.h>
+#include <igl/per_edge_normals.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/centroid.h>
+#include <igl/WindingNumberAABB.h>
+#include <igl/matlab_format.h>
+#include <igl/remove_unreferenced.h>
+
+#include <CGAL/Surface_mesh_default_triangulation_3.h>
+#include <CGAL/Complex_2_in_triangulation_3.h>
+#include <CGAL/make_surface_mesh.h>
+#include <CGAL/Implicit_surface_3.h>
+#include <CGAL/Polyhedron_3.h>
+#include <CGAL/IO/output_surface_facets_to_polyhedron.h>
+// 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>
+#include <vector>
+
+IGL_INLINE bool igl::signed_distance_isosurface(
+  const Eigen::MatrixXd & IV,
+  const Eigen::MatrixXi & IF,
+  const double level,
+  const double angle_bound,
+  const double radius_bound,
+  const double distance_bound,
+  const SignedDistanceType sign_type,
+  Eigen::MatrixXd & V,
+  Eigen::MatrixXi & F)
+{
+  using namespace std;
+
+  // default triangulation for Surface_mesher
+  typedef CGAL::Surface_mesh_default_triangulation_3 Tr;
+  // c2t3
+  typedef CGAL::Complex_2_in_triangulation_3<Tr> C2t3;
+  typedef Tr::Geom_traits GT;//Kernel
+  typedef GT::Sphere_3 Sphere_3;
+  typedef GT::Point_3 Point_3;
+  typedef GT::FT FT;
+  typedef std::function<FT (Point_3)> Function;
+  typedef CGAL::Implicit_surface_3<GT, Function> Surface_3;
+  typedef CGAL::Polyhedron_3<GT> Polyhedron;
+  typedef GT::Kernel Kernel;
+  typedef CGAL::Triangle_3<Kernel> Triangle_3; 
+  typedef typename std::vector<Triangle_3>::iterator Iterator;
+  typedef CGAL::AABB_triangle_primitive<Kernel, Iterator> Primitive;
+  typedef CGAL::AABB_traits<Kernel, Primitive> AABB_triangle_traits;
+  typedef CGAL::AABB_tree<AABB_triangle_traits> Tree;
+  typedef typename Tree::Point_and_primitive_id Point_and_primitive_id;
+
+  Eigen::MatrixXd FN,VN,EN;
+  Eigen::MatrixXi E;
+  Eigen::VectorXi EMAP;
+  // "Signed Distance Computation Using the Angle Weighted Pseudonormal"
+  // [Bærentzen & Aanæs 2005]
+  per_face_normals(IV,IF,FN);
+  per_vertex_normals(IV,IF,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,VN);
+  per_edge_normals(IV,IF,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,EN,E,EMAP);
+  // Prepare distance computation
+  Tree tree;
+  vector<Triangle_3 > T;
+  point_mesh_squared_distance_precompute(IV,IF,tree,T);
+  WindingNumberAABB<Eigen::Vector3d> hier(IV,IF);
+  hier.grow();
+  Tr tr;            // 3D-Delaunay triangulation
+  C2t3 c2t3 (tr);   // 2D-complex in 3D-Delaunay triangulation
+  // defining the surface
+  const auto & IVmax = IV.colwise().maxCoeff();
+  const auto & IVmin = IV.colwise().minCoeff();
+  const double bbd = (IVmax-IVmin).norm();
+  const double r = bbd/2.;
+  const auto & IVmid = 0.5*(IVmax + IVmin);
+  // Supposedly the implict needs to evaluate to <0 at cmid...
+  // http://doc.cgal.org/latest/Surface_mesher/classCGAL_1_1Implicit__surface__3.html
+  Point_3 cmid(IVmid(0),IVmid(1),IVmid(2));
+  Function fun;
+  switch(sign_type)
+  {
+    default:
+      assert(false && "Unknown SignedDistanceType");
+    case SIGNED_DISTANCE_TYPE_DEFAULT:
+    case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
+      fun = 
+        [&tree,&hier,&level](const Point_3 & q) -> FT
+        {
+          return signed_distance_winding_number(tree,hier,q)-level;
+        };
+      break;
+    case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
+      fun = [&tree,&T,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) -> FT
+        {
+          return 
+            igl::signed_distance_pseudonormal(tree,T,IF,FN,VN,EN,EMAP,q) - 
+            level;
+        };
+      break;
+  }
+    //[&tree,&hier,&T,&IF,&FN,&VN,&EN,&EMAP,&level](const Point_3 & q) ->FT
+    //{
+    //  const FT p = signed_distance_pseudonormal(tree,T,IF,FN,VN,EN,EMAP,q);
+    //  const FT w = signed_distance_winding_number(tree,hier,q);
+    //  if(w*p < 0 && (fabs(w) > 0.1 || fabs(p) > 0.1))
+    //  {
+    //    cout<<"q=["<<q.x()<<","<<q.y()<<","<<q.z()<<"];"<<endl;
+    //    cout<<matlab_format(n.transpose().eval(),"n")<<endl;
+    //    cout<<matlab_format(b.transpose().eval(),"b")<<endl;
+    //    cout<<"Sig difference: "<<type<<endl;
+    //    cout<<"w: "<<w<<endl;
+    //    cout<<"p: "<<p<<endl;
+    //    exit(1);
+    //  }
+    //  return w;
+    //},
+  Sphere_3 bounding_sphere(cmid, (r+level)*(r+level));
+  Surface_3 surface(fun,bounding_sphere);
+  CGAL::Surface_mesh_default_criteria_3<Tr> 
+    criteria(angle_bound,radius_bound,distance_bound);
+  // meshing surface
+  CGAL::make_surface_mesh(c2t3, surface, criteria, CGAL::Manifold_tag());
+  // complex to (V,F)
+  return igl::complex_to_mesh(c2t3,V,F);
+}

include/igl/cgal/signed_distance_isosurface.h

@@ -0,0 +1,47 @@
+// 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_SIGNED_DISTANCE_ISOSURFACE_H
+#define IGL_SIGNED_DISTANCE_ISOSURFACE_H
+#include <igl/igl_inline.h>
+#include "signed_distance.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // SIGNED_DISTANCE_ISOSURFACE Compute the contour of an iso-level of the
+  // signed distance field to a given mesh.
+  //
+  // Inputs:
+  //   IV  #IV by 3 list of input mesh vertex positions
+  //   IF  #IF by 3 list of input triangle indices
+  //   level  iso-level to contour (e.g. 0)
+  //   angle_bound  lower bound on triangle angles (mesh quality) (e.g. 28)
+  //   radius_bound  upper bound on triangle size (mesh density?) (e.g. 0.02)
+  //   distance_bound  cgal mysterious parameter (mesh density?) (e.g. 0.01)
+  //   sign_type  method for computing distance _sign_ (see signed_distance.h)
+  // Outputs:
+  //   V  #V by 3 list of input mesh vertex positions
+  //   F  #F by 3 list of input triangle indices
+  //  
+  IGL_INLINE bool signed_distance_isosurface(
+    const Eigen::MatrixXd & IV,
+    const Eigen::MatrixXi & IF,
+    const double level,
+    const double angle_bound,
+    const double radius_bound,
+    const double distance_bound,
+    const SignedDistanceType sign_type,
+    Eigen::MatrixXd & V,
+    Eigen::MatrixXi & F);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "signed_distance_isosurface.cpp"
+#endif
+
+#endif
+

include/igl/deform_skeleton.cpp

@@ -39,8 +39,10 @@ IGL_INLINE void igl::deform_skeleton(
   {
     BET(e,0) = 2*e;
     BET(e,1) = 2*e+1;
+    Matrix4d t;
+    t << T.block(e*4,0,4,3).transpose(), 0,0,0,0;
     Affine3d a;
-    a.matrix() = T.block(e*4,0,4,3).transpose();
+    a.matrix() = t;
     Vector3d c0 = C.row(BE(e,0));
     Vector3d c1 = C.row(BE(e,1));
     CT.row(2*e) =   a * c0;

include/igl/doublearea.cpp

@@ -139,6 +139,7 @@ IGL_INLINE void igl::doublearea(
 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, 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<double, -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> >&);
 template void igl::doublearea<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<int, -1, 1, 0, -1, 1> >&);

include/igl/draw_mesh.cpp

@@ -75,6 +75,7 @@ IGL_INLINE void igl::draw_mesh(
   using namespace std;
   using namespace Eigen;
   const int rF = F.rows();
+  const int cF = F.cols();
   const int cC = C.cols();
   const int rC = C.rows();
   const int cW = W.cols();
@@ -93,6 +94,9 @@ IGL_INLINE void igl::draw_mesh(
     assert(C.cols() == 3 || C.size() == 0);
     assert(N.cols() == 3 || N.size() == 0);
     assert(TC.cols() == 2 || TC.size() == 0);
+    assert(cF == 3 || cF == 4);
+    assert(TF.size() == 0 || TF.cols() == F.cols());
+    assert(NF.size() == 0 || NF.cols() == NF.cols());
   }
   if(W.size()>0)
   {
@@ -101,12 +105,21 @@ IGL_INLINE void igl::draw_mesh(
     assert(W.cols() == WI.cols());
   }
 
-  glBegin(GL_TRIANGLES);
+  switch(F.cols())
+  {
+    default:
+    case 3:
+      glBegin(GL_TRIANGLES);
+      break;
+    case 4:
+      glBegin(GL_QUADS);
+      break;
+  }
   // loop over faces
   for(int i = 0; i<rF;i++)
   {
     // loop over corners of triangle
-    for(int j = 0;j<3;j++)
+    for(int j = 0;j<cF;j++)
     {
 
       int tc = -1;
@@ -119,9 +132,9 @@ IGL_INLINE void igl::draw_mesh(
       }else if(rTC == rV)
       {
         tc = F(i,j);
-      }else if(rTC == rF*2)
+      }else if(rTC == rF*cF)
       {
-        tc = i*2 + j;
+        tc = i*cF + j;
       }else if(rTC == rF)
       {
         tc = i;
@@ -138,9 +151,9 @@ IGL_INLINE void igl::draw_mesh(
       }else if(rC == rV)
       {
         color = C.row(F(i,j));
-      }else if(rC == rF*3)
+      }else if(rC == rF*cF)
       {
-        color = C.row(i*3+j);
+        color = C.row(i*cF+j);
       }else if(rC == rF)
       {
         color = C.row(i);
@@ -152,19 +165,19 @@ IGL_INLINE void igl::draw_mesh(
       int n = -1;
       if(rNF != 0)
       {
-        n = NF(i,j);
-      } else if(rN == 1)
+        n = NF(i,j); // indexed normals
+      } else if(rN == 1) 
       {
-        n = 0;
-      }else if(rN == rV)
+        n = 0; // uniform normals
+      }else if(rN == rF)
       {
-        n = F(i,j);
-      }else if(rN == rF*2)
+        n = i; // face normals
+      }else if(rN == rV)
       {
-        n = i*2 + j;
-      }else if(rN == rF)
+        n = F(i,j); // vertex normals
+      }else if(rN == rF*cF)
       {
-        n = i;
+        n = i*cF + j; // corner normals
       }else
       {
         assert(N.size() == 0);

include/igl/draw_skeleton_3d.cpp

@@ -27,6 +27,10 @@ IGL_INLINE void igl::draw_skeleton_3d(
   // parameter. Further, its joint balls are not rotated with the bones.
   using namespace Eigen;
   using namespace std;
+  if(color.size() == 0)
+  {
+    return draw_skeleton_3d(C,BE,T,MAYA_SEA_GREEN);
+  }
   assert(color.cols() == 4 || color.size() == 4);
   if(T.size() == 0)
   {
@@ -38,7 +42,7 @@ IGL_INLINE void igl::draw_skeleton_3d(
     {
       return;
     }
-    return draw_skeleton_3d(C,BE,T,MAYA_SEA_GREEN);
+    return draw_skeleton_3d(C,BE,T,color);
   }
   assert(T.rows() == BE.rows()*(C.cols()+1));
   assert(T.cols() == C.cols());

include/igl/draw_skeleton_3d.h

@@ -29,19 +29,18 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedBE> & BE,
     const Eigen::PlainObjectBase<DerivedT> & T,
     const Eigen::PlainObjectBase<Derivedcolor> & color);
-  // Wrapper with each T = Identity
+  // Default color
   template <typename DerivedC, typename DerivedBE, typename DerivedT>
   IGL_INLINE void draw_skeleton_3d(
     const Eigen::PlainObjectBase<DerivedC> & C,
     const Eigen::PlainObjectBase<DerivedBE> & BE,
     const Eigen::PlainObjectBase<DerivedT> & T);
-  // Default color
   template <typename DerivedC, typename DerivedBE>
   IGL_INLINE void draw_skeleton_3d(
     const Eigen::PlainObjectBase<DerivedC> & C,
     const Eigen::PlainObjectBase<DerivedBE> & BE);
 };
 #ifndef IGL_STATIC_LIBRARY
-#  include "draw_skeleton_3d.h"
+#  include "draw_skeleton_3d.cpp"
 #endif
 #endif

include/igl/draw_skeleton_vector_graphics.h

@@ -44,6 +44,6 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedT> & T);
 }
 #ifndef IGL_STATIC_LIBRARY
-#  include "draw_skeleton_vector_graphics.h"
+#  include "draw_skeleton_vector_graphics.cpp"
 #endif
 #endif 

include/igl/edge_lengths.cpp

@@ -17,6 +17,15 @@ IGL_INLINE void igl::edge_lengths(
   using namespace std;
   switch(F.cols())
   {
+    case 2:
+    {
+      L.resize(F.rows(),1);
+      for(int i = 0;i<F.rows();i++)
+      {
+        L(i,0) = (V.row(F(i,1))-V.row(F(i,0))).norm();
+      }
+      break;
+    }
     case 3:
     {
       L.resize(F.rows(),3);
@@ -72,4 +81,5 @@ template void igl::edge_lengths<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen:
 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> >&);
 template void igl::edge_lengths<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 4, 0, -1, 4>, Eigen::Matrix<double, -1, 6, 0, -1, 6> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 4, 0, -1, 4> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 6, 0, -1, 6> >&);
+template void igl::edge_lengths<Eigen::Matrix<double, -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<double, -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> >&);
 #endif

include/igl/edge_lengths.h

@@ -22,11 +22,14 @@ namespace igl
   //   DerivedL derived from edge lengths matrix type: i.e. MatrixXd
   // Inputs:
   //   V  eigen matrix #V by 3
+  //   F  #F by 2 list of mesh edges
+  //    or
   //   F  #F by 3 list of mesh faces (must be triangles)
   //    or
   //   T  #T by 4 list of mesh elements (must be tets)
   // Outputs:
-  //   L  #F by {3|6} list of edge lengths 
+  //   L  #F by {1|3|6} list of edge lengths 
+  //     for edges, column of 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]

include/igl/embree/unproject_in_mesh.cpp

@@ -14,8 +14,8 @@
 template <
   typename Derivedobj>
 IGL_INLINE int igl::unproject_in_mesh(
-  const int x,
-  const int y,
+  const double x,
+  const double y,
   const igl::EmbreeIntersector & ei,
   Eigen::PlainObjectBase<Derivedobj> & obj)
 {
@@ -26,8 +26,8 @@ IGL_INLINE int igl::unproject_in_mesh(
 template <
   typename Derivedobj>
 IGL_INLINE int igl::unproject_in_mesh(
-  const int x,
-  const int y,
+  const double x,
+  const double y,
   const igl::EmbreeIntersector & ei,
   Eigen::PlainObjectBase<Derivedobj> & obj,
   std::vector<igl::Hit > & hits)
@@ -115,8 +115,8 @@ IGL_INLINE int igl::unproject_in_mesh(
 
 #ifdef IGL_STATIC_LIBRARY
 #  ifndef IGL_OPENLGL_4
-template int igl::unproject_in_mesh<Eigen::Matrix<double, 3, 1, 0, 3, 1> >(int, int, igl::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
-template int igl::unproject_in_mesh<Eigen::Matrix<double, 3, 1, 0, 3, 1> >(int, int, igl::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
+template int igl::unproject_in_mesh<Eigen::Matrix<double, 3, 1, 0, 3, 1> >(double, double, igl::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
+template int igl::unproject_in_mesh<Eigen::Matrix<double, 3, 1, 0, 3, 1> >(double, double, igl::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 3, 1, 0, 3, 1> >&);
 #  endif
 template int igl::unproject_in_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::Matrix<float, 2, 1, 0, 2, 1> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 4, 0, 4, 4> const&, Eigen::Matrix<float, 4, 1, 0, 4, 1> const&, igl::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, std::vector<igl::Hit, std::allocator<igl::Hit> >&);
 #endif

include/igl/embree/unproject_in_mesh.h

@@ -33,16 +33,16 @@ namespace igl
   template <
     typename Derivedobj>
   IGL_INLINE int unproject_in_mesh(
-    const int x,
-    const int y,
+    const double x,
+    const double y,
     const igl::EmbreeIntersector & ei,
     Eigen::PlainObjectBase<Derivedobj> & obj);
 
   template <
     typename Derivedobj>
   IGL_INLINE int unproject_in_mesh(
-    const int x,
-    const int y,
+    const double x,
+    const double y,
     const igl::EmbreeIntersector & ei,
     Eigen::PlainObjectBase<Derivedobj> & obj,
     std::vector<igl::Hit > & hits);

include/igl/exterior_edges.h

@@ -20,7 +20,7 @@ namespace igl
   Eigen::MatrixXi exterior_edges( const Eigen::MatrixXi & F);
 }
 #ifndef IGL_STATIC_LIBRARY
-#  include "exterior_edges.h"
+#  include "exterior_edges.cpp"
 #endif
 
 #endif

include/igl/forward_kinematics.cpp

@@ -54,6 +54,31 @@ IGL_INLINE void igl::forward_kinematics(
   }
 }
 
+IGL_INLINE void igl::forward_kinematics(
+  const Eigen::MatrixXd & C,
+  const Eigen::MatrixXi & BE,
+  const Eigen::VectorXi & P,
+  const std::vector<
+    Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> > & dQ,
+  Eigen::MatrixXd & T)
+{
+  using namespace Eigen;
+  using namespace std;
+  vector< Quaterniond,aligned_allocator<Quaterniond> > vQ;
+  vector< Vector3d> vT;
+  forward_kinematics(C,BE,P,dQ,vQ,vT);
+  const int dim = C.cols();
+  T.resize(BE.rows()*(dim+1),dim);
+  for(int e = 0;e<BE.rows();e++)
+  {
+    Affine3d a = Affine3d::Identity();
+    a.translate(vT[e]);
+    a.rotate(vQ[e]);
+    T.block(e*(dim+1),0,dim+1,dim) =
+      a.matrix().transpose().block(0,0,dim+1,dim);
+  }
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation
 #endif

include/igl/forward_kinematics.h

@@ -35,6 +35,15 @@ namespace igl
     std::vector<
       Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> > & vQ,
     std::vector<Eigen::Vector3d> & vT);
+  // Outputs:
+  //   T  #BE*(dim+1) by dim stack of transposed transformation matrices
+  IGL_INLINE void forward_kinematics(
+    const Eigen::MatrixXd & C,
+    const Eigen::MatrixXi & BE,
+    const Eigen::VectorXi & P,
+    const std::vector<
+      Eigen::Quaterniond,Eigen::aligned_allocator<Eigen::Quaterniond> > & dQ,
+    Eigen::MatrixXd & T);
 };
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/internal_angles.cpp

@@ -45,4 +45,5 @@ IGL_INLINE void igl::internal_angles(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 template void igl::internal_angles<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::internal_angles<Eigen::Matrix<double, -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<double, -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> >&);
 #endif

include/igl/list_to_matrix.cpp

@@ -53,6 +53,35 @@ IGL_INLINE bool igl::list_to_matrix(const std::vector<std::vector<T > > & V,Mat
   return true;
 }
 
+template <typename T, class Mat>
+IGL_INLINE bool igl::list_to_matrix(
+  const std::vector<std::vector<T > > & V,
+  const int n,
+  const T & padding,
+  Mat & M)
+{
+  const int m = V.size();
+  M.resize(m,n);
+  for(int i = 0;i<m;i++)
+  {
+    const auto & row = V[i];
+    if(row.size()>n)
+    {
+      return false;
+    }
+    int j = 0;
+    for(;j<row.size();j++)
+    {
+      M(i,j) = row[j];
+    }
+    for(;j<n;j++)
+    {
+      M(i,j) = padding;
+    }
+  }
+  return true;
+}
+
 template <typename T, class Mat>
 IGL_INLINE bool igl::list_to_matrix(const std::vector<T > & V,Mat & M)
 {

include/igl/list_to_matrix.h

@@ -26,6 +26,21 @@ namespace igl
   IGL_INLINE bool list_to_matrix(
     const std::vector<std::vector<T > > & V,
     Mat & M);
+  // Convert a list of row vectors of `n` or less to a matrix and pad on
+  // the right with `padding`:
+  //
+  // Inputs:
+  //   V  a m-long list of vectors of size <=n
+  //   n  number of columns
+  //   padding  value to fill in from right for short rows
+  // Outputs:
+  //   M  an m by n matrix
+  template <typename T, class Mat>
+  IGL_INLINE bool list_to_matrix(
+    const std::vector<std::vector<T > > & V,
+    const int n,
+    const T & padding,
+    Mat & M);
   // Vector wrapper
   template <typename T, class Mat>
   IGL_INLINE bool list_to_matrix(const std::vector<T > & V,Mat & M);

include/igl/matlab/mexStream.h

@@ -30,7 +30,6 @@ namespace igl
       inline virtual int overflow(int c = EOF);
   }; 
 }
-#ifndef IGL_STATIC_LIBRARY
-#  include "MexStream.cpp"
-#endif
+// ALWAYS INCLUDE
+#include "MexStream.cpp"
 #endif

include/igl/next_filename.cpp

@@ -0,0 +1,38 @@
+// 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 "next_filename.h"
+#include "STR.h"
+#include "file_exists.h"
+#include <cmath>
+#include <iomanip>
+
+bool igl::next_filename(
+  const std::string & prefix, 
+  const int zeros,
+  const std::string & suffix,
+  std::string & next)
+{
+  using namespace std;
+  // O(n), for huge lists could at least find bounds with exponential search
+  // and then narrow with binary search O(log(n))
+  int i = 0;
+  while(true)
+  {
+    next = STR(prefix << setfill('0') << setw(zeros)<<i<<suffix);
+    if(!file_exists(next.c_str()))
+    {
+      return true;
+    }
+    i++;
+    if(zeros > 0 && i >= pow(10,zeros))
+    {
+      return false;
+    }
+  }
+}
+

include/igl/next_filename.h

@@ -0,0 +1,36 @@
+// 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_NEXT_FILENAME_H
+#define IGL_NEXT_FILENAME_H
+#include "igl_inline.h"
+#include <string>
+namespace igl
+{
+  // Find the file with the first filename of the form
+  // "prefix-%0[zeros]dsuffix"
+  // 
+  // Inputs:
+  //   prefix  path to containing dir and filename prefix
+  //   zeros number of leading zeros as if digit printed with printf
+  //   suffix  suffix of filename and extension
+  // Outputs:
+  //   next  path to next file
+  // Returns true if found, false if exceeding range in zeros
+  IGL_INLINE bool next_filename(
+    const std::string & prefix, 
+    const int zeros,
+    const std::string & suffix,
+    std::string & next);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "next_filename.cpp"
+#endif
+
+#endif
+

include/igl/per_edge_normals.cpp

@@ -14,6 +14,7 @@ template <
 IGL_INLINE void igl::per_edge_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
+  const PerEdgeNormalsWeightingType weighting,
   Eigen::PlainObjectBase<DerivedN> & N,
   Eigen::PlainObjectBase<DerivedE> & E,
   Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
@@ -34,20 +35,51 @@ IGL_INLINE void igl::per_edge_normals(
   MatrixXd FN;
   per_face_normals(V,F,FN);
 
-  VectorXd dblA;
-  doublearea(V,F,dblA);
+  Eigen::VectorXd W(F.rows());
+  switch(weighting)
+  {
+    case PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM:
+      W.setConstant(1.);
+      break;
+    default:
+      assert(false && "Unknown weighting type");
+    case PER_EDGE_NORMALS_WEIGHTING_TYPE_DEFAULT:
+    case PER_EDGE_NORMALS_WEIGHTING_TYPE_AREA:
+    {
+      doublearea(V,F,W);
+      break;
+    }
+  }
+
   N.setConstant(E.rows(),3,0);
   for(int f = 0;f<m;f++)
   {
     for(int c = 0;c<3;c++)
     {
-      N.row(EMAP(f+c*m)) += dblA(f) * FN.row(f);
+      N.row(EMAP(f+c*m)) += W(f) * FN.row(f);
     }
   }
   N.rowwise().normalize();
   
 }
 
+template <
+  typename DerivedV, 
+  typename DerivedF, 
+  typename DerivedN,
+  typename DerivedE,
+  typename DerivedEMAP>
+IGL_INLINE void igl::per_edge_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedN> & N,
+  Eigen::PlainObjectBase<DerivedE> & E,
+  Eigen::PlainObjectBase<DerivedEMAP> & EMAP)
+{
+  return 
+    per_edge_normals(V,F,PER_EDGE_NORMALS_WEIGHTING_TYPE_DEFAULT,N,E,EMAP);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation
 template void igl::per_edge_normals<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<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);

include/igl/per_edge_normals.h

@@ -11,15 +11,39 @@
 #include <Eigen/Core>
 namespace igl
 {
+  enum PerEdgeNormalsWeightingType
+  {
+    // Incident face normals have uniform influence on edge normal
+    PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM = 0,
+    // Incident face normals are averaged weighted by area
+    PER_EDGE_NORMALS_WEIGHTING_TYPE_AREA = 1,
+    // Area weights
+    PER_EDGE_NORMALS_WEIGHTING_TYPE_DEFAULT = 2,
+    NUM_PER_EDGE_NORMALS_WEIGHTING_TYPE = 3
+  };
   // Compute face normals via vertex position list, face list
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F  #F by 3 eigen Matrix of face (triangle) indices
+  //   weighting  weighting type
   // Output:
   //   N  #2 by 3 matrix of mesh edge 3D normals per row
   //   E  #E by 2 matrix of edge indices per row
   //   EMAP  #E by 1 matrix of indices from all edges to E
   //
+  template <
+    typename DerivedV, 
+    typename DerivedF, 
+    typename DerivedN,
+    typename DerivedE,
+    typename DerivedEMAP>
+  IGL_INLINE void per_edge_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const PerEdgeNormalsWeightingType weight,
+    Eigen::PlainObjectBase<DerivedN> & N,
+    Eigen::PlainObjectBase<DerivedE> & E,
+    Eigen::PlainObjectBase<DerivedEMAP> & EMAP);
   template <
     typename DerivedV, 
     typename DerivedF, 

include/igl/per_face_normals.cpp

@@ -19,7 +19,7 @@ IGL_INLINE void igl::per_face_normals(
   N.resize(F.rows(),3);
   // loop over faces
   int Frows = F.rows();
-#pragma omp parallel for
+#pragma omp parallel for if (Frows>10000)
   for(int i = 0; i < Frows;i++)
   {
     const Eigen::Matrix<typename DerivedV::Scalar, 1, 3> v1 = V.row(F(i,1)) - V.row(F(i,0));

include/igl/per_vertex_normals.cpp

@@ -8,28 +8,62 @@
 #include "per_vertex_normals.h"
 
 #include "per_face_normals.h"
+#include "doublearea.h"
+#include "internal_angles.h"
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::per_vertex_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
+  const igl::PerVertexNormalsWeightingType weighting,
   Eigen::PlainObjectBase<DerivedV> & N)
 {
   Eigen::PlainObjectBase<DerivedV> PFN;
   igl::per_face_normals(V,F,PFN);
-  return igl::per_vertex_normals(V,F,PFN,N);
+  return per_vertex_normals(V,F,weighting,PFN,N);
 }
 
 template <typename DerivedV, typename DerivedF>
 IGL_INLINE void igl::per_vertex_normals(
   const Eigen::PlainObjectBase<DerivedV>& V,
   const Eigen::PlainObjectBase<DerivedF>& F,
+  Eigen::PlainObjectBase<DerivedV> & N)
+{
+  return per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT,N);
+}
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::per_vertex_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const igl::PerVertexNormalsWeightingType weighting,
   const Eigen::PlainObjectBase<DerivedV>& FN,
   Eigen::PlainObjectBase<DerivedV> & N)
 {
   // Resize for output
   N.setZero(V.rows(),3);
 
+  Eigen::MatrixXd W(F.rows(),3);
+  switch(weighting)
+  {
+    case PER_VERTEX_NORMALS_WEIGHTING_TYPE_UNIFORM:
+      W.setConstant(1.);
+      break;
+    default:
+      assert(false && "Unknown weighting type");
+    case PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT:
+    case PER_VERTEX_NORMALS_WEIGHTING_TYPE_AREA:
+    {
+      Eigen::VectorXd A;
+      doublearea(V,F,A);
+      W = A.replicate(1,3);
+      break;
+    }
+    case PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE:
+      internal_angles(V,F,W);
+      break;
+  }
+
   // loop over faces
   const int Frows = F.rows();
 //// Minimum number of iterms per openmp thread
@@ -44,15 +78,29 @@ IGL_INLINE void igl::per_vertex_normals(
     {
       // Does this need to be critical?
 //#pragma omp critical
-      N.row(F(i,j)) += FN.row(i);
+      N.row(F(i,j)) += W(i,j)*FN.row(i);
     }
   }
+  // take average via normalization
   N.rowwise().normalize();
 }
 
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE void igl::per_vertex_normals(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedF>& F,
+  const Eigen::PlainObjectBase<DerivedV>& FN,
+  Eigen::PlainObjectBase<DerivedV> & N)
+{
+  return
+    per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT,FN,N);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template void igl::per_vertex_normals<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&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, igl::PerVertexNormalsWeightingType, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template void igl::per_vertex_normals<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&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 template void igl::per_vertex_normals<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<unsigned int, -1, -1, 1, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3> >&);
 template void igl::per_vertex_normals<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&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);

include/igl/per_vertex_normals.h

@@ -9,19 +9,37 @@
 #define IGL_PER_VERTEX_NORMALS_H
 #include "igl_inline.h"
 #include <Eigen/Core>
-// Note: So for this only computes normals per vertex as uniformly weighted
-// averages of incident triangle normals. It would be nice to support more or
-// all of the methods here:
+// Note: It would be nice to support more or all of the methods here:
 // "A comparison of algorithms for vertex normal computation"
 namespace igl
 {
+  enum PerVertexNormalsWeightingType
+  {
+    // Incident face normals have uniform influence on vertex normal
+    PER_VERTEX_NORMALS_WEIGHTING_TYPE_UNIFORM = 0,
+    // Incident face normals are averaged weighted by area
+    PER_VERTEX_NORMALS_WEIGHTING_TYPE_AREA = 1,
+    // Incident face normals are averaged weighted by incident angle of vertex
+    PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE = 2,
+    // Area weights
+    PER_VERTEX_NORMALS_WEIGHTING_TYPE_DEFAULT = 3,
+    NUM_PER_VERTEX_NORMALS_WEIGHTING_TYPE = 4
+  };
   // Compute vertex normals via vertex position list, face list
   // Inputs:
   //   V  #V by 3 eigen Matrix of mesh vertex 3D positions
   //   F  #F by 3 eigne Matrix of face (triangle) indices
+  //   weighting  Weighting type
   // Output:
   //   N  #V by 3 eigen Matrix of mesh vertex 3D normals
   template <typename DerivedV, typename DerivedF>
+  IGL_INLINE void per_vertex_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const igl::PerVertexNormalsWeightingType weighting,
+    Eigen::PlainObjectBase<DerivedV> & N);
+  // Without weighting
+  template <typename DerivedV, typename DerivedF>
   IGL_INLINE void per_vertex_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,
@@ -29,6 +47,14 @@ namespace igl
   // Inputs:
   //   FN  #F by 3 matrix of face (triangle) normals
   template <typename DerivedV, typename DerivedF>
+  IGL_INLINE void per_vertex_normals(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedF>& F,
+    const PerVertexNormalsWeightingType weighting,
+    const Eigen::PlainObjectBase<DerivedV>& FN,
+    Eigen::PlainObjectBase<DerivedV> & N);
+  // Without weighting
+  template <typename DerivedV, typename DerivedF>
   IGL_INLINE void per_vertex_normals(
     const Eigen::PlainObjectBase<DerivedV>& V,
     const Eigen::PlainObjectBase<DerivedF>& F,

include/igl/polygon_mesh_to_triangle_mesh.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 "polygon_mesh_to_triangle_mesh.h"
+#include "matrix_to_list.h"
 
 template <typename Index, typename DerivedF>
 IGL_INLINE void igl::polygon_mesh_to_triangle_mesh(
@@ -54,7 +55,19 @@ IGL_INLINE void igl::polygon_mesh_to_triangle_mesh(
 
 }
 
+template <typename DerivedP, typename DerivedF>
+IGL_INLINE void igl::polygon_mesh_to_triangle_mesh(
+  const Eigen::PlainObjectBase<DerivedP>& P,
+  Eigen::PlainObjectBase<DerivedF>& F)
+{
+  std::vector<std::vector<typename DerivedP::Scalar> > vP;
+  matrix_to_list(P,vP);
+  return polygon_mesh_to_triangle_mesh(vP,F);
+}
+
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template instanciation
+// Explicit template specialization
+// generated by autoexplicit.sh
+template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<int, -1, 3, 0, -1, 3> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
 template void igl::polygon_mesh_to_triangle_mesh<int, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::vector<std::vector<int, std::allocator<int> >, std::allocator<std::vector<int, std::allocator<int> > > > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/polygon_mesh_to_triangle_mesh.h

@@ -35,6 +35,10 @@ namespace igl
   IGL_INLINE void polygon_mesh_to_triangle_mesh(
     const std::vector<std::vector<Index> > & vF,
     Eigen::PlainObjectBase<DerivedF>& F);
+  template <typename DerivedP, typename DerivedF>
+  IGL_INLINE void polygon_mesh_to_triangle_mesh(
+    const Eigen::PlainObjectBase<DerivedP>& P,
+    Eigen::PlainObjectBase<DerivedF>& F);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/principal_curvature.cpp

@@ -307,11 +307,11 @@ public:
 
 IGL_INLINE CurvatureCalculator::CurvatureCalculator()
 {
-  this->localMode=false;
-  this->projectionPlaneCheck=false;
+  this->localMode=true;
+  this->projectionPlaneCheck=true;
   this->sphereRadius=5;
   this->st=SPHERE_SEARCH;
-  this->nt=PROJ_PLANE;
+  this->nt=AVERAGE;
   this->montecarlo=false;
   this->montecarloN=0;
   this->kRing=3;
@@ -418,8 +418,8 @@ IGL_INLINE void CurvatureCalculator::finalEigenStuff (int i, std::vector<Eigen::
   v2global.normalize();
 
   v1global *= c_val(0);
-  v2global *= c_val(1);
-
+  v2global *= c_val(1);  
+  
   if (c_val[0] > c_val[1])
   {
     curv[i]=std::vector<double>(2);
@@ -697,6 +697,18 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
       std::cerr << "Could not compute curvature of radius " << scaledRadius << endl;
       return;
     }
+
+
+    if (projectionPlaneCheck)
+    {
+      vvtmp.reserve (vv.size ());
+      applyProjOnPlane (vertex_normals.row(i), vv, vvtmp);
+      if (vvtmp.size() >= 6 && vvtmp.size()<vv.size())
+	  vv = vvtmp;
+
+    }
+
+
     switch (nt)
     {
       case AVERAGE:
@@ -709,15 +721,6 @@ IGL_INLINE void CurvatureCalculator::computeCurvature()
         fprintf(stderr,"Error: normal type not recognized");
         return;
     }
-
-    if (projectionPlaneCheck)
-    {
-      vvtmp.reserve (vv.size ());
-      applyProjOnPlane (normal, vv, vvtmp);
-      if (vvtmp.size() >= 6)
-        vv = vvtmp;
-    }
-
     if (vv.size()<6)
     {
       std::cerr << "Could not compute curvature of radius " << scaledRadius << endl;

include/igl/project.cpp

@@ -146,6 +146,7 @@ template int igl::project<Eigen::Matrix<float, 3, 1, 0, 3, 1>, Eigen::Matrix<flo
 template Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > igl::project<Eigen::Matrix<float, 3, 1, 0, 3, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<float, 3, 1, 0, 3, 1> > const&);
 template Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > igl::project<Eigen::Matrix<double, 1, -1, 1, 1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, -1, 1, 1, -1> > const&);
 template int igl::project<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> >&);
+template Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > igl::project<Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&);
 #endif
 #endif
 

include/igl/readDMAT.cpp

@@ -59,7 +59,7 @@ template <typename DerivedW>
 IGL_INLINE bool igl::readDMAT(const std::string file_name,
   Eigen::PlainObjectBase<DerivedW> & W)
 {
-  FILE * fp = fopen(file_name.c_str(),"r");
+  FILE * fp = fopen(file_name.c_str(),"rb");
   if(fp == NULL)
   {
     fprintf(stderr,"IOError: readDMAT() could not open %s...\n",file_name.c_str());

include/igl/readMESH.cpp

@@ -467,5 +467,7 @@ IGL_INLINE bool igl::readMESH(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
 // generated by autoexplicit.sh
+template bool igl::readMESH<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&);
+// generated by autoexplicit.sh
 template bool igl::readMESH<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> >, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/readOBJ.cpp

@@ -368,16 +368,17 @@ IGL_INLINE bool igl::readOBJ(
     // igl::list_to_matrix(vF,F) already printed error message to std err
     return false;
   }
-  // Legacy
-  if(F.cols() != 3)
-  {
-    fprintf(stderr,
-            "Error: readOBJ(filename,V,F) is meant for reading triangle-only"
-            " meshes. This mesh has faces all with size %d. See readOBJ.h for other"
-            " options.\n",
-            (int)F.cols());
-    return false;
-  }
+  // THIS IS ANNOYING, I WANT TO READ QUADS, Why not?
+  //// Legacy
+  //if(F.cols() != 3)
+  //{
+  //  fprintf(stderr,
+  //          "Error: readOBJ(filename,V,F) is meant for reading triangle-only"
+  //          " meshes. This mesh has faces all with size %d. See readOBJ.h for other"
+  //          " options.\n",
+  //          (int)F.cols());
+  //  return false;
+  //}
   return true;
 }
 

include/igl/readSTL.cpp

@@ -214,6 +214,9 @@ close_true:
 }
 
 #ifdef IGL_STATIC_LIBRARY
-// Explicit template instanciation
+// Explicit template specialization
+// generated by autoexplicit.sh
+template bool igl::readSTL<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+// generated by autoexplicit.sh
 template bool igl::readSTL<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::basic_string<char, std::char_traits<char>, std::allocator<char> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/read_triangle_mesh.cpp

@@ -7,11 +7,18 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "read_triangle_mesh.h"
 
+#include <igl/list_to_matrix.h>
+#include <igl/readMESH.h>
 #include <igl/readOBJ.h>
 #include <igl/readOFF.h>
+#include <igl/readSTL.h>
+#include <igl/readWRL.h>
 #include <igl/pathinfo.h>
+#include <igl/boundary_facets.h>
+#include <igl/polygon_mesh_to_triangle_mesh.h>
 
 #include <cstdio>
+#include <algorithm>
 #include <iostream>
 
 
@@ -50,25 +57,82 @@ IGL_INLINE bool igl::read_triangle_mesh(
   Eigen::PlainObjectBase<DerivedV>& V,
   Eigen::PlainObjectBase<DerivedF>& F)
 {
-    const char* p;
-    for (p = str.c_str(); *p != '\0'; p++)
-        ;
-    while (*p != '.')
-        p--;
+  std::string _1,_2,_3,_4;
+  return read_triangle_mesh(str,V,F,_1,_2,_3,_4);
+}
+
+template <typename DerivedV, typename DerivedF>
+IGL_INLINE bool igl::read_triangle_mesh(
+  const std::string filename,
+  Eigen::PlainObjectBase<DerivedV>& V,
+  Eigen::PlainObjectBase<DerivedF>& F,
+  std::string & dir,
+  std::string & base,
+  std::string & ext,
+  std::string & name)
+{
+  using namespace std;
+  using namespace Eigen;
 
-    if (!strcmp(p, ".obj") || !strcmp(p, ".OBJ"))
+  // dirname, basename, extension and filename
+  pathinfo(filename,dir,base,ext,name);
+  // Convert extension to lower case
+  transform(ext.begin(), ext.end(), ext.begin(), ::tolower);
+  vector<vector<double > > vV,vN,vTC;
+  vector<vector<int > > vF,vFTC,vFN;
+  if(ext == "mesh")
+  {
+    // Convert extension to lower case
+    MatrixXi T;
+    if(!readMESH(filename,V,T,F))
+    {
+      return 1;
+    }
+    //if(F.size() > T.size() || F.size() == 0)
+    {
+      boundary_facets(T,F);
+    }
+  }else if(ext == "obj")
+  {
+    if(!readOBJ(filename,vV,vTC,vN,vF,vFTC,vFN))
+    {
+      return false;
+    }
+  }else if(ext == "off")
+  {
+    if(!readOFF(filename,vV,vF,vN))
+    {
+      return false;
+    }
+  }else if(ext == "stl")
+  {
+    MatrixXd _;
+    if(!readSTL(filename,V,F,_))
     {
-        return igl::readOBJ(str,V,F);
-    }else if (!strcmp(p, ".off") || !strcmp(p, ".OFF"))
+      return false;
+    }
+  }else if(ext == "wrl")
+  {
+    if(!readWRL(filename,vV,vF))
     {
-        return igl::readOFF(str,V,F);
+      return false;
     }
-    else
+  }else
+  {
+    cerr<<"Error: unknown extension: "<<ext<<endl;
+    return false;
+  }
+  if(vV.size() > 0)
+  {
+    if(!list_to_matrix(vV,V))
     {
-      fprintf(stderr,"read_triangle_mesh() does not recognize extension: %s\n",p);
       return false;
     }
+    polygon_mesh_to_triangle_mesh(vF,F);
+  }
+  return true;
 }
+
 #endif
 
 #ifdef IGL_STATIC_LIBRARY

include/igl/read_triangle_mesh.h

@@ -26,7 +26,6 @@ namespace igl
   //     to Scalar)
   //   Index  type for indices (will be read as int and cast to Index)
   // Inputs:
-  // Inputs:
   //   str  path to .obj/.off file
   // Outputs:
   //   V  eigen double matrix #V by 3
@@ -42,6 +41,20 @@ namespace igl
     const std::string str,
     Eigen::PlainObjectBase<DerivedV>& V,
     Eigen::PlainObjectBase<DerivedF>& F);
+  // Outputs:
+  //  dir  directory path (see pathinfo.h)
+  //  base  base name (see pathinfo.h)
+  //  ext  extension (see pathinfo.h)
+  //  name  filename (see pathinfo.h)
+  template <typename DerivedV, typename DerivedF>
+  IGL_INLINE bool read_triangle_mesh(
+    const std::string str,
+    Eigen::PlainObjectBase<DerivedV>& V,
+    Eigen::PlainObjectBase<DerivedF>& F,
+    std::string & dir,
+    std::string & base,
+    std::string & ext,
+    std::string & name);
 #endif
 }
 

include/igl/remove_unreferenced.cpp

@@ -7,13 +7,18 @@
 // obtain one at http://mozilla.org/MPL/2.0/.
 #include "remove_unreferenced.h"
 
-template <typename Scalar, typename Index>
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedNV,
+  typename DerivedNF,
+  typename DerivedI>
 IGL_INLINE void igl::remove_unreferenced(
-  const Eigen::PlainObjectBase<Scalar> &V,
-  const Eigen::PlainObjectBase<Index> &F,
-  Eigen::PlainObjectBase<Scalar> &NV,
-  Eigen::PlainObjectBase<Index> &NF,
-  Eigen::PlainObjectBase<Index> &I)
+  const Eigen::PlainObjectBase<DerivedV> &V,
+  const Eigen::PlainObjectBase<DerivedF> &F,
+  Eigen::PlainObjectBase<DerivedNV> &NV,
+  Eigen::PlainObjectBase<DerivedNF> &NF,
+  Eigen::PlainObjectBase<DerivedI> &I)
 {
 
   // Mark referenced vertices
@@ -65,4 +70,5 @@ IGL_INLINE void igl::remove_unreferenced(
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template specialization
+template void igl::remove_unreferenced<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<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 #endif

include/igl/remove_unreferenced.h

@@ -27,13 +27,18 @@ namespace igl
   // Output:
   // NV, NF: new mesh without unreferenced vertices
   //
-  template <typename Scalar, typename Index>
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedNV,
+    typename DerivedNF,
+    typename DerivedI>
   IGL_INLINE void remove_unreferenced(
-    const Eigen::PlainObjectBase<Scalar> &V,
-    const Eigen::PlainObjectBase<Index> &F,
-    Eigen::PlainObjectBase<Scalar> &NV,
-    Eigen::PlainObjectBase<Index> &NF,
-    Eigen::PlainObjectBase<Index> &I);
+    const Eigen::PlainObjectBase<DerivedV> &V,
+    const Eigen::PlainObjectBase<DerivedF> &F,
+    Eigen::PlainObjectBase<DerivedNV> &NV,
+    Eigen::PlainObjectBase<DerivedNF> &NF,
+    Eigen::PlainObjectBase<DerivedI> &I);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/tetgen/mesh_with_skeleton.cpp

@@ -11,6 +11,7 @@
 #include <igl/cat.h>
 #include <igl/tetgen/tetrahedralize.h>
 #include <igl/writeOFF.h>
+#include <igl/writeOBJ.h>
 
 #include <iostream>
 // Default settings pq2Y tell tetgen to mesh interior of triangle mesh and
@@ -43,15 +44,13 @@ IGL_INLINE bool igl::mesh_with_skeleton(
   sample_edges(C,BECE,samples_per_bone,S);
   // Vertices we'll constrain tet mesh to meet
   MatrixXd VS = cat(1,V,S);
-  // Boundary faces
-  MatrixXi BF;
   // Use tetgen to mesh the interior of surface, this assumes surface:
   //   * has no holes
   //   * has no non-manifold edges or vertices
   //   * has consistent orientation
   //   * has no self-intersections
   //   * has no 0-volume pieces
-  //writeOFF("mesh_with_skeleton.off",VS,F);
+  //writeOBJ("mesh_with_skeleton.obj",VS,F);
   cerr<<"tetgen begin()"<<endl;
   int status = tetrahedralize( VS,F,eff_tetgen_flags,VV,TT,FF);
   cerr<<"tetgen end()"<<endl;

include/igl/triangle_fan.cpp

@@ -20,7 +20,8 @@ IGL_INLINE void igl::triangle_fan(
   // outgoing (left or right). Thus this will not work.
   assert(E.cols() == 2);
   // Arbitrary starting vertex
-  int s = E(int(((double)rand() / RAND_MAX)*E.rows()),0);
+  //int s = E(int(((double)rand() / RAND_MAX)*E.rows()),0);
+  int s = E(rand()%E.rows(),0);
   vector<vector<int> >  lcap;
   for(int i = 0;i<E.rows();i++)
   {

include/igl/triangle_fan.h

@@ -18,6 +18,6 @@ namespace igl
   IGL_INLINE Eigen::MatrixXi triangle_fan( const Eigen::MatrixXi & E);
 }
 #ifndef IGL_STATIC_LIBRARY
-#  include "triangle_fan.h"
+#  include "triangle_fan.cpp"
 #endif
 #endif

include/igl/viewer/TODOs.txt

@@ -1,4 +1,5 @@
 - `align_and_center_object` continues to zoom out on repeated calls
+- trackball_angle should be a quaternion
 - data.lines, data.points should not concatenate colors with coordinates
 - snap to canonical recenters origin but trackball does not
 - rewrite in libigl style

include/igl/viewer/Viewer.cpp

@@ -410,6 +410,11 @@ namespace igl
                " indices'");
 
     core.init();
+
+    if (callback_init)
+      if (callback_init(*this))
+        return;
+
     init_plugins();
   }
 
@@ -425,6 +430,7 @@ namespace igl
     data.set_face_based(false);
 
     // C-style callbacks
+    callback_init         = 0;
     callback_pre_draw     = 0;
     callback_post_draw    = 0;
     callback_mouse_down   = 0;
@@ -434,14 +440,15 @@ namespace igl
     callback_key_down     = 0;
     callback_key_up       = 0;
 
-    callback_pre_draw_data = 0;
-    callback_post_draw     = 0;
-    callback_mouse_down    = 0;
-    callback_mouse_up      = 0;
-    callback_mouse_move    = 0;
-    callback_mouse_scroll  = 0;
-    callback_key_down      = 0;
-    callback_key_up        = 0;
+    callback_init_data          = 0;
+    callback_pre_draw_data      = 0;
+    callback_post_draw_data     = 0;
+    callback_mouse_down_data    = 0;
+    callback_mouse_up_data      = 0;
+    callback_mouse_move_data    = 0;
+    callback_mouse_scroll_data  = 0;
+    callback_key_down_data      = 0;
+    callback_key_up_data        = 0;
 
   }
 

include/igl/viewer/Viewer.h

@@ -108,6 +108,7 @@ namespace igl
 
 
     // C-style callbacks
+    bool (*callback_init)(Viewer& viewer);
     bool (*callback_pre_draw)(Viewer& viewer);
     bool (*callback_post_draw)(Viewer& viewer);
     bool (*callback_mouse_down)(Viewer& viewer, int button, int modifier);
@@ -118,6 +119,7 @@ namespace igl
     bool (*callback_key_up)(Viewer& viewer, unsigned char key, int modifiers);
 
     // Pointers to per-callback data
+    void* callback_init_data;
     void* callback_pre_draw_data;
     void* callback_post_draw_data;
     void* callback_mouse_down_data;

include/igl/volume.cpp

@@ -22,10 +22,10 @@ IGL_INLINE void igl::volume(
   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));
+    const RowVector3d & a = V.row(T(t,0));
+    const RowVector3d & b = V.row(T(t,1));
+    const RowVector3d & c = V.row(T(t,2));
+    const RowVector3d & d = V.row(T(t,3));
     vol(t) = -(a-d).dot((b-d).cross(c-d))/6.;
   }
 }
@@ -111,4 +111,5 @@ template void igl::volume<Eigen::Matrix<double, -1, 6, 0, -1, 6>, Eigen::Matrix<
 template Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar igl::volume_single<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3> >(Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3> const&);
 template Eigen::Matrix<double, 3, 1, 0, 3, 1>::Scalar igl::volume_single<Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1>, Eigen::Matrix<double, 3, 1, 0, 3, 1> >(Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&, Eigen::Matrix<double, 3, 1, 0, 3, 1> const&);
 template void igl::volume<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, 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> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::volume<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/winding_number.cpp

@@ -36,6 +36,7 @@ IGL_INLINE void igl::winding_number(
         Vector3d p = O.row(o);
         W(o) = hier.winding_number(p);
       }
+      break;
     }
     default: assert(false && "Bad simplex size"); break;
   }

include/igl/winding_number.h

@@ -61,7 +61,7 @@ namespace igl
 }
 
 #ifndef IGL_STATIC_LIBRARY
-#  include "winding_number.h"
+#  include "winding_number.cpp"
 #endif
 
 #endif

include/igl/writeMESH.cpp

@@ -95,7 +95,7 @@ IGL_INLINE bool igl::writeMESH(
   for(int i = 0;i<number_of_tet_vertices;i++)
   {
     // print position of ith tet vertex
-    fprintf(mesh_file,"%lg %lg %lg 1\n",
+    fprintf(mesh_file,"%.17lg %.17lg %.17lg 1\n",
       (double)V(i,0),
       (double)V(i,1),
       (double)V(i,2));

include/igl/writeOBJ.cpp

@@ -8,6 +8,8 @@
 #include "writeOBJ.h"
 
 #include <iostream>
+#include <limits>
+#include <iomanip>
 #include <fstream>
 #include <cstdio>
 
@@ -18,6 +20,7 @@ IGL_INLINE bool igl::writeOBJ(
   const Eigen::PlainObjectBase<DerivedF>& F)
 {
   std::ofstream s(str.c_str());
+  s.precision(std::numeric_limits<double>::digits10 + 1);
 
   if(!s.is_open())
   {
@@ -32,7 +35,16 @@ IGL_INLINE bool igl::writeOBJ(
   
   for(int i=0;i<(int)F.rows();++i)
   {
-    s << "f " << F(i,0)+1 << " " << F(i,1)+1 << " " << F(i,2)+1 << std::endl;
+    s << "f ";
+    for(int c =0;c<(int)F.cols();++c)
+    {
+      if(c>0)
+      {
+        s<<" ";
+      }
+      s<< F(i,c)+1;
+    }
+    s<<std::endl;
   }
   
   s.close();