DGtal  1.5.beta
WindingNumbersShape.h
1 
17 #pragma once
18 
31 #if !defined WindingNumbersShape_h
33 #define WindingNumbersShape_h
34 
35 #ifndef WITH_LIBIGL
36 #error You need to have activated LIBIGL (WITH_LIBIGL flag) to include this file.
37 #endif
38 
39 #include <vector>
40 #include <algorithm>
41 #include <DGtal/base/Common.h>
42 #include <DGtal/base/CountedConstPtrOrConstPtr.h>
43 #include <DGtal/base/ConstAlias.h>
44 #include <DGtal/shapes/CEuclideanOrientedShape.h>
45 
46 
47 //Removing Warnings from libIGL for gcc and clang
48 #pragma GCC system_header // For GCC
49 #pragma clang system_header // For Clang
50 
51 #include <igl/fast_winding_number.h>
52 #include <igl/octree.h>
53 #include <igl/knn.h>
54 #include <igl/copyleft/cgal/point_areas.h>
55 
56 #pragma GCC diagnostic pop // For GCC
57 #pragma clang diagnostic pop // For Clang
58 
59 
60 
61 namespace DGtal
62 {
64 
76  template<typename TSpace>
78  {
80  using Space = TSpace;
81  using RealPoint = typename Space::RealPoint;
82  using RealVector = typename Space::RealVector;
84 
85  //Removing Default constructor
86  WindingNumbersShape() = delete;
87 
99  bool skipPointAreas = false)
100  {
101  myPoints = points;
102  myNormals = normals;
103  myPointAreas = Eigen::VectorXd::Ones(myPoints->rows());
104  // Build octree, from libIGL tutorials
105  igl::octree(*myPoints,myO_PI,myO_CH,myO_CN,myO_W);
106  if (skipPointAreas)
107  trace.warning()<<"[WindingNumberShape] Skipping the CGAL point area estimation. By default, point areas are set to 1.0"<<std::endl;
108  else
109  if (points->rows()> 20)
110  {
111  Eigen::MatrixXi I;
112  igl::knn(*myPoints,(int)points->rows(),myO_PI,myO_CH,myO_CN,myO_W,I);
113  // CGAL is only used to help get point areas
114  igl::copyleft::cgal::point_areas(*myPoints,I,*myNormals,myPointAreas);
115  trace.info()<<"[WindingNumberShape] Min/max point area : "<<myPointAreas.minCoeff()<<" -- "<<myPointAreas.maxCoeff()<<std::endl;
116  }
117  else
118  {
119  trace.warning()<<"[WindingNumberShape] Too few points to use CGAL point_areas. Using the constant area setting."<<std::endl;
120  }
121  }
122 
131  const Eigen::VectorXd &areas)
132  {
133  myPoints = points;
134  myNormals = normals;
135  myPointAreas = areas;
136  igl::octree(*myPoints,myO_PI,myO_CH,myO_CN,myO_W);
137  }
138 
139 
143  {
144  myPointAreas = areas;
145  }
146 
156  const double threshold = 0.3) const
157  {
158  Eigen::MatrixXd queries(1,3);
159  queries << aPoint(0) , aPoint(1) , aPoint(2);
160  auto singlePoint = orientationBatch(queries, threshold);
161  return singlePoint[0];
162  }
163 
170  std::vector<Orientation> orientationBatch(const Eigen::MatrixXd & queries,
171  const double threshold = 0.3) const
172  {
173  Eigen::VectorXd W;
174  std::vector<Orientation> results( queries.rows(), DGtal::OUTSIDE );
175 
176  rawWindingNumberBatch(queries, W);
177 
178  //Reformating the output
179  for(auto i=0u; i < queries.rows(); ++i)
180  {
181  if (std::abs(W(i)) < threshold )
182  results[i] = DGtal::OUTSIDE;
183  else
184  if (std::abs(W(i)) > threshold)
185  results[i] = DGtal::INSIDE;
186  else
187  results[i] = DGtal::ON;
188  }
189  return results;
190  }
191 
192 
197  void rawWindingNumberBatch(const Eigen::MatrixXd & queries,
198  Eigen::VectorXd &W) const
199  {
200  Eigen::MatrixXd O_CM;
201  Eigen::VectorXd O_R;
202  Eigen::MatrixXd O_EC;
203 
204  //Computing the WN values
205  igl::fast_winding_number(*myPoints,*myNormals,myPointAreas,myO_PI,myO_CH,2,O_CM,O_R,O_EC);
206  igl::fast_winding_number(*myPoints,*myNormals,myPointAreas,myO_PI,myO_CH,O_CM,O_R,O_EC,queries,2,W);
207  }
208 
214  Eigen::VectorXd myPointAreas;
215 
217  std::vector<std::vector<int > > myO_PI;
219  Eigen::MatrixXi myO_CH;
221  Eigen::MatrixXd myO_CN;
223  Eigen::VectorXd myO_W;
224 
225 
226  };
227 }
228 
229 #endif // !defined WindingNumbersShape_h
Aim: This class encapsulates its parameter class so that to indicate to the user that the object/poin...
Definition: ConstAlias.h:187
std::ostream & info()
std::ostream & warning()
Space::RealVector RealVector
DGtal is the top-level namespace which contains all DGtal functions and types.
Trace trace
Definition: Common.h:153
Orientation
Definition: Common.h:141
@ INSIDE
Definition: Common.h:141
@ OUTSIDE
Definition: Common.h:141
@ ON
Definition: Common.h:141
Aim: model of a CEuclideanOrientedShape from an implicit function from an oriented point cloud....
Eigen::VectorXd myPointAreas
Const alias to point area measure.
Eigen::MatrixXd myO_CN
libIGL octree for fast queries data structure
Eigen::VectorXd myO_W
libIGL octree for fast queries data structure
Orientation orientation(const RealPoint aPoint, const double threshold=0.3) const
typename Space::RealPoint RealPoint
std::vector< Orientation > orientationBatch(const Eigen::MatrixXd &queries, const double threshold=0.3) const
void setPointAreas(ConstAlias< Eigen::VectorXd > areas)
CountedConstPtrOrConstPtr< Eigen::MatrixXd > myPoints
Const alias to the points.
Eigen::MatrixXi myO_CH
libIGL octree for fast queries data structure
typename Space::RealVector RealVector
WindingNumbersShape(ConstAlias< Eigen::MatrixXd > points, ConstAlias< Eigen::MatrixXd > normals, const Eigen::VectorXd &areas)
WindingNumbersShape(ConstAlias< Eigen::MatrixXd > points, ConstAlias< Eigen::MatrixXd > normals, bool skipPointAreas=false)
std::vector< std::vector< int > > myO_PI
libIGL octree for fast queries data structure
void rawWindingNumberBatch(const Eigen::MatrixXd &queries, Eigen::VectorXd &W) const
CountedConstPtrOrConstPtr< Eigen::MatrixXd > myNormals
Const alias to the normals.
const Point aPoint(3, 4)
PointVector< 3, double > RealPoint