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://oncora1.lyon.fnclcc.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_TXX
20 #define CLITKAFFINEREGISTRATIONGENERICFILTER_TXX
25 //==============================================================================
26 //Creating an observer class that allows us to monitor the registration
27 //================================================================================
28 class CommandIterationUpdate : public itk::Command
31 typedef CommandIterationUpdate Self;
32 typedef itk::Command Superclass;
33 typedef itk::SmartPointer<Self> Pointer;
36 CommandIterationUpdate() {};
38 typedef clitk::GenericOptimizer<args_info_clitkAffineRegistration> OptimizerType;
39 typedef const OptimizerType * OptimizerPointer;
41 // Set the generic optimizer
42 void SetOptimizer(OptimizerPointer o) {
47 void Execute(itk::Object *caller, const itk::EventObject & event) {
48 Execute( (const itk::Object *)caller, event);
51 void Execute(const itk::Object * object, const itk::EventObject & event) {
52 if( !(itk::IterationEvent().CheckEvent( &event )) ) {
56 m_Optimizer->OutputIterationInfo();
59 OptimizerPointer m_Optimizer;
63 //==============================================================================
64 //Creating an observer class that allows us to change parameters at subsequent levels
65 //==============================================================================
66 template <typename TRegistration, class args_info_clitkAffineRegistration>
67 class RegistrationInterfaceCommand : public itk::Command
70 typedef RegistrationInterfaceCommand Self;
71 typedef itk::Command Superclass;
72 typedef itk::SmartPointer<Self> Pointer;
75 RegistrationInterfaceCommand() {};
79 typedef TRegistration RegistrationType;
80 typedef RegistrationType * RegistrationPointer;
83 typedef typename RegistrationType::FixedImageType InternalImageType;
84 typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
86 // Two arguments are passed to the Execute() method: the first
87 // is the pointer to the object which invoked the event and the
88 // second is the event that was invoked.
89 void Execute(itk::Object * object, const itk::EventObject & event) {
90 if( !(itk::IterationEvent().CheckEvent( &event )) ) {
95 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
96 unsigned int numberOfLevels=registration->GetNumberOfLevels();
97 unsigned int currentLevel=registration->GetCurrentLevel()+1;
100 std::cout<<std::endl;
101 std::cout<<"========================================"<<std::endl;
102 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
103 std::cout<<"========================================"<<std::endl;
104 std::cout<<std::endl;
107 if (currentLevel>1) {
108 // Reinitialize the metric (!= number of samples)
109 typename GenericMetricType::Pointer genericMetric= GenericMetricType::New();
110 genericMetric->SetArgsInfo(m_ArgsInfo);
111 genericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
112 if (m_ArgsInfo.referenceMask_given) genericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
113 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
114 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
115 registration->SetMetric(metric);
119 void Execute(const itk::Object * , const itk::EventObject & ) {
123 void SetArgsInfo(args_info_clitkAffineRegistration a) {
126 args_info_clitkAffineRegistration m_ArgsInfo;
130 //==============================================================================
131 // Update with the number of dimensions
132 //==============================================================================
133 template <unsigned int Dimension>
134 void AffineRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
136 if (m_Verbose) std::cout << "Images were detected to be "<< Dimension << "D and " << PixelType << "..." << std::endl;
137 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
138 UpdateWithDimAndPixelType<Dimension, float>();
144 //==============================================================================
145 // Update with the number of dimensions and pixeltype
146 //==============================================================================
147 template <unsigned int Dimension, class PixelType>
149 AffineRegistrationGenericFilter::UpdateWithDimAndPixelType()
151 //=============================================================================
153 //=============================================================================
154 bool threadsGiven=m_ArgsInfo.threads_given;
155 int threads=m_ArgsInfo.threads_arg;
157 //Coordinate Representation
158 typedef double TCoordRep;
160 //The pixeltype of the fixed image will be used for output
161 typedef itk::Image< PixelType, Dimension > FixedImageType;
163 //Whatever the pixel type, internally we work with an image represented in float
164 typedef PixelType InternalPixelType;
165 typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
167 //Read in the reference/fixed image
168 typedef itk::ImageFileReader< InternalImageType > ReaderType;
169 typename ReaderType::Pointer fixedImageReader = ReaderType::New();
170 fixedImageReader->SetFileName( m_ArgsInfo.reference_arg);
173 //Read in the object/moving image
174 typename ReaderType::Pointer movingImageReader = ReaderType::New();
175 movingImageReader->SetFileName( m_ArgsInfo.target_arg );
176 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
177 fixedImageReader->Update();
178 movingImageReader->Update();
180 if (m_Verbose) std::cout << "Reading images... " << std::endl;
182 //we connect pointers to these internal images
183 typename InternalImageType::Pointer fixedImage= fixedImageReader->GetOutput();
184 typename InternalImageType::Pointer movingImage= movingImageReader->GetOutput();
186 //We keep the images used for input for possible output
187 typename InternalImageType::Pointer inputFixedImage= fixedImageReader->GetOutput();
188 typename InternalImageType::Pointer inputMovingImage= movingImageReader->GetOutput();
191 //============================================================================
193 //============================================================================
195 //If given, the intensities of both images are first normalized to a zero mean and SD of 1
196 // (usefull for MI, not necessary for Mattes' MI but performed anyway for the ouput)
197 if ( m_ArgsInfo.normalize_flag ) {
198 typedef itk::NormalizeImageFilter< InternalImageType,InternalImageType > NormalizeFilterType;
200 typename NormalizeFilterType::Pointer fixedNormalizeFilter = NormalizeFilterType::New();
201 typename NormalizeFilterType::Pointer movingNormalizeFilter = NormalizeFilterType::New();
203 fixedNormalizeFilter->SetInput( fixedImage );
204 movingNormalizeFilter->SetInput( movingImage );
206 fixedNormalizeFilter->Update();
207 movingNormalizeFilter->Update();
209 //We keep the images used for input for possible output
210 inputFixedImage= fixedNormalizeFilter->GetOutput();
211 inputMovingImage= movingNormalizeFilter->GetOutput();
213 //the pointers are reconnected for further output
214 fixedImage=fixedNormalizeFilter->GetOutput();
215 movingImage=movingNormalizeFilter->GetOutput();
217 if (m_Verbose) std::cout << "Normalizing image intensities to zero mean and SD of 1..." << std::endl;
221 //If given, the images are blurred before processing
222 if ( m_ArgsInfo.blur_arg!= 0) {
223 typedef itk::DiscreteGaussianImageFilter<InternalImageType,InternalImageType> GaussianFilterType;
224 typename GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
225 typename GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
226 fixedSmoother->SetVariance( m_ArgsInfo.blur_arg );
227 movingSmoother->SetVariance(m_ArgsInfo.blur_arg );
229 fixedSmoother->SetInput( fixedImage );
230 movingSmoother->SetInput( movingImage );
232 fixedSmoother->Update();
233 movingSmoother->Update();
235 fixedImage=fixedSmoother->GetOutput();
236 movingImage=movingSmoother->GetOutput();
238 if (m_Verbose) std::cout << "Blurring images with a Gaussian with standard deviation of " << m_ArgsInfo.blur_arg <<"..." << std::endl;
242 //============================================================================
243 // Setting up the moving image in a reference system
244 //============================================================================
245 const itk::Vector<double, Dimension> movingResolution = movingImage->GetSpacing();
246 typename InternalImageType::RegionType movingRegion = movingImage->GetLargestPossibleRegion();
247 typename InternalImageType::RegionType::SizeType movingSize = movingRegion.GetSize();
249 // Print the parameters of the moving image
251 std::cout << "Object or Moving image:"<<std::endl;
252 std::cout << "Size: " << movingSize[0] << ", " << movingSize[1];
253 if (Dimension==3) std::cout<<", " << movingSize[2];
254 std::cout << std::endl;
256 std::cout<< "Resolution: "<< movingResolution[0] << ", " << movingResolution[1];
257 if (Dimension==3) std::cout<< ", " << movingResolution[2];
258 std::cout << std::endl;
262 //============================================================================
263 // Setting up the fixed image in a reference system
264 //============================================================================
265 const itk::Vector<double, Dimension> fixedResolution = fixedImage->GetSpacing();
266 typename InternalImageType::RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
267 typename InternalImageType::RegionType::SizeType fixedSize = fixedRegion.GetSize();
269 // Print the parameters of the moving image and the transform
271 std::cout << "Target or Moving image:"<<std::endl;
272 std::cout << "Size: " << fixedSize[0] << ", " << fixedSize[1];
273 if (Dimension==3) std::cout<<", " << fixedSize[2];
274 std::cout << std::endl;
276 std::cout<< "Resolution: "<< fixedResolution[0] << ", " << fixedResolution[1];
277 if (Dimension==3) std::cout<< ", " << fixedResolution[2];
278 std::cout << std::endl;
283 //===========================================================================
284 // If given, we connect a mask to reference or target
285 //============================================================================
286 typedef itk::ImageMaskSpatialObject< Dimension > MaskType;
287 typename MaskType::Pointer fixedMask=NULL;
288 if (m_ArgsInfo.referenceMask_given) {
289 fixedMask= MaskType::New();
290 typedef itk::Image< unsigned char, Dimension > ImageMaskType;
291 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
292 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
293 maskReader->SetFileName(m_ArgsInfo.referenceMask_arg);
295 maskReader->Update();
296 } catch( itk::ExceptionObject & err ) {
297 std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
298 std::cerr << err << std::endl;
301 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
303 // Set the image to the spatialObject
304 fixedMask->SetImage( maskReader->GetOutput() );
307 typedef itk::ImageMaskSpatialObject< Dimension > MaskType;
308 typename MaskType::Pointer movingMask=NULL;
309 if (m_ArgsInfo.targetMask_given) {
310 movingMask= MaskType::New();
311 typedef itk::Image< unsigned char, Dimension > ImageMaskType;
312 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
313 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
314 maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
316 maskReader->Update();
317 } catch( itk::ExceptionObject & err ) {
318 std::cerr << "ExceptionObject caught !" << std::endl;
319 std::cerr << err << std::endl;
321 if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
323 movingMask->SetImage( maskReader->GetOutput() );
327 //============================================================================
328 // The image pyramids
329 //============================================================================
330 typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType > FixedImagePyramidType;
331 typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType > MovingImagePyramidType;
332 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
333 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
334 fixedImagePyramid->SetUseShrinkImageFilter(false);
335 fixedImagePyramid->SetInput(fixedImage);
336 fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
337 movingImagePyramid->SetUseShrinkImageFilter(false);
338 movingImagePyramid->SetInput(movingImage);
339 movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
340 if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
341 fixedImagePyramid->Update();
342 movingImagePyramid->Update();
346 //============================================================================
347 // We retrieve the type of metric from the command line
348 //============================================================================
349 typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
350 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
351 genericMetric->SetArgsInfo(m_ArgsInfo);
352 genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
353 if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
354 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
355 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
356 if (movingMask) metric->SetMovingImageMask(movingMask);
358 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
359 if (threadsGiven) metric->SetNumberOfThreads( threads );
361 if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
364 //============================================================================
365 // Initialize using image moments.
366 //============================================================================
367 if (m_ArgsInfo.moment_flag) {
368 typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
369 typename CalculatorType::Pointer fixedCalculator= CalculatorType::New();
371 typename InternalImageType::Pointer fixedThresh;
372 if (m_ArgsInfo.intThreshold_given) {
373 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
374 typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
375 thresholder->SetInput(fixedImage);
376 thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
377 thresholder->Update();
378 fixedThresh=thresholder->GetOutput();
379 } else fixedThresh=fixedImage;
381 fixedCalculator->SetImage(fixedThresh);
382 fixedCalculator->Compute();
383 Vector<double, Dimension> fixedCenter=fixedCalculator->GetCenterOfGravity();
384 if(m_Verbose)std::cout<<"The fixed center of gravity is "<<fixedCenter<<"..."<<std::endl;
386 typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
387 typename CalculatorType::Pointer movingCalculator= CalculatorType::New();
389 typename InternalImageType::Pointer movingThresh;
390 if (m_ArgsInfo.intThreshold_given) {
391 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
392 typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
393 thresholder->SetInput(movingImage);
394 thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
395 thresholder->Update();
396 movingThresh=thresholder->GetOutput();
397 } else movingThresh=movingImage;
399 movingCalculator->SetImage(movingThresh);
400 movingCalculator->Compute();
401 Vector<double, Dimension> movingCenter=movingCalculator->GetCenterOfGravity();
402 if(m_Verbose)std::cout<<"The moving center of gravity is "<<movingCenter<<"..."<<std::endl;
404 Vector<double, Dimension> shift= movingCenter-fixedCenter;
405 if(m_Verbose)std::cout<<"The initial shift applied is "<<shift<<"..."<<std::endl;
407 m_ArgsInfo.transX_arg= shift [0];
408 m_ArgsInfo.transY_arg= shift [1];
409 if (Dimension==3) m_ArgsInfo.transZ_arg=shift [2];
412 //============================================================================
414 //============================================================================
415 typedef clitk::GenericAffineTransform<args_info_clitkAffineRegistration, TCoordRep, Dimension > GenericAffineTransformType;
416 typename GenericAffineTransformType::Pointer genericAffineTransform = GenericAffineTransformType::New();
417 genericAffineTransform->SetArgsInfo(m_ArgsInfo);
418 typedef itk::Transform< double, Dimension, Dimension > TransformType;
419 typename TransformType::Pointer transform = genericAffineTransform->GetTransform();
422 //=======================================================
424 //=======================================================
425 typedef clitk::GenericInterpolator<args_info_clitkAffineRegistration, InternalImageType,TCoordRep > GenericInterpolatorType;
426 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
427 genericInterpolator->SetArgsInfo(m_ArgsInfo);
428 typedef itk::InterpolateImageFunction< InternalImageType, TCoordRep > InterpolatorType;
429 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
432 //============================================================================
434 //============================================================================
435 typedef clitk::GenericOptimizer<args_info_clitkAffineRegistration> GenericOptimizerType;
436 unsigned int nParam = transform->GetNumberOfParameters();
437 GenericOptimizerType::Pointer genericOptimizer=GenericOptimizerType::New();
438 genericOptimizer->SetArgsInfo(m_ArgsInfo);
439 genericOptimizer->SetOutputIteration(m_Verbose);
440 genericOptimizer->SetOutputPosition(m_Verbose);
441 genericOptimizer->SetOutputValue(m_Verbose);
442 genericOptimizer->SetOutputGradient(m_ArgsInfo.gradient_flag);
443 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
444 genericOptimizer->SetNumberOfParameters(nParam);
445 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
446 OptimizerType::Pointer optimizer=genericOptimizer->GetOptimizerPointer();
449 itk::Optimizer::ScalesType scales( nParam );
450 for (unsigned int i=nParam-Dimension; i<nParam; i++) //Translations
451 scales[i] = m_ArgsInfo.tWeight_arg;
452 for (unsigned int i=0; i<nParam-Dimension; i++) //Rest
453 scales[i] = m_ArgsInfo.rWeight_arg*180./M_PI;
454 optimizer->SetScales(scales);
456 //============================================================================
457 // Multiresolution registration
458 //============================================================================
459 typedef itk::MultiResolutionImageRegistrationMethod< InternalImageType,InternalImageType > RegistrationType;
460 typename RegistrationType::Pointer registration = RegistrationType::New();
461 registration->SetFixedImage( fixedImage );
462 registration->SetFixedImageRegion(fixedImage->GetLargestPossibleRegion());
463 registration->SetMovingImage( movingImage );
464 registration->SetFixedImagePyramid( fixedImagePyramid );
465 registration->SetMovingImagePyramid( movingImagePyramid );
466 registration->SetTransform( transform );
467 registration->SetInitialTransformParameters( transform->GetParameters() );
468 registration->SetInterpolator( interpolator );
469 registration->SetMetric(metric);
470 registration->SetOptimizer(optimizer);
471 registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
472 if (m_Verbose) std::cout << "Setting "<< m_ArgsInfo.levels_arg <<" resolution levels..." << std::endl;
473 if (m_Verbose) std::cout << "Initial Transform: "<< registration->GetInitialTransformParameters()<<std::endl;
475 //============================================================================
476 // Connecting the commander to the registration to monitor it
477 //============================================================================
480 // Output iteration info
481 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
482 observer->SetOptimizer(genericOptimizer);
483 optimizer->AddObserver( itk::IterationEvent(), observer );
487 typedef RegistrationInterfaceCommand<RegistrationType, args_info_clitkAffineRegistration> CommandType;
488 typename CommandType::Pointer command = CommandType::New();
489 command->SetArgsInfo(m_ArgsInfo);
490 registration->AddObserver( itk::IterationEvent(), command );
495 //============================================================================
496 // Finally we can start the registration with the given amount of multiresolution levels
497 //============================================================================
498 if (m_Verbose) std::cout << "Starting the registration now..." << std::endl;
501 registration->StartRegistration();
502 } catch( itk::ExceptionObject & err ) {
503 std::cerr << "ExceptionObject caught !" << std::endl;
504 std::cerr << err << std::endl;
508 //============================================================================
509 // Processing the result of the registration
510 //============================================================================
511 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
512 std::cout<< "Result : " <<std::setprecision(12)<<std::endl;
514 for (unsigned int i=nParam-Dimension; i<nParam; i++) //Translations
515 std::cout << " Translation " << i << " = " << finalParameters[i];
516 for (unsigned int i=0; i<nParam-Dimension; i++) //Rest
517 std::cout << " Other parameter " << i << " = " << finalParameters[i];
520 itk::Matrix<double,Dimension+1,Dimension+1> matrix;
521 if(m_ArgsInfo.transform_arg == 3) {
522 for(unsigned int i=0; i<Dimension; i++) {
523 matrix[i][3] = finalParameters[nParam-Dimension+i];
524 for(unsigned int j=0; j<Dimension; j++) {
525 matrix[i][j] = finalParameters[i*3+j];
530 // matrix = clitk::GetBackwardAffineMatrix3D(finalParameters);
533 std::cout << " Affine transform matrix =" << std::endl;
534 std::cout << matrix <<std::setprecision(6)<< std::endl;
536 // Write matrix to a file
537 if(m_ArgsInfo.matrix_given) {
539 mFile.open(m_ArgsInfo.matrix_arg);
540 mFile<<std::setprecision(12)<<matrix<< std::setprecision(6)<<std::endl;
544 //============================================================================
545 // Prepare the resampling filter in order to transform the moving image.
546 //============================================================================
547 if (m_ArgsInfo.output_given || m_ArgsInfo.checker_after_given || m_ArgsInfo.after_given ) {
548 transform->SetParameters( finalParameters );
549 typedef itk::ResampleImageFilter< InternalImageType,InternalImageType > ResampleFilterType;
550 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
552 resampler->SetTransform( transform );
553 resampler->SetInput( movingImageReader->GetOutput() );
554 resampler->SetSize( fixedImageReader ->GetOutput()->GetLargestPossibleRegion().GetSize() );
555 resampler->SetOutputOrigin( fixedImageReader ->GetOutput()->GetOrigin() );
556 resampler->SetOutputSpacing( fixedImageReader ->GetOutput()->GetSpacing() );
557 resampler->SetDefaultPixelValue( 0 );
560 if (m_ArgsInfo.output_given) {
561 //We write an output in the same pixeltype then the input
562 typedef itk::ImageFileWriter< FixedImageType > WriterType;
563 typename WriterType::Pointer outputWriter = WriterType::New();
564 outputWriter->SetFileName(m_ArgsInfo.output_arg );
565 outputWriter->SetInput( resampler->GetOutput() );
566 outputWriter->Update();
570 //============================================================================
572 //============================================================================
573 if (m_ArgsInfo.checker_after_given) {
574 //To display correctly the checkerboard image, the intensities must lie in the same range (normalized)
575 //We write the image in the internal image type
576 typedef itk::ResampleImageFilter< InternalImageType,InternalImageType > ResampleFilterType;
577 typename ResampleFilterType::Pointer internalResampler = ResampleFilterType::New();
578 internalResampler->SetTransform( transform );
579 internalResampler->SetInput( inputMovingImage );
580 internalResampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
581 internalResampler->SetOutputOrigin( fixedImage->GetOrigin() );
582 internalResampler->SetOutputSpacing( fixedImage->GetSpacing() );
583 internalResampler->SetDefaultPixelValue( 0 );
585 //We pass the normalized images to the checker filter
586 typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType;
587 typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New();
589 checkerFilter->SetInput1(inputFixedImage);
590 checkerFilter->SetInput2(internalResampler->GetOutput());
591 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
592 typename InternalWriterType::Pointer checkerWriter = InternalWriterType::New();
593 checkerWriter->SetFileName(m_ArgsInfo.checker_after_arg);
594 checkerWriter->SetInput( checkerFilter->GetOutput() );
595 checkerWriter->Update();
599 //============================================================================
601 //============================================================================
602 if (m_ArgsInfo.checker_before_given) {
603 //To display correctly the checkerboard image, the intensities must lie in the same range (normalized)
604 //We write the image in the internal image type
605 //We pass the normalized images to the checker filter
606 typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType;
607 typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New();
609 checkerFilter->SetInput1(inputFixedImage);
610 checkerFilter->SetInput2(inputMovingImage);
611 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
612 typename InternalWriterType::Pointer checkerWriter = InternalWriterType::New();
613 checkerWriter->SetFileName(m_ArgsInfo.checker_before_arg);
614 checkerWriter->SetInput( checkerFilter->GetOutput() );
615 checkerWriter->Update();
619 //============================================================================
621 //============================================================================
622 if (m_ArgsInfo.after_given) {
623 typedef itk::SubtractImageFilter< InternalImageType, FixedImageType,FixedImageType > DifferenceImageFilterType;
624 typename DifferenceImageFilterType::Pointer differenceAfterFilter= DifferenceImageFilterType::New();
626 differenceAfterFilter->SetInput1(fixedImageReader ->GetOutput());
627 differenceAfterFilter->SetInput2(resampler->GetOutput());
629 // Prepare a writer to write the difference image
630 typedef itk::ImageFileWriter< FixedImageType > WriterType;
631 typename WriterType::Pointer differenceAfterWriter = WriterType::New();
632 differenceAfterWriter->SetFileName(m_ArgsInfo.after_arg );
633 differenceAfterWriter->SetInput( differenceAfterFilter->GetOutput() );
634 differenceAfterWriter->Update();
638 //============================================================================
639 // Difference Before?
640 //============================================================================
641 if (m_ArgsInfo.before_given) {
642 typedef itk::CastImageFilter< InternalImageType,FixedImageType > CastFilterType;
643 typename CastFilterType::Pointer caster = CastFilterType::New();
644 caster->SetInput( movingImageReader->GetOutput() );
646 typedef itk::SubtractImageFilter< InternalImageType, FixedImageType, FixedImageType > DifferenceImageFilterType;
647 typename DifferenceImageFilterType::Pointer differenceBeforeFilter= DifferenceImageFilterType::New();
650 differenceBeforeFilter->SetInput1(fixedImageReader ->GetOutput());
651 differenceBeforeFilter->SetInput2(caster->GetOutput());
653 // Prepare a writer to write the difference image
654 typedef itk::ImageFileWriter< FixedImageType > WriterType;
655 typename WriterType::Pointer differenceBeforeWriter = WriterType::New();
656 differenceBeforeWriter->SetFileName(m_ArgsInfo.before_arg);
657 differenceBeforeWriter->SetInput( differenceBeforeFilter->GetOutput() );
658 differenceBeforeWriter->Update();
663 #endif //#define CLITKAFFINEREGISTRATIONGENERICFILTER_TXX