This documentation describes the tool to construct a binary image from a set of oriented point clouds in \(\mathbb{R}^3\) using the generalized winding number framework of [10]. For short, from a collection of oriented points, the generalized winding number implicit function returns a scalar value at each point accounting for the number of turns of the (unknown) underling surface around that point. From [10], the implicit function for a point cloud \( \{(p_i,n_i)\}_{1..m}\) is approximated by

\[ w(q) = \sum_i^m a_i \frac{(p_i-q)\cdot n_i}{4\pi\|p_i -q\|^3} \]

where \( a_i\) is some area contribution of the point \( p_i\) to the surface (computed from the projection of the k-nearest neighbors to \( p_i\) onto the tangent plane at \( p_i \)). By thresholding the winding function with some values in \([0,1]\), we can define an iso-surface that approximates the exterior envelope of the shape defined by the oriented point cloud. For fast computations, we rely on the libigl of the generalized winding numbers.

Warning: This class requires to have LIBIGL, CGAL and GMP cmake flags set to true (e.g. cmake .. -DWITH_LIBIGL!true -DWITH_CGAL=true -DWITH_GMP=true)

Basic usage

The core of the method is defined in the WindingNumberShape<Space> class. This class is a model of concepts::CEuclideanBoundedShape is the sense that once constructed, an object has an orientation( q ) method that returns the orientation (inside, outside or on) of q to the implicit surface (DGtal::Orientation).

#include <DGtal/helpers/StdDefs.h>
#include <DGtal/shapes/WindingNumberShape.h>
 
....
 
Eigen::MatrixXd &points  = { .... };  //mx3 matrix to encode the position of the m points
Eigen::MatrixXd &normals = { .... };  //mx3 matrix to encode the normal vector for each point
 
WindingNumbersShape<Z3i::Space> winding(points, normals);
 
Z3i::RealPoint q(0.3,0.4,0.0); //a single query point
 
DGtal::Orientation orientation = winding.orientation(q); //returns DGtal::INSIDE, DGtal::ON or DGtal::OUTSIDE

As the WindingNumberShape class is a model of concepts::CEuclideanBoundedShape, the GaussDigitizer could be used (see Shapes, Shapers and Digitizers).

Note: Note that for each query, we have to traverse a hierarchical data structure containing the input points (roughly \(O(\log m)\) expected time for balanced trees).

Example:

Input point cloud	Reconstruction h=2	Reconstruction h=1	Reconstruction h=0.2

Advanced usages: batched queries and modified area measures

First, for multiple queries and to take advantage of the multithreading of the winding number backend:

Eigen::MatrixXd &queries = { .... }; //Qx3 matrix to encode the query points location

std::vector<Orientation> orientations = winding.orientationBatch(queries);

The output vector contains the orientation of the query points in the same order.

At this point, we completely rely on libIGL and CGAL to estimate the \( a_i \) area measures. If the user wants to provide these quantities, the \(\{a_i\}\) can be prescribed at the constructor, or using the setAreaMeasure() method (with a skipPointAreas to true at the constructor).

Note: The orientation methods have a default thresholding parameter set to 0.3. Please check the class documentation for details.

Table of Contents

Basic usage

Advanced usages: batched queries and modified area measures