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 ===========================================================================**/
18 #ifndef __clitkDemonsDeformableRegistrationGenericFilter_txx
19 #define __clitkDemonsDeformableRegistrationGenericFilter_txx
20 #include "clitkDemonsDeformableRegistrationGenericFilter.h"
26 //==============================================================================
27 //Creating an observer class that allows us to change parameters at subsequent levels
28 //==============================================================================
29 template <typename MultiRegistrationFilterType>
30 class CommandResolutionLevelUpdate : public itk::Command
33 typedef CommandResolutionLevelUpdate Self;
34 typedef itk::Command Superclass;
35 typedef itk::SmartPointer<Self> Pointer;
39 itkStaticConstMacro(ImageDimension, unsigned int, MultiRegistrationFilterType::ImageDimension);
40 typedef typename MultiRegistrationFilterType::FixedImageType FixedImageType;
41 typedef typename MultiRegistrationFilterType::MovingImageType MovingImageType;
42 typedef typename MultiRegistrationFilterType::DeformationFieldType DeformationFieldType;
43 typedef typename MultiRegistrationFilterType::RegistrationType RegistrationType;
44 typedef typename RegistrationType::FixedImageType InternalImageType;
45 typedef itk::DiffeomorphicDemonsRegistrationFilter<InternalImageType, InternalImageType, DeformationFieldType> DiffFilterType;
46 typedef itk::FastSymmetricForcesDemonsRegistrationFilter<InternalImageType, InternalImageType, DeformationFieldType> SymFilterType;
49 CommandResolutionLevelUpdate()
60 void SetMaxRMSError(double* m){m_MaxRMSError=m;}
61 void SetSD(double* m){m_SD=m;}
62 void SetMaxStep(double m){m_MaxStep=m;}
63 void SetScaleSD(bool m){m_ScaleSD=m;}
64 void SetScaleStep(bool m){m_ScaleStep=m;}
65 void SetRegistrationType(unsigned int i){m_RegistrationType=i;}
68 void Execute(const itk::Object * caller, const itk::EventObject & event )
70 std::cout << "Warning: The const Execute method shouldn't be called" << std::endl;
72 void Execute(itk::Object *caller, const itk::EventObject & event)
76 if( !(itk::IterationEvent().CheckEvent( &event )) )
82 MultiRegistrationFilterType * filter = dynamic_cast< MultiRegistrationFilterType * >( caller );
85 m_CurrentLevel=filter->GetCurrentLevel();
86 unsigned int numberOfLevels=filter->GetNumberOfLevels();
87 unsigned int expandFactor=1<< (numberOfLevels-m_CurrentLevel-1);
89 double maxStep=m_MaxStep;
94 for (unsigned int i=0 ; i<ImageDimension ; i++) sd[i]*=(double)expandFactor;
95 filter->GetRegistrationFilter()->SetStandardDeviations(sd);
96 filter->GetRegistrationFilter()->SetUpdateFieldStandardDeviations(sd);
99 // Scale max step: 2 & 3 only!
102 maxStep= m_MaxStep*expandFactor;
103 if (m_RegistrationType==2)
105 SymFilterType* symFilter=dynamic_cast< SymFilterType* >(filter->GetRegistrationFilter());
106 symFilter->SetMaximumUpdateStepLength(maxStep);
108 else if (m_RegistrationType==3)
110 DiffFilterType* diffFilter=dynamic_cast< DiffFilterType* >(filter->GetRegistrationFilter());
111 diffFilter->SetMaximumUpdateStepLength(maxStep);
116 filter->GetRegistrationFilter()->SetMaximumRMSError(m_MaxRMSError[m_CurrentLevel]);
119 std::cout << "--------------------------------------------------" << std::endl;
120 std::cout << "Starting resolution level "<<m_CurrentLevel+1<<" of "<<filter->GetNumberOfLevels()<<"..."<<std::endl;
121 std::cout << "Setting the standard deviations to [ "<<sd[0];
122 for (unsigned int i=1; i<3; i++) std::cout<<", "<< sd[i];
123 std::cout <<" ]..."<<std::endl;
124 if ( (m_RegistrationType==2) || (m_RegistrationType==3) ) std::cout <<"Setting the maximum step size to "<< maxStep<<"..."<<std::endl;
125 std::cout <<"Setting the maximum RMS field change to "<< m_MaxRMSError[m_CurrentLevel]<<"..."<<std::endl;
126 std::cout << "--------------------------------------------------" << std::endl;
129 double *m_MaxRMSError,*m_SD;
131 unsigned int m_CurrentLevel, m_RegistrationType;
132 bool m_ScaleSD, m_ScaleStep;
137 //==============================================================================
138 // Creating an observer class that allows output at each iteration
139 //==============================================================================
140 template <typename RegistrationFilterType, typename FixedImageType, typename MovingImageType>
141 class CommandIterationUpdate : public itk::Command
144 typedef CommandIterationUpdate Self;
145 typedef itk::Command Superclass;
146 typedef itk::SmartPointer<Self> Pointer;
147 typedef itk::SmartPointer<const Self> ConstPointer;
151 //find the multiresolution filter
152 // typedef typename RegistrationFilterType::FixedImageType InternalImageType;
153 // typedef typename RegistrationFilterType::MovingImageType MovingImageType;
154 typedef typename RegistrationFilterType::DisplacementFieldType DisplacementFieldType;
155 typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType, MovingImageType, DisplacementFieldType> MultiResolutionRegistrationType;
156 typedef CommandResolutionLevelUpdate<MultiResolutionRegistrationType> LevelObserver;
159 CommandIterationUpdate(){};
164 void SetStop( int* s){m_Stop=s;}
165 void SetLevelObserver(LevelObserver* o ){m_LevelObserver=o;}
169 void Execute(const itk::Object *, const itk::EventObject & )
171 std::cout << "Warning: The const Execute method shouldn't be called" << std::endl;
174 void Execute(itk::Object *caller, const itk::EventObject & event)
176 if( !(itk::IterationEvent().CheckEvent( &event )) )
182 RegistrationFilterType * filter = dynamic_cast< RegistrationFilterType * >( caller );
187 m_Iteration=filter->GetElapsedIterations();
188 m_Metric=filter->GetMetric();
191 std::cout << m_Iteration<<"\t Field Tolerance= "<<filter->GetMaximumRMSError();
192 std::cout <<"\t DVF Change= " << filter->GetRMSChange()<<"\t RMS= "<<m_Metric;
194 // Using stop criterion?
195 if (m_Stop[m_LevelObserver->m_CurrentLevel]>=0)
205 else if(m_Metric<m_Minimum)
211 //Not less then minimum
215 if (m_StopCounter>=m_Stop[m_LevelObserver->m_CurrentLevel])
216 filter->StopRegistration();
220 std::cout <<"\t Stop= "<<m_StopCounter<<" / "<<m_Stop[m_LevelObserver->m_CurrentLevel]<<std::endl;
223 // Not using stop criterion
224 else std::cout <<std::endl;
228 double m_Minimum, m_Metric;
229 int m_StopCounter, m_Iteration;
231 LevelObserver* m_LevelObserver;
235 //==============================================================================
236 // Update with the number of dimensions
237 //==============================================================================
238 template<unsigned int Dimension>
239 void DemonsDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
241 if(PixelType == "short"){
242 if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl;
243 UpdateWithDimAndPixelType<Dimension, signed short>();
245 // else if(PixelType == "unsigned_short"){
246 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl;
247 // UpdateWithDimAndPixelType<Dimension, unsigned short>();
250 // else if (PixelType == "unsigned_char"){
251 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl;
252 // UpdateWithDimAndPixelType<Dimension, unsigned char>();
255 // else if (PixelType == "char"){
256 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl;
257 // UpdateWithDimAndPixelType<Dimension, signed char>();
260 if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and float..." << std::endl;
261 UpdateWithDimAndPixelType<Dimension, float>();
267 //==============================================================================
268 // Update with the number of dimensions and pixeltype
269 //==============================================================================
270 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
272 template<unsigned int ImageDimension, class PixelType>
273 void DemonsDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType()
275 //=======================================================
277 //=======================================================
278 bool threadsGiven=m_ArgsInfo.threads_given;
279 int threads=m_ArgsInfo.threads_arg;
281 typedef itk::Image< PixelType, ImageDimension > FixedImageType;
282 typedef itk::Image< PixelType, ImageDimension > MovingImageType;
283 typedef itk::Image< float, ImageDimension > InternalImageType;
284 typedef double TCoordRep;
287 //=======================================================
289 //=======================================================
290 typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
291 typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
293 typename FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
294 typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
296 fixedImageReader->SetFileName( m_ArgsInfo.reference_arg );
297 movingImageReader->SetFileName( m_ArgsInfo.target_arg );
298 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
299 fixedImageReader->Update();
300 movingImageReader->Update();
302 typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
303 typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
304 typename FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
307 //=======================================================
309 //=======================================================
310 typedef itk::Vector< float, ImageDimension > VectorPixelType;
311 typedef itk::Image< VectorPixelType, ImageDimension > DeformationFieldType;
314 //=======================================================
316 //=======================================================
317 typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, InternalImageType> FixedImagePyramidType;
318 typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, InternalImageType> MovingImagePyramidType;
319 // typedef itk::MultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
320 // typedef itk::MultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
321 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
322 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
325 // =======================================================
326 // Define the registation filters
327 // =======================================================
328 typedef itk::PDEDeformableRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType> RegistrationFilterType;
329 typename RegistrationFilterType::Pointer pdeFilter;
330 typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType,MovingImageType,DeformationFieldType > MultiResolutionRegistrationFilterType;
331 typename MultiResolutionRegistrationFilterType::Pointer multiResolutionFilter = MultiResolutionRegistrationFilterType::New();
332 typedef itk::ESMDemonsRegistrationFunction<InternalImageType, InternalImageType, DeformationFieldType> RegistrationFunctionType;
335 // =======================================================
336 // The multiresolution scheme
337 // =======================================================
339 #if ITK_VERSION_MAJOR <= 4
340 multiResolutionFilter->SetNumberOfThreads(threads);
342 multiResolutionFilter->SetNumberOfWorkUnits(threads);
345 unsigned int nLevels=m_ArgsInfo.levels_arg;
346 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<nLevels<<"..."<<std::endl;
347 multiResolutionFilter->SetFixedImage( fixedImage );
348 multiResolutionFilter->SetMovingImage( movingImage );
349 multiResolutionFilter->SetNumberOfLevels( nLevels );
350 multiResolutionFilter->SetFixedImagePyramid( fixedImagePyramid );
351 multiResolutionFilter->SetMovingImagePyramid( movingImagePyramid );
353 #if ITK_VERSION_MAJOR <= 4
354 multiResolutionFilter->SetNumberOfThreads( threads );
356 multiResolutionFilter->SetNumberOfWorkUnits( threads );
360 //------------------------------------
361 //Set the number of iterations
362 //------------------------------------
363 std::vector<unsigned int> nIterations(nLevels);
364 for (unsigned int i=0 ; i<nLevels; i++)
366 if (m_ArgsInfo.maxIter_given==nLevels)
368 nIterations[i] =m_ArgsInfo.maxIter_arg[i];
371 nIterations[i]=m_ArgsInfo.maxIter_arg[0];
373 multiResolutionFilter->SetNumberOfIterations( &(nIterations[0]) );
375 std::cout<<"Setting the number of iterations to: "<<nIterations[0];
376 for (unsigned int i=1; i<nLevels; i++)
377 std::cout<<", "<<nIterations [i];
378 std::cout<<std::endl;
381 //------------------------------------
382 //Set the max RMS error for the field update
383 //------------------------------------
384 std::vector<double> maxRMSError(nLevels);
385 for (unsigned int i=0 ; i<nLevels; i++)
387 if (m_ArgsInfo.maxRMSError_given==nLevels)
388 maxRMSError[i] =m_ArgsInfo.maxRMSError_arg[i];
390 maxRMSError[i]=m_ArgsInfo.maxRMSError_arg[0];
393 std::cout<<"Setting the max root mean squared error to: "<<maxRMSError[0];
394 for (unsigned int i=1; i<nLevels; i++)
395 std::cout<<", "<<maxRMSError[i];
396 std::cout<<std::endl;
399 //------------------------------------
400 //Get the stop criterion
401 //------------------------------------
402 std::vector<int> stop(nLevels);
403 for (unsigned int i=0; i<nLevels; i++)
404 if (m_ArgsInfo.stop_given==nLevels)
405 stop[i]=m_ArgsInfo.stop_arg[i];
407 stop[i]=m_ArgsInfo.stop_arg[0];
409 std::cout<<"Setting the stop criterion to : "<<stop[0];
410 for (unsigned int i=1; i<nLevels; i++)
411 std::cout<<", "<<stop[i];
412 std::cout<<std::endl;
415 //------------------------------------
417 //------------------------------------
418 typename RegistrationFunctionType::GradientType grad=RegistrationFunctionType::Symmetric;
419 switch (m_ArgsInfo.gradType_arg)
422 grad= RegistrationFunctionType::Symmetric;
425 grad= RegistrationFunctionType::Fixed;
428 grad= RegistrationFunctionType::WarpedMoving;
431 grad= RegistrationFunctionType::MappedMoving;
435 //------------------------------------
436 // Set smoothing standard deviations
437 //------------------------------------
438 double sd [ImageDimension];
439 for (unsigned int i=0; i<ImageDimension; i++)
440 if (m_ArgsInfo.sd_given==ImageDimension)
441 sd[i]=m_ArgsInfo.sd_arg[i];
443 sd[i]=m_ArgsInfo.sd_arg[0];
445 std::cout<<"Setting the final standard deviations for smoothing to: "<<sd[0];
446 for (unsigned int i=1; i<ImageDimension; i++)
447 std::cout<<", "<<sd[i];
448 std::cout<<std::endl;
451 // =======================================================
452 // The level observer: adjust settings at each level
453 // =======================================================
454 typename CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::Pointer levelObserver =
455 CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::New();
456 multiResolutionFilter->AddObserver( itk::IterationEvent(), levelObserver );
457 levelObserver->SetMaxRMSError(&(maxRMSError[0]));
458 levelObserver->SetMaxStep(m_ArgsInfo.maxStep_arg);
459 levelObserver->SetSD(sd);
460 levelObserver->SetScaleStep(m_ArgsInfo.scaleStep_flag);
461 levelObserver->SetScaleSD(m_ArgsInfo.scaleSD_flag);
462 levelObserver->SetRegistrationType(m_ArgsInfo.demons_arg);
465 // =======================================================
466 // The type of filter
467 // =======================================================
468 switch (m_ArgsInfo.demons_arg){
472 typedef itk::DemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
473 typename DemonsFilterType::Pointer m =DemonsFilterType::New();
475 //Set Parameters for this filter
476 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
477 m->SetUseMovingImageGradient( m_ArgsInfo.movGrad_flag);
478 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
479 observer->SetStop(&(stop[0]));
480 observer->SetLevelObserver(levelObserver);
481 m->AddObserver( itk::IterationEvent(), observer );
482 if (m_Verbose) std::cout<<"Using the demons registration filter..."<<std::endl;
489 typedef itk::SymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
490 typename DemonsFilterType::Pointer m =DemonsFilterType::New();
492 //Set Parameters for this filter
493 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
494 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
495 observer->SetStop(&(stop[0]));
496 observer->SetLevelObserver(levelObserver);
497 m->AddObserver( itk::IterationEvent(), observer );
498 if (m_Verbose) std::cout<<"Using the symmetric forces demons registration filter..."<<std::endl;
505 typedef itk::FastSymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
506 typename DemonsFilterType::Pointer m = DemonsFilterType::New();
508 //Set Parameters for this filter
509 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
510 m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
511 m->SetUseGradientType(grad);
512 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
513 observer->SetStop(&(stop[0]));
514 observer->SetLevelObserver(levelObserver);
515 m->AddObserver( itk::IterationEvent(), observer );
516 if (m_Verbose) std::cout<<"Using the fast symmetric forces demons registration filter..."<<std::endl;
523 typedef itk::DiffeomorphicDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
524 typename DemonsFilterType::Pointer m = DemonsFilterType::New();
526 //Set Parameters for this filter
527 m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
528 m->SetUseFirstOrderExp(m_ArgsInfo.firstOrder_flag);
529 m->SetUseGradientType(grad);
530 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
531 observer->SetStop(&(stop[0]));
532 observer->SetLevelObserver(levelObserver);
533 m->AddObserver( itk::IterationEvent(), observer );
534 if (m_Verbose) std::cout<<"Using the diffeomorphic demons registration filter..."<<std::endl;
543 pdeFilter->SetStandardDeviations( sd );
544 pdeFilter->SetUpdateFieldStandardDeviations( sd );
546 // pdeFilter->SetMaximumError(m_ArgsInfo.maxError_arg);
547 // pdeFilter->SetMaximumKernelWidth(m_ArgsInfo.maxError_arg);
548 pdeFilter->SetSmoothDisplacementField(!m_ArgsInfo.fluid_flag);
549 pdeFilter->SetSmoothUpdateField(m_ArgsInfo.fluid_flag);
550 pdeFilter->SetUseImageSpacing( m_ArgsInfo.spacing_flag );
552 //Pass to the multi resolution scheme
553 multiResolutionFilter->SetRegistrationFilter( pdeFilter );
556 // =======================================================
557 // The initial solution
558 // =======================================================
559 if (m_ArgsInfo.init_given)
561 typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
562 typename DeformationFieldReaderType::Pointer defReader=DeformationFieldReaderType::New();
563 defReader->SetFileName(m_ArgsInfo.init_arg);
565 multiResolutionFilter->SetArbitraryInitialDeformationField(defReader->GetOutput());
569 // =======================================================
571 // =======================================================
574 multiResolutionFilter->Update();
576 catch( itk::ExceptionObject & excp )
578 std::cerr <<"Error executing the demons filter: "<< excp << std::endl;
583 //=======================================================
585 //=======================================================
586 typename DeformationFieldType::Pointer deformationField=multiResolutionFilter->GetOutput();
589 //=======================================================
591 //=======================================================
592 typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
593 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
594 fieldWriter->SetInput( deformationField );
595 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
598 fieldWriter->Update();
600 catch( itk::ExceptionObject & excp )
602 std::cerr << "Exception thrown writing the DVF" << std::endl;
603 std::cerr << excp << std::endl;
608 //=======================================================
609 // Warp the moving image
610 //=======================================================
611 typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType;
612 typename WarpFilterType::Pointer warp = WarpFilterType::New();
614 warp->SetDisplacementField( deformationField );
615 warp->SetInput( movingImageReader->GetOutput() );
616 warp->SetOutputOrigin( fixedImage->GetOrigin() );
617 warp->SetOutputSpacing( fixedImage->GetSpacing() );
618 warp->SetOutputDirection( fixedImage->GetDirection() );
619 warp->SetEdgePaddingValue( 0.0 );
623 //=======================================================
624 // Write the warped image
625 //=======================================================
626 typedef itk::ImageFileWriter< FixedImageType > WriterType;
627 typename WriterType::Pointer writer = WriterType::New();
628 writer->SetFileName( m_ArgsInfo.output_arg );
629 writer->SetInput( warp->GetOutput() );
635 catch( itk::ExceptionObject & err )
637 std::cerr << "ExceptionObject caught !" << std::endl;
638 std::cerr << err << std::endl;
643 //=======================================================
644 // Calculate the difference after the deformable transform
645 //=======================================================
646 typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType;
647 if (m_ArgsInfo.after_given)
649 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
650 difference->SetValidInput( fixedImage );
651 difference->SetTestInput( warp->GetOutput() );
655 difference->Update();
657 catch( itk::ExceptionObject & err )
659 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
660 std::cerr << err << std::endl;
664 typename WriterType::Pointer differenceWriter=WriterType::New();
665 differenceWriter->SetInput(difference->GetOutput());
666 differenceWriter->SetFileName(m_ArgsInfo.after_arg);
667 differenceWriter->Update();
672 //=======================================================
673 // Calculate the difference before the deformable transform
674 //=======================================================
675 if( m_ArgsInfo.before_given )
677 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
678 difference->SetValidInput( fixedImage );
679 difference->SetTestInput( movingImage );
683 difference->Update();
685 catch( itk::ExceptionObject & err )
687 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
688 std::cerr << err << std::endl;
692 typename WriterType::Pointer differenceWriter=WriterType::New();
693 writer->SetFileName( m_ArgsInfo.before_arg );
694 writer->SetInput( difference->GetOutput() );
703 #endif // __clitkDemonsDeformableRegistrationGenericFilter_txx