X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=registration%2FclitkAffineRegistrationGenericFilter.cxx;h=2c94379e99e1f8411710e20e968b82c96b3467a5;hb=HEAD;hp=d9bfd2892618e2e8fcb9a0a16baa741a42fb9e7f;hpb=c18059db4f507fd31b5898667f57eced7d48c5f7;p=clitk.git diff --git a/registration/clitkAffineRegistrationGenericFilter.cxx b/registration/clitkAffineRegistrationGenericFilter.cxx index d9bfd28..2c94379 100644 --- a/registration/clitkAffineRegistrationGenericFilter.cxx +++ b/registration/clitkAffineRegistrationGenericFilter.cxx @@ -3,7 +3,7 @@ Authors belong to: - University of LYON http://www.universite-lyon.fr/ - - Léon Bérard cancer center http://oncora1.lyon.fnclcc.fr + - Léon Bérard cancer center http://www.centreleonberard.fr - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr This software is distributed WITHOUT ANY WARRANTY; without even @@ -14,83 +14,736 @@ - BSD See included LICENSE.txt file - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html -======================================================================-====*/ +===========================================================================**/ #ifndef CLITKAFFINEREGISTRATIONGENERICFILTER_CXX #define CLITKAFFINEREGISTRATIONGENERICFILTER_CXX #include "clitkAffineRegistrationGenericFilter.h" +// clitk include +#include "clitkIO.h" +#include "clitkCommon.h" +#include "clitkImageCommon.h" +#include "clitkAffineRegistration_ggo.h" +#include "clitkImageArithm_ggo.h" +#include "clitkCorrelationRatioImageToImageMetric.h" +#include "clitkTransformUtilities.h" +#include "clitkGenericMetric.h" +#include "clitkGenericOptimizer.h" +#include "clitkGenericInterpolator.h" +#include "clitkGenericAffineTransform.h" +#include "clitkImageToImageGenericFilter.h" + + +//itk include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// other includes +#include +#include +#include + namespace clitk { +//============================================================================== +//Creating an observer class that allows us to monitor the registration +//================================================================================ +class CommandIterationUpdate : public itk::Command +{ +public: + typedef CommandIterationUpdate Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer Pointer; + itkNewMacro( Self ); +protected: + CommandIterationUpdate() {}; +public: + typedef clitk::GenericOptimizer OptimizerType; + typedef const OptimizerType * OptimizerPointer; + + // Set the generic optimizer + void SetOptimizer(OptimizerPointer o) { + m_Optimizer=o; + } + + // Execute + void Execute(itk::Object *caller, const itk::EventObject & event) ITK_OVERRIDE { + Execute( (const itk::Object *)caller, event); + } + + void Execute(const itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE { + if ( !(itk::IterationEvent().CheckEvent( &event )) ) { + return; + } + + m_Optimizer->OutputIterationInfo(); + } + + OptimizerPointer m_Optimizer; +}; +//==================================================================================================================================// +//Constructor +//===================================================================================================================================// + +AffineRegistrationGenericFilter::AffineRegistrationGenericFilter(): + ImageToImageGenericFilter("Register") -//==================================================================== -// Constructor -AffineRegistrationGenericFilter::AffineRegistrationGenericFilter() { - m_Verbose=false; + InitializeImageType<2>(); + InitializeImageType<3>(); + m_Verbose=true; } +//==========================================================================================================// +//============================================================================================================// +//-------------------------------------------------------------------- +template +void AffineRegistrationGenericFilter::InitializeImageType() +{ + ADD_DEFAULT_IMAGE_TYPES(Dim); +} +//-------------------------------------------------------------------- -//==================================================================== -// Update -void AffineRegistrationGenericFilter::Update() + +//============================================================================== +//Creating an observer class that allows us to change parameters at subsequent levels +//============================================================================== +template +class RegistrationInterfaceCommand : public itk::Command +{ +public: + typedef RegistrationInterfaceCommand Self; + typedef itk::Command Superclass; + typedef itk::SmartPointer Pointer; + itkNewMacro( Self ); +protected: + RegistrationInterfaceCommand() {}; +public: + + // Registration + typedef TRegistration RegistrationType; + typedef RegistrationType * RegistrationPointer; + + // Metric + typedef typename RegistrationType::FixedImageType InternalImageType; + typedef clitk::GenericMetric GenericMetricType; + + // Two arguments are passed to the Execute() method: the first + // is the pointer to the object which invoked the event and the + // second is the event that was invoked. + void Execute(itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE { + if ( !(itk::IterationEvent().CheckEvent( &event )) ) { + return; + } + + // Get the levels + RegistrationPointer registration = dynamic_cast( object ); + unsigned int numberOfLevels=registration->GetNumberOfLevels(); + unsigned int currentLevel=registration->GetCurrentLevel()+1; + + // Output the levels + std::cout<1) { + // Reinitialize the metric (!= number of samples) + typename GenericMetricType::Pointer genericMetric= GenericMetricType::New(); + genericMetric->SetArgsInfo(m_ArgsInfo); + genericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel())); + if (m_ArgsInfo.referenceMask_given) genericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask()); + typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType; + typename MetricType::Pointer metric=genericMetric->GetMetricPointer(); + registration->SetMetric(metric); + } + } + + void Execute(const itk::Object * , const itk::EventObject & ) ITK_OVERRIDE { + return; + } + + void SetArgsInfo(args_info_clitkAffineRegistration a) { + m_ArgsInfo=a; + } + args_info_clitkAffineRegistration m_ArgsInfo; +}; + +//==============================================================================================// +// ArgsInfo +//==============================================================================================// +void AffineRegistrationGenericFilter::SetArgsInfo(const args_info_clitkAffineRegistration & a) +{ + m_ArgsInfo=a; + if (m_ArgsInfo.reference_given) AddInputFilename(m_ArgsInfo.reference_arg); + + if (m_ArgsInfo.target_given) { + AddInputFilename(m_ArgsInfo.target_arg); + } + + if (m_ArgsInfo.output_given) SetOutputFilename(m_ArgsInfo.output_arg); +} +//============================================================================== +// Update with the number of dimensions and pixeltype +//============================================================================== +template +void AffineRegistrationGenericFilter::UpdateWithInputImageType() { - //Read the PixelType and dimension of the reference image - int Dimension; - std::string PixelType; - clitk::ReadImageDimensionAndPixelType(m_ArgsInfo.reference_arg, Dimension, PixelType); - if (Dimension == 2) UpdateWithDim<2>(PixelType); - else if (Dimension == 3) UpdateWithDim<3>(PixelType); - else { - itkExceptionMacro(<< "Reference Image dimension is " << Dimension - << " but I can only work on 2D and 3D images."); + //============================================================================= + //Input + //============================================================================= + + typedef typename InputImageType::PixelType PixelType; +//typedef typename InputImageType::ImageDimension Dimension; + + bool threadsGiven=m_ArgsInfo.threads_given; + int threads=m_ArgsInfo.threads_arg; + + //Coordinate Representation + typedef double TCoordRep; + + + typename InputImageType::Pointer fixedImage = this->template GetInput(0); + + typename InputImageType::Pointer inputFixedImage = this->template GetInput(0); + + // typedef input2 + typename InputImageType::Pointer movingImage = this->template GetInput(1); + + typename InputImageType::Pointer inputMovingImage = this->template GetInput(1); + + + + //The pixeltype of the fixed image will be used for output + typedef itk::Image< PixelType, InputImageType::ImageDimension > FixedImageType; + + //Whatever the pixel type, internally we work with an image represented in float + typedef typename InputImageType::PixelType InternalPixelType; + typedef itk::Image< PixelType, InputImageType::ImageDimension > InternalImageType; + + + //Read in the reference/fixed image +// typedef itk::ImageFileReader< InternalImageType > ReaderType; +// typename ReaderType::Pointer fixedImageReader = ReaderType::New(); +// fixedImageReader->SetFileName( m_ArgsInfo.reference_arg); + + + //Read in the object/moving image +// typename ReaderType::Pointer movingImageReader = ReaderType::New(); +// movingImageReader->SetFileName( m_ArgsInfo.target_arg ); + if (m_Verbose) std::cout<<"Reading images..."<Update(); +// movingImageReader->Update(); + + if (m_Verbose) std::cout << "Reading images... " << std::endl; + + //we connect pointers to these internal images + // typedef typename fixedImageReader fixedImage; + // typedef typename movingImageReader movingImage; + + //We keep the images used for input for possible output +// typedef typename fixedImageReader inputFixedImage; +// typedef typename movingImageReader inputMovingImage; + + + //============================================================================ + // Preprocessing + //============================================================================ + + //If given, the intensities of both images are first normalized to a zero mean and SD of 1 + // (usefull for MI, not necessary for Mattes' MI but performed anyway for the ouput) + if ( m_ArgsInfo.normalize_flag ) { + typedef itk::NormalizeImageFilter< InternalImageType,InternalImageType > NormalizeFilterType; + + typename NormalizeFilterType::Pointer fixedNormalizeFilter = NormalizeFilterType::New(); + typename NormalizeFilterType::Pointer movingNormalizeFilter = NormalizeFilterType::New(); + + fixedNormalizeFilter->SetInput( fixedImage ); + movingNormalizeFilter->SetInput( movingImage ); + + fixedNormalizeFilter->Update(); + movingNormalizeFilter->Update(); + + //We keep the images used for input for possible output + inputFixedImage= fixedNormalizeFilter->GetOutput(); + inputMovingImage= movingNormalizeFilter->GetOutput(); + + //the pointers are reconnected for further output + fixedImage=fixedNormalizeFilter->GetOutput(); + movingImage=movingNormalizeFilter->GetOutput(); + + if (m_Verbose) std::cout << "Normalizing image intensities to zero mean and SD of 1..." << std::endl; + } + + + //If given, the images are blurred before processing + if ( m_ArgsInfo.blur_arg!= 0) { + typedef itk::DiscreteGaussianImageFilter GaussianFilterType; + typename GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New(); + typename GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New(); + fixedSmoother->SetVariance( m_ArgsInfo.blur_arg ); + movingSmoother->SetVariance(m_ArgsInfo.blur_arg ); + + fixedSmoother->SetInput( fixedImage ); + movingSmoother->SetInput( movingImage ); + + fixedSmoother->Update(); + movingSmoother->Update(); + + fixedImage=fixedSmoother->GetOutput(); + movingImage=movingSmoother->GetOutput(); + + if (m_Verbose) std::cout << "Blurring images with a Gaussian with standard deviation of " << m_ArgsInfo.blur_arg <<"..." << std::endl; + } + + + //============================================================================ + // Setting up the moving image in a reference system + //============================================================================ + const itk::Vector movingResolution = movingImage->GetSpacing(); + typename InternalImageType::RegionType movingRegion = movingImage->GetLargestPossibleRegion(); + typename InternalImageType::RegionType::SizeType movingSize = movingRegion.GetSize(); + + // Print the parameters of the moving image + if (m_Verbose) { + std::cout << "Object or Moving image:"< fixedResolution = fixedImage->GetSpacing(); + typename InternalImageType::RegionType fixedRegion = fixedImage->GetLargestPossibleRegion(); + typename InternalImageType::RegionType::SizeType fixedSize = fixedRegion.GetSize(); + + // Print the parameters of the moving image and the transform + if (m_Verbose) { + std::cout << "Target or Moving image:"< MaskType; + typename MaskType::Pointer fixedMask=ITK_NULLPTR; + if (m_ArgsInfo.referenceMask_given) { + fixedMask= MaskType::New(); + typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType; + typedef itk::ImageFileReader< ImageMaskType > MaskReaderType; + typename MaskReaderType::Pointer maskReader = MaskReaderType::New(); + maskReader->SetFileName(m_ArgsInfo.referenceMask_arg); + try { + maskReader->Update(); + } catch ( itk::ExceptionObject & err ) { + std::cerr << "ExceptionObject caught while reading mask !" << std::endl; + std::cerr << err << std::endl; + return; + } + if (m_Verbose)std::cout <<"Reference image mask was read..." <SetImage( maskReader->GetOutput() ); + } + + typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType; + typename MaskType::Pointer movingMask=ITK_NULLPTR; + if (m_ArgsInfo.targetMask_given) { + movingMask= MaskType::New(); + typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType; + typedef itk::ImageFileReader< ImageMaskType > MaskReaderType; + typename MaskReaderType::Pointer maskReader = MaskReaderType::New(); + maskReader->SetFileName(m_ArgsInfo.targetMask_arg); + try { + maskReader->Update(); + } catch ( itk::ExceptionObject & err ) { + std::cerr << "ExceptionObject caught !" << std::endl; + std::cerr << err << std::endl; + } + if (m_Verbose)std::cout <<"Target image mask was read..." <SetImage( maskReader->GetOutput() ); + } + + + //============================================================================ + // The image pyramids + //============================================================================ + typedef itk::RecursiveMultiResolutionPyramidImageFilter FixedImagePyramidType; + typedef itk::RecursiveMultiResolutionPyramidImageFilter MovingImagePyramidType; + typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New(); + typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New(); + fixedImagePyramid->SetUseShrinkImageFilter(false); + fixedImagePyramid->SetInput(fixedImage); + fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg); + movingImagePyramid->SetUseShrinkImageFilter(false); + movingImagePyramid->SetInput(movingImage); + movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg); + if (m_Verbose) std::cout<<"Creating the image pyramid..."<Update(); + movingImagePyramid->Update(); + + + + //============================================================================ + // We retrieve the type of metric from the command line + //============================================================================ + typedef clitk::GenericMetric GenericMetricType; + typename GenericMetricType::Pointer genericMetric=GenericMetricType::New(); + genericMetric->SetArgsInfo(m_ArgsInfo); + genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0)); + if (fixedMask) genericMetric->SetFixedImageMask(fixedMask); + typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType; + typename MetricType::Pointer metric=genericMetric->GetMetricPointer(); + if (movingMask) metric->SetMovingImageMask(movingMask); + + if (threadsGiven) { +#if ITK_VERSION_MAJOR <= 4 + metric->SetNumberOfThreads( threads ); +#else + metric->SetNumberOfWorkUnits( threads ); +#endif + } + + //============================================================================ + // Initialize using image moments. + //============================================================================ + if (m_ArgsInfo.moment_flag) { + typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType; + typename CalculatorType::Pointer fixedCalculator= CalculatorType::New(); + + typename InternalImageType::Pointer fixedThresh; + if (m_ArgsInfo.intThreshold_given) { + typedef itk::ThresholdImageFilter ThresholdImageFilterType; + typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New(); + thresholder->SetInput(fixedImage); + thresholder->SetLower(m_ArgsInfo.intThreshold_arg); + thresholder->Update(); + fixedThresh=thresholder->GetOutput(); + } else fixedThresh=fixedImage; + + fixedCalculator->SetImage(fixedThresh); + fixedCalculator->Compute(); + Vector fixedCenter=fixedCalculator->GetCenterOfGravity(); + if (m_Verbose)std::cout<<"The fixed center of gravity is "< CalculatorType; + typename CalculatorType::Pointer movingCalculator= CalculatorType::New(); + + typename InternalImageType::Pointer movingThresh; + if (m_ArgsInfo.intThreshold_given) { + typedef itk::ThresholdImageFilter ThresholdImageFilterType; + typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New(); + thresholder->SetInput(movingImage); + thresholder->SetLower(m_ArgsInfo.intThreshold_arg); + thresholder->Update(); + movingThresh=thresholder->GetOutput(); + } else movingThresh=movingImage; + + movingCalculator->SetImage(movingThresh); + movingCalculator->Compute(); + Vector movingCenter=movingCalculator->GetCenterOfGravity(); + if (m_Verbose)std::cout<<"The moving center of gravity is "< shift= movingCenter-fixedCenter; + if (m_Verbose)std::cout<<"The initial shift applied is "< GenericAffineTransformType; + typename GenericAffineTransformType::Pointer genericAffineTransform = GenericAffineTransformType::New(); + genericAffineTransform->SetArgsInfo(m_ArgsInfo); + typedef itk::Transform< double, InputImageType::ImageDimension, InputImageType::ImageDimension > TransformType; + typename TransformType::Pointer transform = genericAffineTransform->GetTransform(); + std::cout< GenericInterpolatorType; + typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New(); + genericInterpolator->SetArgsInfo(m_ArgsInfo); + typedef itk::InterpolateImageFunction< InternalImageType, TCoordRep > InterpolatorType; + typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer(); + std::cout<<"end of interpolator"< GenericOptimizerType; + unsigned int nParam = transform->GetNumberOfParameters(); + typename GenericOptimizerType::Pointer genericOptimizer=GenericOptimizerType::New(); + genericOptimizer->SetArgsInfo(m_ArgsInfo); + genericOptimizer->SetOutputIteration(m_Verbose); + genericOptimizer->SetOutputPosition(m_Verbose); + genericOptimizer->SetOutputValue(m_Verbose); + genericOptimizer->SetOutputGradient(m_ArgsInfo.gradient_flag); + genericOptimizer->SetMaximize(genericMetric->GetMaximize()); + genericOptimizer->SetNumberOfParameters(nParam); + typedef itk::SingleValuedNonLinearOptimizer OptimizerType; + OptimizerType::Pointer optimizer=genericOptimizer->GetOptimizerPointer(); + + // Scales + itk::Optimizer::ScalesType scales( nParam ); + for (unsigned int i=nParam-InputImageType::ImageDimension; iSetScales(scales); + //============================================================================ + // Multiresolution registration + //============================================================================ + std::cout<<"start MultiResolution..."< RegistrationType; + typename RegistrationType::Pointer registration = RegistrationType::New(); + registration->SetFixedImage( fixedImage ); + registration->SetFixedImageRegion(fixedImage->GetLargestPossibleRegion()); + registration->SetMovingImage( movingImage ); + registration->SetFixedImagePyramid( fixedImagePyramid ); + registration->SetMovingImagePyramid( movingImagePyramid ); + registration->SetTransform( transform ); + registration->SetInitialTransformParameters( transform->GetParameters() ); + registration->SetInterpolator( interpolator ); + registration->SetMetric(metric); + registration->SetOptimizer(optimizer); + registration->SetNumberOfLevels( m_ArgsInfo.levels_arg ); + if (m_Verbose) std::cout << "Setting "<< m_ArgsInfo.levels_arg <<" resolution levels..." << std::endl; + if (m_Verbose) std::cout << "Initial Transform: "<< registration->GetInitialTransformParameters()<SetOptimizer(genericOptimizer); + optimizer->AddObserver( itk::IterationEvent(), observer ); + + + // Output level info + typedef RegistrationInterfaceCommand CommandType; + typename CommandType::Pointer command = CommandType::New(); + command->SetArgsInfo(m_ArgsInfo); + registration->AddObserver( itk::IterationEvent(), command ); + } + + + //============================================================================ + // Finally we can start the registration with the given amount of multiresolution levels + //============================================================================ + if (m_Verbose) std::cout << "Starting the registration now..." << std::endl; + + try { + registration->Update(); + } catch ( itk::ExceptionObject & err ) { + std::cerr << "ExceptionObject caught !" << std::endl; + std::cerr << err << std::endl; + } + + + //============================================================================ + // Processing the result of the registration + //============================================================================ + OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters(); + std::cout<< "Result : " <,itk::Image >::Pointer & registration ) -//{ -// typedef itk::Euler2DTransform < double > Transform2DType; -// Transform2DType::Pointer t2 =Transform2DType::New(); -// -// //Initializing the transform -// Transform2DType::OutputVectorType translation; -// translation[0] = m_ArgsInfo.transX_arg; -// translation[1] = m_ArgsInfo.transY_arg; -// t2->SetTranslation(translation); -// t2->SetRotation(m_ArgsInfo.rotX_arg); -// -// //For simplicity we set the center to the top left corner -// Transform2DType::InputPointType center; -// center[0] = 0.; -// center[1] = 0.; -// t2->SetCenter(center); -// registration->SetTransform(t2); -// registration->SetInitialTransformParameters(t2->GetParameters()); -//} - -//void clitk::AffineRegistrationGenericFilter::SetTransfo( itk::MultiResolutionImageRegistrationMethod,itk::Image >::Pointer & registration ) -//{ -// typedef itk::Euler3DTransform < double > Transform3DType; -// Transform3DType::Pointer t3 = Transform3DType::New(); -// -// t3->SetComputeZYX(true); -// //Initializing the transform -// Transform3DType::OutputVectorType translation; -// translation[0] = m_ArgsInfo.transX_arg; -// translation[1] = m_ArgsInfo.transY_arg; -// translation[2] = m_ArgsInfo.transZ_arg; -// t3->SetTranslation(translation); -// t3->SetRotation(m_ArgsInfo.rotX_arg, m_ArgsInfo.rotY_arg, m_ArgsInfo.rotZ_arg); -// -// //For simplicity we set the center to the top left corner -// Transform3DType::InputPointType center; -// center[0] = 0.; -// center[1] = 0.; -// center[2] = 0.; -// t3->SetCenter(center); -// registration->SetTransform(t3); -// registration->SetInitialTransformParameters( t3->GetParameters()); -//} + std::cout << " Affine transform matrix =" << std::endl; + std::cout << matrix <SetParameters( finalParameters ); + typedef itk::ResampleImageFilter< InternalImageType,InternalImageType > ResampleFilterType; + typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); + + resampler->SetTransform( transform ); + resampler->SetInput( movingImage ); + resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() ); + resampler->SetOutputOrigin( fixedImage->GetOrigin() ); + resampler->SetOutputSpacing( fixedImage->GetSpacing() ); + resampler->SetDefaultPixelValue( 0 ); + resampler->Update(); + //Output? + // if (m_ArgsInfo.output_given) { + //We write an output in the same pixeltype then the input + /*typedef itk::ImageFileWriter< FixedImageType > WriterType; + typename WriterType::Pointer outputWriter = WriterType::New(); + outputWriter->SetFileName(m_ArgsInfo.output_arg ); + outputWriter->SetInput( resampler->GetOutput() ); + outputWriter->Update();*/ + typedef InternalImageType OutputImageType; + typename OutputImageType::Pointer outputImage = resampler->GetOutput(); + std::cout<<"Writing Output....."<template SetNextOutput(outputImage); + // } + + + //============================================================================ + // Checker after? + //============================================================================ + if (m_ArgsInfo.checker_after_given) { + //To display correctly the checkerboard image, the intensities must lie in the same range (normalized) + //We write the image in the internal image type + typedef itk::ResampleImageFilter< InternalImageType,InternalImageType > ResampleFilterType; + typename ResampleFilterType::Pointer internalResampler = ResampleFilterType::New(); + internalResampler->SetTransform( transform ); + internalResampler->SetInput( inputMovingImage ); + internalResampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() ); + internalResampler->SetOutputOrigin( fixedImage->GetOrigin() ); + internalResampler->SetOutputSpacing( fixedImage->GetSpacing() ); + internalResampler->SetDefaultPixelValue( 0 ); + + //We pass the normalized images to the checker filter + typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType; + typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New(); + + checkerFilter->SetInput1(inputFixedImage); + checkerFilter->SetInput2(internalResampler->GetOutput()); + typedef itk::ImageFileWriter< InternalImageType > InternalWriterType; + typename InternalWriterType::Pointer checkerWriter = InternalWriterType::New(); + checkerWriter->SetFileName(m_ArgsInfo.checker_after_arg); + checkerWriter->SetInput( checkerFilter->GetOutput() ); + checkerWriter->Update(); + } + + + //============================================================================ + // Checker before? + //============================================================================ + if (m_ArgsInfo.checker_before_given) { + //To display correctly the checkerboard image, the intensities must lie in the same range (normalized) + //We write the image in the internal image type + //We pass the normalized images to the checker filter + typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType; + typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New(); + + checkerFilter->SetInput1(inputFixedImage); + checkerFilter->SetInput2(inputMovingImage); + typedef itk::ImageFileWriter< InternalImageType > InternalWriterType; + typename InternalWriterType::Pointer checkerWriter = InternalWriterType::New(); + checkerWriter->SetFileName(m_ArgsInfo.checker_before_arg); + checkerWriter->SetInput( checkerFilter->GetOutput() ); + checkerWriter->Update(); + } + + + //============================================================================ + // Difference After? + //============================================================================ + if (m_ArgsInfo.after_given) { + typedef itk::SubtractImageFilter< InternalImageType, FixedImageType,FixedImageType > DifferenceImageFilterType; + typename DifferenceImageFilterType::Pointer differenceAfterFilter= DifferenceImageFilterType::New(); + + differenceAfterFilter->SetInput1(fixedImage); + differenceAfterFilter->SetInput2(resampler->GetOutput()); + + // Prepare a writer to write the difference image + typedef itk::ImageFileWriter< FixedImageType > WriterType; + typename WriterType::Pointer differenceAfterWriter = WriterType::New(); + differenceAfterWriter->SetFileName(m_ArgsInfo.after_arg ); + differenceAfterWriter->SetInput( differenceAfterFilter->GetOutput() ); + differenceAfterWriter->Update(); + } +// } + + //============================================================================ + // Difference Before? + //============================================================================ + if (m_ArgsInfo.before_given) { + typedef itk::CastImageFilter< InternalImageType,FixedImageType > CastFilterType; + typename CastFilterType::Pointer caster = CastFilterType::New(); + caster->SetInput( movingImage ); + + typedef itk::SubtractImageFilter< InternalImageType, FixedImageType, FixedImageType > DifferenceImageFilterType; + typename DifferenceImageFilterType::Pointer differenceBeforeFilter= DifferenceImageFilterType::New(); + + + differenceBeforeFilter->SetInput1(fixedImage); + differenceBeforeFilter->SetInput2(caster->GetOutput()); + + // Prepare a writer to write the difference image + typedef itk::ImageFileWriter< FixedImageType > WriterType; + typename WriterType::Pointer differenceBeforeWriter = WriterType::New(); + differenceBeforeWriter->SetFileName(m_ArgsInfo.before_arg); + differenceBeforeWriter->SetInput( differenceBeforeFilter->GetOutput() ); + differenceBeforeWriter->Update(); + } + +} } //====================================================================