accelerated 2D winding number (and sdf in example) (#2522)

* accelerated 2D winding number (and sdf in example)

* line width [ci skip]

* empty

cmake/recipes/external/libigl_tutorial_data.cmake

@@ -8,7 +8,7 @@ include(FetchContent)
 FetchContent_Declare(
     libigl_tutorial_data
     GIT_REPOSITORY https://github.com/libigl/libigl-tutorial-data
-    GIT_TAG        c1f9ede366d02e3531ecbaec5e3769312f31cccd
+    GIT_TAG        6700bf49000cd64199835fb40323e5ca9c7796ab
 )
 FetchContent_MakeAvailable(libigl_tutorial_data)
 

include/igl/box_simplices.cpp

@@ -36,6 +36,8 @@ IGL_INLINE void igl::box_simplices(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::box_simplices<Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 1, -1, 2>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>>&);
+// generated by autoexplicit.sh
 template void igl::box_simplices<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 2, 1, -1, 2>, Eigen::Matrix<double, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 2, 1, -1, 2>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&);
 // generated by autoexplicit.sh
 template void igl::box_simplices<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&);

include/igl/eytzinger_aabb.cpp

@@ -85,6 +85,8 @@ IGL_INLINE void igl::eytzinger_aabb(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::eytzinger_aabb<Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
+// generated by autoexplicit.sh
 template void igl::eytzinger_aabb<Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
 template void igl::eytzinger_aabb<Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
 #endif

include/igl/eytzinger_aabb_sdf.cpp

@@ -94,6 +94,8 @@ IGL_INLINE void igl::eytzinger_aabb_sdf(
 #ifdef IGL_STATIC_LIBRARY
 // Explicit template instantiation
 // generated by autoexplicit.sh
+template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, 1, 2, 1, 1, 2>, std::function<double (int)>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&, std::function<double (int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, 1, 2, 1, 1, 2>::Scalar&);
+// generated by autoexplicit.sh
 template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, -1, 3, 1, -1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)>, Eigen::Matrix<double, -1, 3, 1, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 1, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);
 template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, 1, 3, 1, 1, 3>, std::function<double (int)>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 3, 1, 1, 3>> const&, std::function<double (int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, 1, 3, 1, 1, 3>::Scalar&);
 template void igl::eytzinger_aabb_sdf<Eigen::Matrix<double, -1, 3, 0, -1, 3>, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)>, Eigen::Matrix<double, -1, 3, 0, -1, 3>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<double, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, std::function<double (Eigen::Matrix<double, 1, 3, 1, 1, 3> const&, int)> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 3, 0, -1, 3>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<double, -1, 1, 0, -1, 1>>&);

include/igl/eytzinger_aabb_winding_number.cpp

@@ -0,0 +1,69 @@
+#include "eytzinger_aabb_winding_number.h"
+template <
+  typename Derivedp,
+  typename DerivedV,
+  typename DerivedE,
+  typename DerivedB,
+  typename Derivedleaf,
+  typename DerivedI,
+  typename DerivedC>
+IGL_INLINE void igl::eytzinger_aabb_winding_number(
+  const Eigen::MatrixBase<Derivedp> & p,
+  const Eigen::MatrixBase<DerivedV> & V,
+  const Eigen::MatrixBase<DerivedE> & E,
+  const Eigen::MatrixBase<DerivedB> & B1,
+  const Eigen::MatrixBase<DerivedB> & B2,
+  const Eigen::MatrixBase<Derivedleaf> & leaf,
+  const Eigen::MatrixBase<DerivedI> & I,
+  const Eigen::MatrixBase<DerivedC> & C,
+  typename DerivedV::Scalar & wn)
+{
+  const auto signed_angle = [&V,&p]( const int e0, const int e1)->typename DerivedV::Scalar
+  {
+    const Eigen::Matrix<typename DerivedV::Scalar,1,2> v0 = V.row(e0) - p;
+    const Eigen::Matrix<typename DerivedV::Scalar,1,2> v1 = V.row(e1) - p;
+    const typename DerivedV::Scalar angle = std::atan2(
+      v0(0)*v1(1) - v0(1)*v1(0),
+      v0(0)*v1(0) + v0(1)*v1(1));
+    return angle;
+  };
+
+  wn = 0;
+  if(leaf.size() == 0) { wn = 0; return; }
+  // I don't think stack or queue matters
+  std::vector<int> stack;
+  stack.push_back(0);
+  while(!stack.empty())
+  {
+    int r = stack.back();
+    stack.pop_back();
+    if(leaf(r) >= 0)
+    {
+      wn += signed_angle(E(leaf(r),0),E(leaf(r),1));
+      continue;
+    }
+    // otherwise is p outside the box B1.row(r), B2.row(r)?
+    if( (p(0) < B1(r,0)) || (p(0) > B2(r,0)) ||
+        (p(1) < B1(r,1)) || (p(1) > B2(r,1)) )
+    {
+      // consider the pairs in the sequence at I(C(r)) to I(C(r)+1)
+      assert(C(r+1) - C(r) % 2 == 0);
+      for(int i = C(r); i < C(r+1); i+=2)
+      {
+        wn += signed_angle(I(i),I(i+1));
+      }
+      continue;
+    }
+    // else traverse children
+    const int left_r = 2*r + 1;
+    const int right_r = 2*r + 2;
+    stack.push_back(left_r);
+    stack.push_back(right_r);
+  }
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::eytzinger_aabb_winding_number<Eigen::Matrix<double, 1, 2, 1, 1, 2>, Eigen::Matrix<double, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<double, -1, 2, 1, -1, 2>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<double, 1, 2, 1, 1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<double, -1, 2, 1, -1, 2>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::Matrix<double, -1, -1, 0, -1, -1>::Scalar&);
+#endif
+

include/igl/eytzinger_aabb_winding_number.h

@@ -0,0 +1,51 @@
+#ifndef IGL_EYTZINGER_AABB_WINDING_NUMBER_H
+#define IGL_EYTZINGER_AABB_WINDING_NUMBER_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// Given a query point, compute its (exact) winding number with respect to a
+  /// set of segments stored in an Eytzinger AABB tree with winding number data
+  /// structure. 
+  ///
+  /// @param[in] p  #1 by 2 query point
+  /// @param[in] V  #V by 2 list of vertex positions
+  /// @param[in] E  #E by 2 list of segment endpoint vertex indices
+  /// @param[in] B1  #B by dim list of minimum corners of the Eytzinger AABBs
+  /// @param[in] B2  #B by dim list of maximum corners of the Eytzinger AABBs
+  /// @param[out] leaf #B list of leaf indices, -1 indicates internal node, -2
+  /// indicates empty node
+  /// @param[out] I  #I streamed list of vertex indices I=[(s₀,s₁), (s₂,s₃), ...]
+  /// @param[out] C  #C list of cumulative counts into I for each internal node,
+  /// so that the pairs for internal node i are found in I[C(i):C(i+1)-1]
+  /// @param[out] wn  winding number of p with respect to the segments
+  ///
+  /// \see igl::eytzinger_aabb
+  ///
+  template <
+    typename Derivedp,
+    typename DerivedV,
+    typename DerivedE,
+    typename DerivedB,
+    typename Derivedleaf,
+    typename DerivedI,
+    typename DerivedC>
+  IGL_INLINE void eytzinger_aabb_winding_number(
+    const Eigen::MatrixBase<Derivedp> & p,
+    const Eigen::MatrixBase<DerivedV> & V,
+    const Eigen::MatrixBase<DerivedE> & E,
+    const Eigen::MatrixBase<DerivedB> & B1,
+    const Eigen::MatrixBase<DerivedB> & B2,
+    const Eigen::MatrixBase<Derivedleaf> & leaf,
+    const Eigen::MatrixBase<DerivedI> & I,
+    const Eigen::MatrixBase<DerivedC> & C,
+    typename DerivedV::Scalar & wn);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "eytzinger_aabb_winding_number.cpp"
+#endif
+
+#endif
+

include/igl/eytzinger_aabb_winding_number_tree.cpp

@@ -0,0 +1,118 @@
+#include "eytzinger_aabb_winding_number_tree.h"
+#include "matlab_format.h"
+
+#include <Eigen/Sparse>
+#include <iostream>
+#include <vector>
+#include <queue>
+#include <unordered_map>
+
+template <
+  typename DerivedE,
+  typename Derivedleaf,
+  typename DerivedI,
+  typename DerivedC>
+IGL_INLINE void igl::eytzinger_aabb_winding_number_tree(
+  const Eigen::MatrixBase<DerivedE> & E,
+  const Eigen::MatrixBase<Derivedleaf> & leaf,
+  Eigen::PlainObjectBase<DerivedI> & I,
+  Eigen::PlainObjectBase<DerivedC> & C)
+{
+  const int num_vertices = E.maxCoeff() + 1;
+
+  // helper to increment or decrement count
+  const auto add_into = [](
+      std::unordered_map<int,int> & count_map,
+      const int v, const int delta)
+  {
+    auto it = count_map.find(v);
+    if(it == count_map.end())
+    {
+      count_map[v] = delta;
+    }else
+    {
+      it->second += delta;
+    }
+  };
+
+  // I claim this is still O(n log n) even though it seems like it's 
+  // recomputing each subtree O(log n) times.
+  std::vector<std::vector<int>> vI(leaf.size());
+  for(int r = 0;r<leaf.size();r++)
+  {
+    if(leaf(r) == -2) { continue; }
+    if(leaf(r) >= 0){ continue; }
+
+    // Use unordered map to count contributions
+    std::unordered_map<int,int> count_map;
+
+    std::vector<int> stack;
+    stack.push_back(r);
+    while(!stack.empty())
+    {
+      const int i = stack.back();
+      stack.pop_back();
+      if(leaf(i) >= 0)
+      {
+        // leaf node
+        add_into(count_map,E(leaf(i),0),1);
+        add_into(count_map,E(leaf(i),1),-1);
+      }else
+      {
+        const int left_i = 2*i + 1;
+        const int right_i = 2*i + 2;
+        stack.push_back(right_i);
+        stack.push_back(left_i);
+      }
+    }
+
+    std::vector<int> positive;
+    std::vector<int> negative;
+    for(const auto & it : count_map)
+    {
+      if(it.second != 0)
+      {
+        for(int c = 0;c<std::abs(it.second);c++)
+        {
+          if(it.second > 0)
+          {
+            positive.push_back(it.first);
+          }else
+          {
+            negative.push_back(it.first);
+          }
+        }
+      }
+    }
+    assert(positive.size() == negative.size());
+    std::vector<int> alternating;
+    for(int k = 0;k<positive.size();k++)
+    {
+      alternating.push_back(positive[k]);
+      alternating.push_back(negative[k]);
+    }
+    vI[r] = alternating;
+  }
+ 
+  // first collect cummulative sizes
+  C.resize(leaf.size()+1);
+  C(0) = 0;
+  for(int r = 0;r<leaf.size();r++)
+  {
+    C(r+1) = C(r) + vI[r].size();
+  }
+  I.resize(C(leaf.size()));
+  for(int r = 0;r<leaf.size();r++)
+  {
+    for(int k = 0;k<vI[r].size();k++)
+    {
+      I(C(r)+k) = vI[r][k];
+    }
+  }
+
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instantiation
+template void igl::eytzinger_aabb_winding_number_tree<Eigen::Matrix<int, -1, -1, 0, -1, -1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>, Eigen::Matrix<int, -1, 1, 0, -1, 1>>(Eigen::MatrixBase<Eigen::Matrix<int, -1, -1, 0, -1, -1>> const&, Eigen::MatrixBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>> const&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&, Eigen::PlainObjectBase<Eigen::Matrix<int, -1, 1, 0, -1, 1>>&);
+#endif

include/igl/eytzinger_aabb_winding_number_tree.h

@@ -0,0 +1,42 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2026 Alec Jacobson <[email protected]>
+// 
+// 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_EYTZINGER_AABB_WINDING_NUMBER_TREE_H
+#define IGL_EYTZINGER_AABB_WINDING_NUMBER_TREE_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+namespace igl
+{
+  /// For each internal node of the Eytzinger AABB tree, compute a streamed list
+  /// of pairs of vertices which would "close" the winding number contributions 
+  /// of all segments in the subtree rooted at that node [Jacobson et al. 2013]
+  /// 
+  /// @param[in] E  #E by 2 list of segment endpoint vertex indices
+  /// @param[out] leaf #B list of leaf indices, -1 indicates internal node, -2
+  /// indicates empty node
+  /// @param[out] I  #I streamed list of vertex indices I=[(s₀,s₁), (s₂,s₃), ...]
+  /// @param[out] C  #C list of cumulative counts into I for each internal node,
+  /// so that the pairs for internal node i are found in I[C(i):C(i+1)-1]
+  ///
+  template <
+    typename DerivedE,
+    typename Derivedleaf,
+    typename DerivedI,
+    typename DerivedC>
+  IGL_INLINE void eytzinger_aabb_winding_number_tree(
+    const Eigen::MatrixBase<DerivedE> & E,
+    const Eigen::MatrixBase<Derivedleaf> & leaf,
+    Eigen::PlainObjectBase<DerivedI> & I,
+    Eigen::PlainObjectBase<DerivedC> & C);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#include "eytzinger_aabb_winding_number_tree.cpp"
+#endif
+
+#endif

tutorial/1003_WindingNumber2D/main.cpp

@@ -0,0 +1,192 @@
+#include <igl/opengl/glfw/Viewer.h>
+#include <igl/get_seconds.h>
+#include <igl/readDMAT.h>
+#include <igl/triangulated_grid.h>
+#include <igl/box_simplices.h>
+#include <igl/eytzinger_aabb.h>
+#include <igl/eytzinger_aabb_sdf.h>
+#include <igl/eytzinger_aabb_winding_number_tree.h>
+#include <igl/eytzinger_aabb_winding_number.h>
+double sdCapsule(
+    const Eigen::RowVector2d & p, 
+    const Eigen::RowVector2d & a, 
+    const Eigen::RowVector2d & b, 
+    const double & r )
+{
+  using Scalar = double;
+  Eigen::RowVector2d pa = p - a, ba = b - a;
+  Scalar h = std::max(Scalar(0), std::min(Scalar(1), pa.dot(ba)/ba.dot(ba)));
+  return (pa - ba*h).norm() - r;
+}
+
+int main(int argc, char * argv[])
+{
+  IGL_TICTOC_LAMBDA;
+  Eigen::MatrixXd V;
+  Eigen::MatrixXi E;
+  igl::readDMAT(argc>1?argv[1]: TUTORIAL_SHARED_PATH "/libigl-cursive-V.dmat",V);
+  igl::readDMAT(argc>1?argv[2]: TUTORIAL_SHARED_PATH "/libigl-cursive-E.dmat",E);
+  Eigen::MatrixXi F(E.rows(),3);
+  F<<E,E.col(1);
+
+  Eigen::MatrixXd GV;
+  Eigen::MatrixXi GF;
+  {
+    Eigen::RowVector2d min_corner = V.colwise().minCoeff();
+    Eigen::RowVector2d max_corner = V.colwise().maxCoeff();
+    // samples on x-axis
+    const int sx = 512;
+    const int sy = sx * (max_corner(1)-min_corner(1)) / (max_corner(0)-min_corner(0));
+    Eigen::RowVector2d diag = (max_corner - min_corner).eval();
+    const double bbd = diag.norm();
+    diag /= diag.norm();
+    min_corner -= 0.1*bbd*diag;
+    max_corner += 0.1*bbd*diag;
+
+    igl::triangulated_grid(sx,sy,GV,GF);
+    // Scale and translate grid to fit bounding box
+    GV.col(0) = (GV.col(0).array() * (max_corner(0)-min_corner(0))) + min_corner(0);
+    GV.col(1) = (GV.col(1).array() * (max_corner(1)-min_corner(1))) + min_corner(1);
+  }
+
+  tictoc();
+  Eigen::Matrix<double,Eigen::Dynamic,2,Eigen::RowMajor> EB1,EB2;
+  igl::box_simplices(V,E,EB1,EB2);
+  Eigen::Matrix<double,Eigen::Dynamic,2,Eigen::RowMajor> B1,B2;
+  Eigen::VectorXi leaf;
+  igl::eytzinger_aabb(EB1,EB2,B1,B2,leaf);
+  printf("AABB build: %f s\n",tictoc());
+
+  Eigen::VectorXi I,C;
+  igl::eytzinger_aabb_winding_number_tree(E,leaf,I,C);
+  printf("Winding number tree build: %f s\n",tictoc());
+
+
+  igl::opengl::glfw::Viewer viewer;
+  viewer.data().set_mesh(GV,GF);
+
+  enum Quantity
+  {
+    SIGNED_DISTANCE,
+    UNSIGNED_DISTANCE,
+    WINDING_NUMBER
+  } quantity = SIGNED_DISTANCE;
+  bool naive = false;
+
+  const auto update = [&]()
+  {
+    tictoc();
+    Eigen::VectorXd D(GV.rows());
+    for(int i=0;i<GV.rows();i++)
+    {
+      Eigen::RowVector2d p = GV.row(i);
+      double di,wi;
+      if(naive)
+      {
+        di = std::numeric_limits<double>::infinity();
+        wi = 0;
+        for(int e = 0;e<E.rows();e++)
+        {
+          if(quantity == UNSIGNED_DISTANCE || quantity == SIGNED_DISTANCE)
+          {
+            double d = sdCapsule(
+              p,
+              V.row(E(e,0)),
+              V.row(E(e,1)),
+              0.0);
+            if(d < di){ di = d; }
+          }
+          if(quantity == WINDING_NUMBER || quantity == SIGNED_DISTANCE)
+          {
+            double angle = std::atan2(
+              (V(E(e,0),0)-p(0))*(V(E(e,1),1)-p(1)) - (V(E(e,0),1)-p(1))*(V(E(e,1),0)-p(0)),
+              (V(E(e,0),0)-p(0))*(V(E(e,1),0)-p(0)) + (V(E(e,0),1)-p(1))*(V(E(e,1),1)-p(1))
+            );
+            wi += angle;
+          }
+        }
+      }else
+      {
+        const std::function<double(const int)> primitive_p = [&](const int e)-> double
+        {
+          double d = sdCapsule(
+            p,
+            V.row(E(e,0)),
+            V.row(E(e,1)),
+            0.0);
+          return d;
+        };
+        if(quantity == UNSIGNED_DISTANCE || quantity == SIGNED_DISTANCE)
+        {
+          igl::eytzinger_aabb_sdf(p,primitive_p,B1,B2,leaf,di);
+        }
+        if(quantity == WINDING_NUMBER || quantity == SIGNED_DISTANCE)
+        {
+          igl::eytzinger_aabb_winding_number(p,V,E,B1,B2,leaf,I,C,wi);
+        }
+      }
+      switch(quantity)
+      {
+        case SIGNED_DISTANCE:
+          D(i) = di * (std::abs(wi) > 0.5?-1:1);
+          break;
+        case UNSIGNED_DISTANCE:
+          D(i) = di;
+          break;
+        case WINDING_NUMBER:
+          D(i) = wi;
+          break;
+      }
+    }
+    printf("Compute %s %s: %f s\n",
+      naive?"naive":"aabb",
+      quantity==SIGNED_DISTANCE?"signed distance":
+      quantity==UNSIGNED_DISTANCE?"unsigned distance":"winding number",
+      tictoc());
+    viewer.data().set_data(D);
+  };
+
+  update();
+  viewer.data().show_lines = false;
+  viewer.data().line_width = 2;
+  viewer.core().lighting_factor = 0;
+  viewer.data().set_edges(V,E,Eigen::RowVector3d(0,0,0));
+  // key
+  viewer.callback_key_pressed = [&](igl::opengl::glfw::Viewer &, unsigned int key, int mod)->bool
+  {
+    switch(key)
+    {
+      default:
+        return false;
+      case 'N':
+      case 'n':
+        naive = !naive;
+        break;
+      case 'W':
+      case 'w':
+        quantity = WINDING_NUMBER;
+        break;
+      case 'U':
+      case 'u':
+        quantity = UNSIGNED_DISTANCE;
+        break;
+      case 'S':
+      case 's':
+        quantity = SIGNED_DISTANCE;
+        break;
+    }
+    update();
+    return true;
+  };
+  std::cout<<R"(
+  S,s     Signed distance
+  W,w     Winding number
+  U,u     Unsigned distance
+  N,n     Toggle naive / accelerated
+  )"<<std::endl;
+ 
+
+
+  viewer.launch();
+  return 0;
+}

tutorial/CMakeLists.txt

@@ -136,4 +136,5 @@ endif()
 if(LIBIGL_TUTORIALS_CHAPTER10)
   igl_add_tutorial(1001_LipschitzOctree igl::glfw)
   igl_add_tutorial(1002_EytzingerAABB igl::glfw)
+  igl_add_tutorial(1003_WindingNumber2D igl::glfw)
 endif()