1 /*=========================================================================
2 Program: vv http://www.creatis.insa-lyon.fr/rio/vv
5 - University of LYON http://www.universite-lyon.fr/
6 - Léon Bérard cancer center http://www.centreleonberard.fr
7 - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr
9 This software is distributed WITHOUT ANY WARRANTY; without even
10 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11 PURPOSE. See the copyright notices for more information.
13 It is distributed under dual licence
15 - BSD See included LICENSE.txt file
16 - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
19 #ifndef CLITKAFFINEREGISTRATIONGENERICFILTER_CXX
20 #define CLITKAFFINEREGISTRATIONGENERICFILTER_CXX
22 #include "clitkAffineRegistrationGenericFilter.h"
26 //==============================================================================
27 //Creating an observer class that allows us to monitor the registration
28 //================================================================================
29 class CommandIterationUpdate : public itk::Command
32 typedef CommandIterationUpdate Self;
33 typedef itk::Command Superclass;
34 typedef itk::SmartPointer<Self> Pointer;
37 CommandIterationUpdate() {};
39 typedef clitk::GenericOptimizer<args_info_clitkAffineRegistration> OptimizerType;
40 typedef const OptimizerType * OptimizerPointer;
42 // Set the generic optimizer
43 void SetOptimizer(OptimizerPointer o) {
48 void Execute(itk::Object *caller, const itk::EventObject & event) {
49 Execute( (const itk::Object *)caller, event);
52 void Execute(const itk::Object * object, const itk::EventObject & event) {
53 if ( !(itk::IterationEvent().CheckEvent( &event )) ) {
57 m_Optimizer->OutputIterationInfo();
60 OptimizerPointer m_Optimizer;
62 //==================================================================================================================================//
64 //===================================================================================================================================//
66 AffineRegistrationGenericFilter::AffineRegistrationGenericFilter():
67 ImageToImageGenericFilter<Self>("Register")
70 InitializeImageType<2>();
71 InitializeImageType<3>();
74 //==========================================================================================================//
75 //============================================================================================================//
76 //--------------------------------------------------------------------
77 template<unsigned int Dim>
78 void AffineRegistrationGenericFilter::InitializeImageType()
80 ADD_DEFAULT_IMAGE_TYPES(Dim);
82 //--------------------------------------------------------------------
85 //==============================================================================
86 //Creating an observer class that allows us to change parameters at subsequent levels
87 //==============================================================================
88 template <typename TRegistration>
89 class RegistrationInterfaceCommand : public itk::Command
92 typedef RegistrationInterfaceCommand Self;
93 typedef itk::Command Superclass;
94 typedef itk::SmartPointer<Self> Pointer;
97 RegistrationInterfaceCommand() {};
101 typedef TRegistration RegistrationType;
102 typedef RegistrationType * RegistrationPointer;
105 typedef typename RegistrationType::FixedImageType InternalImageType;
106 typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
108 // Two arguments are passed to the Execute() method: the first
109 // is the pointer to the object which invoked the event and the
110 // second is the event that was invoked.
111 void Execute(itk::Object * object, const itk::EventObject & event) {
112 if ( !(itk::IterationEvent().CheckEvent( &event )) ) {
117 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
118 unsigned int numberOfLevels=registration->GetNumberOfLevels();
119 unsigned int currentLevel=registration->GetCurrentLevel()+1;
122 std::cout<<std::endl;
123 std::cout<<"========================================"<<std::endl;
124 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
125 std::cout<<"========================================"<<std::endl;
126 std::cout<<std::endl;
129 if (currentLevel>1) {
130 // Reinitialize the metric (!= number of samples)
131 typename GenericMetricType::Pointer genericMetric= GenericMetricType::New();
132 genericMetric->SetArgsInfo(m_ArgsInfo);
133 genericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
134 if (m_ArgsInfo.referenceMask_given) genericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
135 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
136 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
137 registration->SetMetric(metric);
141 void Execute(const itk::Object * , const itk::EventObject & ) {
145 void SetArgsInfo(args_info_clitkAffineRegistration a) {
148 args_info_clitkAffineRegistration m_ArgsInfo;
151 //==============================================================================================//
153 //==============================================================================================//
154 void AffineRegistrationGenericFilter::SetArgsInfo(const args_info_clitkAffineRegistration & a)
157 if (m_ArgsInfo.reference_given) AddInputFilename(m_ArgsInfo.reference_arg);
159 if (m_ArgsInfo.target_given) {
160 AddInputFilename(m_ArgsInfo.target_arg);
163 if (m_ArgsInfo.output_given) SetOutputFilename(m_ArgsInfo.output_arg);
165 //==============================================================================
166 // Update with the number of dimensions and pixeltype
167 //==============================================================================
168 template<class InputImageType>
169 void AffineRegistrationGenericFilter::UpdateWithInputImageType()
171 //=============================================================================
173 //=============================================================================
175 typedef typename InputImageType::PixelType PixelType;
176 //typedef typename InputImageType::ImageDimension Dimension;
179 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
180 bool threadsGiven=m_ArgsInfo.threads_given;
181 int threads=m_ArgsInfo.threads_arg;
184 //Coordinate Representation
185 typedef double TCoordRep;
188 typename InputImageType::Pointer fixedImage = this->template GetInput<InputImageType>(0);
190 typename InputImageType::Pointer inputFixedImage = this->template GetInput<InputImageType>(0);
193 typename InputImageType::Pointer movingImage = this->template GetInput<InputImageType>(1);
195 typename InputImageType::Pointer inputMovingImage = this->template GetInput<InputImageType>(1);
199 //The pixeltype of the fixed image will be used for output
200 typedef itk::Image< PixelType, InputImageType::ImageDimension > FixedImageType;
202 //Whatever the pixel type, internally we work with an image represented in float
203 typedef typename InputImageType::PixelType InternalPixelType;
204 typedef itk::Image< PixelType, InputImageType::ImageDimension > InternalImageType;
207 //Read in the reference/fixed image
208 // typedef itk::ImageFileReader< InternalImageType > ReaderType;
209 // typename ReaderType::Pointer fixedImageReader = ReaderType::New();
210 // fixedImageReader->SetFileName( m_ArgsInfo.reference_arg);
213 //Read in the object/moving image
214 // typename ReaderType::Pointer movingImageReader = ReaderType::New();
215 // movingImageReader->SetFileName( m_ArgsInfo.target_arg );
216 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
217 // fixedImageReader->Update();
218 // movingImageReader->Update();
220 if (m_Verbose) std::cout << "Reading images... " << std::endl;
222 //we connect pointers to these internal images
223 // typedef typename fixedImageReader fixedImage;
224 // typedef typename movingImageReader movingImage;
226 //We keep the images used for input for possible output
227 // typedef typename fixedImageReader inputFixedImage;
228 // typedef typename movingImageReader inputMovingImage;
231 //============================================================================
233 //============================================================================
235 //If given, the intensities of both images are first normalized to a zero mean and SD of 1
236 // (usefull for MI, not necessary for Mattes' MI but performed anyway for the ouput)
237 if ( m_ArgsInfo.normalize_flag ) {
238 typedef itk::NormalizeImageFilter< InternalImageType,InternalImageType > NormalizeFilterType;
240 typename NormalizeFilterType::Pointer fixedNormalizeFilter = NormalizeFilterType::New();
241 typename NormalizeFilterType::Pointer movingNormalizeFilter = NormalizeFilterType::New();
243 fixedNormalizeFilter->SetInput( fixedImage );
244 movingNormalizeFilter->SetInput( movingImage );
246 fixedNormalizeFilter->Update();
247 movingNormalizeFilter->Update();
249 //We keep the images used for input for possible output
250 inputFixedImage= fixedNormalizeFilter->GetOutput();
251 inputMovingImage= movingNormalizeFilter->GetOutput();
253 //the pointers are reconnected for further output
254 fixedImage=fixedNormalizeFilter->GetOutput();
255 movingImage=movingNormalizeFilter->GetOutput();
257 if (m_Verbose) std::cout << "Normalizing image intensities to zero mean and SD of 1..." << std::endl;
261 //If given, the images are blurred before processing
262 if ( m_ArgsInfo.blur_arg!= 0) {
263 typedef itk::DiscreteGaussianImageFilter<InternalImageType,InternalImageType> GaussianFilterType;
264 typename GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
265 typename GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
266 fixedSmoother->SetVariance( m_ArgsInfo.blur_arg );
267 movingSmoother->SetVariance(m_ArgsInfo.blur_arg );
269 fixedSmoother->SetInput( fixedImage );
270 movingSmoother->SetInput( movingImage );
272 fixedSmoother->Update();
273 movingSmoother->Update();
275 fixedImage=fixedSmoother->GetOutput();
276 movingImage=movingSmoother->GetOutput();
278 if (m_Verbose) std::cout << "Blurring images with a Gaussian with standard deviation of " << m_ArgsInfo.blur_arg <<"..." << std::endl;
282 //============================================================================
283 // Setting up the moving image in a reference system
284 //============================================================================
285 const itk::Vector<double, InputImageType::ImageDimension> movingResolution = movingImage->GetSpacing();
286 typename InternalImageType::RegionType movingRegion = movingImage->GetLargestPossibleRegion();
287 typename InternalImageType::RegionType::SizeType movingSize = movingRegion.GetSize();
289 // Print the parameters of the moving image
291 std::cout << "Object or Moving image:"<<std::endl;
292 std::cout << "Size: " << movingSize[0] << ", " << movingSize[1];
293 if (InputImageType::ImageDimension==3) std::cout<<", " << movingSize[2];
294 std::cout << std::endl;
296 std::cout<< "Resolution: "<< movingResolution[0] << ", " << movingResolution[1];
297 if (InputImageType::ImageDimension==3) std::cout<< ", " << movingResolution[2];
298 std::cout << std::endl;
302 //============================================================================
303 // Setting up the fixed image in a reference system
304 //============================================================================
305 const itk::Vector<double, InputImageType::ImageDimension> fixedResolution = fixedImage->GetSpacing();
306 typename InternalImageType::RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
307 typename InternalImageType::RegionType::SizeType fixedSize = fixedRegion.GetSize();
309 // Print the parameters of the moving image and the transform
311 std::cout << "Target or Moving image:"<<std::endl;
312 std::cout << "Size: " << fixedSize[0] << ", " << fixedSize[1];
313 if (InputImageType::ImageDimension==3) std::cout<<", " << fixedSize[2];
314 std::cout << std::endl;
316 std::cout<< "Resolution: "<< fixedResolution[0] << ", " << fixedResolution[1];
317 if (InputImageType::ImageDimension==3) std::cout<< ", " << fixedResolution[2];
318 std::cout << std::endl;
323 //===========================================================================
324 // If given, we connect a mask to reference or target
325 //============================================================================
326 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
327 typename MaskType::Pointer fixedMask=NULL;
328 if (m_ArgsInfo.referenceMask_given) {
329 fixedMask= MaskType::New();
330 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
331 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
332 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
333 maskReader->SetFileName(m_ArgsInfo.referenceMask_arg);
335 maskReader->Update();
336 } catch ( itk::ExceptionObject & err ) {
337 std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
338 std::cerr << err << std::endl;
341 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
343 // Set the image to the spatialObject
344 fixedMask->SetImage( maskReader->GetOutput() );
347 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
348 typename MaskType::Pointer movingMask=NULL;
349 if (m_ArgsInfo.targetMask_given) {
350 movingMask= MaskType::New();
351 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
352 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
353 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
354 maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
356 maskReader->Update();
357 } catch ( itk::ExceptionObject & err ) {
358 std::cerr << "ExceptionObject caught !" << std::endl;
359 std::cerr << err << std::endl;
361 if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
363 movingMask->SetImage( maskReader->GetOutput() );
367 //============================================================================
368 // The image pyramids
369 //============================================================================
370 typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType > FixedImagePyramidType;
371 typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType > MovingImagePyramidType;
372 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
373 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
374 fixedImagePyramid->SetUseShrinkImageFilter(false);
375 fixedImagePyramid->SetInput(fixedImage);
376 fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
377 movingImagePyramid->SetUseShrinkImageFilter(false);
378 movingImagePyramid->SetInput(movingImage);
379 movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
380 if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
382 fixedImagePyramid->Update();
383 movingImagePyramid->Update();
387 //============================================================================
388 // We retrieve the type of metric from the command line
389 //============================================================================
390 typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
391 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
392 genericMetric->SetArgsInfo(m_ArgsInfo);
393 genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
394 if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
395 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
396 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
397 if (movingMask) metric->SetMovingImageMask(movingMask);
399 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
400 if (threadsGiven) metric->SetNumberOfThreads( threads );
402 if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
405 //============================================================================
406 // Initialize using image moments.
407 //============================================================================
408 if (m_ArgsInfo.moment_flag) {
409 typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
410 typename CalculatorType::Pointer fixedCalculator= CalculatorType::New();
412 typename InternalImageType::Pointer fixedThresh;
413 if (m_ArgsInfo.intThreshold_given) {
414 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
415 typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
416 thresholder->SetInput(fixedImage);
417 thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
418 thresholder->Update();
419 fixedThresh=thresholder->GetOutput();
420 } else fixedThresh=fixedImage;
422 fixedCalculator->SetImage(fixedThresh);
423 fixedCalculator->Compute();
424 Vector<double, InputImageType::ImageDimension> fixedCenter=fixedCalculator->GetCenterOfGravity();
425 if (m_Verbose)std::cout<<"The fixed center of gravity is "<<fixedCenter<<"..."<<std::endl;
427 typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
428 typename CalculatorType::Pointer movingCalculator= CalculatorType::New();
430 typename InternalImageType::Pointer movingThresh;
431 if (m_ArgsInfo.intThreshold_given) {
432 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
433 typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
434 thresholder->SetInput(movingImage);
435 thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
436 thresholder->Update();
437 movingThresh=thresholder->GetOutput();
438 } else movingThresh=movingImage;
440 movingCalculator->SetImage(movingThresh);
441 movingCalculator->Compute();
442 Vector<double, InputImageType::ImageDimension> movingCenter=movingCalculator->GetCenterOfGravity();
443 if (m_Verbose)std::cout<<"The moving center of gravity is "<<movingCenter<<"..."<<std::endl;
445 Vector<double, InputImageType::ImageDimension> shift= movingCenter-fixedCenter;
446 if (m_Verbose)std::cout<<"The initial shift applied is "<<shift<<"..."<<std::endl;
448 m_ArgsInfo.transX_arg= shift [0];
449 m_ArgsInfo.transY_arg= shift [1];
450 if (InputImageType::ImageDimension==3) m_ArgsInfo.transZ_arg=shift [2];
453 //============================================================================
455 //============================================================================
456 typedef clitk::GenericAffineTransform<args_info_clitkAffineRegistration, TCoordRep, InputImageType::ImageDimension > GenericAffineTransformType;
457 typename GenericAffineTransformType::Pointer genericAffineTransform = GenericAffineTransformType::New();
458 genericAffineTransform->SetArgsInfo(m_ArgsInfo);
459 typedef itk::Transform< double, InputImageType::ImageDimension, InputImageType::ImageDimension > TransformType;
460 typename TransformType::Pointer transform = genericAffineTransform->GetTransform();
461 std::cout<<m_ArgsInfo.transform_arg<<std::endl;
464 //=======================================================
466 //=======================================================
467 std::cout<<"setting Interpolator..."<<std::endl;
468 typedef clitk::GenericInterpolator<args_info_clitkAffineRegistration, InternalImageType,TCoordRep > GenericInterpolatorType;
469 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
470 genericInterpolator->SetArgsInfo(m_ArgsInfo);
471 typedef itk::InterpolateImageFunction< InternalImageType, TCoordRep > InterpolatorType;
472 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
473 std::cout<<"end of interpolator"<<std::endl;
475 //============================================================================
477 //============================================================================
478 typedef clitk::GenericOptimizer<args_info_clitkAffineRegistration> GenericOptimizerType;
479 unsigned int nParam = transform->GetNumberOfParameters();
480 typename GenericOptimizerType::Pointer genericOptimizer=GenericOptimizerType::New();
481 genericOptimizer->SetArgsInfo(m_ArgsInfo);
482 genericOptimizer->SetOutputIteration(m_Verbose);
483 genericOptimizer->SetOutputPosition(m_Verbose);
484 genericOptimizer->SetOutputValue(m_Verbose);
485 genericOptimizer->SetOutputGradient(m_ArgsInfo.gradient_flag);
486 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
487 genericOptimizer->SetNumberOfParameters(nParam);
488 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
489 OptimizerType::Pointer optimizer=genericOptimizer->GetOptimizerPointer();
492 itk::Optimizer::ScalesType scales( nParam );
493 for (unsigned int i=nParam-InputImageType::ImageDimension; i<nParam; i++) //Translations
494 scales[i] = m_ArgsInfo.tWeight_arg;
495 for (unsigned int i=0; i<nParam-InputImageType::ImageDimension; i++) //Rest
496 scales[i] = m_ArgsInfo.rWeight_arg*180./M_PI;
497 optimizer->SetScales(scales);
498 //============================================================================
499 // Multiresolution registration
500 //============================================================================
501 std::cout<<"start MultiResolution..."<<std::endl;
502 typedef itk::MultiResolutionImageRegistrationMethod< InternalImageType,InternalImageType > RegistrationType;
503 typename RegistrationType::Pointer registration = RegistrationType::New();
504 registration->SetFixedImage( fixedImage );
505 registration->SetFixedImageRegion(fixedImage->GetLargestPossibleRegion());
506 registration->SetMovingImage( movingImage );
507 registration->SetFixedImagePyramid( fixedImagePyramid );
508 registration->SetMovingImagePyramid( movingImagePyramid );
509 registration->SetTransform( transform );
510 registration->SetInitialTransformParameters( transform->GetParameters() );
511 registration->SetInterpolator( interpolator );
512 registration->SetMetric(metric);
513 registration->SetOptimizer(optimizer);
514 registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
515 if (m_Verbose) std::cout << "Setting "<< m_ArgsInfo.levels_arg <<" resolution levels..." << std::endl;
516 if (m_Verbose) std::cout << "Initial Transform: "<< registration->GetInitialTransformParameters()<<std::endl;
518 //============================================================================
519 // Connecting the commander to the registration to monitor it
520 //============================================================================
523 // Output iteration info
524 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
525 observer->SetOptimizer(genericOptimizer);
526 optimizer->AddObserver( itk::IterationEvent(), observer );
530 typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
531 typename CommandType::Pointer command = CommandType::New();
532 command->SetArgsInfo(m_ArgsInfo);
533 registration->AddObserver( itk::IterationEvent(), command );
538 //============================================================================
539 // Finally we can start the registration with the given amount of multiresolution levels
540 //============================================================================
541 if (m_Verbose) std::cout << "Starting the registration now..." << std::endl;
544 #if ITK_VERSION_MAJOR < 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR <= 2)
545 registration->StartRegistration();
547 registration->Update();
549 } catch ( itk::ExceptionObject & err ) {
550 std::cerr << "ExceptionObject caught !" << std::endl;
551 std::cerr << err << std::endl;
555 //============================================================================
556 // Processing the result of the registration
557 //============================================================================
558 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
559 std::cout<< "Result : " <<std::setprecision(12)<<std::endl;
561 for (unsigned int i=nParam-InputImageType::ImageDimension; i<nParam; i++) //Translations
562 std::cout << " Translation " << i << " = " << finalParameters[i];
563 for (unsigned int i=0; i<nParam-InputImageType::ImageDimension; i++) //Rest
564 std::cout << " Other parameter " << i << " = " << finalParameters[i];
566 itk::Matrix<double,InputImageType::ImageDimension+1,InputImageType::ImageDimension+1> matrix;
567 if (m_ArgsInfo.transform_arg == 3) {
568 for (unsigned int i=0; i<InputImageType::ImageDimension; i++) {
569 matrix[i][3] = finalParameters[nParam-InputImageType::ImageDimension+i];
570 for (unsigned int j=0; j<InputImageType::ImageDimension; j++) {
571 matrix[i][j] = finalParameters[i*3+j];
576 matrix = clitk::GetBackwardAffineMatrix<InputImageType::ImageDimension>(finalParameters);
577 std::cout<<"outside GetBackWardAffineMatrix...."<<std::endl;
580 std::cout << " Affine transform matrix =" << std::endl;
581 std::cout << matrix <<std::setprecision(6)<< std::endl;
582 std::cout << " End of Registration" << std::endl;
583 // Write matrix to a file
584 if (m_ArgsInfo.matrix_given) {
586 mFile.open(m_ArgsInfo.matrix_arg);
587 mFile<<std::setprecision(12)<<matrix<< std::setprecision(6)<<std::endl;
591 //============================================================================
592 // Prepare the resampling filter in order to transform the moving image.
593 //============================================================================
594 // if (m_ArgsInfo.output_given || m_ArgsInfo.checker_after_given || m_ArgsInfo.after_given ) {
595 transform->SetParameters( finalParameters );
596 typedef itk::ResampleImageFilter< InternalImageType,InternalImageType > ResampleFilterType;
597 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
599 resampler->SetTransform( transform );
600 resampler->SetInput( movingImage );
601 resampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
602 resampler->SetOutputOrigin( fixedImage->GetOrigin() );
603 resampler->SetOutputSpacing( fixedImage->GetSpacing() );
604 resampler->SetDefaultPixelValue( 0 );
607 // if (m_ArgsInfo.output_given) {
608 //We write an output in the same pixeltype then the input
609 /*typedef itk::ImageFileWriter< FixedImageType > WriterType;
610 typename WriterType::Pointer outputWriter = WriterType::New();
611 outputWriter->SetFileName(m_ArgsInfo.output_arg );
612 outputWriter->SetInput( resampler->GetOutput() );
613 outputWriter->Update();*/
614 typedef InternalImageType OutputImageType;
615 typename OutputImageType::Pointer outputImage = resampler->GetOutput();
616 std::cout<<"Writing Output....."<<std::endl;
617 this->template SetNextOutput<OutputImageType>(outputImage);
621 //============================================================================
623 //============================================================================
624 if (m_ArgsInfo.checker_after_given) {
625 //To display correctly the checkerboard image, the intensities must lie in the same range (normalized)
626 //We write the image in the internal image type
627 typedef itk::ResampleImageFilter< InternalImageType,InternalImageType > ResampleFilterType;
628 typename ResampleFilterType::Pointer internalResampler = ResampleFilterType::New();
629 internalResampler->SetTransform( transform );
630 internalResampler->SetInput( inputMovingImage );
631 internalResampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
632 internalResampler->SetOutputOrigin( fixedImage->GetOrigin() );
633 internalResampler->SetOutputSpacing( fixedImage->GetSpacing() );
634 internalResampler->SetDefaultPixelValue( 0 );
636 //We pass the normalized images to the checker filter
637 typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType;
638 typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New();
640 checkerFilter->SetInput1(inputFixedImage);
641 checkerFilter->SetInput2(internalResampler->GetOutput());
642 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
643 typename InternalWriterType::Pointer checkerWriter = InternalWriterType::New();
644 checkerWriter->SetFileName(m_ArgsInfo.checker_after_arg);
645 checkerWriter->SetInput( checkerFilter->GetOutput() );
646 checkerWriter->Update();
650 //============================================================================
652 //============================================================================
653 if (m_ArgsInfo.checker_before_given) {
654 //To display correctly the checkerboard image, the intensities must lie in the same range (normalized)
655 //We write the image in the internal image type
656 //We pass the normalized images to the checker filter
657 typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType;
658 typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New();
660 checkerFilter->SetInput1(inputFixedImage);
661 checkerFilter->SetInput2(inputMovingImage);
662 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
663 typename InternalWriterType::Pointer checkerWriter = InternalWriterType::New();
664 checkerWriter->SetFileName(m_ArgsInfo.checker_before_arg);
665 checkerWriter->SetInput( checkerFilter->GetOutput() );
666 checkerWriter->Update();
670 //============================================================================
672 //============================================================================
673 if (m_ArgsInfo.after_given) {
674 typedef itk::SubtractImageFilter< InternalImageType, FixedImageType,FixedImageType > DifferenceImageFilterType;
675 typename DifferenceImageFilterType::Pointer differenceAfterFilter= DifferenceImageFilterType::New();
677 differenceAfterFilter->SetInput1(fixedImage);
678 differenceAfterFilter->SetInput2(resampler->GetOutput());
680 // Prepare a writer to write the difference image
681 typedef itk::ImageFileWriter< FixedImageType > WriterType;
682 typename WriterType::Pointer differenceAfterWriter = WriterType::New();
683 differenceAfterWriter->SetFileName(m_ArgsInfo.after_arg );
684 differenceAfterWriter->SetInput( differenceAfterFilter->GetOutput() );
685 differenceAfterWriter->Update();
689 //============================================================================
690 // Difference Before?
691 //============================================================================
692 if (m_ArgsInfo.before_given) {
693 typedef itk::CastImageFilter< InternalImageType,FixedImageType > CastFilterType;
694 typename CastFilterType::Pointer caster = CastFilterType::New();
695 caster->SetInput( movingImage );
697 typedef itk::SubtractImageFilter< InternalImageType, FixedImageType, FixedImageType > DifferenceImageFilterType;
698 typename DifferenceImageFilterType::Pointer differenceBeforeFilter= DifferenceImageFilterType::New();
701 differenceBeforeFilter->SetInput1(fixedImage);
702 differenceBeforeFilter->SetInput2(caster->GetOutput());
704 // Prepare a writer to write the difference image
705 typedef itk::ImageFileWriter< FixedImageType > WriterType;
706 typename WriterType::Pointer differenceBeforeWriter = WriterType::New();
707 differenceBeforeWriter->SetFileName(m_ArgsInfo.before_arg);
708 differenceBeforeWriter->SetInput( differenceBeforeFilter->GetOutput() );
709 differenceBeforeWriter->Update();
715 //====================================================================
717 #endif //#define CLITKAFFINEREGISTRATIONCGENERICFILTER_CXX