Voxelization from an oriented point cloud

Author

David Coeurjolly

Since

1.4

Part of Shapes package

See also

examples/polyscope_examples/windingNumberShape.cpp, tests/shapes/testWindingNumberShape.cpp

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 cmake flags set to true (e.g. cmake .. -DDGTAL_WITH_LIBIGL!true - DDGTAL_WITH_CGAL=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.