[Insight-users] All zeros displacement field when
using DeformableRegistration1
Luis Ibanez
luis.ibanez at kitware.com
Tue Nov 9 17:07:45 EST 2004
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>::ImplementImageMetricLoad;
> DispatcherType::RegisterVisitor((ImageLoadType*)0,fp);
> // Software Guide : EndCodeSnippet
> }
> {
> ElementType2::LoadImplementationFunctionPointer fp =
> &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
> 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