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

Conflicts:
	include/igl/viewer/Viewer.cpp

Former-commit-id: a60ebfa8010c8933f59e3946a326690cdeb0c36a

include/igl/any_of.h

@@ -15,6 +15,8 @@ namespace igl
   // Inputs:
   //   S  matrix
   // Returns whether any entries are true
+  //
+  // Seems that Eigen (now) implements this for `Eigen::Array` 
   template <typename Mat>
   IGL_INLINE bool any_of(const Mat & S);
 }

include/igl/parula.cpp

@@ -35,7 +35,7 @@ IGL_INLINE void igl::parula(
   Eigen::PlainObjectBase<DerivedC> & C)
 {
   const double min_z = (normalize?Z.minCoeff():0);
-  const double max_z = (normalize?Z.maxCoeff():-1);
+  const double max_z = (normalize?Z.maxCoeff():1);
   return parula(Z,min_z,max_z,C);
 }
 template <typename DerivedZ, typename DerivedC>
@@ -55,4 +55,5 @@ IGL_INLINE void igl::parula(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instanciation
 template void igl::parula<double>(double, double*);
+template void igl::parula<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&, bool, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
 #endif

include/igl/slice_tets.cpp

@@ -0,0 +1,221 @@
+#include "slice_tets.h"
+#include <igl/sort.h>
+#include <igl/cat.h>
+#include <igl/per_face_normals.h>
+#include <cassert>
+#include <algorithm>
+#include <vector>
+
+template <
+  typename DerivedV, 
+  typename DerivedT, 
+  typename Derivedplane,
+  typename DerivedU,
+  typename DerivedG,
+  typename DerivedJ,
+  typename BCType>
+IGL_INLINE void igl::slice_tets(
+  const Eigen::PlainObjectBase<DerivedV>& V,
+  const Eigen::PlainObjectBase<DerivedT>& T,
+  const Eigen::PlainObjectBase<Derivedplane> & plane,
+  Eigen::PlainObjectBase<DerivedU>& U,
+  Eigen::PlainObjectBase<DerivedG>& G,
+  Eigen::PlainObjectBase<DerivedJ>& J,
+  Eigen::SparseMatrix<BCType> & BC)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(V.cols() == 3 && "V should be #V by 3");
+  assert(T.cols() == 4 && "T should be #T by 4");
+  assert(plane.size() == 4 && "Plane equation should be 4 coefficients");
+
+  // number of tets
+  const size_t m = T.rows();
+
+  typedef typename DerivedV::Scalar Scalar;
+  typedef typename DerivedT::Scalar Index;
+  typedef Matrix<Scalar,Dynamic,1> VectorXS;
+  typedef Matrix<Scalar,Dynamic,4> MatrixX4S;
+  typedef Matrix<Scalar,Dynamic,3> MatrixX3S;
+  typedef Matrix<Scalar,Dynamic,2> MatrixX2S;
+  typedef Matrix<Index,Dynamic,4> MatrixX4I;
+  typedef Matrix<Index,Dynamic,3> MatrixX3I;
+  typedef Matrix<Index,Dynamic,1> VectorXI;
+  typedef Matrix<bool,Dynamic,1> VectorXb;
+  
+  // Value of plane's implicit function at all vertices
+  VectorXS IV = 
+    (V.col(0)*plane(0) + 
+     V.col(1)*plane(1) + 
+     V.col(2)*plane(2)).array()
+    + plane(3);
+  MatrixX4S IT(m,4);
+  for(size_t t = 0;t<m;t++)
+  {
+    for(size_t c = 0;c<4;c++)
+    {
+      IT(t,c) = IV(T(t,c));
+    }
+  }
+
+  const auto & extract_rows = [](
+    const PlainObjectBase<DerivedT> & T,
+    const MatrixX4S & IT,
+    const VectorXb & I,
+    MatrixX4I  & TI,
+    MatrixX4S & ITI,
+    VectorXI & JI)
+  {
+    const Index num_I = std::count(I.data(),I.data()+I.size(),true);
+    TI.resize(num_I,4);
+    ITI.resize(num_I,4);
+    JI.resize(num_I,1);
+    {
+      size_t k = 0;
+      for(size_t t = 0;t<(size_t)T.rows();t++)
+      {
+        if(I(t))
+        {
+          TI.row(k) = T.row(t);
+          ITI.row(k) = IT.row(t);
+          JI(k) = t;
+          k++;
+        }
+      }
+      assert(k == num_I);
+    }
+  };
+
+  VectorXb I13 = (IT.array()<0).rowwise().count()==1;
+  VectorXb I31 = (IT.array()>0).rowwise().count()==1;
+  VectorXb I22 = (IT.array()<0).rowwise().count()==2;
+  MatrixX4I T13,T31,T22;
+  MatrixX4S IT13,IT31,IT22;
+  VectorXI J13,J31,J22;
+  extract_rows(T,IT,I13,T13,IT13,J13);
+  extract_rows(T,IT,I31,T31,IT31,J31);
+  extract_rows(T,IT,I22,T22,IT22,J22);
+
+  const auto & apply_sort = [] (
+     const MatrixX4I & T, 
+     const MatrixX4I & sJ, 
+     MatrixX4I & sT)
+  {
+    sT.resize(T.rows(),4);
+    for(size_t t = 0;t<(size_t)T.rows();t++)
+    {
+      for(size_t c = 0;c<4;c++)
+      {
+        sT(t,c) = T(t,sJ(t,c));
+      }
+    }
+  };
+
+  const auto & one_below = [&V,&apply_sort](
+    const MatrixX4I & T,
+    const MatrixX4S & IT,
+    MatrixX3I & G,
+    SparseMatrix<BCType> & BC)
+  {
+    // Number of tets
+    const size_t m = T.rows();
+    MatrixX4S sIT;
+    MatrixX4I sJ;
+    sort(IT,2,true,sIT,sJ);
+    MatrixX4I sT;
+    apply_sort(T,sJ,sT);
+    MatrixX3S lambda = 
+      sIT.rightCols(3).array() /
+      (sIT.rightCols(3).colwise()-sIT.col(0)).array();
+    vector<Triplet<BCType> > IJV;
+    IJV.reserve(m*3*2);
+    for(size_t c = 0;c<3;c++)
+    {
+      for(size_t t = 0;t<(size_t)m;t++)
+      {
+        IJV.push_back(Triplet<BCType>(c*m+t,  sT(t,0),  lambda(t,c)));
+        IJV.push_back(Triplet<BCType>(c*m+t,sT(t,c+1),1-lambda(t,c)));
+      }
+    }
+    BC.resize(m*3,V.rows());
+    BC.reserve(m*3*2);
+    BC.setFromTriplets(IJV.begin(),IJV.end());
+    G.resize(m,3);
+    for(size_t c = 0;c<3;c++)
+    {
+      G.col(c).setLinSpaced(m,0+c*m,(m-1)+c*m);
+    }
+  };
+
+  const auto & two_below = [&V,&apply_sort](
+    const MatrixX4I & T,
+    const MatrixX4S & IT,
+    MatrixX3I & G,
+    SparseMatrix<BCType> & BC)
+  {
+    // Number of tets
+    const size_t m = T.rows();
+    MatrixX4S sIT;
+    MatrixX4I sJ;
+    sort(IT,2,true,sIT,sJ);
+    MatrixX4I sT;
+    apply_sort(T,sJ,sT);
+    MatrixX2S lambda = 
+      sIT.rightCols(2).array() /
+      (sIT.rightCols(2).colwise()-sIT.col(0)).array();
+    MatrixX2S gamma = 
+      sIT.rightCols(2).array() /
+      (sIT.rightCols(2).colwise()-sIT.col(1)).array();
+    vector<Triplet<BCType> > IJV;
+    IJV.reserve(m*4*2);
+    for(size_t c = 0;c<2;c++)
+    {
+      for(size_t t = 0;t<(size_t)m;t++)
+      {
+        IJV.push_back(Triplet<BCType>(0*2*m+c*m+t,  sT(t,0),  lambda(t,c)));
+        IJV.push_back(Triplet<BCType>(0*2*m+c*m+t,sT(t,c+2),1-lambda(t,c)));
+        IJV.push_back(Triplet<BCType>(1*2*m+c*m+t,  sT(t,1),   gamma(t,c)));
+        IJV.push_back(Triplet<BCType>(1*2*m+c*m+t,sT(t,c+2),1- gamma(t,c)));
+      }
+    }
+    BC.resize(m*4,V.rows());
+    BC.reserve(m*4*2);
+    BC.setFromTriplets(IJV.begin(),IJV.end());
+    G.resize(2*m,3);
+    G.block(0,0,m,1) = VectorXI::LinSpaced(m,0+0*m,(m-1)+0*m);
+    G.block(0,1,m,1) = VectorXI::LinSpaced(m,0+1*m,(m-1)+1*m);
+    G.block(0,2,m,1) = VectorXI::LinSpaced(m,0+3*m,(m-1)+3*m);
+    G.block(m,0,m,1) = VectorXI::LinSpaced(m,0+0*m,(m-1)+0*m);
+    G.block(m,1,m,1) = VectorXI::LinSpaced(m,0+3*m,(m-1)+3*m);
+    G.block(m,2,m,1) = VectorXI::LinSpaced(m,0+2*m,(m-1)+2*m);
+  };
+
+  MatrixX3I G13,G31,G22;
+  SparseMatrix<BCType> BC13,BC31,BC22;
+  one_below(T13,IT13,G13,BC13);
+  one_below(T31,-IT31,G31,BC31);
+  two_below(T22,IT22,G22,BC22);
+
+  BC = cat(1,cat(1,BC13,BC31),BC22);
+  U = BC*V;
+  G.resize(G13.rows()+G31.rows()+G22.rows(),3);
+  G<<G13,(G31.array()+BC13.rows()),(G22.array()+BC13.rows()+BC31.rows());
+  MatrixX3S N;
+  per_face_normals(U,G,N);
+  Matrix<Scalar,1,3> planeN(plane(0),plane(1),plane(2));
+  VectorXb flip = (N.array().rowwise() * planeN.array()).rowwise().sum()<0;
+  for(size_t g = 0;g<(size_t)G.rows();g++)
+  {
+    if(flip(g))
+    {
+      G.row(g) = G.row(g).reverse().eval();
+    }
+  }
+
+  J.resize(G.rows());
+  J<<J13,J31,J22,J22;
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template specialization
+#endif

include/igl/slice_tets.h

@@ -0,0 +1,60 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2015 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_SLICE_TETS_H
+#define IGL_SLICE_TETS_H
+#include "igl_inline.h"
+
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+
+#include <vector>
+
+namespace igl
+{
+  // SLICE_TETS Slice through a tet mesh (V,T) along a given plane (via its
+  // implicit equation).
+  //
+  // Inputs:
+  //   V  #V by 3 list of tet mesh vertices
+  //   T  #T by 4 list of tet indices into V 
+  //   plane  list of 4 coefficients in the plane equation: [x y z 1]'*plane = 0
+  //   Optional:
+  //     'Manifold' followed by whether to stitch together triangles into a
+  //       manifold mesh {true}: results in more compact U but slightly slower.
+  // Outputs:
+  //   U  #U by 3 list of triangle mesh vertices along slice
+  //   G  #G by 3 list of triangles indices into U
+  //   J  #G list of indices into T revealing from which tet each faces comes
+  //   BC  #U by #V list of barycentric coordinates (or more generally: linear
+  //     interpolation coordinates) so that U = BC*V
+  // 
+  template <
+    typename DerivedV, 
+    typename DerivedT, 
+    typename Derivedplane,
+    typename DerivedU,
+    typename DerivedG,
+    typename DerivedJ,
+    typename BCType>
+  IGL_INLINE void slice_tets(
+    const Eigen::PlainObjectBase<DerivedV>& V,
+    const Eigen::PlainObjectBase<DerivedT>& T,
+    const Eigen::PlainObjectBase<Derivedplane> & plane,
+    Eigen::PlainObjectBase<DerivedU>& U,
+    Eigen::PlainObjectBase<DerivedG>& G,
+    Eigen::PlainObjectBase<DerivedJ>& J,
+    Eigen::SparseMatrix<BCType> & BC);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "slice_tets.cpp"
+#endif
+
+#endif
+
+

include/igl/viewer/TODOs.txt

@@ -16,6 +16,7 @@
 - zoom with pan rather than scaling
 - refresh draw while resizing
 - use constructor initializer list rather than complicated constructor
+- support per-element alpha values
 + snap to canonical view key shortcut is not working
 + resize TwBar with window
 + trackball should be able to drag over TwBar

include/igl/viewer/Viewer.cpp

@@ -468,8 +468,12 @@ namespace igl
 
     // first try to load it with a plugin
     for (unsigned int i = 0; i<plugins.size(); ++i)
+    {
       if (plugins[i]->load(mesh_file_name_string))
+      {
         return true;
+      }
+    }
 
     data.clear();
 
@@ -982,9 +986,15 @@ namespace igl
 
     opengl.init();
 
+    // Alec: It seems silly to overload launch to take a filename as an
+    // argument. load_mesh_from_file has many side effects so it makes
+    // debugging launch difficult.
+
     // Load the mesh passed as input
     if (filename.size() > 0)
+    {
       load_mesh_from_file(filename.c_str());
+    }
 
     core.align_camera_center(data.V,data.F);
 

include/igl/viewer/ViewerCore.cpp

@@ -10,6 +10,7 @@
 #include <igl/quat_to_mat.h>
 #include <igl/massmatrix.h>
 #include <Eigen/Geometry>
+#include <iostream>
 
 
 Eigen::Matrix4f lookAt (
@@ -153,7 +154,11 @@ IGL_INLINE void igl::ViewerCore::align_camera_center(
     return;
 
   get_scale_and_shift_to_fit_mesh(V,F,model_zoom,model_translation);
-  object_scale = (V.colwise().maxCoeff() - V.colwise().minCoeff()).norm();
+  // Rather than crash on empty mesh...
+  if(V.size() > 0)
+  {
+    object_scale = (V.colwise().maxCoeff() - V.colwise().minCoeff()).norm();
+  }
 }
 
 IGL_INLINE void igl::ViewerCore::get_scale_and_shift_to_fit_mesh(

tutorial/605_Tetgen/main.cpp

@@ -65,7 +65,7 @@ int main(int argc, char *argv[])
   igl::readOFF("../shared/fertility.off",V,F);
 
   // Tetrahedralize the interior
-  igl::tetrahedralize(V,F,"pq1.414", TV,TT,TF);
+  igl::tetrahedralize(V,F,"pq1.414Y", TV,TT,TF);
 
   // Compute barycenters
   igl::barycenter(TV,TT,B);

tutorial/702_WindingNumber/CMakeLists.txt

@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.6)
+project(702_WindingNumber)
+
+include("../CMakeLists.shared")
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${PROJECT_NAME}_bin ${SOURCES} ${SHARED_SOURCES})
+target_link_libraries(${PROJECT_NAME}_bin ${SHARED_LIBRARIES})

tutorial/702_WindingNumber/main.cpp

@@ -0,0 +1,139 @@
+#include <igl/readMESH.h>
+#include <igl/winding_number.h>
+#include <igl/barycenter.h>
+#include <igl/boundary_facets.h>
+#include <igl/parula.h>
+#include <igl/slice_tets.h>
+#include <igl/slice.h>
+#include <igl/viewer/Viewer.h>
+#include <Eigen/Sparse>
+#include <iostream>
+
+Eigen::MatrixXd V,BC;
+Eigen::VectorXd W;
+Eigen::MatrixXi T,F,G;
+double slice_z = 0.5;
+enum OverLayType
+{
+  OVERLAY_NONE = 0,
+  OVERLAY_INPUT = 1,
+  OVERLAY_OUTPUT = 2,
+  NUM_OVERLAY = 3,
+} overlay = OVERLAY_NONE;
+
+void update_visualization(igl::Viewer & viewer)
+{
+  using namespace Eigen;
+  using namespace std;
+  Eigen::Vector4d plane(
+    0,0,1,-((1-slice_z)*V.col(2).minCoeff()+slice_z*V.col(2).maxCoeff()));
+  MatrixXd V_vis;
+  MatrixXi F_vis;
+  VectorXi J;
+  SparseMatrix<double> bary;
+  igl::slice_tets(V,T,plane,V_vis,F_vis,J,bary);
+  VectorXd W_vis;
+  igl::slice(W,J,W_vis);
+  MatrixXd C_vis;
+  // color without normalizing
+  igl::parula(W_vis,false,C_vis);
+
+
+  const auto & append_mesh = [&C_vis,&F_vis,&V_vis](
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const RowVector3d & color)
+  {
+    F_vis.conservativeResize(F_vis.rows()+F.rows(),3);
+    F_vis.bottomRows(F.rows()) = F.array()+V_vis.rows();
+    V_vis.conservativeResize(V_vis.rows()+V.rows(),3);
+    V_vis.bottomRows(V.rows()) = V;
+    C_vis.conservativeResize(C_vis.rows()+F.rows(),3);
+    C_vis.bottomRows(F.rows()).rowwise() = color;
+  };
+  switch(overlay)
+  {
+    case OVERLAY_INPUT:
+      append_mesh(V,F,RowVector3d(1.,0.894,0.227));
+      break;
+    case OVERLAY_OUTPUT:
+      append_mesh(V,G,RowVector3d(0.8,0.8,0.8));
+      break;
+    default:
+      break;
+  }
+  viewer.data.clear();
+  viewer.data.set_mesh(V_vis,F_vis);
+  viewer.data.set_colors(C_vis);
+  viewer.data.set_face_based(true);
+}
+
+bool key_down(igl::Viewer& viewer, unsigned char key, int mod)
+{
+  switch(key)
+  {
+    default:
+      return false;
+    case ' ':
+      overlay = (OverLayType)((1+(int)overlay)%NUM_OVERLAY);
+      break;
+    case '.':
+      slice_z = std::min(slice_z+0.01,0.99);
+      break;
+    case ',':
+      slice_z = std::max(slice_z-0.01,0.01);
+      break;
+  }
+  update_visualization(viewer);
+  return true;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace Eigen;
+  using namespace std;
+
+  cout<<"Usage:"<<endl;
+  cout<<"[space]  toggle showing input mesh, output mesh or slice "<<endl;
+  cout<<"         through tet-mesh of convex hull."<<endl;
+  cout<<"'.'/','  push back/pull forward slicing plane."<<endl;
+  cout<<endl;
+
+  // Load mesh: (V,T) tet-mesh of convex hull, F contains facets of input
+  // surface mesh _after_ self-intersection resolution
+  igl::readMESH("../shared/big-sigcat.mesh",V,T,F);
+
+  // Compute barycenters of all tets
+  igl::barycenter(V,T,BC);
+
+  // Compute generalized winding number at all barycenters
+  cout<<"Computing winding number over all "<<T.rows()<<" tets..."<<endl;
+  igl::winding_number(V,F,BC,W);
+
+  // Extract interior tets
+  MatrixXi CT((W.array()>0.5).count(),4);
+  {
+    size_t k = 0;
+    for(size_t t = 0;t<T.rows();t++)
+    {
+      if(W(t)>0.5)
+      {
+        CT.row(k) = T.row(t);
+        k++;
+      }
+    }
+  }
+  // find bounary facets of interior tets
+  igl::boundary_facets(CT,G);
+  // boundary_facets seems to be reversed...
+  G = G.rowwise().reverse().eval();
+
+  // normalize
+  W = (W.array() - W.minCoeff())/(W.maxCoeff()-W.minCoeff());
+
+  // Plot the generated mesh
+  igl::Viewer viewer;
+  update_visualization(viewer);
+  viewer.callback_key_down = &key_down;
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -78,3 +78,4 @@ if(CGAL_FOUND)
 add_subdirectory("609_Boolean")
 endif()
 add_subdirectory("701_Statistics")
+add_subdirectory("702_WindingNumber")

tutorial/images/big-sigcat-winding-number.gif.REMOVED.git-id

@@ -0,0 +1 @@
+848269f54baa921cd2b5799d8cc77e3024bae9d5

tutorial/shared/big-sigcat.mesh.REMOVED.git-id

@@ -0,0 +1 @@
+0ca1d1e69f9b544b53246a4645a4a711c40de98c

tutorial/tutorial.html.REMOVED.git-id

@@ -1 +1 @@
-03ab4587d9ff6fd12ba10a0bd4df804713ded055
+2e5cd62765c38c507481d1880609faab4f97a9a9

tutorial/tutorial.md.REMOVED.git-id

@@ -1 +1 @@
-40f385aefcf815b8fa5d92f6333da0f5e793633c
+65f03843aee5fa051933483e51bc0a045058efeb