<div>Hi,</div> <div> </div> <div>I am trying to register 2 3D labelled volumes, using ITK<FONT color=#008000 size=2></div> <div>BSplineDeformableTransform , metric is set <FONT size=2></div> <div>itkMatchCardinalityImageToImageMetric and <FONT size=2></div> <div>itkLBFGSBOptimizer as it is recommended for registration of labelled volume...</div> <div> </div> <div>I set <FONT size=2>AffineTransform as bulkTransform in the code. </FONT></div> <div><FONT size=2></FONT> </div> <div><FONT size=2><FONT size=2>I have experimented with various on-image grid sizes of ( 8x8x8, 12x12x3 , 27x27x3 ). </FONT></FONT></div> <div> </div> <div>The border grid size was kept fixed.</div> <div> </div> <div>The results were exactly the same as using only affine transforms.</div> <div> </div> <div>Am I doing anything wrong? please advise</div> <div> </div>
<div>*************************************************************************************************</div><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> "itkImageRegistrationMethod.h"</FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> "itkMattesMutualInformationImageToImageMetric.h"</FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> "itkLinearInterpolateImageFunction.h"</FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> "itkImage.h"</FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> "itkTimeProbesCollectorBase.h"</FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> "itkMatchCardinalityImageToImageMetric.h"</FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT
size=2><FONT color=#000000> "itkNearestNeighborInterpolateImageFunction.h"</FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> <conio.h></FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> "itkCenteredTransformInitializer.h"</FONT></div></FONT><FONT color=#0000ff size=2> <div>#include</FONT><FONT size=2><FONT color=#000000> "itkAffineTransform.h"</FONT></div></FONT> <div> </div> <div>#include "itkBSplineDeformableTransform.h"<BR>#include "itkLBFGSBOptimizer.h"<BR>// Software Guide : EndCodeSnippet</div> <div>// Software Guide : BeginLatex<BR>// <BR>// The parameter space of the \code{BSplineDeformableTransform} is composed by<BR>// the set of all the deformations associated with the nodes of the BSpline<BR>// grid. This large number of parameters makes possible to represent a wide<BR>// variety of deformations, but it
also has the price of requiring a<BR>// significant amount of computation time.<BR>//<BR>// \index{itk::BSplineDeformableTransform!header}<BR>// <BR>// Software Guide : EndLatex </div> <div>#include "itkImageFileReader.h"<BR>#include "itkImageFileWriter.h"</div> <div>#include "itkResampleImageFilter.h"<BR>#include "itkCastImageFilter.h"<BR>#include "itkSquaredDifferenceImageFilter.h"</div> <div><BR>// The following section of code implements a Command observer<BR>// used to monitor the evolution of the registration process.<BR>//<BR>#include "itkCommand.h"<BR>class CommandIterationUpdate : public itk::Command <BR>{<BR>public:<BR> typedef CommandIterationUpdate Self;<BR> typedef itk::Command Superclass;<BR> typedef itk::SmartPointer<Self> Pointer;<BR> itkNewMacro( Self );<BR>protected:<BR>
CommandIterationUpdate() {};<BR>public:<BR> typedef itk::LBFGSBOptimizer OptimizerType;<BR> typedef const OptimizerType * OptimizerPointer;</div> <div> void Execute(itk::Object *caller, const itk::EventObject & event)<BR> {<BR> Execute( (const itk::Object *)caller, event);<BR> }</div> <div> void Execute(const itk::Object * object, const itk::EventObject & event)<BR> {<BR> OptimizerPointer optimizer = <BR> dynamic_cast< OptimizerPointer >( object );<BR> if( typeid( event ) != typeid( itk::IterationEvent ) )<BR> {<BR> return;<BR> }<BR>
std::cout << optimizer->GetCurrentIteration() << " ";<BR> std::cout << optimizer->GetValue() << " ";<BR> std::cout << optimizer->GetInfinityNormOfProjectedGradient() << std::endl;<BR> }<BR>};</div> <div><BR>int main( int argc, char *argv[] )<BR>{<BR> if( EXEC_MODE )<BR> {<BR> if( argc < 4 )<BR> {<BR> std::cerr << "Missing Parameters " << std::endl;<BR> std::cerr << "Usage: " << argv[0];<BR> std::cerr << " fixedImageFile movingImageFile outputImagefile ";<BR> std::cerr << " [differenceOutputfile] [differenceBeforeRegistration] ";<BR> std::cerr << "[xy-size] [z-size] "; <BR> std::cerr << " [deformationField] ";<BR> return 1;<BR> }<BR> }
</div> <div><BR> const unsigned int ImageDimension = 3;<BR> typedef unsigned char PixelType;</div> <div> typedef itk::Image< PixelType, ImageDimension > FixedImageType;<BR> typedef itk::Image< PixelType, ImageDimension > MovingImageType;</div> <div><BR> // Software Guide : BeginLatex<BR> //<BR> // We instantiate now the type of the \code{BSplineDeformableTransform} using<BR> // as template parameters the type for coordinates representation, the<BR> // dimension of the space, and the order of the BSpline. <BR> // <BR> // \index{BSplineDeformableTransform|New}<BR> // \index{BSplineDeformableTransform|Instantiation}<BR> //<BR> // Software Guide : EndLatex </div> <div> // Software Guide : BeginCodeSnippet<BR> const unsigned int SpaceDimension =
ImageDimension;<BR> const unsigned int SplineOrder = 3;<BR> typedef double CoordinateRepType;</div> <div> typedef itk::BSplineDeformableTransform<<BR> CoordinateRepType,<BR> SpaceDimension,<BR> SplineOrder > TransformType;</div> <div> typedef itk::AffineTransform< double,
<BR> ImageDimension > BulkTransformType;</div> <div> // Software Guide : EndCodeSnippet</div> <div><BR> typedef itk::LBFGSBOptimizer OptimizerType;</div> <div><BR> /*<BR> typedef itk::MattesMutualInformationImageToImageMetric< <BR> FixedImageType, <BR> MovingImageType
> MetricType;<BR> */<BR> typedef itk::MatchCardinalityImageToImageMetric< <BR> FixedImageType, <BR> MovingImageType > MetricType;<BR> <BR> /*<BR> typedef itk:: LinearInterpolateImageFunction< <BR>
MovingImageType,<BR> double > InterpolatorType;<BR> */<BR> typedef itk::NearestNeighborInterpolateImageFunction< <BR> MovingImageType,<BR> double > InterpolatorType;<BR> <BR>
typedef itk::ImageRegistrationMethod< <BR> FixedImageType, <BR> MovingImageType > RegistrationType;</div> <div> MetricType::Pointer metric = MetricType::New();<BR> metric->MeasureMatchesOff();</div> <div><BR> OptimizerType::Pointer optimizer = OptimizerType::New();<BR> InterpolatorType::Pointer interpolator = InterpolatorType::New();<BR>
RegistrationType::Pointer registration = RegistrationType::New();<BR> </div> <div> registration->SetMetric( metric );<BR> registration->SetOptimizer( optimizer );<BR> registration->SetInterpolator( interpolator );</div> <div><BR> // Software Guide : BeginLatex<BR> //<BR> // The transform object is constructed below and passed to the registration<BR> // method.<BR> // \index{itk::RegistrationMethod!SetTransform()}<BR> //<BR> // Software Guide : EndLatex </div> <div> // Software Guide : BeginCodeSnippet<BR> TransformType::Pointer transform = TransformType::New();<BR> BulkTransformType::Pointer bulkTransform = BulkTransformType::New();</div> <div>
registration->SetTransform( transform );<BR> // Software Guide : EndCodeSnippet</div> <div> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;<BR> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;</div> <div> FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();<BR> MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();</div> <div> if( EXEC_MODE )<BR> fixedImageReader->SetFileName( argv[1] );<BR> else<BR> fixedImageReader->SetFileName( STR_FIXED_IMAGE );</div> <div> fixedImageReader->Update();</div> <div> if( EXEC_MODE )<BR> movingImageReader->SetFileName( argv[2] );<BR> else<BR> movingImageReader->SetFileName( STR_MOVING_IMAGE );</div> <div> movingImageReader->Update();</div> <div>
FixedImageType::ConstPointer fixedImage = fixedImageReader->GetOutput();</div> <div> registration->SetFixedImage( fixedImage );<BR> registration->SetMovingImage( movingImageReader->GetOutput() );</div> <div> FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();<BR> <BR> registration->SetFixedImageRegion( fixedRegion );</div> <div> // Software Guide : BeginLatex<BR> //<BR> // Here we define the parameters of the BSplineDeformableTransform grid. We<BR> // arbitrarily decide to use a grid with $5 \times 5$ nodes within the image. <BR> // The reader should note that the BSpline computation requires a<BR> // finite support region ( 1 grid node at the lower borders and 2<BR> // grid nodes at upper borders). Therefore in this example, we set<BR> // the grid size to be $8
\times 8$ and place the grid origin such that<BR> // grid node (1,1) coinicides with the first pixel in the fixed image.<BR> // <BR> // \index{BSplineDeformableTransform}<BR> //<BR> // Software Guide : EndLatex </div> <div> // setup affine transform<BR> typedef itk::CenteredTransformInitializer< <BR> BulkTransformType, <BR> FixedImageType,
<BR> MovingImageType > TransformInitializerType;<BR> TransformInitializerType::Pointer initializer = TransformInitializerType::New();<BR> initializer->SetTransform( bulkTransform );<BR> initializer->SetFixedImage( fixedImageReader->GetOutput() );<BR> initializer->SetMovingImage( movingImageReader->GetOutput() );<BR> initializer->MomentsOn();<BR> initializer->InitializeTransform();</div> <div> <BR> // Software Guide : BeginCodeSnippet<BR> typedef TransformType::RegionType RegionType;<BR> RegionType bsplineRegion;<BR> RegionType::SizeType gridSizeOnImage;<BR> RegionType::SizeType gridBorderSize;<BR>
RegionType::SizeType totalGridSize;</div> <div> gridSizeOnImage.Fill( atoi( argv[6] ) );<BR> gridSizeOnImage[2] = atoi( argv[7] );<BR> gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper )<BR> totalGridSize = gridSizeOnImage + gridBorderSize;</div> <div> bsplineRegion.SetSize( totalGridSize );</div> <div> typedef TransformType::SpacingType SpacingType;<BR> SpacingType spacing = fixedImage->GetSpacing();</div> <div> typedef TransformType::OriginType OriginType;<BR> OriginType origin = fixedImage->GetOrigin();;</div> <div> FixedImageType::SizeType fixedImageSize = fixedRegion.GetSize();</div> <div> for(unsigned int r=0; r<ImageDimension; r++)<BR> {<BR> spacing[r] *= floor( static_cast<double>(fixedImageSize[r] - 1) /
<BR> static_cast<double>(gridSizeOnImage[r] - 1) );<BR> origin[r] -= spacing[r]; <BR> }</div> <div> transform->SetGridSpacing( spacing );<BR> transform->SetGridOrigin( origin );<BR> transform->SetGridRegion( bsplineRegion );<BR> </div> <div> typedef TransformType::ParametersType ParametersType;</div> <div> const unsigned int numberOfParameters =<BR> transform->GetNumberOfParameters();<BR> <BR> ParametersType parameters( numberOfParameters );</div> <div> parameters.Fill( 0.0 );</div> <div> transform->SetParameters( parameters );</div> <div> transform->SetBulkTransform( bulkTransform );</div> <div> // Software
Guide : EndCodeSnippet</div> <div> </div> <div> // Software Guide : BeginLatex<BR> // <BR> // We now pass the parameters of the current transform as the initial<BR> // parameters to be used when the registration process starts. <BR> //<BR> // Software Guide : EndLatex </div> <div> // Software Guide : BeginCodeSnippet<BR> registration->SetInitialTransformParameters( transform->GetParameters() );<BR> // Software Guide : EndCodeSnippet</div> <div> // std::cout << "Intial Parameters = " << std::endl;<BR> // std::cout << transform->GetParameters() << std::endl;</div> <div> // Software Guide : BeginLatex<BR> // <BR> // Next we set the parameters of the LBFGSB Optimizer. <BR> //<BR> // Software Guide : EndLatex </div> <div><BR> // Software Guide : BeginCodeSnippet<BR>
OptimizerType::BoundSelectionType boundSelect( transform->GetNumberOfParameters() );<BR> OptimizerType::BoundValueType upperBound( transform->GetNumberOfParameters() );<BR> OptimizerType::BoundValueType lowerBound( transform->GetNumberOfParameters() );</div> <div> boundSelect.Fill( 0 );<BR> upperBound.Fill( 0.0 );<BR> lowerBound.Fill( 0.0 );</div> <div> optimizer->SetBoundSelection( boundSelect );<BR> optimizer->SetUpperBound( upperBound );<BR> optimizer->SetLowerBound( lowerBound );</div> <div> optimizer->SetCostFunctionConvergenceFactor( 1e+7 );<BR> optimizer->SetProjectedGradientTolerance( 1e-4 );<BR> optimizer->SetMaximumNumberOfIterations( 500 );<BR> optimizer->SetMaximumNumberOfEvaluations( 500 );<BR> optimizer->SetMaximumNumberOfCorrections( 12 );</div> <div> // Software Guide : EndCodeSnippet</div> <div> // Create the Command observer and
register it with the optimizer.<BR> //<BR> CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();<BR> optimizer->AddObserver( itk::IterationEvent(), observer );</div> <div><BR> // Software Guide : BeginLatex<BR> // <BR> // Next we set the parameters of the Mattes Mutual Information Metric. <BR> //<BR> // Software Guide : EndLatex </div> <div> // Software Guide : BeginCodeSnippet<BR> //metric->SetNumberOfHistogramBins( 50 );<BR> <BR> //const unsigned int numberOfSamples = fixedRegion.GetNumberOfPixels() / 10;</div> <div> //metric->SetNumberOfSpatialSamples( numberOfSamples );<BR> // Software Guide : EndCodeSnippet<BR> </div> <div> // Software Guide : BeginLatex<BR> // <BR> // Given that the Mattes Mutual Information metric uses a random iterator in<BR> // order to collect the samples
from the images, it is usually convenient to<BR> // initialize the seed of the random number generator.<BR> //<BR> // Software Guide : EndLatex </div> <div> // Software Guide : BeginCodeSnippet<BR> vnl_sample_reseed( 76926294 );<BR> // Software Guide : EndCodeSnippet</div> <div><BR> // Add a time probe<BR> itk::TimeProbesCollectorBase collector;</div> <div> std::cout << std::endl << "Starting Registration" << std::endl;</div> <div> try <BR> { <BR> collector.Start( "Registration" );<BR> registration->StartRegistration(); <BR> collector.Stop( "Registration" );<BR> } <BR> catch( itk::ExceptionObject & err ) <BR> { <BR> std::cerr << "ExceptionObject caught !" << std::endl; <BR> std::cerr << err <<
std::endl; <BR> return -1;<BR> } <BR> <BR> OptimizerType::ParametersType finalParameters = <BR> registration->GetLastTransformParameters();</div> <div>// std::cout << "Last Transform Parameters" << std::endl;<BR>// std::cout << finalParameters << std::endl;</div> <div><BR> // Report the time taken by the registration<BR> collector.Report();</div> <div> // Software Guide : BeginCodeSnippet<BR> transform->SetParameters( finalParameters );<BR> // Software Guide : EndCodeSnippet</div> <div><BR> typedef itk::ResampleImageFilter< <BR> MovingImageType,
<BR> FixedImageType > ResampleFilterType;</div> <div> ResampleFilterType::Pointer resample = ResampleFilterType::New();</div> <div> resample->SetTransform( transform );<BR> resample->SetInput( movingImageReader->GetOutput() );</div> <div> resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );<BR> resample->SetOutputOrigin( fixedImage->GetOrigin() );<BR> resample->SetOutputSpacing( fixedImage->GetSpacing() );<BR> resample->SetDefaultPixelValue( 0 );<BR> resample->SetInterpolator( interpolator );<BR> <BR> typedef unsigned char OutputPixelType;</div> <div> typedef itk::Image< OutputPixelType, ImageDimension > OutputImageType;<BR>
<BR> typedef itk::CastImageFilter< <BR> FixedImageType,<BR> OutputImageType > CastFilterType;<BR> <BR> typedef itk::ImageFileWriter< OutputImageType > WriterType;</div> <div><BR> WriterType::Pointer writer = WriterType::New();<BR> CastFilterType::Pointer caster = CastFilterType::New();</div> <div><BR> // write aligned image<BR> if( EXEC_MODE )<BR> writer->SetFileName( argv[3] );<BR> else<BR> writer->SetFileName( STR_OUTPUT_IMAGE );</div> <div><BR>
caster->SetInput( resample->GetOutput() );<BR> writer->SetInput( caster->GetOutput() );</div> <div><BR> try<BR> {<BR> writer->Update();<BR> }<BR> catch( itk::ExceptionObject & err ) <BR> { <BR> std::cerr << "ExceptionObject caught !" << std::endl; <BR> std::cerr << err << std::endl; <BR> return -1;<BR> } <BR> </div> <div><BR> typedef itk::SquaredDifferenceImageFilter< <BR> FixedImageType,
<BR> FixedImageType, <BR> OutputImageType > DifferenceFilterType;</div> <div> DifferenceFilterType::Pointer difference = DifferenceFilterType::New();</div> <div> WriterType::Pointer writer2 = WriterType::New();<BR> writer2->SetInput( difference->GetOutput() ); <BR> </div> <div> // Compute the difference image between the <BR> // fixed and resampled moving image.<BR> difference->SetInput1( fixedImageReader->GetOutput() );<BR> difference->SetInput2( resample->GetOutput() );</div> <div> if(
EXEC_MODE )<BR> writer2->SetFileName( argv[5] );<BR> else<BR> writer2->SetFileName( STR_DIFF_AFTER_IMAGE );</div> <div> try<BR> {<BR> writer2->Update();<BR> }<BR> catch( itk::ExceptionObject & err ) <BR> { <BR> std::cerr << "ExceptionObject caught !" << std::endl; <BR> std::cerr << err << std::endl; <BR> return -1;<BR> } </div> <div> // Compute the difference image between the <BR> // fixed and moving image before registration.</div> <div> <BR> if( EXEC_MODE )<BR> writer2->SetFileName( argv[4] );<BR> else<BR> writer2->SetFileName( STR_DIFF_BEFORE_IMAGE );</div> <div>
difference->SetInput1( fixedImageReader->GetOutput() );<BR> difference->SetInput2( movingImageReader->GetOutput() );<BR> try<BR> {<BR> writer2->Update();<BR> }<BR> catch( itk::ExceptionObject & err ) <BR> { <BR> std::cerr << "ExceptionObject caught !" << std::endl; <BR> std::cerr << err << std::endl; <BR> return -1;<BR> } </div> <div><BR> // Generate the explicit deformation field resulting from <BR> // the registration.<BR> if( EXEC_MODE && argc >= 8 )<BR> {</div> <div> typedef itk::Vector< float, ImageDimension > VectorType;<BR> typedef
itk::Image< VectorType, ImageDimension > DeformationFieldType;</div> <div> DeformationFieldType::Pointer field = DeformationFieldType::New();<BR> field->SetRegions( fixedRegion );<BR> field->SetOrigin( fixedImage->GetOrigin() );<BR> field->SetSpacing( fixedImage->GetSpacing() );<BR> field->Allocate();</div> <div> typedef itk::ImageRegionIterator< DeformationFieldType > FieldIterator;<BR> FieldIterator fi( field, fixedRegion );</div> <div> fi.GoToBegin();</div> <div> TransformType::InputPointType fixedPoint;<BR> TransformType::OutputPointType movingPoint;<BR> DeformationFieldType::IndexType index;</div> <div> VectorType displacement;</div> <div> while( ! fi.IsAtEnd()
)<BR> {<BR> index = fi.GetIndex();<BR> field->TransformIndexToPhysicalPoint( index, fixedPoint );<BR> movingPoint = transform->TransformPoint( fixedPoint );<BR> displacement[0] = movingPoint[0] - fixedPoint[0];<BR> displacement[1] = movingPoint[1] - fixedPoint[1];<BR> fi.Set( displacement );<BR> ++fi;<BR> }</div> <div> </div> <div> typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;<BR> FieldWriterType::Pointer fieldWriter = FieldWriterType::New();</div> <div> fieldWriter->SetInput( field );</div> <div> fieldWriter->SetFileName( argv[8] );<BR>
try<BR> {<BR> fieldWriter->Update();<BR> }<BR> catch( itk::ExceptionObject & excp )<BR> {<BR> std::cerr << "Exception thrown " << std::endl;<BR> std::cerr << excp << std::endl;<BR> return EXIT_FAILURE;<BR> }<BR> }</div> <div> return 0;<BR>}</div> <div> </div> <div><FONT size=2><FONT size=2></FONT></FONT> </div> <div> </div> <div><FONT size=2><FONT size=2> </div></FONT></FONT></FONT></FONT></FONT><p> 
<hr size=1>Get your email and see which of your friends are online - Right on the <a href="http://us.rd.yahoo.com/evt=42973/*http://www.yahoo.com/preview"> new Yahoo.com</a>