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>
 
 
//Removing Warnings from libIGL for gcc and clang
#pragma GCC system_header  // For GCC
#pragma clang system_header  // For Clang
 
#include <igl/fast_winding_number.h>
#include <igl/octree.h>
#include <igl/knn.h>
#include <igl/copyleft/cgal/point_areas.h>
 
#pragma GCC diagnostic pop  // For GCC
#pragma clang diagnostic pop  // For Clang
 
 
 
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,
                        bool skipPointAreas = false)
    {
      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 (skipPointAreas)
        trace.warning()<<"[WindingNumberShape] Skipping the CGAL point area estimation. By default, point areas are set to 1.0"<<std::endl;
      else
        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()<<"[WindingNumberShape] 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(), DGtal::OUTSIDE );
 
      rawWindingNumberBatch(queries, 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;
    }
    
 
    void rawWindingNumberBatch(const Eigen::MatrixXd & queries,
                               Eigen::VectorXd &W) const
    {
      Eigen::MatrixXd O_CM;
      Eigen::VectorXd O_R;
      Eigen::MatrixXd O_EC;
      
      //Computing the WN values
      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);
    }
    
    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