Parametric Surface Creation in Python

Question

Is there a Python module for handling parametric (u-v) surfaces? I'm looking for something that's the 3D analogue to scipy.interpolate's spline functions, where I can create parametric splines through a set of 2D points thusly:

xypts = [[0., 1., 5., 2.], [4., 3., 6., 7.]]
tck, u = scipy.interpolate.splprep(xypts, s=0, k=3)

and then get a point at any t-value on the spline like this:

t = 0.5
intxypt = scipy.interpolate.splev(t, tck)

So, what I'd like is something that works like this:

# xyzpts is a 3 x m x n matrix, with m and n >= 4 for a cubic surface
tck, s, t = srfprep(xyzpts, s=0, k=3)
u, v = 0.5, 0.5
intxyzpt = srfev(u, v, tck)

I wrote my own code to do just this some time ago, but frankly it kinda sucks (slow and fragile, especially at the surface edges) and I'm looking for something more standard and optimized.

Answer 1

This is probably obvious, but if you are able to guess the u-v coordinates corresponding to each data point (simplest case u=x, v=y if the surface is a graph), parametric interpolation (u,v) -> (x,y,z) is essentially 2-D interpolation of 3 separate datasets (the x, y, and z coordinates), so you can use any usual 2-D interpolation method.

splprep in fact works this way, by assuming the points are ordered and assigning the u coordinates according to u[i] = u[i-1] + dist(p[i], p[j]) using the euclidean distance. This generalizes to 2-D if you know which points are "next to each other". For example, if the x,y,z data come as 2-D arrays, you can do

from scipy import interpolate
import numpy as np

# example dataset (wavy cylinder)

def surf(u, v):
    x = np.cos(v*np.pi*2) * (1 + 0.3*np.cos(30*u))
    y = np.sin(v*np.pi*2) * (1 + 0.3*np.cos(30*u))
    z = 2*u
    return x, y, z

ux, vx = np.meshgrid(np.linspace(0, 1, 20),
                     np.linspace(0, 1, 20))
x, y, z = surf(ux, vx)

# reconstruct (u, v) using the existing (!) neighbourhood information
du = np.sqrt(np.diff(x, axis=0)**2 + np.diff(y, axis=0)**2 + np.diff(z, axis=0)**2)
dv = np.sqrt(np.diff(x, axis=1)**2 + np.diff(y, axis=1)**2 + np.diff(z, axis=1)**2)
u = np.zeros_like(x)
v = np.zeros_like(x)
u[1:,:] = np.cumsum(du, axis=0)
v[:,1:] = np.cumsum(dv, axis=1)

u /= u.max(axis=0)[None,:] # hmm..., or maybe skip this scaling step -- may distort the result
v /= v.max(axis=1)[:,None]

# construct interpolant (unstructured grid)
ip_surf = interpolate.CloughTocher2DInterpolator(
        (u.ravel(), v.ravel()), 
        np.c_[x.ravel(), y.ravel(), z.ravel()])

# the BivariateSpline classes might also work here, but the above is more robust

# plot projections
import matplotlib.pyplot as plt

u = np.random.rand(2000)
v = np.random.rand(2000)

plt.subplot(131)
plt.plot(ip_surf(u, v)[:,0], ip_surf(u, v)[:,1], '.')
plt.title('xy')
plt.subplot(132)
plt.plot(ip_surf(u, v)[:,1], ip_surf(u, v)[:,2], '.')
plt.title('yz')
plt.subplot(133)
plt.plot(ip_surf(u, v)[:,2], ip_surf(u, v)[:,0], '.')
plt.title('zx')
plt.show()

EDIT: Ok, I'm not fully sure how robust the above computed u,v in practice are, as it seems there's room for distortion. However, the LocallyLinearEmbedding below may work better in this respect.

If you are not able to guess the u,v values, for instance you have just a bunch of points and no neighborhood information, the problem becomes more difficult. The appropriate keywords here seem to be "surface reconstruction" and "manifold learning".

I didn't try, but it seems to me that you easily can get suitable u,v coordinates using LocallyLinearEmbedding from scikits-learn, see this example. They have a bunch of different algorithms, and this seems solid enough. The resulting u = Y[:,0]; v = Y[:,1] you can then use in unstructured 2-D interpolation methods, as shown above.

Maybe googling more will reveal more packages.