[Insight-users] deformable registration in binary image
Luis Ibanez
luis.ibanez at kitware.com
Fri Jan 20 10:59:52 EST 2006
Hi Marc,
It seems that the metric is being evaluated with such a small
change in the parameters of the deformable transform, that the
values Metric gradient appears to be zero.
One of the important characteristics of image metrics on binary
image is that they do not change progressively.
That is, a small change in the transform parameters, may not
produce a small change in the metric value. Simply because the
fixed image pixels may still map on top of the same pixels of the
moving image.
You can easily verify this by adding a print out in line 289 of
the file
Insight/
Code/
Numerics/
itkLBFGSBOptimizer.cxx
with code like:
std::cout << "Gradient = " << gradient << " std::endl;
If our assumption is correct, then you should see here an
array of zeros in the print out.
If the problem is the null gradient, then you may want
to try one of the two following options:
1) Blur your images with a Smoothing filter such as a Gaussian,
and then use a normal Image metric suited for registering
grayscale images.
OR
2) Keep the images as binary images, and increase the mobility
of the transform, in order to make that the metric actually
change values, and provide you with non null derivatives.
Please let us know what you observe.
Thanks
Luis
--------------------
Marc Ruiz wrote:
> Hi!
>
> We are stopped and we do not how to continue... Let me explain with
> precision the problem.
>
> We are registering 2 segmented images. Our goal is registrate them the
> most closely as it was possible. Then we would try with rigid, affine
> and deformable registration. Rigid and affine work well but deformable
> does not. This is the code of deformable registration. You will see it
> is deformable registration 7 or 8 changing the metric
> (matchCardinality) and the interpolator (NN).
>
> The problem is that deformable registration stops without doing any
> iteration. It stops because of these values:
> GetValue()= 0.0709776 and GetInfinityNormOfProjectedGradient()=0
>
> Luis suggests that a possible solution was the method
> SetBulkTransform(), because the images are close one from each other,
> but I have tried with to images not registered (so far from each
> other..) at I get the same problem.
>
> So, is there any mistake in the code?? Or any advice??
>
> I would be very grateful of your indications. Thank you!
>
> On 1/16/06, *Luis Ibanez* <luis.ibanez at kitware.com
> <mailto:luis.ibanez at kitware.com>> wrote:
>
>
> Hi Marc,
>
>
> I'm not sure I understand your description.
>
> I would assume that if you initialize the transform
> properly then the images would be close to each other...
>
>
>
> In any case....
>
>
> You seem to be going in the right direction, since you are first
> solving the Rigid registration, then the Affine registration and
> finally moving into the Deformable registration.
>
> Please not that if you are using the BSplineDeformableTransform
> http://www.itk.org/Insight/Doxygen/html/d7/d3b/classitk_1_1BSplineDeformableTransform.html
>
> this Transform let you define a "Bulk" Transform. If you
> provide this bulk transform, the BsplineDeformable transform
> will use it as a pre-composed transform.
>
> (see the SetBulkTransform() method...)
>
>
> What it is usually done, is to take the Affine transform
> resulting from your initial registration, and plug it in
> as the bulk transform for the deformable registration.
>
> In this way, you can initialize the deformable transform
> as an identity transform (e.g. no deformations initially),
> and have the bulk transform represent the compensation for
> the general misalignment of the two images.
>
>
>
> Please give it a try,
> and let us know if you find any problems.
>
>
>
> Thanks
>
>
>
> Luis
>
>
>
> ----------------
> Marc Ruiz wrote:
> > Hello!
> >
> > The 2 images are first rigid and then affine registered. If they are
> > initialized they are far one from the other.
> >
> > So... I have to launch the registration without doing these two
> steps??
> > But then, I read deformable registration is not capable of doing big
> > translations or rotations.
> >
> > Any suggest??
> >
> >
> > On 1/16/06, *Luis Ibanez* < luis.ibanez at kitware.com
> <mailto:luis.ibanez at kitware.com>
> > <mailto:luis.ibanez at kitware.com
> <mailto:luis.ibanez at kitware.com>>> wrote:
> >
> >
> > Hi Roger,
> >
> >
> > Your modifications to the code seem to be reasonable.
> >
> > The main difference between the examples
> >
> >
> > DeformableRegistration7
> >
> > and
> >
> > DeformableRegistration8
> >
> >
> > is that the first one is intended for images of the
> > same modality, while the second is intended for images
> > of different modality.
> >
> > Given that you replaced the Metric with the MatchCardinality
> > metric, then there is no major difference left between the
> > two examples.
> >
> >
> > I'm not sure about the problem you find with the
> > infinity norm of the projected gradient,.. but if
> > you are using the LBFGSB optimizer, it is known
> > that it has some issues if you start it on a local
> > minimum (maximum), in other words, if you start the
> > optimization problem in a place where the derivative
> > of the metric is null, then the optimizer can not
> > complete the first iteration.
> >
> >
> > Please give us more details on how you are initializing
> > the Transform, and how similar the Fixed and Moving
> > images are.
> >
> >
> > Thanks
> >
> >
> >
> > Luis
> >
> >
> >
> >
> >
> >
> > --------------------
> > Roger Alvaredo wrote:
> > > Hi!
> > >
> > > I tested deformableRegistration7.cxx or
> > deformableRegistration8.cxx in
> > > MRI's and they work correctly. Now I would like to re
> adapt def8 for
> > > working with binary images.
> > >
> > > I change the metric to MatcthCardinality and the
> interpolator to
> > > NearestNeighbor.
> > >
> > > In the case of the metric I add the
> > method: metric->MeasureMatchesOff()
> > > because the optimizer minimizes the metric, and I remove
> the lines of
> > > MattesMutualInformation as the number of histograms and so
> on. It
> > is ok??
> > >
> > > So, I would like to know...:
> > >
> > > 1. Advantatges and disvantatges of each one (def7 and def8).
> > > 2. I have a problem because the change of metric...:
> > > The initial values are: GetValue()= 0.0709776
> > > and GetInfinityNormOfProjectedGradient()=0
> > > why is 0 the norm? it stops before doing an iteration for
> this...
> > >
> > > Any advice?
> > >
> > > Thanks!
> > >
> > >
> > >
> > >
> > >
> >
> ------------------------------------------------------------------------
> > >
> > > _______________________________________________
> > > Insight-users mailing list
> > > Insight-users at itk.org <mailto:Insight-users at itk.org>
> <mailto:Insight-users at itk.org <mailto:Insight-users at itk.org>>
> > > http://www.itk.org/mailman/listinfo/insight-users
> > <http://www.itk.org/mailman/listinfo/insight-users
> <http://www.itk.org/mailman/listinfo/insight-users>>
> >
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users at itk.org <mailto:Insight-users at itk.org> <mailto:
> Insight-users at itk.org <mailto:Insight-users at itk.org>>
> > http://www.itk.org/mailman/listinfo/insight-users
> >
> >
> >
> >
> > --
> >
> > MaRC//
>
>
>
>
> --
>
> MaRC//
>
>
> ------------------------------------------------------------------------
>
> /*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: MaskDeformableRegistration.cxx,v $
> Language: C++
> =========================================================================*/
>
>
> #include "itkImageRegistrationMethod.h"
> #include "itkMatchCardinalityImageToImageMetric.h"
> #include "itkNearestNeighborInterpolateImageFunction.h"
> #include "itkImage.h"
>
> #include "itkTimeProbesCollectorBase.h"
>
> #include "itkBSplineDeformableTransform.h"
> #include "itkLBFGSBOptimizer.h"
>
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
>
> #include "itkResampleImageFilter.h"
> #include "itkCastImageFilter.h"
>
> #include "itkCommand.h"
> class CommandIterationUpdate : public itk::Command
> {
> public:
> typedef CommandIterationUpdate Self;
> typedef itk::Command Superclass;
> typedef itk::SmartPointer<Self> Pointer;
> itkNewMacro( Self );
> protected:
> CommandIterationUpdate() {};
> public:
> typedef itk::LBFGSBOptimizer OptimizerType;
> typedef const OptimizerType * OptimizerPointer;
>
> void Execute(itk::Object *caller, const itk::EventObject & event)
> {
> Execute( (const itk::Object *)caller, event);
> }
>
> void Execute(const itk::Object * object, const itk::EventObject & event)
> {
> OptimizerPointer optimizer =
> dynamic_cast< OptimizerPointer >( object );
> if( typeid( event ) != typeid( itk::IterationEvent ) )
> {
> return;
> }
> std::cout << optimizer->GetCurrentIteration() << " ";
> std::cout << optimizer->GetValue() << " ";
> std::cout << optimizer->GetInfinityNormOfProjectedGradient() << std::endl;
> }
> };
>
>
> int main( int argc, char *argv[] )
> {
> if( argc < 4 )
> {
> std::cerr << "Missing Parameters " << std::endl;
> std::cerr << "Usage: " << argv[0];
> std::cerr << " fixedImageFile movingImageFile outputImagefile ";
> return 1;
> }
>
> const unsigned int ImageDimension = 3;
> typedef signed short PixelType;
>
> typedef itk::Image< PixelType, ImageDimension > FixedImageType;
> typedef itk::Image< PixelType, ImageDimension > MovingImageType;
>
>
> const unsigned int SpaceDimension = ImageDimension;
> const unsigned int SplineOrder = 3;
> typedef double CoordinateRepType;
>
> typedef itk::BSplineDeformableTransform<
> CoordinateRepType,
> SpaceDimension,
> SplineOrder > TransformType;
>
>
> typedef itk::LBFGSBOptimizer OptimizerType;
>
> typedef itk::MatchCardinalityImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;
>
> typedef itk:: NearestNeighborInterpolateImageFunction<
> MovingImageType,
> double > InterpolatorType;
>
> typedef itk::ImageRegistrationMethod<
> FixedImageType,
> MovingImageType > RegistrationType;
>
> MetricType::Pointer metric = MetricType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> RegistrationType::Pointer registration = RegistrationType::New();
>
>
> registration->SetMetric( metric );
> registration->SetOptimizer( optimizer );
> registration->SetInterpolator( interpolator );
>
>
> TransformType::Pointer transform = TransformType::New();
> registration->SetTransform( transform );
>
> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
>
> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
>
> fixedImageReader->SetFileName( argv[1] );
> movingImageReader->SetFileName( argv[2] );
>
> FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
>
> registration->SetFixedImage( fixedImage );
> registration->SetMovingImage( movingImageReader->GetOutput() );
>
> fixedImageReader->Update();
>
> FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
>
> registration->SetFixedImageRegion( fixedRegion );
>
>
> typedef TransformType::RegionType RegionType;
> RegionType bsplineRegion;
> RegionType::SizeType gridSizeOnImage;
> RegionType::SizeType gridBorderSize;
> RegionType::SizeType totalGridSize;
>
> gridSizeOnImage.Fill( 12 );
> gridSizeOnImage[2] = 3;
> gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper )
> totalGridSize = gridSizeOnImage + gridBorderSize;
>
> bsplineRegion.SetSize( totalGridSize );
>
> typedef TransformType::SpacingType SpacingType;
> SpacingType spacing = fixedImage->GetSpacing();
>
> typedef TransformType::OriginType OriginType;
> OriginType origin = fixedImage->GetOrigin();;
>
> FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();
>
> for(unsigned int r=0; r<ImageDimension; r++)
> {
> spacing[r] *= floor( static_cast<double>(fixedImageSize[r] - 1) /
> static_cast<double>(gridSizeOnImage[r] - 1) );
> origin[r] -= spacing[r];
> }
>
> transform->SetGridSpacing( spacing );
> transform->SetGridOrigin( origin );
> transform->SetGridRegion( bsplineRegion );
>
>
> typedef TransformType::ParametersType ParametersType;
>
> const unsigned int numberOfParameters =
> transform->GetNumberOfParameters();
>
> ParametersType parameters( numberOfParameters );
>
> parameters.Fill( 0.0 );
>
> transform->SetParameters( parameters );
>
>
> registration->SetInitialTransformParameters( transform->GetParameters() );
>
> metric->MeasureMatchesOff();
>
> OptimizerType::BoundSelectionType boundSelect( transform->GetNumberOfParameters() );
> OptimizerType::BoundValueType upperBound( transform->GetNumberOfParameters() );
> OptimizerType::BoundValueType lowerBound( transform->GetNumberOfParameters() );
>
> boundSelect.Fill( 0 );
> upperBound.Fill( 0.0 );
> lowerBound.Fill( 0.0 );
>
> optimizer->SetBoundSelection( boundSelect );
> optimizer->SetUpperBound( upperBound );
> optimizer->SetLowerBound( lowerBound );
>
> optimizer->SetCostFunctionConvergenceFactor( 1e1 );
> optimizer->SetProjectedGradientTolerance( 3.5e-7);
> optimizer->SetMaximumNumberOfIterations( 500 );
> optimizer->SetMaximumNumberOfEvaluations( 500 );
> optimizer->SetMaximumNumberOfCorrections( 12 );
>
>
> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
> optimizer->AddObserver( itk::IterationEvent(), observer );
>
>
> // Add a time probe
> itk::TimeProbesCollectorBase collector;
>
> std::cout << std::endl << "Starting Registration" << std::endl;
>
> try
> {
> collector.Start( "Registration" );
> registration->StartRegistration();
> collector.Stop( "Registration" );
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
>
> OptimizerType::ParametersType finalParameters =
> registration->GetLastTransformParameters();
>
> // std::cout << "Last Transform Parameters" << std::endl;
> // std::cout << finalParameters << std::endl;
>
>
> // Report the time taken by the registration
> collector.Report();
>
>
> transform->SetParameters( finalParameters );
>
>
> std::cout << optimizer->GetCurrentIteration() << " ";
> std::cout << optimizer->GetValue() << " ";
> std::cout << optimizer->GetInfinityNormOfProjectedGradient() << std::endl;
>
>
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType > ResampleFilterType;
>
> ResampleFilterType::Pointer resample = ResampleFilterType::New();
>
> resample->SetTransform( transform );
> resample->SetInput( movingImageReader->GetOutput() );
>
> resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
> resample->SetOutputOrigin( fixedImage->GetOrigin() );
> resample->SetOutputSpacing( fixedImage->GetSpacing() );
> resample->SetDefaultPixelValue( 100 );
>
> typedef signed short OutputPixelType;
>
> typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;
>
> typedef itk::CastImageFilter<
> FixedImageType,
> OutputImageType > CastFilterType;
>
> typedef itk::ImageFileWriter< OutputImageType > WriterType;
>
>
> WriterType::Pointer writer = WriterType::New();
> CastFilterType::Pointer caster = CastFilterType::New();
>
>
> writer->SetFileName( argv[3] );
>
>
> caster->SetInput( resample->GetOutput() );
> writer->SetInput( caster->GetOutput() );
>
>
> try
> {
> writer->Update();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
>
>
> return 0;
> }
>
More information about the Insight-users
mailing list