[Insight-users] All zeros displacement field when using
DeformableRegistration1
Barbara
bgarita at ufl.edu
Tue Nov 9 14:09:38 EST 2004
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
>
>
>
>
>
>
>
-------------- next part --------------
/*=========================================================================
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;
}
More information about the Insight-users
mailing list