Recipe for implicit surface reconstruction with HRBF

In this entry I'll explain and give code for the easiest method I know to reconstruct a surface represented by a scalar field \(f: \mathbb R^3 \rightarrow \mathbb R\) by interpolating a set of point with their normals. The method is called Hermite Radial Basis Function (HRBF) interpolation. The audience aimed is anyone with some basic knowledge about implicit surfaces (which people often know as metaballs and render with the marching cube algorithm)

Quick links:


An example

The goal is to build a surface \( f \) going through \( N \) points \( \mathbf p_i \), moreover the normals of \( f \) must match the user defined normals at each point \( \mathbf p_i \):

In these figures we can see two configurations. Here, points coordinates are the same, only a normal changes its direction (they are both represented in green). As we can see it plays an important role in the final shape of the surface. Using normals in addition to points has several advantages: the underlying algorithm will be more robust and easier to implement unlike other methods who only uses points (like standard RBF reconstruction). Also normals gives more control to the user like tangents used for 2D curves.

Other examples: look at the introductory figures, we see that when a point is too far away from the other points, a piece of surface gets separated from the main object producing a second piece of surface shaped like a bubble.

The maths

I'll split this section into two parts. First part (Conventions) I'll describe the conventions and basic notions I'm using, it's the easy part. The second part (Building the scalar field \( f \)) is more difficult but not mandatory to use the C++ code available below, however it will be usefull to anyone who want to re-implement the method.


So, given a set of \( N \) points and normals \( (\mathbf p_i, \mathbf n_i) \) we seek a scalar field \(f: \mathbb R^3 \rightarrow \mathbb R\). By convention we want \( f  \) to be equal to zero everywhere the surface lies. Since the surface must go through every points \( \mathbf p_i \) then \( f(\mathbf p_i) = 0\). Inside the shape \( f \) is negative and outside \( f \) is positive. The function \( f \) returns values ranging from \( [-\infty; +\infty ] \) and has increasing values from the inside to the outside. We say that \( f \) has a global support or is globally supported. It means the function varies everywhere in the ambiant space \( \mathbf R^3 \).
It's worth noticing Metaballs are not defined with globally supported functions but compactly supported functions. Indeed the scalar field of a metaball decreases from the center untill it vanishes to zero at the metaball's radius. So the function of a metaball stops varying outside its radius which means we can bound the variation inside a bounding box hence the name compactly supported. If you want to blend objects by summing them as you usualy do with metaballs they have to be compactly supported. In this other entry I explain how to convert a globally supported scalar field to a compact support.

Here is a picture to help you visualize our function \( f \):

We see the surface for \(f(\mathbf p) \) = 0. Curves inside are red with \( f(\mathbf p) < 0\) and outside blue with \( f(\mathbf p) > 0\). What I did is cut the 3D space with a 2D plane and draw at regular interval curves. Every points that lies on the same curve  has the same scalar value of \( f \). For instance points at the frontier between blue and red are equal to zero. The next curve farther from the surface represent points for which \( f(\mathbf p) = 1 \) and so on. These curves are analog to contour lines in hiking maps but instead of representing the terrain elevation it tells the distance between a point and the implicit surface. (N.B: blue curves should continue indefinitely as the distance increases but I did not render all of them)

There is one thing I haven't defined yet which is how to compute the normals to the implicit surface. We have to compute what we call the gradient of the implicit surface which is written \( \nabla f \). The gradient return a vector which contains in each component the partial derivatives of \( f \):

$$ \nabla f = \begin{pmatrix}  \frac {\partial f(x, y, z)}{\partial x}, && \frac {\partial f(x, y, z)}{\partial y}, && \frac {\partial f(x, y, z)}{\partial z}
\end{pmatrix} $$

The gradient vector gives at any point of \( f \) the direction of steepest ascent i.e. the direction towards which \(f\) increases most. If you look back to the iso-curves you'll understand that the gradient is always orthogonal to the curves and hence to the implicit surface.  Indeed \( f \) doesn't vary if we follow a curve but increases the most when we go right through them.

Now we know how to compute the normal with \( \nabla f \) at any point of the implicit surface. In order to the implicit surface normals to match normals \( n_i \), we have to set \( \nabla f(\mathbf p_i) = \mathbf n_i \). This is what we have defined so far:

Building the scalar field \( f \)

Note: reading this section is not mandatory to use the C++ code.

Now comes the hard part: finding \( f \) that interpolates the set of points and normals \( (\mathbf p_i, \mathbf n_i) \). First we have to define what \( f \) looks like. How to find the form of \( f \) is outside the topic of this tutorial I'll just give the formula directly:

$$ f(\mathbf{x}) = \sum_i^N \mathbf{\alpha}_i \phi(\| \mathbf{x} - \mathbf{p}_i\|) +
\boldsymbol{\beta}_i . \nabla \left [\phi(\| \mathbf{x} - \mathbf{p}_i \|) \right ]$$

So with have \( \mathbf x \in \mathbb R^3 \) the evaluated position and \( \mathbf{p}_i \) the samples/points to be interpolated. The function \( \phi : \mathbb R \rightarrow \mathbb R \) is a called a radial basis function which I advice to define as \( \phi(x) = x^3 \). there are other possibilities but in my experience the method works best with \( x^3\). In order to be able to compute the value of \( f \) we will have to find the \( N \) values of \( \alpha_i \in \mathbb R \) and \( \boldsymbol{\beta}_i \in \mathbb R^3 \). Notice that \( \boldsymbol{\beta}_i \) and \( \nabla \phi \) being vectors \( \boldsymbol{\beta}_i . \nabla \phi \) is the dot product. In order to find the unknown \( \alpha_i \) and \( \boldsymbol{\beta}_i \) we have to solve this linear system of equations:

$$ \begin{pmatrix}
 f( \mathbf{p}_i ) \\
 \nabla f( \mathbf{p}_i )
  0 \\
\end{pmatrix} $$

As we said before points on the surface are associated to the value zero and normals and gradient must match. For this system I'm using a rather compact notation because there is actually \( 4N \) lines with \( 4N \) unknown. As it would take too much space I have put all the [  devellopements details in a pdf ]. If you're not interested in the devellopments I've made a [  summary with only the detailed matrix ] to solve the system.

Interpreting the formula of \( f \)

It may be a bit abstract to understand how summing the expression \( \mathbf{\alpha}_i \phi(\| \mathbf{x} - \mathbf{p}_i\|) +
\boldsymbol{\beta}_i . \nabla \left [\phi(\| \mathbf{x} - \mathbf{p}_i \|) \right ] \) with appropriate weights \( \mathbf{\alpha}_i \) and \( \boldsymbol{\beta}_i \) gives a correct interpolation of the implicit surface. Lets take an example and draw the implicit surface:

$$ g(\mathbf x) = 1 \phi(\| \mathbf{x} \|) + (1,0,0) . \nabla \left [\phi(\| \mathbf{x} \|) \right ]  $$

It's a single HRBF centered at \( (0, 0, 0) \) with \( \mathbf{\alpha}_i = 1 \) and \( \boldsymbol{\beta}_i = (1, 0, 0) \) which looks like this:

The point position in red is \( (0,0,0) \) and the normal direction is \( (1,0,0) \). So we can observe that each point correspond to a blob like function. However the blob is not centered at the point position but its surface goes through the point. Moreover the normal controls the orientation of this blob. So the intuition is that by blending all the blobs with a sum we generate a final object whose shape goes through all the points.


The essential

The C++ code for HRBF interpolation uses the librarie Eigen to solve the linear system. Fortunately it's a header library. It means you won't have to compile/install it because it is only composed of headers. They are readily usable within any projects using the #include directive. The [  HRBF C++ code ] should not be hard to use. The class takes in input a vector of points and normals to compute the surface as shown here:

#include "hrbf_phi_funcs.h"
#include "hrbf_core.h"

typedef HRBF_fit<float, 3, Rbf_pow3<float> > HRBF;

    HRBF fit;
    // Define samples (position, normal);
    HRBF::Vector pts[] = { HRBF::Vector( 0.f, 0.f, 0.f), HRBF::Vector( 1.f, 0.f, 0.f), HRBF::Vector( 0.f, 0.f, 2.f) };
    HRBF::Vector nrs[] = { HRBF::Vector(-1.f, 0.f, 0.f), HRBF::Vector( 1.f, 0.f, 0.f), HRBF::Vector( 0.f, 0.f, 1.f) };

    const int size = sizeof(pts) / sizeof(HRBF::Vector);
    std::vector p(pts, pts + size );
    std::vector n(nrs, nrs + size );
    // Solve linear system;
    fit.hermite_fit(p, n);

Then the implicit surface potential can be easily evaluated at any point of the ambiant space:

HRBF::Vector x( 0.f, 0.f, 0.f);
float potential = fit.eval( x );
HRBF::Vector gf = fit.grad( x );

Note also that the code works for any dimensions, for instance in 2D you can easily generate the function of an implicit curve \( f: \mathbb R^2 \rightarrow \mathbb R \)

typedef HRBF_fit<float, 2, Rbf_pow3<float> > HRBF_curve;


Be aware compiling in Debug mode is a lot slower than Release mode to solve linear systems.

The toy application

In order to do some testing rapidly I provide a toy application. It enables the edition of HRBF points and the vizualisation with a marching cube algorithm. There is also an automatic bounding box fitting around the HRBF. Feel free to play and modify the program:


Hermite Interpolation of Implicit Surfaces with Radial Basis Functions; Ives Macêdo & Al

five comments

An impressive work, keep going :D

URIEN Loic - 18-06-’13 19:46

Great job, learning… :)

flyiran - 22-10-’13 11:00

By any chance are you developing a Blender plug-in?
Rodolphe: Nop sorry just the standalone app.

H. Rhynedahll - 29-11-’13 19:23

Nice write up. :)

Ives M. - 16-04-’15 19:41

Nice, thanks. Have you any more info, about selection of basis function as x^3 ?
Rodolphe: not really it comes from my personal experience, I tried others (polynomials, thin plates etc.) but among the ones I tried x^3 performed the best.

Perry - 13-11-’15 11:33
All html tags except <b> and <i> will be removed from your comment. You can make links by just typing the url or mail-address.
Spam bot question: