[Insight-users] All zeros displacement field when using DeformableRegistration1

Barbara Garita bgarita at ufl.edu
Wed Nov 10 11:48:13 EST 2004


Hello Luis:

Thanks for you help, I have some question about your response:

1) What is the relation between deformation field and the image pixel
spacing.  How can this information  eliminate the extreme deformations as
the possible cause of the black resampled image?

2) I used a  histogram count.

I am trying to use the code for simpler images and try to go from there.  I
will try the Demons algorithm as well.
Thanks again, Barbara


----- Original Message ----- 
From: "Luis Ibanez" <luis.ibanez at kitware.com>
To: "Barbara" <bgarita at ufl.edu>
Cc: <insight-users at itk.org>
Sent: Tuesday, November 09, 2004 5:07 PM
Subject: Re: [Insight-users] All zeros displacement field when using
DeformableRegistration1


>
> Hi Barbara,
>
> Thanks for posting your additional information.
>
> First of all, Note that the FEM model that is used here for
> image registration is *totally* disconnected from any physical
> reality.
>
> The coefficients that you set for the elastic model *are not*
> related at all to the physical properties of the tissue present
> in the images (e.g. bone, muscle...).
>
> The elastic model is used *just* as a paradigm for managing a
> deformation field under a regularization constraint.
>
> You are free to set the FEM parameters in any way that lead to
> the proper level of flexibility in the model grid and with the
> proper level of regularization.
>
> -----------
>
>  From the image data that you posted,
>
> 1.1 - 1.3) The values of the deformation field seem to be well
>             in the expected range for an image of pixel spacing
>             0.32 ^ 3.
>
> This eliminates extreme deformations as the possible cause of
> your black resampled image.
>
> The next suspects in line are:
>
>
> A) The resampling process
> B) The file writing process
> C) The viewer that you are using to conclude that the image is black.
>
> since from your code, it seems that (A) and (B) are just the
> same that we use in the nightly testing.... the remaining main
> suspect is (C)...
>
> So... please let us know what kind of visualization or measuring
> program are you using that lead you to conclude that your image
> is filled with zeros...
>
> (e.g. it is a common mistake to view images of 16bits with viewers
> that do not normalize intensities, and therefore produce the impression
> that the image is black.)
>
>
> Please let us know what you are using,
>
>
>
> Thanks
>
>
>
>     Luis
>
>
>
>
> -----------------
> Barbara wrote:
>
> > Hello Luis,
> >
> > 1.1)  See attachement. 1.2) pixel spacing is 0.32 x 0.32 x 0.32
> > 1.3) from my readings of the scalar images (vector components in x, y
> > and z) there were a few values in the -100 for both x, and y.  The
> > majority of the values around 90% remaind in and itnerval between -10 to
> > 10 (mm).  These histogram measures are for the registration of two
> > images which their origens were a distrance of 6 mm apart in the x -y
> > plane.
> >
> > Additional comments:
> > a1) When doing this nonrigid registration I get huge energy values, in
> > the order of :
> > ....energy 54939155425.463
> >     interpolating vector field of size [38, 90, 106] interpolation done
> >     min E 0.000 delt E -54939155425.463 iter 1
> >     load sizes [38, 90, 106]  image [38, 90, 106]
> >     energy 69884204337.691
> >     interpolating vector field of size [38, 90, 106] interpolation done
> >     min E 54939155425.463 delt E -14945048912.228 iter 2
> >     load sizes [38, 90, 106]  image [38, 90, 106]
> >     energy 70060092642.231 .....
> >
> > a2) I haven't teaked with the material properties.  I wonder if I should
> > use a E (modulus of elasticity), v (poisson ratio), density  close to
> > real bone values.
> > a3) I was looking at my warperfilter code, were is the DefaultPixel
> > Value() set.
> >
> > Thanks, Barbara
> >
> >
> >
> > Luis Ibanez wrote:
> >
> >>
> >> Hi Barbara,
> >>
> >> 1) One possibility is that there may be a mismatch in
> >>    the magnitude of the displacement vectors relative
> >>    to the pixel spacing in the image.
> >>
> >>
> >>    1.1)   Could you please post your code to the list
> >>    1.2)   Please indicate the pixel spacing of your image
> >>    1.3)   Please indicate the max values of the deformation field
> >>
> >>
> >>    You may want to try setting the DefaultPixelValue() to something
> >>    different from zero, in order to verify if what you are seeing
> >>    is your image being mapped outside of the field of view or
> >>    if there is an intensity corruption on the resampling.
> >>
> >>
> >> 2) FEM is fine for what you are doing.
> >>    You could try BSplines too, both techniques are in principle ok
> >>    for this application.
> >>
> >>    Note that if the Vertebrae are really close in shape you may
> >>    also want to try the Demons deformable registration and the
> >>    LevelSet deformable registration.
> >>
> >>
> >>
> >>
> >> Please let us know about the questions above,
> >>
> >>
> >>   Thanks
> >>
> >>
> >>     Luis
> >>
> >>
> >>
> >> ----------------------
> >> Barbara Garita wrote:
> >>
> >>> Hello All:
> >>>
> >>> I will like some advice concerning two issues about non-rigid
> >>> registration:
> >>> 1) I am using DeformableRegistration to register two 3D images (data
> >>> sets of human vertebrae). After I made the necessary accommodation
> >>> for 3D to the parameterlist and the code, I ran the
> >>> DeformableRegistration1, I checked for the correctness of the
> >>> vectordeformationfield scalar components  by extracting them and
> >>> analyzing their values.  However, when I use the Warperfilter and the
> >>> vectordeformationfield to get back to a scalar image I get a image of
> >>> the right size but it is all zeros.  2) I am trying to non-rigidly
> >>> register two image data set of human vertebrae (I am interested in
> >>> having a displacement field that gives me correspondence between the
> >>> surfaces (edges, contours...) of both vertebrae, I am not really
> >>> interested in the information inside the vertebrae--), should I use
> >>> another filter, like affine or B spline....
> >>>
> >>> Thanks for you advice, Barbara
> >>>
> >>>
> >>>
>
>>> ------------------------------------------------------------------------
> >>>
> >>> _______________________________________________
> >>> Insight-users mailing list
> >>> Insight-users at itk.org
> >>> http://www.itk.org/mailman/listinfo/insight-users
> >>
> >>
> >>
> >>
> >>
> >>
> >>
> >>
> >
> >
> > ------------------------------------------------------------------------
> >
> >
/*=========================================================================
> >
> >   Program:   Insight Segmentation & Registration Toolkit
> >   Module:    $RCSfile: DeformableRegistration1.cxx,v $
> >   Language:  C++
> >   Date:      $Date: 2004/05/12 23:21:36 $
> >   Version:   $Revision: 1.25 $
> >
> >   Copyright (c) Insight Software Consortium. All rights reserved.
> >   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for
details.
> >
> >      This software is distributed WITHOUT ANY WARRANTY; without even
> >      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
> >      PURPOSE.  See the above copyright notices for more information.
> >
> >
=========================================================================*/
> >
> >
> > #include "itkImageFileReader.h"
> > #include "itkImageFileWriter.h"
> > //#include "itkMetaImageIO.h"
> >
> > #include "itkRescaleIntensityImageFilter.h"
> > #include "itkHistogramMatchingImageFilter.h"
> >
> > //  Software Guide : BeginLatex
> > //
> > // The finite element (FEM) library within the Insight Toolkit can be
> > // used to solve deformable image registration problems.  The first step
in
> > // implementing a FEM-based registration is to include the appropriate
> > // header files.
> > //
> > //  \index{Registration!Finite Element-Based}
> > //
> > //  Software Guide : EndLatex
> >
> >
> > // Software Guide : BeginCodeSnippet
> > #include "itkFEM.h"
> > #include "itkFEMRegistrationFilter.h"
> >
> > // Software Guide : EndCodeSnippet
> >
> > //#include "itkFEMFiniteDifferenceFunctionLoad.h"
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  Next, we use \code{typedef}s to instantiate all necessary classes.
We
> > //  define the image and element types we plan to use to solve a
> > //  two-dimensional registration problem.  We define multiple element
> > //  types so that they can be used without recompiling the code.
> > //
> > //  Software Guide : EndLatex
> >
> >
> > //  Software Guide : BeginCodeSnippet
> > typedef itk::Image<unsigned char, 2>
fileImageType;
> > typedef itk::Image<float, 2>                               ImageType;
> > typedef itk::fem::Element2DC0LinearQuadrilateralMembrane   ElementType;
> > typedef itk::fem::Element2DC0LinearTriangularMembrane      ElementType2;
> > //  Software Guide : EndCodeSnippet
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  Note that in order to solve a three-dimensional registration
> > //  problem, we would simply define 3D image and element types in lieu
> > //  of those above.  The following declarations could be used for a 3D
> > //  problem:
> > //
> > //  SoftwareGuide : EndLatex
> >
> >
> > //  SoftwareGuide : BeginCodeSnippet
> > typedef itk::Image<unsigned char, 3>                    fileImage3DType;
> > typedef itk::Image<float, 3>                            Image3DType;
> > typedef itk::fem::Element3DC0LinearHexahedronMembrane   Element3DType;
> > typedef itk::fem::Element3DC0LinearTetrahedronMembrane  Element3DType2;
> > //  Software Guide : EndCodeSnippet
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  Here, we instantiate the load types and explicitly template the
> > //  load implementation type.  We also define visitors that allow the
> > //  elements and loads to communicate with one another.
> > //
> > //  Software Guide : EndLatex
> >
> >
> > //typedef itk::fem::ImageMetricLoad<ImageType,ImageType>
ImageLoadType;
> >
> > //  Software Guide : BeginCodeSnippet
> >
> > typedef itk::fem::FiniteDifferenceFunctionLoad<ImageType,ImageType>
ImageLoadType;
> > template class itk::fem::ImageMetricLoadImplementation<ImageLoadType>;
> >
> > typedef ElementType::LoadImplementationFunctionPointer     LoadImpFP;
> > typedef ElementType::LoadType
ElementLoadType;
> >
> > typedef ElementType2::LoadImplementationFunctionPointer    LoadImpFP2;
> > typedef ElementType2::LoadType
ElementLoadType2;
> >
> > typedef itk::fem::VisitorDispatcher<ElementType,ElementLoadType,
LoadImpFP>
> >
DispatcherType;
> >
> > typedef itk::fem::VisitorDispatcher<ElementType2,ElementLoadType2,
LoadImpFP2>
> >
DispatcherType2;
> > //  Software Guide : EndCodeSnippet
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  Once all the necessary components have been instantianted, we can
> > //  instantiate the \doxygen{FEMRegistrationFilter}, which depends on
the
> > //  image input and output types.
> > //
> > //  Software Guide : EndLatex
> >
> >
> > //  Software Guide : BeginCodeSnippet
> > typedef itk::fem::FEMRegistrationFilter<ImageType,ImageType>
RegistrationType;
> > //  Software Guide : EndCodeSnippet
> >
> >
> > int main(int argc, char *argv[])
> > {
> >   char *paramname;
> >   if ( argc < 2 )
> >     {
> >     std::cout << "Parameter file name missing" << std::endl;
> >     std::cout << "Usage: " << argv[0] << " param.file" << std::endl;
> >     return -1;
> >     }
> >   else
> >     {
> >     paramname=argv[1];
> >     }
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  The \doxygen{fem::ImageMetricLoad} must be registered before it
> > //  can be used correctly with a particular element type.  An example
> > //  of this is shown below for ElementType.  Similar
> > //  definitions are required for all other defined element types.
> > //
> > //  Software Guide : EndLatex
> >
> >   // Register the correct load implementation with the element-typed
visitor dispatcher.
> >   {
> > //  Software Guide : BeginCodeSnippet
> >   ElementType::LoadImplementationFunctionPointer fp =
> >
&itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetri
cLoad;
> >   DispatcherType::RegisterVisitor((ImageLoadType*)0,fp);
> > //  Software Guide : EndCodeSnippet
> >   }
> >   {
> >   ElementType2::LoadImplementationFunctionPointer fp =
> >
&itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetri
cLoad;
> >   DispatcherType2::RegisterVisitor((ImageLoadType*)0,fp);
> >   }
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  In order to begin the registration, we declare an instance of the
> > //  FEMRegistrationFilter.  For simplicity, we will call
> > //  it \code{registrationFilter}.
> > //
> > //  Software Guide : EndLatex
> >
> > //  Software Guide : BeginCodeSnippet
> >   RegistrationType::Pointer registrationFilter =
RegistrationType::New();
> > //  Software Guide : EndCodeSnippet
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  Next, we call \code{registrationFilter->SetConfigFileName()} to read
the parameter
> > //  file containing information we need to set up the registration
> > //  filter (image files, image sizes, etc.).  A sample parameter file is
shown at the end of this
> > //  section, and the individual components are labeled.
> > //
> > //  Software Guide : EndLatex
> >
> >
> >   // Attempt to read the parameter file, and exit if an error occurs
> >   registrationFilter->SetConfigFileName(paramname);
> >   if ( !registrationFilter->ReadConfigFile(
> >            (registrationFilter->GetConfigFileName()).c_str() ) )
> >     {
> >     return -1;
> >     }
> >
> >   // Read the image files
> >   typedef itk::ImageFileReader< fileImageType >      FileSourceType;
> >   typedef fileImageType::PixelType PixType;
> >
> >   FileSourceType::Pointer movingfilter = FileSourceType::New();
> >   movingfilter->SetFileName(
(registrationFilter->GetMovingFile()).c_str() );
> >   FileSourceType::Pointer fixedfilter = FileSourceType::New();
> >   fixedfilter->SetFileName(
(registrationFilter->GetFixedFile()).c_str() );
> >   std::cout << " reading moving " << registrationFilter->GetMovingFile()
<< std::endl;
> >   std::cout << " reading fixed " << registrationFilter->GetFixedFile()
<< std::endl;
> >
> >
> >   try
> >     {
> >     movingfilter->Update();
> >     }
> >   catch( itk::ExceptionObject & e )
> >     {
> >     std::cerr << "Exception caught during reference file reading " <<
std::endl;
> >     std::cerr << e << std::endl;
> >     return -1;
> >     }
> >   try
> >     {
> >     fixedfilter->Update();
> >     }
> >   catch( itk::ExceptionObject & e )
> >     {
> >     std::cerr << "Exception caught during target file reading " <<
std::endl;
> >     std::cerr << e << std::endl;
> >     return -1;
> >     }
> >
> >
> >   // Rescale the image intensities so that they fall between 0 and 255
> >   typedef itk::RescaleIntensityImageFilter<fileImageType,ImageType>
FilterType;
> >   FilterType::Pointer movingrescalefilter = FilterType::New();
> >   FilterType::Pointer fixedrescalefilter = FilterType::New();
> >
> >   movingrescalefilter->SetInput(movingfilter->GetOutput());
> >   fixedrescalefilter->SetInput(fixedfilter->GetOutput());
> >
> >   const double desiredMinimum =  0.0;
> >   const double desiredMaximum =  255.0;
> >
> >   movingrescalefilter->SetOutputMinimum( desiredMinimum );
> >   movingrescalefilter->SetOutputMaximum( desiredMaximum );
> >   movingrescalefilter->UpdateLargestPossibleRegion();
> >   fixedrescalefilter->SetOutputMinimum( desiredMinimum );
> >   fixedrescalefilter->SetOutputMaximum( desiredMaximum );
> >   fixedrescalefilter->UpdateLargestPossibleRegion();
> >
> >
> >   // Histogram match the images
> >   typedef itk::HistogramMatchingImageFilter<ImageType,ImageType>
HEFilterType;
> >   HEFilterType::Pointer IntensityEqualizeFilter = HEFilterType::New();
> >
> >   IntensityEqualizeFilter->SetReferenceImage(
fixedrescalefilter->GetOutput() );
> >   IntensityEqualizeFilter->SetInput( movingrescalefilter->GetOutput() );
> >   IntensityEqualizeFilter->SetNumberOfHistogramLevels( 100);
> >   IntensityEqualizeFilter->SetNumberOfMatchPoints( 15);
> >   IntensityEqualizeFilter->ThresholdAtMeanIntensityOn();
> >   IntensityEqualizeFilter->Update();
> >
> >   registrationFilter->SetFixedImage(fixedrescalefilter->GetOutput());
> >
registrationFilter->SetMovingImage(IntensityEqualizeFilter->GetOutput());
> >
> >
> >   itk::ImageFileWriter<ImageType>::Pointer writer;
> >   writer = itk::ImageFileWriter<ImageType>::New();
> >   std::string ofn="fixed.hdr";
> >   writer->SetFileName(ofn.c_str());
> >   writer->SetInput(registrationFilter->GetFixedImage() );
> >   writer->Write();
> >
> >   ofn="moving.hdr";
> >   itk::ImageFileWriter<ImageType>::Pointer writer2;
> >   writer2 =  itk::ImageFileWriter<ImageType>::New();
> >   writer2->SetFileName(ofn.c_str());
> >   writer2->SetInput(registrationFilter->GetMovingImage() );
> >   writer2->Write();
> >
> >
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  In order to initialize the mesh of elements, we must first create
> > //  ``dummy'' material and element objects and assign them to the
> > //  registration filter.  These objects are subsequently used to
> > //  either read a predefined mesh from a file or generate a mesh using
> > //  the software.  The values assigned to the fields within the
> > //  material object are arbitrary since they will be replaced with
> > //  those specified in the parameter file.  Similarly, the element
> > //  object will be replaced with those from the desired mesh.
> > //
> > //  Software Guide : EndLatex
> >
> > //  Software Guide : BeginCodeSnippet
> >   // Create the material properties
> >   itk::fem::MaterialLinearElasticity::Pointer m;
> >   m = itk::fem::MaterialLinearElasticity::New();
> >   m->GN = 0;                  // Global number of the material
> >   m->E = registrationFilter->GetElasticity();  // Young's modulus -- 
used in the membrane
> >   m->A = 1.0;                 // Cross-sectional area
> >   m->h = 1.0;                 // Thickness
> >   m->I = 1.0;                 // Moment of inertia
> >   m->nu = 0.;                 // Poisson's ratio -- DONT CHOOSE 1.0!!
> >   m->RhoC = 1.0;              // Density
> >
> >   // Create the element type
> >   ElementType::Pointer e1=ElementType::New();
> >   e1->m_mat=dynamic_cast<itk::fem::MaterialLinearElasticity*>( m );
> >   registrationFilter->SetElement(e1);
> >   registrationFilter->SetMaterial(m);
> > //  Software Guide : EndCodeSnippet
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  Now we are ready to run the registration:
> > //
> > //  Software Guide : EndLatex
> >
> > //  Software Guide : BeginCodeSnippet
> >   registrationFilter->RunRegistration();
> > //  Software Guide : EndCodeSnippet
> >
> >
> > //  Software Guide : BeginLatex
> > //
> > //  To output the image resulting from the registration, we can call
> > //  \code{WriteWarpedImage()}.  The image is written in floating point
> > //  format.
> > //
> > //  Software Guide : EndLatex
> >
> > //  Software Guide : BeginCodeSnippet
> >   registrationFilter->WriteWarpedImage(
> >         (registrationFilter->GetResultsFileName()).c_str());
> > //  Software Guide : EndCodeSnippet
> >
> > //  Software Guide : BeginLatex
> > //
> > //  We can also output the displacement fields resulting from the
> > //  registration, we can call \code{WriteDisplacementField()} with the
> > //  desired vector component as an argument.  For a $2D$ registration,
> > //  you would want to write out both the $x$ and $y$ displacements, and
> > //  this requires two calls to the aforementioned function.
> > //
> > //  Software Guide : EndLatex
> >
> > //  Software Guide : BeginCodeSnippet
> >   if (registrationFilter->GetWriteDisplacements())
> >     {
> >     registrationFilter->WriteDisplacementField(0);
> >     registrationFilter->WriteDisplacementField(1);
> >     // If this were a 3D example, you might also want to call this line:
> >     // registrationFilter->WriteDisplacementField(2);
> >
> >     // We can also write it as a multicomponent vector field
> >     registrationFilter->WriteDisplacementFieldMultiComponent();
> >     }
> > //  Software Guide : EndCodeSnippet
> >
> >   //  This is a documented sample parameter file that can be used with
> >   //  this deformable registration example.
> >   //
> >   //  ../Data/FiniteElementRegistrationParameters1.txt
> >   //
> >
> >   return 0;
> > }
> >
> >
> >
> >
> > ------------------------------------------------------------------------
> >
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users at itk.org
> > http://www.itk.org/mailman/listinfo/insight-users
>
>
>
>
>



More information about the Insight-users mailing list