basic shape diameter function implementation

Former-commit-id: 44ec6b46df31036cd861e5db18820430505ca2c2

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/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
+