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 ======================================================================-====*/
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::DeformationFieldType DeformationFieldType;
155 typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType, MovingImageType, DeformationFieldType> 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 // =======================================================
338 if (threadsGiven) multiResolutionFilter->SetNumberOfThreads(threads);
339 unsigned int nLevels=m_ArgsInfo.levels_arg;
340 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<nLevels<<"..."<<std::endl;
341 multiResolutionFilter->SetFixedImage( fixedImage );
342 multiResolutionFilter->SetMovingImage( movingImage );
343 multiResolutionFilter->SetNumberOfLevels( nLevels );
344 multiResolutionFilter->SetFixedImagePyramid( fixedImagePyramid );
345 multiResolutionFilter->SetMovingImagePyramid( movingImagePyramid );
346 if (threadsGiven) multiResolutionFilter->SetNumberOfThreads( threads );
348 //------------------------------------
349 //Set the number of iterations
350 //------------------------------------
351 std::vector<unsigned int> nIterations(nLevels);
352 for (unsigned int i=0 ; i<nLevels; i++)
354 if (m_ArgsInfo.maxIter_given==nLevels)
356 nIterations[i] =m_ArgsInfo.maxIter_arg[i];
359 nIterations[i]=m_ArgsInfo.maxIter_arg[0];
361 multiResolutionFilter->SetNumberOfIterations( &(nIterations[0]) );
363 std::cout<<"Setting the number of iterations to: "<<nIterations[0];
364 for (unsigned int i=1; i<nLevels; i++)
365 std::cout<<", "<<nIterations [i];
366 std::cout<<std::endl;
369 //------------------------------------
370 //Set the max RMS error for the field update
371 //------------------------------------
372 std::vector<double> maxRMSError(nLevels);
373 for (unsigned int i=0 ; i<nLevels; i++)
375 if (m_ArgsInfo.maxRMSError_given==nLevels)
376 maxRMSError[i] =m_ArgsInfo.maxRMSError_arg[i];
378 maxRMSError[i]=m_ArgsInfo.maxRMSError_arg[0];
381 std::cout<<"Setting the max root mean squared error to: "<<maxRMSError[0];
382 for (unsigned int i=1; i<nLevels; i++)
383 std::cout<<", "<<maxRMSError[i];
384 std::cout<<std::endl;
387 //------------------------------------
388 //Get the stop criterion
389 //------------------------------------
390 std::vector<int> stop(nLevels);
391 for (unsigned int i=0; i<nLevels; i++)
392 if (m_ArgsInfo.stop_given==nLevels)
393 stop[i]=m_ArgsInfo.stop_arg[i];
395 stop[i]=m_ArgsInfo.stop_arg[0];
397 std::cout<<"Setting the stop criterion to : "<<stop[0];
398 for (unsigned int i=1; i<nLevels; i++)
399 std::cout<<", "<<stop[i];
400 std::cout<<std::endl;
403 //------------------------------------
405 //------------------------------------
406 typename RegistrationFunctionType::GradientType grad=RegistrationFunctionType::Symmetric;
407 switch (m_ArgsInfo.gradType_arg)
410 grad= RegistrationFunctionType::Symmetric;
413 grad= RegistrationFunctionType::Fixed;
416 grad= RegistrationFunctionType::WarpedMoving;
419 grad= RegistrationFunctionType::MappedMoving;
423 //------------------------------------
424 // Set smoothing standard deviations
425 //------------------------------------
426 double sd [ImageDimension];
427 for (unsigned int i=0; i<ImageDimension; i++)
428 if (m_ArgsInfo.sd_given==ImageDimension)
429 sd[i]=m_ArgsInfo.sd_arg[i];
431 sd[i]=m_ArgsInfo.sd_arg[0];
433 std::cout<<"Setting the final standard deviations for smoothing to: "<<sd[0];
434 for (unsigned int i=1; i<ImageDimension; i++)
435 std::cout<<", "<<sd[i];
436 std::cout<<std::endl;
439 // =======================================================
440 // The level observer: adjust settings at each level
441 // =======================================================
442 typename CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::Pointer levelObserver =
443 CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::New();
444 multiResolutionFilter->AddObserver( itk::IterationEvent(), levelObserver );
445 levelObserver->SetMaxRMSError(&(maxRMSError[0]));
446 levelObserver->SetMaxStep(m_ArgsInfo.maxStep_arg);
447 levelObserver->SetSD(sd);
448 levelObserver->SetScaleStep(m_ArgsInfo.scaleStep_flag);
449 levelObserver->SetScaleSD(m_ArgsInfo.scaleSD_flag);
450 levelObserver->SetRegistrationType(m_ArgsInfo.demons_arg);
453 // =======================================================
454 // The type of filter
455 // =======================================================
456 switch (m_ArgsInfo.demons_arg){
460 typedef itk::DemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
461 typename DemonsFilterType::Pointer m =DemonsFilterType::New();
463 //Set Parameters for this filter
464 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
465 m->SetUseMovingImageGradient( m_ArgsInfo.movGrad_flag);
466 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
467 observer->SetStop(&(stop[0]));
468 observer->SetLevelObserver(levelObserver);
469 m->AddObserver( itk::IterationEvent(), observer );
470 if (m_Verbose) std::cout<<"Using the demons registration filter..."<<std::endl;
477 typedef itk::SymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
478 typename DemonsFilterType::Pointer m =DemonsFilterType::New();
480 //Set Parameters for this filter
481 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
482 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
483 observer->SetStop(&(stop[0]));
484 observer->SetLevelObserver(levelObserver);
485 m->AddObserver( itk::IterationEvent(), observer );
486 if (m_Verbose) std::cout<<"Using the symmetric forces demons registration filter..."<<std::endl;
493 typedef itk::FastSymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
494 typename DemonsFilterType::Pointer m = DemonsFilterType::New();
496 //Set Parameters for this filter
497 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
498 m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
499 m->SetUseGradientType(grad);
500 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
501 observer->SetStop(&(stop[0]));
502 observer->SetLevelObserver(levelObserver);
503 m->AddObserver( itk::IterationEvent(), observer );
504 if (m_Verbose) std::cout<<"Using the fast symmetric forces demons registration filter..."<<std::endl;
511 typedef itk::DiffeomorphicDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
512 typename DemonsFilterType::Pointer m = DemonsFilterType::New();
514 //Set Parameters for this filter
515 m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
516 m->SetUseFirstOrderExp(m_ArgsInfo.firstOrder_flag);
517 m->SetUseGradientType(grad);
518 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
519 observer->SetStop(&(stop[0]));
520 observer->SetLevelObserver(levelObserver);
521 m->AddObserver( itk::IterationEvent(), observer );
522 if (m_Verbose) std::cout<<"Using the diffeomorphic demons registration filter..."<<std::endl;
531 pdeFilter->SetStandardDeviations( sd );
532 pdeFilter->SetUpdateFieldStandardDeviations( sd );
534 // pdeFilter->SetMaximumError(m_ArgsInfo.maxError_arg);
535 // pdeFilter->SetMaximumKernelWidth(m_ArgsInfo.maxError_arg);
536 pdeFilter->SetSmoothDeformationField(!m_ArgsInfo.fluid_flag);
537 pdeFilter->SetSmoothUpdateField(m_ArgsInfo.fluid_flag);
538 pdeFilter->SetUseImageSpacing( m_ArgsInfo.spacing_flag );
540 //Pass to the multi resolution scheme
541 multiResolutionFilter->SetRegistrationFilter( pdeFilter );
544 // =======================================================
545 // The initial solution
546 // =======================================================
547 if (m_ArgsInfo.init_given)
549 typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
550 typename DeformationFieldReaderType::Pointer defReader=DeformationFieldReaderType::New();
551 defReader->SetFileName(m_ArgsInfo.init_arg);
553 multiResolutionFilter->SetArbitraryInitialDeformationField(defReader->GetOutput());
557 // =======================================================
559 // =======================================================
562 multiResolutionFilter->Update();
564 catch( itk::ExceptionObject & excp )
566 std::cerr <<"Error executing the demons filter: "<< excp << std::endl;
571 //=======================================================
573 //=======================================================
574 typename DeformationFieldType::Pointer deformationField=multiResolutionFilter->GetOutput();
577 //=======================================================
579 //=======================================================
580 typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
581 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
582 fieldWriter->SetInput( deformationField );
583 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
586 fieldWriter->Update();
588 catch( itk::ExceptionObject & excp )
590 std::cerr << "Exception thrown writing the DVF" << std::endl;
591 std::cerr << excp << std::endl;
596 //=======================================================
597 // Warp the moving image
598 //=======================================================
599 typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType;
600 typename WarpFilterType::Pointer warp = WarpFilterType::New();
602 warp->SetDeformationField( deformationField );
603 warp->SetInput( movingImageReader->GetOutput() );
604 warp->SetOutputOrigin( fixedImage->GetOrigin() );
605 warp->SetOutputSpacing( fixedImage->GetSpacing() );
606 warp->SetOutputDirection( fixedImage->GetDirection() );
607 warp->SetEdgePaddingValue( 0.0 );
611 //=======================================================
612 // Write the warped image
613 //=======================================================
614 typedef itk::ImageFileWriter< FixedImageType > WriterType;
615 typename WriterType::Pointer writer = WriterType::New();
616 writer->SetFileName( m_ArgsInfo.output_arg );
617 writer->SetInput( warp->GetOutput() );
623 catch( itk::ExceptionObject & err )
625 std::cerr << "ExceptionObject caught !" << std::endl;
626 std::cerr << err << std::endl;
631 //=======================================================
632 // Calculate the difference after the deformable transform
633 //=======================================================
634 typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType;
635 if (m_ArgsInfo.after_given)
637 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
638 difference->SetValidInput( fixedImage );
639 difference->SetTestInput( warp->GetOutput() );
643 difference->Update();
645 catch( itk::ExceptionObject & err )
647 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
648 std::cerr << err << std::endl;
652 typename WriterType::Pointer differenceWriter=WriterType::New();
653 differenceWriter->SetInput(difference->GetOutput());
654 differenceWriter->SetFileName(m_ArgsInfo.after_arg);
655 differenceWriter->Update();
660 //=======================================================
661 // Calculate the difference before the deformable transform
662 //=======================================================
663 if( m_ArgsInfo.before_given )
665 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
666 difference->SetValidInput( fixedImage );
667 difference->SetTestInput( movingImage );
671 difference->Update();
673 catch( itk::ExceptionObject & err )
675 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
676 std::cerr << err << std::endl;
680 typename WriterType::Pointer differenceWriter=WriterType::New();
681 writer->SetFileName( m_ArgsInfo.before_arg );
682 writer->SetInput( difference->GetOutput() );
691 #endif // __clitkDemonsDeformableRegistrationGenericFilter_txx