[Insight-users] Follow up: About the up-sampling of deformation field from low resolution to high in Examples/DeformableRegistration6

[Luis Ibanez]

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
>
>