[Insight-users] Voxel Value Ratio Registration Metric

[Geoff Topping]

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=9666047
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090709/e5f4546d/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: RatioVarianceImageToImageMetricExcerpts.txx
Type: application/octet-stream
Size: 11158 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090709/e5f4546d/attachment.obj>