[Insight-users] About the accuracy of up-sampling of deformation field from low resolution to high

Mon Apr 5 15:35:52 EDT 2010

Hi Luis,

   Your detailed explanation makes a lot more sense now. Thank you so much!

   I have another question about the accuracy of the upsampling (fitting the
lower resolution deformation field with higher resolution coefficients).
I noticed the discrepancy (magnitude of difference of displacement vector
between lower and higher resolution) at the grid locations of higher
resolution
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).
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).
Please see an attached histogram of these differences of all grid locations
(higher resolution).

  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
achieve the same level of accuracy compared to conventional matrix solving
approach (as your mentioned)? Why indeed the higher resolution cannot
achieve
a perfect accuracy to fit the lower resolution?

[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 bigwww.epfl.ch by Philippe Thevenaz

Thanks,
Mengda

2010/4/3 Luis Ibanez <luis.ibanez at kitware.com>

> Hi Mengda,
>
>
>               Thanks for your detailed comments.
>
>
> Let's look at the basic math to see
> if we can make sense of the code.
>
>
> The BSpline deformable transform class is representing
> a deformation field D at every point "p" in space.
>
> The field is computed as a linear combination of the
> BSpline basis (B):
>
>
>                  D(p) =  Sum_i  {  Bi(p-qi) * Vi  }
>
> where
>
> B(p) are the Bplines Basis
>
> p is a position vector (x,y,z).
>
> Vi is the value to be interpolated, in this case Vi is
>     a displacement vector located at the "i-th" node
>
> qi are the locations of each one of the grid nodes
>    in (x,y,z) coordinates.
>
> and "i" goes from 1 to N, where N is the number
> of nodes in the Bspline grid.
>
>
> What we are attempting to do in this code is to
> compute a finer grid of nodes that approximate
> well the same deformation field D(p).
>
>
> So we are looking for Wj elements that are each
> one a displacement vector (x,y,z), located at
> positions rj, (the location so the new grid nodes)
> such that
>
> (Eq1)      D(p) =  Sum_i  {  Bi(p-qi) * Vi  }
> (Eq2)              =  Sum_j  {  B'j(p-rj) * Wj }
>
>
> Where
>
> B'(p) are the basis of the finer BSpline grid.
>
> and j goes from 1 to M, where M is the number
> of grid points of the new BSpline grid.
>
> --
>
> The process is easier to understand by looking
> backwards first.  So, in order to find the coefficients
> "Wj" we need to solve a linear system
>
>         D(pk) =  Sum_j { B'j( pk - rj ) * Wj
>
> Where we know the value of the deformation field
> in a set of positions pk (which are coordinates in
> space.
>
> In order to be able to solve this linear system
> we need at least the same number of equations
> as the number of variables, therefore we need
> to know the deformation field D in M points, and
> the ideal locations of those points are just below
> the same grid nodes of the new grid.
>
> The Decomposition filter in lines 358-364 takes
> care of solving this part of the inverse problem,
> (the one in Eq2). once we provide a number M
> of values of the D(pk) deformation field.
>
> Now, in order to compute those values of D(pk),
> the ones that connect Eq1 with Eq2, we use the
> resample image filer lines 340-356. Note that
> the key here is that we use as interpolator a
>
>        BSplineResampleImageFunction
>
> This interpolator does the equivalent of the
> multiplication by B() functions and the sum
> over "i", as described in Eq1.
>
> Therefore, the outcome of the resample filter is
> not an image of coefficients (as it look at first
> sight), but indeed an image of deformations
> (D(p)) where every pixel matches the location
> of the new higher resolution BSpline grid.
>
>
> Please read the documentation of the class:
>
> http://www.itk.org/Doxygen/html/classitk_1_1BSplineResampleImageFunction.html
>
> and its superclass:
>
> http://www.itk.org/Doxygen/html/classitk_1_1BSplineInterpolateImageFunction.html
>
>
> and let us know if you still have any questions,
>
>
>     Thanks
>
>
>           Luis
>
>
>
>
> -----------------------------------------------------------------------------------------------
> On Thu, Apr 1, 2010 at 5:13 PM, Mengda Wu <wumengda at gmail.com> wrote:
> > Hi all,
> >
> >     I just would like to follow up with my last email. I think the
> > up-sampling code (line 340 to 375) is wrong.
> >
> >    Since the ResampleImageFilter is operated on the coefficients (not the
> > deformation field mentioned in my last email), what its output should
> > be already coefficients (not the deformation field) in higher resolution.
> > Using BSplineDecompositionImageFilter to get the coefficients again
> > does not make sense.
> >
> >    For reference, a method for up-sampling of the deformation field is
> > mentioned in Mattes's 2003 paper, in particular equations (17) and (18).
> >
> > Thanks,
> > Mengda
> >
> > ---------- Forwarded message ----------
> > From: Mengda Wu <wumengda at gmail.com>
> > Date: 2010/4/1
> > Subject: About the up-sampling of deformation field from low resolution
> to
> > high in Examples/DeformableRegistration6
> > To: Insight-users at itk.org
> >
> >
> > Hi all,
> >
> >    I have a question about conversion of deformation field from low
> > resolution to high in Examples/DeformableRegistration6. I do not
> understand
> > why we need both ResampleImageFilter and BSplineDecompositionImageFilter.
> I
> > think we only need to up sample the deformation field (per dimension)
> > using ResampleImageFilter with BSplineResampleImageFunction.
> >
> >   Then, why is the output of ResampleImageFilter fed into
> > BSplineDecompositionImageFilter and then put into the parameters? What
> does
> > BSplineDecompositionImageFilter do?
> > I cannot get this from the document of BSplineDecompositionImageFilter.
> >
> > Thanks,
> > Mengda
> >
> >
> > _____________________________________
> > Powered by www.kitware.com
> >
> > Visit other Kitware open-source projects at
> > http://www.kitware.com/opensource/opensource.html
> >
> > Kitware offers ITK Training Courses, for more information visit:
> > http://www.kitware.com/products/protraining.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/20100405/9931ab57/attachment-0001.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: fitting-histogram.png
Type: image/png
Size: 6287 bytes
Desc: not available
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100405/9931ab57/attachment-0001.png>