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 __clitkBSplineDeformableRegistrationGenericFilter_txx
19 #define __clitkBSplineDeformableRegistrationGenericFilter_txx
20 #include "clitkBSplineDeformableRegistrationGenericFilter.h"
26 //==============================================================================
27 //Creating an observer class that allows us to change parameters at subsequent levels
28 //==============================================================================
29 template <typename TRegistration>
30 class RegistrationInterfaceCommand : public itk::Command
33 typedef RegistrationInterfaceCommand Self;
36 typedef itk::Command Superclass;
37 typedef itk::SmartPointer<Self> Pointer;
40 RegistrationInterfaceCommand() {};
42 typedef TRegistration RegistrationType;
43 typedef RegistrationType * RegistrationPointer;
45 // Two arguments are passed to the Execute() method: the first
46 // is the pointer to the object which invoked the event and the
47 // second is the event that was invoked.
48 void Execute(itk::Object * object, const itk::EventObject & event)
50 if( !(itk::IterationEvent().CheckEvent( &event )) )
54 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
55 unsigned int numberOfLevels=registration->GetNumberOfLevels();
56 unsigned int currentLevel=registration->GetCurrentLevel()+1;
58 std::cout<<"========================================"<<std::endl;
59 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
60 std::cout<<"========================================"<<std::endl;
64 void Execute(const itk::Object * , const itk::EventObject & )
70 //==============================================================================
71 // Creating an observer class that allows output at each iteration
72 //==============================================================================
73 class CommandIterationUpdate : public itk::Command
76 typedef CommandIterationUpdate Self;
77 typedef itk::Command Superclass;
78 typedef itk::SmartPointer<Self> Pointer;
81 CommandIterationUpdate() {};
83 typedef clitk::GenericOptimizer<args_info_clitkBSplineDeformableRegistration> OptimizerType;
84 typedef const OptimizerType * OptimizerPointer;
86 // We set the generic optimizer
87 void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
90 void Execute(itk::Object *caller, const itk::EventObject & event)
92 Execute( (const itk::Object *)caller, event);
95 void Execute(const itk::Object * object, const itk::EventObject & event)
97 if( !(itk::IterationEvent().CheckEvent( &event )) )
102 m_Optimizer->OutputIterationInfo();
105 OptimizerPointer m_Optimizer;
109 //==============================================================================
110 // Update with the number of dimensions
111 //==============================================================================
112 template<unsigned int Dimension>
113 void BSplineDeformableRegistrationGenericFilter::UpdateWithDim(std::string PixelType)
116 if (m_Verbose) std::cout << "Images were detected to be "<< Dimension << "D and " << PixelType << "..." << std::endl;
118 if(PixelType == "short"){
119 if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed short..." << std::endl;
120 UpdateWithDimAndPixelType<Dimension, signed short>();
122 // else if(PixelType == "unsigned_short"){
123 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_short..." << std::endl;
124 // UpdateWithDimAndPixelType<Dimension, unsigned short>();
127 // else if (PixelType == "unsigned_char"){
128 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and unsigned_char..." << std::endl;
129 // UpdateWithDimAndPixelType<Dimension, unsigned char>();
132 // else if (PixelType == "char"){
133 // if (m_Verbose) std::cout << "Launching warp in "<< Dimension <<"D and signed_char..." << std::endl;
134 // UpdateWithDimAndPixelType<Dimension, signed char>();
137 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
138 UpdateWithDimAndPixelType<Dimension, float>();
144 //==============================================================================
145 // Update with the number of dimensions and pixeltype
146 //==============================================================================
147 template<unsigned int ImageDimension, class PixelType>
148 void BSplineDeformableRegistrationGenericFilter::UpdateWithDimAndPixelType()
151 //=======================================================
153 //=======================================================
154 bool threadsGiven=m_ArgsInfo.threads_given;
155 int threads=m_ArgsInfo.threads_arg;
157 typedef itk::Image< PixelType, ImageDimension > FixedImageType;
158 typedef itk::Image< PixelType, ImageDimension > MovingImageType;
159 const unsigned int SpaceDimension = ImageDimension;
160 typedef double TCoordRep;
163 //=======================================================
165 //=======================================================
166 typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
167 typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
169 typename FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
170 typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
172 fixedImageReader->SetFileName( m_ArgsInfo.reference_arg );
173 movingImageReader->SetFileName( m_ArgsInfo.target_arg );
174 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
175 fixedImageReader->Update();
176 movingImageReader->Update();
178 typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
179 typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
180 typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
182 // The metric region: where should the metric be CALCULATED (depends on mask)
183 typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
184 typename FixedImageType::RegionType::SizeType metricRegionSize=metricRegion.GetSize();
185 typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
186 typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
188 // The transform region: where should the transform be DEFINED (depends on mask)
189 typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
190 typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
191 typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
192 typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
195 //=======================================================
196 // If given, we connect a mask to the fixed image
197 //======================================================
198 typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
199 typename MaskType::Pointer spatialObjectMask=NULL;
201 if (m_ArgsInfo.mask_given)
203 typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
204 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
205 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
206 maskReader->SetFileName(m_ArgsInfo.mask_arg);
210 maskReader->Update();
212 catch( itk::ExceptionObject & err )
214 std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
215 std::cerr << err << std::endl;
218 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
221 // Set the image to the spatialObject
222 spatialObjectMask = MaskType::New();
223 spatialObjectMask->SetImage( maskReader->GetOutput() );
225 // Find the bounding box of the "inside" label
226 typedef itk::LabelStatisticsImageFilter<ImageMaskType, ImageMaskType> StatisticsImageFilterType;
227 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
228 statisticsImageFilter->SetInput(maskReader->GetOutput());
229 statisticsImageFilter->SetLabelInput(maskReader->GetOutput());
230 statisticsImageFilter->Update();
231 typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
233 // Limit the transform region to the mask
234 for (unsigned int i=0; i<ImageDimension; i++)
236 transformRegionIndex[i]=boundingBox[2*i];
237 transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
239 transformRegion.SetSize(transformRegionSize);
240 transformRegion.SetIndex(transformRegionIndex);
241 fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
243 // Limit the metric region to the mask
244 metricRegion=transformRegion;
245 fixedImage->TransformIndexToPhysicalPoint(metricRegion.GetIndex(), metricRegionOrigin);
251 //=======================================================
253 //=====================================================
256 // Fixed image region
257 std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
258 <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
259 <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
262 std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
263 <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
264 <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
267 std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
268 <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
269 <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
274 //=======================================================
276 //=======================================================
277 typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
278 typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
279 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
280 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
283 // //=======================================================
284 // // Rigid Transform
285 // //=======================================================
286 // typedef itk::Euler3DTransform <double> RigidTransformType;
287 // RigidTransformType::Pointer rigidTransform=RigidTransformType::New();
289 // if (m_ArgsInfo.rigid_given)
291 // itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.rigid_arg);
293 // //Set the rotation
294 // itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
295 // rigidTransform->SetMatrix(finalRotation);
297 // //Set the translation
298 // itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
299 // rigidTransform->SetTranslation(finalTranslation);
304 //=======================================================
306 //=======================================================
307 typename FixedImageType::RegionType::SizeType splineOrders ;
309 //Default is cubic splines
310 splineOrders.Fill(3);
311 if (m_ArgsInfo.order_given)
312 for(unsigned int i=0; i<ImageDimension;i++)
313 splineOrders[i]=m_ArgsInfo.order_arg[i];
316 typedef itk::Transform<TCoordRep, ImageDimension, SpaceDimension> TransformType;
317 typename TransformType::Pointer transform;
318 typedef itk::BSplineDeformableTransform<TCoordRep,SpaceDimension, 3> BSplineTransformType;
319 typedef BSplineTransformType* BSplineTransformPointer;
320 typedef clitk::BSplineDeformableTransform<TCoordRep,ImageDimension, SpaceDimension > BLUTTransformType;
321 typedef BLUTTransformType* BLUTTransformPointer;
323 // JV parameter array is passed by reference, create outside context so it exists afterwards!!!!!
324 typedef typename TransformType::ParametersType ParametersType;
325 ParametersType parameters;
328 // CLITK BLUT transform
329 if(m_ArgsInfo.wlut_flag)
331 typename BLUTTransformType::Pointer bsplineTransform = BLUTTransformType::New();
332 if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
333 bsplineTransform->SetSplineOrders(splineOrders);
335 //-------------------------------------------------------------------------
336 // Define the region: Either the spacing or the number of CP should be given
337 //-------------------------------------------------------------------------
340 typedef typename BSplineTransformType::RegionType RegionType;
341 RegionType bsplineRegion;
342 typename RegionType::SizeType gridSizeOnImage;
343 typename RegionType::SizeType gridBorderSize;
344 typename RegionType::SizeType totalGridSize;
347 typedef typename BSplineTransformType::SpacingType SpacingType;
348 SpacingType fixedImageSpacing, chosenSpacing, adaptedSpacing;
349 fixedImageSpacing = fixedImage->GetSpacing();
351 // Only spacing given: adjust if necessary
352 if (m_ArgsInfo.spacing_given && !m_ArgsInfo.control_given)
354 for(unsigned int r=0; r<ImageDimension; r++)
356 chosenSpacing[r]= m_ArgsInfo.spacing_arg[r];
357 gridSizeOnImage[r] = ceil( (double) transformRegion.GetSize()[r] / ( round(chosenSpacing[r]/fixedImageSpacing[r]) ) );
358 adaptedSpacing[r]= ( round(chosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
360 if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
361 if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
362 if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
365 // Only number of CP given: adjust if necessary
366 else if (m_ArgsInfo.control_given && !m_ArgsInfo.spacing_given)
368 for(unsigned int r=0; r<ImageDimension; r++)
370 gridSizeOnImage[r]= m_ArgsInfo.control_arg[r];
371 chosenSpacing[r]=fixedImageSpacing[r]*( (double)(transformRegion.GetSize()[r]) /
372 (double)(gridSizeOnImage[r]) );
373 adaptedSpacing[r]= fixedImageSpacing[r]* ceil( (double)(transformRegion.GetSize()[r] - 1) /
374 (double)(gridSizeOnImage[r] - 1) );
376 if (m_Verbose) std::cout<<"The chosen control point spacing "<<chosenSpacing<<"..."<<std::endl;
377 if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<adaptedSpacing<<"..."<<std::endl;
378 if (m_Verbose) std::cout<<"The number of (internal) control points is "<<gridSizeOnImage<<"..."<<std::endl;
381 // Spacing and number of CP given: no adjustment adjust, just warnings
382 else if (m_ArgsInfo.control_given && m_ArgsInfo.spacing_given)
384 for(unsigned int r=0; r<ImageDimension; r++)
386 adaptedSpacing[r]= m_ArgsInfo.spacing_arg[r];
387 gridSizeOnImage[r] = m_ArgsInfo.control_arg[r];
388 if (gridSizeOnImage[r]*adaptedSpacing[r]< transformRegion.GetSize()[r]*fixedImageSpacing[r])
390 std::cout<<"WARNING: Specified control point region ("<<gridSizeOnImage[r]*adaptedSpacing[r]
391 <<"mm) does not cover the transform region ("<< transformRegion.GetSize()[r]*fixedImageSpacing[r]
392 <<"mm) for dimension "<<r<<"!" <<std::endl
393 <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
395 if ( fmod(adaptedSpacing[r], fixedImageSpacing[r]) )
397 std::cout<<"WARNING: Specified control point spacing for dimension "<<r
398 <<" does not allow exact representation of BLUT FFD!"<<std::endl
399 <<"Spacing ratio is non-integer: "<<adaptedSpacing[r]/ fixedImageSpacing[r]<<std::endl
400 <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
403 if (m_Verbose) std::cout<<"The control points spacing was set to "<<adaptedSpacing<<"..."<<std::endl;
404 if (m_Verbose) std::cout<<"The number of (internal) control points spacing is "<<gridSizeOnImage<<"..."<<std::endl;
407 //JV border size should depend on spline order
408 for(unsigned int r=0; r<ImageDimension; r++) gridBorderSize[r]=splineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
409 totalGridSize = gridSizeOnImage + gridBorderSize;
410 bsplineRegion.SetSize( totalGridSize );
411 if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
414 typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
415 SpacingType gridOriginOffset = gridDirection * adaptedSpacing;
417 // Origin: 1 CP border for spatial dimensions
418 typedef typename BSplineTransformType::OriginType OriginType;
419 OriginType gridOrigin = transformRegionOrigin - gridOriginOffset;
420 if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;
423 bsplineTransform->SetGridSpacing( adaptedSpacing );
424 bsplineTransform->SetGridOrigin( gridOrigin );
425 bsplineTransform->SetGridRegion( bsplineRegion );
426 bsplineTransform->SetGridDirection( gridDirection );
429 //if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
430 //bsplineTransform->SetBulkTransform( rigidTransform );
432 //Vector BSpline interpolator
433 //bsplineTransform->SetOutputSpacing(fixedImage->GetSpacing());
434 typename RegionType::SizeType samplingFactors;
435 for (unsigned int i=0; i< ImageDimension; i++)
437 if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
438 << (double)adaptedSpacing[i]/ (double) fixedImageSpacing[i]<<"..."<<std::endl;
439 if (m_ArgsInfo.samplingFactor_given) samplingFactors[i]=m_ArgsInfo.samplingFactor_arg[i];
440 else samplingFactors[i]=(int) ((double)adaptedSpacing[i]/ (double) movingImage->GetSpacing()[i]);
441 if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<samplingFactors[i]<<"..."<<std::endl;
443 bsplineTransform->SetLUTSamplingFactors(samplingFactors);
446 if (m_ArgsInfo.init_given)
448 typedef itk::ImageFileReader<typename BLUTTransformType::CoefficientImageType> CoefficientReaderType;
449 typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
450 coeffReader->SetFileName(m_ArgsInfo.init_arg[0]);
451 coeffReader->Update();
452 bsplineTransform->SetCoefficientImage(coeffReader->GetOutput());
456 //typedef typename TransformType::ParametersType ParametersType;
457 const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters();
458 parameters=ParametersType( numberOfParameters );
459 parameters.Fill( 0.0 );
460 bsplineTransform->SetParameters( parameters );
464 if (spatialObjectMask) bsplineTransform->SetMask( spatialObjectMask );
467 transform=bsplineTransform;
470 //ITK BSpline transform
473 typename BSplineTransformType::Pointer bsplineTransform = BSplineTransformType::New();
476 typedef typename BSplineTransformType::RegionType RegionType;
477 RegionType bsplineRegion;
478 typename RegionType::SizeType gridSizeOnImage;
479 typename RegionType::SizeType gridBorderSize;
480 typename RegionType::SizeType totalGridSize;
482 //Set the number of control points
483 for(unsigned int r=0; r<ImageDimension; r++) gridSizeOnImage[r]=m_ArgsInfo.control_arg[r];
484 if (m_Verbose) std::cout<<"Setting the number of internal control points "<<gridSizeOnImage<<"..."<<std::endl;
485 gridBorderSize.Fill( 3 ); // Border for spline order = 3 ( 1 lower, 2 upper )
486 totalGridSize = gridSizeOnImage + gridBorderSize;
487 bsplineRegion.SetSize( totalGridSize );
490 typedef typename BSplineTransformType::SpacingType SpacingType;
491 SpacingType spacing = fixedImage->GetSpacing();
492 typename FixedImageType::SizeType fixedImageSize = fixedImageRegion.GetSize();
493 if (m_ArgsInfo.spacing_given)
496 for(unsigned int r=0; r<ImageDimension; r++)
498 spacing[r] =m_ArgsInfo.spacing_arg[r];
503 for(unsigned int r=0; r<ImageDimension; r++)
505 spacing[r] *= static_cast<double>(fixedImageSize[r] - 1) /
506 static_cast<double>(gridSizeOnImage[r] - 1);
509 if (m_Verbose) std::cout<<"The control points spacing was set to "<<spacing<<"..."<<std::endl;
512 typename FixedImageType::DirectionType gridDirection = fixedImage->GetDirection();
513 SpacingType gridOriginOffset = gridDirection * spacing;
516 typedef typename BSplineTransformType::OriginType OriginType;
517 OriginType origin = fixedImage->GetOrigin();
518 OriginType gridOrigin = origin - gridOriginOffset;
521 bsplineTransform->SetGridSpacing( spacing );
522 bsplineTransform->SetGridOrigin( gridOrigin );
523 bsplineTransform->SetGridRegion( bsplineRegion );
524 bsplineTransform->SetGridDirection( gridDirection );
527 // if (m_Verbose) std::cout<<"Setting rigid transform..."<<std::endl;
528 // bsplineTransform->SetBulkTransform( rigidTransform );
530 // Initial parameters
531 if (m_ArgsInfo.init_given)
533 typedef itk::ImageFileReader<typename BSplineTransformType::ImageType> CoefficientReaderType;
534 typename BSplineTransformType::ImageType::Pointer coeffImages[SpaceDimension];
535 for(unsigned int i=0; i<SpaceDimension; i++)
537 typename CoefficientReaderType::Pointer coeffReader=CoefficientReaderType::New();
538 coeffReader->SetFileName(m_ArgsInfo.init_arg[i]);
539 coeffReader->Update();
540 coeffImages[i]=coeffReader->GetOutput();
542 bsplineTransform->SetCoefficientImage(coeffImages);
546 const unsigned int numberOfParameters = bsplineTransform->GetNumberOfParameters();
547 parameters=ParametersType( numberOfParameters );
548 parameters.Fill( 0.0 );
549 bsplineTransform->SetParameters( parameters );
553 transform=bsplineTransform;
558 //=======================================================
560 //=======================================================
561 typedef clitk::GenericInterpolator<args_info_clitkBSplineDeformableRegistration, FixedImageType, TCoordRep > GenericInterpolatorType;
562 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
563 genericInterpolator->SetArgsInfo(m_ArgsInfo);
564 typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
565 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
568 //=======================================================
570 //=======================================================
571 typedef clitk::GenericMetric<args_info_clitkBSplineDeformableRegistration, FixedImageType,MovingImageType > GenericMetricType;
572 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
573 genericMetric->SetArgsInfo(m_ArgsInfo);
574 genericMetric->SetFixedImageRegion(metricRegion);
575 typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
576 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
577 if (spatialObjectMask) metric->SetFixedImageMask( spatialObjectMask );
579 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
580 if (threadsGiven) metric->SetNumberOfThreads( threads );
582 if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
586 //=======================================================
588 //=======================================================
589 typedef clitk::GenericOptimizer<args_info_clitkBSplineDeformableRegistration> GenericOptimizerType;
590 GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
591 genericOptimizer->SetArgsInfo(m_ArgsInfo);
592 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
593 genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters());
594 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
595 OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
598 //=======================================================
600 //=======================================================
601 typedef itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;
602 typename RegistrationType::Pointer registration = RegistrationType::New();
603 registration->SetMetric( metric );
604 registration->SetOptimizer( optimizer );
605 registration->SetInterpolator( interpolator );
606 registration->SetTransform (transform);
607 if(threadsGiven) registration->SetNumberOfThreads(threads);
608 registration->SetFixedImage( fixedImage );
609 registration->SetMovingImage( movingImage );
610 registration->SetFixedImageRegion( metricRegion );
611 registration->SetFixedImagePyramid( fixedImagePyramid );
612 registration->SetMovingImagePyramid( movingImagePyramid );
613 registration->SetInitialTransformParameters( transform->GetParameters() );
614 registration->SetNumberOfLevels(m_ArgsInfo.levels_arg);
615 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
618 //================================================================================================
620 //================================================================================================
623 // Output iteration info
624 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
625 observer->SetOptimizer(genericOptimizer);
626 optimizer->AddObserver( itk::IterationEvent(), observer );
630 typedef RegistrationInterfaceCommand<RegistrationType> CommandType;
631 typename CommandType::Pointer command = CommandType::New();
632 registration->AddObserver( itk::IterationEvent(), command );
636 //=======================================================
638 //=======================================================
639 if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
643 registration->StartRegistration();
645 catch( itk::ExceptionObject & err )
647 std::cerr << "ExceptionObject caught while registering!" << std::endl;
648 std::cerr << err << std::endl;
653 //=======================================================
655 //=======================================================
656 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
657 transform->SetParameters( finalParameters );
660 //=======================================================
661 // Get the BSpline coefficient images and write them
662 //=======================================================
663 if (m_ArgsInfo.coeff_given)
665 if(m_ArgsInfo.wlut_flag)
667 BLUTTransformPointer bsplineTransform=dynamic_cast<BLUTTransformPointer>(registration->GetTransform());
668 typedef itk::Image<itk::Vector<TCoordRep, SpaceDimension>, ImageDimension> CoefficientImageType;
669 typename CoefficientImageType::Pointer coefficientImage =bsplineTransform->GetCoefficientImage();
670 typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
671 typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
672 coeffWriter->SetInput(coefficientImage);
673 coeffWriter->SetFileName(m_ArgsInfo.coeff_arg[0]);
674 coeffWriter->Update();
678 BSplineTransformPointer bsplineTransform=dynamic_cast<BSplineTransformPointer>(registration->GetTransform());
679 typedef itk::Image<TCoordRep, ImageDimension> CoefficientImageType;
680 #if ITK_VERSION_MAJOR > 3
681 typename BSplineTransformType::CoefficientImageArray coefficientImages = bsplineTransform->GetCoefficientImage();
683 typename CoefficientImageType::Pointer *coefficientImages =bsplineTransform->GetCoefficientImage();
685 typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
686 for (unsigned int i=0;i<SpaceDimension; i ++)
688 typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
689 coeffWriter->SetInput(coefficientImages[i]);
690 coeffWriter->SetFileName(m_ArgsInfo.coeff_arg[i]);
691 coeffWriter->Update();
697 //=======================================================
699 //=======================================================
700 typedef itk::Vector< float, SpaceDimension > DisplacementType;
701 typedef itk::Image< DisplacementType, ImageDimension > DeformationFieldType;
703 typename DeformationFieldType::Pointer field = DeformationFieldType::New();
704 field->SetRegions( fixedImageRegion );
705 field->SetOrigin( fixedImage->GetOrigin() );
706 field->SetSpacing( fixedImage->GetSpacing() );
707 field->SetDirection( fixedImage->GetDirection() );
710 typedef itk::ImageRegionIteratorWithIndex< DeformationFieldType > FieldIterator;
711 FieldIterator fi( field, fixedImageRegion );
714 typename TransformType::InputPointType fixedPoint;
715 typename TransformType::OutputPointType movingPoint;
716 typename DeformationFieldType::IndexType index;
718 DisplacementType displacement;
719 while( ! fi.IsAtEnd() )
721 index = fi.GetIndex();
722 field->TransformIndexToPhysicalPoint( index, fixedPoint );
723 movingPoint = transform->TransformPoint( fixedPoint );
724 displacement = movingPoint - fixedPoint;
725 fi.Set( displacement );
730 //=======================================================
732 //=======================================================
733 typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
734 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
735 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
736 fieldWriter->SetInput( field );
739 fieldWriter->Update();
741 catch( itk::ExceptionObject & excp )
743 std::cerr << "Exception thrown writing the DVF" << std::endl;
744 std::cerr << excp << std::endl;
749 //=======================================================
750 // Resample the moving image
751 //=======================================================
752 typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType;
753 typename WarpFilterType::Pointer warp = WarpFilterType::New();
755 warp->SetDeformationField( field );
756 warp->SetInput( movingImageReader->GetOutput() );
757 warp->SetOutputOrigin( fixedImage->GetOrigin() );
758 warp->SetOutputSpacing( fixedImage->GetSpacing() );
759 warp->SetOutputDirection( fixedImage->GetDirection() );
760 warp->SetEdgePaddingValue( 0.0 );
764 //=======================================================
765 // Write the warped image
766 //=======================================================
767 typedef itk::ImageFileWriter< FixedImageType > WriterType;
768 typename WriterType::Pointer writer = WriterType::New();
769 writer->SetFileName( m_ArgsInfo.output_arg );
770 writer->SetInput( warp->GetOutput() );
776 catch( itk::ExceptionObject & err )
778 std::cerr << "ExceptionObject caught !" << std::endl;
779 std::cerr << err << std::endl;
784 //=======================================================
785 // Calculate the difference after the deformable transform
786 //=======================================================
787 typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType;
788 if (m_ArgsInfo.after_given)
790 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
791 difference->SetValidInput( fixedImage );
792 difference->SetTestInput( warp->GetOutput() );
796 difference->Update();
798 catch( itk::ExceptionObject & err )
800 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
801 std::cerr << err << std::endl;
805 typename WriterType::Pointer differenceWriter=WriterType::New();
806 differenceWriter->SetInput(difference->GetOutput());
807 differenceWriter->SetFileName(m_ArgsInfo.after_arg);
808 differenceWriter->Update();
813 //=======================================================
814 // Calculate the difference before the deformable transform
815 //=======================================================
816 if( m_ArgsInfo.before_given )
819 typename FixedImageType::Pointer moving=FixedImageType::New();
820 if (m_ArgsInfo.rigid_given)
822 typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
823 typename ResamplerType::Pointer resampler=ResamplerType::New();
824 resampler->SetInput(movingImage);
825 resampler->SetOutputOrigin(fixedImage->GetOrigin());
826 resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
827 resampler->SetOutputSpacing(fixedImage->GetSpacing());
828 resampler->SetDefaultPixelValue( 0. );
829 //resampler->SetTransform(rigidTransform);
831 moving=resampler->GetOutput();
836 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
837 difference->SetValidInput( fixedImage );
838 difference->SetTestInput( moving );
842 difference->Update();
844 catch( itk::ExceptionObject & err )
846 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
847 std::cerr << err << std::endl;
851 typename WriterType::Pointer differenceWriter=WriterType::New();
852 writer->SetFileName( m_ArgsInfo.before_arg );
853 writer->SetInput( difference->GetOutput() );
861 #endif // __clitkBSplineDeformableRegistrationGenericFilter_txx