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 #if ITK_VERSION_MAJOR >= 4
155 typedef typename RegistrationFilterType::DisplacementFieldType DisplacementFieldType;
157 typedef typename RegistrationFilterType::DeformationFieldType DisplacementFieldType;
159 typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType, MovingImageType, DisplacementFieldType> MultiResolutionRegistrationType;
160 typedef CommandResolutionLevelUpdate<MultiResolutionRegistrationType> LevelObserver;
163 CommandIterationUpdate(){};
168 void SetStop( int* s){m_Stop=s;}
169 void SetLevelObserver(LevelObserver* o ){m_LevelObserver=o;}
173 void Execute(const itk::Object *, const itk::EventObject & )
175 std::cout << "Warning: The const Execute method shouldn't be called" << std::endl;
178 void Execute(itk::Object *caller, const itk::EventObject & event)
180 if( !(itk::IterationEvent().CheckEvent( &event )) )
186 RegistrationFilterType * filter = dynamic_cast< RegistrationFilterType * >( caller );
191 m_Iteration=filter->GetElapsedIterations();
192 m_Metric=filter->GetMetric();
195 std::cout << m_Iteration<<"\t Field Tolerance= "<<filter->GetMaximumRMSError();
196 std::cout <<"\t DVF Change= " << filter->GetRMSChange()<<"\t RMS= "<<m_Metric;
198 // Using stop criterion?
199 if (m_Stop[m_LevelObserver->m_CurrentLevel]>=0)
209 else if(m_Metric<m_Minimum)
215 //Not less then minimum
219 if (m_StopCounter>=m_Stop[m_LevelObserver->m_CurrentLevel])
220 filter->StopRegistration();
224 std::cout <<"\t Stop= "<<m_StopCounter<<" / "<<m_Stop[m_LevelObserver->m_CurrentLevel]<<std::endl;
227 // Not using stop criterion
228 else std::cout <<std::endl;
232 double m_Minimum, m_Metric;
233 int m_StopCounter, m_Iteration;
235 LevelObserver* m_LevelObserver;
239 //==============================================================================
240 // Update with the number of dimensions
241 //==============================================================================
242 template<unsigned int Dimension>
243 void DemonsDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
245 if(PixelType == "short"){
246 if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl;
247 UpdateWithDimAndPixelType<Dimension, signed short>();
249 // else if(PixelType == "unsigned_short"){
250 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl;
251 // UpdateWithDimAndPixelType<Dimension, unsigned short>();
254 // else if (PixelType == "unsigned_char"){
255 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl;
256 // UpdateWithDimAndPixelType<Dimension, unsigned char>();
259 // else if (PixelType == "char"){
260 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl;
261 // UpdateWithDimAndPixelType<Dimension, signed char>();
264 if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and float..." << std::endl;
265 UpdateWithDimAndPixelType<Dimension, float>();
271 //==============================================================================
272 // Update with the number of dimensions and pixeltype
273 //==============================================================================
274 //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
276 template<unsigned int ImageDimension, class PixelType>
277 void DemonsDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType()
279 //=======================================================
281 //=======================================================
282 bool threadsGiven=m_ArgsInfo.threads_given;
283 int threads=m_ArgsInfo.threads_arg;
285 typedef itk::Image< PixelType, ImageDimension > FixedImageType;
286 typedef itk::Image< PixelType, ImageDimension > MovingImageType;
287 typedef itk::Image< float, ImageDimension > InternalImageType;
288 typedef double TCoordRep;
291 //=======================================================
293 //=======================================================
294 typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
295 typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
297 typename FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
298 typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
300 fixedImageReader->SetFileName( m_ArgsInfo.reference_arg );
301 movingImageReader->SetFileName( m_ArgsInfo.target_arg );
302 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
303 fixedImageReader->Update();
304 movingImageReader->Update();
306 typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
307 typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
308 typename FixedImageType::RegionType fixedRegion = fixedImage->GetBufferedRegion();
311 //=======================================================
313 //=======================================================
314 typedef itk::Vector< float, ImageDimension > VectorPixelType;
315 typedef itk::Image< VectorPixelType, ImageDimension > DeformationFieldType;
318 //=======================================================
320 //=======================================================
321 typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, InternalImageType> FixedImagePyramidType;
322 typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, InternalImageType> MovingImagePyramidType;
323 // typedef itk::MultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
324 // typedef itk::MultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
325 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
326 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
329 // =======================================================
330 // Define the registation filters
331 // =======================================================
332 typedef itk::PDEDeformableRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType> RegistrationFilterType;
333 typename RegistrationFilterType::Pointer pdeFilter;
334 typedef clitk::MultiResolutionPDEDeformableRegistration<FixedImageType,MovingImageType,DeformationFieldType > MultiResolutionRegistrationFilterType;
335 typename MultiResolutionRegistrationFilterType::Pointer multiResolutionFilter = MultiResolutionRegistrationFilterType::New();
336 typedef itk::ESMDemonsRegistrationFunction<InternalImageType, InternalImageType, DeformationFieldType> RegistrationFunctionType;
339 // =======================================================
340 // The multiresolution scheme
341 // =======================================================
342 if (threadsGiven) multiResolutionFilter->SetNumberOfThreads(threads);
343 unsigned int nLevels=m_ArgsInfo.levels_arg;
344 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<nLevels<<"..."<<std::endl;
345 multiResolutionFilter->SetFixedImage( fixedImage );
346 multiResolutionFilter->SetMovingImage( movingImage );
347 multiResolutionFilter->SetNumberOfLevels( nLevels );
348 multiResolutionFilter->SetFixedImagePyramid( fixedImagePyramid );
349 multiResolutionFilter->SetMovingImagePyramid( movingImagePyramid );
350 if (threadsGiven) multiResolutionFilter->SetNumberOfThreads( threads );
352 //------------------------------------
353 //Set the number of iterations
354 //------------------------------------
355 std::vector<unsigned int> nIterations(nLevels);
356 for (unsigned int i=0 ; i<nLevels; i++)
358 if (m_ArgsInfo.maxIter_given==nLevels)
360 nIterations[i] =m_ArgsInfo.maxIter_arg[i];
363 nIterations[i]=m_ArgsInfo.maxIter_arg[0];
365 multiResolutionFilter->SetNumberOfIterations( &(nIterations[0]) );
367 std::cout<<"Setting the number of iterations to: "<<nIterations[0];
368 for (unsigned int i=1; i<nLevels; i++)
369 std::cout<<", "<<nIterations [i];
370 std::cout<<std::endl;
373 //------------------------------------
374 //Set the max RMS error for the field update
375 //------------------------------------
376 std::vector<double> maxRMSError(nLevels);
377 for (unsigned int i=0 ; i<nLevels; i++)
379 if (m_ArgsInfo.maxRMSError_given==nLevels)
380 maxRMSError[i] =m_ArgsInfo.maxRMSError_arg[i];
382 maxRMSError[i]=m_ArgsInfo.maxRMSError_arg[0];
385 std::cout<<"Setting the max root mean squared error to: "<<maxRMSError[0];
386 for (unsigned int i=1; i<nLevels; i++)
387 std::cout<<", "<<maxRMSError[i];
388 std::cout<<std::endl;
391 //------------------------------------
392 //Get the stop criterion
393 //------------------------------------
394 std::vector<int> stop(nLevels);
395 for (unsigned int i=0; i<nLevels; i++)
396 if (m_ArgsInfo.stop_given==nLevels)
397 stop[i]=m_ArgsInfo.stop_arg[i];
399 stop[i]=m_ArgsInfo.stop_arg[0];
401 std::cout<<"Setting the stop criterion to : "<<stop[0];
402 for (unsigned int i=1; i<nLevels; i++)
403 std::cout<<", "<<stop[i];
404 std::cout<<std::endl;
407 //------------------------------------
409 //------------------------------------
410 typename RegistrationFunctionType::GradientType grad=RegistrationFunctionType::Symmetric;
411 switch (m_ArgsInfo.gradType_arg)
414 grad= RegistrationFunctionType::Symmetric;
417 grad= RegistrationFunctionType::Fixed;
420 grad= RegistrationFunctionType::WarpedMoving;
423 grad= RegistrationFunctionType::MappedMoving;
427 //------------------------------------
428 // Set smoothing standard deviations
429 //------------------------------------
430 double sd [ImageDimension];
431 for (unsigned int i=0; i<ImageDimension; i++)
432 if (m_ArgsInfo.sd_given==ImageDimension)
433 sd[i]=m_ArgsInfo.sd_arg[i];
435 sd[i]=m_ArgsInfo.sd_arg[0];
437 std::cout<<"Setting the final standard deviations for smoothing to: "<<sd[0];
438 for (unsigned int i=1; i<ImageDimension; i++)
439 std::cout<<", "<<sd[i];
440 std::cout<<std::endl;
443 // =======================================================
444 // The level observer: adjust settings at each level
445 // =======================================================
446 typename CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::Pointer levelObserver =
447 CommandResolutionLevelUpdate<MultiResolutionRegistrationFilterType>::New();
448 multiResolutionFilter->AddObserver( itk::IterationEvent(), levelObserver );
449 levelObserver->SetMaxRMSError(&(maxRMSError[0]));
450 levelObserver->SetMaxStep(m_ArgsInfo.maxStep_arg);
451 levelObserver->SetSD(sd);
452 levelObserver->SetScaleStep(m_ArgsInfo.scaleStep_flag);
453 levelObserver->SetScaleSD(m_ArgsInfo.scaleSD_flag);
454 levelObserver->SetRegistrationType(m_ArgsInfo.demons_arg);
457 // =======================================================
458 // The type of filter
459 // =======================================================
460 switch (m_ArgsInfo.demons_arg){
464 typedef itk::DemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
465 typename DemonsFilterType::Pointer m =DemonsFilterType::New();
467 //Set Parameters for this filter
468 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
469 m->SetUseMovingImageGradient( m_ArgsInfo.movGrad_flag);
470 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
471 observer->SetStop(&(stop[0]));
472 observer->SetLevelObserver(levelObserver);
473 m->AddObserver( itk::IterationEvent(), observer );
474 if (m_Verbose) std::cout<<"Using the demons registration filter..."<<std::endl;
481 typedef itk::SymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
482 typename DemonsFilterType::Pointer m =DemonsFilterType::New();
484 //Set Parameters for this filter
485 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
486 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
487 observer->SetStop(&(stop[0]));
488 observer->SetLevelObserver(levelObserver);
489 m->AddObserver( itk::IterationEvent(), observer );
490 if (m_Verbose) std::cout<<"Using the symmetric forces demons registration filter..."<<std::endl;
497 typedef itk::FastSymmetricForcesDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
498 typename DemonsFilterType::Pointer m = DemonsFilterType::New();
500 //Set Parameters for this filter
501 m->SetIntensityDifferenceThreshold( m_ArgsInfo.intThreshold_arg);
502 m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
503 m->SetUseGradientType(grad);
504 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
505 observer->SetStop(&(stop[0]));
506 observer->SetLevelObserver(levelObserver);
507 m->AddObserver( itk::IterationEvent(), observer );
508 if (m_Verbose) std::cout<<"Using the fast symmetric forces demons registration filter..."<<std::endl;
515 typedef itk::DiffeomorphicDemonsRegistrationFilter< InternalImageType, InternalImageType, DeformationFieldType > DemonsFilterType;
516 typename DemonsFilterType::Pointer m = DemonsFilterType::New();
518 //Set Parameters for this filter
519 m->SetMaximumUpdateStepLength( m_ArgsInfo.maxStep_arg);
520 m->SetUseFirstOrderExp(m_ArgsInfo.firstOrder_flag);
521 m->SetUseGradientType(grad);
522 typename CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::Pointer observer = CommandIterationUpdate<DemonsFilterType, FixedImageType, MovingImageType>::New();
523 observer->SetStop(&(stop[0]));
524 observer->SetLevelObserver(levelObserver);
525 m->AddObserver( itk::IterationEvent(), observer );
526 if (m_Verbose) std::cout<<"Using the diffeomorphic demons registration filter..."<<std::endl;
535 pdeFilter->SetStandardDeviations( sd );
536 pdeFilter->SetUpdateFieldStandardDeviations( sd );
538 // pdeFilter->SetMaximumError(m_ArgsInfo.maxError_arg);
539 // pdeFilter->SetMaximumKernelWidth(m_ArgsInfo.maxError_arg);
540 #if ITK_VERSION_MAJOR >= 4
541 pdeFilter->SetSmoothDisplacementField(!m_ArgsInfo.fluid_flag);
543 pdeFilter->SetSmoothDeformationField(!m_ArgsInfo.fluid_flag);
545 pdeFilter->SetSmoothUpdateField(m_ArgsInfo.fluid_flag);
546 pdeFilter->SetUseImageSpacing( m_ArgsInfo.spacing_flag );
548 //Pass to the multi resolution scheme
549 multiResolutionFilter->SetRegistrationFilter( pdeFilter );
552 // =======================================================
553 // The initial solution
554 // =======================================================
555 if (m_ArgsInfo.init_given)
557 typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
558 typename DeformationFieldReaderType::Pointer defReader=DeformationFieldReaderType::New();
559 defReader->SetFileName(m_ArgsInfo.init_arg);
561 multiResolutionFilter->SetArbitraryInitialDeformationField(defReader->GetOutput());
565 // =======================================================
567 // =======================================================
570 multiResolutionFilter->Update();
572 catch( itk::ExceptionObject & excp )
574 std::cerr <<"Error executing the demons filter: "<< excp << std::endl;
579 //=======================================================
581 //=======================================================
582 typename DeformationFieldType::Pointer deformationField=multiResolutionFilter->GetOutput();
585 //=======================================================
587 //=======================================================
588 typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
589 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
590 fieldWriter->SetInput( deformationField );
591 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
594 fieldWriter->Update();
596 catch( itk::ExceptionObject & excp )
598 std::cerr << "Exception thrown writing the DVF" << std::endl;
599 std::cerr << excp << std::endl;
604 //=======================================================
605 // Warp the moving image
606 //=======================================================
607 typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType;
608 typename WarpFilterType::Pointer warp = WarpFilterType::New();
610 #if ITK_VERSION_MAJOR >= 4
611 warp->SetDisplacementField( deformationField );
613 warp->SetDeformationField( 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