top of page
    • Nov 1, 2022

RBF Interpolation in Houdini




Recently I set out to understand interpolation more as a concept in the context of computer graphics and vfx work. I don't think it's very important to know everything about interpolation to create great-looking vfx, but for me it can help to understand underlying processes to implement more complicated and interesting techniques in


the future. Additionally, computer graphics and Houdini heavily rely on the concept of interpolation in for most of its operators and processes.

For example, the resample SOP uses the idea of interpolation to generate new point positions that lie between surrounding input points. The key defining factor of interpolation is that it creates new data (new points in this case) based off an arbitrary set of input data (input points/ polyline). Subdivision, smooth, poly-reduce, point deform, divide, scatter – all of these use interpolation to achieve their function.




Here's a list of some different types of interpolation

  1. Linear interpolation (bilinear and trilinear) – creates a linear function that describes the relationship between two points of data, repeated for the amount of dimensions needed. For example in houdini this is lerp(), slerp(), etc vex functions.

  2. Polynomial interpolation (spline) – creates a polynomial function that describes a line that moves through all input data points. Input values can be run through a basis function that can adjust curve profiles (linear, cubic, etc).

  3. ​IDW approximation – inverse distance weighting; this averages surrounding weighted values relative to input data’s proximity


to surrounding points


As a challenge while learning more about interpolation, I set out to create a point deformer using radial basis function (RBF) interpolation. This method interpolates data based solely on the distances between input data by creating weights associated with values.


It is generally more computationally heavy than other methods for large data-sets because it involves solving a system of equations and also creating large matrices. The benefit is that it produces extremely smooth and more visually consistent results and will affect all points differently regardless of distance from the deformed data. In this way it deforms the whole space if used in the context of (xyz) coordinates and moving points. An example of how smooth it is when interpolating a large data set ( volumes and geometry across a large area) below:


Implementation

Let's create the tool from scratch in Houdini. First we store the difference between the variable in an attribute called "@_def". In the case of position, this would represent the difference between the static and deformed mesh.

Next we need to use a python node to with the _def geometry plugged into input 1 ( 2nd input on the python sop). We have to switch to python to build the interpolant ( a matrix that describes the interpolation ) because vex doesn't support more than a 4x4 matrix. I'll break down the individual python components then share the full snippet at the end.


Here we prep the coordinates into arrays:

# create arrays of x, y, z coords

def organize(geo):
    x = []
    y = []
    z = []
    
    for pt in geo.points():
        pos = list(pt.position())
        x.append(pos[0])
        y.append(pos[1])
        z.append(pos[2])
    return (x, y, z)
    
o = organize(data)
x = o[0]
y = o[1]
z = o[2]

Then similarly we need to assemble the deformation values into arrays 

xd = []
yd = []
zd = []

for point in data.points():
    xd.append(point.attribValue('_def')[0])
    yd.append(point.attribValue('_def')[1])
    zd.append(point.attribValue('_def')[2])

Then I build a function that takes a value (r) and runs that through the selected basis function (this is where basis comes from in "radial basis function"):

def phi(r, constant=10):
    c = node.evalParm('constant')
    if node.evalParm('basis_function') == 0:
        # linear
        return r
    if node.evalParm('basis_function') == 1:
        # cubic
        return r*r*r
    if node.evalParm('basis_function') == 2:
        # thin plate
        return r * r * np.log(r)
    if node.evalParm('basis_function') == 3:
        # multi quad
        return np.sqrt(pow(r*c, 2)+1.0)
    if node.evalParm('basis_function') == 4:
        # gauss
        return np.exp(c*pow(r, 2))
    if node.evalParm('basis_function') > 4:
        # linear if off list 
        return r

This post is not finished - 11/4/2022

bottom of page