Python-VTK 3D Spline Regression through STL model of vascular tree

I need to create a spline or polyline representation of a vascular tree model (see below).

The model is in a STL format, thus I have the x-y-z coordinates of all vertices. The lines should run through the center of the vessel mesh thus I thought that the best approach would be a spline regression through the vertex cloud. In addition it would be great if I can have the radius of the vessel at given points, e.g. the coordinates of the polyline.

I looked through this forum and the VTK website (assuming they have a straightforward implementation for this sort of thing) but so far I haven't found something I can use. Does anyone know of a Python module or VTK class (which I would call from Python) that can do this? The python modules I found on this are all for 2D data.

Thanks a lot!

EDIT: I came across this library called VMTK that deals almost exclusively with vessel segmentation and has functionality for what they call 'centerline calculation'. However, they usually require the vessels to be 'cut' at their ends and 'source points' to be defined. In the case of my model, however, one can see that the end points are 'capped' which makes matters more complicated. If I find a solution I'll post here

I don't know any software or python classes exactly on your problem. Maybe python interpolate.splev will help you with a single vessel. You may try the following code as an example:

from scipy import interpolate
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np

# 3D example
total_rad = 10
z_factor = 3
noise = 0.1

num_true_pts = 200
s_true = np.linspace(0, total_rad, num_true_pts)
x_true = np.cos(s_true)
y_true = np.sin(s_true)
z_true = s_true/z_factor

num_sample_pts = 100
s_sample = np.linspace(0, total_rad, num_sample_pts)
x_sample = np.cos(s_sample) + noise * np.random.randn(num_sample_pts)
y_sample = np.sin(s_sample) + noise * np.random.randn(num_sample_pts)
z_sample = s_sample/z_factor + noise * np.random.randn(num_sample_pts)

tck, u = interpolate.splprep([x_sample,y_sample,z_sample], s=2)
x_knots, y_knots, z_knots = interpolate.splev(tck[0], tck)
u_fine = np.linspace(0,1,num_true_pts)
x_fine, y_fine, z_fine = interpolate.splev(u_fine, tck)

fig2 = plt.figure(2)
ax3d = fig2.add_subplot(111, projection='3d')
# blue line shows true helix
ax3d.plot(x_true, y_true, z_true, 'b')
# red stars show distorted sample around a blue line
ax3d.plot(x_sample, y_sample, z_sample, 'r*')
# green line and dots show fitted curve
ax3d.plot(x_knots, y_knots, z_knots, 'go')
ax3d.plot(x_fine, y_fine, z_fine, 'g')
plt.show()

This code uses noisy centerline path of a single vessel and fit it with a smooth curve (see the result below):

interolation result

Usually, two user seeds are used to mark centerline ends, in the case of centerline representation as in VMTK.

The other way to get centerlines automatically is to voxelize your stl mesh, costruct a voxel skeleton, and separate skeletal segment to represent each vessel. Then you can interpolate each centerline to get the smooth curves. Unprocessed skeletal segments usualy have zigzags.