skeleton example

Former-commit-id: 00597b8f23225ed8e47e04622e3f462e361e504d

examples/skeleton/Makefile

@@ -0,0 +1,72 @@
+
+.PHONY: all
+
+# Shared flags etc.
+include ../../Makefile.conf
+
+all: obj example
+
+.PHONY: example
+
+LIBIGL=../../
+LIBIGL_INC=-I$(LIBIGL)/include
+LIBIGL_LIB=-L$(LIBIGL)/lib -ligl -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)
+
+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=../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
+
+INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(MATLAB_INC) $(GLUT_INC) ${TETGEN_INC} $(MOSEK_INC)
+LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(MATLAB_LIB) $(CARBON_LIB) $(TETGEN_LIB) $(MOSEK_LIB)
+
+CPP_FILES=$(wildcard ./*.cpp)
+OBJ_FILES=$(addprefix obj/,$(notdir $(CPP_FILES:.cpp=.o))) 
+
+CFLAGS+=-std=c++11 
+
+example: obj $(OBJ_FILES)
+	g++ $(OPENMP) $(AFLAGS) $(CFLAGS) -o example $(OBJ_FILES) $(LIB)
+
+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/draw_skeleton.h

@@ -0,0 +1,119 @@
+#include <Eigen/Core>
+// Draw a skeleton
+//
+// Inputs:
+//   C  #C by dim List of joint rest positions
+//   BE  #BE by 2 list of bone edge indices into C
+template <typename DerivedC, typename DerivedBE>
+void draw_skeleton(
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::PlainObjectBase<DerivedBE> & BE);
+template <typename DerivedC, typename DerivedBE, typename DerivedT>
+void draw_skeleton(
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::PlainObjectBase<DerivedBE> & BE,
+  const Eigen::PlainObjectBase<DerivedT> & T);
+
+// Implementation
+#include <igl/PI.h>
+#include <igl/OpenGL_convenience.h>
+#include <Eigen/Geometry>
+#include <iostream>
+template <typename DerivedC, typename DerivedBE>
+void draw_skeleton(
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::PlainObjectBase<DerivedBE> & BE)
+{
+  using namespace Eigen;
+  typedef Eigen::Matrix<typename DerivedC::Scalar,Dynamic,Dynamic> Mat;
+  Mat I = Mat::Identity(C.cols()+1,C.cols());
+  Mat T = I.replicate(BE.rows(),1);
+  return draw_skeleton(C,BE,T);
+}
+
+template <typename DerivedC, typename DerivedBE, typename DerivedT>
+void draw_skeleton(
+  const Eigen::PlainObjectBase<DerivedC> & C,
+  const Eigen::PlainObjectBase<DerivedBE> & BE,
+  const Eigen::PlainObjectBase<DerivedT> & T)
+{
+  using namespace Eigen;
+  using namespace std;
+  glDisable(GL_LIGHTING);
+  glColor3d(70./255.,252./255.,167./255.);
+  glLineWidth(2.0);
+
+  auto draw_sphere = [](const double r)
+  {
+    auto draw_circle = []()
+    {
+      glBegin(GL_LINE_STRIP);
+      for(double th = 0;th<2.0*igl::PI;th+=(2.0*igl::PI/30.0))
+      {
+        glVertex3d(cos(th),sin(th),0.0);
+      }
+      glVertex3d(cos(0),sin(0),0.0);
+      glEnd();
+    };
+    glPushMatrix();
+    glScaled(r,r,r);
+    draw_circle();
+    glRotated(90.0,1.0,0.0,0.0);
+    draw_circle();
+    glRotated(90.0,0.0,1.0,0.0);
+    draw_circle();
+    glPopMatrix();
+  };
+  auto draw_pyramid = [](const double r)
+  {
+    glBegin(GL_LINE_STRIP);
+    glVertex3d(0, 1,-1);
+    glVertex3d(0,-1,-1);
+    glVertex3d(0,-1, 1);
+    glVertex3d(0, 1, 1);
+    glVertex3d(0, 1,-1);
+    glEnd();
+    glBegin(GL_LINES);
+    glVertex3d(0, 1,-1);
+    glVertex3d(1,0,0);
+    glVertex3d(0,-1,-1);
+    glVertex3d(1,0,0);
+    glVertex3d(0,-1, 1);
+    glVertex3d(1,0,0);
+    glVertex3d(0, 1, 1);
+    glVertex3d(1,0,0);
+    glEnd();
+  };
+
+  // Loop over bones
+  for(int e = 0;e < BE.rows();e++)
+  {
+    // Draw a sphere
+    auto s = C.row(BE(e,0));
+    auto d = C.row(BE(e,1));
+    auto b = (d-s).transpose().eval();
+    double r = 0.02;
+    Matrix4d Te = Matrix4d::Identity();
+    Te.block(0,0,3,4) = T.block(e*4,0,4,3).transpose();
+    Quaterniond q;
+    q.setFromTwoVectors(Vector3d(1,0,0),b);
+    glPushMatrix();
+    glMultMatrixd(Te.data());
+    glTranslated(s(0),s(1),s(2));
+    draw_sphere(r);
+    const double len = b.norm()-2.*r;
+    if(len>=0)
+    {
+      auto u = b.normalized()*r;
+      glPushMatrix();
+      glTranslated(u(0),u(1),u(2));
+      glMultMatrixd(Affine3d(q).matrix().data());
+      glScaled(b.norm()-2.*r,r,r);
+      draw_pyramid(r);
+      glPopMatrix();
+    }
+    glTranslated(b(0),b(1),b(2));
+    draw_sphere(r);
+    glPopMatrix();
+  }
+}

examples/skeleton/example.cpp

@@ -0,0 +1,676 @@
+#include "draw_skeleton.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/triangulate.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/tetgen/mesh_with_skeleton.h>
+#include <igl/boundary_conditions.h>
+#include <igl/normalize_row_sums.h>
+#include <igl/lbs_matrix.h>
+#include <igl/bbw/bbw.h>
+
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+#include <GLUT/glut.h>
+
+#include <string>
+#include <vector>
+#include <stack>
+#include <iostream>
+
+
+Eigen::MatrixXd V,N,C,W,M;
+Eigen::VectorXd Vmid,Vmin,Vmax;
+double bbd = 1.0;
+Eigen::MatrixXi F,BE;
+Eigen::VectorXi P;
+struct State
+{
+  igl::Camera camera;
+} s;
+
+bool wireframe = false;
+bool skeleton_on_top = 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;
+
+std::stack<State> undo_stack;
+std::stack<State> redo_stack;
+
+bool is_rotating = false;
+int down_x,down_y;
+igl::Camera down_camera;
+
+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;
+
+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 = s.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);
+  s.camera.m_aspect = (double)width/(double)height;
+}
+
+void push_scene()
+{
+  using namespace igl;
+  using namespace std;
+  glMatrixMode(GL_PROJECTION);
+  glPushMatrix();
+  glLoadIdentity();
+  auto & camera = s.camera;
+  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] = {30.0/255.0,30.0/255.0,50.0/255.0,0};
+  glClearColor(back[0],back[1],back[2],0);
+  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
+
+  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();
+    auto & camera = s.camera;
+    camera.orbit(q.conjugate());
+  }
+
+  glEnable(GL_DEPTH_TEST);
+  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);
+  draw_floor(GREY,DARK_GREY);
+  glPopMatrix();
+
+  push_object();
+
+  // Set material properties
+  glDisable(GL_COLOR_MATERIAL);
+  glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT,  SILVER_AMBIENT);
+  glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE,  SILVER_DIFFUSE  );
+  glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, SILVER_SPECULAR);
+  glMaterialf (GL_FRONT_AND_BACK, GL_SHININESS, 128);
+
+  vector<Quaterniond> dQ(BE.rows(),Quaterniond::Identity()),vQ;
+  vector<Vector3d> vT;
+  Matrix3d A = Matrix3d::Identity();
+  for(int e = 0;e<BE.rows();e++)
+  {
+    dQ[e] = AngleAxisd((sin(get_seconds()+e))*0.06*PI,A.col(e%3));
+  }
+  forward_kinematics(C,BE,P,dQ,vQ,vT);
+  const int dim = C.cols();
+  MatrixXd T(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);
+  }
+  
+  if(wireframe)
+  {
+    glPolygonMode(GL_FRONT_AND_BACK,GL_LINE);
+  }
+  glLineWidth(1.0);
+  MatrixXd U = M*T;
+  per_face_normals(U,F,N);
+  draw_mesh(U,F,N);
+  glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
+
+  if(skeleton_on_top)
+  {
+    glDisable(GL_DEPTH_TEST);
+  }
+
+  draw_skeleton(C,BE,T);
+
+  pop_object();
+
+  pop_scene();
+
+  report_gl_error();
+
+  TwDraw();
+  glutSwapBuffers();
+  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();
+
+  auto & camera = s.camera;
+  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);
+  }
+}
+
+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);
+  switch(glutButton)
+  {
+    case GLUT_RIGHT_BUTTON:
+    case GLUT_LEFT_BUTTON:
+    {
+      switch(glutState)
+      {
+        case 1:
+          // up
+          glutSetCursor(GLUT_CURSOR_INHERIT);
+          is_rotating = false;
+          break;
+        case 0:
+          if(!tw_using)
+          {
+            push_undo();
+            glutSetCursor(GLUT_CURSOR_CYCLE);
+            // collect information for trackball
+            is_rotating = true;
+            down_camera = s.camera;
+            down_x = mouse_x;
+            down_y = mouse_y;
+          }
+        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;
+    }
+  }
+}
+
+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;
+    auto & camera = s.camera;
+    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());
+  }
+}
+
+void init_relative()
+{
+  using namespace Eigen;
+  using namespace igl;
+  per_face_normals(V,F,N);
+  Vmax = V.colwise().maxCoeff();
+  Vmin = V.colwise().minCoeff();
+  Vmid = 0.5*(Vmax + Vmin);
+  bbd = (Vmax-Vmin).norm();
+}
+
+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 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 'z':
+    case 'Z':
+      if(command_down)
+      {
+        if(shift_down)
+        {
+          redo();
+        }else
+        {
+          undo();
+        }
+        break;
+      }else
+      {
+        push_undo();
+        Quaterniond q;
+        snap_to_canonical_view_quat(s.camera.m_rotation_conj,1.0,q);
+        s.camera.orbit(q.conjugate());
+      }
+    default:
+      if(!TwEventKeyboardGLUT(key,mouse_x,mouse_y))
+      {
+        cout<<"Unknown key command: "<<key<<" "<<int(key)<<endl;
+      }
+  }
+  
+}
+
+bool init_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;
+  // Mesh with samples on skeleton
+  // New vertices of tet mesh, V prefaces VV
+  MatrixXd VV;
+  MatrixXi TT,FF,CE;
+  VectorXi P;
+  if(!mesh_with_skeleton(V,F,C,P,BE,CE,10,VV,TT,FF))
+  {
+    return false;
+  }
+  // List of boundary indices (aka fixed value indices into VV)
+  VectorXi b;
+  // List of boundary conditions of each weight function
+  MatrixXd bc;
+  if(!boundary_conditions(VV,TT,C,P,BE,CE,b,bc))
+  {
+    return false;
+  }
+
+  // compute BBW 
+  // Default bbw data and flags
+  BBWData bbw_data;
+  bbw_data.active_set_params.max_iter = 4;
+  // Weights matrix
+  if(!bbw(VV,TT,b,bc,bbw_data,W))
+  {
+    return false;
+  }
+  // Normalize weights to sum to one
+  normalize_row_sums(W,W);
+  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.obj";
+  string skel_filename = "../shared/cheburashka.tgf";
+  if(argc < 3)
+  {
+    cerr<<"Usage:"<<endl<<"    ./example input.obj"<<endl;
+    cout<<endl<<"Opening default mesh..."<<endl;
+  }else
+  {
+    // Read and prepare mesh
+    filename = argv[1];
+    skel_filename = argv[2];
+  }
+
+  // print key commands
+  cout<<"[Click] and [drag]  Rotate model using trackball."<<endl;
+  cout<<"[Z,z]               Snap rotation to canonical view."<<endl;
+  cout<<"[⌘ Z]               Undo."<<endl;
+  cout<<"[⇧ ⌘ Z]             Redo."<<endl;
+  cout<<"[^C,ESC]            Exit."<<endl;
+
+  // dirname, basename, extension and filename
+  string d,b,ext,f;
+  pathinfo(filename,d,b,ext,f);
+  // 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_faces(T,F);
+  //  }
+  }
+  if(vV.size() > 0)
+  {
+    if(!list_to_matrix(vV,V))
+    {
+      return 1;
+    }
+    triangulate(vF,F);
+  }
+
+  readTGF(skel_filename,C,BE);
+  // Recover parent indices because (C,BE) is crappy format for a tree.
+  P.resize(BE.rows(),1);
+  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;
+      }
+    }
+  }
+
+  init_weights(V,F,C,BE,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, 
+    s.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");
+  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("upright");
+  glutDisplayFunc(display);
+  glutReshapeFunc(reshape);
+  glutKeyboardFunc(key);
+  glutMouseFunc(mouse);
+  glutMotionFunc(mouse_drag);
+  glutPassiveMotionFunc((GLUTmousemotionfun)TwEventMouseMotionGLUT);
+  glutMainLoop();
+
+  return 0;
+}