Merge pull request #810 from GavinBarill/master

Fast Winding Numbers for points

Former-commit-id: 382c3a911a786d122a6ff7a8d7da4570c1966c79

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_COPYLEFT_CGAL_FAST_WINDING_NUMBER
+#define IGL_COPYLEFT_CGAL_FAST_WINDING_NUMBER
+#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.cpp

@@ -0,0 +1,181 @@
+#include "point_areas.h"
+#include "delaunay_triangulation.h"
+
+#include "../../colon.h"
+#include "../../slice.h"
+#include "../../slice_mask.h"
+#include "../../parallel_for.h"
+
+#include "CGAL/Exact_predicates_inexact_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 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<DerivedN>& N,
+                                        Eigen::PlainObjectBase<DerivedA> & A,
+                                        Eigen::PlainObjectBase<DerivedT> & T)
+      {
+        typedef typename DerivedP::Scalar real;
+        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<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;
+        
+        
+        
+        const int n = P.rows();
+        
+        assert(P.cols() == 3 && "P must have exactly 3 columns");
+        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);
+        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(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);
+            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;
+          } 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();
+            
+            T.row(i) = svd.matrixV().col(2).transpose();
+            if(T.row(i).dot(N.row(i)) < 0){
+              T.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.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/fast_winding_number.cpp

@@ -0,0 +1,323 @@
+#include "fast_winding_number.h"
+#include "octree.h"
+#include "knn.h"
+#include "parallel_for.h"
+#include <vector>
+
+namespace igl {
+  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 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 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 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;
+    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;
+        }
+    };
+  
+    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::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::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;
+              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 =
+                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,1);
+  
+    std::function< real_wn(const RowVec, const std::vector<int>) > helper;
+    helper = [&helper,
+              &P,&N,&A,
+              &point_indices,&CH,
+              &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;
+      real_wn wn = 0;
+      for(int i = 0; i < near_indices.size(); i++){
+        int index = near_indices.at(i);
+        //Leaf Case, Brute force
+        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,
+                              N.row(curr_row)*A(curr_row));
+          }
+        }
+        //Non-Leaf Case
+        else {
+          for(int child = 0; child < 8; 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(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,
+                                        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(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);
+    }
+  }
+  
+  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;
+    
+    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;
+  
+    octree(P,point_indices,CH,CN,W);
+  
+    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,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,
+    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);
+  }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

include/igl/fast_winding_number.h

@@ -0,0 +1,122 @@
+#ifndef IGL_FAST_WINDING_NUMBER
+#define IGL_FAST_WINDING_NUMBER
+#include "igl_inline.h"
+#include <Eigen/Core>
+#include <vector>
+namespace igl
+{
+  // 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::octree, and the areas
+  // (if not obtained at scan time), may be calculated using
+  // igl::copyleft::cgal::point_areas.
+  //
+  // 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
+  //   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
+  //   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 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 Eigen::MatrixBase<DerivedCH>& CH,
+            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
+  //   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
+  //   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 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 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);
+  
+  // 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
+                                      );
+}
+#ifndef IGL_STATIC_LIBRARY
+#  include "fast_winding_number.cpp"
+#endif
+
+#endif
+

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
+