[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