fast winding number computation of single point

include/igl/fast_winding_number.cpp

@@ -387,28 +387,45 @@ template <
   typename DerivedW>
 IGL_INLINE void igl::fast_winding_number(
   const FastWindingNumberBVH & fwn_bvh,
-  const float _accuracy_scale,
+  const float accuracy_scale,
   const Eigen::MatrixBase<DerivedQ> & Q,
   Eigen::PlainObjectBase<DerivedW> & W)
 {
-  const double accuracy_scale = 2;
   assert(Q.cols() == 3 && "Q should be 3D");
   W.resize(Q.rows(),1);
   igl::parallel_for(Q.rows(),[&](int p)
-    {
+  {
     FastWindingNumber::HDK_Sample::UT_Vector3T<float>Qp;
-          Qp[0] = Q(p,0);
-          Qp[1] = Q(p,1);
-          Qp[2] = Q(p,2);
-      W(p) = fwn_bvh.ut_solid_angle.computeSolidAngle(
-        Qp,
-        accuracy_scale) 
-        / (4.0*igl::PI);
-    },1000);
+    Qp[0] = Q(p,0);
+    Qp[1] = Q(p,1);
+    Qp[2] = Q(p,2);
+    W(p) = fwn_bvh.ut_solid_angle.computeSolidAngle(Qp,accuracy_scale) / (4.0*igl::PI);
+  },1000);
+}
+
+//TODO: there should be a function that given row id, we return solid angle
+//        however its weird without that function being private... so keeping like this.
+template <typename Derivedp>
+IGL_INLINE typename Derivedp::Scalar igl::fast_winding_number(
+  const FastWindingNumberBVH & fwn_bvh,
+  const float accuracy_scale,
+  const Eigen::MatrixBase<Derivedp> & p)
+{
+  assert(p.cols() == 3 && "p should be 3D");
+
+  FastWindingNumber::HDK_Sample::UT_Vector3T<float>Qp;
+  Qp[0] = p(0,0);
+  Qp[1] = p(0,1);
+  Qp[2] = p(0,2);
+
+  typename Derivedp::Scalar w = fwn_bvh.ut_solid_angle.computeSolidAngle(Qp,accuracy_scale) / (4.0*igl::PI);
+
+  return w;
 }
 
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
+template Eigen::Matrix<float, -1, -1, 0, -1, -1>::Scalar igl::fast_winding_number<Eigen::Matrix<float, -1, -1, 0, -1, -1> >(igl::FastWindingNumberBVH const&, float, Eigen::MatrixBase<Eigen::Matrix<float, -1, -1, 0, -1, -1> > const&);
 // generated by autoexplicit.sh
 template void igl::fast_winding_number<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, -1, -1, 0, -1, -1> const>, Eigen::Matrix<int, -1, -1, 0, -1, -1> >(Eigen::MatrixBase<Eigen::CwiseUnaryOp<Eigen::internal::scalar_cast_op<double, float>, Eigen::Matrix<double, -1, -1, 0, -1, -1> const> > const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1> > const&, int, igl::FastWindingNumberBVH&);
 // generated by autoexplicit.sh

include/igl/fast_winding_number.h

@@ -206,6 +206,19 @@ namespace igl
     const float accuracy_scale,
     const Eigen::MatrixBase<DerivedQ> & Q,
     Eigen::PlainObjectBase<DerivedW> & W);
+  // After precomputation, compute winding number at a a single point
+  //
+  // Inputs:
+  //   fwn_bvh  Precomputed bounding volume hierarchy
+  //   accuracy_scale  parameter controlling accuracy (e.g., 2)
+  //   p single position
+  // Outputs:
+  //   w  winding number of this point
+  template <typename Derivedp>
+  IGL_INLINE typename Derivedp::Scalar fast_winding_number(
+    const FastWindingNumberBVH & fwn_bvh,
+    const float accuracy_scale,
+    const Eigen::MatrixBase<Derivedp> & p);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "fast_winding_number.cpp"