Get Jacobian wrt to Position for BSplineTransform

As far as I can tell, the function ComputeJacobianWithRespectToPosition is not implemented for BSplineTransform (as it is for DisplacementFieldTransform) but only ComputeJacobianWithRespectToParameters, for example using this snippet:

using BSplineTransformType = itk::BSplineTransform<double, 3, 3>;
using JType = BSplineTransformType::JacobianType;
using PType = BSplineTransformType::InputPointType;

//  [...] read transform into tfm from file
const BSplineTransformType::Pointer bspline = static_cast<BSplineTransformType *>(tfm.GetPointer());

auto J = JType();
// some point p
auto p = PType({4, 5, 6});
// get Jacobian at point {4, 5, 6}
bspline->ComputeJacobianWithRespectToParameters(p, J);

Is there a way to get the Jacobian matrix from the Jacobian of the Parameters or is the easiest way to use a DisplacementFieldTransform for that? I guess I can apply the chain rule to get the correct derivatives, but how to get to the missing derivates… I have no idea :sweat_smile:
I have searched for a while for an example or some usage in the code but could not find any - or I’m overlooking something…

One option is to convert the bspline transform to displacement field transform and use its ComputeJacobianWithRespectToParameters method.

Yep, I did that so far. It requires me to build the displacement field first, which takes some amount of memory and time. Calculating it directly from the parametrization would be much more efficient.

Maybe it isn’t that easy and thus the ComputeJacobianWithRespectToPosition isn’t implemented?

If you want to compute the derivatives at random points then you can use vtkGridTransform. It is pretty amazing, it can compute derivatives of both forward and inversre transform at any point, very quickly. We use VTK transform classes in 3D Slicer exactly because ITK requires computing a dense displacement field for many operations, which is just way too slow.

That one sounds extremely interesting! However, why does it not compute the derivative of the inverse:

Calculate the inverse transform as well as the derivative of the forward transform (that’s correct: the derivative of the forward transform, not of the inverse transform)

VTK: vtkGridTransform Class Reference

This is actually a good thing! It means when you implement a new transform type then you only need to implement the derivative of the forward transform. The inverse transform and its derivative are automatically computed for you! To get the inverse transform, simply call the GetInverse() method and then you can get the derivative from that.

It works for all important transform types (linear, b-spline, thin-plate-spline, and grid), even for arbitrary composition and any combination of inversion of any components. Transformation is available for the entire 3D space (not limited to the area of b-spine grid, etc) by smooth extrapolation, allowing robust inverse computation near the edges.

This is all really amazing, pretty much solves all spatial transformation related tasks in medical image computing. Credits go to David Gobbi.

1 Like

Could you hint me to an example where this is used in combination with BSpline to get the derivatives of the forward and inverse transform?
Is it possible to translate an existing ITK BSplineTransform to VTK?