Hi Luis,<br><br> Your detailed explanation makes a lot more sense now. Thank you so much!<br><br> I have another question about the accuracy of the upsampling (fitting the lower resolution deformation field with higher resolution coefficients). <br>
I noticed the discrepancy (magnitude of difference of displacement vector between lower and higher resolution) at the grid locations of higher resolution <br>could be quite large. I have a 3D volume (with isotropic pixel spacing 0.78 mm) and the b-spline grid spacing (higher resolution) is about 15 pixels (that is 0.78x15 = 11.7 mm). <br>
This difference can be as large as 45 mm in some locations. The b-spline grid spacing (lower resolution) is about 30 pixels (that is 0.78x30 = 23.4 mm). <br>Please see an attached histogram of these differences of all grid locations (higher resolution).<br>
<br> I am wondering if you can provide any clue of that. I tried to understand how this fitting works. Is the filtering approach (as in the papers below) able to <br>achieve the same level of accuracy compared to conventional matrix solving approach (as your mentioned)? Why indeed the higher resolution cannot achieve<br>
a perfect accuracy to fit the lower resolution?<br><br>[1] M. Unser, "Splines: A Perfect Fit for Signal and Image Processing,"
IEEE Signal Processing Magazine, vol. 16, no. 6, pp. 22-38, November
1999. [2] M. Unser, A. Aldroubi and M. Eden, "B-Spline Signal
Processing: Part I--Theory," IEEE Transactions on Signal Processing,
vol. 41, no. 2, pp. 821-832, February 1993. [3] M. Unser, A. Aldroubi
and M. Eden, "B-Spline Signal Processing: Part II--Efficient Design and
Applications," IEEE Transactions on Signal Processing, vol. 41, no. 2,
pp. 834-848, February 1993. And code obtained from <a href="http://bigwww.epfl.ch">bigwww.epfl.ch</a> by
Philippe Thevenaz<br><br>Thanks,<br>Mengda<br><br><div class="gmail_quote">2010/4/3 Luis Ibanez <span dir="ltr"><<a href="mailto:luis.ibanez@kitware.com">luis.ibanez@kitware.com</a>></span><br><blockquote class="gmail_quote" style="margin: 0pt 0pt 0pt 0.8ex; border-left: 1px solid rgb(204, 204, 204); padding-left: 1ex;">
Hi Mengda,<br>
<br>
<br>
Thanks for your detailed comments.<br>
<br>
<br>
Let's look at the basic math to see<br>
if we can make sense of the code.<br>
<br>
<br>
The BSpline deformable transform class is representing<br>
a deformation field D at every point "p" in space.<br>
<br>
The field is computed as a linear combination of the<br>
BSpline basis (B):<br>
<br>
<br>
D(p) = Sum_i { Bi(p-qi) * Vi }<br>
<br>
where<br>
<br>
B(p) are the Bplines Basis<br>
<br>
p is a position vector (x,y,z).<br>
<br>
Vi is the value to be interpolated, in this case Vi is<br>
a displacement vector located at the "i-th" node<br>
<br>
qi are the locations of each one of the grid nodes<br>
in (x,y,z) coordinates.<br>
<br>
and "i" goes from 1 to N, where N is the number<br>
of nodes in the Bspline grid.<br>
<br>
<br>
What we are attempting to do in this code is to<br>
compute a finer grid of nodes that approximate<br>
well the same deformation field D(p).<br>
<br>
<br>
So we are looking for Wj elements that are each<br>
one a displacement vector (x,y,z), located at<br>
positions rj, (the location so the new grid nodes)<br>
such that<br>
<br>
(Eq1) D(p) = Sum_i { Bi(p-qi) * Vi }<br>
(Eq2) = Sum_j { B'j(p-rj) * Wj }<br>
<br>
<br>
Where<br>
<br>
B'(p) are the basis of the finer BSpline grid.<br>
<br>
and j goes from 1 to M, where M is the number<br>
of grid points of the new BSpline grid.<br>
<br>
--<br>
<br>
The process is easier to understand by looking<br>
backwards first. So, in order to find the coefficients<br>
"Wj" we need to solve a linear system<br>
<br>
D(pk) = Sum_j { B'j( pk - rj ) * Wj<br>
<br>
Where we know the value of the deformation field<br>
in a set of positions pk (which are coordinates in<br>
space.<br>
<br>
In order to be able to solve this linear system<br>
we need at least the same number of equations<br>
as the number of variables, therefore we need<br>
to know the deformation field D in M points, and<br>
the ideal locations of those points are just below<br>
the same grid nodes of the new grid.<br>
<br>
The Decomposition filter in lines 358-364 takes<br>
care of solving this part of the inverse problem,<br>
(the one in Eq2). once we provide a number M<br>
of values of the D(pk) deformation field.<br>
<br>
Now, in order to compute those values of D(pk),<br>
the ones that connect Eq1 with Eq2, we use the<br>
resample image filer lines 340-356. Note that<br>
the key here is that we use as interpolator a<br>
<br>
BSplineResampleImageFunction<br>
<br>
This interpolator does the equivalent of the<br>
multiplication by B() functions and the sum<br>
over "i", as described in Eq1.<br>
<br>
Therefore, the outcome of the resample filter is<br>
not an image of coefficients (as it look at first<br>
sight), but indeed an image of deformations<br>
(D(p)) where every pixel matches the location<br>
of the new higher resolution BSpline grid.<br>
<br>
<br>
Please read the documentation of the class:<br>
<a href="http://www.itk.org/Doxygen/html/classitk_1_1BSplineResampleImageFunction.html" target="_blank">http://www.itk.org/Doxygen/html/classitk_1_1BSplineResampleImageFunction.html</a><br>
<br>
and its superclass:<br>
<a href="http://www.itk.org/Doxygen/html/classitk_1_1BSplineInterpolateImageFunction.html" target="_blank">http://www.itk.org/Doxygen/html/classitk_1_1BSplineInterpolateImageFunction.html</a><br>
<br>
<br>
and let us know if you still have any questions,<br>
<br>
<br>
Thanks<br>
<br>
<br>
Luis<br>
<br>
<br>
<br>
-----------------------------------------------------------------------------------------------<br>
<div><div></div><div class="h5">On Thu, Apr 1, 2010 at 5:13 PM, Mengda Wu <<a href="mailto:wumengda@gmail.com">wumengda@gmail.com</a>> wrote:<br>
> Hi all,<br>
><br>
> I just would like to follow up with my last email. I think the<br>
> up-sampling code (line 340 to 375) is wrong.<br>
><br>
> Since the ResampleImageFilter is operated on the coefficients (not the<br>
> deformation field mentioned in my last email), what its output should<br>
> be already coefficients (not the deformation field) in higher resolution.<br>
> Using BSplineDecompositionImageFilter to get the coefficients again<br>
> does not make sense.<br>
><br>
> For reference, a method for up-sampling of the deformation field is<br>
> mentioned in Mattes's 2003 paper, in particular equations (17) and (18).<br>
><br>
> Thanks,<br>
> Mengda<br>
><br>
> ---------- Forwarded message ----------<br>
> From: Mengda Wu <<a href="mailto:wumengda@gmail.com">wumengda@gmail.com</a>><br>
> Date: 2010/4/1<br>
> Subject: About the up-sampling of deformation field from low resolution to<br>
> high in Examples/DeformableRegistration6<br>
> To: <a href="mailto:Insight-users@itk.org">Insight-users@itk.org</a><br>
><br>
><br>
> Hi all,<br>
><br>
> I have a question about conversion of deformation field from low<br>
> resolution to high in Examples/DeformableRegistration6. I do not understand<br>
> why we need both ResampleImageFilter and BSplineDecompositionImageFilter. I<br>
> think we only need to up sample the deformation field (per dimension)<br>
> using ResampleImageFilter with BSplineResampleImageFunction.<br>
><br>
> Then, why is the output of ResampleImageFilter fed into<br>
> BSplineDecompositionImageFilter and then put into the parameters? What does<br>
> BSplineDecompositionImageFilter do?<br>
> I cannot get this from the document of BSplineDecompositionImageFilter.<br>
><br>
> Thanks,<br>
> Mengda<br>
><br>
><br>
</div></div><div class="im">> _____________________________________<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>
</div><div class="im">> Kitware offers ITK Training Courses, for more information visit:<br>
> <a href="http://www.kitware.com/products/protraining.html" target="_blank">http://www.kitware.com/products/protraining.html</a><br>
><br>
> Please keep messages on-topic and check the ITK FAQ at:<br>
</div>> <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
<div class="im">><br>
> Follow this link to subscribe/unsubscribe:<br>
</div><div><div></div><div class="h5">> <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
><br>
><br>
</div></div></blockquote></div><br>