Colleagues of mine use a spline library called splinepy. It has some nice features, I wanted to try out, using some ITK BSplineTransforms I have obtained from image registration. However, I wasn’t able to convert them into the representation which splinepy accepts and also my colleagues struggled to understand how a BSplineTransform is actually created in ITK.
I think the fundamental difference is, that in splinepy, I’m required to give a knot vector and the positions of the control points, while in ITK I just give displacements of control points, which are on a regular grid. Finding the CP coordinates and converting displacements to positions is not a big deal, however, we were not able to produce the correct knot vector, which would give the same spline. I tested my conversion by transforming points in ITK and by requesting the point in splinepy. However, I can only make it work when using a spline of degree 1.
I tried to lookup the original literature, on which the ITK implementation is based on.
What I could find is, that the ITK implementation is apparently based upon the work of Unser et al. - however, I wasn’t able to deduce from skimming the paper how the spline is actually implemented and what the knot vector might be, but I may just have missed important information in this paper (perhaps because at first glance they used yet another notation than in The NURBS Book or in the paper by Lee et al)… Reading in the ITK source code wasn’t a success either, but again, I may just have not found the important bits.
I believe that the issue is in the way the spline is setup in combination with the transformation domain. I know that ITK will place control points outside the transformation domain, depending on the order (btw, I believe that “order” in ITK actually means “degree”?! A p-degree spline has order p+1, but I think that ITK uses the term order instead of degree?)
If I remember correctly, if the spline is written down in the way as Lee et al does it (see Eq.1), then the spline is not interpolating until the edge of the control point grid - hence it should not have knot multiplicity p+1 at the ends of the knot vector. But messing around with different knot vectors has not worked…
To show what is already working, here is an example:
import SimpleITK as sitk
import splinepy as sp
import numpy as np
u = np.array([
0.0, 0.0,
0.0, 0.0,
0.0, 0.0,
0.0, 0.0,
0.1, 0.2,
0.0, 0.0,
0.0, 0.0,
0.0, 0.0,
-0.1, 0.3,
]).reshape(9, 2)
# splinepy:
s = sp.BSpline(degrees=[1, 1], knot_vectors=[[0, 0, 0.5, 1, 1], [0, 0, 0.5, 1, 1]], control_points=sp.utils.data.cartesian_product(np.array([[0, 0.5, 1], [0, 0.5, 1]])) + u)
# SimpleITK
bs = sitk.BSplineTransform(2, 1) # Create a new BSpline, dimension=2, order=1
# 3x3 grid, origin=0,0, spacing=0.5,0.5, direction=unity
bs.SetFixedParameters((3.0, 3.0, 0.0, 0.0, 0.5, 0.5, 1.0, 0.0, 0.0, 1.0))
bs.SetParameters(u.flatten(order='F'))
test_points = [[0, 0], [0.5, 0.5], [0.3, 0.75], [1, 0], [0.25, 0.33]]
for point in test_points:
print(f"p={str(point) :15s} splinepy: {str(s.evaluate([point])[0].tolist()):30s} sitk: {str(list(bs.TransformPoint(point))):30s}")
with the output:
p=[0, 0] splinepy: [0.0, 0.0] sitk: [0.0, 0.0]
p=[0.5, 0.5] splinepy: [0.6, 0.7] sitk: [0.6, 0.7]
p=[0.3, 0.75] splinepy: [0.32999999999999996, 0.81] sitk: [0.32999999999999996, 0.81]
p=[1, 0] splinepy: [1.0, 0.0] sitk: [1.0, 0.0]
p=[0.25, 0.33] splinepy: [0.28300000000000003, 0.396] sitk: [0.28300000000000003, 0.396]
which looks really good! But here is with one degree higher:
# Assuming I need tripple knot multiplicity at the ends...
# Assuming I need to have the parametric coordinate running from -0.5 to 1.5, as the
# ITK BSpline will have this as the control point grid.
s = sp.BSpline(degrees=[2, 2], knot_vectors=[[-0.5, -0.5, -0.5, 1.5, 1.5, 1.5], [-0.5, -0.5, -0.5, 1.5, 1.5, 1.5]], control_points=sp.utils.data.cartesian_product(np.array([[-0.5, 0.5, 1.5], [-0.5, 0.5, 1.5]])) + u)
bs = sitk.BSplineTransform(2, 2)
bs.SetParameters(u.flatten(order='F'))
test_points = [[0, 0], [0.5, 0.5], [0.3, 0.75], [1, 0], [0.25, 0.33]]
for point in test_points:
print(f"p={str(point) :15s} splinepy: {str(s.evaluate([point])[0].tolist()):30s} sitk: {str(list(bs.TransformPoint(point))):30s}")
Here, it fails:
p=[0, 0] splinepy: [0.013671874999999991, 0.029296874999999997] sitk: [0.025, 0.05]
p=[0.5, 0.5] splinepy: [0.51875, 0.56875] sitk: [0.5546875, 0.6171875]
p=[0.3, 0.75] splinepy: [0.31625000000000003, 0.81375] sitk: [0.347546875, 0.851421875]
p=[1, 0] splinepy: [1.010546875, 0.03867187500000001] sitk: [1.0250000000000001, 0.050000000000000135]
p=[0.25, 0.33] splinepy: [0.27033824218749997, 0.38278605468750004] sitk: [0.29940546875, 0.42966171875000003]
Maybe someone with more BSpline knowledge can give me a hint how to do that?
).
