buggy bounds on signed distance

Former-commit-id: 4a3c70671767ff95bdd5d7d9a94ab1c32bb37039

include/igl/AABB.cpp

@@ -370,6 +370,7 @@ igl::AABB<DerivedV,DIM>::squared_distance(
   int & i,
   Eigen::PlainObjectBase<RowVectorDIMS> & c) const
 {
+  assert(max_sqr_d <= min_sqr_d);
   using namespace Eigen;
   using namespace std;
   if(max_sqr_d > min_sqr_d)
@@ -1015,8 +1016,19 @@ namespace igl
     Eigen::Matrix<float, 1, 2, 1, 1, 2> const&, 
     int&, 
     Eigen::PlainObjectBase<Eigen::Matrix<float, 1, 2, 1, 1, 2> >&) const { return -1;};
+  template <>
+  template <>
+  IGL_INLINE float igl::AABB<Eigen::Matrix<float, -1, 3, 1, -1, 3>, 2>::squared_distance(
+    Eigen::MatrixBase<Eigen::Matrix<float, -1, 3, 1, -1, 3> > const&, 
+    Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3> > const&, 
+    Eigen::Matrix<float, 1, 2, 1, 1, 2> const&, 
+    float, 
+    float, 
+    int&, 
+    Eigen::PlainObjectBase<Eigen::Matrix<float, 1, 2, 1, 1, 2> >&) const { return -1;};
 }
 
+
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 template void igl::AABB<Eigen::Matrix<double, -1, -1, 0, -1, -1>, 3>::serialize<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1> >(Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, -1, 0, -1, -1> >&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1> >&, int) const;

include/igl/copyleft/offset_surface.cpp

@@ -1,6 +1,8 @@
 #include "offset_surface.h"
 #include "marching_cubes.h"
 #include "../voxel_grid.h"
+#include "../signed_distance.h"
+#include "../flood_fill.h"
 #include <cassert>
 #include <iostream>
 
@@ -26,8 +28,9 @@ void igl::copyleft::offset_surface(
   Eigen::PlainObjectBase<Derivedside> & side,
   Eigen::PlainObjectBase<DerivedS> & S)
 {
+  typedef typename DerivedV::Scalar Scalar;
+  typedef typename DerivedF::Scalar Index;
   {
-    typedef typename DerivedV::Scalar Scalar;
     Eigen::AlignedBox<Scalar,3> box;
     typedef Eigen::Matrix<Scalar,1,3> RowVector3S;
     assert(V.cols() == 3 && "V must contain positions in 3D");
@@ -37,12 +40,19 @@ void igl::copyleft::offset_surface(
     box.extend(max_ext.transpose());
     igl::voxel_grid(box,s,1,GV,side);
   }
+
+  const Scalar h = 
+    (GV.col(0).maxCoeff()-GV.col(0).minCoeff())/((Scalar)(side(0)-1));
+  const Scalar lower_bound = isolevel-sqrt(3.0)*h;
+  const Scalar upper_bound = isolevel+sqrt(3.0)*h;
   {
-    Eigen::VectorXi I;
+    Eigen::Matrix<Index,Eigen::Dynamic,1> I;
     Eigen::Matrix<typename DerivedV::Scalar,Eigen::Dynamic,3> C,N;
     igl::signed_distance(
-      GV,V,F,signed_distance_type,S,I,C,N);
+      GV,V,F,signed_distance_type,lower_bound,upper_bound,S,I,C,N);
   }
+  igl::flood_fill(side,S);
+  
   DerivedS SS = S.array()-isolevel;
   igl::copyleft::marching_cubes(SS,GV,side(0),side(1),side(2),SV,SF);
 }

include/igl/signed_distance.cpp

@@ -28,6 +28,8 @@ IGL_INLINE void igl::signed_distance(
   const Eigen::MatrixBase<DerivedV> & V,
   const Eigen::MatrixBase<DerivedF> & F,
   const SignedDistanceType sign_type,
+  const typename DerivedV::Scalar lower_bound,
+  const typename DerivedV::Scalar upper_bound,
   Eigen::PlainObjectBase<DerivedS> & S,
   Eigen::PlainObjectBase<DerivedI> & I,
   Eigen::PlainObjectBase<DerivedC> & C,
@@ -36,6 +38,7 @@ IGL_INLINE void igl::signed_distance(
   using namespace Eigen;
   using namespace std;
   const int dim = V.cols();
+  assert(lower_bound <= upper_bound);
   assert((V.cols() == 3||V.cols() == 2) && "V should have 3d or 2d positions");
   assert((P.cols() == 3||P.cols() == 2) && "P should have 3d or 2d positions");
   assert(V.cols() == P.cols() && "V should have same dimension as P");
@@ -121,6 +124,14 @@ IGL_INLINE void igl::signed_distance(
   S.resize(P.rows(),1);
   I.resize(P.rows(),1);
   C.resize(P.rows(),dim);
+
+  // convert to bounds on (unsiged) squared distances
+  typedef typename DerivedV::Scalar Scalar; 
+  const Scalar max_abs = std::max(std::abs(lower_bound),std::abs(upper_bound));
+  const Scalar up_sqr_d = std::pow(max_abs,2.0);
+  const Scalar low_sqr_d = 
+    std::pow(std::max(max_abs-(upper_bound-lower_bound),(Scalar)0.0),2.0);
+
   parallel_for(P.rows(),[&](const int p)
   //for(int p = 0;p<P.rows();p++)
   {
@@ -148,13 +159,15 @@ IGL_INLINE void igl::signed_distance(
       case SIGNED_DISTANCE_TYPE_UNSIGNED:
         s = 1.;
         sqrd = dim==3?
-          tree3.squared_distance(V,F,q3,i,c3):
-          tree2.squared_distance(V,F,q2,i,c2);
+          tree3.squared_distance(V,F,q3,low_sqr_d,up_sqr_d,i,c3):
+          tree2.squared_distance(V,F,q2,low_sqr_d,up_sqr_d,i,c2);
         break;
       case SIGNED_DISTANCE_TYPE_DEFAULT:
       case SIGNED_DISTANCE_TYPE_WINDING_NUMBER:
         dim==3 ? 
-          signed_distance_winding_number(tree3,V,F,hier3,q3,s,sqrd,i,c3):
+          signed_distance_winding_number(
+            tree3,V,F,hier3,q3,low_sqr_d,up_sqr_d,s,sqrd,i,c3):
+          // Too lazy to add the bounds here... no reason not to
           signed_distance_winding_number(tree2,V,F,q2,s,sqrd,i,c2);
         break;
       case SIGNED_DISTANCE_TYPE_PSEUDONORMAL:
@@ -162,8 +175,10 @@ IGL_INLINE void igl::signed_distance(
         RowVector3S n3;
         Eigen::Matrix<typename DerivedV::Scalar,1,2>  n2;
         dim==3 ?
-          signed_distance_pseudonormal(tree3,V,F,FN,VN,EN,EMAP,q3,s,sqrd,i,c3,n3):
-          signed_distance_pseudonormal(tree2,V,F,FN,VN,q2,s,sqrd,i,c2,n2);
+          signed_distance_pseudonormal(
+            tree3,V,F,FN,VN,EN,EMAP,q3,low_sqr_d,up_sqr_d,s,sqrd,i,c3,n3):
+          signed_distance_pseudonormal(
+            tree2,V,F,FN,VN,q2,        low_sqr_d,up_sqr_d,s,sqrd,i,c2,n2);
         Eigen::Matrix<typename DerivedV::Scalar,1,Eigen::Dynamic>  n;
         (dim==3 ? n = n3 : n = n2);
         N.row(p) = n;
@@ -177,6 +192,30 @@ IGL_INLINE void igl::signed_distance(
   ,10000);
 }
 
+template <
+  typename DerivedP,
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedS,
+  typename DerivedI,
+  typename DerivedC,
+  typename DerivedN>
+IGL_INLINE void igl::signed_distance(
+  const Eigen::MatrixBase<DerivedP> & P,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const SignedDistanceType sign_type,
+  Eigen::PlainObjectBase<DerivedS> & S,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C,
+  Eigen::PlainObjectBase<DerivedN> & N)
+{
+  typedef typename DerivedV::Scalar Scalar;
+  Scalar lower = std::numeric_limits<Scalar>::min();
+  Scalar upper = std::numeric_limits<Scalar>::max();
+  return signed_distance(P,V,F,sign_type,lower,upper,S,I,C,N);
+}
+
 
 template <
   typename DerivedV,
@@ -271,6 +310,8 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
   const Eigen::MatrixBase<DerivedEN> & EN,
   const Eigen::MatrixBase<DerivedEMAP> & EMAP,
   const Eigen::MatrixBase<Derivedq> & q,
+  const Scalar low_sqr_d,
+  const Scalar up_sqr_d,
   Scalar & s,
   Scalar & sqrd,
   int & i,
@@ -284,10 +325,42 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
   //sqrd = tree.squared_distance(V,F,RowVector3S(q),i,(RowVector3S&)c);
   // Alec: Why was this constructor around c necessary?
   //sqrd = tree.squared_distance(V,F,q,i,(RowVector3S&)c);
-  sqrd = tree.squared_distance(V,F,q,i,c);
+  sqrd = tree.squared_distance(V,F,q,low_sqr_d,up_sqr_d,i,c);
   pseudonormal_test(V,F,FN,VN,EN,EMAP,q,i,c,s,n);
 }
 
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename DerivedFN,
+  typename DerivedVN,
+  typename DerivedEN,
+  typename DerivedEMAP,
+  typename Derivedq,
+  typename Scalar,
+  typename Derivedc,
+  typename Derivedn>
+IGL_INLINE void igl::signed_distance_pseudonormal(
+  const AABB<DerivedV,3> & tree,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const Eigen::MatrixBase<DerivedFN> & FN,
+  const Eigen::MatrixBase<DerivedVN> & VN,
+  const Eigen::MatrixBase<DerivedEN> & EN,
+  const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+  const Eigen::MatrixBase<Derivedq> & q,
+  Scalar & s,
+  Scalar & sqrd,
+  int & i,
+  Eigen::PlainObjectBase<Derivedc> & c,
+  Eigen::PlainObjectBase<Derivedn> & n)
+{
+  Scalar low_sqr_d = 0;
+  Scalar up_sqr_d = std::numeric_limits<Scalar>::max();
+  return signed_distance_pseudonormal(
+    tree,V,F,FN,VN,EN,EMAP,q,low_sqr_d,up_sqr_d,s,sqrd,i,c,n);
+}
+
 template <
   typename DerivedV,
   typename DerivedE,
@@ -304,6 +377,8 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
   const Eigen::MatrixBase<DerivedEN> & EN,
   const Eigen::MatrixBase<DerivedVN> & VN,
   const Eigen::MatrixBase<Derivedq> & q,
+  const Scalar low_sqr_d,
+  const Scalar up_sqr_d,
   Scalar & s,
   Scalar & sqrd,
   int & i,
@@ -313,10 +388,39 @@ IGL_INLINE void igl::signed_distance_pseudonormal(
   using namespace Eigen;
   using namespace std;
   typedef Eigen::Matrix<typename DerivedV::Scalar,1,2> RowVector2S;
-  sqrd = tree.squared_distance(V,E,RowVector2S(q),i,(RowVector2S&)c);
+  sqrd = tree.squared_distance(
+    V,E,RowVector2S(q),low_sqr_d,up_sqr_d,i,(RowVector2S&)c);
   pseudonormal_test(V,E,EN,VN,q,i,c,s,n);
 }
 
+template <
+  typename DerivedV,
+  typename DerivedE,
+  typename DerivedEN,
+  typename DerivedVN,
+  typename Derivedq,
+  typename Scalar,
+  typename Derivedc,
+  typename Derivedn>
+IGL_INLINE void igl::signed_distance_pseudonormal(
+  const AABB<DerivedV,2> & tree,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedE> & E,
+  const Eigen::MatrixBase<DerivedEN> & EN,
+  const Eigen::MatrixBase<DerivedVN> & VN,
+  const Eigen::MatrixBase<Derivedq> & q,
+  Scalar & s,
+  Scalar & sqrd,
+  int & i,
+  Eigen::PlainObjectBase<Derivedc> & c,
+  Eigen::PlainObjectBase<Derivedn> & n)
+{
+  Scalar low_sqr_d = 0;
+  Scalar up_sqr_d = std::numeric_limits<Scalar>::max();
+  return signed_distance_pseudonormal(
+    tree,V,E,EN,VN,q,low_sqr_d,up_sqr_d,s,sqrd,i,c,n);
+}
+
 template <
   typename DerivedV,
   typename DerivedF,
@@ -349,6 +453,8 @@ IGL_INLINE void igl::signed_distance_winding_number(
   const Eigen::MatrixBase<DerivedF> & F,
   const igl::WindingNumberAABB<Derivedq,DerivedV,DerivedF> & hier,
   const Eigen::MatrixBase<Derivedq> & q,
+  const Scalar low_sqr_d,
+  const Scalar up_sqr_d,
   Scalar & s,
   Scalar & sqrd,
   int & i,
@@ -357,11 +463,35 @@ IGL_INLINE void igl::signed_distance_winding_number(
   using namespace Eigen;
   using namespace std;
   typedef Eigen::Matrix<typename DerivedV::Scalar,1,3> RowVector3S;
-  sqrd = tree.squared_distance(V,F,RowVector3S(q),i,(RowVector3S&)c);
+  sqrd = tree.squared_distance(
+      V,F,RowVector3S(q),low_sqr_d,up_sqr_d,i,(RowVector3S&)c);
   const Scalar w = hier.winding_number(q.transpose());
   s = 1.-2.*w;
 }
 
+template <
+  typename DerivedV,
+  typename DerivedF,
+  typename Derivedq,
+  typename Scalar,
+  typename Derivedc>
+IGL_INLINE void igl::signed_distance_winding_number(
+  const AABB<DerivedV,3> & tree,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedF> & F,
+  const igl::WindingNumberAABB<Derivedq,DerivedV,DerivedF> & hier,
+  const Eigen::MatrixBase<Derivedq> & q,
+  Scalar & s,
+  Scalar & sqrd,
+  int & i,
+  Eigen::PlainObjectBase<Derivedc> & c)
+{
+  Scalar low_sqr_d = 0;
+  Scalar up_sqr_d = std::numeric_limits<Scalar>::max();
+  return signed_distance_winding_number(
+    tree,V,F,hier,q,low_sqr_d,up_sqr_d,s,sqrd,i,c);
+}
+
 template <
   typename DerivedV,
   typename DerivedF,

include/igl/signed_distance.h

@@ -32,6 +32,8 @@ namespace igl
   //   F  #F by ss list of triangle indices, ss should be 3 unless sign_type ==
   //     SIGNED_DISTANCE_TYPE_UNSIGNED
   //   sign_type  method for computing distance _sign_ S
+  //   lower_bound  lower bound of distances needed {std::numeric_limits::min}
+  //   upper_bound  lower bound of distances needed {std::numeric_limits::max}
   // Outputs:
   //   S  #P list of smallest signed distances
   //   I  #P list of facet indices corresponding to smallest distances
@@ -41,6 +43,26 @@ namespace igl
   //
   // Known bugs: This only computes distances to triangles. So unreferenced
   // vertices and degenerate triangles are ignored.
+  template <
+    typename DerivedP,
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedS,
+    typename DerivedI,
+    typename DerivedC,
+    typename DerivedN>
+  IGL_INLINE void signed_distance(
+    const Eigen::MatrixBase<DerivedP> & P,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const SignedDistanceType sign_type,
+    const typename DerivedV::Scalar lower_bound,
+    const typename DerivedV::Scalar upper_bound,
+    Eigen::PlainObjectBase<DerivedS> & S,
+    Eigen::PlainObjectBase<DerivedI> & I,
+    Eigen::PlainObjectBase<DerivedC> & C,
+    Eigen::PlainObjectBase<DerivedN> & N);
+  // Default bounds
   template <
     typename DerivedP,
     typename DerivedV,
@@ -138,11 +160,39 @@ namespace igl
     const Eigen::MatrixBase<DerivedEN> & EN,
     const Eigen::MatrixBase<DerivedEMAP> & EMAP,
     const Eigen::MatrixBase<Derivedq> & q,
+    const Scalar low_sqr_d,
+    const Scalar up_sqr_d,
     Scalar & s,
     Scalar & sqrd,
     int & i,
     Eigen::PlainObjectBase<Derivedc> & c,
     Eigen::PlainObjectBase<Derivedn> & n);
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename DerivedFN,
+    typename DerivedVN,
+    typename DerivedEN,
+    typename DerivedEMAP,
+    typename Derivedq,
+    typename Scalar,
+    typename Derivedc,
+    typename Derivedn>
+  IGL_INLINE void signed_distance_pseudonormal(
+    const AABB<DerivedV,3> & tree,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const Eigen::MatrixBase<DerivedFN> & FN,
+    const Eigen::MatrixBase<DerivedVN> & VN,
+    const Eigen::MatrixBase<DerivedEN> & EN,
+    const Eigen::MatrixBase<DerivedEMAP> & EMAP,
+    const Eigen::MatrixBase<Derivedq> & q,
+    Scalar & s,
+    Scalar & sqrd,
+    int & i,
+    Eigen::PlainObjectBase<Derivedc> & c,
+    Eigen::PlainObjectBase<Derivedn> & n);
+
   template <
     typename DerivedV,
     typename DerivedE,
@@ -159,11 +209,36 @@ namespace igl
     const Eigen::MatrixBase<DerivedEN> & EN,
     const Eigen::MatrixBase<DerivedVN> & VN,
     const Eigen::MatrixBase<Derivedq> & q,
+    const Scalar low_sqr_d,
+    const Scalar up_sqr_d,
     Scalar & s,
     Scalar & sqrd,
     int & i,
     Eigen::PlainObjectBase<Derivedc> & c,
     Eigen::PlainObjectBase<Derivedn> & n);
+
+  template <
+    typename DerivedV,
+    typename DerivedE,
+    typename DerivedEN,
+    typename DerivedVN,
+    typename Derivedq,
+    typename Scalar,
+    typename Derivedc,
+    typename Derivedn>
+  IGL_INLINE void signed_distance_pseudonormal(
+    const AABB<DerivedV,2> & tree,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedE> & E,
+    const Eigen::MatrixBase<DerivedEN> & EN,
+    const Eigen::MatrixBase<DerivedVN> & VN,
+    const Eigen::MatrixBase<Derivedq> & q,
+    Scalar & s,
+    Scalar & sqrd,
+    int & i,
+    Eigen::PlainObjectBase<Derivedc> & c,
+    Eigen::PlainObjectBase<Derivedn> & n);
+
   // Inputs:
   //   tree  AABB acceleration tree (see cgal/point_mesh_squared_distance.h)
   //   hier  Winding number evaluation hierarchy
@@ -183,6 +258,24 @@ namespace igl
   //   s  sign
   //   sqrd  squared distance
   //   pp  closest point and primitve
+  template <
+    typename DerivedV,
+    typename DerivedF,
+    typename Derivedq,
+    typename Scalar,
+    typename Derivedc>
+  IGL_INLINE void signed_distance_winding_number(
+    const AABB<DerivedV,3> & tree,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedF> & F,
+    const igl::WindingNumberAABB<Derivedq,DerivedV,DerivedF> & hier,
+    const Eigen::MatrixBase<Derivedq> & q,
+    const Scalar low_sqr_d,
+    const Scalar up_sqr_d,
+    Scalar & s,
+    Scalar & sqrd,
+    int & i,
+    Eigen::PlainObjectBase<Derivedc> & c);
   template <
     typename DerivedV,
     typename DerivedF,
@@ -199,6 +292,7 @@ namespace igl
     Scalar & sqrd,
     int & i,
     Eigen::PlainObjectBase<Derivedc> & c);
+
   template <
     typename DerivedV,
     typename DerivedF,