User Manual

Authors

Nicolas Mellado, Gautier Ciaudo, Gael Guennebaud, Pascal Barla

Introduction

This module is dedicated to the smooth fitting of point clouds and extraction of useful geometric properties. Figure 1(a) shows a typical example in 2D: we reconstruct a potential function (shown in fake colors) from a few 2D points equipped with normals; then the 0-isoline of this potential gives us an implicit curve (in orange) from which we can readily extract further properties like curvature. A great benefit of this implicit technique [3] is that it works in arbitrary dimensions: Figures 1(b-c) show how we reconstruct an implicit 3D surface with the same approach, starting from a 3D point cloud. Working with meshes then simply consists of only considering their vertices.

This is just the tip of the iceberg though, as we also provide methods for dealing with points equipped with non-oriented normals [2], techniques to analyze points clouds in scale-space to discover salient structures [4], methods to compute multi-scale principal curvatures [5] and methods to compute surface variation using a plane instead of a sphere for the fitting [7]. See the table below:

Primitive	Supported Input	Fitting techniques	Analysis/Tools	Other usages
Plane	Points	CovariancePlaneFit	Surface Variation[7]
Algebraic Sphere	Oriented points	OrientedSphereFit [3]	GLS[4]	Ray Traced Curvature[5]
Algebraic Sphere	Non-oriented points	UnorientedSphereFit [2]

In the following, we focus on a basic use of the module, using 3D points equipped with oriented normals fitted by a Grenaille::AlgebraicSphere. First, one must set up data samples that interface with an external code; then run the fitting process in itself, and finally collect outputs. We also show how to compute curvatures with the very same approach, as this is a common requirement of many geometric processing algorithms.

Figure 1: (a) An implicit 2D curve fit to a 2D point cloud. (b) A 3D point cloud shown with splats. (c) An implicit 3D surface reconstructed from (b).

Code structure

We structured Grenaille as follow: a "core" folder defining operators that rely on no data structure and work both with CUDA and C++; and an "algorithms" folder with methods that may interface with user-provided data structures, with optimized solutions in C++, CUDA or both. The rationale behind this separation is to provide two levels of complexity: core operators implement atomic scientific contributions that are agnostic of the host application; algorithms implement step-by-step procedures that will likely rely on core operators and/or structure traversal.

In its current state, we do not have algorithms for Grenaille implemented: it's work in progress (we are open to contributions). If you want to use Grenaille, just include its main header:

#include <Patate/grenaille.h>

Cuda

Grenaille can be used directly on GPU, thanks to several mechanisms:

• Eigen Cuda capabilities, see Eigen documentation for more details. You need to use a consistent Eigen::Index on both CPU and GPU if you plan to transfer memory between the computing units. That's why we recommend to set the following preprocessor variable when compiling your project:

  -DEIGEN_DEFAULT_DENSE_INDEX_TYPE=int

• C++03 Compliant: right now we do not use C++11 in Patate, to maximize compatibility across platforms and Cuda versions.
• Automatic CPU/GPU compilation qualifiers. We use the macro

  MULTIARCH void function();

  to use the same code for C++ and CUDA. It has no effect when the code is compiled with GCC or Clang, but it will force the compilation for both host and device architectures when compiling with nvcc. A similar macro system is provided for mathematical functions, to switch between STL and CUDA versions.

Check the C++/Cuda and Python/Cuda (using PyCuda) examples for more details and how-to.

You are now ready to read the next sections, and discover how to interface Grenaille with your data-structures, how to fit various primitives and compute geometrical and differential properties.

Data Samples

The first step needed to use Grenaille is to define how samples are represented inside the module. A strength of Grenaille is to define all these structures at compile time to generate optimized code for fast evaluation at runtime.

The class Grenaille::Concept::PointConcept defines the interface that has to be implemented to represent a sample. Observe that there is no need for data conversion: all you need to do is to indicate how to access existing data.

Warning

You should avoid data of low magnitude (i.e., 1 should be a significant value) to get good results; thus rescaling might be necessary.

As an example, let's fit a Grenaille::AlgebraicSphere onto points equipped with normals using the Grenaille::OrientedSphereFit. To this end, we must define a structure MyPoint containing a normal vector and its associated accessors. Depending on the fitting procedure we will use, we may need to define a Matrix Type. This is for instance required to compute Principal Curvatures. This leads to the following class:

using namespace Grenaille;

// This class defines the input data format
class MyPoint
{
public:
  enum {Dim = 3};
  typedef double Scalar;
  typedef Eigen::Matrix<Scalar, Dim, 1>   VectorType;
  typedef Eigen::Matrix<Scalar, Dim, Dim> MatrixType; // Type needed by Grenaille::GLSCurvatureHelper

  MULTIARCH inline MyPoint(const VectorType &pos    = VectorType::Zero(), 
         const VectorType& normal = VectorType::Zero())
    : _pos(pos), _normal(normal) {}
    
  MULTIARCH inline const VectorType& pos()    const { return _pos; }  
  MULTIARCH inline const VectorType& normal() const { return _normal; }

  MULTIARCH inline VectorType& pos()    { return _pos; }  
  MULTIARCH inline VectorType& normal() { return _normal; }

private:
  VectorType _pos, _normal;
};

To make the code more readable, we also define typedef helpers for Scalar and Vector types outside of the MyPoint class:

typedef MyPoint::Scalar Scalar;
typedef MyPoint::VectorType VectorType;

Fitting Process

Two template classes must be specialized to indicate how a fit will be applied.

The first step consists in identifying a weighting function: it defines how neighbor samples will contribute to the fit, as illustrated in Figure 2(a) in 2D. In this example, we choose a weight based on the Euclidean distance using Grenaille::DistWeightFunc, remapped through a bisquare kernel defined in Grenaille::SmoothWeightKernel:

typedef DistWeightFunc<MyPoint,SmoothWeightKernel<Scalar> > WeightFunc; 

The second step identifies a complete fitting procedure through the specialization of a Grenaille::Basket. In our example, we want to apply an Grenaille::OrientedSphereFit to input data points, which outputs an Grenaille::AlgebraicSphere by default. We further require such a sphere to be reparametrized using a Grenaille::GLSParam extension, which provides more intuitive parameters. This leads to the following specialization:

typedef Basket<MyPoint,WeightFunc,OrientedSphereFit, GLSParam> Fit1;  

At this point, most of the hard job has already been performed. All we have to do now is to provide an instance of the weight function, where t refers to the neighborhood size, and iniate the fit at an arbitrary position p.

// Create the previously defined fitting procedure
Fit1 fit;

// Set a weighting function instance
fit.setWeightFunc(WeightFunc(t));

// Set the evaluation position
fit.init(p);

Then neighbors are added sequentially: in this example, we traverse a simple array, and samples outside of the neighborhood are automatically ignored by the weighting function. Once all neighbors have been incorporated, the fit is performed and results stored in the specialized Basket object. STL-like iterators can be used directly for the fit by calling

fit.compute(vecs.begin(), vecs.end());

Internally, the container is traversed and the method finalize is called at the end:

// Iterate over samples and fit the primitive
for(vector<MyPoint>::iterator it = vecs.begin(); it != vecs.end(); it++)
  fit.addNeighbor(*it); 
  
//finalize fitting 
fit.finalize();

After calling finalize or compute, it is better to test the return state of the fitting before using it. There are diffent states but the most important one is STABLE.

FIT_RESULT eResult = fit.compute(vecs.begin(), vecs.end()); // or = fit.finalize();
if(eResult == STABLE) 
{
  //do things...
}

You may also use accessors to perform the same check:

fit.compute(vecs.begin(), vecs.end());

if(fit.isStable()) //You can also check the function isReady()
{
  //do things;
}

Figure 2. (a) Fitting a 2D point cloud of positions and normals at a point p (in red) requires to define a weight function of size t (in green). (b) This results in an implicit scalar field (in fake colors), from which parameters of a local spherical surface can be extracted: an offset tau, a normal eta and a curvature kappa.

Basic Outputs

Now that you have performed fitting, you may use its outputs in a number of ways (see Figure 2(b) for an illustration in 2D).

Scalar field

You may directly access properties of the fitted scalar-field, as defined in Grenaille::AlgebraicSphere :

cout << "Value of the scalar field at the initial point: " 
     << p.transpose() 
     << " is equal to " << fit.potential(p)
     << endl;
     
cout << "Its gradient is equal to: "
     << fit.primitiveGradient(p).transpose()
     << endl;

This generates the following output:

Value of the scalar field at the initial point: 0 0 0 is equal to -0.501162
Its gradient is equal to:   0.00016028  0.000178782 -0.000384989

Sphere

You may rather access properties of the fitted sphere (the 0-isosurface of the fitted scalar field), as defined in Grenaille::AlgebraicSphere :

cout << "Center: [" << fit.center().transpose() << "] ;  radius: " << fit.radius() << endl;

You will obtain:

Center: [-0.000160652 -0.000179197  0.000385884] ;  radius: 1.00232

Alternatively, you may prefer accessing parameters provided by the Grenaille::GLSParam extension:

cout << "Fitted Sphere: " << endl
     << "\t Tau  : "      << fit.tau()             << endl
     << "\t Eta  : "      << fit.eta().transpose() << endl
     << "\t Kappa: "      << fit.kappa()           << endl;

You will obtain:

Fitted Sphere: 
  Tau  : -0.501162
  Eta  :   0.35325  0.394028 -0.848502
  Kappa: 0.997682 

Other methods

Thanks to the Grenaille::AlgebraicSphere, you may also project an arbitrary point onto the fitted sphere via:

cout << "The initial point " << p.transpose()              << endl
     << "Is projected at   " << fit.project(p).transpose() << endl;

You will then obtain:

The initial point 0 0 0
Is projected at    0.353911  0.394765 -0.850088

Computing Curvatures

This part presents how to compute curvature values from the various fitting procedures and extensions available in Grenaille.

Mean Curvature

A fast approximation of the mean curvature is provided by GLSParam::kappa() and illustrated at two scales on a 3D surface in Figure 3:

typedef Basket<MyPoint, WeightFunc, OrientedSphereFit, // Sphere fitting
                                    GLSParam> Fit;     // GLS reparametrization

Fit fit;
// Fit primitive
// ... 

// Get mean curvature
Scalar curvature = fit.kappa();

It can also be extracted by analysing the spatial derivatives of fitted normals, provided by GLSDer::deta(). Note that the computational cost is more important in this case.

typedef Basket<MyPoint, WeightFunc, OrientedSphereFit,    // Sphere fitting
                                  GLSParam,               // GLS reparametrization
                                  OrientedSphereSpaceDer, // Spatial derivatives 
                                  GLSDer > Fit;           // GLS differentiation                                    

Fit fit;
// Fit primitive
// ... 

// Mean curvature values is given by half of the trace of the jacobian matrix
MatrixType jacobian = fit.deta();  
cout << " Mean curvature: " << 0.5 * jacobian.trace() << endl;

Note

In this example, we build the Jacobian from fit.deta(), using an implicit cast from VectorArray to MatrixType. It will require a more careful conversion when the basket contains different extensions. More details here.

Figure 3. Mean curvature computed at a fine (left) and a coarse (right) scale, and rendered with a simple color map (orange for concavities, blue for convexities).

Principal Curvatures

Principal curvatures and their associated directions can be computed by an eigen decomposition of the spatial derivatives of eta, provided by the extension GLSCurvatureHelper (based on the analysis of GLSDer::deta()):

typedef Basket<MyPoint, WeightFunc, OrientedSphereFit,    // Sphere fitting
                                  GLSParam,               // GLS reparametrization
                                  OrientedSphereSpaceDer, // Spatial derivatives 
                                  GLSDer,                 // GLS differentiation   
                                  GLSCurvatureHelper > Fit;           

Fit fit;
// Fit primitive
// ... 
// The eigen decomposition is called during the fit finalize
fit.finalize();

if(fit.isStable())
{
  // Get principal curvatures
  Scalar k1 = fit.GLSk1();
  Scalar k2 = fit.GLSk2();

  // Get associated directions
  VectorType d1 = fit.GLSk1Direction();
  VectorType d2 = fit.GLSk2Direction();

  // Get gaussian curvature
  Scalar K = fit.GLSGaussianCurvature();
}

Going Further

This page was intended to show you a standard use of the Grenaille module. However, there is more to it than Grenaille::OrientedSphereFit, as shown in its reference. For example, you can create a totally different Basket with a Grenaille::CompactPlane and its extension Grenaille::CovariancePlaneFit using plane instead of algebraic sphere to procede the fit in order to compute Grenaille::CovariancePlaneFit::surfaceVariation. We encourage you to make your own tests and have a look at our examples. It will give you ideas about how to use this module inside your own applications.