[Insight-users] demons deformable filter

[Luis Ibanez]

Hi Corinne,

I assume that you are making reference to the documentation in

http://www.itk.org/Insight/Doxygen/html/classitk_1_1DemonsRegistrationFunction.html


  "This class encapsulate the PDE which drives the demons registration
   algorithm. It is used by DemonsRegistrationFilter to compute the
   output deformation field which will map a moving image onto a a
   fixed image."


The deformation field computed by the DemonsRegistrationFilter is
the vector field that can be used by the WrapImageFilter in order
to deform (resample) the moving image.

This deformation field is read using a point P the fixed image
coordinates and it tells what displacement you have to apply
to P in order to find its corresponding point Q in the moving
image coordinates.

            Q = P + DeformationField( P )

Where P is a point in the fixed  image coordinate system
and   Q is a point in the moving image coordinate system


so... the deformation field maps position from the fixed image
coordinate system into the moving image coordinate system, and
by doing this, it is the field which can be used to "map a moving
image onto a fixed image".



Please let us know if you have further questions.



Thanks


   Luis



------------------------
Corinne Mattmann wrote:

> Hi Luis,
> 
> Thanks very much for your answer. At the moment there's just one thing I
> don't understand:
> In itkDemonsRegistrationFunction I have the variable "index" which
> contains the position in the fixed image I am looking at that moment. To
> get the intensity value of the moving image at the surrounding points I
> have to know which points of the original moving image are mapped to 
> that surrounding points.
> In point B) you wrote that I just have to take those surrounding points
> and map them through the deformation field. And here is what I don't
> understand: In the description of the DemonsRegistrationFilter I've read
> that the "deformation field maps the moving image onto the fixed image".
> If I understand this statement correctly, this means that the
> deformation field contains a vector in each point, telling me where the
> intensity value of the moving image is positioned after the
> registration.
> Therefore to get the positions in the moving image of the mentioned
> surrounding points I would somehow have to go backwards through the
> deformation field. Is this right or do I misunderstand something? And if
> I understand it right how can I then get the mapped moving values
> without calculating tho whole image in InitializeIteration()?
> 
> Thanks very much,
> Corinne
> 
> 
> -----Original Message-----
> From: Luis Ibanez [mailto:luis.ibanez at kitware.com] 
> Sent: Tuesday, January 13, 2004 4:37 PM
> To: Corinne Mattmann
> Cc: insight-users at itk.org
> Subject: Re: [Insight-users] demons deformable filter
> 
> 
> 
> Hi Corinne,
> 
> 
> You probably want to start by creating new classes
> 
> 
>      itkDemons5RegistrationFunction
>      itkDemons5RegistrationFilter
> 
> 
> You can simply copy & rename them from the
> currently existing DemonsRegistration classes,
> and then do the following modifications on them:
> 
> 
> 
> 
> A) itk::Demons5RegistrationFunction still derives
>     from itk::PDEDeformableRegistrationFunction so
>     it has access to the protected member variable:
> 
>          m_DeformationField
> 
> 
>     This variable is not set by default, so you have
>     to modify the itkDemons5RegistrationFunction and
>     add in the InitializeIteration() method the following
>     lines
> 
> 
>>  // update variables in the equation object
>>  DemonsRegistrationFunctionType *f = 
>>    dynamic_cast<DemonsRegistrationFunctionType *>
>>    (this->GetDifferenceFunction().GetPointer());
>>
>>  if ( !f )
>>    {
>>    itkExceptionMacro(<<"FiniteDifferenceFunction not of type
> 
> DemonsRegistrationFunctionType");
> 
>>    }
>>
>>  f->SetDeformationField( this->GetDeformationField() );
> 
>  >
>      This will set the pointer on the function.
> 
>      Make sure that this new lines of code are BEFORE
>      the line invoking
> 
>             Superclass::InitializeIteration();
> 
> 
> B) Now, you can modify itkDemons5RegistrationFunction
>     for computing an estimation of the gradient by using
>     finite differences. For this you take the point on
>     the fixed image coordinate frame where you want to
>     compute the gradient of the deformed moving image.
> 
>     From this point you generate 2 x N points
>     (N = image dimension).  Those points are computed
>     as P +- delta[i].
> 
>     For example in 2D you will get points like:
> 
>                  Py+
>                  |
>                  |
>          Px- ----P----Px+
>                  |
>                  |
>                  Py-
> 
>     Then you take those points and map them through
>     the deformation field, so you will get points
>     Qx+,Qx-,Qy+,Qy-  (and so on if you are in N>2).
> 
>     With these points you can use the interpolator
>     and get intensity values from the moving image
>     and finally compute the gradient vector as
> 
>        Gx = (Qx+ - Qx-) / (2.0 * xspacing)
>        Gy = (Qy+ - Qy-) / (2.0 * yspacing)
>         (and so on if you are in N > 2).
> 
> 
> 
> 
> 
> This may not be very elegant but should give you
> a reasonable estimation of the gradient of the
> deformed moving image.
> 
> 
> BTW If you get it working and you are willing to
> contribute this code back to ITK we will be happy
> to include it in the toolkit. Probably with better
> names for the files.
> 
> 
> 
> Regards,
> 
> 
>    Luis
> 
> 
> --------------------------
> Corinne Mattmann wrote:
> 
> 
>>Hi,
>>
>>I am using the demons deformable filter but as I have images with just
> 
> 
>>one "material" - therefore basically black-white-images with smoothed 
>>transitions - the registration is not perfect. In fact, the quality of
> 
> 
>>the registration depends on the similarity of the two initial images. 
>>That's why I would like to try another demons-formula taken out of 
>>Jean-Philippe Thirion's paper "Fast Non-Rigid Matching of 3D Medical 
>>Images" (formula 5). For that formula you need to calculate the 
>>gradient of the reconstructed moving image at each iteration. Can you 
>>tell me how I can get to this gradient in the function
>>ComputeUpdate(...) in itkDemonsRegistrationFunction.txx?
>>
>>Thanks very much,
>>Corinne
>>
>>_______________________________________________
>>Insight-users mailing list
>>Insight-users at itk.org 
>>http://www.itk.org/mailman/listinfo/insight-users
>>
> 
> 
> 
> 
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>