Last year I was involved in the development of a Pose Space Deformation system. One interesting key-component of this system is the global interpolation mechanism between vertices. In this article I will quickly go over the fundamentals of interpolation and present the RBF interpolation algorithm.
Interpolation 101
Interpolation is a cornerstone of Computer Graphics as it enables to treat discrete samples as continuous data. The idea is to evaluate an unknown function f at an arbitrary point x. The value of f is known for a number of sample points x0, x1, .., xn-1. A sample point p is a n-dimensional point associated with a n-dimensional value f(p). For example, you can consider the pixels of a RGBA texture like samples which are 2D points associated to 4D values. Interpolation can be either global or local. In global interpolation each sample is taken into account to compute any interpolated value, whereas in local interpolation only the neighbor points are used to compute interpolation.
Linear interpolation or lerp
Linear interpolation is used to approximate 1-variable functions, or in more simple terms, to trace a straight line between two points. Representing lerp with weight between samples is the approach used in GLSL via the mix function. Linear interpolation on multiple samples will be computed using two samples at a time and is obviously local.
y = f(t, y0, y1) = t * y0 + (1-t) * y1 where 0 <= t <=1 and t = (x1 – x0) / (x1 – x)
Bilinear interpolation
This one is used in many applications, notability in hardware texture mapping. Bilinear interpolation is the extension of linear interpolation for 2-variable functions e.g. 2D points. You will use it to smoothly interpolate the value of samples on a 2D surface. In 3D this interpolation will be called trilinear interpolation. It takes place on a regular grid. For a given point p = (x,y), you’ll need the four surrounding samples. The main idea is to perform a cheap quadratic interpolation using only two lerps, one on each axis.
We define a 2-variables function f, and 4 2D samples S00 S01 S10 S11 associated with a 1D value (be that a hue, a height or something else). You can easily compute f(x,y) using the lerp formula (meaning 0 <= x <=1 and 0<= y <= 1) with :
If you plot this using the value of f(x,y) either as a hue or a height you got something like this :
Polynomial interpolation
Polynomial interpolation is global and consists in finding for a set of n samples, a polynomial p of degree at most n-1 , where p(xi) = yi (in the case of a single variable polynomial). In practice, polynomial interpolation can be used to trace a curve passing by all points. Linear interpolation is a polynomial interpolation of degree 1.
To build the polynomial, it is necessary to build the given set of linear equations, each line representing a sample and its value. Such a system can be resolved with either Gaussian elimination or SVD, but more on that later. Once the system is resolved the coefficients S represent the coefficients of the polynomial.
Polynomial interpolation is impractical for large values of n. The much more used spline interpolation, uses lower-degrees polynomials to interpolate between a subset of the samples.
Interpolation of scattered data in n-dimensions
Samples are not always nicely ordered along a regular grid. RBF interpolation is used to interpolate globally between samples scattered in the n-dimensional space. As the interpolation system is built upon distances between samples, the method can be used with any set of n-dimensional samples, such as 3D points, 2D points or quaternions.
Radial basis functions
RBF are functions whose values depend only on the distance d of a point p to the center c. Most of RBFs take one additional parameter, for example the Gaussian function takes the width of the curve as a parameter. The huge number of RBF available and their various parameters, make the RBF interpolation very flexible.
- Gaussian function (orange)
- Multiquadric (red)
- Inverse Quadratic (magenta)
- Thin plate (brown)
The algorithm
The algorithm is quite simple, it consists again in building a set of linear equations and solving it.
- Build a (n,n) matrix. Each cell contains the value of the RBF of the distance between two samples
- Solve the system, store the coefficients S
- To find the interpolated value at p, take the sum of the RBF of the distance between p and each sample multiplied by the coefficients.
A bit of code now.
First the Gaussian RBF :
1 | double RBFGauss(double d, double e) { return exp(-(d*e)*(d*e)); } |
Second, building the system and solving it :
Numerical Recipes recommends using the LU decomposition to inverse the matrix and solve the system. I ended up using a SVD instead, mainly because it is faster and it handles a lot of degenerated cases. Since the matrix is symmetrical, you could use something like Cholesky decomposition (which is a special case of LU). For simplicity, I used 2D samples and 1D values.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | bool buildWeights(int numSamples, const Vector2d * s, const double * v, double * w) { MatrixXd m(numSamples, numSamples); // The system VectorXd values(numSamples); VectorXd weights(numSamples); for (int i = 0; i < numSamples; ++i) { values[i] = v[i]; for (int j = 0; j < numSamples; ++j) { m(i,j) = RBFGauss((s[i]-s[j]).norm() ,EPSILON); } } SVD svd(m); bool status = svd.solve(values, &weights); memcpy(w, weights.data(), numSamples * sizeof(double)); return status; } |
Third, the interpolation itself, as you can see it is quite straightforward :
1 2 3 4 5 6 7 8 9 | double interpolate(int numSamples, const Vector2d & p, const Vector2d * s, const double * w) { double res = 0.0; for (int i = 0; i < numSamples; ++i) { res += w[i] * RBFGauss((p-s[i]).norm(), EPSILON); } return res; } |
Keep in mind that the code in the snippets is clarity-oriented, there are tons of optimizations you can do, notably if your values are multidimensional.
A few 3D plots to finish (click on the images, the thumbnail engine ate my samples). Here the samples (in red) are 2D (x,z) and are associated to a height (y). The height is interpolated along a 2D surface between the samples. You can notice the behaviour when the interpolation takes place further away from the samples.
When using 3D points both as positions and values one can generate interesting shapes. Here (x,y,z) values are interpolated along (x,y).
In production : Position Space Deformation
I won’t elaborate too much on it since I would like to write about it more in-depth later. Position space deformation, in opposition to skeleton space deformation, is a skinning/rigging/animation technique used a lot in VFX. The goal is to reduce loss of volume in shapes during torsion, like in the elbow or the shoulder.
Artists will sculpt key-poses for shapes. In these poses they manually correct loss of volume, like in the picture below. RBF interpolation is then used to interpolate vertices positions between poses, each vertex being a sample. To interpolate, we use the rigidly deformed vertices as base positions. The interpolation then corrects all the positions according to the sculpted poses. Obviously naive approach to RBF interpolation won’t work since for a shape with 10k vertices you will need to inverse (10k * number of poses) matrices. Weta Digital addressed this issues in a 2010 Siggraph talk for the PSD system on Avatar.
References
I used benchmark . Though it is a bit overkill if you only use it for mat4/vec3, and probably to much template oriented, at least the syntax doesn’t get in the way.
obsolete to say the least. However it is generally well-written and contains a lot of useful info. I would advise the curious to take a look at it one way or another.
Pose space deformation: a unified approach to shape interpolation and skeleton-driven deformation, by Lewis, J.P. and Cordner, M. and Fong, N.