[Insight-users] demons deformable filter
Corinne Mattmann
mattmaco at ee.ethz.ch
Thu, 15 Jan 2004 16:49:36 -0700
Hi Luis,
I have two more questions:
1) In itkDemonsRegistrationFunction the fixedGradientSquaredMagnitude
and update-value are calculated like this:
fixedGradientSquaredMagnitude += vnl_math_sqr( fixedGradient[j] ) *
m_FixedImageSpacing[j];
update[j] = speedValue * fixedGradient[j] *
vnl_math_sqr(m_FixedImageSpacing[j]) / denominator;
Why do this two formulas contain a multiplication with the ImageSpacing?
2) With it.GetCenterPixel() I get the deformation vector at the current
position - the same as when I go like
"m_DeformationField->GetPixel(index)". Do this two variables point to
the same deformation field or does one of them contain the updated
values (half actual iteration, half previous iteration)?
Thanks very much for all your answers,
Corinne
ps: The itkDemons5RegistrationFunction is working so far and with my
images I get better results than with the other demons function :) . I
just would like to make sure that my formula is right (at the moment I
don't multiply with the ImageSpacing -> see 1). If I delete this
multiplication in itkDemonsRegistrationFunction as well, the actual
registration gets better with my images but the reconstruction of the
image (apply the warp filter to the moving image) gets noisy.
-----Original Message-----
From: Luis Ibanez [mailto:luis.ibanez at kitware.com]
Sent: Wednesday, January 14, 2004 11:55 AM
To: Corinne Mattmann
Cc: insight-users at itk.org
Subject: Re: [Insight-users] demons deformable filter
Hi Corinne,
I assume that you are making reference to the documentation in
http://www.itk.org/Insight/Doxygen/html/classitk_1_1DemonsRegistrationFu
nction.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
>