WindingNumbersShape.h

 
#pragma once
 
#if !defined WindingNumbersShape_h
#define WindingNumbersShape_h
 
#ifndef WITH_LIBIGL
#error You need to have activated LIBIGL (WITH_LIBIGL flag) to include this file.
#endif
 
#include <vector>
#include <algorithm>
#include <DGtal/base/Common.h>
#include <DGtal/base/CountedConstPtrOrConstPtr.h>
#include <DGtal/base/ConstAlias.h>
 
#include <DGtal/shapes/CEuclideanOrientedShape.h>
#include <igl/fast_winding_number.h>
#include <igl/octree.h>
#include <igl/knn.h>
#include <igl/copyleft/cgal/point_areas.h>
namespace DGtal
{
 
  template<typename TSpace>
  struct WindingNumbersShape
  {
    using Space       = TSpace;
    using RealPoint   = typename Space::RealPoint;
    using RealVector  = typename Space::RealVector;
    using Orientation = DGtal::Orientation;
    
    //Removing Default constructor
    WindingNumbersShape() = delete;
    
    WindingNumbersShape(ConstAlias<Eigen::MatrixXd> points, 
                        ConstAlias<Eigen::MatrixXd> normals)
    {
      myPoints  = points;
      myNormals = normals;
      myPointAreas = Eigen::VectorXd::Ones(myPoints->rows());
      // Build octree, from libIGL tutorials
      igl::octree(*myPoints,myO_PI,myO_CH,myO_CN,myO_W);
      if (points->rows()> 20)
      {
        Eigen::MatrixXi I;
        igl::knn(*myPoints,(int)points->rows(),myO_PI,myO_CH,myO_CN,myO_W,I);
        // CGAL is only used to help get point areas
        igl::copyleft::cgal::point_areas(*myPoints,I,*myNormals,myPointAreas);
        trace.info()<<" Min/max point area : "<<myPointAreas.minCoeff()<<" -- "<<myPointAreas.maxCoeff()<<std::endl;
      }
      else
      {
        trace.warning()<<"[WindingNumberShape] Too few points to use CGAL point_areas. Using the constant area setting."<<std::endl;
      }
    }
    
    WindingNumbersShape(ConstAlias<Eigen::MatrixXd> points,
                        ConstAlias<Eigen::MatrixXd> normals,
                        const Eigen::VectorXd &areas)
    {
      myPoints  = points;
      myNormals = normals;
      myPointAreas = areas;
      igl::octree(*myPoints,myO_PI,myO_CH,myO_CN,myO_W);
    }
    
    
    void setPointAreas(ConstAlias<Eigen::VectorXd> areas)
    {
      myPointAreas = areas;
    }
    
    Orientation orientation(const RealPoint aPoint,
                            const double threshold = 0.3) const
    {
      Eigen::MatrixXd queries(1,3);
      queries << aPoint(0) , aPoint(1) , aPoint(2);
      auto singlePoint = orientationBatch(queries, threshold);
      return singlePoint[0];
    }
    
    std::vector<Orientation> orientationBatch(const Eigen::MatrixXd & queries,
                                              const double threshold = 0.3) const
    {
      Eigen::VectorXd W;
      std::vector<Orientation> results( queries.rows() );
      Eigen::MatrixXd O_CM;
      Eigen::VectorXd O_R;
      Eigen::MatrixXd O_EC;
      
      //Checking if the areas
      igl::fast_winding_number(*myPoints,*myNormals,myPointAreas,myO_PI,myO_CH,2,O_CM,O_R,O_EC);
      igl::fast_winding_number(*myPoints,*myNormals,myPointAreas,myO_PI,myO_CH,O_CM,O_R,O_EC,queries,2,W);
      
      //Reformating the output
      for(auto i=0u; i < queries.rows(); ++i)
      {
        if (std::abs(W(i)) < threshold )
          results[i] = DGtal::OUTSIDE;
        else
          if (std::abs(W(i)) > threshold)
            results[i] = DGtal::INSIDE;
          else
            results[i] = DGtal::ON;
      }
      return results;
    }
    
    CountedConstPtrOrConstPtr<Eigen::MatrixXd> myPoints;
    CountedConstPtrOrConstPtr<Eigen::MatrixXd> myNormals;
    Eigen::VectorXd myPointAreas;
    
    
    std::vector<std::vector<int > > myO_PI;
    Eigen::MatrixXi myO_CH;
    Eigen::MatrixXd myO_CN;
    Eigen::VectorXd myO_W;
   
    
  };
}
 
#endif // !defined WindingNumbersShape_h