[Insight-users] All zeros displacement field when using
DeformableRegistration1
Luis Ibanez
luis.ibanez at kitware.com
Wed Nov 10 17:44:57 EST 2004
Hi Barbara,
About your questions:
1) When you compute a deformation field in Image Registration,
the pixels in the deformation field are vectors that indicate
how you should move in physical space in order to go from one
pixel in one of the images to the corresponding matching pixel
in the other image.
The magnitude of the vectors in the deformation field should
therefore be interpreted in physical units of length
(e.g. millimeters).
The reason for my question is that if you end up with a deformation
field whose pixel vectors have values in the range of the thousands,
while you are registering images whose extent is in the range of
the hundreds, then that will indicate that the deformation field
if displacing the image *so far* that there is no overlap between
the mapped image and the reference image. This will result in a
resampled image that is full of zero pixels. That's the situation
I was referring to as "extreme deformation".
After your response... that doesn't seem to be the case of your data.
since you have pixel spacings of 0.32 millimeters, and displacement
vectors in the range of 10, so most pixels are being mapped to
positions that are 3 pixels away from its current site.
2) An Histogram count is *suspicious* given that you are dealing with
float images. Computing an histogram of a float image is not as
trivial as it may sound. You have to be careful with the setup of
the minimum and maximum values of the histogram, the number of bins
and so on.
It is possible for example that your image is plenty of information
in the range 0.1 to 0.9, and that you are attempting to compute an
histogram for the classical 8 bit range of 0-255. In that case all
the values of the image will appear as Zeros, despite the fact that
they contain useful and valid information.
Please let us know how did you computed this histogram.
You could double check those values by connecting the MaximumMinimum
image calculator to your output image and printing the result of the
Maximum and Minimum pixel values.
Something that you may want to try rapidly is to save the "zero" image
using a MetaImage file format, and load it in ParaView, where you will
be able to see the entire range of the intensities form the suspect
image.
ParaView is an open source and free application that can be downloaded
from
http://www.paraview.org
you will find binary build for most major platforms.
Please let us know what you find,
Thanks
Luis
------------------------
Barbara Garita wrote:
> 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