Radial Basis Function Interpolation

Python Solve
Lets start by dropping down a Python SOP and initializing our variables in typical fashion. The array lpos will hold the positions of all points in order and can be empty for now. Instances for matrices A and v are created at the proper sizes, assuming a float3 attribute type for v.

import numpy as np
import math

node = hou.pwd()
mode = node.parm("mode").eval()
geo = node.geometry(0)
pts = geo.points()
npts = len(pts)
lpos = []
A = np.empty([npts, npts], dtype=float)
v = np.empty([npts, 3], dtype=float)

Next we’ll iterate over all the source points in a cocatenated loop to fill the rows and columns of matrices A and v and store the source point positions in lpos.

for p in pts:
pP = p.position();
lpos.extend([pP.x(), pP.y(), pP.z()])
v[p.number()] = p.attribValue("def")
for l in pts:
    A[p.number(),l.number()]= basis(mag(l.position() - pP), mode)

You probaly spotted the two functions in there we still haven’t declared yet, lets talk about them. The basis function defines how the point distances stored in A are weighted. We’ll just implement linear d, cubic d^3 and thin plate spline d^2\ln(d) but theres quite a few more to be found on wikipedia (also some utilizing a shaping parameter). Variable m is linked to a parameter (mode) to allow switching between interpolation weighting modes.

def basis(d, m):
    if m==1:
        return(d*d*d)
    elif m==2:
         return(d*d*math.log(d))
    else:
        return(d)

To calculate the length of the distance vector we can define a function as follows. NumPy comes with a builtin function to handle this (np.linalg.norm) but oddly its quite a bit slower than the user-defined variant in this case.

def mag(x):
    return(math.sqrt(sum(i*i for i in x)))

Now that we actually know A and v we can finally solve for x.

x = np.linalg.solve(A, v)

To pass the data collected here on to the deformation stage we’ll return them as detail attributes. The array x is a cocatenated list of lists [[x_1,y_1,z_1],[x_2,y_2,z_2],\dots,[x_n,y_n,z_n]] so we’ll have to flatten it [x_1,y_1,z_1,\dots,x_n,y_n,z_n] to be recognized as an array attribute when binding. Since we will use the OpenCL SOP later, which doesnt have a secondary input, let’s also export the list of source point positions. And the number of source points to avoid checking the array’s length in the OpenCL kernel for every deformation point.

coef = []
for i in x.tolist():
    coef.extend(i)

geo.addAttrib(hou.attribType.Global, "lpts", len(pts))
geo.addArrayAttrib(hou.attribType.Global, "lpos", hou.attribData.Float, 3)
geo.addArrayAttrib(hou.attribType.Global, "coef", hou.attribData.Float, 3)
geo.setGlobalAttribValue("lpos", lpos)
geo.setGlobalAttribValue("coef", coef)

With those sweet detail attributes set-up we can finally do some deformation.