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;
61 //==================================================================================================================================//
63 //===================================================================================================================================//
65 template<class args_info_clitkAffineRegistration>
66 AffineRegistrationGenericFilter<args_info_clitkAffineRegistration>::AffineRegistrationGenericFilter():
67 ImageToImageGenericFilter<Self>("Register")
70 InitializeImageType<2>();
71 InitializeImageType<3>();
74 //==========================================================================================================//
75 //============================================================================================================//
76 //--------------------------------------------------------------------
77 template<class args_info_type>
78 template<unsigned int Dim>
79 void AffineRegistrationGenericFilter<args_info_type>::InitializeImageType()
81 ADD_DEFAULT_IMAGE_TYPES(Dim);
83 //--------------------------------------------------------------------
86 //==============================================================================
87 //Creating an observer class that allows us to change parameters at subsequent levels
88 //==============================================================================
89 template <typename TRegistration, class args_info_clitkAffineRegistration>
90 class RegistrationInterfaceCommand : public itk::Command
93 typedef RegistrationInterfaceCommand Self;
94 typedef itk::Command Superclass;
95 typedef itk::SmartPointer<Self> Pointer;
98 RegistrationInterfaceCommand() {};
102 typedef TRegistration RegistrationType;
103 typedef RegistrationType * RegistrationPointer;
106 typedef typename RegistrationType::FixedImageType InternalImageType;
107 typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
109 // Two arguments are passed to the Execute() method: the first
110 // is the pointer to the object which invoked the event and the
111 // second is the event that was invoked.
112 void Execute(itk::Object * object, const itk::EventObject & event) {
113 if ( !(itk::IterationEvent().CheckEvent( &event )) ) {
118 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
119 unsigned int numberOfLevels=registration->GetNumberOfLevels();
120 unsigned int currentLevel=registration->GetCurrentLevel()+1;
123 std::cout<<std::endl;
124 std::cout<<"========================================"<<std::endl;
125 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
126 std::cout<<"========================================"<<std::endl;
127 std::cout<<std::endl;
130 if (currentLevel>1) {
131 // Reinitialize the metric (!= number of samples)
132 typename GenericMetricType::Pointer genericMetric= GenericMetricType::New();
133 genericMetric->SetArgsInfo(m_ArgsInfo);
134 genericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
135 if (m_ArgsInfo.referenceMask_given) genericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
136 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
137 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
138 registration->SetMetric(metric);
142 void Execute(const itk::Object * , const itk::EventObject & ) {
146 void SetArgsInfo(args_info_clitkAffineRegistration a) {
149 args_info_clitkAffineRegistration m_ArgsInfo;
152 //==============================================================================
153 // Update with the number of dimensions and pixeltype
154 //==============================================================================
155 template<class args_info_clitkAffineRegistration>
156 template<class InputImageType>
157 void AffineRegistrationGenericFilter<args_info_clitkAffineRegistration>::UpdateWithInputImageType()
159 //=============================================================================
161 //=============================================================================
163 typedef typename InputImageType::PixelType PixelType;
164 //typedef typename InputImageType::ImageDimension Dimension;
167 bool threadsGiven=m_ArgsInfo.threads_given;
168 int threads=m_ArgsInfo.threads_arg;
170 //Coordinate Representation
171 typedef double TCoordRep;
173 //The pixeltype of the fixed image will be used for output
174 typedef itk::Image< PixelType, InputImageType::ImageDimension > FixedImageType;
176 //Whatever the pixel type, internally we work with an image represented in float
177 typedef typename InputImageType::PixelType InternalPixelType;
178 typedef itk::Image< InternalPixelType, InputImageType::ImageDimension > InternalImageType;
180 //Read in the reference/fixed image
181 typedef itk::ImageFileReader< InternalImageType > ReaderType;
182 typename ReaderType::Pointer fixedImageReader = ReaderType::New();
183 fixedImageReader->SetFileName( m_ArgsInfo.reference_arg);
186 //Read in the object/moving image
187 typename ReaderType::Pointer movingImageReader = ReaderType::New();
188 movingImageReader->SetFileName( m_ArgsInfo.target_arg );
189 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
190 fixedImageReader->Update();
191 movingImageReader->Update();
193 if (m_Verbose) std::cout << "Reading images... " << std::endl;
195 //we connect pointers to these internal images
196 typename InternalImageType::Pointer fixedImage= fixedImageReader->GetOutput();
197 typename InternalImageType::Pointer movingImage= movingImageReader->GetOutput();
199 //We keep the images used for input for possible output
200 typename InternalImageType::Pointer inputFixedImage= fixedImageReader->GetOutput();
201 typename InternalImageType::Pointer inputMovingImage= movingImageReader->GetOutput();
204 //============================================================================
206 //============================================================================
208 //If given, the intensities of both images are first normalized to a zero mean and SD of 1
209 // (usefull for MI, not necessary for Mattes' MI but performed anyway for the ouput)
210 if ( m_ArgsInfo.normalize_flag ) {
211 typedef itk::NormalizeImageFilter< InternalImageType,InternalImageType > NormalizeFilterType;
213 typename NormalizeFilterType::Pointer fixedNormalizeFilter = NormalizeFilterType::New();
214 typename NormalizeFilterType::Pointer movingNormalizeFilter = NormalizeFilterType::New();
216 fixedNormalizeFilter->SetInput( fixedImage );
217 movingNormalizeFilter->SetInput( movingImage );
219 fixedNormalizeFilter->Update();
220 movingNormalizeFilter->Update();
222 //We keep the images used for input for possible output
223 inputFixedImage= fixedNormalizeFilter->GetOutput();
224 inputMovingImage= movingNormalizeFilter->GetOutput();
226 //the pointers are reconnected for further output
227 fixedImage=fixedNormalizeFilter->GetOutput();
228 movingImage=movingNormalizeFilter->GetOutput();
230 if (m_Verbose) std::cout << "Normalizing image intensities to zero mean and SD of 1..." << std::endl;
234 //If given, the images are blurred before processing
235 if ( m_ArgsInfo.blur_arg!= 0) {
236 typedef itk::DiscreteGaussianImageFilter<InternalImageType,InternalImageType> GaussianFilterType;
237 typename GaussianFilterType::Pointer fixedSmoother = GaussianFilterType::New();
238 typename GaussianFilterType::Pointer movingSmoother = GaussianFilterType::New();
239 fixedSmoother->SetVariance( m_ArgsInfo.blur_arg );
240 movingSmoother->SetVariance(m_ArgsInfo.blur_arg );
242 fixedSmoother->SetInput( fixedImage );
243 movingSmoother->SetInput( movingImage );
245 fixedSmoother->Update();
246 movingSmoother->Update();
248 fixedImage=fixedSmoother->GetOutput();
249 movingImage=movingSmoother->GetOutput();
251 if (m_Verbose) std::cout << "Blurring images with a Gaussian with standard deviation of " << m_ArgsInfo.blur_arg <<"..." << std::endl;
255 //============================================================================
256 // Setting up the moving image in a reference system
257 //============================================================================
258 const itk::Vector<double, InputImageType::ImageDimension> movingResolution = movingImage->GetSpacing();
259 typename InternalImageType::RegionType movingRegion = movingImage->GetLargestPossibleRegion();
260 typename InternalImageType::RegionType::SizeType movingSize = movingRegion.GetSize();
262 // Print the parameters of the moving image
264 std::cout << "Object or Moving image:"<<std::endl;
265 std::cout << "Size: " << movingSize[0] << ", " << movingSize[1];
266 if (InputImageType::ImageDimension==3) std::cout<<", " << movingSize[2];
267 std::cout << std::endl;
269 std::cout<< "Resolution: "<< movingResolution[0] << ", " << movingResolution[1];
270 if (InputImageType::ImageDimension==3) std::cout<< ", " << movingResolution[2];
271 std::cout << std::endl;
275 //============================================================================
276 // Setting up the fixed image in a reference system
277 //============================================================================
278 const itk::Vector<double, InputImageType::ImageDimension> fixedResolution = fixedImage->GetSpacing();
279 typename InternalImageType::RegionType fixedRegion = fixedImage->GetLargestPossibleRegion();
280 typename InternalImageType::RegionType::SizeType fixedSize = fixedRegion.GetSize();
282 // Print the parameters of the moving image and the transform
284 std::cout << "Target or Moving image:"<<std::endl;
285 std::cout << "Size: " << fixedSize[0] << ", " << fixedSize[1];
286 if (InputImageType::ImageDimension==3) std::cout<<", " << fixedSize[2];
287 std::cout << std::endl;
289 std::cout<< "Resolution: "<< fixedResolution[0] << ", " << fixedResolution[1];
290 if (InputImageType::ImageDimension==3) std::cout<< ", " << fixedResolution[2];
291 std::cout << std::endl;
296 //===========================================================================
297 // If given, we connect a mask to reference or target
298 //============================================================================
299 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
300 typename MaskType::Pointer fixedMask=NULL;
301 if (m_ArgsInfo.referenceMask_given) {
302 fixedMask= MaskType::New();
303 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
304 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
305 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
306 maskReader->SetFileName(m_ArgsInfo.referenceMask_arg);
308 maskReader->Update();
309 } catch ( itk::ExceptionObject & err ) {
310 std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
311 std::cerr << err << std::endl;
314 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
316 // Set the image to the spatialObject
317 fixedMask->SetImage( maskReader->GetOutput() );
320 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
321 typename MaskType::Pointer movingMask=NULL;
322 if (m_ArgsInfo.targetMask_given) {
323 movingMask= MaskType::New();
324 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
325 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
326 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
327 maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
329 maskReader->Update();
330 } catch ( itk::ExceptionObject & err ) {
331 std::cerr << "ExceptionObject caught !" << std::endl;
332 std::cerr << err << std::endl;
334 if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
336 movingMask->SetImage( maskReader->GetOutput() );
340 //============================================================================
341 // The image pyramids
342 //============================================================================
343 typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType > FixedImagePyramidType;
344 typedef itk::RecursiveMultiResolutionPyramidImageFilter<InternalImageType,InternalImageType > MovingImagePyramidType;
345 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
346 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
347 fixedImagePyramid->SetUseShrinkImageFilter(false);
348 fixedImagePyramid->SetInput(fixedImage);
349 fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
350 movingImagePyramid->SetUseShrinkImageFilter(false);
351 movingImagePyramid->SetInput(movingImage);
352 movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
353 if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
354 fixedImagePyramid->Update();
355 movingImagePyramid->Update();
359 //============================================================================
360 // We retrieve the type of metric from the command line
361 //============================================================================
362 typedef clitk::GenericMetric<args_info_clitkAffineRegistration, InternalImageType, InternalImageType> GenericMetricType;
363 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
364 genericMetric->SetArgsInfo(m_ArgsInfo);
365 genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
366 if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
367 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
368 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
369 if (movingMask) metric->SetMovingImageMask(movingMask);
371 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
372 if (threadsGiven) metric->SetNumberOfThreads( threads );
374 if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
377 //============================================================================
378 // Initialize using image moments.
379 //============================================================================
380 if (m_ArgsInfo.moment_flag) {
381 typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
382 typename CalculatorType::Pointer fixedCalculator= CalculatorType::New();
384 typename InternalImageType::Pointer fixedThresh;
385 if (m_ArgsInfo.intThreshold_given) {
386 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
387 typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
388 thresholder->SetInput(fixedImage);
389 thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
390 thresholder->Update();
391 fixedThresh=thresholder->GetOutput();
392 } else fixedThresh=fixedImage;
394 fixedCalculator->SetImage(fixedThresh);
395 fixedCalculator->Compute();
396 Vector<double, InputImageType::ImageDimension> fixedCenter=fixedCalculator->GetCenterOfGravity();
397 if (m_Verbose)std::cout<<"The fixed center of gravity is "<<fixedCenter<<"..."<<std::endl;
399 typedef itk::ImageMomentsCalculator< InternalImageType > CalculatorType;
400 typename CalculatorType::Pointer movingCalculator= CalculatorType::New();
402 typename InternalImageType::Pointer movingThresh;
403 if (m_ArgsInfo.intThreshold_given) {
404 typedef itk::ThresholdImageFilter<InternalImageType> ThresholdImageFilterType;
405 typename ThresholdImageFilterType::Pointer thresholder = ThresholdImageFilterType::New();
406 thresholder->SetInput(movingImage);
407 thresholder->SetLower(m_ArgsInfo.intThreshold_arg);
408 thresholder->Update();
409 movingThresh=thresholder->GetOutput();
410 } else movingThresh=movingImage;
412 movingCalculator->SetImage(movingThresh);
413 movingCalculator->Compute();
414 Vector<double, InputImageType::ImageDimension> movingCenter=movingCalculator->GetCenterOfGravity();
415 if (m_Verbose)std::cout<<"The moving center of gravity is "<<movingCenter<<"..."<<std::endl;
417 Vector<double, InputImageType::ImageDimension> shift= movingCenter-fixedCenter;
418 if (m_Verbose)std::cout<<"The initial shift applied is "<<shift<<"..."<<std::endl;
420 m_ArgsInfo.transX_arg= shift [0];
421 m_ArgsInfo.transY_arg= shift [1];
422 if (InputImageType::ImageDimension==3) m_ArgsInfo.transZ_arg=shift [2];
425 //============================================================================
427 //============================================================================
428 typedef clitk::GenericAffineTransform<args_info_clitkAffineRegistration, TCoordRep, InputImageType::ImageDimension > GenericAffineTransformType;
429 typename GenericAffineTransformType::Pointer genericAffineTransform = GenericAffineTransformType::New();
430 genericAffineTransform->SetArgsInfo(m_ArgsInfo);
431 typedef itk::Transform< double, InputImageType::ImageDimension, InputImageType::ImageDimension > TransformType;
432 typename TransformType::Pointer transform = genericAffineTransform->GetTransform();
435 //=======================================================
437 //=======================================================
438 typedef clitk::GenericInterpolator<args_info_clitkAffineRegistration, InternalImageType,TCoordRep > GenericInterpolatorType;
439 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
440 genericInterpolator->SetArgsInfo(m_ArgsInfo);
441 typedef itk::InterpolateImageFunction< InternalImageType, TCoordRep > InterpolatorType;
442 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
445 //============================================================================
447 //============================================================================
448 typedef clitk::GenericOptimizer<args_info_clitkAffineRegistration> GenericOptimizerType;
449 unsigned int nParam = transform->GetNumberOfParameters();
450 typename GenericOptimizerType::Pointer genericOptimizer=GenericOptimizerType::New();
451 genericOptimizer->SetArgsInfo(m_ArgsInfo);
452 genericOptimizer->SetOutputIteration(m_Verbose);
453 genericOptimizer->SetOutputPosition(m_Verbose);
454 genericOptimizer->SetOutputValue(m_Verbose);
455 genericOptimizer->SetOutputGradient(m_ArgsInfo.gradient_flag);
456 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
457 genericOptimizer->SetNumberOfParameters(nParam);
458 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
459 OptimizerType::Pointer optimizer=genericOptimizer->GetOptimizerPointer();
462 itk::Optimizer::ScalesType scales( nParam );
463 for (unsigned int i=nParam-InputImageType::ImageDimension; i<nParam; i++) //Translations
464 scales[i] = m_ArgsInfo.tWeight_arg;
465 for (unsigned int i=0; i<nParam-InputImageType::ImageDimension; i++) //Rest
466 scales[i] = m_ArgsInfo.rWeight_arg*180./M_PI;
467 optimizer->SetScales(scales);
468 //============================================================================
469 // Multiresolution registration
470 //============================================================================
471 typedef itk::MultiResolutionImageRegistrationMethod< InternalImageType,InternalImageType > RegistrationType;
472 typename RegistrationType::Pointer registration = RegistrationType::New();
473 registration->SetFixedImage( fixedImage );
474 registration->SetFixedImageRegion(fixedImage->GetLargestPossibleRegion());
475 registration->SetMovingImage( movingImage );
476 registration->SetFixedImagePyramid( fixedImagePyramid );
477 registration->SetMovingImagePyramid( movingImagePyramid );
478 registration->SetTransform( transform );
479 registration->SetInitialTransformParameters( transform->GetParameters() );
480 registration->SetInterpolator( interpolator );
481 registration->SetMetric(metric);
482 registration->SetOptimizer(optimizer);
483 registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
484 if (m_Verbose) std::cout << "Setting "<< m_ArgsInfo.levels_arg <<" resolution levels..." << std::endl;
485 if (m_Verbose) std::cout << "Initial Transform: "<< registration->GetInitialTransformParameters()<<std::endl;
487 //============================================================================
488 // Connecting the commander to the registration to monitor it
489 //============================================================================
492 // Output iteration info
493 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
494 observer->SetOptimizer(genericOptimizer);
495 optimizer->AddObserver( itk::IterationEvent(), observer );
499 typedef RegistrationInterfaceCommand<RegistrationType, args_info_clitkAffineRegistration> CommandType;
500 typename CommandType::Pointer command = CommandType::New();
501 command->SetArgsInfo(m_ArgsInfo);
502 registration->AddObserver( itk::IterationEvent(), command );
507 //============================================================================
508 // Finally we can start the registration with the given amount of multiresolution levels
509 //============================================================================
510 if (m_Verbose) std::cout << "Starting the registration now..." << std::endl;
513 registration->StartRegistration();
514 } catch ( itk::ExceptionObject & err ) {
515 std::cerr << "ExceptionObject caught !" << std::endl;
516 std::cerr << err << std::endl;
520 //============================================================================
521 // Processing the result of the registration
522 //============================================================================
523 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
524 std::cout<< "Result : " <<std::setprecision(12)<<std::endl;
526 for (unsigned int i=nParam-InputImageType::ImageDimension; i<nParam; i++) //Translations
527 std::cout << " Translation " << i << " = " << finalParameters[i];
528 for (unsigned int i=0; i<nParam-InputImageType::ImageDimension; i++) //Rest
529 std::cout << " Other parameter " << i << " = " << finalParameters[i];
532 itk::Matrix<double,InputImageType::ImageDimension+1,InputImageType::ImageDimension+1> matrix;
533 if (m_ArgsInfo.transform_arg == 3) {
534 for (unsigned int i=0; i<InputImageType::ImageDimension; i++) {
535 matrix[i][3] = finalParameters[nParam-InputImageType::ImageDimension+i];
536 for (unsigned int j=0; j<InputImageType::ImageDimension; j++) {
537 matrix[i][j] = finalParameters[i*3+j];
542 matrix = clitk::GetBackwardAffineMatrix<InputImageType::ImageDimension>(finalParameters);
545 std::cout << " Affine transform matrix =" << std::endl;
546 std::cout << matrix <<std::setprecision(6)<< std::endl;
548 // Write matrix to a file
549 if (m_ArgsInfo.matrix_given) {
551 mFile.open(m_ArgsInfo.matrix_arg);
552 mFile<<std::setprecision(12)<<matrix<< std::setprecision(6)<<std::endl;
556 //============================================================================
557 // Prepare the resampling filter in order to transform the moving image.
558 //============================================================================
559 if (m_ArgsInfo.output_given || m_ArgsInfo.checker_after_given || m_ArgsInfo.after_given ) {
560 transform->SetParameters( finalParameters );
561 typedef itk::ResampleImageFilter< InternalImageType,InternalImageType > ResampleFilterType;
562 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
564 resampler->SetTransform( transform );
565 resampler->SetInput( movingImageReader->GetOutput() );
566 resampler->SetSize( fixedImageReader ->GetOutput()->GetLargestPossibleRegion().GetSize() );
567 resampler->SetOutputOrigin( fixedImageReader ->GetOutput()->GetOrigin() );
568 resampler->SetOutputSpacing( fixedImageReader ->GetOutput()->GetSpacing() );
569 resampler->SetDefaultPixelValue( 0 );
572 if (m_ArgsInfo.output_given) {
573 //We write an output in the same pixeltype then the input
574 typedef itk::ImageFileWriter< FixedImageType > WriterType;
575 typename WriterType::Pointer outputWriter = WriterType::New();
576 outputWriter->SetFileName(m_ArgsInfo.output_arg );
577 outputWriter->SetInput( resampler->GetOutput() );
578 outputWriter->Update();
582 //============================================================================
584 //============================================================================
585 if (m_ArgsInfo.checker_after_given) {
586 //To display correctly the checkerboard image, the intensities must lie in the same range (normalized)
587 //We write the image in the internal image type
588 typedef itk::ResampleImageFilter< InternalImageType,InternalImageType > ResampleFilterType;
589 typename ResampleFilterType::Pointer internalResampler = ResampleFilterType::New();
590 internalResampler->SetTransform( transform );
591 internalResampler->SetInput( inputMovingImage );
592 internalResampler->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
593 internalResampler->SetOutputOrigin( fixedImage->GetOrigin() );
594 internalResampler->SetOutputSpacing( fixedImage->GetSpacing() );
595 internalResampler->SetDefaultPixelValue( 0 );
597 //We pass the normalized images to the checker filter
598 typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType;
599 typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New();
601 checkerFilter->SetInput1(inputFixedImage);
602 checkerFilter->SetInput2(internalResampler->GetOutput());
603 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
604 typename InternalWriterType::Pointer checkerWriter = InternalWriterType::New();
605 checkerWriter->SetFileName(m_ArgsInfo.checker_after_arg);
606 checkerWriter->SetInput( checkerFilter->GetOutput() );
607 checkerWriter->Update();
611 //============================================================================
613 //============================================================================
614 if (m_ArgsInfo.checker_before_given) {
615 //To display correctly the checkerboard image, the intensities must lie in the same range (normalized)
616 //We write the image in the internal image type
617 //We pass the normalized images to the checker filter
618 typedef itk::CheckerBoardImageFilter< InternalImageType > CheckerBoardFilterType;
619 typename CheckerBoardFilterType::Pointer checkerFilter= CheckerBoardFilterType::New();
621 checkerFilter->SetInput1(inputFixedImage);
622 checkerFilter->SetInput2(inputMovingImage);
623 typedef itk::ImageFileWriter< InternalImageType > InternalWriterType;
624 typename InternalWriterType::Pointer checkerWriter = InternalWriterType::New();
625 checkerWriter->SetFileName(m_ArgsInfo.checker_before_arg);
626 checkerWriter->SetInput( checkerFilter->GetOutput() );
627 checkerWriter->Update();
631 //============================================================================
633 //============================================================================
634 if (m_ArgsInfo.after_given) {
635 typedef itk::SubtractImageFilter< InternalImageType, FixedImageType,FixedImageType > DifferenceImageFilterType;
636 typename DifferenceImageFilterType::Pointer differenceAfterFilter= DifferenceImageFilterType::New();
638 differenceAfterFilter->SetInput1(fixedImageReader ->GetOutput());
639 differenceAfterFilter->SetInput2(resampler->GetOutput());
641 // Prepare a writer to write the difference image
642 typedef itk::ImageFileWriter< FixedImageType > WriterType;
643 typename WriterType::Pointer differenceAfterWriter = WriterType::New();
644 differenceAfterWriter->SetFileName(m_ArgsInfo.after_arg );
645 differenceAfterWriter->SetInput( differenceAfterFilter->GetOutput() );
646 differenceAfterWriter->Update();
650 //============================================================================
651 // Difference Before?
652 //============================================================================
653 if (m_ArgsInfo.before_given) {
654 typedef itk::CastImageFilter< InternalImageType,FixedImageType > CastFilterType;
655 typename CastFilterType::Pointer caster = CastFilterType::New();
656 caster->SetInput( movingImageReader->GetOutput() );
658 typedef itk::SubtractImageFilter< InternalImageType, FixedImageType, FixedImageType > DifferenceImageFilterType;
659 typename DifferenceImageFilterType::Pointer differenceBeforeFilter= DifferenceImageFilterType::New();
662 differenceBeforeFilter->SetInput1(fixedImageReader ->GetOutput());
663 differenceBeforeFilter->SetInput2(caster->GetOutput());
665 // Prepare a writer to write the difference image
666 typedef itk::ImageFileWriter< FixedImageType > WriterType;
667 typename WriterType::Pointer differenceBeforeWriter = WriterType::New();
668 differenceBeforeWriter->SetFileName(m_ArgsInfo.before_arg);
669 differenceBeforeWriter->SetInput( differenceBeforeFilter->GetOutput() );
670 differenceBeforeWriter->Update();
675 #endif //#define CLITKAFFINEREGISTRATIONGENERICFILTER_TXX