Made the changes Alec asked for

Former-commit-id: 2bf99911002a6b3fbb6c75027b1864317b7f161d

include/igl/copyleft/cgal/fast_winding_number.cpp

@@ -0,0 +1,82 @@
+#include "../../fast_winding_number.h"
+#include "../../octree.h"
+#include "../../knn.h"
+#include "../../parallel_for.h"
+#include "point_areas.h"
+#include <vector>
+
+namespace igl {
+  namespace copyleft{
+    namespace cgal{
+      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,Eigen::Dynamic,Eigen::Dynamic>
+                                                                  RealMatrix;
+        
+        std::vector<std::vector<int> > point_indices;
+        Eigen::Matrix<int,Eigen::Dynamic,8> CH;
+        Eigen::Matrix<real,Eigen::Dynamic,3> CN;
+        Eigen::Matrix<real,Eigen::Dynamic,1> W;
+        Eigen::MatrixXi I;
+        Eigen::Matrix<real,Eigen::Dynamic,1> A;
+        
+        octree(P,point_indices,CH,CN,W);
+        knn(P,21,point_indices,CH,CN,W,I);
+        point_areas(P,I,N,A);
+        
+        Eigen::Matrix<real,Eigen::Dynamic,Eigen::Dynamic> EC;
+        Eigen::Matrix<real,Eigen::Dynamic,3> CM;
+        Eigen::Matrix<real,Eigen::Dynamic,1> R;
+        
+        igl::fast_winding_number(P,N,A,point_indices,CH,
+                                 expansion_order,CM,R,EC);
+        igl::fast_winding_number(P,N,A,point_indices,CH,
+                                 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);
+      }
+    }
+  }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

include/igl/copyleft/cgal/fast_winding_number.h

@@ -0,0 +1,66 @@
+#ifndef IGL_FAST_WINDING_NUMBER_CGAL
+#define IGL_FAST_WINDING_NUMBER_CGAL
+#include "../../igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // Evaluate the fast winding number for point data, without known areas. The
+  // areas are calculated using igl::knn and igl::copyleft::cgal::point_areas.
+  //
+  // 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, which are the first two
+  // functions see in igl/fast_winding_number.h
+  //
+  // 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 and
+  // igl::point_areas. 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"
+#endif
+
+#endif
+

include/igl/copyleft/cgal/point_areas_and_normals.cpp → include/igl/copyleft/cgal/point_areas.cpp

@@ -1,4 +1,4 @@
-#include "point_areas_and_normals.h"
+#include "point_areas.h"
 #include "delaunay_triangulation.h"
 
 #include "../../colon.h"
@@ -22,24 +22,38 @@ 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(
+      
+      template <typename DerivedP, typename DerivedI, typename DerivedN,
+      typename DerivedA>
+      IGL_INLINE void point_areas(
+                                  const Eigen::MatrixBase<DerivedP>& P,
+                                  const Eigen::MatrixBase<DerivedI>& I,
+                                  const Eigen::MatrixBase<DerivedN>& N,
+                                  Eigen::PlainObjectBase<DerivedA> & A)
+      {
+        Eigen::MatrixXd T;
+        point_areas(P,I,N,A,T);
+      }
+      
+      
+      template <typename DerivedP, typename DerivedI, typename DerivedN,
+      typename DerivedA, typename DerivedT>
+      IGL_INLINE void point_areas(
                                         const Eigen::MatrixBase<DerivedP>& P,
                                         const Eigen::MatrixBase<DerivedI>& I,
-                                        const Eigen::MatrixBase<DerivedO>& O,
+                                        const Eigen::MatrixBase<DerivedN>& N,
                                         Eigen::PlainObjectBase<DerivedA> & A,
-                                        Eigen::PlainObjectBase<DerivedN> & N)
+                                        Eigen::PlainObjectBase<DerivedT> & T)
       {
         typedef typename DerivedP::Scalar real;
-        typedef typename DerivedO::Scalar scalarO;
+        typedef typename DerivedN::Scalar scalarN;
         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,
+        typedef Eigen::Matrix<scalarN, Eigen::Dynamic, Eigen::Dynamic> MatrixN;
+        typedef Eigen::Matrix<typename DerivedN::Scalar,
                   Eigen::Dynamic, Eigen::Dynamic> VecotorO;
         typedef Eigen::Matrix<typename DerivedI::Scalar,
                   Eigen::Dynamic, Eigen::Dynamic> MatrixI;
@@ -49,23 +63,23 @@ namespace igl {
         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() == N.rows()
+               && "P and N 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);
+        T.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 =
+          if(N.rows() && neighbors.rows() > 1){
+            MatrixN neighbor_normals;
+            igl::slice(N,neighbor_index,1,neighbor_normals);
+            Eigen::Matrix<scalarN,1,3> poi_normal = neighbor_normals.row(0);
+            Eigen::Matrix<scalarN,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);
@@ -74,7 +88,6 @@ namespace igl {
           }
           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
@@ -84,9 +97,9 @@ namespace igl {
                                     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;
+            T.row(i) = svd.matrixV().col(2).transpose();
+            if(T.row(i).dot(N.row(i)) < 0){
+              T.row(i) *= -1;
             }
             
             MatrixP plane;

include/igl/copyleft/cgal/point_areas.h

@@ -0,0 +1,78 @@
+// 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_H
+#define IGL_POINT_AREAS_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::octree. We reccomend using a k value
+    // between 15 and 20 inclusive for accurate area estimation.
+    //
+    // N is used filter the neighbours, 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 the tangent
+    // plane normal.
+    //
+    // Note: This function *should* be implemented by pre-filtering I, rather
+    // than filtering in this function using N. In this case, the function
+    // would only take P and I as input.
+    //
+    // Inputs:
+    //   P  #P by 3 list of point locations
+    //   I  #P by k list of k-nearest-neighbor indices into P
+    //   N  #P by 3 list of point normals
+    // Outputs:
+    //   A  #P list of estimated areas
+    template <typename DerivedP, typename DerivedI, typename DerivedN,
+      typename DerivedA>
+    IGL_INLINE void point_areas(
+                                        const Eigen::MatrixBase<DerivedP>& P,
+                                        const Eigen::MatrixBase<DerivedI>& I,
+                                        const Eigen::MatrixBase<DerivedN>& N,
+                                        Eigen::PlainObjectBase<DerivedA> & A);
+      
+    // This version can be used to output the tangent plane normal at each
+    // point. Since we area already fitting a plane to each point's neighbour
+    // set, the tangent plane normals come "for free"
+    //
+    // Inputs:
+    //   P  #P by 3 list of point locations
+    //   I  #P by k list of k-nearest-neighbor indices into P
+    //   N  #P by 3 list of point normals
+    // Outputs:
+    //   A  #P list of estimated areas
+    //   T  #P by 3 list of tangent plane normals for each point
+    template <typename DerivedP, typename DerivedI, typename DerivedN,
+    typename DerivedA, typename DerivedT>
+    IGL_INLINE void point_areas(
+                                const Eigen::MatrixBase<DerivedP>& P,
+                                const Eigen::MatrixBase<DerivedI>& I,
+                                const Eigen::MatrixBase<DerivedN>& N,
+                                Eigen::PlainObjectBase<DerivedA> & A,
+                                Eigen::PlainObjectBase<DerivedT> & T);
+      
+    }
+  }
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "point_areas.cpp"
+#endif
+
+#endif
+

include/igl/copyleft/cgal/point_areas_and_normals.h

@@ -1,59 +0,0 @@
-// 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,133 +1,131 @@
 #include "fast_winding_number.h"
-#include "build_octree.h"
-#include "knn_octree.h"
+#include "octree.h"
+#include "knn.h"
 #include "parallel_for.h"
-#include "copyleft/cgal/point_areas_and_normals.h"
 #include <vector>
 
 namespace igl {
-template <typename DerivedP, typename DerivedA, typename DerivedN,
-    typename Index, typename DerivedCM, typename DerivedR, typename DerivedEC>
+  template <typename DerivedP, typename DerivedA, typename DerivedN,
+  typename Index, typename DerivedCH, 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;
-    
-      typedef Eigen::Matrix<real_p,1,3> RowVec3p;
-    
-      int m = children.size();
-      int num_terms;
-    
-      assert(expansion_order < 3 && expansion_order >= 0 && "m must be less than n");
-      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;
-      }
-    
-      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;
-          zeroth_expansion << 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);
-          }
-        
-          masscenter = masscenter/areatotal;
-          CM.row(index) = masscenter;
-          EC.block(index,0,1,3) = zeroth_expansion;
-        
-          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::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::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(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);
-              }
-          }
-      };
-      helper(0);
+                                      const Eigen::MatrixBase<DerivedN>& N,
+                                      const Eigen::MatrixBase<DerivedA>& A,
+                                      const std::vector<std::vector<Index> > & point_indices,
+                                      const Eigen::MatrixBase<DerivedCH>& CH,
+                                      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;
+  
+    typedef Eigen::Matrix<real_p,1,3> RowVec3p;
+  
+    int m = CH.size();
+    int num_terms;
+  
+    assert(expansion_order < 3 && expansion_order >= 0 && "m must be less than n");
+    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;
+    }
+  
+    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,&CH,&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;
+        zeroth_expansion << 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);
+        }
+      
+        masscenter = masscenter/areatotal;
+        CM.row(index) = masscenter;
+        EC.block(index,0,1,3) = zeroth_expansion;
+      
+        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::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::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(index,12+9*k,1,9) += Eigen::Map<
+                      Eigen::Matrix<real_ec,1,9> >(TempCoeffs.data(),
+                                                   TempCoeffs.size());
+                }
+            }
+        }
+      
+        R(index) = max_norm;
+        if(CH(index,0) != -1)
+        {
+            for(int i = 0; i < 8; i++){
+                int child = CH(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>
+  typename Index, typename DerivedCH, 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
-            ){
+                        const Eigen::MatrixBase<DerivedN>& N,
+                        const Eigen::MatrixBase<DerivedA>& A,
+                        const std::vector<std::vector<Index> > & point_indices,
+                        const Eigen::MatrixBase<DerivedCH>& CH,
+                        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;
@@ -177,7 +175,8 @@ template <typename DerivedP, typename DerivedA, typename DerivedN,
               Eigen::Matrix<real_ec,3,3> RowCol;
               RowCol.setZero(3,3);
               RowCol.row(i) = loc;
-              RowCol = RowCol + RowCol.transpose();
+              Eigen::Matrix<real_ec,3,3> RowColT = RowCol.transpose();
+              RowCol = RowCol + RowColT;
               ThirdDerivative +=
                   -3.0/(4.0*M_PI*std::pow(r,5))*(RowCol+Diagonal);
               Eigen::Matrix<real_ec,1,9> derivative_vector =
@@ -196,7 +195,7 @@ template <typename DerivedP, typename DerivedA, typename DerivedN,
     std::function< real_wn(const RowVec, const std::vector<int>) > helper;
     helper = [&helper,
               &P,&N,&A,
-              &point_indices,&children,
+              &point_indices,&CH,
               &CM,&R,&EC,&beta,
               &direct_eval,&expansion_eval]
     (const RowVec query, const std::vector<int> near_indices)-> real_wn
@@ -206,7 +205,7 @@ template <typename DerivedP, typename DerivedA, typename DerivedN,
       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){
+        if(CH(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,
@@ -216,10 +215,10 @@ template <typename DerivedP, typename DerivedA, typename DerivedN,
         //Non-Leaf Case
         else {
           for(int child = 0; child < 8; child++){
-              int child_index = children.at(index)(child);
+              int child_index = CH(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){
+                  if(CH(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,
@@ -270,24 +269,20 @@ template <typename DerivedP, typename DerivedA, typename DerivedN,
                                       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;
+    Eigen::Matrix<int,Eigen::Dynamic,8> CH;
+    Eigen::Matrix<real,Eigen::Dynamic,3> CN;
+    Eigen::Matrix<real,Eigen::Dynamic,1> W;
   
-    build_octree(P,point_indices,children,centers,widths);
+    octree(P,point_indices,CH,CN,W);
   
-    Eigen::MatrixXd EC, CM;
-    Eigen::VectorXd R;
+    Eigen::Matrix<real,Eigen::Dynamic,Eigen::Dynamic> EC;
+    Eigen::Matrix<real,Eigen::Dynamic,3> CM;
+    Eigen::Matrix<real,Eigen::Dynamic,1> 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);
+    fast_winding_number(P,N,A,point_indices,CH,expansion_order,CM,R,EC);
+    fast_winding_number(P,N,A,point_indices,CH,CM,R,EC,Q,beta,WN);
   }
   
   template <typename DerivedP, typename DerivedA, typename DerivedN,
@@ -300,49 +295,6 @@ template <typename DerivedP, typename DerivedA, typename DerivedN,
                                       ){
     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);
-  }
 }
 
 

include/igl/fast_winding_number.h

@@ -12,9 +12,9 @@ namespace igl
   // 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
+  // The octree data is designed to come from igl::octree, and the areas
   // (if not obtained at scan time), may be calculated using
-  // igl::point_areas_and_normals.
+  // igl::copyleft::cgal::point_areas.
   //
   // Inputs:
   //   P  #P by 3 list of point locations
@@ -22,8 +22,8 @@ namespace igl
   //   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
+  //   CH             #OctreeCells by 8, where the ith row is the indices of
+  //                  the ith octree cell's 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
@@ -33,12 +33,13 @@ namespace igl
   //       (Note that #TaylorCoefficients = ∑_{i=1}^{expansion_order} 3^i)
   //
   template <typename DerivedP, typename DerivedA, typename DerivedN,
-    typename Index, typename DerivedCM, typename DerivedR, typename DerivedEC>
+    typename Index, typename DerivedCH, 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 Eigen::MatrixBase<DerivedCH>& CH,
             const int exansion_order,
             Eigen::PlainObjectBase<DerivedCM>& CM,
             Eigen::PlainObjectBase<DerivedR>& R,
@@ -53,8 +54,8 @@ namespace igl
   //   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
+  //   CH  #OctreeCells by 8, where the ith row is the indices of
+  //       the ith octree cell's 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
@@ -70,13 +71,14 @@ namespace igl
   //   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>
+    typename Index, typename DerivedCH, 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<DerivedCH>& CH,
             const Eigen::MatrixBase<DerivedCM>& CM,
             const Eigen::MatrixBase<DerivedR>& R,
             const Eigen::MatrixBase<DerivedEC>& EC,
@@ -84,35 +86,6 @@ namespace igl
             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).
   //
@@ -140,59 +113,6 @@ namespace igl
                                       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.cpp

@@ -0,0 +1,105 @@
+#include "knn.h"
+#include "parallel_for.h"
+
+#include <cmath>
+#include <queue>
+
+namespace igl {
+  template <typename DerivedP, typename KType, typename IndexType,
+  typename DerivedCH, typename DerivedCN, typename DerivedW,
+  typename DerivedI>
+  IGL_INLINE void knn(const Eigen::MatrixBase<DerivedP>& P,
+                      const KType & k,
+                      const std::vector<std::vector<IndexType> > & point_indices,
+                      const Eigen::MatrixBase<DerivedCH>& CH,
+                      const Eigen::MatrixBase<DerivedCN>& CN,
+                      const Eigen::MatrixBase<DerivedW>& W,
+                      Eigen::PlainObjectBase<DerivedI> & I)
+  {
+    typedef typename DerivedCN::Scalar CentersType;
+    typedef typename DerivedW::Scalar WidthsType;
+    
+    typedef Eigen::Matrix<typename DerivedP::Scalar, 1, 3> RowVector3PType;
+    
+    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)
+    {
+      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, &CN, &W,
+                  &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,
+                                            CN.row(left-n),
+                                            W(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,
+                                             CN.row(right-n),
+                                             W(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(CH(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(CH(curr_cell,j)+n);
+            }
+          }
+        }
+      }
+    },1000);
+  }
+}
+
+
+
+

include/igl/knn.h

@@ -0,0 +1,52 @@
+// 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
+#define IGL_KNN
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+namespace igl
+{
+  // 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::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
+  //   CH     #OctreeCells by 8, where the ith row is the indices of
+  //          the ith octree cell's children
+  //   CN     #OctreeCells by 3, where the ith row is a 3d row vector
+  //          representing the position of the ith cell's center
+  //   W      #OctreeCells, 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 DerivedCH, typename DerivedCN, typename DerivedW,
+    typename DerivedI>
+  IGL_INLINE void knn(const Eigen::MatrixBase<DerivedP>& P,
+    const KType & k,
+    const std::vector<std::vector<IndexType> > & point_indices,
+    const Eigen::MatrixBase<DerivedCH>& CH,
+    const Eigen::MatrixBase<DerivedCN>& CN,
+    const Eigen::MatrixBase<DerivedW>& W,
+    Eigen::PlainObjectBase<DerivedI> & I);
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "knn.cpp"
+#endif
+#endif
+

include/igl/octree.cpp

@@ -0,0 +1,176 @@
+#include "octree.h"
+#include <vector>
+#include <queue>
+
+namespace igl {
+  template <typename DerivedP, typename IndexType, typename DerivedCH,
+    typename DerivedCN, typename DerivedW>
+  IGL_INLINE void octree(const Eigen::MatrixBase<DerivedP>& P,
+                         std::vector<std::vector<IndexType> > & point_indices,
+                         Eigen::PlainObjectBase<DerivedCH>& CH,
+                         Eigen::PlainObjectBase<DerivedCN>& CN,
+                         Eigen::PlainObjectBase<DerivedW>& W)
+  {
+    
+    
+    
+    const int MAX_DEPTH = 30000;
+
+    typedef typename DerivedCH::Scalar ChildrenType;
+    typedef typename DerivedCN::Scalar CentersType;
+    typedef typename DerivedW::Scalar WidthsType;
+    typedef Eigen::Matrix<ChildrenType,8,1> Vector8i;
+    typedef Eigen::Matrix<typename DerivedP::Scalar, 1, 3> RowVector3PType;
+    typedef Eigen::Matrix<CentersType, 1, 3>       RowVector3CentersType;
+    
+    std::vector<Eigen::Matrix<ChildrenType,8,1>,
+        Eigen::aligned_allocator<Eigen::Matrix<ChildrenType,8,1> > > children;
+    std::vector<Eigen::Matrix<CentersType,1,3>,
+        Eigen::aligned_allocator<Eigen::Matrix<CentersType,1,3> > > centers;
+    std::vector<WidthsType> widths;
+    
+    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)){
+        index = index + 2;
+      }
+      if( location(2) >= center(2)){
+        index = index + 4;
+      }
+      return index;
+    };
+
+    
+    std::function< RowVector3CentersType(const RowVector3CentersType,
+                                         const CentersType,
+                                         const ChildrenType) >
+    translate_center =
+        [](const RowVector3CentersType & parent_center,
+           const CentersType h,
+           const ChildrenType 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;
+      }
+      RowVector3CentersType output = parent_center + change_vector;
+      return output;
+    };
+  
+    // 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();
+  
+    std::function< void(const ChildrenType, const int) > helper;
+    helper = [&helper,&translate_center,&get_octant,&m,
+              &zero_to_seven,&neg_ones,&P,
+              &point_indices,&children,&centers,&widths]
+    (const ChildrenType 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
+        CentersType h = widths.at(index)/2;
+        RowVector3CentersType curr_center = centers.at(index);
+        
+
+        for(ChildrenType 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<IndexType> all(P.rows());
+      for(IndexType i = 0;i<all.size();i++) all[i]=i;
+      point_indices.emplace_back(all);
+    }
+    children.emplace_back(neg_ones);
+  
+    //Get the minimum AABB for the points
+    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));
+    centers.emplace_back( aabb_center );
+  
+    //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);
+    
+    //Now convert from vectors to Eigen matricies:
+    CH.resize(children.size(),8);
+    CN.resize(centers.size(),3);
+    W.resize(widths.size(),1);
+    
+    for(int i = 0; i < children.size(); i++){
+      CH.row(i) = children.at(i);
+    }
+    for(int i = 0; i < centers.size(); i++){
+      CN.row(i) = centers.at(i);
+    }
+    for(int i = 0; i < widths.size(); i++){
+      W(i) = widths.at(i);
+    }
+  }
+}
+

include/igl/octree.h

@@ -0,0 +1,62 @@
+// 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_OCTREE
+#define IGL_OCTREE
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+
+
+
+
+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
+  //   CH     #OctreeCells by 8, where the ith row is the indices of
+  //          the ith octree cell's children
+  //   CN     #OctreeCells by 3, where the ith row is a 3d row vector
+  //          representing the position of the ith cell's center
+  //   W      #OctreeCells, a vector where the ith entry is the width
+  //          of the ith octree cell
+  //
+  template <typename DerivedP, typename IndexType, typename DerivedCH,
+  typename DerivedCN, typename DerivedW>
+  IGL_INLINE void octree(const Eigen::MatrixBase<DerivedP>& P,
+    std::vector<std::vector<IndexType> > & point_indices,
+    Eigen::PlainObjectBase<DerivedCH>& CH,
+    Eigen::PlainObjectBase<DerivedCN>& CN,
+    Eigen::PlainObjectBase<DerivedW>& W);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "octree.cpp"
+#endif
+
+#endif
+