first draft ready for libigl

Former-commit-id: 9056ae18d515054227208686cfc4cd6ea41e13d7

include/igl/build_octree.cpp

@@ -2,135 +2,141 @@
 #include <vector>
 #include <queue>
 
-const int MAX_DEPTH = 1000;
 
-int get_octant(const Eigen::RowVector3d & location, const Eigen::RowVector3d & center){
-    //Binary numbering of children:
-    //if we treat location as the origin,
-    //first bit is 1 if x is positive, 0 if not
-    //second bit is 1 if y is positive, 0 if not
-    //third bit is 1 if z is positive, 0 if not
-    //for example, the quadrant with negative x, positive y, positive z is:
-    //110 binary = 6 decimal
-    int index = 0;
-    if( location(0) >= center(0)){
+namespace igl {
+  template <typename DerivedP, typename IndexType, typename CentersType,
+  	typename WidthsType>
+  IGL_INLINE void build_octree(const Eigen::MatrixBase<DerivedP>& P,
+    std::vector<std::vector<IndexType> > & point_indices,
+    std::vector<Eigen::Matrix<IndexType,8,1>,
+    	Eigen::aligned_allocator<Eigen::Matrix<IndexType,8,1> > > & children,
+    std::vector<Eigen::Matrix<CentersType,1,3>,
+    	Eigen::aligned_allocator<Eigen::Matrix<CentersType,1,3> > > & centers,
+    std::vector<WidthsType> & widths)
+  {
+    const int MAX_DEPTH = 30000;
+
+    typedef Eigen::Matrix<IndexType,8,1> Vector8i;
+    typedef Eigen::Matrix<typename DerivedP::Scalar, 1, 3> RowVector3PType;
+    typedef Eigen::Matrix<CentersType, 1, 3> RowVector3CentersType;
+    
+    auto get_octant = [](RowVector3PType location,
+                         RowVector3CentersType center){
+      // We use a binary numbering of children. Treating the parent cell's
+      // center as the origin, we number the octants in the following manner:
+      // The first bit is 1 iff the octant's x coordinate is positive
+      // The second bit is 1 iff the octant's y coordinate is positive
+      // The third bit is 1 iff the octant's z coordinate is positive
+      //
+      // For example, the octant with negative x, positive y, positive z is:
+      // 110 binary = 6 decimal
+      IndexType index = 0;
+      if( location(0) >= center(0)){
         index = index + 1;
-    }
-    if( location(1) >= center(1)){
+      }
+      if( location(1) >= center(1)){
         index = index + 2;
-    }
-    if( location(2) >= center(2)){
+      }
+      if( location(2) >= center(2)){
         index = index + 4;
-    }
-    return index;
-}
-
-
-Eigen::RowVector3d translate_center(const Eigen::RowVector3d & parent_center, double h, int child_index){
-    Eigen::RowVector3d change_vector;
-    change_vector << h,h,h;
-    if(child_index % 2){ //positive x chilren are 1,3,4,7
-        change_vector(0) = -h;
-    }
-    if(child_index == 2 || child_index == 3 ||
-       child_index == 6 || child_index == 7){ //positive y children are 2,3,6,7
-        change_vector(1) = -h;
-    }
-    if(child_index > 3){ //positive z children are 4,5,6,7
-        change_vector(2) = -h;
-    }
-    return parent_center - change_vector;
-}
-
-void build_octree(const Eigen::MatrixXd & P,
-                  std::vector<std::vector<int> > & point_indices,
-                  std::vector<Eigen::Matrix<int,8,1>, Eigen::aligned_allocator<Eigen::Matrix<int,8,1>>> & children,
-                  std::vector<Eigen::RowVector3d, Eigen::aligned_allocator<Eigen::RowVector3d>> & centers,
-                  std::vector<double> & widths
-                  ){
-    typedef Eigen::Matrix<int,8,1> Vector8i;
-    typedef Eigen::Matrix<double,8,1> Vector8d;
+      }
+      return index;
+    };
 
-    // How many cells do we have so far?
-    int m = 0;
     
+    auto translate_center = [](const RowVector3CentersType & parent_center,
+                               const CentersType h,
+                               const IndexType child_index){
+      RowVector3CentersType change_vector;
+      change_vector << -h,-h,-h;
+      //positive x chilren are 1,3,4,7
+      if(child_index % 2){
+        change_vector(0) = h;
+      }
+      //positive y children are 2,3,6,7
+      if(child_index == 2 || child_index == 3 ||
+        child_index == 6 || child_index == 7){
+        change_vector(1) = h;
+      }
+      //positive z children are 4,5,6,7
+      if(child_index > 3){
+        change_vector(2) = h;
+      }
+      return parent_center + change_vector;
+    };
+  
+    // How many cells do we have so far?
+    IndexType m = 0;
+  
     // Useful list of number 0..7
     const Vector8i zero_to_seven = (Vector8i()<<0,1,2,3,4,5,6,7).finished();
     const Vector8i neg_ones = (Vector8i()<<-1,-1,-1,-1,-1,-1,-1,-1).finished();
-    
-    // This function variable needs to be declared before it is defined so that
-    // you can capture it and call it recursively. Annoyingly this means you need
-    // to fill in the function prototype here, too.
+  
     std::function< void(const int, const int) > helper;
-    helper = [
-              // These variables from the parent scope are "captured" you can read and
-              // write to them
-              &helper,
-              &zero_to_seven,&neg_ones,&P,&point_indices,&children,&centers,&widths,&m]
-    // The "-> bool" means that this function will return bool (I don't think
-    // you need this, but in case you do.)
-    (const int index, const int depth)-> void
+    helper = [&helper,&translate_center,&get_octant,&m,
+              &zero_to_seven,&neg_ones,&P,
+              &point_indices,&children,&centers,&widths]
+    (const IndexType index, const int depth)-> void
     {
-        if(point_indices.at(index).size() > 1 && depth < MAX_DEPTH){
-            //give the parent access to the children
-            children.at(index) = zero_to_seven.array() + m;
-            //make the children's data in our arrays
-            
-            //Add the children to the lists, as default children
-            double h = widths.at(index)/2;
-            Eigen::RowVector3d curr_center = centers.at(index);
-            for(int i = 0; i < 8; i++){
-                children.emplace_back(-1 * Eigen::MatrixXi::Ones(8,1));
-                point_indices.emplace_back(std::vector<int>());
-                centers.emplace_back(translate_center(curr_center,h,i));
-                widths.emplace_back(h);
-            }
-            
-            //Split up the points into the corresponding children
-            for(int j = 0; j < point_indices.at(index).size(); j++){
-                int curr_point_index = point_indices.at(index).at(j);
-                int cell_of_curr_point = get_octant(P.row(curr_point_index),curr_center)+m;
-                point_indices.at(cell_of_curr_point).emplace_back(curr_point_index);
-            }
-            
-            //Now increase m
-            m += 8;
-            
-            // Look ma, I'm calling myself.
-            for(int i = 0; i < 8; i++){
-                helper(children.at(index)(i),depth+1);
-            }
+      if(point_indices.at(index).size() > 1 && depth < MAX_DEPTH){
+        //give the parent access to the children
+        children.at(index) = zero_to_seven.array() + m;
+        //make the children's data in our arrays
+      
+        //Add the children to the lists, as default children
+        WidthsType h = widths.at(index)/2;
+        RowVector3CentersType curr_center = centers.at(index);
+        for(int i = 0; i < 8; i++){
+          children.emplace_back(neg_ones);
+          point_indices.emplace_back(std::vector<IndexType>());
+          centers.emplace_back(translate_center(curr_center,h,i));
+          widths.emplace_back(h);
         }
+      
+        //Split up the points into the corresponding children
+        for(int j = 0; j < point_indices.at(index).size(); j++){
+          IndexType curr_point_index = point_indices.at(index).at(j);
+          IndexType cell_of_curr_point =
+          	get_octant(P.row(curr_point_index),curr_center)+m;
+          point_indices.at(cell_of_curr_point).emplace_back(curr_point_index);
+        }
+      
+        //Now increase m
+        m += 8;
+      
+        // Look ma, I'm calling myself.
+        for(int i = 0; i < 8; i++){
+          helper(children.at(index)(i),depth+1);
+        }
+      }
     };
-    
+  
     {
-        std::vector<int> all(P.rows());
-        for(int i = 0;i<all.size();i++) all[i]=i;
-        point_indices.emplace_back(all);
+      std::vector<IndexType> all(P.rows());
+      for(IndexType i = 0;i<all.size();i++) all[i]=i;
+      point_indices.emplace_back(all);
     }
-    children.emplace_back(-1 * Eigen::MatrixXi::Ones(8,1));
-    
+    children.emplace_back(neg_ones);
+  
     //Get the minimum AABB for the points
-    Eigen::RowVector3d backleftbottom(P.col(0).minCoeff(),P.col(1).minCoeff(),P.col(2).minCoeff());
-    Eigen::RowVector3d frontrighttop(P.col(0).maxCoeff(),P.col(1).maxCoeff(),P.col(2).maxCoeff());
-    Eigen::RowVector3d aabb_center = (backleftbottom+frontrighttop)/2.0;
-    double aabb_width = std::max(std::max(frontrighttop(0) - backleftbottom(0),
+    RowVector3PType backleftbottom(P.col(0).minCoeff(),
+                                   P.col(1).minCoeff(),
+                                   P.col(2).minCoeff());
+    RowVector3PType frontrighttop(P.col(0).maxCoeff(),
+                                  P.col(1).maxCoeff(),
+                                  P.col(2).maxCoeff());
+    RowVector3CentersType aabb_center = (backleftbottom+frontrighttop)/2.0;
+    WidthsType aabb_width = std::max(std::max(
+                                          frontrighttop(0) - backleftbottom(0),
                                           frontrighttop(1) - backleftbottom(1)),
-                                 frontrighttop(2) - backleftbottom(2));
+                                 					frontrighttop(2) - backleftbottom(2));
     centers.emplace_back( aabb_center );
-    widths.emplace_back( aabb_width ); //Widths are the side length of the cube, (not half the side length)
+  
+    //Widths are the side length of the cube, (not half the side length):
+    widths.emplace_back( aabb_width );
     m++;
     // then you have to actually call the function
     helper(0,0);
-    
-    assert(point_indices.size() == m);
-    assert(children.size() == m);
-    assert(centers.size() == m);
-    assert(widths.size() == m);
+  }
 }
-//}
-//#ifndef IGL_STATIC_LIBRARY
-//#  include "point_areas_and_normals.cpp"
-//#endif
-//#endif
 

include/igl/build_octree.h

@@ -1,15 +1,57 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Gavin Barill <gavinpcb@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_BUILD_OCTREE
 #define IGL_BUILD_OCTREE
 #include <Eigen/Core>
-//namespace igl
-//{
-    void build_octree(const Eigen::MatrixXd & P,
-                    std::vector<std::vector<int> > & point_indices,
-                    std::vector<Eigen::Matrix<int,8,1>, Eigen::aligned_allocator<Eigen::Matrix<int,8,1>>> & children,
-                    std::vector<Eigen::RowVector3d, Eigen::aligned_allocator<Eigen::RowVector3d>> & centers,
-                    std::vector<double> & widths
-                    );
-//}
+#include "igl_inline.h"
+
+namespace igl
+{
+  // Given a set of 3D points P, generate data structures for a pointerless
+  // octree. Each cell stores its points, children, center location and width.
+  // Our octree is not dense. We use the following rule: if the current cell
+  // has any number of points, it will have all 8 children. A leaf cell will
+  // have -1's as its list of child indices.
+  //
+  // We use a binary numbering of children. Treating the parent cell's center
+  // as the origin, we number the octants in the following manner:
+  // The first bit is 1 iff the octant's x coordinate is positive
+  // The second bit is 1 iff the octant's y coordinate is positive
+  // The third bit is 1 iff the octant's z coordinate is positive
+  //
+  // For example, the octant with negative x, positive y, positive z is:
+  // 110 binary = 6 decimal
+  //
+  // Inputs:
+  //   P  #P by 3 list of point locations
+  //
+  // Outputs:
+  //   point_indices  a vector of vectors, where the ith entry is a vector of
+  //                  the indices into P that are the ith octree cell's points
+  //   children       a vector of vectors, where the ith entry is a vector of
+  //                  the ith octree cell's of octree children
+  //   centers        a vector where the ith entry is a 3d row vector
+  //                  representing the position of the ith cell's center
+  //   widths          a vector where the ith entry is the width of the ith
+  //                  octree cell
+  //
+  template <typename DerivedP, typename IndexType, typename CentersType,
+  	typename WidthsType>
+  IGL_INLINE void build_octree(const Eigen::MatrixBase<DerivedP>& P,
+    std::vector<std::vector<IndexType> > & point_indices,
+    std::vector<Eigen::Matrix<IndexType,8,1>,
+    	Eigen::aligned_allocator<Eigen::Matrix<IndexType,8,1> > > & children,
+    std::vector<Eigen::Matrix<CentersType,1,3>,
+      Eigen::aligned_allocator<Eigen::Matrix<CentersType,1,3> > > & centers,
+  	std::vector<WidthsType> & widths);
+}
+
 #ifndef IGL_STATIC_LIBRARY
 #  include "build_octree.cpp"
 #endif

include/igl/copyleft/cgal/point_areas_and_normals.cpp

@@ -0,0 +1,168 @@
+#include "point_areas_and_normals.h"
+#include "delaunay_triangulation.h"
+
+#include "../../colon.h"
+#include "../../slice.h"
+#include "../../slice_mask.h"
+#include "../../parallel_for.h"
+
+#include "CGAL/Exact_predicates_exact_constructions_kernel.h"
+#include "CGAL/Triangulation_vertex_base_with_info_2.h"
+#include "CGAL/Triangulation_data_structure_2.h"
+#include "CGAL/Delaunay_triangulation_2.h"
+
+
+
+typedef CGAL::Exact_predicates_inexact_constructions_kernel           Kernel;
+typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned int, Kernel> Vb;
+typedef CGAL::Triangulation_data_structure_2<Vb>                      Tds;
+typedef CGAL::Delaunay_triangulation_2<Kernel, Tds>                   Delaunay;
+typedef Kernel::Point_2                                               Point;
+
+namespace igl {
+	namespace copyleft{
+		namespace cgal{
+      template <typename DerivedP, typename DerivedI, typename DerivedO,
+      typename DerivedA, typename DerivedN>
+      IGL_INLINE void point_areas_and_normals(
+                                        const Eigen::MatrixBase<DerivedP>& P,
+                                        const Eigen::MatrixBase<DerivedI>& I,
+                                      	const Eigen::MatrixBase<DerivedO>& O,
+                                        Eigen::PlainObjectBase<DerivedA> & A,
+                                        Eigen::PlainObjectBase<DerivedN> & N)
+      {
+        typedef typename DerivedP::Scalar real;
+        typedef typename DerivedO::Scalar scalarO;
+        typedef typename DerivedA::Scalar scalarA;
+        typedef Eigen::Matrix<real,1,3> RowVec3;
+        typedef Eigen::Matrix<real,1,2> RowVec2;
+        
+        typedef Eigen::Matrix<real, Eigen::Dynamic, Eigen::Dynamic> MatrixP;
+        typedef Eigen::Matrix<scalarO, Eigen::Dynamic, Eigen::Dynamic> MatrixO;
+        typedef Eigen::Matrix<typename DerivedO::Scalar,
+        					Eigen::Dynamic, Eigen::Dynamic> VecotorO;
+        typedef Eigen::Matrix<typename DerivedI::Scalar,
+        					Eigen::Dynamic, Eigen::Dynamic> MatrixI;
+        
+        
+        
+        const int n = P.rows();
+        
+        assert(P.cols() == 3 && "P must have exactly 3 columns");
+        assert(P.rows() == O.rows()
+               && "P and O must have the same number of rows");
+        assert(P.rows() == I.rows()
+               && "P and I must have the same number of rows");
+        
+        A.setZero(n,1);
+        N.setZero(n,3);
+        igl::parallel_for(P.rows(),[&](int i)
+        {
+          MatrixI neighbor_index = I.row(i);
+          MatrixP neighbors;
+          igl::slice(P,neighbor_index,1,neighbors);
+          if(O.rows() && neighbors.rows() > 1){
+            MatrixO neighbor_normals;
+            igl::slice(O,neighbor_index,1,neighbor_normals);
+            Eigen::Matrix<scalarO,1,3> poi_normal = neighbor_normals.row(0);
+            Eigen::Matrix<scalarO,Eigen::Dynamic,1> dotprod =
+            								poi_normal(0)*neighbor_normals.col(0)
+            + poi_normal(1)*neighbor_normals.col(1)
+            + poi_normal(2)*neighbor_normals.col(2);
+            Eigen::Array<bool,Eigen::Dynamic,1> keep = dotprod.array() > 0;
+            igl::slice_mask(Eigen::MatrixXd(neighbors),keep,1,neighbors);
+          }
+          if(neighbors.rows() <= 2){
+            A(i) = 0;
+            N.row(i) = Eigen::RowVector3d::Zero();
+          } else {
+            //subtract the mean from neighbors, then take svd,
+            //the scores will be U*S. This is our pca plane fitting
+            RowVec3 mean = neighbors.colwise().mean();
+            MatrixP mean_centered = neighbors.rowwise() - mean;
+            Eigen::JacobiSVD<MatrixP> svd(mean_centered,
+                                    Eigen::ComputeThinU | Eigen::ComputeThinV);
+            MatrixP scores = svd.matrixU() * svd.singularValues().asDiagonal();
+            
+            N.row(i) = svd.matrixV().col(2).transpose();
+            if(N.row(i).dot(O.row(i)) < 0){
+              N.row(i) *= -1;
+            }
+            
+            MatrixP plane;
+            igl::slice(scores,igl::colon<int>(0,scores.rows()-1),
+                     igl::colon<int>(0,1),plane);
+            
+            std::vector< std::pair<Point,unsigned> > points;
+            //This is where we obtain a delaunay triangulation of the points
+            for(unsigned iter = 0; iter < plane.rows(); iter++){
+              points.push_back( std::make_pair(
+                      Point(plane(iter,0),plane(iter,1)), iter ) );
+            }
+            Delaunay triangulation;
+            triangulation.insert(points.begin(),points.end());
+            Eigen::MatrixXi F(triangulation.number_of_faces(),3);
+            int f_row = 0;
+            for(Delaunay::Finite_faces_iterator fit =
+                triangulation.finite_faces_begin();
+              	fit != triangulation.finite_faces_end(); ++fit) {
+              Delaunay::Face_handle face = fit;
+              F.row(f_row) = Eigen::RowVector3i((int)face->vertex(0)->info(),
+                                                (int)face->vertex(1)->info(),
+                                                (int)face->vertex(2)->info());
+              f_row++;
+            }
+            
+            //Here we calculate the voronoi area of the point
+            scalarA area_accumulator = 0;
+            for(int f = 0; f < F.rows(); f++){
+              int X = -1;
+              for(int face_iter = 0; face_iter < 3; face_iter++){
+                if(F(f,face_iter) == 0){
+                  X = face_iter;
+                }
+              }
+              if(X >= 0){
+              //Triangle XYZ with X being the point we want the area of
+                int Y = (X+1)%3;
+                int Z = (X+2)%3;
+                scalarA x = (plane.row(F(f,Y))-plane.row(F(f,Z))).norm();
+                scalarA y = (plane.row(F(f,X))-plane.row(F(f,Z))).norm();
+                scalarA z = (plane.row(F(f,Y))-plane.row(F(f,X))).norm();
+                scalarA cosX = (z*z + y*y - x*x)/(2*y*z);
+                scalarA cosY = (z*z + x*x - y*y)/(2*x*z);
+                scalarA cosZ = (x*x + y*y - z*z)/(2*y*x);
+                Eigen::Matrix<scalarA,1,3> barycentric;
+                barycentric << x*cosX, y*cosY, z*cosZ;
+                barycentric /= (barycentric(0)+barycentric(1)+barycentric(2));
+                
+                //TODO: to make numerically stable, reorder so that x≥y≥z:
+                scalarA full_area = 0.25*std::sqrt(
+                                    (x+(y+z))*(z-(x-y))*(z+(x-y))*(x+(y-z)));
+                Eigen::Matrix<scalarA,1,3> partial_area =
+                																		barycentric * full_area;
+                if(cosX < 0){
+                  area_accumulator += 0.5*full_area;
+                } else if (cosY < 0 || cosZ < 0){
+                  area_accumulator += 0.25*full_area;
+                } else {
+                  area_accumulator += (partial_area(1) + partial_area(2))/2;
+                }
+              }
+            }
+            
+            if(std::isfinite(area_accumulator)){
+              A(i) = area_accumulator;
+            } else {
+              A(i) = 0;
+            }
+          }
+        },1000);
+      }
+		}
+	}
+}
+
+
+
+

include/igl/copyleft/cgal/point_areas_and_normals.h

@@ -0,0 +1,59 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Gavin Barill <gavinpcb@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_POINT_AREAS_AND_NORMALS_H
+#define IGL_POINT_AREAS_AND_NORMALS_H
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+namespace igl
+{
+	namespace copyleft
+  {
+  	namespace cgal
+    {
+    // Given a 3D set of points P, each with a list of k-nearest-neighbours,
+    // estimate the geodesic voronoi area associated with each point.
+    //
+  	// The k nearest neighbours may be known from running igl::knn_octree on
+    // the output data from igl::build_octree. We reccomend using a k value
+    // between 15 and 20 inclusive for accurate area estimation.
+    //
+    // If given, O represents an input oritentaiton to each point. This could
+    // be an initial guess at the normals, a vector to a camera position (for
+    // scanned data), or ground truth normals if you're only interested in the
+    // area estimation. It's used to ensure area estimation only occurs using
+    // neighbors that are on the same side of the surface (ie for thin
+    // sheets), as well as to solve the orientation ambiguity of normal
+    // estimation.
+    //
+    // Inputs:
+    //   P  #P by 3 list of point locations
+    //   I  #P by k list of k-nearest-neighbor indices into P
+    //   O  #P by 3 list of point orientation vectors (optional)
+    // Outputs:
+    //   A  #P list of estimated areas
+    //   N  #P by 3 list of estimated normals
+		template <typename DerivedP, typename DerivedI, typename DerivedO,
+      typename DerivedA, typename DerivedN>
+    IGL_INLINE void point_areas_and_normals(
+                                        const Eigen::MatrixBase<DerivedP>& P,
+                                        const Eigen::MatrixBase<DerivedI>& I,
+                                        const Eigen::MatrixBase<DerivedO>& O,
+                                        Eigen::PlainObjectBase<DerivedA> & A,
+                                        Eigen::PlainObjectBase<DerivedN> & N);
+      
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_areas_and_normals.cpp"
+#endif
+
+#endif
+

include/igl/fast_winding_number.cpp

@@ -1,207 +1,348 @@
 #include "fast_winding_number.h"
+#include "build_octree.h"
+#include "knn_octree.h"
+#include "parallel_for.h"
+#include "copyleft/cgal/point_areas_and_normals.h"
 #include <vector>
-#include <iostream>
-#include <igl/parallel_for.h>
-
-void fast_winding_number_precompute(const Eigen::MatrixXd & P,
-                                        const Eigen::MatrixXd & N,
-                                        const Eigen::VectorXd & A,
-                                        const std::vector<std::vector<int> > & point_indices,
-                                        const std::vector<Eigen::Matrix<int,8,1>, Eigen::aligned_allocator<Eigen::Matrix<int,8,1>>> & children,
-                                        const int expansion_order,
-                                        Eigen::MatrixXd & CM,
-                                        Eigen::VectorXd & R,
-                                        Eigen::MatrixXd & EC
-                                    ){
-    int m = children.size();
-    int num_terms;
-    if(expansion_order == 0){
-        num_terms = 3;
-    } else if(expansion_order ==1){
-        num_terms = 3 + 9;
-    } else if(expansion_order == 2){
-        num_terms = 3 + 9 + 27;
-    } else {
-        assert(false);
-    }
+
+namespace igl {
+template <typename DerivedP, typename DerivedA, typename DerivedN,
+  	typename Index, typename DerivedCM, typename DerivedR, typename DerivedEC>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+              const Eigen::MatrixBase<DerivedN>& N,
+              const Eigen::MatrixBase<DerivedA>& A,
+              const std::vector<std::vector<Index> > & point_indices,
+              const std::vector<Eigen::Matrix<Index,8,1>,
+              	Eigen::aligned_allocator<Eigen::Matrix<Index,8,1>>> & children,
+              const int expansion_order,
+              Eigen::PlainObjectBase<DerivedCM>& CM,
+              Eigen::PlainObjectBase<DerivedR>& R,
+              Eigen::PlainObjectBase<DerivedEC>& EC
+              ){
+      typedef typename DerivedP::Scalar real_p;
+      typedef typename DerivedN::Scalar real_n;
+      typedef typename DerivedA::Scalar real_a;
+      typedef typename DerivedCM::Scalar real_cm;
+      typedef typename DerivedR::Scalar real_r;
+      typedef typename DerivedEC::Scalar real_ec;
     
-    R.resize(m);
-    CM.resize(m,3);
-    EC.resize(m,num_terms);
-    EC = Eigen::MatrixXd::Zero(m,num_terms);
-    std::function< void(const int) > helper;
-    helper = [&helper,
-              &P,&N,&A,&expansion_order,&point_indices,&children,&EC,&R,&CM]
-    (const int index)-> void
-    {
-        double sum_area = 0;
-        
-        Eigen::RowVector3d masscenter = Eigen::RowVector3d::Zero();
-        Eigen::RowVector3d zeroth_expansion = Eigen::RowVector3d::Zero();
-        double areatotal = 0.0;
-        for(int j = 0; j < point_indices.at(index).size(); j++){
-            int curr_point_index = point_indices.at(index).at(j);
+      typedef Eigen::Matrix<real_p,1,3> RowVec3p;
+    
+      int m = children.size();
+      int num_terms;
+      if(expansion_order == 0){
+          num_terms = 3;
+      } else if(expansion_order ==1){
+          num_terms = 3 + 9;
+      } else if(expansion_order == 2){
+          num_terms = 3 + 9 + 27;
+      } else {
+          assert(false);
+      }
+    
+      R.resize(m);
+      CM.resize(m,3);
+      EC.resize(m,num_terms);
+      EC.setZero(m,num_terms);
+      std::function< void(const int) > helper;
+      helper = [&helper,
+                &P,&N,&A,&expansion_order,&point_indices,&children,&EC,&R,&CM]
+      (const int index)-> void
+      {
+          Eigen::Matrix<real_cm,1,3> masscenter;
+          masscenter << 0,0,0;
+          Eigen::Matrix<real_ec,1,3> zeroth_expansion;
+          masscenter << 0,0,0;
+          real_p areatotal = 0.0;
+          for(int j = 0; j < point_indices.at(index).size(); j++){
+              int curr_point_index = point_indices.at(index).at(j);
             
-            areatotal += A(curr_point_index);
-            masscenter += A(curr_point_index)*P.row(curr_point_index);
-            zeroth_expansion += A(curr_point_index)*N.row(curr_point_index);
-        }
+              areatotal += A(curr_point_index);
+              masscenter += A(curr_point_index)*P.row(curr_point_index);
+              zeroth_expansion += A(curr_point_index)*N.row(curr_point_index);
+          }
         
-        masscenter = masscenter/areatotal;
-        CM.row(index) = masscenter;
-        EC.block<1,3>(index,0) = zeroth_expansion;
+          masscenter = masscenter/areatotal;
+          CM.row(index) = masscenter;
+          EC.block(index,0,1,3) = zeroth_expansion;
         
-        double max_norm = 0;
-        double curr_norm;
+          real_r max_norm = 0;
+          real_r curr_norm;
         
-        for(int i = 0; i < point_indices.at(index).size(); i++){
-            //Get max distance from center of mass:
-            int curr_point_index = point_indices.at(index).at(i);
-            Eigen::RowVector3d point = P.row(curr_point_index)-masscenter;
-            curr_norm = point.norm();
-            if(curr_norm > max_norm){
-                max_norm = curr_norm;
-            }
+          for(int i = 0; i < point_indices.at(index).size(); i++){
+              //Get max distance from center of mass:
+              int curr_point_index = point_indices.at(index).at(i);
+              Eigen::Matrix<real_r,1,3> point =
+            			P.row(curr_point_index)-masscenter;
+              curr_norm = point.norm();
+              if(curr_norm > max_norm){
+                  max_norm = curr_norm;
+              }
             
-            //Calculate higher order terms if necessary
-            Eigen::Matrix3d TempCoeffs;
-            if(EC.cols() >= (3+9)){
-                TempCoeffs = A(curr_point_index) * point.transpose() * N.row(curr_point_index);
-                EC.block<1,9>(index,3) += Eigen::Map<Eigen::RowVectorXd>(TempCoeffs.data(), TempCoeffs.size());
-            }
+              //Calculate higher order terms if necessary
+              Eigen::Matrix<real_ec,3,3> TempCoeffs;
+              if(EC.cols() >= (3+9)){
+                  TempCoeffs = A(curr_point_index)*point.transpose()*
+                									N.row(curr_point_index);
+                  EC.block(index,3,1,9) +=
+                  Eigen::Map<Eigen::Matrix<real_ec,1,9> >(TempCoeffs.data(),
+                                                          TempCoeffs.size());
+              }
             
-            if(EC.cols() == (3+9+27)){
-                for(int k = 0; k < 3; k++){
-                    TempCoeffs = 0.5 * point(k) * (A(curr_point_index) * point.transpose() * N.row(curr_point_index)) ;
-                    EC.block<1,9>(index,12+9*k) += Eigen::Map<Eigen::RowVectorXd>(TempCoeffs.data(), TempCoeffs.size());
-                }
-            }
-        }
+              if(EC.cols() == (3+9+27)){
+                  for(int k = 0; k < 3; k++){
+                      TempCoeffs = 0.5 * point(k) * (A(curr_point_index)*
+                                    point.transpose()*N.row(curr_point_index));
+                      EC.block(index,12+9*k,1,9) += Eigen::Map<
+                    		Eigen::Matrix<real_ec,1,9> >(TempCoeffs.data(),
+                                                     TempCoeffs.size());
+                  }
+              }
+          }
         
-        R(index) = max_norm;
-        if(children.at(index)(0) != -1)
-        {
-            for(int i = 0; i < 8; i++){
-                int child = children.at(index)(i);
-                helper(child);
-            }
+          R(index) = max_norm;
+          if(children.at(index)(0) != -1)
+          {
+              for(int i = 0; i < 8; i++){
+                  int child = children.at(index)(i);
+                  helper(child);
+              }
+          }
+      };
+      helper(0);
+  }
+  
+  template <typename DerivedP, typename DerivedA, typename DerivedN,
+  	typename Index, typename DerivedCM, typename DerivedR, typename DerivedEC,
+  	typename DerivedQ, typename BetaType, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+              const Eigen::MatrixBase<DerivedN>& N,
+              const Eigen::MatrixBase<DerivedA>& A,
+              const std::vector<std::vector<Index> > & point_indices,
+              const std::vector<Eigen::Matrix<Index,8,1>,
+                Eigen::aligned_allocator<Eigen::Matrix<Index,8,1>>> & children,
+              const Eigen::MatrixBase<DerivedCM>& CM,
+              const Eigen::MatrixBase<DerivedR>& R,
+              const Eigen::MatrixBase<DerivedEC>& EC,
+              const Eigen::MatrixBase<DerivedQ>& Q,
+              const BetaType beta,
+              Eigen::PlainObjectBase<DerivedWN>& WN
+              ){
+  
+    typedef typename DerivedP::Scalar real_p;
+    typedef typename DerivedN::Scalar real_n;
+    typedef typename DerivedA::Scalar real_a;
+    typedef typename DerivedCM::Scalar real_cm;
+    typedef typename DerivedR::Scalar real_r;
+    typedef typename DerivedEC::Scalar real_ec;
+    typedef typename DerivedQ::Scalar real_q;
+    typedef typename DerivedWN::Scalar real_wn;
+  
+    typedef Eigen::Matrix<real_q,1,3> RowVec;
+    typedef Eigen::Matrix<real_ec,3,3> EC_3by3;
+  
+    auto direct_eval = [](const RowVec & loc,
+                          const Eigen::Matrix<real_ec,1,3> & anorm){
+        real_wn wn = (loc(0)*anorm(0)+loc(1)*anorm(1)+loc(2)*anorm(2))
+      															/(4.0*M_PI*std::pow(loc.norm(),3));
+        if(std::isnan(wn)){
+            return 0.5;
+        }else{
+            return wn;
         }
     };
-    helper(0);
-}
-
-double direct_eval(const Eigen::RowVector3d & loc, const Eigen::RowVector3d & anorm){
-    double wn = (loc(0)*anorm(0) + loc(1)*anorm(1) + loc(2)*anorm(2))/(4.0*M_PI*std::pow(loc.norm(),3));
-    if(std::isnan(wn)){
-        return 0.5;
-    }else{
-        return wn;
-    }
-}
-
-double expansion_eval(const Eigen::RowVector3d & loc, const Eigen::RowVectorXd & EC){
-    double wn = direct_eval(loc,EC.head<3>());
-    double r = loc.norm();
-    if(EC.size()>3){
-        Eigen::Matrix3d SecondDerivative = Eigen::Matrix3d::Identity()/(4.0*M_PI*std::pow(r,3));
+  
+    auto expansion_eval = [&direct_eval](const RowVec & loc,
+                                         const Eigen::RowVectorXd & EC){
+      real_wn wn = direct_eval(loc,EC.head<3>());
+      double r = loc.norm();
+      if(EC.size()>3){
+        Eigen::Matrix<real_ec,3,3> SecondDerivative =
+          	Eigen::Matrix<real_ec,3,3>::Identity()/(4.0*M_PI*std::pow(r,3));
         SecondDerivative += -3.0*loc.transpose()*loc/(4.0*M_PI*std::pow(r,5));
-        Eigen::RowVectorXd derivative_vector = Eigen::Map<Eigen::RowVectorXd>(SecondDerivative.data(), SecondDerivative.size());
+        Eigen::Matrix<real_ec,1,9> derivative_vector =
+        	Eigen::Map<Eigen::Matrix<real_ec,1,9> >(SecondDerivative.data(),
+                                                  SecondDerivative.size());
         wn += derivative_vector.cwiseProduct(EC.segment<9>(3)).sum();
-    }
-    if(EC.size()>3+9){
-        Eigen::Matrix3d ThirdDerivative;
-        for(int i = 0; i < 3; i++){
-            ThirdDerivative = 15.0*loc(i)*loc.transpose()*loc/(4.0*M_PI*std::pow(r,7));
-            Eigen::Matrix3d Diagonal;
-            Diagonal << loc(i), 0, 0,
-                        0, loc(i), 0,
-                        0, 0, loc(i);
-            Eigen::Matrix3d RowCol = Eigen::Matrix3d::Zero();
-            RowCol.row(i) = loc;
-            RowCol = RowCol + RowCol.transpose();
-            ThirdDerivative += -3.0/(4.0*M_PI*std::pow(r,5)) * (RowCol + Diagonal);
-            Eigen::RowVectorXd derivative_vector = Eigen::Map<Eigen::RowVectorXd>(ThirdDerivative.data(), ThirdDerivative.size());
-            wn += derivative_vector.cwiseProduct(EC.segment<9>(12 + i*9)).sum();
-        }
-    }
-    return wn;
-}
-
-void fast_winding_number(const Eigen::MatrixXd & P,
-                         const Eigen::MatrixXd & N,
-                         const Eigen::VectorXd & A,
-                         const std::vector<std::vector<int> > & point_indices,
-                         const std::vector<Eigen::Matrix<int,8,1>,               Eigen::aligned_allocator<Eigen::Matrix<int,8,1>>> & children,
-                         const Eigen::MatrixXd & CM,
-                         const Eigen::VectorXd & R,
-                         const Eigen::MatrixXd & EC,
-                         const Eigen::MatrixXd & Q,
-                         const double & beta,
-                         Eigen::VectorXd & WN
-                         ){
+      }
+      if(EC.size()>3+9){
+          Eigen::Matrix<real_ec,3,3> ThirdDerivative;
+          for(int i = 0; i < 3; i++){
+              ThirdDerivative =
+            			15.0*loc(i)*loc.transpose()*loc/(4.0*M_PI*std::pow(r,7));
+              Eigen::Matrix<real_ec,3,3> Diagonal;
+              Diagonal << loc(i), 0, 0,
+              0, loc(i), 0,
+              0, 0, loc(i);
+              Eigen::Matrix<real_ec,3,3> RowCol;
+              RowCol.setZero(3,3);
+              RowCol.row(i) = loc;
+              RowCol = RowCol + RowCol.transpose();
+              ThirdDerivative +=
+            			-3.0/(4.0*M_PI*std::pow(r,5))*(RowCol+Diagonal);
+              Eigen::Matrix<real_ec,1,9> derivative_vector =
+            		Eigen::Map<Eigen::Matrix<real_ec,1,9> >(ThirdDerivative.data(),
+                                                        ThirdDerivative.size());
+              wn += derivative_vector.cwiseProduct(
+                                              EC.segment<9>(12 + i*9)).sum();
+          }
+      }
+      return wn;
+    };
+  
     int m = Q.rows();
-    WN.resize(m);
-    
-    std::function< double(const Eigen::RowVector3d, const std::vector<int>) > helper;
+    WN.resize(m,1);
+  
+    std::function< real_wn(const RowVec, const std::vector<int>) > helper;
     helper = [&helper,
               &P,&N,&A,
               &point_indices,&children,
-              &CM,&R,&EC,&beta]
-    (const Eigen::RowVector3d query, const std::vector<int> near_indices)-> double
+              &CM,&R,&EC,&beta,
+              &direct_eval,&expansion_eval]
+    (const RowVec query, const std::vector<int> near_indices)-> real_wn
     {
-        std::vector<int> new_near_indices;
-        double wn = 0;
-        for(int i = 0; i < near_indices.size(); i++){
-            int index = near_indices.at(i);
-            //Leaf Case, Brute force
-            if(children.at(index)(0) == -1){
-                for(int j = 0; j < point_indices.at(index).size(); j++){
-                    int curr_row = point_indices.at(index).at(j);
-                    wn += direct_eval(P.row(curr_row)-query,N.row(curr_row)*A(curr_row));
-                }
-            }
-            //Non-Leaf Case
-            else {
-                for(int child = 0; child < 8; child++){
-                    int child_index = children.at(index)(child);
-                    if(point_indices.at(child_index).size() > 0){
-                        if((CM.row(child_index)-query).norm() > beta*R(child_index)){
-                            if(children.at(child_index)(0) == -1){
-                                for(int j = 0; j < point_indices.at(child_index).size(); j++){
-                                    int curr_row = point_indices.at(child_index).at(j);
-                                    wn += direct_eval(P.row(curr_row)-query,N.row(curr_row)*A(curr_row));
-                                }
-                            }else{
-                                wn += expansion_eval(CM.row(child_index)-query,EC.row(child_index));
-                            }
-                        }else {
-                            new_near_indices.emplace_back(child_index);
-                        }
+      std::vector<int> new_near_indices;
+      real_wn wn = 0;
+      for(int i = 0; i < near_indices.size(); i++){
+        int index = near_indices.at(i);
+        //Leaf Case, Brute force
+        if(children.at(index)(0) == -1){
+          for(int j = 0; j < point_indices.at(index).size(); j++){
+            int curr_row = point_indices.at(index).at(j);
+            wn += direct_eval(P.row(curr_row)-query,
+                              N.row(curr_row)*A(curr_row));
+          }
+        }
+        //Non-Leaf Case
+        else {
+          for(int child = 0; child < 8; child++){
+              int child_index = children.at(index)(child);
+              if(point_indices.at(child_index).size() > 0){
+                if((CM.row(child_index)-query).norm() > beta*R(child_index)){
+                  if(children.at(child_index)(0) == -1){
+                    for(int j=0;j<point_indices.at(child_index).size();j++){
+                      int curr_row = point_indices.at(child_index).at(j);
+                      wn += direct_eval(P.row(curr_row)-query,
+                                        N.row(curr_row)*A(curr_row));
                     }
-                }
+                  }else{
+                    wn += expansion_eval(CM.row(child_index)-query,
+                                         EC.row(child_index));
+                  }
+                }else {
+                  new_near_indices.emplace_back(child_index);
+              }
             }
+          }
         }
-        if(new_near_indices.size() > 0){
-            wn += helper(query,new_near_indices);
-        }
-        return wn;
+      }
+      if(new_near_indices.size() > 0){
+          wn += helper(query,new_near_indices);
+      }
+      return wn;
     };
-    
-    
-    if(beta >= 0){
-        std::vector<int> near_indices_start = {0};
-        igl::parallel_for(m,[&](int iter){
-            WN(iter) = helper(Q.row(iter),near_indices_start);
-        },1000);
+  
+  
+    if(beta > 0){
+      std::vector<int> near_indices_start = {0};
+      igl::parallel_for(m,[&](int iter){
+        WN(iter) = helper(Q.row(iter),near_indices_start);
+      },1000);
     } else {
-        igl::parallel_for(m,[&](int iter){
-            double wn = 0;
-            for(int j = 0; j <P.rows(); j++){
-                wn += direct_eval(P.row(j)-Q.row(iter),N.row(j)*A(j));
-            }
-            WN(iter) = wn;
-        },1000);
+      igl::parallel_for(m,[&](int iter){
+        double wn = 0;
+        for(int j = 0; j <P.rows(); j++){
+          wn += direct_eval(P.row(j)-Q.row(iter),N.row(j)*A(j));
+        }
+        WN(iter) = wn;
+      },1000);
     }
+  }
+  
+  template <typename DerivedP, typename DerivedA, typename DerivedN,
+  	typename DerivedQ, typename BetaType, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedA>& A,
+                                      const Eigen::MatrixBase<DerivedQ>& Q,
+                                      const int expansion_order,
+                                      const BetaType beta,
+                                      Eigen::PlainObjectBase<DerivedWN>& WN
+                                      ){
+    typedef typename DerivedWN::Scalar real;
+    typedef typename Eigen::Matrix<real,1,3> RowVec;
+  
+    std::vector<std::vector<int>> point_indices;
+    std::vector<Eigen::Matrix<int,8,1>,
+      Eigen::aligned_allocator<Eigen::Matrix<int,8,1> > > children;
+    std::vector<RowVec, Eigen::aligned_allocator<RowVec> > centers;
+    std::vector<real> widths;
+    Eigen::MatrixXi I;
+    Eigen::MatrixXd NotUsed;
+  
+    build_octree(P,point_indices,children,centers,widths);
+  
+    Eigen::MatrixXd EC, CM;
+    Eigen::VectorXd R;
+  
+    fast_winding_number(P,N,A,point_indices,children,
+                          expansion_order,CM,R,EC);
+    fast_winding_number(P,N,A,point_indices,children,CM,R,EC,Q,beta,WN);
+  }
+  
+  template <typename DerivedP, typename DerivedA, typename DerivedN,
+  	typename DerivedQ, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedA>& A,
+                                      const Eigen::MatrixBase<DerivedQ>& Q,
+                                      Eigen::PlainObjectBase<DerivedWN>& WN
+                                      ){
+    fast_winding_number(P,N,A,Q,2,2.0,WN);
+  }
+  
+  template <typename DerivedP, typename DerivedN, typename DerivedQ,
+  	typename BetaType, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedQ>& Q,
+                                      const int expansion_order,
+                                      const BetaType beta,
+                                      Eigen::PlainObjectBase<DerivedWN>& WN
+                                      ){
+    typedef typename DerivedWN::Scalar real;
+    typedef typename Eigen::Matrix<real,1,3> RowVec;
+  
+    std::vector<std::vector<int>> point_indices;
+    std::vector<Eigen::Matrix<int,8,1>,
+    	Eigen::aligned_allocator<Eigen::Matrix<int,8,1> > > children;
+    std::vector<RowVec, Eigen::aligned_allocator<RowVec> > centers;
+    std::vector<real> widths;
+    Eigen::MatrixXi I;
+    Eigen::Matrix<real,Eigen::Dynamic,Eigen::Dynamic> NotUsed;
+    Eigen::Matrix<real,Eigen::Dynamic,1> A;
+  
+    build_octree(P,point_indices,children,centers,widths);
+    knn_octree(P,21,point_indices,children,centers,widths,I);
+    copyleft::cgal::point_areas_and_normals(P,I,N,A,NotUsed);
+  
+    Eigen::MatrixXd EC, CM;
+    Eigen::VectorXd R;
+
+    fast_winding_number(P,N,A,point_indices,children,expansion_order,CM,R,EC);
+    fast_winding_number(P,N,A,point_indices,children,CM,R,EC,Q,beta,WN);
+  }
+  
+  template <typename DerivedP, typename DerivedN,
+  	typename DerivedQ, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(
+                            const Eigen::MatrixBase<DerivedP>& P,
+                            const Eigen::MatrixBase<DerivedN>& N,
+                            const Eigen::MatrixBase<DerivedQ>& Q,
+                            Eigen::PlainObjectBase<DerivedWN>& WN
+                            ){
+    fast_winding_number(P,N,Q,2,2.0,WN);
+  }
 }
 
 
@@ -223,6 +364,8 @@ void fast_winding_number(const Eigen::MatrixXd & P,
 
 
 
+
+
 
 
 

include/igl/fast_winding_number.h

@@ -1,29 +1,197 @@
 #ifndef IGL_FAST_WINDING_NUMBER
 #define IGL_FAST_WINDING_NUMBER
+#include "igl_inline.h"
 #include <Eigen/Core>
 namespace igl
 {
-    void fast_winding_number_precompute(const Eigen::MatrixXd & P,
-                                        const Eigen::MatrixXd & N,
-                                        const Eigen::MatrixXd & A,
-                                        const std::vector<std::vector<int> > & point_indices,
-                                        const std::vector<Eigen::Matrix<int,8,1>, Eigen::aligned_allocator<Eigen::Matrix<int,8,1>>> & children,
-                                        Eigen::MatrixXd & CM,
-                                        Eigen::VectorXd & R,
-                                        Eigen::MatrixXd & EC
-                    );
-    
-    void fast_winding_number(const Eigen::MatrixXd & P,
-                             const Eigen::MatrixXd & N,
-                             const Eigen::MatrixXd & A,
-                             const std::vector<std::vector<int> > & point_indices,
-                             const std::vector<Eigen::Matrix<int,8,1>,               Eigen::aligned_allocator<Eigen::Matrix<int,8,1>>> & children,
-                             const Eigen::MatrixXd & CM,
-                             const Eigen::VectorXd & R,
-                             const Eigen::MatrixXd & EC,
-                             const Eigen::MatrixXd & Q,
-                             Eigen::VectorXd & WN
-                             );
+  // Generate the precomputation for the fast winding number for point data
+  // [Barill et. al 2018].
+  //
+  // Given a set of 3D points P, with normals N, areas A, along with octree
+  // data, and an expansion order, we define a taylor series expansion at each
+  // octree cell.
+  //
+  // The octree data is designed to come from igl::build_octree, and the areas
+  // (if not obtained at scan time), may be calculated using
+  // igl::point_areas_and_normals.
+  //
+  // Inputs:
+  //   P  #P by 3 list of point locations
+  //   N  #P by 3 list of point normals
+  //   A  #P by 1 list of point areas
+  //   point_indices  a vector of vectors, where the ith entry is a vector of
+  //                  the indices into P that are the ith octree cell's points
+  //   children       a vector of vectors, where the ith entry is a vector of
+  //                  the ith octree cell's of octree children
+  //   expansion_order    the order of the taylor expansion. We support 0,1,2.
+  // Outputs:
+  //   CM  #OctreeCells by 3 list of each cell's center of mass
+  //   R	 #OctreeCells by 1 list of each cell's maximum distance of any point
+  //			 to the center of mass
+  //   EC  #OctreeCells by #TaylorCoefficients list of expansion coefficients.
+  //	 		 (Note that #TaylorCoefficients = ∑_{i=1}^{expansion_order} 3^i)
+  //
+  template <typename DerivedP, typename DerivedA, typename DerivedN,
+  	typename Index, typename DerivedCM, typename DerivedR, typename DerivedEC>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedA>& A,
+            const std::vector<std::vector<Index> > & point_indices,
+            const std::vector<Eigen::Matrix<Index,8,1>, Eigen::aligned_allocator<Eigen::Matrix<Index,8,1>>> & children,
+            const int exansion_order,
+            Eigen::PlainObjectBase<DerivedCM>& CM,
+            Eigen::PlainObjectBase<DerivedR>& R,
+            Eigen::PlainObjectBase<DerivedEC>& EC);
+  
+  // Evaluate the fast winding number for point data, having already done the
+  // the precomputation
+  //
+  // Inputs:
+  //   P  #P by 3 list of point locations
+  //   N  #P by 3 list of point normals
+  //   A  #P by 1 list of point areas
+  //   point_indices  a vector of vectors, where the ith entry is a vector of
+  //                  the indices into P that are the ith octree cell's points
+  //   children       a vector of vectors, where the ith entry is a vector of
+  //                  the ith octree cell's of octree children
+  //   CM  #OctreeCells by 3 list of each cell's center of mass
+  //   R   #OctreeCells by 1 list of each cell's maximum distance of any point
+  //       to the center of mass
+  //   EC  #OctreeCells by #TaylorCoefficients list of expansion coefficients.
+  //        (Note that #TaylorCoefficients = ∑_{i=1}^{expansion_order} 3^i)
+  //   Q  #Q by 3 list of query points for the winding number
+  //   beta  This is a Barnes-Hut style accuracy term that separates near feild
+  //				 from far field. The higher the beta, the more accurate and slower
+  //         the evaluation. We reccommend using a beta value of 2. Note that
+  //				 for a beta value ≤ 0, we use the direct evaluation, rather than
+  //         the fast approximation
+  // Outputs:
+  //	 WN  #Q by 1 list of windinng number values at each query point
+  //
+  template <typename DerivedP, typename DerivedA, typename DerivedN,
+  	typename Index, typename DerivedCM, typename DerivedR, typename DerivedEC,
+  	typename DerivedQ, typename BetaType, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedA>& A,
+              const std::vector<std::vector<Index> > & point_indices,
+              const std::vector<Eigen::Matrix<Index,8,1>, Eigen::aligned_allocator<Eigen::Matrix<Index,8,1>>> & children,
+              const Eigen::MatrixBase<DerivedCM>& CM,
+              const Eigen::MatrixBase<DerivedR>& R,
+              const Eigen::MatrixBase<DerivedEC>& EC,
+              const Eigen::MatrixBase<DerivedQ>& Q,
+              const BetaType beta,
+              Eigen::PlainObjectBase<DerivedWN>& WN);
+  
+  // Evaluate the fast winding number for point data.
+  //
+  // This function performes the precomputation and evaluation all in one.
+  // If you need to acess the precomuptation for repeated evaluations, use the
+  // two functions designed for exposed precomputation (described above).
+  
+  // Inputs:
+  //   P  #P by 3 list of point locations
+  //   N  #P by 3 list of point normals
+  //   A  #P by 1 list of point areas
+  //   Q  #Q by 3 list of query points for the winding number
+  //   beta  This is a Barnes-Hut style accuracy term that separates near feild
+  //         from far field. The higher the beta, the more accurate and slower
+  //         the evaluation. We reccommend using a beta value of 2.
+  //   expansion_order    the order of the taylor expansion. We support 0,1,2.
+  // Outputs:
+  //   WN  #Q by 1 list of windinng number values at each query point
+  //
+  template <typename DerivedP, typename DerivedA, typename DerivedN,
+  	typename DerivedQ, typename BetaType, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedA>& A,
+                                      const Eigen::MatrixBase<DerivedQ>& Q,
+                                      const int expansion_order,
+                                      const BetaType beta,
+                                      Eigen::PlainObjectBase<DerivedWN>& WN
+                                      );
+  
+  // Evaluate the fast winding number for point data, with default expansion
+  // order and beta (both are set to 2).
+  //
+  // This function performes the precomputation and evaluation all in one.
+  // If you need to acess the precomuptation for repeated evaluations, use the
+  // two functions designed for exposed precomputation (described above).
+  //
+  // Inputs:
+  //   P  #P by 3 list of point locations
+  //   N  #P by 3 list of point normals
+  //   A  #P by 1 list of point areas
+  //   Q  #Q by 3 list of query points for the winding number
+  //   beta  This is a Barnes-Hut style accuracy term that separates near feild
+  //         from far field. The higher the beta, the more accurate and slower
+  //         the evaluation. We reccommend using a beta value of 2.
+  //   expansion_order    the order of the taylor expansion. We support 0,1,2.
+  // Outputs:
+  //   WN  #Q by 1 list of windinng number values at each query point
+  //
+  template <typename DerivedP, typename DerivedA, typename DerivedN,
+  	typename DerivedQ, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedA>& A,
+                                      const Eigen::MatrixBase<DerivedQ>& Q,
+                                      Eigen::PlainObjectBase<DerivedWN>& WN
+                                      );
+  
+  // Evaluate the fast winding number for point data, without known areas. The
+  // areas are calculated using igl::knn_search and
+  // igl::point_areas_and_normals.
+  //
+  // This function performes the precomputation and evaluation all in one.
+  // If you need to acess the precomuptation for repeated evaluations, use the
+  // two functions designed for exposed precomputation (described above).
+  
+  // Inputs:
+  //   P  #P by 3 list of point locations
+  //   N  #P by 3 list of point normals
+  //   Q  #Q by 3 list of query points for the winding number
+  //   beta  This is a Barnes-Hut style accuracy term that separates near feild
+  //         from far field. The higher the beta, the more accurate and slower
+  //         the evaluation. We reccommend using a beta value of 2.
+  //   expansion_order    the order of the taylor expansion. We support 0,1,2.
+  // Outputs:
+  //   WN  #Q by 1 list of windinng number values at each query point
+  //
+  template <typename DerivedP, typename DerivedN, typename DerivedQ,
+  	typename BetaType, typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedQ>& Q,
+                                      const int expansion_order,
+                                      const BetaType beta,
+                                      Eigen::PlainObjectBase<DerivedWN>& WN
+                                      );
+  
+  // Evaluate the fast winding number for point data, without known areas. The
+  // areas are calculated using igl::knn_search and
+  // igl::point_areas_and_normals. This function uses the default expansion
+  // order and beta (both are set to 2).
+  //
+  // This function performes the precomputation and evaluation all in one.
+  // If you need to acess the precomuptation for repeated evaluations, use the
+  // two functions designed for exposed precomputation (described above).
+  
+  // Inputs:
+  //   P  #P by 3 list of point locations
+  //   N  #P by 3 list of point normals
+  //   Q  #Q by 3 list of query points for the winding number
+  // Outputs:
+  //   WN  #Q by 1 list of windinng number values at each query point
+  //
+  template <typename DerivedP, typename DerivedN, typename DerivedQ,
+  	typename DerivedWN>
+  IGL_INLINE void fast_winding_number(const Eigen::MatrixBase<DerivedP>& P,
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedQ>& Q,
+                                      Eigen::PlainObjectBase<DerivedWN>& WN
+                                      );
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "fast_winding_number.cpp"

include/igl/knn_octree.cpp

@@ -1,112 +1,102 @@
-#include "point_areas_and_normals.h"
-#include <cmath>
-#include <igl/parallel_for.h>
-
-double distance_to_width_one_cube(Eigen::RowVector3d point){
-    return std::sqrt(std::pow(std::max(std::abs(point(0))-1,0.0),2)
-                     + std::pow(std::max(std::abs(point(1))-1,0.0),2)
-                     + std::pow(std::max(std::abs(point(2))-1,0.0),2));
-}
+#include "knn_octree.h"
+#include "parallel_for.h"
 
-double distance_to_cube(Eigen::RowVector3d point, Eigen::RowVector3d cube_center, double cube_width){
-    Eigen::RowVector3d transformed_point = (point-cube_center)/cube_width;
-    return cube_width*distance_to_width_one_cube(transformed_point);
-}
-
-//    template <typename DerivedP, typename DerivedI, typename DerivedO>
-//    IGL_INLINE void point_areas_and_normals(
-//                               const Eigen::MatrixBase<DerivedP>& P,
-//                               const Eigen::MatrixBase<DerivedI>& I,
-//                               const Eigen::MatrixBase<DerivedO>& O,
-//                               Eigen::MatrixBase<DerivedA> & A,
-//                               Eigen::MatrixBase<DerivedN> & N);
+#include <cmath>
 
 namespace igl {
-    void knn_octree(const Eigen::MatrixXd & P,
-                    const int & k,
-                    const std::vector<std::vector<int> > & point_indices,
-                    const std::vector<Eigen::Matrix<int,8,1>, Eigen::aligned_allocator<Eigen::Matrix<int,8,1>>> & children,
-                    const std::vector<Eigen::RowVector3d, Eigen::aligned_allocator<Eigen::RowVector3d>> & centers,
-                    const std::vector<double> & widths,
-                    Eigen::MatrixXi & I
-                    )
+  template <typename DerivedP, typename KType, typename IndexType,
+    typename CentersType, typename WidthsType, typename DerivedI>
+  IGL_INLINE void knn_octree(
+    const Eigen::MatrixBase<DerivedP>& P,
+    const KType & k,
+    const std::vector<std::vector<IndexType> > & point_indices,
+    const std::vector<Eigen::Matrix<IndexType,8,1>,
+        Eigen::aligned_allocator<Eigen::Matrix<IndexType,8,1> > > & children,
+    const std::vector<Eigen::Matrix<CentersType,1,3>,
+        Eigen::aligned_allocator<Eigen::Matrix<CentersType,1,3> > > & centers,
+    const std::vector<WidthsType> & widths,
+    Eigen::PlainObjectBase<DerivedI> & I)
+  {
+    
+    typedef Eigen::Matrix<typename DerivedP::Scalar, 1, 3> RowVector3PType;
+    
+    const int n = P.rows();
+    const KType real_k = std::min(n,k);
+    
+    auto distance_to_width_one_cube = [](RowVector3PType point){
+      return std::sqrt(std::pow(std::max(std::abs(point(0))-1,0.0),2)
+                       + std::pow(std::max(std::abs(point(1))-1,0.0),2)
+                       + std::pow(std::max(std::abs(point(2))-1,0.0),2));
+    };
+    
+    auto distance_to_cube = [&distance_to_width_one_cube]
+              (RowVector3PType point,
+               Eigen::Matrix<CentersType,1,3> cube_center,
+               WidthsType cube_width){
+      RowVector3PType transformed_point = (point-cube_center)/cube_width;
+      return cube_width*distance_to_width_one_cube(transformed_point);
+    };
+    
+    I.resize(n,real_k);
+    
+    igl::parallel_for(n,[&](int i)
     {
-        const int n = P.rows();
-        const int real_k = std::min(n,k);
-        I.resize(n,real_k);
-        
-      igl::parallel_for(n,[&](int i)
-      {
-          int points_found = 0;
-          Eigen::RowVector3d point_of_interest = P.row(i);
-          
-          //To make my priority queue take both points and octree cells, I use the indices 0 to n-1
-          //for the n points, and the indices n to n+m-1 for the m octree cells
-          
-          // Using lambda to compare elements.
-          auto cmp = [&point_of_interest, &P, &centers, &widths, &n](int left, int right) {
-              double leftdistance, rightdistance;
-              if(left < n){ //left is a point index
-                  leftdistance = (P.row(left) - point_of_interest).norm();
-              } else { //left is an octree cell
-                  leftdistance = distance_to_cube(point_of_interest,centers.at(left-n),widths.at(left-n));
-              }
-              
-              if(right < n){ //left is a point index
-                  rightdistance = (P.row(right) - point_of_interest).norm();
-              } else { //left is an octree cell
-                  rightdistance = distance_to_cube(point_of_interest,centers.at(right-n),widths.at(right-n));
-              }
-              return leftdistance >= rightdistance;
-          };
-          std::priority_queue<int, std::vector<int>, decltype(cmp)> queue(cmp);
-          queue.push(n); //This is the 0th octree cell (ie the root)
-          while(points_found < real_k){
-              int curr_cell_or_point = queue.top();
-              queue.pop();
-              if(curr_cell_or_point < n){ //current index is for is a point
-                  I(i,points_found) = curr_cell_or_point;
-                  points_found++;
-              } else {
-                  int curr_cell = curr_cell_or_point - n;
-                  if(children.at(curr_cell)(0) == -1){ //In the case of a leaf
-                      if(point_indices.at(curr_cell).size() > 0){ //Assumption: Leaves either have one point, or none
-                          queue.push(point_indices.at(curr_cell).at(0)); //push the point (pardon the pun)
-                      }
-                  } else { //Not a leaf
-                      for(int j = 0; j < 8; j++){
-                          queue.push(children.at(curr_cell)(j)+n); //+n to adjust for the octree cells
-                      }
-                  }
-              }
+      int points_found = 0;
+      RowVector3PType point_of_interest = P.row(i);
+      
+      //To make my priority queue take both points and octree cells,
+      //I use the indices 0 to n-1 for the n points,
+      // and the indices n to n+m-1 for the m octree cells
+      
+      // Using lambda to compare elements.
+      auto cmp = [&point_of_interest, &P, &centers, &widths,
+                  &n, &distance_to_cube](int left, int right) {
+        double leftdistance, rightdistance;
+        if(left < n){ //left is a point index
+          leftdistance = (P.row(left) - point_of_interest).norm();
+        } else { //left is an octree cell
+          leftdistance = distance_to_cube(point_of_interest,
+                                            centers.at(left-n),
+                                            widths.at(left-n));
+        }
+      
+        if(right < n){ //left is a point index
+          rightdistance = (P.row(right) - point_of_interest).norm();
+        } else { //left is an octree cell
+          rightdistance = distance_to_cube(point_of_interest,
+                                             centers.at(right-n),
+                                             widths.at(right-n));
+        }
+        return leftdistance >= rightdistance;
+      };
+      
+      std::priority_queue<IndexType, std::vector<IndexType>,
+        decltype(cmp)> queue(cmp);
+      
+      queue.push(n); //This is the 0th octree cell (ie the root)
+      while(points_found < real_k){
+        IndexType curr_cell_or_point = queue.top();
+        queue.pop();
+        if(curr_cell_or_point < n){ //current index is for is a point
+          I(i,points_found) = curr_cell_or_point;
+          points_found++;
+        } else {
+          IndexType curr_cell = curr_cell_or_point - n;
+          if(children.at(curr_cell)(0) == -1){ //In the case of a leaf
+            if(point_indices.at(curr_cell).size() > 0){
+              //Assumption: Leaves either have one point, or none
+              queue.push(point_indices.at(curr_cell).at(0));
+            }
+          } else { //Not a leaf
+            for(int j = 0; j < 8; j++){
+              //+n to adjust for the octree cells
+              queue.push(children.at(curr_cell)(j)+n);
+            }
           }
-//          points_found = 0;
-//          std::list<int> l = {n}; //This is the 0th octree cell (ie the root)
-//          while(points_found < real_k){
-//              int curr_cell_or_point = l.front();
-//              l.pop_front();
-//              if(curr_cell_or_point < n){ //current index is for is a point
-//                  I(i,points_found) = curr_cell_or_point;
-//                  points_found++;
-//              } else {
-//                  int curr_cell = curr_cell_or_point - n;
-//                  if(children.at(curr_cell)(0) == -1){ //In the case of a leaf
-//                      if(point_indices.at(curr_cell).size() > 0){ //Assumption: Leaves either have one point, or none
-//                          l.emplace_back(point_indices.at(curr_cell).at(0)); //push the point (pardon the pun)
-//                      }
-//                  } else { //Not a leaf
-//                      for(int j = 0; j < 8; j++){
-//                          l.emplace_back(children.at(curr_cell)(j)+n); //+n to adjust for the octree cells
-//                      }
-//                  }
-//                  l.sort(cmp);
-//                  int new_size = std::min<int>(l.size(),real_k-points_found);
-//                  l.resize(new_size);
-//              }
-//          }
-//
-      },1000);
-    }
+        }
+      }
+    },1000);
+  }
 }
 
 

include/igl/knn_octree.h

@@ -1,16 +1,47 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2018 Gavin Barill <gavinpcb@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_KNN_OCTREE
 #define IGL_KNN_OCTREE
+#include "igl_inline.h"
 #include <Eigen/Core>
+
 namespace igl
 {
-void knn_octree(const Eigen::MatrixXd & P,
-                const int & k,
-                const std::vector<std::vector<int> > & point_indices,
-                const std::vector<Eigen::Matrix<int,8,1>, Eigen::aligned_allocator<Eigen::Matrix<int,8,1>>> & children,
-                const std::vector<Eigen::RowVector3d, Eigen::aligned_allocator<Eigen::RowVector3d>> & centers,
-                const std::vector<double> & widths,
-                Eigen::MatrixXi & I
-                );
+  // Given a 3D set of points P, an whole number k, and an octree
+  // find the indicies of the k nearest neighbors for each point in P.
+  // Note that each point is its own neighbor.
+  //
+  // The octree data structures used in this function are intended to be the
+  // same ones output from igl::build_octree
+  //
+  // Inputs:
+  //   P  #P by 3 list of point locations
+  //	 k  number of neighbors to find
+  //   point_indices  a vector of vectors, where the ith entry is a vector of
+  //                  the indices into P that are the ith octree cell's points
+  //   children 			a vector of vectors, where the ith entry is a vector of
+  //									the ith octree cell's of octree children
+  //	 centers				a vector where the ith entry is a 3d row vector
+  //									representing the position of the ith cell's center
+  //   widths					a vector where the ith entry is the width of the ith
+  //									octree cell
+  // Outputs:
+  //   I  #P by k list of k-nearest-neighbor indices into P
+  template <typename DerivedP, typename KType, typename IndexType,
+    typename CentersType, typename WidthsType, typename DerivedI>
+  IGL_INLINE void knn_octree(const Eigen::MatrixBase<DerivedP>& P,
+    const KType & k,
+    const std::vector<std::vector<IndexType> > & point_indices,
+    const std::vector<Eigen::Matrix<IndexType,8,1>, Eigen::aligned_allocator<Eigen::Matrix<IndexType,8,1> > > & children,
+    const std::vector<Eigen::Matrix<CentersType,1,3>, Eigen::aligned_allocator<Eigen::Matrix<CentersType,1,3> > > & centers,
+    const std::vector<WidthsType> & widths,
+    Eigen::PlainObjectBase<DerivedI> & I);
 }
 #ifndef IGL_STATIC_LIBRARY
 #  include "knn_octree.cpp"