Browse Source

skeleton building and posing examples

Former-commit-id: 98c3146acce92f56b5311e9d329210f4d882d647
Alec Jacobson 10 years ago
parent
commit
5e41228487

+ 8 - 0
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');
+  %

+ 2 - 2
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);

+ 76 - 0
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

+ 1211 - 0
examples/skeleton-builder/example.cpp

@@ -0,0 +1,1211 @@
+#include <igl/draw_skeleton_3d.h>
+#include <igl/draw_skeleton_vector_graphics.h>
+#include <igl/two_axis_valuator_fixed_up.h>
+#include <igl/readOBJ.h>
+#include <igl/readTGF.h>
+#include <igl/writeOBJ.h>
+#include <igl/writeOFF.h>
+#include <igl/readWRL.h>
+#include <igl/report_gl_error.h>
+#include <igl/polygon_mesh_to_triangle_mesh.h>
+#include <igl/readOFF.h>
+#include <igl/readMESH.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;
+
+  // dirname, basename, extension and filename
+  string dir,base,ext,name;
+  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 == "obj")
+  {
+    // Convert extension to lower case
+    if(!igl::readOBJ(filename,vV,vTC,vN,vF,vFTC,vFN))
+    {
+      return 1;
+    }
+  }else if(ext == "off")
+  {
+    // Convert extension to lower case
+    if(!igl::readOFF(filename,vV,vF,vN))
+    {
+      return 1;
+    }
+  }else if(ext == "wrl")
+  {
+    // Convert extension to lower case
+    if(!igl::readWRL(filename,vV,vF))
+    {
+      return 1;
+    }
+  //}else
+  //{
+  //  // Convert extension to lower case
+  //  MatrixXi T;
+  //  if(!igl::readMESH(filename,V,T,F))
+  //  {
+  //    return 1;
+  //  }
+  //  //if(F.size() > T.size() || F.size() == 0)
+  //  {
+  //    boundary_facets(T,F);
+  //  }
+  }
+  if(vV.size() > 0)
+  {
+    if(!list_to_matrix(vV,V))
+    {
+      return 1;
+    }
+    polygon_mesh_to_triangle_mesh(vF,F);
+  }
+  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;
+}

+ 83 - 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

+ 982 - 0
examples/skeleton-posing/example.cpp

@@ -0,0 +1,982 @@
+#include <igl/draw_skeleton_3d.h>
+#include <igl/draw_skeleton_vector_graphics.h>
+#include <igl/two_axis_valuator_fixed_up.h>
+#include <igl/readOBJ.h>
+#include <igl/readTGF.h>
+#include <igl/writeOBJ.h>
+#include <igl/writeOFF.h>
+#include <igl/writeDMAT.h>
+#include <igl/readDMAT.h>
+#include <igl/readWRL.h>
+#include <igl/report_gl_error.h>
+#include <igl/polygon_mesh_to_triangle_mesh.h>
+#include <igl/readOFF.h>
+#include <igl/readMESH.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/lbs_matrix.h>
+#include <igl/colon.h>
+#include <igl/writeTGF.h>
+#include <igl/file_exists.h>
+#include <igl/MouseController.h>
+#include <igl/STR.h>
+#include <igl/bone_parents.h>
+#include <igl/cgal/remesh_self_intersections.h>
+#include <igl/tetgen/tetrahedralize.h>
+#include <igl/barycenter.h>
+#include <igl/remove_unreferenced.h>
+#include <igl/boundary_facets.h>
+#include <igl/winding_number.h>
+#include <igl/bbw/bbw.h>
+#include <igl/normalize_row_sums.h>
+#include <igl/centroid.h>
+#include <igl/tetgen/mesh_with_skeleton.h>
+#include <igl/draw_beach_ball.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>
+
+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;
+}
+
+std::string next_filename(
+  const std::string & prefix, 
+  const int zeros,
+  const std::string & suffix)
+{
+  using namespace std;
+  using namespace igl;
+  // O(n)
+  int i = 0;
+  while(true)
+  {
+    string next = STR(prefix << setfill('0') << setw(zeros)<<i<<suffix);
+    if(!file_exists(next.c_str()))
+    {
+      return next;
+    }
+    i++;
+  }
+}
+
+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<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());
+  }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");
+  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,_3;
+    remesh_self_intersections(V,F,{},CV,CF,_1,_2,_3);
+  }
+  MatrixXd TV;
+  MatrixXi TT;
+  {
+    MatrixXi _1;
+    // c  convex hull
+    // Y  no boundary steiners
+    // p  polygon input
+    if(tetrahedralize(V,F,"cYp",TV,TT,_1) != 0)
+    {
+      cout<<REDRUM("CDT failed.")<<endl;
+      return false;
+    }
+  }
+  {
+    MatrixXd BC;
+    barycenter(TV,TT,BC);
+    VectorXd W;
+    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;
+    if(!mesh_with_skeleton(V,F,C,{},BE,{},10,TV,TT,_1))
+    {
+      cout<<REDRUM("tetgen failed.")<<endl;
+      return false;
+    }
+  }
+  // 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;
+
+  // dirname, basename, extension and filename
+  string dir,base,ext,name;
+  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 == "obj")
+  {
+    // Convert extension to lower case
+    if(!igl::readOBJ(filename,vV,vTC,vN,vF,vFTC,vFN))
+    {
+      return 1;
+    }
+  }else if(ext == "off")
+  {
+    // Convert extension to lower case
+    if(!igl::readOFF(filename,vV,vF,vN))
+    {
+      return 1;
+    }
+  }else if(ext == "wrl")
+  {
+    // Convert extension to lower case
+    if(!igl::readWRL(filename,vV,vF))
+    {
+      return 1;
+    }
+  //}else
+  //{
+  //  // Convert extension to lower case
+  //  MatrixXi T;
+  //  if(!igl::readMESH(filename,V,T,F))
+  //  {
+  //    return 1;
+  //  }
+  //  //if(F.size() > T.size() || F.size() == 0)
+  //  {
+  //    boundary_facets(T,F);
+  //  }
+  }
+  if(vV.size() > 0)
+  {
+    if(!list_to_matrix(vV,V))
+    {
+      return 1;
+    }
+    polygon_mesh_to_triangle_mesh(vF,F);
+  }
+  if(output_prefix.size() == 0)
+  {
+    output_prefix = dir+"/"+name+"-pose-";
+  }
+
+  {
+    string output_filename = next_filename(output_prefix,4,".dmat");
+    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;
+}
+

+ 1 - 1
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)/

+ 10 - 0
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);

+ 14 - 0
include/igl/MouseController.h

@@ -122,6 +122,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 +501,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())

+ 25 - 0
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

+ 32 - 0
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
+

+ 6 - 0
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)

+ 9 - 5
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,

+ 63 - 0
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

+ 49 - 0
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
+

+ 5 - 1
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());

+ 1 - 2
include/igl/draw_skeleton_3d.h

@@ -29,13 +29,12 @@ 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,

+ 6 - 6
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

+ 4 - 4
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);

+ 25 - 0
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

+ 9 - 0
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

+ 1 - 0
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
 

+ 1 - 1
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());

+ 12 - 6
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

+ 11 - 6
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

+ 1 - 0
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