Patate Lib  0.5
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Pages
User Manual
Nicolas Mellado, Gautier Ciaudo, Gael Guennebaud, Pascal Barla


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>


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:
  • 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.

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
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; }
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
// Set the evaluation position

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++)
//finalize fitting

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


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

cout << "Center: [" << << "] ; 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;
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
// 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.