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 clitkShapedBLUTSpatioTemporalDIRGenericFilter_txx
19 #define clitkShapedBLUTSpatioTemporalDIRGenericFilter_txx
21 /* =================================================
22 * @file clitkShapedBLUTSpatioTemporalDIRGenericFilter.txx
28 ===================================================*/
30 #include "clitkShapedBLUTSpatioTemporalDIRGenericFilter.h"
35 //==============================================================================
36 // Creating an observer class that allows output at each iteration
37 //==============================================================================
38 template <typename TRegistration, class args_info_clitkShapedBLUTSpatioTemporalDIR>
39 class CommandIterationUpdate : public itk::Command
42 typedef CommandIterationUpdate Self;
43 typedef itk::Command Superclass;
44 typedef itk::SmartPointer<Self> Pointer;
48 typedef TRegistration RegistrationType;
49 typedef RegistrationType * RegistrationPointer;
52 typedef typename RegistrationType::FixedImageType FixedImageType;
53 typedef typename FixedImageType::RegionType RegionType;
54 itkStaticConstMacro(ImageDimension, unsigned int,FixedImageType::ImageDimension);
55 typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform<double, ImageDimension, ImageDimension> TransformType;
56 typedef clitk::ShapedBLUTSpatioTemporalDeformableTransformInitializer<TransformType, FixedImageType> InitializerType;
57 typedef typename InitializerType::CoefficientImageType CoefficientImageType;
58 typedef itk::CastImageFilter<CoefficientImageType, CoefficientImageType> CastImageFilterType;
59 typedef typename TransformType::ParametersType ParametersType;
60 typedef typename InitializerType::Pointer InitializerPointer;
63 CommandIterationUpdate() { m_IterationCounter=0;}
65 typedef clitk::GenericOptimizer<args_info_clitkShapedBLUTSpatioTemporalDIR> OptimizerType;
66 typedef const OptimizerType * OptimizerPointer;
69 void Execute(itk::Object *caller, const itk::EventObject & event)
71 Execute( (const itk::Object *)caller, event);
74 void Execute(const itk::Object * object, const itk::EventObject & event)
76 if( !(itk::IterationEvent().CheckEvent( &event )) )
82 m_Optimizer->OutputIterationInfo();
85 // Write intermediate result
86 if (m_ArgsInfo.current_given && (m_IterationCounter> m_ArgsInfo.current_arg) )
89 writeImage<CoefficientImageType>(m_Initializer->GetTransform()->GetCoefficientImage(), m_ArgsInfo.coeff_arg, m_ArgsInfo.verbose_flag );
95 void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
96 OptimizerPointer m_Optimizer;
98 void SetInitializer(InitializerPointer i){m_Initializer=i;}
99 InitializerPointer m_Initializer;
101 void SetArgsInfo(args_info_clitkShapedBLUTSpatioTemporalDIR a){m_ArgsInfo=a;}
102 args_info_clitkShapedBLUTSpatioTemporalDIR m_ArgsInfo;
104 int m_IterationCounter;
108 //==============================================================================
109 //Creating an observer class that allows us to change parameters at subsequent levels
110 //==============================================================================
111 template <typename TRegistration, class args_info_clitkShapedBLUTSpatioTemporalDIR>
112 class RegistrationInterfaceCommand : public itk::Command
115 typedef RegistrationInterfaceCommand Self;
116 typedef itk::Command Superclass;
117 typedef itk::SmartPointer<Self> Pointer;
120 RegistrationInterfaceCommand() {};
124 typedef TRegistration RegistrationType;
125 typedef RegistrationType * RegistrationPointer;
128 typedef typename RegistrationType::FixedImageType FixedImageType;
129 typedef typename FixedImageType::RegionType RegionType;
130 itkStaticConstMacro(ImageDimension, unsigned int,FixedImageType::ImageDimension);
131 typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform<double, ImageDimension, ImageDimension> TransformType;
132 typedef clitk::ShapedBLUTSpatioTemporalDeformableTransformInitializer<TransformType, FixedImageType> InitializerType;
133 typedef typename InitializerType::CoefficientImageType CoefficientImageType;
134 typedef itk::CastImageFilter<CoefficientImageType, CoefficientImageType> CastImageFilterType;
135 typedef typename TransformType::ParametersType ParametersType;
136 typedef typename InitializerType::Pointer InitializerPointer;
137 typedef typename CommandIterationUpdate< TRegistration, args_info_clitkShapedBLUTSpatioTemporalDIR>::Pointer CommandIterationUpdatePointer;
140 typedef clitk::GenericOptimizer<args_info_clitkShapedBLUTSpatioTemporalDIR> GenericOptimizerType;
141 typedef typename GenericOptimizerType::Pointer GenericOptimizerPointer;
144 typedef typename RegistrationType::FixedImageType InternalImageType;
145 typedef clitk::GenericMetric<args_info_clitkShapedBLUTSpatioTemporalDIR, InternalImageType, InternalImageType> GenericMetricType;
146 typedef typename GenericMetricType::Pointer GenericMetricPointer;
148 // Two arguments are passed to the Execute() method: the first
149 // is the pointer to the object which invoked the event and the
150 // second is the event that was invoked.
151 void Execute(itk::Object * object, const itk::EventObject & event)
153 if( !(itk::IterationEvent().CheckEvent( &event )) )
159 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
160 unsigned int numberOfLevels=registration->GetNumberOfLevels();
161 unsigned int currentLevel=registration->GetCurrentLevel()+1;
164 std::cout<<std::endl;
165 std::cout<<"========================================"<<std::endl;
166 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
167 std::cout<<"========================================"<<std::endl;
168 std::cout<<std::endl;
173 // fixed image region pyramid
174 typedef clitk::MultiResolutionPyramidRegionFilter<InternalImageType> FixedImageRegionPyramidType;
175 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
176 fixedImageRegionPyramid->SetRegion(m_MetricRegion);
177 fixedImageRegionPyramid->SetSchedule(registration->GetFixedImagePyramid()->GetSchedule());
179 // Reinitialize the metric (!= number of samples)
180 m_GenericMetric= GenericMetricType::New();
181 m_GenericMetric->SetArgsInfo(m_ArgsInfo);
182 m_GenericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
183 if (m_ArgsInfo.referenceMask_given) m_GenericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
184 m_GenericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(registration->GetCurrentLevel()));
185 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
186 typename MetricType::Pointer metric=m_GenericMetric->GetMetricPointer();
187 registration->SetMetric(metric);
189 // Get the current coefficient image and make a COPY
190 typename itk::ImageDuplicator<CoefficientImageType>::Pointer caster=itk::ImageDuplicator<CoefficientImageType>::New();
191 caster->SetInputImage(m_Initializer->GetTransform()->GetCoefficientImage());
193 typename CoefficientImageType::Pointer currentCoefficientImage=caster->GetOutput();
195 // Write the intermediate result?
196 if (m_ArgsInfo.intermediate_given>=numberOfLevels)
197 writeImage<CoefficientImageType>(currentCoefficientImage, m_ArgsInfo.intermediate_arg[currentLevel-2], m_ArgsInfo.verbose_flag);
199 // Set the new transform properties
200 m_Initializer->SetImage(registration->GetFixedImagePyramid()->GetOutput(currentLevel-1));
201 if( m_Initializer->m_ControlPointSpacingIsGiven)
202 m_Initializer->SetControlPointSpacing(m_Initializer->m_ControlPointSpacingArray[registration->GetCurrentLevel()]);
203 if( m_Initializer->m_NumberOfControlPointsIsGiven)
204 m_Initializer->SetNumberOfControlPointsInsideTheImage(m_Initializer->m_NumberOfControlPointsInsideTheImageArray[registration->GetCurrentLevel()]);
206 // Reinitialize the transform
207 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
208 m_Initializer->InitializeTransform();
209 ParametersType* newParameters= new typename TransformType::ParametersType(m_Initializer->GetTransform()->GetNumberOfParameters());
211 // Reinitialize an Optimizer (!= number of parameters)
212 m_GenericOptimizer = GenericOptimizerType::New();
213 m_GenericOptimizer->SetArgsInfo(m_ArgsInfo);
214 m_GenericOptimizer->SetMaximize(m_Maximize);
215 m_GenericOptimizer->SetNumberOfParameters(m_Initializer->GetTransform()->GetNumberOfParameters());
216 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
217 OptimizerType::Pointer optimizer = m_GenericOptimizer->GetOptimizerPointer();
218 optimizer->AddObserver( itk::IterationEvent(), m_CommandIterationUpdate);
219 registration->SetOptimizer(optimizer);
220 m_CommandIterationUpdate->SetOptimizer(m_GenericOptimizer);
222 // Set the previous transform parameters to the registration
223 // if(m_Initializer->m_Parameters!=NULL )delete m_Initializer->m_Parameters;
224 m_Initializer->SetInitialParameters(currentCoefficientImage,*newParameters);
225 registration->SetInitialTransformParametersOfNextLevel(*newParameters);
229 void Execute(const itk::Object * , const itk::EventObject & )
234 void SetInitializer(InitializerPointer i){m_Initializer=i;}
235 InitializerPointer m_Initializer;
237 void SetArgsInfo(args_info_clitkShapedBLUTSpatioTemporalDIR a){m_ArgsInfo=a;}
238 args_info_clitkShapedBLUTSpatioTemporalDIR m_ArgsInfo;
240 void SetCommandIterationUpdate(CommandIterationUpdatePointer c){m_CommandIterationUpdate=c;};
241 CommandIterationUpdatePointer m_CommandIterationUpdate;
243 GenericOptimizerPointer m_GenericOptimizer;
244 void SetMaximize(bool b){m_Maximize=b;}
247 GenericMetricPointer m_GenericMetric;
248 void SetMetricRegion(RegionType i){m_MetricRegion=i;}
249 RegionType m_MetricRegion;
255 //==============================================================================
256 // Update with the number of dimensions
257 //==============================================================================
258 template<unsigned int Dimension>
260 ShapedBLUTSpatioTemporalDIRGenericFilter::UpdateWithDim(std::string PixelType)
262 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
264 if(PixelType == "short"){
265 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
266 UpdateWithDimAndPixelType<Dimension, signed short>();
268 // else if(PixelType == "unsigned_short"){
269 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
270 // UpdateWithDimAndPixelType<Dimension, unsigned short>();
273 // else if (PixelType == "unsigned_char"){
274 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
275 // UpdateWithDimAndPixelType<Dimension, unsigned char>();
278 // else if (PixelType == "char"){
279 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
280 // UpdateWithDimAndPixelType<Dimension, signed char>();
283 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
284 UpdateWithDimAndPixelType<Dimension, float>();
289 //==============================================================================
290 // Update with the number of dimensions and pixeltype
291 //==============================================================================
292 template<unsigned int ImageDimension, class PixelType>
294 ShapedBLUTSpatioTemporalDIRGenericFilter::UpdateWithDimAndPixelType()
298 //=============================================================================
300 //=============================================================================
301 bool threadsGiven=m_ArgsInfo.threads_given;
302 int threads=m_ArgsInfo.threads_arg;
304 typedef itk::Image< PixelType, ImageDimension > FixedImageType;
305 typedef itk::Image< PixelType, ImageDimension > MovingImageType;
306 const unsigned int SpaceDimension = ImageDimension-1;
307 typedef double TCoordRep;
310 //=======================================================
312 //=======================================================
313 typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
314 typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
316 typename FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
317 typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
319 fixedImageReader->SetFileName( m_ArgsInfo.reference_arg );
320 movingImageReader->SetFileName( m_ArgsInfo.target_arg );
321 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
322 fixedImageReader->Update();
323 movingImageReader->Update();
325 typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
326 typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
327 typename FixedImageType::Pointer croppedFixedImage=fixedImage;
330 //=======================================================
332 //=======================================================
334 // The original input region
335 typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
337 // The transform region with respect to the input region:
338 // where should the transform be DEFINED (depends on mask)
339 typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
340 typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
341 typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
342 typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
344 // The metric region with respect to the extracted transform region:
345 // where should the metric be CALCULATED (depends on transform)
346 typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
347 typename FixedImageType::RegionType::SizeType metricRegionSize=metricRegion.GetSize();
348 typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
349 typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
351 //===========================================================================
352 // If given, we connect a mask to reference or target
353 //============================================================================
354 typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
355 typename MaskType::Pointer fixedMask=NULL;
356 if (m_ArgsInfo.referenceMask_given)
358 fixedMask= MaskType::New();
359 typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
360 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
361 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
362 maskReader->SetFileName(m_ArgsInfo.referenceMask_arg);
365 maskReader->Update();
367 catch( itk::ExceptionObject & err )
369 std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
370 std::cerr << err << std::endl;
373 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
375 // Set the image to the spatialObject
376 fixedMask->SetImage( maskReader->GetOutput() );
378 // Find the bounding box of the "inside" label
379 typedef itk::LabelStatisticsImageFilter<ImageMaskType, ImageMaskType> StatisticsImageFilterType;
380 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
381 statisticsImageFilter->SetInput(maskReader->GetOutput());
382 statisticsImageFilter->SetLabelInput(maskReader->GetOutput());
383 statisticsImageFilter->Update();
384 typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
386 // Limit the transform region to the mask
387 for (unsigned int i=0; i<ImageDimension; i++)
389 transformRegionIndex[i]=boundingBox[2*i];
390 transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
392 transformRegion.SetSize(transformRegionSize);
393 transformRegion.SetIndex(transformRegionIndex);
394 fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
396 // Crop the fixedImage to the bounding box to facilitate multi-resolution
397 typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
398 typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
399 extractImageFilter->SetInput(fixedImage);
400 extractImageFilter->SetExtractionRegion(transformRegion);
401 extractImageFilter->Update();
402 croppedFixedImage=extractImageFilter->GetOutput();
404 // Update the metric region
405 metricRegion = croppedFixedImage->GetLargestPossibleRegion();
406 metricRegionIndex=metricRegion.GetIndex();
407 metricRegionSize=metricRegion.GetSize();
408 croppedFixedImage->TransformIndexToPhysicalPoint(metricRegionIndex, metricRegionOrigin);
410 // Set start index to zero (with respect to croppedFixedImage/transform region)
411 metricRegionIndex.Fill(0);
412 metricRegion.SetIndex(metricRegionIndex);
413 croppedFixedImage->SetRegions(metricRegion);
414 croppedFixedImage->SetOrigin(metricRegionOrigin);
418 typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
419 typename MaskType::Pointer movingMask=NULL;
420 if (m_ArgsInfo.targetMask_given)
422 movingMask= MaskType::New();
423 typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
424 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
425 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
426 maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
429 maskReader->Update();
431 catch( itk::ExceptionObject & err )
433 std::cerr << "ExceptionObject caught !" << std::endl;
434 std::cerr << err << std::endl;
436 if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
438 movingMask->SetImage( maskReader->GetOutput() );
442 //=======================================================
444 //=======================================================
448 // Fixed image region
449 std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
450 <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
451 <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
454 std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
455 <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
456 <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
459 std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
460 <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
461 <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
466 //=======================================================
467 // Pyramids (update them for transform initializer)
468 //=======================================================
469 typedef clitk::RecursiveSpatioTemporalMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
470 typedef clitk::RecursiveSpatioTemporalMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
471 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
472 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
473 fixedImagePyramid->SetUseShrinkImageFilter(false);
474 fixedImagePyramid->SetInput(croppedFixedImage);
475 fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
476 movingImagePyramid->SetUseShrinkImageFilter(false);
477 movingImagePyramid->SetInput(movingImage);
478 movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
479 if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
480 fixedImagePyramid->Update();
481 movingImagePyramid->Update();
482 typedef clitk::MultiResolutionPyramidRegionFilter<FixedImageType> FixedImageRegionPyramidType;
483 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
484 fixedImageRegionPyramid->SetRegion(metricRegion);
485 fixedImageRegionPyramid->SetSchedule(fixedImagePyramid->GetSchedule());
490 // //=======================================================
491 // // Rigid Transform
492 // //=======================================================
493 // typedef itk::Euler3DTransform <double> RigidTransformType;
494 // RigidTransformType::Pointer rigidTransform;
495 // if (m_ArgsInfo.rigid_given)
497 // rigidTransform=RigidTransformType::New();
498 // itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.rigid_arg);
500 // //Set the rotation
501 // itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
502 // rigidTransform->SetMatrix(finalRotation);
504 // //Set the translation
505 // itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
506 // rigidTransform->SetTranslation(finalTranslation);
511 //=======================================================
512 // ShapedBSplineSpatioTemporal Transform
513 //=======================================================
514 typedef clitk::ShapedBLUTSpatioTemporalDeformableTransform<TCoordRep,ImageDimension, ImageDimension > TransformType;
515 typename TransformType::Pointer transform= TransformType::New();
516 transform->SetTransformShape(m_ArgsInfo.shape_arg);
517 if (fixedMask) transform->SetMask( fixedMask );
518 // if (rigidTransform) transform->SetBulkTransform( rigidTransform );
520 //-------------------------------------------------------------------------
521 // The transform initializer
522 //-------------------------------------------------------------------------
523 typedef clitk::ShapedBLUTSpatioTemporalDeformableTransformInitializer< TransformType,FixedImageType> InitializerType;
524 typename InitializerType::Pointer initializer = InitializerType::New();
525 initializer->SetVerbose(m_Verbose);
526 initializer->SetImage(fixedImagePyramid->GetOutput(0));
527 initializer->SetTransform(transform);
529 //-------------------------------------------------------------------------
531 //-------------------------------------------------------------------------
532 typename FixedImageType::RegionType::SizeType splineOrders ;
533 splineOrders.Fill(3);
534 if (m_ArgsInfo.order_given)
535 for(unsigned int i=0; i<ImageDimension;i++)
536 splineOrders[i]=m_ArgsInfo.order_arg[i];
537 if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
538 initializer->SetSplineOrders(splineOrders);
540 //-------------------------------------------------------------------------
542 //-------------------------------------------------------------------------
545 if (m_ArgsInfo.spacing_given)
547 initializer->m_ControlPointSpacingArray.resize(m_ArgsInfo.levels_arg);
548 initializer->SetControlPointSpacing(m_ArgsInfo.spacing_arg);
549 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]=initializer->m_ControlPointSpacing;
550 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]
551 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
553 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
555 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]=initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-i]*2;
556 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i][ImageDimension-1]=2;
557 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]
558 <<" at level "<<m_ArgsInfo.levels_arg-i<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
564 if (m_ArgsInfo.control_given)
566 initializer->m_NumberOfControlPointsInsideTheImageArray.resize(m_ArgsInfo.levels_arg);
567 initializer->SetNumberOfControlPointsInsideTheImage(m_ArgsInfo.control_arg);
568 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]=initializer->m_NumberOfControlPointsInsideTheImage;
569 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]<<"control points inside the image"
570 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
572 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
574 for(unsigned int j=0;j<ImageDimension-1;j++)
575 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i][j]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i][j]/2.);
576 // initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i]/2.);
577 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]<<"control points inside the image"
578 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
583 // Inialize on the first level
584 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level 1 of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
585 if (m_ArgsInfo.spacing_given) initializer->SetControlPointSpacing( initializer->m_ControlPointSpacingArray[0]);
586 if (m_ArgsInfo.control_given) initializer->SetNumberOfControlPointsInsideTheImage(initializer->m_NumberOfControlPointsInsideTheImageArray[0]);
587 if (m_ArgsInfo.samplingFactor_given) initializer->SetSamplingFactors(m_ArgsInfo.samplingFactor_arg);
590 initializer->InitializeTransform();
592 //-------------------------------------------------------------------------
593 // Initial parameters (passed by reference)
594 //-------------------------------------------------------------------------
595 typedef typename TransformType::ParametersType ParametersType;
596 const unsigned int numberOfParameters = transform->GetNumberOfParameters();
597 ParametersType parameters(numberOfParameters);
598 parameters.Fill( 0.0 );
599 transform->SetParameters( parameters );
600 if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters);
603 //=======================================================
605 //=======================================================
606 typedef clitk::GenericInterpolator<args_info_clitkShapedBLUTSpatioTemporalDIR, FixedImageType,TCoordRep > GenericInterpolatorType;
607 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
608 genericInterpolator->SetArgsInfo(m_ArgsInfo);
609 typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
610 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
613 //=======================================================
615 //=======================================================
616 typedef clitk::GenericMetric< args_info_clitkShapedBLUTSpatioTemporalDIR, FixedImageType,MovingImageType > GenericMetricType;
617 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
618 genericMetric->SetArgsInfo(m_ArgsInfo);
619 genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
620 if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
621 genericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(0));
622 typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
623 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
624 if (movingMask) metric->SetMovingImageMask(movingMask);
626 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
627 if (threadsGiven) metric->SetNumberOfThreads( threads );
629 if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
633 //=======================================================
635 //=======================================================
636 typedef clitk::GenericOptimizer<args_info_clitkShapedBLUTSpatioTemporalDIR> GenericOptimizerType;
637 GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
638 genericOptimizer->SetArgsInfo(m_ArgsInfo);
639 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
640 genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters());
641 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
642 OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
645 //=======================================================
647 //=======================================================
648 typedef clitk::SpatioTemporalMultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;
649 typename RegistrationType::Pointer registration = RegistrationType::New();
650 registration->SetMetric( metric );
651 registration->SetOptimizer( optimizer );
652 registration->SetInterpolator( interpolator );
653 registration->SetTransform (transform);
654 if(threadsGiven) registration->SetNumberOfThreads(threads);
655 registration->SetFixedImage( croppedFixedImage );
656 registration->SetMovingImage( movingImage );
657 registration->SetFixedImageRegion( metricRegion );
658 registration->SetFixedImagePyramid( fixedImagePyramid );
659 registration->SetMovingImagePyramid( movingImagePyramid );
660 registration->SetInitialTransformParameters( transform->GetParameters() );
661 registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
662 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
665 //================================================================================================
667 //================================================================================================
670 // Output iteration info
671 typedef CommandIterationUpdate <RegistrationType,args_info_clitkShapedBLUTSpatioTemporalDIR> CommandIterationUpdateType;
672 typename CommandIterationUpdateType::Pointer observer = CommandIterationUpdateType::New();
673 observer->SetInitializer(initializer);
674 observer->SetArgsInfo(m_ArgsInfo);
675 observer->SetOptimizer(genericOptimizer);
676 optimizer->AddObserver( itk::IterationEvent(), observer );
679 typedef RegistrationInterfaceCommand<RegistrationType,args_info_clitkShapedBLUTSpatioTemporalDIR> CommandType;
680 typename CommandType::Pointer command = CommandType::New();
681 command->SetInitializer(initializer);
682 command->SetArgsInfo(m_ArgsInfo);
683 command->SetCommandIterationUpdate(observer);
684 command->SetMaximize(genericMetric->GetMaximize());
685 command->SetMetricRegion(metricRegion);
686 registration->AddObserver( itk::IterationEvent(), command );
690 //=======================================================
692 //=======================================================
693 if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
697 registration->StartRegistration();
699 catch( itk::ExceptionObject & err )
701 std::cerr << "ExceptionObject caught while registering!" << std::endl;
702 std::cerr << err << std::endl;
707 //=======================================================
709 //=======================================================
710 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
711 transform->SetParameters( finalParameters );
714 std::cout<<"Stop condition description: "
715 <<registration->GetOptimizer()->GetStopConditionDescription()<<std::endl;
719 //=======================================================
720 // Get the BSplineSpatioTemporal coefficient images and write them
721 //=======================================================
722 if (m_ArgsInfo.coeff_given)
724 typedef typename TransformType::CoefficientImageType CoefficientImageType;
725 typename CoefficientImageType::Pointer coefficientImage =transform->GetCoefficientImage();
726 typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
727 typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
728 coeffWriter->SetInput(coefficientImage);
729 coeffWriter->SetFileName(m_ArgsInfo.coeff_arg);
730 coeffWriter->Update();
733 if (m_ArgsInfo.padCoeff_given)
735 typedef itk::Image<itk::Vector<TCoordRep, SpaceDimension>, ImageDimension> CoefficientImageType;
736 typename CoefficientImageType::Pointer coefficientImage =transform->GetPaddedCoefficientImage();
737 typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
738 typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
739 coeffWriter->SetInput(coefficientImage);
740 coeffWriter->SetFileName(m_ArgsInfo.padCoeff_arg);
741 coeffWriter->Update();
745 //=======================================================
746 // Generate the DVF (4D with 3 comps and 4D with 4 comps)
747 //=======================================================
748 typedef itk::Vector< float, SpaceDimension > DisplacementType;
749 typedef itk::Vector< float, ImageDimension > Displacement4DType;
750 typedef itk::Image< DisplacementType, ImageDimension > DeformationFieldType;
751 typedef itk::Image< Displacement4DType, ImageDimension > DeformationField4DType;
753 typename DeformationFieldType::Pointer field = DeformationFieldType::New();
754 typename DeformationField4DType::Pointer field4D = DeformationField4DType::New();
755 field->SetRegions( fixedImageRegion );
756 field->SetOrigin( fixedImage->GetOrigin() );
757 field->SetSpacing( fixedImage->GetSpacing() );
758 field->SetDirection( fixedImage->GetDirection() );
760 field4D->SetRegions( fixedImageRegion );
761 field4D->SetOrigin( fixedImage->GetOrigin() );
762 field4D->SetSpacing( fixedImage->GetSpacing() );
763 field4D->SetDirection( fixedImage->GetDirection() );
766 typedef itk::ImageRegionIteratorWithIndex< DeformationFieldType > FieldIterator;
767 typedef itk::ImageRegionIteratorWithIndex< DeformationField4DType > Field4DIterator;
768 FieldIterator fi( field, fixedImageRegion );
769 Field4DIterator fi4D( field4D, fixedImageRegion );
773 typename TransformType::InputPointType fixedPoint;
774 typename TransformType::OutputPointType movingPoint;
775 typename DeformationFieldType::IndexType index;
777 DisplacementType displacement;
778 Displacement4DType displacement4D;
779 displacement4D[ImageDimension-1]=0;
780 while( ! fi.IsAtEnd() )
782 index = fi.GetIndex();
783 field->TransformIndexToPhysicalPoint( index, fixedPoint );
784 movingPoint = transform->TransformPoint( fixedPoint );
785 for (unsigned int i=0; i<SpaceDimension; i++)
787 displacement[i] = movingPoint[i] - fixedPoint[i];
788 displacement4D[i] = movingPoint[i] - fixedPoint[i];
790 fi.Set( displacement );
791 fi4D.Set(displacement4D);
797 //=======================================================
799 //=======================================================
800 typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
801 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
802 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
803 fieldWriter->SetInput( field );
806 fieldWriter->Update();
808 catch( itk::ExceptionObject & excp )
810 std::cerr << "Exception thrown writing the DVF" << std::endl;
811 std::cerr << excp << std::endl;
816 //=======================================================
817 // Resample the moving image
818 //=======================================================
819 typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationField4DType > WarpFilterType;
820 typename WarpFilterType::Pointer warp = WarpFilterType::New();
822 #if ITK_VERSION_MAJOR >= 4
823 warp->SetDisplacementField( field4D );
825 warp->SetDeformationField( field4D );
827 warp->SetInput( movingImageReader->GetOutput() );
828 warp->SetOutputOrigin( fixedImage->GetOrigin() );
829 warp->SetOutputSpacing( fixedImage->GetSpacing() );
830 warp->SetOutputDirection( fixedImage->GetDirection() );
831 warp->SetEdgePaddingValue( 0.0 );
835 //=======================================================
836 // Write the warped image
837 //=======================================================
838 typedef itk::ImageFileWriter< FixedImageType > WriterType;
839 typename WriterType::Pointer writer = WriterType::New();
840 writer->SetFileName( m_ArgsInfo.output_arg );
841 writer->SetInput( warp->GetOutput() );
847 catch( itk::ExceptionObject & err )
849 std::cerr << "ExceptionObject caught !" << std::endl;
850 std::cerr << err << std::endl;
855 //=======================================================
856 // Calculate the difference after the deformable transform
857 //=======================================================
858 typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType;
859 if (m_ArgsInfo.after_given)
861 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
862 difference->SetValidInput( fixedImage );
863 difference->SetTestInput( warp->GetOutput() );
867 difference->Update();
869 catch( itk::ExceptionObject & err )
871 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
872 std::cerr << err << std::endl;
876 typename WriterType::Pointer differenceWriter=WriterType::New();
877 differenceWriter->SetInput(difference->GetOutput());
878 differenceWriter->SetFileName(m_ArgsInfo.after_arg);
879 differenceWriter->Update();
884 //=======================================================
885 // Calculate the difference before the deformable transform
886 //=======================================================
887 if( m_ArgsInfo.before_given )
890 typename FixedImageType::Pointer moving=FixedImageType::New();
891 if (m_ArgsInfo.rigid_given)
893 typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
894 typename ResamplerType::Pointer resampler=ResamplerType::New();
895 resampler->SetInput(movingImage);
896 resampler->SetOutputOrigin(fixedImage->GetOrigin());
897 resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
898 resampler->SetOutputSpacing(fixedImage->GetSpacing());
899 resampler->SetDefaultPixelValue( 0. );
900 //resampler->SetTransform(rigidTransform);
902 moving=resampler->GetOutput();
907 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
908 difference->SetValidInput( fixedImage );
909 difference->SetTestInput( moving );
913 difference->Update();
915 catch( itk::ExceptionObject & err )
917 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
918 std::cerr << err << std::endl;
922 typename WriterType::Pointer differenceWriter=WriterType::New();
923 writer->SetFileName( m_ArgsInfo.before_arg );
924 writer->SetInput( difference->GetOutput() );
933 #endif //#define clitkShapedBLUTSpatioTemporalDIRGenericFilter_txx