Hi Geof,<br><br>Thanks for your detailed email.<br><br>Here is how you may want to compute the derivatives of this metric.<br><br>Your metric is:<br><br>M = (1 / (n - 1)) * ( sum(x(i)*x(i)) - (1/n) * sum(x(i)) * sum(x(i)) ) [Eq.1]<br>
<br>where <br><br>x(i) = m(i) / f(i) [Eq.2]<br><br>f(i) = i-th pixel of the fixed image<br>m(i) = moving image pixel corresponding to the f(i) pixel under the<br> current Transform parameters.<br><br>strictly speaking m(i) should be written as<br>
<br>m( q(i),T) where<br><br>p(i) = physical coordinate of pixel f(i) and<br><br>q(i) = T ( p(i) ) = point in the moving image to <br> which the Transform has mapped point<br> p(i).<br><br>T = array of Transform parameters<br>
<br>That is, p(i) is a point in the coordinate system of the<br>fixed image, while q(i) is its corresponding point in the<br>moving image, under the Transformation T.<br><br><br>The Metric derivative is computed against the array of <br>
Transform parameters. Therefore we are looking for <br><br> dM / dT<br><br>and therefore it follows<br><br>dM/dT = d (1 / (n - 1)) * ( sum(x(i)*x(i)) - (1/n) * sum(x(i)) * sum(x(i)) ) / dT [Eq.3]<br><br>which is<br><br>
(1 / (n - 1)) * ( 2 * sum( x(i) * dx(i)/dT ) - (1/n) * sum(x(i)) * sum( d x(i) /d T ) ) [Eq.4]<br><br>the terms d x(i) / d T are<br><br>d x(i) / dT = d ( m(i)/f(i) ) / d T [Eq. 5]<br><br>f(i) is independent of T, therefore<br>
<br>d x(i) / dT = ( 1 / f(i) ) * d ( m(i) ) / dT [Eq.6]<br><br>where<br><br>d ( m(i) ) / dT is "qualitatively": how much the intensity of the corresponding<br>pixel m(i) will change, with respect to each one of the components of the<br>
Transform parameter array T.<br><br>This term is the one that you can compute (using the chain-rule) as a product<br>of the moving image gradient and the Transform Jacobian.<br><br>Strictly speaking:<br><br>m(i) is m( q(i),T) (see above)<br>
<br>d (m(i) ) / dT = d ( m( q(i), T ) / dT = GM( q(i) ) * J( p(i), T) [Eq.7]<br>
<br>Where <br><br>GM( q(i) ) is the Gradient of the moving image evaluated at point q(i)<br><br>where, again, q(i) is the point in the moving image that corresponds<br>to point p(i) of the fixed image.<br><br>J( p(i), T ) is the Jacobian matrix of the Transform that indicates<br>
how much the coordinates of the image point q(i) will change, with<br>respect to variations of each of the components of the Transform<br>parameters array.<br><br>Once you compute d( m(i) ) / dT, you can put this value back in <br>
equation 6., to compute the term d x(i) / dT<br> <br> d x(i) / dT = ( 1 / f(i) ) * d( m(i) ) / dT<br><br>and with this value, you can compute Equation 4. which will give you<br>the value of d M / d T.<br><br><br>To quickly answer your question: <br>
<br>YES, you have to multiply d m(i) / dT by ( 1 / f(i) ),<br>in order to get d x(i) / dT.<br><br><br>If you are looking for a detailed description of the calculus used<br>for computing Metric derivatives, you will find a discussion on this <br>
computation in the book:<br><br><br>"Insight into Images"<br><a href="http://www.amazon.com/Insight-into-Images-Segmentation-Registration/dp/1568812175">http://www.amazon.com/Insight-into-Images-Segmentation-Registration/dp/1568812175</a><br>
<br><br><br>BTW: As usual, <br>Once you get this metric to work, we strongly encourage you<br>to share it with the ITK community by submitting it as a paper<br>to the Insight Journal :-)<br><br><br><br> Regards,<br>
<br><br> Luis<br><br><br><br>-------------------------------------------------------------------------------------------------<br><div class="gmail_quote">On Thu, Jul 9, 2009 at 9:02 AM, Geoff Topping <span dir="ltr"><<a href="mailto:g_topping@hotmail.com">g_topping@hotmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
<div>
Hi all,<br><br>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.<br>
<br>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,<br>
<br>x(i) = m(i) / f(i)<br><br>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:<br><br>M = (1 / (n - 1)) * ( sum(x(i)*x(i)) - (1/n) * sum(x(i)) * sum(x(i)) )<br>
<br>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 )<br><br>where mean(x) is estimated from the sample data.<br><br>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.<br>
<br>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.<br><br>For each transform parameter s, I took the derivative of the metric M to get:<br>
<br>dM/ds = (2 / (n - 1)) * ( sum( x(i) * dx(i)/ds ) - (1/n) * sum(x(i)) * sum(dx(i)/ds) )<br><br>where dx(i)/ds is somewhat compliated to determine.<br><br>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.<br>
<br>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.<br>
<br>If someone could clarify that, or whether I'm understanding this at all, or perhaps that my derivation is right, that'd be nice.<br><br>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.<br>
<br>Thanks,<br><font color="#888888">Geoff Topping<br><hr>We are your photos. Share us now with <a href="http://go.microsoft.com/?linkid=9666045" target="_blank">Windows Live Photos.</a></font></div>
<br>_____________________________________<br>
Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
<br>
Visit other Kitware open-source projects at<br>
<a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
<br>
Please keep messages on-topic and check the ITK FAQ at: <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<br>
Follow this link to subscribe/unsubscribe:<br>
<a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
<br></blockquote></div><br>