mv swept volume from reconfigurables

Former-commit-id: d0fdbacae95b2326711ac0fd8aff640ea93697ff

include/igl/copyleft/cgal/swept_volume.cpp

@@ -0,0 +1,58 @@
+#include "swept_volume.h"
+//#include "tictoc.h"
+#include "../../swept_volume_bounding_box.h"
+#include "../../swept_volume_signed_distance.h"
+#include "../../voxel_grid.h"
+#include "../../get_seconds.h"
+#include "../marching_cubes.h"
+#include <iostream>
+
+IGL_INLINE void igl::copyleft::cgal::swept_volume(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const std::function<Eigen::Affine3d(const double t)> & transform,
+  const size_t steps,
+  const size_t grid_res,
+  const size_t isolevel_grid,
+  Eigen::MatrixXd & SV,
+  Eigen::MatrixXi & SF)
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  using namespace igl::copyleft;
+
+  const auto & Vtransform = 
+    [&V,&transform](const size_t vi,const double t)->RowVector3d
+  {
+    Vector3d Vvi = V.row(vi).transpose();
+    return (transform(t)*Vvi).transpose();
+  };
+  AlignedBox3d Mbox;
+  swept_volume_bounding_box(V.rows(),Vtransform,steps,Mbox);
+
+  // Amount of padding: pad*h should be >= isolevel
+  const int pad = isolevel_grid+1;
+  // number of vertices on the largest side
+  const int s = grid_res+2*pad;
+  const double h = Mbox.diagonal().maxCoeff()/(double)(s-2.*pad-1.);
+  const double isolevel = isolevel_grid*h;
+
+  // create grid
+  //cerr<<"Creating grid "<<s<<" ..."<<endl;
+  //tictoc();
+  RowVector3i res;
+  MatrixXd GV;
+  voxel_grid(Mbox,s,pad,GV,res);
+  //cerr<<tictoc()<<" seconds."<<endl;
+
+  // compute values
+  VectorXd S;
+  //cerr<<"Signed Distance..."<<endl;
+  swept_volume_signed_distance(V,F,transform,steps,GV,res,h,isolevel,S);
+  //cerr<<tictoc()<<" seconds."<<endl;
+  //cerr<<"Marching cubes..."<<endl;
+  S.array()-=isolevel;
+  marching_cubes(S,GV,res(0),res(1),res(2),SV,SF);
+  //cerr<<tictoc()<<" seconds."<<endl;
+}

include/igl/copyleft/cgal/swept_volume.h

@@ -0,0 +1,44 @@
+#ifndef IGL_COPYLEFT_CGAL_SWEPT_VOLUME_H
+#define IGL_COPYLEFT_CGAL_SWEPT_VOLUME_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+namespace igl
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+    // Compute the surface of the swept volume of a solid object with surface
+    // (V,F) mesh under going rigid motion.
+    // 
+    // Inputs:
+    //   V  #V by 3 list of mesh positions in reference pose
+    //   F  #F by 3 list of mesh indices into V
+    //   transform  function handle so that transform(t) returns the rigid
+    //     transformation at time t∈[0,1]
+    //   steps  number of time steps: steps=3 --> t∈{0,0.5,1}
+    //   grid_res  number of grid cells on the longest side containing the
+    //     motion (isolevel+1 cells will also be added on each side as padding)
+    //   isolevel  distance level to be contoured as swept volume
+    // Outputs:
+    //   SV  #SV by 3 list of mesh positions of the swept surface
+    //   SF  #SF by 3 list of mesh faces into SV
+    IGL_INLINE void swept_volume(
+      const Eigen::MatrixXd & V,
+      const Eigen::MatrixXi & F,
+      const std::function<Eigen::Affine3d(const double t)> & transform,
+      const size_t steps,
+      const size_t grid_res,
+      const size_t isolevel,
+      Eigen::MatrixXd & SV,
+      Eigen::MatrixXi & SF);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "swept_volume.cpp"
+#endif
+
+#endif

include/igl/flood_fill.cpp

@@ -0,0 +1,81 @@
+#include "flood_fill.h"
+#include <limits>
+
+IGL_INLINE void igl::flood_fill(
+  const Eigen::RowVector3i & res,
+  Eigen::VectorXd & S)
+{
+  using namespace Eigen;
+  using namespace std;
+  const auto flood = [&res,&S] (
+     const int xi, 
+     const int yi, 
+     const int zi,
+     const int signed_xi, 
+     const int signed_yi, 
+     const int signed_zi,
+     const double s)
+    {
+      // flood fill this value back on this row
+      for(int bxi = xi;signed_xi<--bxi;)
+      {
+        S(bxi+res(0)*(yi + res(1)*zi)) = s;
+      }
+      // flood fill this value back on any previous rows
+      for(int byi = yi;signed_yi<--byi;)
+      {
+        for(int xi = 0;xi<res(0);xi++)
+        {
+          S(xi+res(0)*(byi + res(1)*zi)) = s;
+        }
+      }
+      // flood fill this value back on any previous "sheets"
+      for(int bzi = zi;signed_zi<--bzi;)
+      {
+        for(int yi = 0;yi<res(1);yi++)
+        {
+          for(int xi = 0;xi<res(0);xi++)
+          {
+            S(xi+res(0)*(yi + res(1)*bzi)) = s;
+          }
+        }
+      }
+    };
+  int signed_zi = -1;
+  double s = numeric_limits<double>::quiet_NaN();
+  for(int zi = 0;zi<res(2);zi++)
+  {
+    int signed_yi = -1;
+    if(zi != 0)
+    {
+      s = S(0+res(0)*(0 + res(1)*(zi-1)));
+    }
+    for(int yi = 0;yi<res(1);yi++)
+    {
+      // index of last signed item on this row
+      int signed_xi = -1;
+      if(yi != 0)
+      {
+        s = S(0+res(0)*(yi-1 + res(1)*zi));
+      }
+      for(int xi = 0;xi<res(0);xi++)
+      {
+        int i = xi+res(0)*(yi + res(1)*zi);
+        if(S(i)!=S(i))
+        {
+          if(s == s)
+          {
+            S(i) = s;
+          }
+          continue;
+        }
+        s = S(i);
+        flood(xi,yi,zi,signed_xi,signed_yi,signed_zi,s);
+        // remember this position
+        signed_xi = xi;
+        signed_yi = yi;
+        signed_zi = zi;
+      }
+    }
+  }
+}

include/igl/flood_fill.h

@@ -0,0 +1,23 @@
+#ifndef IGL_FLOOD_FILL_H
+#define IGL_FLOOD_FILL_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Given a 3D array with sparse non-nan (number?) data fill in the NaNs via
+  // flood fill. This should ensure that, e.g., if data near 0 always has a band
+  // (surface) of numbered and signed data, then components of nans will be
+  // correctly signed.
+  //
+  // Inputs:
+  //   res  3-long list of dimensions of grid
+  //   S  res(0)*res(1)*res(2)  list of scalar values (with (many) nans), see
+  //     output
+  // Outputs:
+  //   S  flood fill data in place
+  IGL_INLINE void flood_fill(const Eigen::RowVector3i& res, Eigen::VectorXd& S);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "flood_fill.cpp"
+#endif
+#endif

include/igl/grid.cpp

@@ -0,0 +1,23 @@
+#include "grid.h"
+
+IGL_INLINE void igl::grid(const Eigen::RowVector3i & res, Eigen::MatrixXd & GV)
+{
+  using namespace Eigen;
+  GV.resize(res(0)*res(1)*res(2),3);
+  for(int zi = 0;zi<res(2);zi++)
+  {
+    const auto lerp = 
+      [&](const double di, const int d)->double{return di/(double)(res(d)-1);};
+    const double z = lerp(zi,2);
+    for(int yi = 0;yi<res(1);yi++)
+    {
+      const double y = lerp(yi,1);
+      for(int xi = 0;xi<res(0);xi++)
+      {
+        const double x = lerp(xi,0);
+        GV.row(xi+res(0)*(yi + res(1)*zi)) = RowVector3d(x,y,z);
+      }
+    }
+  }
+}
+

include/igl/grid.h

@@ -0,0 +1,20 @@
+#ifndef IGL_GRID_H
+#define IGL_GRID_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  // Construct vertices of a regular grid, suitable for input to
+  // `igl::marching_cubes`
+  //
+  // Inputs:
+  //   res  3-long list of number of vertices along the x y and z dimensions
+  // Outputs:
+  //   GV  res(0)*res(1)*res(2) by 3 list of mesh vertex positions.
+  //   
+  IGL_INLINE void grid(const Eigen::RowVector3i & res, Eigen::MatrixXd & GV);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "grid.cpp"
+#endif
+#endif 

include/igl/swept_volume_bounding_box.cpp

@@ -0,0 +1,20 @@
+#include "swept_volume_bounding_box.h"
+
+IGL_INLINE void igl::swept_volume_bounding_box(
+  const size_t & n,
+  const std::function<Eigen::RowVector3d(const size_t vi, const double t)> & V,
+  const size_t & steps,
+  Eigen::AlignedBox3d & box)
+{
+  using namespace Eigen;
+  box.setEmpty();
+  const VectorXd t = VectorXd::LinSpaced(steps,0,1);
+  // Find extent over all time steps
+  for(int ti = 0;ti<t.size();ti++)
+  {
+    for(size_t vi = 0;vi<n;vi++)
+    {
+      box.extend(V(vi,t(ti)).transpose());
+    }
+  }
+}

include/igl/swept_volume_bounding_box.h

@@ -0,0 +1,31 @@
+#ifndef IGL_SWEPT_VOLUME_BOUNDING_BOX_H
+#define IGL_SWEPT_VOLUME_BOUNDING_BOX_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <functional>
+namespace igl
+{
+  // Construct an axis-aligned bounding box containing a shape undergoing a
+  // motion sampled at `steps` discrete momements.
+  //
+  // Inputs:
+  //   n  number of mesh vertices
+  //   V  function handle so that V(i,t) returns the 3d position of vertex
+  //     i at time t, for t∈[0,1]
+  //   steps  number of time steps: steps=3 --> t∈{0,0.5,1}
+  // Outputs:
+  //   box  box containing mesh under motion
+  IGL_INLINE void swept_volume_bounding_box(
+    const size_t & n,
+    const std::function<
+      Eigen::RowVector3d(const size_t vi, const double t)> & V,
+    const size_t & steps,
+    Eigen::AlignedBox3d & box);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "swept_volume_bounding_box.cpp"
+#endif
+
+#endif 

include/igl/swept_volume_signed_distance.cpp

@@ -0,0 +1,112 @@
+#include "swept_volume_signed_distance.h"
+#include "flood_fill.h"
+#include "signed_distance.h"
+#include "AABB.h"
+#include "pseudonormal_test.h"
+#include "per_face_normals.h"
+#include "per_vertex_normals.h"
+#include "per_edge_normals.h"
+#include <Eigen/Geometry>
+#include <cmath>
+
+IGL_INLINE void igl::swept_volume_signed_distance(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const std::function<Eigen::Affine3d(const double t)> & transform,
+  const size_t & steps,
+  const Eigen::MatrixXd & GV,
+  const Eigen::RowVector3i & res,
+  const double h,
+  const double isolevel,
+  const Eigen::VectorXd & S0,
+  Eigen::VectorXd & S)
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  S = S0;
+  const VectorXd t = VectorXd::LinSpaced(steps,0,1);
+  const bool finite_iso = isfinite(isolevel);
+  const double extension = (finite_iso ? isolevel : 0) + sqrt(3.0)*h;
+  Eigen::AlignedBox3d box(
+    V.colwise().minCoeff().array()-extension,
+    V.colwise().maxCoeff().array()+extension);
+  // Precomputation
+  Eigen::MatrixXd FN,VN,EN;
+  Eigen::MatrixXi E;
+  Eigen::VectorXi EMAP;
+  per_face_normals(V,F,FN);
+  per_vertex_normals(V,F,PER_VERTEX_NORMALS_WEIGHTING_TYPE_ANGLE,FN,VN);
+  per_edge_normals(
+    V,F,PER_EDGE_NORMALS_WEIGHTING_TYPE_UNIFORM,FN,EN,E,EMAP);
+  AABB<MatrixXd,3> tree;
+  tree.init(V,F);
+  for(int ti = 0;ti<t.size();ti++)
+  {
+    const Affine3d At = transform(t(ti));
+    for(int g = 0;g<GV.rows();g++)
+    {
+      // Don't bother finding out how deep inside points are.
+      if(finite_iso && S(g)==S(g) && S(g)<isolevel-sqrt(3.0)*h)
+      {
+        continue;
+      }
+      const RowVector3d gv = 
+        (GV.row(g) - At.translation().transpose())*At.linear();
+      // If outside of extended box, then consider it "far away enough"
+      if(finite_iso && !box.contains(gv.transpose()))
+      {
+        continue;
+      }
+      RowVector3d c,n;
+      int i;
+      double sqrd,s;
+      //signed_distance_pseudonormal(tree,V,F,FN,VN,EN,EMAP,gv,s,sqrd,i,c,n);
+      const double min_sqrd = 
+        finite_iso ? 
+        pow(sqrt(3.)*h+isolevel,2) : 
+        numeric_limits<double>::infinity();
+      sqrd = tree.squared_distance(V,F,gv,min_sqrd,i,c);
+      if(sqrd<min_sqrd)
+      {
+        pseudonormal_test(V,F,FN,VN,EN,EMAP,gv,i,c,s,n);
+        if(S(g) == S(g))
+        {
+          S(g) = min(S(g),s*sqrt(sqrd));
+        }else
+        {
+          S(g) = s*sqrt(sqrd);
+        }
+      }
+    }
+  }
+
+  if(finite_iso)
+  {
+    flood_fill(res,S);
+  }else
+  {
+#ifndef NDEBUG
+    for_each(S.data(),S.data()+S.size(),[](const double s){assert(s==s);});
+#endif
+  }
+}
+
+IGL_INLINE void igl::swept_volume_signed_distance(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const std::function<Eigen::Affine3d(const double t)> & transform,
+  const size_t & steps,
+  const Eigen::MatrixXd & GV,
+  const Eigen::RowVector3i & res,
+  const double h,
+  const double isolevel,
+  Eigen::VectorXd & S)
+{
+  using namespace std;
+  using namespace igl;
+  using namespace Eigen;
+  S = VectorXd::Constant(GV.rows(),1,numeric_limits<double>::quiet_NaN());
+  return 
+    swept_volume_signed_distance(V,F,transform,steps,GV,res,h,isolevel,S,S);
+}

include/igl/swept_volume_signed_distance.h

@@ -0,0 +1,57 @@
+#ifndef IGL_SWEPT_VOLUME_SIGNED_DISTANCE_H
+#define IGL_SWEPT_VOLUME_SIGNED_DISTANCE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+#include <functional>
+namespace igl
+{
+  // Compute the signed distance to a sweep surface of a mesh under-going
+  // an arbitrary motion V(t) discretely sampled at `steps`-many moments in
+  // time at a grid.
+  //
+  // Inputs:
+  //   V  #V by 3 list of mesh positions in reference pose
+  //   F  #F by 3 list of triangle indices [0,n)
+  //   transform  function handle so that transform(t) returns the rigid
+  //     transformation at time t∈[0,1]
+  //   steps  number of time steps: steps=3 --> t∈{0,0.5,1}
+  //   GV  #GV by 3 list of evaluation point grid positions
+  //   res  3-long resolution of GV grid
+  //   h  edge-length of grid
+  //   isolevel  isolevel to "focus" on; grid positions far enough away from
+  //     isolevel (based on h) will get approximate values). Set
+  //     isolevel=infinity to get good values everywhere (slow and
+  //     unnecessary if just trying to extract isolevel-level set).
+  //   S0  #GV initial values (will take minimum with these), can be same
+  //     as S)
+  // Outputs:
+  //   S  #GV list of signed distances
+  IGL_INLINE void swept_volume_signed_distance(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const std::function<Eigen::Affine3d(const double t)> & transform,
+    const size_t & steps,
+    const Eigen::MatrixXd & GV,
+    const Eigen::RowVector3i & res,
+    const double h,
+    const double isolevel,
+    const Eigen::VectorXd & S0,
+    Eigen::VectorXd & S);
+  IGL_INLINE void swept_volume_signed_distance(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const std::function<Eigen::Affine3d(const double t)> & transform,
+    const size_t & steps,
+    const Eigen::MatrixXd & GV,
+    const Eigen::RowVector3i & res,
+    const double h,
+    const double isolevel,
+    Eigen::VectorXd & S);
+  }
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "swept_volume_signed_distance.cpp"
+#endif 
+
+#endif

include/igl/voxel_grid.cpp

@@ -0,0 +1,51 @@
+#include "voxel_grid.h"
+#include "grid.h"
+
+IGL_INLINE void igl::voxel_grid(
+  const Eigen::AlignedBox3d & box, 
+  const int in_s,
+  const int pad_count,
+  Eigen::MatrixXd & GV,
+  Eigen::RowVector3i & side)
+{
+  using namespace Eigen;
+  using namespace std;
+  MatrixXd::Index si = -1;
+  box.diagonal().maxCoeff(&si);
+  //MatrixXd::Index si = 0;
+  //assert(si>=0);
+  const double s_len = box.diagonal()(si);
+  assert(in_s>(pad_count*2+1) && "s should be > 2*pad_count+1");
+  const double s = in_s - 2*pad_count;
+  side(si) = s;
+  for(int i = 0;i<3;i++)
+  {
+    if(i!=si)
+    {
+      side(i) = ceil(s * (box.max()(i)-box.min()(i))/s_len);
+    }
+  }
+  side.array() += 2*pad_count;
+  grid(side,GV);
+  // A *    p/s  + B = min
+  // A * (1-p/s) + B = max
+  // B = min - A * p/s
+  // A * (1-p/s) + min - A * p/s = max
+  // A * (1-p/s) - A * p/s = max-min
+  // A * (1-2p/s) = max-min
+  // A  = (max-min)/(1-2p/s)
+  const Array<double,3,1> ps= 
+    (double)(pad_count)/(side.transpose().cast<double>().array()-1.);
+  const Array<double,3,1> A = box.diagonal().array()/(1.0-2.*ps);
+  //// This would result in an "anamorphic", but perfectly fit grid:
+  //const Array<double,3,1> B = box.min().array() - A.array()*ps;
+  //GV.array().rowwise() *= A.transpose();
+  //GV.array().rowwise() += B.transpose();
+  // Instead scale by largest factor and move to match center
+  Array<double,3,1>::Index ai = -1;
+  double a = A.maxCoeff(&ai);
+  const Array<double,1,3> ratio = 
+    a*(side.cast<double>().array()-1.0)/(double)(side(ai)-1.0);
+  GV.array().rowwise() *= ratio;
+  GV.rowwise() += (box.center().transpose()-GV.colwise().mean()).eval();
+}

include/igl/voxel_grid.h

@@ -0,0 +1,28 @@
+#ifndef IGL_VOXEL_GRID_H
+#define IGL_VOXEL_GRID_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <Eigen/Geometry>
+namespace igl
+{
+  // Construct the cell center positions of a regular voxel grid (lattice) made
+  // of perfectly square voxels.
+  // 
+  // Inputs:
+  //   box  bounding box to enclose by grid
+  //   s  number of cell centers on largest side (including 2*pad_count)
+  //   pad_count  number of cells beyond box
+  // Outputs:
+  //   GV  res(0)*res(1)*res(2) by 3 list of cell center positions
+  //   res  3-long list of dimension of voxel grid
+  IGL_INLINE void voxel_grid(
+    const Eigen::AlignedBox3d & box, 
+    const int s,
+    const int pad_count,
+    Eigen::MatrixXd & GV,
+    Eigen::RowVector3i & side);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "voxel_grid.cpp"
+#endif
+#endif