sample mesh, cumulative sum, histogram

Former-commit-id: 481aacbf65d95d9b42ab25384a93d849f58294c3

examples/ambient-occlusion-mex/example.m

@@ -7,6 +7,8 @@ M = massmatrix(V,F,'voronoi');
 Z = EV(:,7);
 qZ = round(255*(Z - min(Z))/(max(Z)-min(Z)));
 C = ind2rgb(qZ,jet(256));
+%Z = connected_components(F);
+%C = ind2rgb(Z(:),jet(max(Z)));
 
 nsp = 3;
 fs = 20;

examples/randomly-sample-mesh/Makefile

@@ -0,0 +1,54 @@
+
+.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 -liglmatlab -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
+
+MATLAB_INC=-I$(MATLAB)/extern/include/
+MATLAB_LIB=-L$(MATLAB)/bin/maci64 -lmx -lmat -lmex -lstdc++
+
+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=
+
+INC=$(LIBIGL_INC) $(ANTTWEAKBAR_INC) $(EIGEN3_INC) $(MATLAB_INC) $(GLUT_INC)
+LIB=$(OPENGL_LIB) $(GLUT_LIB) $(ANTTWEAKBAR_LIB) $(LIBIGL_LIB) $(MATLAB_LIB) $(CARBON_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/randomly-sample-mesh/example.cpp

@@ -0,0 +1,602 @@
+#include <igl/two_axis_valuator_fixed_up.h>
+#include <igl/readOBJ.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/sample_mesh.h>
+
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+
+#include <GLUT/glut.h>
+
+#include <Carbon/Carbon.h>
+
+#include <string>
+#include <vector>
+#include <stack>
+#include <iostream>
+
+
+Eigen::MatrixXd V,N,S;
+Eigen::VectorXd Vmid,Vmin,Vmax;
+double bbd = 1.0;
+Eigen::MatrixXi F;
+struct State
+{
+  igl::Camera camera;
+} s;
+
+// 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;
+  glClearColor(1,1,1,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();
+  push_object();
+
+  // 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);
+
+  
+  draw_mesh(V,F,N);
+
+  glPointSize(3.);
+  glEnable(GL_POINT_SMOOTH);
+  glColor4f(0,0.5,1,1);
+  glDisable(GL_LIGHTING);
+  glBegin(GL_POINTS);
+  for(int s = 0;s<S.rows();s++)
+  {
+    glVertex3d(S(s,0),S(s,1),S(s,2));
+  }
+  glEnd();
+
+  pop_object();
+
+  // Draw a nice floor
+  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};
+  draw_floor(GREY,DARK_GREY);
+  glPopMatrix();
+
+  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 init_samples()
+{
+  using namespace Eigen;
+  using namespace igl;
+  const int n = V.rows()*10;
+  SparseMatrix<double> B;
+  VectorXi FI;
+  sample_mesh(n,V,F,B,FI);
+  S = B*V;
+}
+
+
+KeyMap keyStates ;
+bool IS_KEYDOWN( uint16_t vKey )
+{
+  uint8_t index = vKey / 32 ;
+  uint8_t shift = vKey % 32 ;
+  return keyStates[index].bigEndianValue & (1 << shift) ;
+}
+
+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;
+      }
+  }
+  
+}
+
+int main(int argc, char * argv[])
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  string filename = "../shared/cheburashka.obj";
+  if(argc < 2)
+  {
+    cerr<<"Usage:"<<endl<<"    ./example input.obj"<<endl;
+    cout<<endl<<"Opening default mesh..."<<endl;
+  }else
+  {
+    // Read and prepare mesh
+    filename = argv[1];
+  }
+
+  // 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);
+  }
+
+  init_relative();
+  init_samples();
+
+  // 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.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;
+}

include/igl/cumsum.cpp

@@ -0,0 +1,64 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "cumsum.h"
+#include <numeric>
+#include <iostream>
+
+template <typename DerivedX, typename DerivedY>
+IGL_INLINE void igl::cumsum(
+  const Eigen::PlainObjectBase<DerivedX > & X,
+  const int dim,
+  Eigen::PlainObjectBase<DerivedY > & Y)
+{
+  using namespace Eigen;
+  using namespace std;
+  Y.resize(X.rows(),X.cols());
+  // get number of columns (or rows)
+  int num_outer = (dim == 1 ? X.cols() : X.rows() );
+  // get number of rows (or columns)
+  int num_inner = (dim == 1 ? X.rows() : X.cols() );
+  // This has been optimized so that dim = 1 or 2 is roughly the same cost.
+  // (Optimizations assume ColMajor order)
+  if(dim == 1)
+  {
+#pragma omp parallel for
+    for(int o = 0;o<num_outer;o++)
+    {
+      typename DerivedX::Scalar sum = 0;
+      for(int i = 0;i<num_inner;i++)
+      {
+        sum += X(i,o);
+        Y(i,o) = sum;
+      }
+    }
+  }else
+  {
+    for(int i = 0;i<num_inner;i++)
+    {
+      // Notice that it is *not* OK to put this above the inner loop
+      // Though here it doesn't seem to pay off...
+//#pragma omp parallel for
+      for(int o = 0;o<num_outer;o++)
+      {
+        if(i == 0)
+        {
+          Y(o,i) = X(o,i);
+        }else
+        {
+          Y(o,i) = Y(o,i-1) + X(o,i);
+        }
+      }
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::cumsum<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+template void igl::cumsum<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/cumsum.h

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_CUMSUM_H
+#define IGL_CUMSUM_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Computes a cumulative sum of the columns of X, like matlab's `cumsum`.
+  //
+  // Templates:
+  //   DerivedX  Type of matrix X
+  //   DerivedY  Type of matrix Y
+  // Inputs:
+  //   X  m by n Matrix to be cumulatively summed.
+  //   dim  dimension to take cumlative sum (1 or 2)
+  // Output:
+  //   Y  m by n Matrix containing cumulative sum.
+  //
+  template <typename DerivedX, typename DerivedY>
+  IGL_INLINE void cumsum(
+    const Eigen::PlainObjectBase<DerivedX > & X,
+    const int dim,
+    Eigen::PlainObjectBase<DerivedY > & Y);
+  //template <typename DerivedX, typename DerivedY>
+  //IGL_INLINE void cumsum(
+  //  const Eigen::PlainObjectBase<DerivedX > & X,
+  //  Eigen::PlainObjectBase<DerivedY > & Y);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "cumsum.cpp"
+#endif
+
+#endif
+

include/igl/histc.cpp

@@ -0,0 +1,66 @@
+#include "histc.h"
+#include <cassert>
+
+template <typename DerivedX, typename DerivedE, typename DerivedN, typename DerivedB>
+IGL_INLINE void igl::histc(
+  const Eigen::PlainObjectBase<DerivedX > & X,
+  const Eigen::PlainObjectBase<DerivedE > & E,
+  Eigen::PlainObjectBase<DerivedN > & N,
+  Eigen::PlainObjectBase<DerivedB > & B)
+{
+  const int m = X.size();
+  const int n = E.size();
+  assert( 
+    (E.topLeftCorner(E.size()-1) - 
+      E.bottomRightCorner(E.size()-1)).maxCoeff() >= 0 && 
+    "E should be monotonically increasing");
+  N.resize(n,1);
+  B.resize(m,1);
+  N.setConstant(0);
+#pragma omp parallel for
+  for(int j = 0;j<m;j++)
+  {
+    const double x = X(j);
+    // Boring one-offs
+    if(x < E(0) || x > E(E.size()-1))
+    {
+      B(j) = -1;
+      continue;
+    }
+    // Find x in E
+    int l = 0;
+    int h = E.size()-1;
+    int k = l;
+    while((h-l)>1)
+    {
+      assert(x >= E(l));
+      assert(x <= E(h));
+      k = (h+l)/2;
+      if(x < E(k))
+      {
+        h = k;
+      }else
+      {
+        l = k;
+      }
+    }
+    if(x == E(h))
+    {
+      k = h;
+    }else
+    {
+      k = l;
+    }
+    B(j) = k;
+    // Not sure if this pragma is needed
+#pragma omp atomic
+    N(k)++;
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::histc<Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::histc<Eigen::Matrix<double, -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<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::histc<Eigen::Matrix<double, -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<double, -1, 1, 0, -1, 1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/histc.h

@@ -0,0 +1,43 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_histc_H
+#define IGL_histc_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Like matlab's histc. Count occurances of values in X between consecutive
+  // entries in E
+  //
+  // Inputs:
+  //   X  m-long Vector of values
+  //   E  n-long Monotonically increasing vector of edges
+  // Outputs:
+  //   N  n-long vector where N(k) reveals how many values in X fall between
+  //     E(k) <= X < E(k+1)
+  //   B  m-long vector of bin ids so that B(j) = k if E(k) <= X(j) < E(k+1).
+  //     B(j) = -1 if X(j) is outside of E.
+  //
+  template <typename DerivedX, typename DerivedE, typename DerivedN, typename DerivedB>
+  IGL_INLINE void histc(
+    const Eigen::PlainObjectBase<DerivedX > & X,
+    const Eigen::PlainObjectBase<DerivedE > & E,
+    Eigen::PlainObjectBase<DerivedN > & N,
+    Eigen::PlainObjectBase<DerivedB > & B);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "histc.cpp"
+#endif
+
+#endif
+
+
+

include/igl/sample_mesh.cpp

@@ -0,0 +1,74 @@
+#include "sample_mesh.h"
+#include "doublearea.h"
+#include "cumsum.h"
+#include "histc.h"
+#include <cassert>
+
+template <typename DerivedV, typename DerivedF, typename DerivedB, typename DerivedFI>
+IGL_INLINE void igl::sample_mesh(
+  const int n,
+  const Eigen::PlainObjectBase<DerivedV > & V,
+  const Eigen::PlainObjectBase<DerivedF > & F,
+  Eigen::PlainObjectBase<DerivedB > & B,
+  Eigen::PlainObjectBase<DerivedFI > & FI)
+{
+  using namespace Eigen;
+  typedef typename DerivedV::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,1> VectorXs;
+  VectorXs A;
+  doublearea(V,F,A);
+  A /= A.array().sum();
+  // Should be traingle mesh. Although Turk's method 1 generalizes...
+  assert(F.cols() == 3);
+  VectorXs C;
+  VectorXs A0(A.size()+1);
+  A0(0) = 0;
+  A0.bottomRightCorner(A.size(),1) = A;
+  cumsum(A0,1,C);
+  const VectorXs R = (VectorXs::Random(n,1).array() + 1.)/2.;
+  assert(R.minCoeff() >= 0);
+  assert(R.maxCoeff() <= 1);
+  VectorXi _;
+  histc(R,C,_,FI);
+  const VectorXs S = (VectorXs::Random(n,1).array() + 1.)/2.;
+  const VectorXs T = (VectorXs::Random(n,1).array() + 1.)/2.;
+  B.resize(n,3);
+  B.col(0) = 1.-T.array().sqrt();
+  B.col(1) = (1.-S.array()) * T.array().sqrt();
+  B.col(2) = S.array() * T.array().sqrt();
+}
+
+template <typename DerivedV, typename DerivedF, typename ScalarB, typename DerivedFI>
+IGL_INLINE void igl::sample_mesh(
+  const int n,
+  const Eigen::PlainObjectBase<DerivedV > & V,
+  const Eigen::PlainObjectBase<DerivedF > & F,
+  Eigen::SparseMatrix<ScalarB > & B,
+  Eigen::PlainObjectBase<DerivedFI > & FI)
+{
+  using namespace Eigen;
+  using namespace std;
+  Matrix<ScalarB,Dynamic,3> BC;
+  sample_mesh(n,V,F,BC,FI);
+  vector<Triplet<ScalarB> > BIJV;
+  BIJV.reserve(n*3);
+  for(int s = 0;s<n;s++)
+  {
+    for(int c = 0;c<3;c++)
+    {
+      assert(FI(s) < F.rows());
+      assert(FI(s) >= 0);
+      const int v = F(FI(s),c);
+      BIJV.push_back(Triplet<ScalarB>(s,v,BC(s,c)));
+    }
+  }
+  B.resize(n,V.rows());
+  B.reserve(n*3);
+  B.setFromTriplets(BIJV.begin(),BIJV.end());
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::sample_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);
+template void igl::sample_mesh<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, Eigen::SparseMatrix<double, 0, int>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
+#endif

include/igl/sample_mesh.h

@@ -0,0 +1,50 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_SAMPLE_MESH_H
+#define IGL_SAMPLE_MESH_H
+
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Sparse>
+
+namespace igl
+{
+  // SAMPLE_MESH Randomly sample a mesh (V,F) n times.
+  //
+  // Inputs:
+  //   n  number of samples
+  //   V  #V by dim list of mesh vertex positions
+  //   F  #F by 3 list of mesh triangle indices
+  // Outputs:
+  //   B  n by #V list of barycentric coordinates such that each row i has
+  //     three nonzero entries hoding barycentric corridinates of the point.
+  //   FI  n list of indices into F 
+  //
+  template <typename DerivedV, typename DerivedF, typename DerivedB, typename DerivedFI>
+  IGL_INLINE void sample_mesh(
+    const int n,
+    const Eigen::PlainObjectBase<DerivedV > & V,
+    const Eigen::PlainObjectBase<DerivedF > & F,
+    Eigen::PlainObjectBase<DerivedB > & B,
+    Eigen::PlainObjectBase<DerivedFI > & FI);
+  template <typename DerivedV, typename DerivedF, typename ScalarB, typename DerivedFI>
+  IGL_INLINE void sample_mesh(
+    const int n,
+    const Eigen::PlainObjectBase<DerivedV > & V,
+    const Eigen::PlainObjectBase<DerivedF > & F,
+    Eigen::SparseMatrix<ScalarB > & B,
+    Eigen::PlainObjectBase<DerivedFI > & FI);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "sample_mesh.cpp"
+#endif
+
+#endif
+
+