[Insight-users] Voxel Value Ratio Registration Metric

[Luis Ibanez]

Hi Geof,

Thanks for your detailed email.

Here is how you may want to compute the derivatives of this metric.

Your metric is:

M = (1 / (n - 1)) * ( sum(x(i)*x(i)) - (1/n) * sum(x(i)) * sum(x(i)) )
[Eq.1]

where

x(i) = m(i) / f(i)      [Eq.2]

f(i) = i-th pixel of the fixed image
m(i) = moving image pixel corresponding to the f(i) pixel under the
          current Transform parameters.

strictly speaking m(i) should be written as

m( q(i),T)   where

p(i) = physical coordinate of pixel f(i)   and

q(i) = T ( p(i) ) = point in the moving image to
         which the Transform has mapped point
         p(i).

T = array of Transform parameters

That is, p(i) is a point in the  coordinate system of the
fixed image, while q(i) is its corresponding point in the
moving image, under the Transformation T.


The Metric derivative is computed against the array of
Transform parameters. Therefore we are looking for

  dM / dT

and therefore it follows

dM/dT = d (1 / (n - 1)) * ( sum(x(i)*x(i)) - (1/n) * sum(x(i)) * sum(x(i)) )
/ dT    [Eq.3]

which is

(1 / (n - 1)) * ( 2 * sum( x(i) * dx(i)/dT ) - (1/n) * sum(x(i)) * sum( d
x(i) /d T ) )   [Eq.4]

the terms d x(i) / d T  are

d x(i) / dT = d ( m(i)/f(i) ) / d T     [Eq. 5]

f(i) is independent of T, therefore

d x(i) / dT  = ( 1 / f(i) ) * d ( m(i) ) / dT          [Eq.6]

where

d ( m(i) ) / dT  is "qualitatively": how much the intensity of the
corresponding
pixel m(i) will change, with respect to each one of the components of the
Transform parameter array T.

This term is the one that you can compute (using the chain-rule) as a
product
of the moving image gradient and the Transform Jacobian.

Strictly speaking:

m(i)  is   m( q(i),T)   (see above)

d (m(i) ) / dT  = d ( m( q(i), T ) / dT =  GM( q(i) )  *  J( p(i), T)
[Eq.7]

Where

GM( q(i) ) is the Gradient of the moving image evaluated at point q(i)

where, again, q(i) is the point in the moving image that corresponds
to point p(i) of the fixed image.

J( p(i), T ) is the Jacobian matrix of the Transform that indicates
how much the coordinates of the image point q(i) will change, with
respect to variations of each of the components of the Transform
parameters array.

Once you compute d( m(i) ) / dT, you can put this value back in
equation 6., to compute the term  d x(i) / dT

 d x(i) / dT = ( 1 / f(i) ) * d( m(i) ) / dT

and with this value, you can compute Equation 4. which will give you
the value of     d M / d T.


To quickly answer your question:

YES, you have to multiply   d m(i) / dT     by    ( 1 / f(i) ),
in order to get d x(i)  / dT.


If you are looking for a detailed description of the calculus used
for computing Metric derivatives, you will find a discussion on this
computation in the book:


"Insight into Images"
http://www.amazon.com/Insight-into-Images-Segmentation-Registration/dp/1568812175



BTW: As usual,
Once you get this metric to work, we strongly encourage you
to share it with the ITK community by submitting it as a paper
to the Insight Journal   :-)



     Regards,


           Luis



-------------------------------------------------------------------------------------------------
On Thu, Jul 9, 2009 at 9:02 AM, Geoff Topping <g_topping at hotmail.com> wrote:

>  Hi all,
>
> I've been working on making a new ITK image registration metric based on
> the variance ratio of voxel values, similar to the automated image
> registration (AIR) program by Roger Woods.  I've got some initial code and
> theory I'd like to run by anyone who'll have a look.
>
> The metric measure is given by the variance of the included voxels' value
> ratio between the moving and fixed image.  If f(i) is the fixed image voxel
> value with index i, and m(i) is the corresponding moving image voxel value,
> and x(i) is their ratio,
>
> x(i) = m(i) / f(i)
>
> Then the measure M is given by the variance of these values, or, for all
> accepted voxels (which must have nonzero f(i) and be in the accepted
> region), with n voxels included:
>
> M = (1 / (n - 1)) * ( sum(x(i)*x(i)) - (1/n) * sum(x(i)) * sum(x(i)) )
>
> This form of the variance calculation can be derived from the formula for
> sample variance formula which is (1 / (n - 1)) * sum( (x(i) - mean(x))^2 )
>
> where mean(x) is estimated from the sample data.
>
> M is relatively easy to implement as a measure, as the sums of x(i)^2 and
> x(i) can be calculated with a single loop over the image data.
>
> The derivative is a bit more complicated to calculate, and is where I'd
> like some additional attention from anyone experienced in writing this sort
> of function for an ITK metric.
>
> For each transform parameter s, I took the derivative of the metric M to
> get:
>
> dM/ds = (2 / (n - 1)) * ( sum( x(i) * dx(i)/ds ) - (1/n) * sum(x(i)) *
> sum(dx(i)/ds) )
>
> where dx(i)/ds is somewhat compliated to determine.
>
> Looking at MeanSquaresImageToImageMetric and
> NormalizedCorrelationImageToImageMetric, the dx(i)/ds term seems to be
> determined using a jacobian matrix and a gradient image.  I'm not sure
> exactly how the product of (the jacobian for the appropriate transform
> parameter and fixed image dimension) and (the gradient image for the
> appropriate dimension at the appropriate pixel) should be used to get my
> dx(i)/ds.  I couldn't find any theoretical explanation of this code, so am
> working backwards to try to understand how the examples I have work.
>
> I'm currently just using the product of the gradient and jacobian to give
> me the dx(i)/ds value, but I'm wondering if this needs to be divided by the
> fixed image voxel value as well?  I'm not sure whether the jacobian and
> gradient product is giving me the differential of the moving image value, or
> the differential of the ratio of the moving and fixed image values.  Since
> x(i) = m(i) / f(i), this uncertainty would determine whether I need an extra
> division.
>
> If someone could clarify that, or whether I'm understanding this at all, or
> perhaps that my derivation is right, that'd be nice.
>
> I've also attached my current GetValue and GetDerivative functions code.
> My tests suggest they're at least somewhat functional, although in some
> cases the measure will improve for a few registration iterations, and then
> get progressively worse for the remainder, suggesting something might be
> wrong with the derivative.  If somsone could look that over to see whether
> I'm approaching this correctly, that'd be nice.  This isn't meant to be a
> completed insight journal submission yet, though.
>
> Thanks,
> Geoff Topping
> ------------------------------
> We are your photos. Share us now with Windows Live Photos.<http://go.microsoft.com/?linkid=9666045>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090714/3ace14be/attachment.htm>