Merge branch 'alecjacobson' of https://github.com/libigl/libigl into alecjacobson

Former-commit-id: 3887cdd5ecd3f1f35654236f17c2336c0224f1cc

include/igl/ambient_occlusion.cpp

@@ -34,7 +34,6 @@ IGL_INLINE void igl::ambient_occlusion(
   const int n = P.rows();
   // Resize output
   S.resize(n,1);
-  VectorXi hits = VectorXi::Zero(n,1);
   // Embree seems to be parallel when constructing but not when tracing rays
   const MatrixXf D = random_dir_stratified(num_samples).cast<float>();
 

include/igl/copyleft/cgal/hausdorff.cpp

@@ -0,0 +1,39 @@
+#include "hausdorff.h"
+#include "../../hausdorff.h"
+#include <functional>
+
+template <
+  typename DerivedV,
+  typename Kernel,
+  typename Scalar>
+IGL_INLINE void igl::copyleft::cgal::hausdorff(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const CGAL::AABB_tree<
+    CGAL::AABB_traits<Kernel, 
+      CGAL::AABB_triangle_primitive<Kernel, 
+        typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+      >
+    >
+  > & treeB,
+  const std::vector<CGAL::Triangle_3<Kernel> > & /*TB*/,
+  Scalar & l,
+  Scalar & u)
+{
+  // Not sure why using `auto` here doesn't work with the `hausdorff` function
+  // parameter but explicitly naming the type does...
+  const std::function<double(const double &,const double &,const double &)> 
+    dist_to_B = [&treeB](
+    const double & x, const double & y, const double & z)->double
+  {
+    CGAL::Point_3<Kernel> query(x,y,z);
+    typename CGAL::AABB_tree<
+      CGAL::AABB_traits<Kernel, 
+        CGAL::AABB_triangle_primitive<Kernel, 
+          typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+        >
+      >
+    >::Point_and_primitive_id pp = treeB.closest_point_and_primitive(query);
+    return std::sqrt((query-pp.first).squared_length());
+  };
+  return igl::hausdorff(V,dist_to_B,l,u);
+}

include/igl/copyleft/cgal/hausdorff.h

@@ -0,0 +1,62 @@
+// 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_COPYLEFT_CGAL_HAUSDORFF_H
+#define IGL_COPYLEFT_CGAL_HAUSDORFF_H
+#include "../../igl_inline.h"
+
+#include <Eigen/Dense>
+#include "CGAL_includes.hpp"
+#include <vector>
+
+namespace igl 
+{
+  namespace copyleft
+  {
+    namespace cgal
+    {
+      // Compute lower and upper bounds (l,u) on the Hausdorff distance between a triangle
+      // (V) and a pointset (e.g., mesh, triangle soup) given by a distance function
+      // handle (dist_to_B).
+      //
+      // Inputs:
+      //   V   3 by 3 list of corner positions so that V.row(i) is the position of the
+      //     ith corner
+      //   treeB  CGAL's AABB tree containing triangle soup (VB,FB)
+      //   TB  list of CGAL triangles in order of FB (for determining which was found
+      //     in computation)
+      // Outputs:
+      //   l  lower bound on Hausdorff distance 
+      //   u  upper bound on Hausdorff distance
+      //
+      template <
+        typename DerivedV,
+        typename Kernel,
+        typename Scalar>
+      IGL_INLINE void hausdorff(
+        const Eigen::MatrixBase<DerivedV>& V,
+        const CGAL::AABB_tree<
+          CGAL::AABB_traits<Kernel, 
+            CGAL::AABB_triangle_primitive<Kernel, 
+              typename std::vector<CGAL::Triangle_3<Kernel> >::iterator
+            >
+          >
+        > & treeB,
+        const std::vector<CGAL::Triangle_3<Kernel> > & TB,
+        Scalar & l,
+        Scalar & u);
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "hausdorff.cpp"
+#endif
+
+#endif
+
+

include/igl/embree/shape_diameter_function.cpp

@@ -0,0 +1,69 @@
+// 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 "shape_diameter_function.h"
+#include "../shape_diameter_function.h"
+#include "EmbreeIntersector.h"
+#include "../Hit.h"
+
+template <
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::embree::shape_diameter_function(
+  const igl::embree::EmbreeIntersector & ei,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  const auto & shoot_ray = [&ei](
+    const Eigen::Vector3f& s,
+    const Eigen::Vector3f& dir)->double
+  {
+    igl::Hit hit;
+    const float tnear = 1e-4f;
+    if(ei.intersectRay(s,dir,hit,tnear))
+    {
+      return hit.t;
+    }else
+    {
+      return std::numeric_limits<double>::infinity();
+    }
+  };
+  return igl::shape_diameter_function(shoot_ray,P,N,num_samples,S);
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::embree::shape_diameter_function(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  using namespace Eigen;
+  EmbreeIntersector ei;
+  ei.init(V.template cast<float>(),F.template cast<int>());
+  shape_diameter_function(ei,P,N,num_samples,S);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::embree::shape_diameter_function<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::embree::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, 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::embree::shape_diameter_function<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(igl::embree::EmbreeIntersector const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::embree::shape_diameter_function<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<int, -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> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::embree::shape_diameter_function<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::embree::shape_diameter_function<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -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<int, -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> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&);
+#endif
+

include/igl/embree/shape_diameter_function.h

@@ -0,0 +1,60 @@
+// 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_EMBREE_SHAPE_DIAMETER_FUNCTION_H
+#define IGL_EMBREE_SHAPE_DIAMETER_FUNCTION_H
+#include "../igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+  namespace embree
+  {
+    // Forward define
+    class EmbreeIntersector;
+    // Compute shape diamter function per given point
+    //
+    // Inputs:
+    //    ei  EmbreeIntersector containing (V,F)
+    //    P  #P by 3 list of origin points
+    //    N  #P by 3 list of origin normals
+    // Outputs:
+    //    S  #P list of shape diamater function values between bounding box
+    //    diagonal (perfect sphere) and 0 (perfect needle hook)
+    //
+    template <
+      typename DerivedP,
+      typename DerivedN,
+      typename DerivedS >
+    IGL_INLINE void shape_diameter_function(
+      const EmbreeIntersector & ei,
+      const Eigen::PlainObjectBase<DerivedP> & P,
+      const Eigen::PlainObjectBase<DerivedN> & N,
+      const int num_samples,
+      Eigen::PlainObjectBase<DerivedS> & S);
+    // Wrapper which builds new EmbreeIntersector for (V,F). That's expensive so
+    // avoid this if repeatedly calling.
+    template <
+      typename DerivedV,
+      typename DerivedF,
+      typename DerivedP,
+      typename DerivedN,
+      typename DerivedS >
+    IGL_INLINE void shape_diameter_function(
+      const Eigen::PlainObjectBase<DerivedV> & V,
+      const Eigen::PlainObjectBase<DerivedF> & F,
+      const Eigen::PlainObjectBase<DerivedP> & P,
+      const Eigen::PlainObjectBase<DerivedN> & N,
+      const int num_samples,
+      Eigen::PlainObjectBase<DerivedS> & S);
+  }
+};
+#ifndef IGL_STATIC_LIBRARY
+#  include "shape_diameter_function.cpp"
+#endif
+
+#endif
+

include/igl/hausdorff.cpp

@@ -36,6 +36,54 @@ IGL_INLINE void igl::hausdorff(
   d = sqrt(std::max(dba,dab));
 }
 
+template <
+  typename DerivedV,
+  typename Scalar>
+IGL_INLINE void igl::hausdorff(
+  const Eigen::MatrixBase<DerivedV>& V,
+  const std::function<Scalar(const Scalar &,const Scalar &, const Scalar &)> & dist_to_B,
+  Scalar & l,
+  Scalar & u)
+{
+  // e  3-long vector of opposite edge lengths
+  Eigen::Matrix<typename DerivedV::Scalar,1,3> e;
+  // Maximum edge length
+  Scalar e_max = 0;
+  for(int i=0;i<3;i++)
+  {
+    e(i) = (V.row((i+1)%3)-V.row((i+2)%3)).norm();
+    e_max = std::max(e_max,e(i));
+  }
+  // Semiperimeter
+  const Scalar s = (e(0)+e(1)+e(2))*0.5;
+  // Area
+  const Scalar A = sqrt(s*(s-e(0))*(s-e(1))*(s-e(2)));
+  // Circumradius
+  const Scalar R = e(0)*e(1)*e(2)/(4.*A);
+  // inradius
+  const Scalar r = A/s;
+  // Initialize lower bound to ∞
+  l = std::numeric_limits<Scalar>::infinity();
+  // d  3-long vector of distance from each corner to B
+  Eigen::Matrix<typename DerivedV::Scalar,1,3> d;
+  Scalar u1 = std::numeric_limits<Scalar>::infinity();
+  Scalar u2 = 0;
+  for(int i=0;i<3;i++)
+  {
+    d(i) = dist_to_B(V(i,0),V(i,1),V(i,2));
+    // Lower bound is simply the max over vertex distances
+    l = std::max(d(i),l);
+    // u1 is the minimum of corner distances + maximum adjacent edge 
+    u1 = std::min(u1,d(i) + std::max(e((i+1)%3),e((i+2)%3)));
+    // u2 first takes the maximum over corner distances
+    u2 = std::max(u2,d(i));
+  }
+  // u2 is the distance from the circumcenter/midpoint of obtuse edge plus the
+  // largest corner distance 
+  u2 += (s-r>2.*R ? R : 0.5*e_max);
+  u = std::min(u1,u2);
+}
+
 #ifdef IGL_STATIC_LIBRARY
 template void igl::hausdorff<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, double>(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -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> > const&, double&);
 #endif

include/igl/hausdorff.h

@@ -10,6 +10,7 @@
 #include "igl_inline.h"
 
 #include <Eigen/Dense>
+#include <functional>
 
 namespace igl 
 {
@@ -50,6 +51,28 @@ namespace igl
     const Eigen::PlainObjectBase<DerivedVB> & VB, 
     const Eigen::PlainObjectBase<DerivedFB> & FB,
     Scalar & d);
+  // Compute lower and upper bounds (l,u) on the Hausdorff distance between a triangle
+  // (V) and a pointset (e.g., mesh, triangle soup) given by a distance function
+  // handle (dist_to_B).
+  //
+  // Inputs:
+  //   V   3 by 3 list of corner positions so that V.row(i) is the position of the
+  //     ith corner
+  //   dist_to_B  function taking the x,y,z coordinate of a query position and
+  //     outputing the closest-point distance to some point-set B
+  // Outputs:
+  //   l  lower bound on Hausdorff distance 
+  //   u  upper bound on Hausdorff distance
+  //
+  template <
+    typename DerivedV,
+    typename Scalar>
+  IGL_INLINE void hausdorff(
+    const Eigen::MatrixBase<DerivedV>& V,
+    const std::function<
+      Scalar(const Scalar &,const Scalar &, const Scalar &)> & dist_to_B,
+    Scalar & l,
+    Scalar & u);
 }
 
 #ifndef IGL_STATIC_LIBRARY

include/igl/matlab/MatlabWorkspace.h

@@ -476,15 +476,15 @@ inline bool igl::matlab::MatlabWorkspace::find(
   //cout<<name<<": "<<mxGetM(mx_data)<<" "<<mxGetN(mx_data)<<endl;
   const int m = mxGetM(mx_data);
   const int n = mxGetN(mx_data);
-
+  // TODO: It should be possible to directly load the data into the sparse
+  // matrix without going through the triplets
   // Copy data immediately
   double * pr = mxGetPr(mx_data);
   mwIndex * ir = mxGetIr(mx_data);
   mwIndex * jc = mxGetJc(mx_data);
-
   vector<Triplet<MT> > MIJV;
-  MIJV.reserve(mxGetNumberOfElements(mx_data));
-
+  const int nnz = mxGetNzmax(mx_data);
+  MIJV.reserve(nnz);
   // Iterate over outside
   int k = 0;
   for(int j=0; j<n;j++)
@@ -501,8 +501,8 @@ inline bool igl::matlab::MatlabWorkspace::find(
   }
   M.resize(m,n);
   M.setFromTriplets(MIJV.begin(),MIJV.end());
-  return true;
 
+  return true;
 }
 
 inline bool igl::matlab::MatlabWorkspace::find( 

include/igl/matlab/parse_rhs.cpp

@@ -29,6 +29,51 @@ IGL_INLINE void igl::matlab::parse_rhs_index(
   V.array() -= 1;
 }
 
+template <typename MT>
+IGL_INLINE void igl::matlab::parse_rhs(
+  const mxArray *prhs[], 
+  Eigen::SparseMatrix<MT> & M)
+{
+  using namespace Eigen;
+  using namespace std;
+  const mxArray * mx_data = prhs[0];
+  // Handle boring case where matrix is actually an empty dense matrix
+  if(mxGetNumberOfElements(mx_data) == 0)
+  {
+    M.resize(0,0);
+    return;
+  }
+  assert(mxIsSparse(mx_data));
+  assert(mxGetNumberOfDimensions(mx_data) == 2);
+  //cout<<name<<": "<<mxGetM(mx_data)<<" "<<mxGetN(mx_data)<<endl;
+  const int m = mxGetM(mx_data);
+  const int n = mxGetN(mx_data);
+  // TODO: It should be possible to directly load the data into the sparse
+  // matrix without going through the triplets
+  // Copy data immediately
+  double * pr = mxGetPr(mx_data);
+  mwIndex * ir = mxGetIr(mx_data);
+  mwIndex * jc = mxGetJc(mx_data);
+  vector<Triplet<MT> > MIJV;
+  MIJV.reserve(mxGetNumberOfElements(mx_data));
+  // Iterate over outside
+  int k = 0;
+  for(int j=0; j<n;j++)
+  {
+    // Iterate over inside
+    while(k<(int)jc[j+1])
+    {
+      //cout<<ir[k]<<" "<<j<<" "<<pr[k]<<endl;
+      assert((int)ir[k]<m);
+      assert((int)j<n);
+      MIJV.push_back(Triplet<MT >(ir[k],j,pr[k]));
+      k++;
+    }
+  }
+  M.resize(m,n);
+  M.setFromTriplets(MIJV.begin(),MIJV.end());
+}
+
 #ifdef IGL_STATIC_LIBRARY
 template void igl::matlab::parse_rhs_index<Eigen::Matrix<int, -1, 1, 0, -1, 1> >(mxArray_tag const**, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&);
 template void igl::matlab::parse_rhs_index<Eigen::Matrix<int, -1, -1, 0, -1, -1> >(mxArray_tag const**, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> >&);

include/igl/matlab/parse_rhs.h

@@ -10,6 +10,7 @@
 #include <igl/igl_inline.h>
 #include <mex.h>
 #include <Eigen/Dense>
+#include <Eigen/Sparse>
 namespace igl
 {
   namespace matlab
@@ -29,6 +30,10 @@ namespace igl
     IGL_INLINE void parse_rhs_index(
       const mxArray *prhs[], 
       Eigen::PlainObjectBase<DerivedV> & V);
+    template <typename VType>
+    IGL_INLINE void parse_rhs(
+      const mxArray *prhs[], 
+      Eigen::SparseMatrix<VType> & M);
   }
 };
 #ifndef IGL_STATIC_LIBRARY

include/igl/shape_diameter_function.cpp

@@ -0,0 +1,151 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 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 "shape_diameter_function.h"
+#include "random_dir.h"
+#include "ray_mesh_intersect.h"
+#include "EPS.h"
+#include "Hit.h"
+#include "parallel_for.h"
+#include <functional>
+#include <vector>
+#include <algorithm>
+
+template <
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::shape_diameter_function(
+  const std::function<
+    double(
+      const Eigen::Vector3f&,
+      const Eigen::Vector3f&)
+      > & shoot_ray,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  using namespace Eigen;
+  const int n = P.rows();
+  // Resize output
+  S.resize(n,1);
+  // Embree seems to be parallel when constructing but not when tracing rays
+  const MatrixXf D = random_dir_stratified(num_samples).cast<float>();
+
+  const auto & inner = [&P,&N,&num_samples,&D,&S,&shoot_ray](const int p)
+  {
+    const Vector3f origin = P.row(p).template cast<float>();
+    const Vector3f normal = N.row(p).template cast<float>();
+    int num_hits = 0;
+    double total_distance = 0;
+    for(int s = 0;s<num_samples;s++)
+    {
+      Vector3f d = D.row(s);
+      // Shoot _inward_
+      if(d.dot(normal) > 0)
+      {
+        // reverse ray
+        d *= -1;
+      }
+      const double dist = shoot_ray(origin,d);
+      if(std::isfinite(dist))
+      {
+        total_distance += dist;
+        num_hits++;
+      }
+    }
+    S(p) = total_distance/(double)num_hits;
+  };
+  parallel_for(n,inner,1000);
+}
+
+template <
+  typename DerivedV,
+  int DIM,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::shape_diameter_function(
+  const igl::AABB<DerivedV,DIM> & aabb,
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  const auto & shoot_ray = [&aabb,&V,&F](
+    const Eigen::Vector3f& _s,
+    const Eigen::Vector3f& dir)->double
+  {
+    Eigen::Vector3f s = _s+1e-4*dir;
+    igl::Hit hit;
+    if(aabb.intersect_ray(
+      V,
+      F,
+      s  .cast<typename DerivedV::Scalar>().eval(),
+      dir.cast<typename DerivedV::Scalar>().eval(),
+      hit))
+    {
+      return hit.t;
+    }else
+    {
+      return std::numeric_limits<double>::infinity();
+    }
+  };
+  return shape_diameter_function(shoot_ray,P,N,num_samples,S);
+
+}
+
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedP,
+  typename DerivedN,
+  typename DerivedS >
+IGL_INLINE void igl::shape_diameter_function(
+  const Eigen::PlainObjectBase<DerivedV> & V,
+  const Eigen::PlainObjectBase<DerivedF> & F,
+  const Eigen::PlainObjectBase<DerivedP> & P,
+  const Eigen::PlainObjectBase<DerivedN> & N,
+  const int num_samples,
+  Eigen::PlainObjectBase<DerivedS> & S)
+{
+  if(F.rows() < 100)
+  {
+    // Super naive
+    const auto & shoot_ray = [&V,&F](
+      const Eigen::Vector3f& _s,
+      const Eigen::Vector3f& dir)->double
+    {
+      Eigen::Vector3f s = _s+1e-4*dir;
+      igl::Hit hit;
+      if(ray_mesh_intersect(s,dir,V,F,hit))
+      {
+        return hit.t;
+      }else
+      {
+        return std::numeric_limits<double>::infinity();
+      }
+    };
+    return shape_diameter_function(shoot_ray,P,N,num_samples,S);
+  }
+  AABB<DerivedV,3> aabb;
+  aabb.init(V,F);
+  return shape_diameter_function(aabb,V,F,P,N,num_samples,S);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::shape_diameter_function<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<double (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, 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::shape_diameter_function<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<double (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::shape_diameter_function<Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, 1, 3, 1, 1, 3>, Eigen::Matrix<double, -1, 1, 0, -1, 1> >(std::function<double (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, Eigen::PlainObjectBase<Eigen::Matrix<double, 1, 3, 1, 1, 3> > const&, int, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1> >&);
+template void igl::shape_diameter_function<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1> >(std::function<double (Eigen::Matrix<float, 3, 1, 0, 3, 1> const&, Eigen::Matrix<float, 3, 1, 0, 3, 1> const&)> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> > const&, 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/shape_diameter_function.h

@@ -0,0 +1,95 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2017 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_SHAPE_DIAMETER_FUNCTION_H
+#define IGL_SHAPE_DIAMETER_FUNCTION_H
+#include "igl_inline.h"
+#include "AABB.h"
+#include <Eigen/Core>
+#include <functional>
+namespace igl
+{
+  // Compute shape diamater function per given point. In the parlence of the
+  // paper "Consistent Mesh Partitioning and Skeletonisation using the Shape
+  // Diameter Function" [Shapiro et al. 2008], this implementation uses a 180°
+  // cone and a _uniform_ average (_not_ a average weighted by inverse angles).
+  //
+  // Inputs:
+  //    shoot_ray  function handle that outputs hits of a given ray against a
+  //      mesh (embedded in function handles as captured variable/data)
+  //    P  #P by 3 list of origin points
+  //    N  #P by 3 list of origin normals
+  // Outputs:
+  //    S  #P list of shape diamater function values between bounding box
+  //    diagonal (perfect sphere) and 0 (perfect needle hook)
+  //
+  template <
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void shape_diameter_function(
+    const std::function<
+      double(
+        const Eigen::Vector3f&,
+        const Eigen::Vector3f&)
+        > & shoot_ray,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+  // Inputs:
+  //   AABB  axis-aligned bounding box hierarchy around (V,F)
+  template <
+    typename DerivedV,
+    int DIM,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void shape_diameter_function(
+    const igl::AABB<DerivedV,DIM> & aabb,
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+  // Inputs:
+  //    V  #V by 3 list of mesh vertex positions
+  //    F  #F by 3 list of mesh face indices into V
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedP,
+    typename DerivedN,
+    typename DerivedS >
+  IGL_INLINE void shape_diameter_function(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const Eigen::PlainObjectBase<DerivedP> & P,
+    const Eigen::PlainObjectBase<DerivedN> & N,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+  //   per_face  whether to compute per face (S is #F by 1) or per vertex (S is
+  //     #V by 1)
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedS>
+  IGL_INLINE void shape_diameter_function(
+    const Eigen::PlainObjectBase<DerivedV> & V,
+    const Eigen::PlainObjectBase<DerivedF> & F,
+    const bool per_face,
+    const int num_samples,
+    Eigen::PlainObjectBase<DerivedS> & S);
+};
+#ifndef IGL_STATIC_LIBRARY
+#  include "shape_diameter_function.cpp"
+#endif
+
+#endif
+