Image registration metric derivative

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!

The following resources may be useful:

But, the source code is the ultimate authority.

1 Like

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:

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


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

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.

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.

/* 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.

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

Thanks, this is really helpful.