slice tets and tempaltes

Former-commit-id: 35c90027fb7bdb6410b18e5b117ab578ef117abc

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/Viewer.cpp

@@ -475,8 +475,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();
 
@@ -492,7 +496,9 @@ namespace igl
     if (extension == "off" || extension =="OFF")
     {
       if (!igl::readOFF(mesh_file_name_string, data.V, data.F))
+      {
         return false;
+      }
     }
     else if (extension == "obj" || extension =="OBJ")
     {
@@ -503,7 +509,9 @@ namespace igl
       Eigen::MatrixXi UV_F;
 
       if (!(igl::readOBJ(mesh_file_name_string, data.V, data.F, corner_normals, fNormIndices, UV_V, UV_F)))
+      {
         return false;
+      }
     }
     else
     {
@@ -968,9 +976,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 (
@@ -143,7 +144,11 @@ IGL_INLINE void igl::ViewerCore::align_camera_center(
   const Eigen::MatrixXi& F)
 {
   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(