mv most winding number code from project

Former-commit-id: 1177af4015d294d9e16227698dcc4c9b86504c39

include/igl/WindingNumberAABB.h

@@ -0,0 +1,340 @@
+#ifndef IGL_WINDINGNUMBERAABB_H
+#define IGL_WINDINGNUMBERAABB_H
+#include "WindingNumberTree.h"
+
+namespace igl
+{
+  template <typename Point>
+  class WindingNumberAABB : public WindingNumberTree<Point>
+  {
+    protected:
+      Point min_corner;
+      Point max_corner;
+      double total_positive_area;
+    public: 
+      enum SplitMethod
+      {
+        CENTER_ON_LONGEST_AXIS = 0,
+        MEDIAN_ON_LONGEST_AXIS = 1,
+        NUM_SPLIT_METHODS = 3
+      } split_method;
+    public:
+      inline WindingNumberAABB(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F);
+      inline WindingNumberAABB(
+        const WindingNumberTree<Point> & parent,
+        const Eigen::MatrixXi & F);
+      // Initialize some things
+      inline void init();
+      inline bool inside(const Point & p) const;
+      inline virtual void grow();
+      // Compute min and max corners
+      inline void compute_min_max_corners();
+      inline double max_abs_winding_number(const Point & p) const;
+      inline double max_simple_abs_winding_number(const Point & p) const;
+  };
+}
+
+// Implementation
+
+#include "winding_number.h"
+
+#include "barycenter.h"
+#include "median.h"
+#include "doublearea.h"
+#include "per_face_normals.h"
+
+#include <limits>
+#include <vector>
+#include <iostream>
+
+// Minimum number of faces in a hierarchy element (this is probably dependent
+// on speed of machine and compiler optimization)
+#ifndef WindingNumberAABB_MIN_F
+#  define WindingNumberAABB_MIN_F 100
+#endif
+
+template <typename Point>
+inline void igl::WindingNumberAABB<Point>::init()
+{
+  using namespace igl;
+  using namespace Eigen;
+  assert(max_corner.size() == 3);
+  assert(min_corner.size() == 3);
+  compute_min_max_corners();
+  VectorXd dblA;
+  doublearea(this->getV(),this->getF(),dblA);
+  total_positive_area = dblA.sum()/2.0;
+}
+
+template <typename Point>
+inline igl::WindingNumberAABB<Point>::WindingNumberAABB(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F):
+  WindingNumberTree<Point>(V,F),
+  min_corner(),
+  max_corner(),
+  total_positive_area(std::numeric_limits<double>::infinity()),
+  split_method(MEDIAN_ON_LONGEST_AXIS)
+{
+  init();
+}
+
+template <typename Point>
+inline igl::WindingNumberAABB<Point>::WindingNumberAABB(
+  const WindingNumberTree<Point> & parent,
+  const Eigen::MatrixXi & F):
+  WindingNumberTree<Point>(parent,F),
+  min_corner(),
+  max_corner(),
+  total_positive_area(std::numeric_limits<double>::infinity()),
+  split_method(MEDIAN_ON_LONGEST_AXIS)
+{
+  init();
+}
+
+template <typename Point>
+inline void igl::WindingNumberAABB<Point>::grow()
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  //cout<<"cap.rows(): "<<this->getcap().rows()<<endl;
+  //cout<<"F.rows(): "<<this->getF().rows()<<endl;
+
+  // Base cases
+  if(
+    this->getF().rows() <= (WindingNumberAABB_MIN_F>0?WindingNumberAABB_MIN_F:0) ||
+    (this->getcap().rows() - 2) >= this->getF().rows())
+  {
+    // Don't grow
+    return;
+  }
+
+  // Compute longest direction
+  int max_d = -1;
+  double max_len = -numeric_limits<double>::infinity();
+  for(int d = 0;d<min_corner.size();d++)
+  {
+    if( (max_corner[d] - min_corner[d]) > max_len )
+    {
+      max_len = (max_corner[d] - min_corner[d]);
+      max_d = d;
+    }
+  }
+  // Compute facet barycenters
+  MatrixXd BC;
+  barycenter(this->getV(),this->getF(),BC);
+
+
+  // Blerg, why is selecting rows so difficult
+
+  vector<int> id( this->getF().rows());
+  double split_value;
+  // Split in longest direction
+  switch(split_method)
+  {
+    case MEDIAN_ON_LONGEST_AXIS:
+      // Determine median
+      median(BC.col(max_d),split_value);
+      break;
+    default:
+      assert(false);
+    case CENTER_ON_LONGEST_AXIS:
+      split_value = 0.5*(max_corner[max_d] + min_corner[max_d]);
+      break;
+  }
+  //cout<<"c: "<<0.5*(max_corner[max_d] + min_corner[max_d])<<" "<<
+  //  "m: "<<split_value<<endl;;
+
+  for(int i = 0;i<this->getF().rows();i++)
+  {
+    if(BC(i,max_d) <= split_value)
+    {
+      id[i] = 0; //left
+    }else
+    {
+      id[i] = 1; //right
+    }
+  }
+
+  const int lefts = (int) count(id.begin(),id.end(),0);
+  const int rights = (int) count(id.begin(),id.end(),1);
+  if(lefts == 0 || rights == 0)
+  {
+    // badly balanced base case (could try to recut)
+    return;
+  }
+  assert(lefts+rights == this->getF().rows());
+  MatrixXi leftF(lefts,  this->getF().cols());
+  MatrixXi rightF(rights,this->getF().cols());
+  int left_i = 0;
+  int right_i = 0;
+  for(int i = 0;i<this->getF().rows();i++)
+  {
+    if(id[i] == 0)
+    {
+      leftF.row(left_i++) = this->getF().row(i);
+    }else if(id[i] == 1)
+    {
+      rightF.row(right_i++) = this->getF().row(i);
+    }else
+    {
+      assert(false);
+    }
+  }
+  assert(right_i == rightF.rows());
+  assert(left_i == leftF.rows());
+  // Finally actually grow children and Recursively grow
+  WindingNumberAABB<Point> * leftWindingNumberAABB = new WindingNumberAABB<Point>(*this,leftF);
+  leftWindingNumberAABB->grow();
+  this->children.push_back(leftWindingNumberAABB);
+  WindingNumberAABB<Point> * rightWindingNumberAABB = new WindingNumberAABB<Point>(*this,rightF);
+  rightWindingNumberAABB->grow();
+  this->children.push_back(rightWindingNumberAABB);
+}
+
+template <typename Point>
+inline bool igl::WindingNumberAABB<Point>::inside(const Point & p) const
+{
+  assert(p.size() == max_corner.size());
+  assert(p.size() == min_corner.size());
+  for(int i = 0;i<p.size();i++)
+  {
+    if( p(i) < min_corner(i) || p(i) >= max_corner(i))
+    {
+      return false;
+    }
+  }
+  return true;
+}
+
+template <typename Point>
+inline void igl::WindingNumberAABB<Point>::compute_min_max_corners()
+{
+  using namespace std;
+  // initialize corners
+  for(int d = 0;d<min_corner.size();d++)
+  {
+    min_corner[d] =  numeric_limits<double>::infinity();
+    max_corner[d] = -numeric_limits<double>::infinity();
+  }
+
+  this->center = Point(0,0,0);
+  // Loop over facets
+  for(int i = 0;i<this->getF().rows();i++)
+  {
+    for(int j = 0;j<this->getF().cols();j++)
+    {
+      for(int d = 0;d<min_corner.size();d++)
+      {
+        min_corner[d] = 
+          this->getV()(this->getF()(i,j),d) < min_corner[d] ?  
+            this->getV()(this->getF()(i,j),d) : min_corner[d];
+        max_corner[d] = 
+          this->getV()(this->getF()(i,j),d) > max_corner[d] ?  
+            this->getV()(this->getF()(i,j),d) : max_corner[d];
+      }
+      // This is biased toward vertices incident on more than one face, but
+      // perhaps that's good
+      this->center += this->getV().row(this->getF()(i,j));
+    }
+  }
+  // Average
+  this->center.array() /= this->getF().size();
+
+  //cout<<"min_corner: "<<this->min_corner.transpose()<<endl;
+  //cout<<"Center: "<<this->center.transpose()<<endl;
+  //cout<<"max_corner: "<<this->max_corner.transpose()<<endl;
+  //cout<<"Diag center: "<<((this->max_corner + this->min_corner)*0.5).transpose()<<endl;
+  //cout<<endl;
+
+  this->radius = (max_corner-min_corner).norm()/2.0;
+}
+
+template <typename Point>
+inline double igl::WindingNumberAABB<Point>::max_abs_winding_number(const Point & p) const
+{
+  using namespace std;
+  // Only valid if not inside
+  if(inside(p))
+  {
+    return numeric_limits<double>::infinity();
+  }
+  // Q: we know the total positive area so what's the most this could project
+  // to? Remember it could be layered in the same direction.
+  return numeric_limits<double>::infinity();
+}
+
+template <typename Point>
+inline double igl::WindingNumberAABB<Point>::max_simple_abs_winding_number(const Point & p) const
+{
+  using namespace std;
+  using namespace Eigen;
+  using namespace igl;
+  // Only valid if not inside
+  if(inside(p))
+  {
+    return numeric_limits<double>::infinity();
+  }
+  // Max simple is the same as sum of positive winding number contributions of
+  // bounding box
+
+  // begin precomputation
+  //MatrixXd BV((int)pow(2,3),3);
+  MatrixXd BV((int)(1<<3),3);
+  BV <<
+    min_corner[0],min_corner[1],min_corner[2],
+    min_corner[0],min_corner[1],max_corner[2],
+    min_corner[0],max_corner[1],min_corner[2],
+    min_corner[0],max_corner[1],max_corner[2],
+    max_corner[0],min_corner[1],min_corner[2],
+    max_corner[0],min_corner[1],max_corner[2],
+    max_corner[0],max_corner[1],min_corner[2],
+    max_corner[0],max_corner[1],max_corner[2];
+  MatrixXi BF(2*2*3,3);
+  BF <<
+    0,6,4,
+    0,2,6,
+    0,3,2,
+    0,1,3,
+    2,7,6,
+    2,3,7,
+    4,6,7,
+    4,7,5,
+    0,4,5,
+    0,5,1,
+    1,5,7,
+    1,7,3;
+  MatrixXd BFN;
+  per_face_normals(BV,BF,BFN);
+  // end of precomputation
+
+  // Only keep those with positive dot products
+  MatrixXi PBF(BF.rows(),BF.cols());
+  int pbfi = 0;
+  Point p2c = 0.5*(min_corner+max_corner)-p;
+  for(int i = 0;i<BFN.rows();i++)
+  {
+    if(p2c.dot(BFN.row(i)) > 0)
+    {
+      PBF.row(pbfi++) = BF.row(i);
+    }
+  }
+  PBF.conservativeResize(pbfi,PBF.cols());
+  double w = numeric_limits<double>::infinity();
+  igl::winding_number_3(
+    BV.data(),
+    BV.rows(),
+    PBF.data(),
+    PBF.rows(),
+    p.data(),
+    1,
+    &w);
+  return w;
+}
+
+//// Explicit instanciation
+//template class igl::WindingNumberAABB<Eigen::Vector3d >;
+#endif

include/igl/WindingNumberMethod.h

@@ -0,0 +1,13 @@
+#ifndef IGL_WINDINGNUMBERMETHOD_H
+#define IGL_WINDINGNUMBERMETHOD_H
+namespace igl
+{
+  enum WindingNumberMethod
+  {
+    EXACT_WINDING_NUMBER_METHOD = 0, // Exact hierarchical evaluation
+    APPROX_SIMPLE_WINDING_NUMBER_METHOD = 1,
+    APPROX_CACHE_WINDING_NUMBER_METHOD = 2, // Approximate hierarchical evaluation
+    NUM_WINDING_NUMBER_METHODS = 3
+  };
+}
+#endif

include/igl/WindingNumberTree.h

@@ -0,0 +1,454 @@
+#ifndef IGL_BOUNDINGTREE_H
+#define IGL_BOUNDINGTREE_H
+#include <list>
+#include <map>
+#include <Eigen/Dense>
+#include "WindingNumberMethod.h"
+
+namespace igl
+{
+  // Templates:
+  //   Point  type for points in space, e.g. Eigen::Vector3d
+  template <typename Point>
+  class WindingNumberTree
+  {
+    public:
+      // Method to use (see enum above)
+      //static double min_max_w;
+      static std::map< 
+        std::pair<const WindingNumberTree*,const WindingNumberTree*>, double>
+          cached;
+    protected:
+      WindingNumberMethod method;
+      const WindingNumberTree * parent;
+      std::list<WindingNumberTree * > children;
+      //// List of boundary edges (recall edges are vertices in 2d)
+      //const Eigen::MatrixXi boundary;
+      // Base mesh vertices
+      Eigen::MatrixXd & V;
+      // Base mesh vertices with duplicates removed
+      Eigen::MatrixXd SV;
+      // Facets in this bounding volume
+      Eigen::MatrixXi F;
+      // Tesselated boundary curve
+      Eigen::MatrixXi cap;
+      // Upper Bound on radius of enclosing ball
+      double radius;
+      // (Approximate) center (of mass)
+      Point center;
+    public:
+      // For root
+      inline WindingNumberTree(
+        const Eigen::MatrixXd & V,
+        const Eigen::MatrixXi & F);
+      // For chilluns 
+      inline WindingNumberTree(
+        const WindingNumberTree<Point> & parent,
+        const Eigen::MatrixXi & F);
+      inline virtual ~WindingNumberTree();
+      // Set method
+      inline void set_method( const WindingNumberMethod & m);
+    public:
+      inline const Eigen::MatrixXd & getV() const;
+      inline const Eigen::MatrixXi & getF() const;
+      inline const Eigen::MatrixXi & getcap() const;
+      // Grow the Tree recursively
+      inline virtual void grow();
+      // Determine whether a given point is inside the bounding 
+      //
+      // Inputs:
+      //   p  query point 
+      // Returns true if the point p is inside this bounding volume
+      inline virtual bool inside(const Point & p) const;
+      // Compute the (partial) winding number of a given point p
+      // According to method
+      //  
+      // Inputs:
+      //   p  query point 
+      // Returns winding number 
+      inline double winding_number(const Point & p) const;
+      // Same as above, but always computes winding number using exact method
+      // (sum over every facet)
+      inline double winding_number_all(const Point & p) const;
+      // Same as above, but always computes using sum over tesslated boundary
+      inline double winding_number_boundary(const Point & p) const;
+      //// Same as winding_number above, but if max_simple_abs_winding_number is
+      //// less than some threshold min_max_w just return 0 (colloquially the "fast
+      //// multipole method)
+      ////
+      ////
+      //// Inputs:
+      ////   p  query point 
+      ////   min_max_w  minimum max simple w to be processed
+      //// Returns approximate winding number
+      //double winding_number_approx_simple(
+      //  const Point & p, 
+      //  const double min_max_w);
+      // Print contents of Tree
+      //
+      // Optional input:
+      //   tab  tab to show depth
+      inline void print(const char * tab="");
+      // Determine max absolute winding number
+      //
+      // Inputs:
+      //   p  query point 
+      // Returns max winding number of 
+      inline virtual double max_abs_winding_number(const Point & p) const; 
+      // Same as above, but stronger assumptions on (V,F). Assumes (V,F) is a
+      // simple polyhedron
+      inline virtual double max_simple_abs_winding_number(const Point & p) const;
+      // Compute or read cached winding number for point p with respect to mesh
+      // in bounding box, recursing according to approximation criteria
+      //
+      // Inputs:
+      //   p  query point 
+      //   that  WindingNumberTree containing mesh w.r.t. which we're computing w.n.
+      // Returns cached winding number
+      inline virtual double cached_winding_number(const WindingNumberTree & that, const Point & p) const;
+  };
+}
+
+// Implementation
+
+#include "WindingNumberTree.h"
+#include "winding_number.h"
+#include "triangle_fan.h"
+#include "exterior_edges.h"
+
+#include <igl/PI.h>
+#include <igl/remove_duplicate_vertices.h>
+
+#include <iostream>
+#include <limits>
+
+//template <typename Point>
+//WindingNumberMethod WindingNumberTree<Point>::method = EXACT_WINDING_NUMBER_METHOD;
+//template <typename Point>
+//double WindingNumberTree<Point>::min_max_w = 0;
+template <typename Point>
+std::map< std::pair<const igl::WindingNumberTree<Point>*,const igl::WindingNumberTree<Point>*>, double>
+  igl::WindingNumberTree<Point>::cached;
+
+static Eigen::MatrixXd dummyV;
+template <typename Point>
+inline igl::WindingNumberTree<Point>::WindingNumberTree(
+  const Eigen::MatrixXd & _V,
+  const Eigen::MatrixXi & _F):
+  method(EXACT_WINDING_NUMBER_METHOD),
+  parent(NULL),
+  V(dummyV),
+  SV(),
+  F(),
+  //boundary(igl::boundary_faces<Eigen::MatrixXi,Eigen::MatrixXi>(F))
+  cap(),
+  radius(std::numeric_limits<double>::infinity()),
+  center(0,0,0)
+{
+  // Remove any exactly duplicate vertices
+  // Q: Can this ever increase the complexity of the boundary?
+  // Q: Would we gain even more by remove almost exactly duplicate vertices?
+  Eigen::MatrixXi SF,SVI,SVJ;
+  igl::remove_duplicate_vertices(_V,_F,0.0,SV,SVI,SVJ,F);
+  triangle_fan(exterior_edges(F),cap);
+  V = SV;
+}
+
+template <typename Point>
+inline igl::WindingNumberTree<Point>::WindingNumberTree(
+  const igl::WindingNumberTree<Point> & parent,
+  const Eigen::MatrixXi & _F):
+  method(parent.method),
+  parent(&parent),
+  V(parent.V),
+  SV(),
+  F(_F),
+  cap(triangle_fan(exterior_edges(_F)))
+{
+}
+
+template <typename Point>
+inline igl::WindingNumberTree<Point>::~WindingNumberTree()
+{
+  using namespace std;
+  // Delete children
+  typename list<WindingNumberTree<Point>* >::iterator cit = children.begin();
+  while(cit != children.end())
+  {
+    // clear the memory of this item
+    delete (* cit);
+    // erase from list, returns next element in iterator
+    cit = children.erase(cit);
+  }
+}
+      
+template <typename Point>
+inline void igl::WindingNumberTree<Point>::set_method(const WindingNumberMethod & m)
+{
+  this->method = m;
+  for(auto child : children)
+  {
+    child->set_method(m);
+  }
+}
+
+template <typename Point>
+inline const Eigen::MatrixXd & igl::WindingNumberTree<Point>::getV() const
+{
+  return V;
+}
+
+template <typename Point>
+inline const Eigen::MatrixXi & igl::WindingNumberTree<Point>::getF() const
+{
+  return F;
+}
+
+template <typename Point>
+inline const Eigen::MatrixXi & igl::WindingNumberTree<Point>::getcap() const
+{
+  return cap;
+}
+
+template <typename Point>
+inline void igl::WindingNumberTree<Point>::grow()
+{
+  // Don't grow
+  return;
+}
+
+template <typename Point>
+inline bool igl::WindingNumberTree<Point>::inside(const Point & /*p*/) const
+{
+  return true;
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::winding_number(const Point & p) const
+{
+  using namespace std;
+  //cout<<"+"<<boundary.rows();
+  // If inside then we need to be careful
+  if(inside(p))
+  {
+    // If not a leaf then recurse
+    if(children.size()>0)
+    {
+      // Recurse on each child and accumulate
+      double sum = 0;
+      for(
+        typename list<WindingNumberTree<Point>* >::const_iterator cit = children.begin();
+        cit != children.end();
+        cit++)
+      {
+        switch(method)
+        {
+          case EXACT_WINDING_NUMBER_METHOD:
+            sum += (*cit)->winding_number(p);
+            break;
+          case APPROX_SIMPLE_WINDING_NUMBER_METHOD:
+          case APPROX_CACHE_WINDING_NUMBER_METHOD:
+            //if((*cit)->max_simple_abs_winding_number(p) > min_max_w)
+            //{
+              sum += (*cit)->winding_number(p);
+            //}
+            break;
+          default:
+            assert(false);
+            break;
+        }
+      }
+      return sum;
+    }else
+    {
+      return winding_number_all(p);
+    }
+  }else{
+    // Otherwise we can just consider boundary
+    // Q: If we using the "multipole" method should we also subdivide the
+    // boundary case?
+    if((cap.rows() - 2) < F.rows())
+    {
+      switch(method)
+      {
+        case EXACT_WINDING_NUMBER_METHOD:
+          return winding_number_boundary(p);
+        case APPROX_SIMPLE_WINDING_NUMBER_METHOD:
+        {
+          double dist = (p-center).norm();
+          // Radius is already an overestimate of inside
+          if(dist>1.0*radius)
+          {
+            return 0;
+          }else
+          {
+            return winding_number_boundary(p);
+          }
+        }
+        case APPROX_CACHE_WINDING_NUMBER_METHOD:
+        {
+          return parent->cached_winding_number(*this,p);
+        }
+        default: assert(false);break;
+      }
+    }else
+    {
+      // doesn't pay off to use boundary
+      return winding_number_all(p);
+    }
+  }
+  return 0;
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::winding_number_all(const Point & p) const
+{
+  double w = 0;
+  igl::winding_number_3(
+    V.data(),
+    V.rows(),
+    F.data(),
+    F.rows(),
+    p.data(),
+    1,
+    &w);
+  return w;
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::winding_number_boundary(const Point & p) const
+{
+  using namespace Eigen;
+  using namespace std;
+
+  double w = 0;
+  igl::winding_number_3(
+    V.data(),
+    V.rows(),
+    cap.data(),
+    cap.rows(),
+    &p[0],
+    1,
+    &w);
+  return w;
+}
+
+//template <typename Point>
+//inline double igl::WindingNumberTree<Point>::winding_number_approx_simple(
+//  const Point & p, 
+//  const double min_max_w)
+//{
+//  using namespace std;
+//  if(max_simple_abs_winding_number(p) > min_max_w)
+//  {
+//    return winding_number(p);
+//  }else
+//  {
+//    cout<<"Skipped! "<<max_simple_abs_winding_number(p)<<"<"<<min_max_w<<endl;
+//    return 0;
+//  }
+//}
+
+template <typename Point>
+inline void igl::WindingNumberTree<Point>::print(const char * tab)
+{
+  using namespace std;
+  // Print all facets
+  cout<<tab<<"["<<endl<<F<<endl<<"]";
+  // Print children
+  for(
+      typename list<WindingNumberTree<Point>* >::iterator cit = children.begin();
+      cit != children.end();
+      cit++)
+  {
+    cout<<","<<endl;
+    (*cit)->print((string(tab)+"").c_str());
+  }
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::max_abs_winding_number(const Point & p) const
+{
+  return std::numeric_limits<double>::infinity();
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::max_simple_abs_winding_number(const Point & p) const
+{
+  using namespace std;
+  return numeric_limits<double>::infinity();
+}
+
+template <typename Point>
+inline double igl::WindingNumberTree<Point>::cached_winding_number(
+  const igl::WindingNumberTree<Point> & that,
+  const Point & p) const
+{
+  using namespace std;
+  using namespace igl;
+  // Simple metric for "far".
+  //   this             that
+  //                   --------
+  //   -----          /   |    \ .
+  //  /  r  \        /    R     \ .
+  // | p !   |      |     !      |
+  //  \_____/        \          /
+  //                  \________/
+  //
+  // 
+  // a = angle formed by trapazoid formed by raising sides with lengths r and R
+  // at respective centers.
+  //
+  // a = atan2(R-r,d), where d is the distance between centers
+
+  // That should be bigger (what about parent? what about sister?)
+  bool far = this->radius<that.radius;
+  if(far)
+  {
+    double a = atan2(
+      that.radius - this->radius,
+      (that.center - this->center).norm());
+    assert(a>0);
+    far = (a<PI/8.0);
+  }
+
+  if(far)
+  {
+    // Not implemented yet
+    pair<const WindingNumberTree*,const WindingNumberTree*> this_that(this,&that);
+    // Need to compute it for first time?
+    if(cached.count(this_that)==0)
+    {
+      cached[this_that] = 
+        that.winding_number_boundary(this->center);
+    }
+    return cached[this_that];
+  }else if(children.size() == 0)
+  {
+    // not far and hierarchy ended too soon: can't use cache
+    return that.winding_number_boundary(p);
+  }else
+  {
+    for(
+      typename list<WindingNumberTree<Point>* >::const_iterator cit = children.begin();
+      cit != children.end();
+      cit++)
+    {
+      if((*cit)->inside(p))
+      {
+        return (*cit)->cached_winding_number(that,p);
+      }
+    }
+    // Not inside any children? This can totally happen because bounding boxes
+    // are set to bound contained facets. So sibilings may overlap and their
+    // union may not contain their parent (though, their union is certainly a
+    // subset of their parent).
+    assert(false);
+  }
+  return 0;
+}
+
+// Explicit instanciation
+//template class igl::WindingNumberTree<Eigen::Vector3d >;
+
+#endif

include/igl/exterior_edges.cpp

@@ -0,0 +1,104 @@
+#include "exterior_edges.h"
+#include "all_edges.h"
+
+#include <cassert>
+#include <map>
+#include <utility>
+
+//template <typename T> inline int sgn(T val) {
+//      return (T(0) < val) - (val < T(0));
+//}
+
+//static void mod2(std::pair<const std::pair<const int, const int>, int>& p)
+//{
+//  using namespace std;
+//  // Be sure that sign of mod matches sign of argument
+//  p.second = p.second%2 ? sgn(p.second) : 0;
+//}
+
+//// http://stackoverflow.com/a/5517869/148668
+//struct Compare
+//{
+//   int i;
+//   Compare(const int& i) : i(i) {}
+//};
+//bool operator==(const std::pair<std::pair<const int, const int>,int>&p, const Compare& c)
+//{
+//  return c.i == p.second;
+//}
+//bool operator==(const Compare& c, const std::pair<std::pair<const int, const int>, int> &p)
+//{
+//  return c.i == p.second;
+//}
+
+IGL_INLINE void igl::exterior_edges(
+  const Eigen::MatrixXi & F,
+  Eigen::MatrixXi & E)
+{
+  using namespace Eigen;
+  using namespace std;
+  assert(F.cols() == 3);
+  MatrixXi all_E;
+  all_edges(F,all_E);
+  // Count occurances looking only at pairs (i,j) where i<j, so we count and
+  // edge i-->j as +1 and j-->i as -1
+  map<pair<const int,const int>,int> C;
+  // Loop over edges
+  for(int e = 0;e<all_E.rows();e++)
+  {
+    int i = all_E(e,0);
+    int j = all_E(e,1);
+    if(i<j)
+    {
+      // Forward direction --> +1
+      C[pair<const int,const int>(i,j)]++;
+    }else
+    {
+      // Backward direction --> -1
+      C[pair<const int,const int>(j,i)]--;
+    }
+  }
+  // Q: Why mod off factors of 2? +1/-1 already takes care of interior edges?
+  //// Mod off any factors of 2
+  //for_each(C.begin(),C.end(),mod2);
+  //int zeros = (int) count(C.begin(),C.end(),Compare(0));
+  //E.resize(C.size() - zeros,2);
+  E.resize(all_E.rows(),all_E.cols());
+  int e = 0;
+  // Find all edges with -1 or 1 occurances, flipping direction for -1
+  for(
+    map<pair<const int, const int>,int>::const_iterator cit=C.begin();
+    cit!=C.end();
+    cit++)
+  {
+    int i,j;
+    if(cit->second > 0)
+    {
+      i = cit->first.first;
+      j = cit->first.second;
+    } else if(cit->second < 0)
+    {
+      i = cit->first.second;
+      j = cit->first.first;
+    } else if(cit->second == 0)
+    {
+      continue;
+    }
+    for(int k = 0;k<abs(cit->second);k++)
+    {
+      E(e,0) = i;
+      E(e,1) = j;
+      e++;
+    }
+  }
+  E.conservativeResize(e,2);
+  assert(e == E.rows());
+}
+
+IGL_INLINE Eigen::MatrixXi igl::exterior_edges( const Eigen::MatrixXi & F)
+{
+  using namespace Eigen;
+  MatrixXi E;
+  exterior_edges(F,E);
+  return E;
+}

include/igl/exterior_edges.h

@@ -0,0 +1,26 @@
+#ifndef IGL_EXTERIOR_EDGES_H
+#define IGL_EXTERIOR_EDGES_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+namespace igl
+{
+  // EXTERIOR_EDGES Determines boundary "edges" and also edges with an
+  // odd number of occurances where seeing edge (i,j) counts as +1 and seeing
+  // the opposite edge (j,i) counts as -1
+  //
+  // Inputs:
+  //   F  #F by simplex_size list of "faces"
+  // Outputs:
+  //   E  #E by simplex_size-1  list of exterior edges
+  //
+  IGL_INLINE void exterior_edges(
+    const Eigen::MatrixXi & F,
+    Eigen::MatrixXi & E);
+  // Inline version
+  Eigen::MatrixXi exterior_edges( const Eigen::MatrixXi & F);
+}
+#ifdef IGL_HEADER_ONLY
+#  include "exterior_edges.h"
+#endif
+
+#endif

include/igl/triangle_fan.cpp

@@ -0,0 +1,46 @@
+#include "triangle_fan.h"
+#include "exterior_edges.h"
+#include "list_to_matrix.h"
+
+IGL_INLINE void igl::triangle_fan(
+  const Eigen::MatrixXi & E,
+  Eigen::MatrixXi & cap)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  // Handle lame base case
+  if(E.size() == 0)
+  {
+    cap.resize(0,E.cols()+1);
+    return;
+  }
+  // "Triangulate" aka "close" the E trivially with facets
+  // Note: in 2D we need to know if E endpoints are incoming or
+  // outgoing (left or right). Thus this will not work.
+  assert(E.cols() == 2);
+  // Arbitrary starting vertex
+  int s = E(int(((double)rand() / RAND_MAX)*E.rows()),0);
+  vector<vector<int> >  lcap;
+  for(int i = 0;i<E.rows();i++)
+  {
+    // Skip edges incident on s (they would be zero-area)
+    if(E(i,0) == s || E(i,1) == s)
+    {
+      continue;
+    }
+    vector<int> e(3);
+    e[0] = s;
+    e[1] = E(i,0);
+    e[2] = E(i,1);
+    lcap.push_back(e);
+  }
+  list_to_matrix(lcap,cap);
+}
+
+IGL_INLINE Eigen::MatrixXi igl::triangle_fan( const Eigen::MatrixXi & E)
+{
+  Eigen::MatrixXi cap;
+  triangle_fan(E,cap);
+  return cap;
+}

include/igl/triangle_fan.h

@@ -0,0 +1,23 @@
+#ifndef IGL_TRIANGLE_FAN_H
+#define IGL_TRIANGLE_FAN_H
+#include "igl_inline.h"
+#include <Eigen/Dense>
+namespace igl
+{
+  // Given a list of faces tesselate all of the "exterior" edges forming another
+  // list of 
+  //
+  // Inputs:
+  //   E  #E by simplex_size-1  list of exterior edges (see exterior_edges.h)
+  // Outputs:
+  //   cap  #cap by simplex_size  list of "faces" tessleating the boundary edges
+  IGL_INLINE void triangle_fan(
+    const Eigen::MatrixXi & E,
+    Eigen::MatrixXi & cap);
+  // In-line version
+  IGL_INLINE Eigen::MatrixXi triangle_fan( const Eigen::MatrixXi & E);
+}
+#ifdef IGL_HEADER_ONLY
+#  include "triangle_fan.h"
+#endif
+#endif

include/igl/winding_number.cpp

@@ -0,0 +1,215 @@
+#include "winding_number.h"
+#include "WindingNumberAABB.h"
+
+#include <igl/PI.h>
+#include <cmath>
+
+IGL_INLINE void igl::winding_number(
+  const Eigen::MatrixXd & V,
+  const Eigen::MatrixXi & F,
+  const Eigen::MatrixXd & O,
+  Eigen::VectorXd & W)
+{
+  using namespace Eigen;
+  // make room for output
+  W.resize(O.rows(),1);
+  switch(F.cols())
+  {
+    case 2:
+      return winding_number_2(
+        V.data(),
+        V.rows(),
+        F.data(),
+        F.rows(),
+        O.data(),
+        O.rows(),
+        W.data());
+    case 3:
+    {
+      WindingNumberAABB<Vector3d> hier(V,F);
+      hier.grow();
+      // loop over origins
+      const int no = O.rows();
+#   pragma omp parallel for if (no>IGL_WINDING_NUMBER_OMP_MIN_VALUE)
+      for(int o = 0;o<no;o++)
+      {
+        Vector3d p = O.row(o);
+        W(o) = hier.winding_number(p);
+      }
+    }
+    default: assert(false && "Bad simplex size"); break;
+  }
+}
+
+template <typename DerivedF>
+IGL_INLINE void igl::winding_number_3(
+  const double * V,
+  const int n,
+  const DerivedF * F,
+  const int m,
+  const double * O,
+  const int no,
+  double * S)
+{
+  // Initialize output
+  // loop over origins
+#pragma omp parallel for if (no>IGL_WINDING_NUMBER_OMP_MIN_VALUE)
+  for(int o = 0;o<no;o++)
+  {
+    S[o] = 0;
+  }
+  // Only use parallel for if there are many facets and more than one origin.
+  // Assumes that if there is exactly one origin then this is being called
+  // within an outer for loop which may be parallel
+#pragma omp parallel for if (m>IGL_WINDING_NUMBER_OMP_MIN_VALUE && no>1)
+  // loop over faces
+  for(int f = 0;f<m;f++)
+  {
+    // Gather corners 
+    double C[3][3];
+    // loop around triangle
+    for(int t=0;t<3;t++)
+    {
+      // loop over dimensions
+      for(int d = 0;d<3;d++)
+      {
+        // Indices are offset by 1
+        int Ff = F[m*t + f];
+        C[t][d] = V[d*n + Ff];
+      }
+    }
+    // loop over origins
+    for(int o = 0;o<no;o++)
+    {
+      // Gather vectors to corners
+      double v[3][3];
+      double vl[3];
+      // loop around triangle
+      for(int t=0;t<3;t++)
+      {
+        vl[t] = 0;
+        // loop over dimensions
+        for(int d = 0;d<3;d++)
+        {
+          v[t][d] = C[t][d] - O[d*no + o];
+          // compute edge length contribution
+          vl[t] += v[t][d]*v[t][d];
+        }
+        // finish edge length computation
+        // Matlab crashes on NaN
+        if(vl[t]!=0)
+        {
+          vl[t] /= sqrt(vl[t]);
+        }
+      }
+      //printf("\n");
+      // Compute determinant
+      double detf = 
+        v[0][0]*v[1][1]*v[2][2]+
+        v[1][0]*v[2][1]*v[0][2]+
+        v[2][0]*v[0][1]*v[1][2]-
+        v[2][0]*v[1][1]*v[0][2]-
+        v[1][0]*v[0][1]*v[2][2]-
+        v[0][0]*v[2][1]*v[1][2];
+      // Compute pairwise dotproducts
+      double dp[3];
+      dp[0] = v[1][0]*v[2][0];
+      dp[0] += v[1][1]*v[2][1];
+      dp[0] += v[1][2]*v[2][2];
+      dp[1] = v[2][0]*v[0][0];
+      dp[1] += v[2][1]*v[0][1];
+      dp[1] += v[2][2]*v[0][2];
+      dp[2] = v[0][0]*v[1][0];
+      dp[2] += v[0][1]*v[1][1];
+      dp[2] += v[0][2]*v[1][2];
+      // Compute winding number
+      // Only divide by TWO_PI instead of 4*pi because there was a 2 out front
+      double val = atan2(detf,
+        vl[0]*vl[1]*vl[2] + 
+        dp[0]*vl[0] +
+        dp[1]*vl[1] +
+        dp[2]*vl[2]) / (2.*igl::PI);
+#pragma omp atomic
+      S[o] += val;
+    }
+  }
+}
+
+template <typename DerivedF>
+IGL_INLINE void igl::winding_number_2(
+  const double * V,
+  const int n,
+  const DerivedF * F,
+  const int m,
+  const double * O,
+  const int no,
+  double * S)
+{
+  // Initialize output
+  // loop over origins
+#pragma omp parallel for if (no>IGL_WINDING_NUMBER_OMP_MIN_VALUE)
+  for(int o = 0;o<no;o++)
+  {
+    S[o] = 0;
+  }
+#pragma omp parallel for if (m>IGL_WINDING_NUMBER_OMP_MIN_VALUE && no>1)
+  // loop over faces
+  for(int f = 0;f<m;f++)
+  {
+    // Index of source and destination
+    int s = F[m*0 + f];
+    int d = F[m*1 + f];
+    // Positions of source and destination
+    double vs[2];
+    double vd[2];
+    vs[0] = V[0*n + s];
+    vs[1] = V[1*n + s];
+    vd[0] = V[0*n + d];
+    vd[1] = V[1*n + d];
+    
+    // loop over origins
+    for(int o = 0;o<no;o++)
+    {
+      // Gather vectors to source and destination
+      double o2vs[2];
+      double o2vd[2];
+      // and lengths
+      double o2vsl = 0;
+      double o2vdl = 0;
+      for(int i = 0;i<2;i++)
+      {
+        o2vs[i] = O[i*no + o] - vs[i];
+        o2vd[i] = O[i*no + o] - vd[i];
+        o2vsl += o2vs[i]*o2vs[i];
+        o2vdl += o2vd[i]*o2vd[i];
+      }
+      o2vsl = sqrt(o2vsl);
+      o2vdl = sqrt(o2vdl);
+      // Normalize
+      for(int i = 0;i<2;i++)
+      {
+        // Matlab crashes on NaN
+        if(o2vsl!=0)
+        {
+          o2vs[i] /= o2vsl;
+        }
+        if(o2vdl!=0)
+        {
+          o2vd[i] /= o2vdl;
+        }
+      }
+      double val =
+        atan2(o2vd[0]*o2vs[1]-o2vd[1]*o2vs[0],o2vd[0]*o2vs[0]+o2vd[1]*o2vs[1])/
+        (2.*igl::PI);
+#pragma omp atomic
+      S[o] += val;
+    }
+  }
+}
+
+#ifndef IGL_HEADER_ONLY
+// Explicit template specialization
+template void igl::winding_number_3<int>(double const*, int, int const*, int, double const*, int, double*);
+template void igl::winding_number_2<double>(double const*, int, double const*, int, double const*, int, double*);
+template void igl::winding_number_3<double>(double const*, int, double const*, int, double const*, int, double*);
+#endif

include/igl/winding_number.h

@@ -0,0 +1,67 @@
+#ifndef IGL_WINDING_NUMBER_H
+#define IGL_WINDING_NUMBER_H
+#include "igl_inline.h"
+#include <Eigen/Core>
+
+// Minimum number of iterms per openmp thread
+#ifndef IGL_WINDING_NUMBER_OMP_MIN_VALUE
+#  define IGL_WINDING_NUMBER_OMP_MIN_VALUE 1000
+#endif
+namespace igl
+{
+  // WINDING_NUMBER Compute the sum of solid angles of a triangle/tetrahedron
+  // described by points (vectors) V
+  //
+  // Templates:
+  //   dim  dimension of input
+  // Inputs:
+  //  V  n by 3 list of vertex positions
+  //  n  number of mesh vertices
+  //  F  #F by 3 list of triangle indices, minimum index is 0
+  //  m  number of faces
+  //  O  no by 3 list of origin positions
+  //  no  number of origins
+  // Outputs:
+  //  S  no by 1 list of winding numbers
+  //
+  // 3d
+  IGL_INLINE void winding_number(
+    const Eigen::MatrixXd & V,
+    const Eigen::MatrixXi & F,
+    const Eigen::MatrixXd & O,
+    Eigen::VectorXd & W);
+  template <typename DerivedF>
+  IGL_INLINE void winding_number_3(
+    const double * V,
+    const int n,
+    const DerivedF * F,
+    const int m,
+    const double * O,
+    const int no,
+    double * S);
+  // Only one evaluation origin
+  template <typename DerivedF>
+  IGL_INLINE void winding_number_3(
+    const double * V,
+    const int n,
+    const DerivedF * F,
+    const int m,
+    const double * O,
+    double * S);
+  // 2d
+  template <typename DerivedF>
+  IGL_INLINE void winding_number_2(
+    const double * V,
+    const int n,
+    const DerivedF * F,
+    const int m,
+    const double * O,
+    const int no,
+    double * S);
+}
+
+#ifdef IGL_HEADER_ONLY
+#  include "winding_number.h"
+#endif
+
+#endif

todos.txt

@@ -1,3 +1,5 @@
+- extras should be in namespaces, i.e. igl::tetgen::tetrahedralize
+- extra header blocks should be verbose, i.e. IGL_TETGEN_TETRAHEDRALIZE_H
 - reorient_facets_raycast depends on boost and embree
 - sort out `grad.*` vs `gradMat.*`
 - clean up externals