Interpolate BSpline transform to region of interest

1

I would like to fine-tune a BSpline transform in a region of interest. Since I already have a coarse bspline for the whole MRI, I would like to interpolate this transform to the region of interest (I don’t really care about what happens outside this region, but the alignment in the region should improve).

I have following helper functions

def interpolate_bspline(
    original_bspline: sitk.BSplineTransform, new_bspline: sitk.BSplineTransform
) -> sitk.BSplineTransform:
    """Interpolate some original BSpline transform to a target transform"""

    dim = new_bspline.GetDimension()
    sp_order = new_bspline.GetOrder()

    img_params = []
    for i in range(dim):
        # resampling the image field to the desired BSpline
        resample = sitk.ResampleImageFilter()
        resample.SetInterpolator(sitk.sitkBSpline)
        resample.SetDefaultPixelValue(0)
        # by default the Identity is used as transform
        resample.SetSize(new_bspline.GetCoefficientImages()[i].GetSize())
        resample.SetOutputSpacing(new_bspline.GetCoefficientImages()[i].GetSpacing())
        resample.SetOutputOrigin(new_bspline.GetCoefficientImages()[i].GetOrigin())
        resample.SetOutputDirection(new_bspline.GetTransformDomainDirection())
        out = resample.Execute(original_bspline.GetCoefficientImages()[i])

        decompose = sitk.BSplineDecompositionImageFilter()
        decompose.SetSplineOrder(sp_order)

        img_params.append(decompose.Execute(out))

    return sitk.BSplineTransform(img_params, sp_order)

def create_bspline(
    fixed_image: sitk.Image, grid_spacing: Sequence[float], order: int = 3
) -> sitk.BSplineTransform:
    image_physical_size = [
        size * spacing
        for size, spacing in zip(fixed_image.GetSize(), fixed_image.GetSpacing())
    ]
    mesh_size = [
        int(image_size / grid_spacing + 0.5)
        for image_size, grid_spacing in zip(image_physical_size, grid_spacing)
    ]
    tx = sitk.BSplineTransform(fixed_image.GetDimension(), order)
    tx.SetTransformDomainDirection(fixed_image.GetDirection())
    tx.SetTransformDomainOrigin(fixed_image.GetOrigin())
    tx.SetTransformDomainPhysicalDimensions(image_physical_size)
    tx.SetTransformDomainMeshSize(mesh_size)
    return tx

now, try to use the interpolation like this:

            tx_roi = create_bspline(fixed_image, grid_spacing=[15.0] * 3)
            tx_roi = interpolate_bspline(tx_existing , tx_roi)

To test if the interpolation works I applied it to the moving image, but the result looks wrong.

I got the idea to use BSplineDecompositionImageFilter from this example: ITK: Examples/RegistrationITKv4/DeformableRegistration15.cxx

What am I doing wrong?

2

Not sure if interpolating the coefficients makes sense.

Meanwhile I found this post. I will give the LandmarkBasedTransformInitializerFilter a try (by sampling my bspline displacements at fixed/moving points).

3

using that approach I get the same issue as reported there

ITK ERROR: BSplineScatteredDataPointSetToImageFilter(000002B8FEDF8530): The reparameterized point component 5.09091 is outside the corresponding parametric domain of [0, 5)

For your reference:

def interpolate_bspline2(
    original_bspline: sitk.BSplineTransform, fixed_image: sitk.Image
) -> sitk.BSplineTransform:
    
    dim = fixed_image.GetDimension()
    num_control_pts = 8

    nx, ny, nz = tuple(fixed_image.GetSize())
    ix = np.linspace(0, nx, num=num_control_pts)
    iy = np.linspace(0, ny, num=num_control_pts)
    iz = np.linspace(0, nz, num=num_control_pts)
    ixv, iyv, izv = np.meshgrid(ix, iy, iz, indexing="ij")
    indices = np.vstack([ixv.ravel(), iyv.ravel(), izv.ravel()]).T

    fixed_pts = []
    moving_pts = []
    for ip in indices:
        fp = fixed_image.TransformContinuousIndexToPhysicalPoint(ip)
        mp = original_bspline.TransformPoint(fp)
        fixed_pts += list(fp)
        moving_pts += list(mp)

    landmark_tx = sitk.LandmarkBasedTransformInitializer(
        transform=sitk.BSplineTransform(dim),
        fixedLandmarks=fixed_pts,
        movingLandmarks=moving_pts,
        referenceImage=fixed_image,
        numberOfControlPoints=num_control_pts)
    return typing.cast(sitk.BSplineTransform, landmark_tx)
4

That error means that the displacement origin for one or more points is outside the defined B-spline domain. I wrote the underlying B-spline filter but not the landmark initializer so I can’t say where the error might be. One thing you can try is the ANTs tool ants.fit_transform_to_paired_points for which several deformable transforms are available and for which a handful are based on the ITK B-spline filter. A small, self-contained example is here.

2 Likes