[Insight-users] Voxel Value Ratio Registration Metric
Geoff Topping
g_topping at hotmail.com
Wed Jul 15 05:41:20 EDT 2009
Thanks for the reply Luis,
I have two further questions:
1) In your equation 4, shouldn't there be a factor of 2 from taking the derivative of sum(x(i)) * sum(x(i)) ? I only see the 2 for the derivative of sum(x(i)*x(i)).
2) Why is ImageToImageMetric::GetGradientImage named that, and not GetMovingImageGradientImage? My confusion about whether I needed to divide by the fixed image voxel value was partly due to not being sure what image the gradient was being produced from. It would seem clearer to specify in the function name. (Although now that I've seen your explanation, it's clear enough what is meant and why, it could be clearer...)
Thanks,
Geoff Topping
Date: Tue, 14 Jul 2009 11:40:32 -0400
Subject: Re: [Insight-users] Voxel Value Ratio Registration Metric
From: luis.ibanez at kitware.com
To: g_topping at hotmail.com
CC: insight-users at itk.org
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.
_____________________________________
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
_________________________________________________________________
Internet explorer 8 lets you browse the web faster.
http://go.microsoft.com/?linkid=9655582
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090715/9ee28700/attachment.htm>
More information about the Insight-users
mailing list