Image registration metric derivative

ibr · August 14, 2018, 3:32pm

I am using image registration for some time now but I am still trying to understand things in more deep. It seems tracing ITK code is not that easy.

What is the math behind the metric derivative and where can I find the implementation part i.e the update function of the transform parameters? it would be nice if someone provide explanation or some simple-to-read references.

Assuming a 2D rigid transformation with only 3 parameters translation in x and y, and a rotation angle; assuming also a mean squared error metric, how is the metric derivative and the update equation look like?

is extend the concept to 3D is straight forward or involving more math and code?

I hope the questions are clear!
Thanks!

matt.mccormick · August 15, 2018, 2:21pm

The following resources may be useful:

The ITKv4 Image Registration Framework article
The Registration Chapter of the ITK Software Guide
John Galeotti’s Course’s ITK Registration Chapter
Materials from the many university courses based on ITK

But, the source code is the ultimate authority.

ibr · August 15, 2018, 3:25pm

Thanks Matt for your concern. What I am looking is more details about the metric derivative. I already went through the links you provided but couldn’t find a clear answer.
If you take a look at the paper in the first link, details about equation 2 is what I need (by the way, I found the same equation in this paper as well) .

I like to know:

How this derivative computed in more details.
How to implement this equation e.g. pseudo code or the part in itk that does it.

Thanks!

matt.mccormick · August 15, 2018, 3:57pm

The code that implements Equation 2 for the first link for the mean squares metric is here:

github.com

InsightSoftwareConsortium/ITK/blob/6b41fe1cbf7a71a6003d4824e02a69b41a508cda/Modules/Registration/Metricsv4/include/itkMeanSquaresImageToImageMetricv4GetValueAndDerivativeThreader.hxx#L42-L83


/** Only the voxelwise contribution given the point pairs. */
FixedImagePixelType diff = fixedImageValue - movingImageValue;
const unsigned int nComponents = NumericTraits<FixedImagePixelType>::GetLength( diff );
metricValueReturn = NumericTraits<MeasureType>::ZeroValue();


for ( unsigned int nc = 0; nc < nComponents; nc++ )
  {
  MeasureType diffC = DefaultConvertPixelTraits<FixedImagePixelType>::GetNthComponent(nc, diff);
  metricValueReturn += diffC*diffC;
  }


if( ! this->GetComputeDerivative() )
  {
  return true;
  }


/* Use a pre-allocated jacobian object for efficiency */
using JacobianReferenceType = typename TImageToImageMetric::JacobianType &;
JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian;
JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobianPositional;

This file has been truncated. show original

It is a bit complicated because common computation is spread across the class hierarchy. To inspect the hierarchy, see the class’s Doxyen page.

We could make these classes easier to follow and more performant if we replaced the DomainThreader classes with the ITK 5 multi-threading API, described in the Multi-Threading section of the ITK Software Guide.

ibr · August 15, 2018, 4:04pm

Thanks a lot, I will give it a try.

It is a bit complicated

Indeed , probably a dedicated book/course on ITK registration would be a good idea to explain the theory and the implementation.

ibr · October 11, 2018, 1:00pm

/* Use a pre-allocated jacobian object for efficiency */ using JacobianReferenceType = typename TImageToImageMetric::JacobianType &; JacobianReferenceType jacobian = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobian; JacobianReferenceType jacobianPositional = this->m_GetValueAndDerivativePerThreadVariables[threadId].MovingTransformJacobianPositional;

I kept tracing it but I could not find where the computation of Jacobian implemented. Could you please add a link to it.

dzenanz · October 11, 2018, 2:41pm

Each transform type can provide its own way to compute a Jacobian. The base implementation for “4x4 transform matrix” is here:

github.com

InsightSoftwareConsortium/ITK/blob/5142ed7bb5e4e5042734fc47bf556a2342be6cec/Modules/Core/Transform/include/itkMatrixOffsetTransformBase.hxx#L647-L673


::ComputeJacobianWithRespectToParameters(const InputPointType & p, JacobianType & jacobian) const
{
// This will not reallocate memory if the dimensions are equal
// to the matrix's current dimensions.
jacobian.SetSize( NOutputDimensions, this->GetNumberOfLocalParameters() );
jacobian.Fill(0.0);


// The Jacobian of the affine transform is composed of
// subblocks of diagonal matrices, each one of them having
// a constant value in the diagonal.
const InputVectorType v = p - this->GetCenter();


unsigned int blockOffset = 0;
for( unsigned int block = 0; block < NInputDimensions; block++ )
  {
  for( unsigned int dim = 0; dim < NOutputDimensions; dim++ )
    {
    jacobian(block, blockOffset + dim) = v[dim];
    }

This file has been truncated. show original

ibr · October 11, 2018, 2:49pm

Thanks, this is really helpful.