[Insight-users] Re: some questions on BSpline Registration
Luis Ibanez
luis.ibanez at kitware.com
Tue Jul 20 00:30:18 EDT 2004
Hi He Smth,
YOU DO NOT NEED TO CHANGE THE PIXEL SPACING
OR THE IMAGE ORIGIN !!!
YOU SHOULD NOT MODIFY THE PIXEL SPACING OR
THE IMAGE ORIGIN OF AN IMAGE THAT YOU READ
FROM A DICOM FILE.
PLEASE REMOVE THE LINES
spacing.Fill(1);
fixedImage->SetSpacing(spacing);
movingImage->SetSpacing(spacing);
THESE LINES ARE WRONG AND ARE DANGEROUS
..........
ITK read the origin and spacing of your
image directly from the DICOM files.
YOU SHOULD NOT CHANGE IT !
The Example DeformableRegistration4 is using
the origin and spacing of the Input image in
order to set to origin and spacing of the GRID
that support the BSpline. This is done in order
to ensure that the GRID has one node past the
image in the lower range of the coordinates, and
two nodes past the image in the higher range of
the coordinates.
Regards,
Luis
----------------
he smth wrote:
> hi Luis and all,
> this is the code of BSplineRegistration using DICOM files
> I has set all the argv parameters for debug
> I am still confused whith the origin and spacing
> thanks for any help
>
> btw, it seems not support long file name. so I set fixedImage,
> movingImage as Image1 and Image2.
>
>
>
> /*=========================================================================
>
> Program: Insight Segmentation & Registration Toolkit
> Module: $RCSfile: DeformableRegistration4.cxx,v $
> Language: C++
> Date: $Date: 2004/05/22 01:11:13 $
> Version: $Revision: 1.7 $
>
> 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.
>
> =========================================================================*/
>
> // Software Guide : BeginLatex
> //
> // This example illustrates the use of the \doxygen{BSplineDeformableTransform}
> // class for performing registration of two $2D$ images. The example code is
> // for the most part identical to the code presented in
> // Section~\ref{sec:RigidRegistrationIn2D}. The major difference is that this
> // example we replace the Transform for a more generic one endowed with a large
> // number of degrees of freedom. Due to the large number of parameters, we will
> // also replace the simple steepest descent optimizer with the
> // \doxygen{LBFGSOptimizer}.
> //
> //
> // \index{itk::BSplineDeformableTransform}
> // \index{itk::BSplineDeformableTransform!DeformableRegistration}
> // \index{itk::LBFGSOptimizer}
> //
> //
> // Software Guide : EndLatex
>
>
> #include "itkImageRegistrationMethod.h"
> #include "itkMeanSquaresImageToImageMetric.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkImage.h"
> #include "itkDICOMImageIO2Factory.h"
>
> #include "itkTimeProbesCollectorBase.h"
>
> // Software Guide : BeginLatex
> //
> // The following are the most relevant headers to this example.
> //
> // \index{itk::BSplineDeformableTransform!header}
> // \index{itk::LBFGSOptimizer!header}
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> #include "itkBSplineDeformableTransform.h"
> #include "itkLBFGSOptimizer.h"
> // Software Guide : EndCodeSnippet
>
> // Software Guide : BeginLatex
> //
> // The parameter space of the \code{BSplineDeformableTransform} is composed by
> // the set of all the deformations associated with the nodes of the BSpline
> // grid. This large number of parameters makes possible to represent a wide
> // variety of deformations, but it also has the price of requiring a
> // significant amount of computation time.
> //
> // \index{itk::BSplineDeformableTransform!header}
> //
> // Software Guide : EndLatex
>
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
>
> #include "itkResampleImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkSquaredDifferenceImageFilter.h"
>
>
> // NOTE: the LBFGSOptimizer does not invoke events
>
>
> 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 ";
> std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";
> std::cerr << " [deformationField] ";
> return 1;
> }
> */
> const unsigned int ImageDimension = 2;
> typedef short PixelType;
>
> typedef itk::Image< PixelType, ImageDimension > FixedImageType;
> typedef itk::Image< PixelType, ImageDimension > MovingImageType;
>
>
> // Software Guide : BeginLatex
> //
> // We instantiate now the type of the \code{BSplineDeformableTransform} using
> // as template parameters the type for coordinates representation, the
> // dimension of the space, and the order of the BSpline.
> //
> // \index{BSplineDeformableTransform|New}
> // \index{BSplineDeformableTransform|Instantiation}
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> const unsigned int SpaceDimension = ImageDimension;
> const unsigned int SplineOrder = 3;
> typedef double CoordinateRepType;
>
> typedef itk::BSplineDeformableTransform<
> CoordinateRepType,
> 2,
> SplineOrder > TransformType;
> // Software Guide : EndCodeSnippet
>
>
> typedef itk::LBFGSOptimizer OptimizerType;
>
>
> typedef itk::MeanSquaresImageToImageMetric<
> FixedImageType,
> MovingImageType > MetricType;
>
> typedef itk:: LinearInterpolateImageFunction<
> 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 );
>
>
> // Software Guide : BeginLatex
> //
> // The transform object is constructed below and passed to the registration
> // method.
> // \index{itk::RegistrationMethod!SetTransform()}
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> TransformType::Pointer transform = TransformType::New();
> registration->SetTransform( transform );
> // Software Guide : EndCodeSnippet
>
> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
>
> FixedImageReaderType::Pointer fixedImageReader =
> FixedImageReaderType::New();
> MovingImageReaderType::Pointer movingImageReader =
> MovingImageReaderType::New();
>
> itk::DICOMImageIO2Factory::RegisterOneFactory();
>
> fixedImageReader->SetFileName( "Image1.dcm" );
> movingImageReader->SetFileName( "Image2.dcm" );
>
> FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();
> MovingImageType::Pointer movingImage = movingImageReader->GetOutput();
>
>
>
>
> registration->SetFixedImage( fixedImage );
> registration->SetMovingImage( movingImageReader->GetOutput() );
>
> fixedImageReader->Update();
>
> FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
>
> registration->SetFixedImageRegion( fixedRegion );
>
> movingImageReader->Update();
> std::cout << fixedImage->GetBufferedRegion()<<std::endl;
> std::cout << movingImage->GetBufferedRegion()<<std::endl;
>
> // Software Guide : BeginLatex
> //
> // Here we define the parameters of the BSplineDeformableTransform grid. We
> // arbitrarily decide to use a grid with $5 \times 5$ nodes within
> the image.
> // The reader should note that the BSpline computation requires a
> // finite support region ( 1 grid node at the lower borders and 2
> // grid nodes at upper borders). Therefore in this example, we set
> // the grid size to be $8 \times 8$ and place the grid origin such that
> // grid node (1,1) coinicides with the first pixel in the fixed image.
> //
> // \index{BSplineDeformableTransform}
> //
> // Software Guide : EndLatex
>
>
> // Software Guide : BeginCodeSnippet
> typedef TransformType::RegionType RegionType;
> RegionType bsplineRegion;
> RegionType::SizeType gridSizeOnImage;
> RegionType::SizeType gridBorderSize;
> RegionType::SizeType totalGridSize;
>
> gridSizeOnImage.Fill( 5 );
> gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1
> lower, 2 upper )
> totalGridSize = gridSizeOnImage + gridBorderSize;
> std::cout << "totalGridSize: " << totalGridSize << std::endl;
>
> bsplineRegion.SetSize( totalGridSize );
>
> typedef TransformType::SpacingType SpacingType;
> SpacingType spacing =
> fixedImage->GetSpacing();//spacing[0]=1.5625;spacing[1]=1.5625;
> std::cout << "spacing: "<< spacing<<std::endl;
> spacing.Fill(1);
> // fixedImage->SetSpacing(spacing);
> // movingImage->SetSpacing(spacing);
>
> typedef TransformType::OriginType OriginType;
> OriginType origin =
> fixedImage->GetOrigin();//origin[0]=-9.0569601058959961;origin[1]=-191.00799560546875;
> std::cout << "origin: "<< origin <<std::endl;
> origin.Fill(0);
> // fixedImage->SetOrigin(origin);
> // movingImage->SetOrigin(origin);
>
> 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();
>
> std::cout << "numberOfParameters:" << numberOfParameters << std::endl;
>
> ParametersType parameters( numberOfParameters );
>
> parameters.Fill( 0.0 );
>
> transform->SetParameters( parameters );
> // Software Guide : EndCodeSnippet
>
>
>
> // Software Guide : BeginLatex
> //
> // We now pass the parameters of the current transform as the initial
> // parameters to be used when the registration process starts.
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> registration->SetInitialTransformParameters( transform->GetParameters() );
> // Software Guide : EndCodeSnippet
>
> std::cout << "Intial Parameters = " << std::endl;
> std::cout << transform->GetParameters() << std::endl;
>
> // Software Guide : BeginLatex
> //
> // Next we set the parameters of the LBFGS Optimizer.
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> optimizer->SetGradientConvergenceTolerance( 0.05 );
> optimizer->SetLineSearchAccuracy( 0.9 );
> optimizer->SetDefaultStepLength( 1.5 );
> optimizer->TraceOn();
> // optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
> //I think it is too big to waste time
> optimizer->SetMaximumNumberOfFunctionEvaluations( 1000 );
> // Software Guide : EndCodeSnippet
>
>
>
> // 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();
>
> // Software Guide : BeginLatex
> //
> // Let's execute this example using the rat lung images from the
> previous examples.
> //
> // \begin{itemize}
> // \item \code{RatLungSlice1.mha}
> // \item \code{RatLungSlice2.mha}
> // \end{itemize}
> //
> //
> // Software Guide : EndLatex
>
> // Software Guide : BeginCodeSnippet
> transform->SetParameters( finalParameters );
> // Software Guide : EndCodeSnippet
>
>
> typedef itk::ResampleImageFilter<
> MovingImageType,
> FixedImageType > ResampleFilterType;
>
> ResampleFilterType::Pointer resample = ResampleFilterType::New();
>
> resample->SetTransform( transform );
> resample->SetInput( movingImageReader->GetOutput() );
>
> resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
>
>
> // typedef TransformType::SpacingType SpacingType;
> SpacingType spacing1 = fixedImage->GetSpacing(); //not used
> // spacing1.Fill(1);
>
> // typedef TransformType::OriginType OriginType;
> OriginType origin1 = fixedImage->GetOrigin(); //not used
> // origin1.Fill(0);
> // resample->SetOutputOrigin( origin1 );
> // resample->SetOutputSpacing( spacing1 );
>
> resample->SetOutputOrigin( fixedImage->GetOrigin() );
> resample->SetOutputSpacing( fixedImage->GetSpacing() );
> resample->SetDefaultPixelValue( 100 );
>
> typedef unsigned char 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( "output.png" );
>
>
> 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;
> }
>
>
>
> typedef itk::SquaredDifferenceImageFilter<
> FixedImageType,
> FixedImageType,
> OutputImageType > DifferenceFilterType;
>
> DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
>
> WriterType::Pointer writer2 = WriterType::New();
> writer2->SetInput( difference->GetOutput() );
>
>
> // Compute the difference image between the
> // fixed and resampled moving image.
> if( 1 )
> {
> difference->SetInput1( fixedImageReader->GetOutput() );
> difference->SetInput2( resample->GetOutput() );
> writer2->SetFileName( "Diff_fixed_resampled.png" );
> try
> {
> writer2->Update();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
> }
>
>
> // Compute the difference image between the
> // fixed and moving image before registration.
> if( 1)
> {
> writer2->SetFileName( "Diff_fixed_moving.png" );
> difference->SetInput1( fixedImageReader->GetOutput() );
> difference->SetInput2( movingImageReader->GetOutput() );
> try
> {
> writer2->Update();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return -1;
> }
> }
>
>
>
> // Generate the explicit deformation field resulting from
> // the registration.
>
> typedef itk::Vector< float, ImageDimension > VectorType;
> typedef itk::Image< VectorType, ImageDimension > DeformationFieldType;
>
> DeformationFieldType::Pointer field = DeformationFieldType::New();
> field->SetRegions( fixedRegion );
> field->SetOrigin( fixedImage->GetOrigin() );
> field->SetSpacing( fixedImage->GetSpacing() );
> std::cout << "fixedRegion " << fixedRegion << std::endl;
> std::cout << "fixedImage->GetOrigin(): " << fixedImage->GetOrigin()<<
> std::endl;
> std::cout << "fixedImage->GetSpacing(): " << fixedImage->GetSpacing()
> << std::endl;
> // origin1.Fill(0);
> // spacing1.Fill(1);
> // field->SetOrigin(origin);
> // field->SetSpacing(spacing);
> field->Allocate();
>
> typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;
> FieldIterator fi( field, fixedRegion );
>
> fi.GoToBegin();
>
> TransformType::InputPointType fixedPoint;
> TransformType::OutputPointType movingPoint;
> DeformationFieldType::IndexType index;
>
> VectorType displacement;
> int n = 0;
> while( ! fi.IsAtEnd() )
> {
> index = fi.GetIndex();
> field->TransformIndexToPhysicalPoint( index, fixedPoint );
> movingPoint = transform->TransformPoint( fixedPoint );
> displacement[0] = movingPoint[0] - fixedPoint[0];
> displacement[1] = movingPoint[1] - fixedPoint[1];
> fi.Set( displacement );
> ++fi;
> n++;
> }
> std::cout << "Iterator Number: " << n << std::endl;
>
>
>
> typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
> FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
>
> fieldWriter->SetInput( field );
>
> if( 1 )
> {
> fieldWriter->SetFileName( "DeformationField.mhd" );
> try
> {
> fieldWriter->Update();
> }
> catch( itk::ExceptionObject & excp )
> {
> std::cerr << "Exception thrown " << std::endl;
> std::cerr << excp << std::endl;
> return EXIT_FAILURE;
> }
> }
>
> return 0;
> }
>
>
>
> On Sun, 11 Jul 2004 18:21:26 +0800, he smth <kingaza at gmail.com> wrote:
>
>>hi all,
>>I have some questions in using DeformableRegistration4.
>>
>>1. I run the programe DeformableRegistration4.exe and get the
>>deformable field. but then I have a question when viewing the result
>>and the two input file(RatLungSlice1.mha and RatLungSlice2.mha) in
>>paraview: can the result be accepted? I have save it to a .png file in
>>the attachment. I hope that it is my fault and please tell me why.
>>
>>2. when I make it in use of dicom file registration, I am confused
>>with the parameters of origin and spacing. please tell me how to set
>>these values if any successful experience.
>>
>>3. what the parameters means? is it the motion field of each grid
>>node? I have red the code for itkBSplineDeformableTransform, and found
>>that the number of parameters is defined by the size of grid node.
>>
>>4. how to define the parameters of the BSplineDeformableTransform
>>grid? In this example, we set gridSizeOnImage as 5 and gridBorderSize
>>as 3, and so the size of grid node is 8*8.
>>
>>thanks for any help!
>>
>>
>>
>
> _______________________________________________
> 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