[Insight-users] binary deformable registration : THERE IS NO
SCIENCE WITHOUT REPRODUCIBILITY
Luis Ibanez
luis.ibanez at kitware.com
Sun Sep 17 20:38:05 EDT 2006
Hi Siddharth
Thanks a lot for posting your images, source code and parameters.
This information now makes possible to track the problem that you
are reporting.
We used a modified version of the source code in
Insight/Examples/Registration/
DeformableRegistration1.cxx
The code was modified to work in 3D images.
After playing with your images and the parameters file,
We got to a relatively acceptable value of Elasticity
would be in the range of 6e4.
You probably already know this, but just in case, the values
of the Elasticity modulus are a bit counterintuitive. The larger
values correspond to the harder materials.
http://en.wikipedia.org/wiki/Young%27s_modulus#Approximate_values
The way to formally find the desired value of the modulus to use
is to consider that the forces are computed as derivative of the
image metrics, and the deformations are measured in the units of
spacing.
The derivatives of image metrics will have virtual units such as
ImageIntensity / mm
while the deformations will have units of mm (assuming that you
use an image modality such as CT and MRI.
In your binary images, the kind of Modulus that you want is the
one that will produce a displacement of 1 pixel when you shift
the intensities of your images by 1 pixel. Therefore, you want
it to be in the order of "255". Which is the variation of a
MeanSquares metric when you shift an element by one pixel.
Note however that this may result in a too soft and mushy material,
that will easily shift even in conditions where it shouldn't.
Therefore you want to increase the elasticity (hardness) of the
material, as a mechanism for regularizing the deformation field.
That is, the harder the material, the more uniform the deformation
field will be.
Please find attached the source code and parameters file that we
found to be appropriate for registering your two binary images.
Please run it and let us know what you find.
Thanks
Luis
=========================
Siddharth Vikal wrote:
> Thanks Luis for detailed email.
>
> Yes I agree that there is no science without reproducibility, and that
> my statements were incomplete without specifiying parameters and images
> that I tried upon.
>
> Find in this email the two binary files to be registered, FEM parameters
> file as attachments, and code snippet inline.
>
> But thanks for clarifying about FEM toolkit that choosing Young's
> modulus and elasticity values from tissue properties is not going to be
> useful and it is images dependent.
>
> Pardon my lack of understanding from scanty documentation on FEM
> registration that i found and from looking at .h file, I'm still not
> clear about how to set up the registration parameters for my problem.
> You say:
>
> the magnitude of this forces are *NOT* measured in Newtons, or any other
> physically consistent unit, but in the units of the Image Metric that
> you are using divided by units of pixel spacing; The values that you
> should explore for the Young modulus and Poisson ratio are based on the
> magnitude of the image metric derivatives.
>
> So say if my image metric derivative is x, how will i infer/explore
> suitable Young's modulus?
> Perhaps you can throw some more light on the relationship more specifically.
>
> And after looking at images, you can provide me some advice for setting
> up parameters for my problem.
>
>
> The files attached to be deformable registered have already been rigid
> registered and made isotropic.
>
> Looking forward to hear from you.
>
> regards
> siddharth
>
>
>
===============================================
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: DeformableRegistration1.cxx,v $
Language: C++
Date: $Date: 2005/11/19 16:31:50 $
Version: $Revision: 1.32 $
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.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkRescaleIntensityImageFilter.h"
#include "itkHistogramMatchingImageFilter.h"
#include "itkFEM.h"
#include "itkFEMRegistrationFilter.h"
/* Example of FEM-base deformable registration in 3D */
const unsigned int Dimension = 3;
typedef itk::Image<unsigned char, Dimension> FileImageType;
typedef itk::Image<float, Dimension> ImageType;
typedef itk::fem::Element3DC0LinearHexahedronMembrane ElementType;
typedef itk::fem::Element3DC0LinearTetrahedronMembrane ElementType2;
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;
typedef itk::fem::FEMRegistrationFilter<ImageType,ImageType> RegistrationType;
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];
}
// Register the correct load implementation with the element-typed visitor
// dispatcher.
typedef itk::fem::ImageMetricLoadImplementation<
ImageLoadType> LoadImplementationType;
{
ElementType::LoadImplementationFunctionPointer fp =
&LoadImplementationType::ImplementImageMetricLoad;
DispatcherType::RegisterVisitor((ImageLoadType*)0,fp);
}
{
ElementType2::LoadImplementationFunctionPointer fp =
&LoadImplementationType::ImplementImageMetricLoad;
DispatcherType2::RegisterVisitor((ImageLoadType*)0,fp);
}
RegistrationType::Pointer registrationFilter = RegistrationType::New();
// 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 ";
std::cout << registrationFilter->GetMovingFile() << std::endl;
std::cout << " reading fixed ";
std::cout << registrationFilter->GetFixedFile() << std::endl;
try
{
movingfilter->Update();
}
catch( itk::ExceptionObject & e )
{
std::cerr << "Exception caught during reference file reading ";
std::cerr << std::endl << e << std::endl;
return -1;
}
try
{
fixedfilter->Update();
}
catch( itk::ExceptionObject & e )
{
std::cerr << "Exception caught during target file reading ";
std::cerr << std::endl << 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();
writer->SetFileName("fixed.mhd");
writer->SetInput(registrationFilter->GetFixedImage() );
writer->Write();
itk::ImageFileWriter<ImageType>::Pointer writer2;
writer2 = itk::ImageFileWriter<ImageType>::New();
writer2->SetFileName("moving.mhd");
writer2->SetInput(registrationFilter->GetMovingImage() );
writer2->Write();
// 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
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);
registrationFilter->RunRegistration();
registrationFilter->WarpImage( registrationFilter->GetMovingImage() );
// Warp the moving image and write it to a file.
writer->SetFileName("warpedMovingImage.mhd");
writer->SetInput( registrationFilter->GetWarpedImage() );
writer->Update();
registrationFilter->WriteDisplacementFieldMultiComponent();
// This is a documented sample parameter file that can be used with
// this deformable registration example.
//
// ../Data/FiniteElementRegistrationParameters3.txt
//
return 0;
}
-------------- next part --------------
Reading config file...FEM_RegistrationParameters.txt
Example configured. E 6000 rho 3000
reading moving ResampledVolume_RigidRegistered.mha
reading fixed ResampledVolume_Intraop.mha
beginning level 0
scaling 1
scaling 1
scaling 1
ElementsPerDim 27.25 27.25 14.5
generating regular mesh
generating regular mesh done
DO NOT init interpolation grid : im sz 110 110 59 MeshSize 109 109 58
applying loads
which node 1
edge coord 4.03704 0 0
which node 1
edge coord 4.03704 0 0
which node 1
edge coord 4.03704 0 0
which node 3
edge coord 0 4.03704 0
which node 3
edge coord 0 4.03704 0
which node 3
edge coord 0 4.03704 0
which node 4
edge coord 0 0 4.14286
which node 4
edge coord 0 0 4.14286
which node 4
edge coord 0 0 4.14286
node 1
allocating deformation field
load sizes [109, 109, 58] image [109, 109, 58]
load sizes [109, 109, 58] image [109, 109, 58]
energy 3.3748e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 0 delt E -3.3748e+07 iter 1
load sizes [109, 109, 58] image [109, 109, 58]
energy 4.8111e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 3.3748e+07 delt E -1.43631e+07 iter 2
load sizes [109, 109, 58] image [109, 109, 58]
energy 4.47308e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 4.8111e+07 delt E 3.38028e+06 iter 3
load sizes [109, 109, 58] image [109, 109, 58]
energy 4.10104e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 4.47308e+07 delt E 3.7204e+06 iter 4
load sizes [109, 109, 58] image [109, 109, 58]
energy 3.80697e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 4.10104e+07 delt E 2.94068e+06 iter 5
load sizes [109, 109, 58] image [109, 109, 58]
energy 3.53508e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 3.80697e+07 delt E 2.71886e+06 iter 6
load sizes [109, 109, 58] image [109, 109, 58]
energy 3.29253e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 3.53508e+07 delt E 2.42557e+06 iter 7
load sizes [109, 109, 58] image [109, 109, 58]
energy 3.05455e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 3.29253e+07 delt E 2.37979e+06 iter 8
load sizes [109, 109, 58] image [109, 109, 58]
energy 2.78391e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 3.05455e+07 delt E 2.70641e+06 iter 9
load sizes [109, 109, 58] image [109, 109, 58]
energy 2.56444e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 2.78391e+07 delt E 2.19472e+06 iter 10
load sizes [109, 109, 58] image [109, 109, 58]
energy 2.38429e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 2.56444e+07 delt E 1.80141e+06 iter 11
load sizes [109, 109, 58] image [109, 109, 58]
energy 2.22448e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 2.38429e+07 delt E 1.59812e+06 iter 12
load sizes [109, 109, 58] image [109, 109, 58]
energy 2.08255e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 2.22448e+07 delt E 1.4193e+06 iter 13
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.9299e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 2.08255e+07 delt E 1.52647e+06 iter 14
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.78411e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.9299e+07 delt E 1.45797e+06 iter 15
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.65446e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.78411e+07 delt E 1.29643e+06 iter 16
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.55862e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.65446e+07 delt E 958488 iter 17
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.48775e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.55862e+07 delt E 708685 iter 18
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.43914e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.48775e+07 delt E 486120 iter 19
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.39865e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.43914e+07 delt E 404843 iter 20
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.36418e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.39865e+07 delt E 344714 iter 21
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.33111e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.36418e+07 delt E 330711 iter 22
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.2984e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.33111e+07 delt E 327115 iter 23
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.2747e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.2984e+07 delt E 236955 iter 24
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.25828e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.2747e+07 delt E 164223 iter 25
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.23904e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.25828e+07 delt E 192366 iter 26
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.21319e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.23904e+07 delt E 258563 iter 27
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.18503e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.21319e+07 delt E 281567 iter 28
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.16941e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.18503e+07 delt E 156169 iter 29
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.16216e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.16941e+07 delt E 72541.9 iter 30
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.15762e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.16216e+07 delt E 45340 iter 31
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.15781e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.15762e+07 delt E -1827.12 iter 32
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.16178e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.15781e+07 delt E -39709.5 iter 33
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.16829e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.16178e+07 delt E -65162.4 iter 34
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.17962e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.16829e+07 delt E -113207 iter 35
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.19609e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.17962e+07 delt E -164757 iter 36
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.21704e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.19609e+07 delt E -209509 iter 37
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.23546e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.21704e+07 delt E -184190 iter 38
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.24554e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.23546e+07 delt E -100827 iter 39
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.25305e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.24554e+07 delt E -75102.3 iter 40
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.26003e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.25305e+07 delt E -69775.6 iter 41
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.26793e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.26003e+07 delt E -79034.3 iter 42
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.27771e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.26793e+07 delt E -97764.5 iter 43
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.28793e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.27771e+07 delt E -102180 iter 44
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.30376e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.28793e+07 delt E -158351 iter 45
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.32156e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.30376e+07 delt E -177949 iter 46
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.33962e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.32156e+07 delt E -180622 iter 47
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.35825e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.33962e+07 delt E -186281 iter 48
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.37654e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.35825e+07 delt E -182886 iter 49
load sizes [109, 109, 58] image [109, 109, 58]
energy 1.39565e+07
interpolating vector field of size [109, 109, 58] interpolation done
min E 1.37654e+07 delt E -191086 iter 50
field pix [7.47107e-08, -5.71859e-07, -2.86566e-06]
max vec 38.4134
field size [109, 109, 58]
end level 0 allocating m_FloatImage
get jacobian
min Jacobian 0
Warping image
Warping image
Writing multi-component displacement vector field...done
-------------- next part --------------
% ---------------------------------------------------------
% Parameters for the single- or multi-resolution techniques
% ---------------------------------------------------------
1 % Number of levels in the multi-res pyramid (1 = single-res)
1 % Highest level to use in the pyramid
1 1 1 % Scaling at lowest level of pyramid
4 % Number of pixels per element
6.e3 % Elasticity (E)
3.e3 % Density x capacity (RhoC)
1 % Image energy scaling (gamma) - sets gradient step size
2 % NumberOfIntegrationPoints
1 % WidthOfMetricRegion
50 % MaximumIterations
% -------------------------------
% Parameters for the registration
% -------------------------------
0 0.99 % Similarity metric (0=mean sq, 1 = ncc, 2=pattern int, 3=MI, 5=demons)
1.0 % Alpha
0 % DescentDirection (1 = max, 0 = min)
0 % DoLineSearch (0=never, 1=always, 2=if needed)
1.e1 % TimeStep
0.5 % Landmark variance
0 % Employ regridding / enforce diffeomorphism ( >= 1 -> true)
% ----------------------------------
% Information about the image inputs
% ----------------------------------
109 % Nx (image x dimension)
109 % Ny (image y dimension)
58 % Nz (image z dimension - not used if 2D)
ResampledVolume_RigidRegistered.mha % ReferenceFileName
ResampledVolume_Intraop.mha % TargetFileName
% -------------------------------------------------------------------
% The actions below depend on the values of the flags preceding them.
% For example, to write out the displacement fields, you have to set
% the value of WriteDisplacementField to 1.
% -------------------------------------------------------------------
0 % UseLandmarks? - read the file name below if this is true
FiniteElementRegistrationLandmarks.txt % LandmarkFileName
DeformableRegistered.mha % ResultsFileName (prefix only)
0 % WriteDisplacementField?
./RatLung_disp % DisplacementsFileName (prefix only)
0 % ReadMeshFile?
- % MeshFileName
END
More information about the Insight-users
mailing list