[Insight-users] Problems with 3D DeformableRegistrationFilter
Luis Ibanez
luis.ibanez@kitware.com
Tue May 4 03:30:52 EDT 2004
Hi Jayme,
1) Yes, the FEM Registration filter is a bit unusual in the sense
that it groups a lot of functionality. It is almost an application
in itself.
The reason why it cannot use the basic registration framework is
that its approach doesn't quite match with the structure of an
optimizer, metric, ... etc.
You may want to add a couple of "Feature requests" in the
bug tracker, in particular regarding the emision of Events.
http://www.itk.org/Bug/
Note that refreshing the deformation field at each iteration
maybe too costly. This filter represents the deformation using
the displacements of the mesh nodes. The computation of the
full deformation field requires to sample pixels over the mesh
and will produce a performace hit if performed at each iteration.
You could however request to have a callable method that computes
the resample deformation field. In this case, your observer could
force the computation only when needed.
2) That's an excellent question. :-)
To give you a Mathematical answer:
"Yes,
sure there is an intelligent way of selecting parameters."
The elasticity and rhoc parameters depend on the intensity of your
images, their similarity and the type of image metric that you use.
If you think in physics, you will find a more natural approach to
parameter selection. Set a minimal FEM grid, for example 3x3x3 nodes,
and make the metric print out all the derivative values computed at
the nodes. Those are the forces excerted over the nodes. Start with
a low elasticity parameter, which should result in large
deformations. Then restart the registration with progressively higher
elasticity values until you find displacement that are in the range
of your expected deformations.
I'm affraid that you will need to read at least an introduction to FEM
concepts in order to make good use of this filter. Otherwise it will
quite difficult for you to follow the relationships between the
parameters.
Regards,
Luis
----------------------
Jayme Kosior wrote:
> Hi Luis,
>
> I apologize for the lateness of my response. My inbox was cluttered and
> I somehow missed your reply and just found it now as I was cleaning it
> out. :)
>
> My problem with the pixel type was that I was using a reader that was
> setup to read 8-bit images and yet my data was 32-bit. I was wrong about
> the Rat Lung example and it is setup correctly. I do have 3D FEM
> deformable registration running (attached is the code that I am using
> along with one of the configuration files), but I have been having a
> *very* difficult time determining how to drive the parameters. Here are
> some questions that I have:
>
> 1. FEMRegistrationFilter Design Question
>
> I am unclear as to why the FEMRegistrationFilter does not use the
> ImageRegistrationMethod framework and instead inherits from
> ImageToImageFilter. Since it is implemented as a regular
> ImageToImageFilter, there are no methods available to supply an observer
> (part of the Command/Observer pattern used by other registration
> methods). For example, this means that I cannot change interpolators
> between multi-resolution levels or get feedback (besides the
> command-line feedback) For example, I would like to get the current
> warped image to visually monitor registration between the long
> iterations periods. From what I can tell, the multi-resolution scheme is
> handled inside the FEMRegistrationFilter and there is no way to access
> it. Also, the FEMRegistrationFilter interface does not provide the
> ability to choose other metrics (in particular Mattes mutual information
> metric), nor the ability to modify the the metrics themselves (i.e the
> number of spatial samples for Mutual information metric). Perhaps I am
> missing something, but I haven't been able to answer this question by
> reading the documentation or the source code.
>
>
> 2. Is there an intelligent way to determine the optimal parameters for a
> particular application?
>
> The parameter space is quite large. In order to understand the use of
> each parameter, I have been trying to balance the theory behind FEM
> deformable registration with the ITK implementation (i.e. how to set
> rhoC to 1e3, elasticity to 1e2 etc. for a given application). So far, I
> have been using the good old brute force approach and have been working
> the computers overtime, but I have not had much success. Currently, I am
> trying to register two FLAIR volumes from different patients. First, I
> use an Affine registration and then apply the FEM registration process.
>
> I have read previous threads and the source code to get an idea of how
> to set them. The main parameters I am experimenting with are:
> - elasticity (1e1 - 1e4)
> - density/rhoc (1e1 - 1e4)
> - time step (1e0 - 1e3)
>
> The other parameters I set as suggested in previous threads (see
> configuration file below). Qualitatively by visual inspection and using
> normalized mean squared difference images, I observed that there was
> very little deformation when the elasticity/density ratio was ~1 as
> recommended in the threads (in fact, the differences I observed could
> quite possibly be from the random sampling of the mutual information
> metric). I found that density needed to be less than 1e3, elasticity
> around 1e1-1e2, and the time step around 1e2 to get a noticeable warp.
> However, the images were not registered very well.
>
> Is there a better way to find the optimal parameters for a particular
> application? For example, with the anisotropic diffusion filter you can
> select the conductance parameter K to be around the gradient magnitude
> produced by the noise level that you wish to remove. Are there any
> starting points for the FEM registration methods?
>
> The difficulty that I am having with FEM registration is in applying FEM
> principles (which are mainly based on physical principles) with
> intensity based registration. I don't have a FEM background.
>
> I appreciate you taking the time to help me. Thank you.
>
> Sincerely,
>
> Jayme Kosior
>
> NOTE: The FEM3D file is one that I modified slightly from the CVS
> repository. Namely, I just changed it to use 3D elements and work with
> 32-bit volumes only.
>
>
>
> ------------------------------------------------------------------------
>
> % Configuration file #1 for DeformableRegistration1.cxx
> %
> % This example demonstrates the setup of a basic registration
> % problem that does NOT use multi-resolution strategies. As a
> % result, only one value for the parameters between
> % (# of pixels per element) and (maximum iterations) is necessary.
> % If you were using multi-resolution, you would have to specify
> % values for those parameters at each level of the pyramid.
> %
> % Note: the paths in the parameters assume you have the traditional
> % ITK file hierarchy as shown below:
> %
> % ITK/Insight/Examples/Registration/DeformableRegistration1.cxx
> % ITK/Insight/Examples/Data/RatLungSlice*
> % ITK/Insight-Bin/bin/DeformableRegistration1
> %
> % ---------------------------------------------------------
> % Parameters for the single- or multi-resolution techniques
> % ---------------------------------------------------------
> 3 % Number of levels in the multi-res pyramid (1 = single-res)
> 3 % Highest level to use in the pyramid
> 1 1 1 % Scaling at lowest level of pyramid
> 10 10 10 % Number of pixels per element
> 1.e1 1.e1 1.e1 % Elasticity (E)
> 1.e2 1.e2 1.e2 % Density x capacity (RhoC)
> 1 1 1 % Image energy scaling (gamma) - sets gradient step size
> 3 3 3 % NumberOfIntegrationPoints
> 3 3 3 % WidthOfMetricRegion
> 30 30 30 % MaximumIterations
> % -------------------------------
> % Parameters for the registration
> % -------------------------------
> 3 0.99 % Similarity metric (0=mean sq, 1 = ncc, 2=pattern int, 3=MI, 5=demons) 2nd num changes NormalizeGradient
> 1.0 % Alpha
> 1 % DescentDirection (1 = max, 0 = min)
> 2 % DoLineSearch (0=never, 1=always, 2=if needed)
> 1.e2 % TimeStep
> 0.5 % Landmark variance
> 0 % Employ regridding / enforce diffeomorphism ( >= 1 -> true)
> % ----------------------------------
> % Information about the image inputs
> % ----------------------------------
> 256 % Nx (image x dimension)
> 256 % Ny (image y dimension)
> 12 % Nz (image z dimension - not used if 2D)
> ../CISGResults/Rigid/KFmToHFm_Rigid_NoZ/32-bit/CISG-KtoH-sinc-Rigid.hdr % moving image
> ../images3D/images/HFLAIR30m.hdr % fixed image
> % -------------------------------------------------------------------
> % The actions below depend on the values of the flags preceding them.
> % For example, to write out the displacement fields, you have to set
> % the value of WriteDisplacementField to 1.
> % -------------------------------------------------------------------
> 0 % UseLandmarks? - read the file name below if this is true
> - % LandmarkFileName
> ./finalResults/MI/test2/KtoH-ITK.hdr % ResultsFileName (prefix only)
> 1 % WriteDisplacementField?
> ./finalResults/MI/test2/KtoH-ITK-VecFld % DisplacementsFileName (prefix only)
> 0 % ReadMeshFile?
> - % MeshFileName
> END
>
>
>
> ------------------------------------------------------------------------
>
>
>
>
>
>
> On Apr 10, 2004, at 10:25, Luis Ibanez wrote:
>
>
> Hi Jayme,
>
> Thanks for sending your code.
>
> It is not clear why the pixel type should make
> any difference, specially given the fact that
> after the readers there is a RescaleImageFilter.
>
> >From your message it is not yet clear to me if
> you got the registration working or not.
> If you did, it will be nice to add your file
> as an example into
>
> Insight/Examples/Registration
>
> along with the configuration input file.
>
>
> Please let us know about the status of your
> registration process.
>
>
>
> Thanks,
>
>
> Luis
>
>
> -------------------
> Jayme Kosior wrote:
>
> Hi Luis,
> I believe I figured out what was causing the registration to
> fail. The data that I was using was in 8-bit format. I converted
> it into 32-bit format and I changed the line:
> //typedef itk::Image<unsigned char, 3> fileImage3DType;
> to:
> typedef itk::Image<float, 3> fileImage3DType;
> since the image pixel type was defined as:
> typedef itk::Image<float, 3> Image3DType;
> The RatLungSlice1.mha and RatLungSlice2.mha provided with the
> ITK examples were both 8-bits, and so I had assumed that there
> was a type cast/conversion after the RescaleIntensityImageFilter
> ran, since it was setup as:
> typedef
> itk::RescaleIntensityImageFilter<fileImage3DType,Image3DType>
> FilterType;
> One thing that is strange is that 2D deformable registration
> worked (i.e. produced a noticeable deformation) using 8-bit
> images as input, but the same is not true for 3D images. As I
> said before, all that happened is that the moving images would
> be transformed into the fixed image space, but no deformation
> would occur.
> In order to test to see if the registration is working, I have
> been registering 3D MR FLAIR to 3D Perfusion data sets. They are
> both at slightly different angles (~8 degress x-y plane, ~20
> degrees z-plane, plus EPI distortion), but I did not register
> them using a rigid body method (i.e. VersorRigid3DTransform -->
> in fact, that's what I'm working on right now). I realize that
> rigid registration should be performed first, but I just wanted
> to see a noticeable deformation in the images before I started
> to pursue the FEMRegistration method any further (just to be
> sure I could get it working).
> Anyways, here is the code that I was using. All that has really
> changed is that I am using 3D elements. I left the original
> comments.
> Thank you.
> Jayme
> On Mar 22, 2004, at 16:36, Luis Ibanez wrote:
>
>
> Hi Jayme,
>
> At first sight this looks like you forgot to switch
> some of the FEM components from 2D to 3D.
>
> Could you please post your modified version of the
> DeformableRegistration1.cxx example.
>
> Once we debug it, we will actually like to add it
> as an extra example in the same directory.
>
>
> Thanks
>
>
> Luis
>
>
> --------------------------
> Jayme Kosior wrote:
>
> Hello.
> I can successfully run DeformableRegistration1 to
> register 2D images. I have modified
> DeformableRegistration1 and the parameter file
> accordingly to perform 3D deformable registration. After
> performing the registration, only the first image in the
> 3D dataset is deformed. The rest of the images are not
> modified. The only change in the rest of the registered
> data set is that it is now in the fixed image coordinate
> space.
> I'm not sure if the registration is not warping the
> images properly, or if maybe it is not writing them
> correctly. Of course, maybe I'm not doing something
> correctly. :)
> Here is the parameter file that I am using.
> Thank you in advance for your help.
> Jayme
> Setup Info: 700 Mhz Athlon, Red Hat Linux 9.0, ITK
> version 1.6 release, 2 8-bit 3D data sets in Analyze
> format at input.
> % Configuration file #1 for DeformableRegistration1.cxx
> %
> % This example demonstrates the setup of a basic
> registration
> % problem that does NOT use multi-resolution strategies.
> As a
> % result, only one value for the parameters between
> % (# of pixels per element) and (maximum iterations) is
> necessary.
> % If you were using multi-resolution, you would have to
> specify
> % values for those parameters at each level of the pyramid.
> %
> % Note: the paths in the parameters assume you have the
> traditional
> % ITK file hierarchy as shown below:
> %
> %
> ITK/Insight/Examples/Registration/DeformableRegistration1.cxx
>
> % ITK/Insight/Examples/Data/RatLungSlice*
> % ITK/Insight-Bin/bin/DeformableRegistration1
> %
> % ---------------------------------------------------------
> % Parameters for the single- or multi-resolution techniques
> % ---------------------------------------------------------
> 1 % Number of levels in the multi-res pyramid (1 =
> single-res)
> 1 % Highest level to use in the pyramid
> 1 1 1 % Scaling at lowest level of pyramid
> 8 % Number of pixels per element
> 1.e4 % Elasticity (E)
> 1.e3 % Density x capacity (RhoC)
> 1 % Image energy scaling (gamma) - sets gradient step size
> 6 % NumberOfIntegrationPoints
> 4 % WidthOfMetricRegion
> 50 % MaximumIterations
> % -------------------------------
> % Parameters for the registration
> % -------------------------------
> 4 1.0 % Similarity metric (0=mean sq, 1 = ncc, 2=pattern
> int, 3=MI, 5=demons) 2nd num changes NormalizeGradient
> 0.5 % Alpha
> 0 % DescentDirection (1 = max, 0 = min)
> 2 % DoLineSearch (0=never, 1=always, 2=if needed)
> 1.e1 % TimeStep
> 0.5 % Landmark variance
> 0 % Employ regridding / enforce diffeomorphism ( >= 1 ->
> true)
> % ----------------------------------
> % Information about the image inputs
> % ----------------------------------
> 256 % Nx (image x dimension)
> 256 % Ny (image y dimension)
> 12 % Nz (image z dimension - not used if 2D)
> ./images-8bit/CFLAIR30m.hdr % moving image
> ./images-8bit/CPWI0d_brain3.hdr % fixed image
> %
> -------------------------------------------------------------------
>
> % The actions below depend on the values of the flags
> preceding them.
> % For example, to write out the displacement fields, you
> have to set
> % the value of WriteDisplacementField to 1.
> %
> -------------------------------------------------------------------
>
> 0 % UseLandmarks? - read the file name below if this is
> true
> - % LandmarkFileName
> ./results/PatientC-FtoP-8bit % ResultsFileName (prefix
> only)
> 1 % WriteDisplacementField?
> ./results/OtherPatientC-FtoP_disp %
> DisplacementsFileName (prefix only)
> 0 % ReadMeshFile?
> - % MeshFileName
> END
> _______________________________________________
> Insight-users mailing list
> Insight-users@itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
>
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users@itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users@itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list