happysplines.py

import numpy as np, matplotlib.pyplot as pp, mpl_toolkits.mplot3d as mp
import scipy.interpolate as si, scipy.optimize as so, scipy.spatial.distance as ssd, scipy.integrate

data = (1,2,3,4,5,6,7,8),(1,2.5,4,5.5,6.5,8.5,10,12.5),(1,2.5,4.5,6.5,8,8.5,10,12.5)
p = 6.5,9.5,9

# Fit a spline to the data - s is the amount of smoothing, tck is the parameters of the resulting spline
(tck, uu) = si.splprep(data, s=0)

# Return distance from 3d point p to a point on the spline at spline parameter u
def distToP(u):
	s = si.splev(u, tck)
	return ssd.euclidean(p, s)

# Return the magnitude of the gradient of the spline at u
def slopeMag(u):
	return ssd.euclidean(0, si.splev(u, tck, 1))

# Return the distance from u to v along the spline
def distAlong(u, v):
	return scipy.integrate.quad(slopeMag, u, v)[0]

# Return the distance from u along the spline to the halfway point
def distAlongToHalf(u):
	return abs(distAlong(0.0, u) - (distAlong(0.0, 1.0)/2))

# Find the closest point on the spline to our 3d point p
# We do this by finding a value for the spline parameter u which
# gives the minimum distance in 3d to p
closestu = so.fmin(distToP, 0.5)
closest = si.splev(closestu, tck)

# Find out where on the spline the halfway point is
# Spline parameter u = 0.5 isn't necessarily halfway along the curve
halfwayu = so.fmin(distAlongToHalf, 0.5)
halfway = si.splev(halfwayu, tck)
s = "distance along spline to halfway point is %f" % distAlong(halfwayu, closestu)

# Build a list of 3d points along the spline
spline = si.splev(np.linspace(0,1,100), tck)

# Plot things!
ax = mp.Axes3D(pp.figure(figsize=(8,8)))
ax.plot((p[0],), (p[1],), (p[2],), 'o', label='input point')
ax.plot(closest[0], closest[1], closest[2], 'o', label='closest point on spline')
ax.plot(halfway[0], halfway[1], halfway[2], 'o', label='halfway along spline')
ax.plot(data[0], data[1], data[2], label='raw data points')
ax.plot(spline[0], spline[1], spline[2], label='spline fit to data')
pp.legend(loc='lower right')
pp.title(s)
pp.show()