Jacobian of B-spline coefficients

For a motion-estimation algorithm, we’d like to be able to calculate the Jacboian of b-spline coefficients. Alternatively (and less computationally expensive), we’d like to be able to apply this Jacboian onto a pre-existing transformation.

We would like to do this in VTK or ITK, whichever is easier. Matching ITK post is here.

Presumably classes like vtkImageBSplineCoefficients and vtkBSplineTransform might do what we want, but I couldn’t find methods with names I was expecting when looking through their documentation.

Could anyone point us to the relevant classes/methods (dare I say examples) for this?

The vtkBSplineTransform::InternalTransformDerivative() method does this, it is declared in the vtkAbstractTransform base class:

  /**
   * This will transform a point and, at the same time, calculate a
   * 3x3 Jacobian matrix that provides the partial derivatives of the
   * transformation at that point.  This method does not call Update.
   * Meant for use only within other VTK classes.
   */
  virtual void InternalTransformDerivative(
    const float in[3], float out[3], float derivative[3][3]) = 0;
  virtual void InternalTransformDerivative(
    const double in[3], double out[3], double derivative[3][3]) = 0;

This method computes the partial derivatives of the output coords with respect to the input coords. It’s called “Internal” just because its primary purpose is to be used internally by other vtkAbstractTransform methods, but it’s a public method so it’s fine to call it directly. Of course this method only computes the Jacobian at a single point.

Most VTK users do not grab the Jacobian directly. A more typical usage pattern, e.g. if someone has a finite element mesh that includes vector data, is to use vtkTransformFilter to apply a transform to the mesh. In this case, VTK takes care of the details of applying the Jacobian to the vectors.

As a follow-up, if you are specifically interested in concatenating grid transforms, this can be achieved with the vtkGeneralTransform class. This class has a similar interface to the vtkTransform class:

auto multi_transform = vtkSmartPointer<vtkGeneralTransform>::New();
multi_transform->SetInput(transform_A);
// PostMultiply() adjusts the order of concatenation for the Concatenate
// method, so that "this = transform_B*transform_A",
// i.e. B is applied after A, rather than vice-versa.
multi_transform->PostMultiply();
multi_transform->Concatenate(transform_B);
multi_transform->Update();

The vtkGeneralTransform has vtkAbstractTransform as its base, so it can be used in VTK wherever a vtkAbstractTransform can be used (e.g. it can be used by vtkTransformToGrid). It automatically applies the Jacobian as needed while performing transformations. To generate a new grid, use this as the input of vtkTransformToGrid.

This is great info, thanks. We’ll need to think a little bit more about the best way to do it, but it’s good to know the functionality is there (don’t think we’ll be concatenating transforms, but useful info too).

3D Slicer uses VTK transform’s dynamic inverse computation when an inverse is needed. We only compute an inverse displacement field when you choose to export the segmentation into a displacement field.

Of course an inverse does not always exist (e.g., when the forward transform maps multiple input positions into the same output position). You can detect validity of an inverse with quite good confidence by transforming a point with the inverse transform and then applying the forward transform and check if you get the original position.

In Slicer, we implemented variants of VTK transforms that compute the Jacobian for arbitrarily oriented volumes (in vtkAddon library). Now that VTK supports oriented image data natively, these computations could be moved into VTK core.