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 clitkBLUTDIRGenericFilter_cxx
19 #define clitkBLUTDIRGenericFilter_cxx
21 /* =================================================
22 * @file clitkBLUTDIRGenericFilter.cxx
28 ===================================================*/
30 #include "clitkBLUTDIRGenericFilter.h"
31 #include "clitkBLUTDIRCommandIterationUpdateDVF.h"
32 #include "itkCenteredTransformInitializer.h"
33 #include "itkLabelStatisticsImageFilter.h"
34 #if (ITK_VERSION_MAJOR == 4) && (ITK_VERSION_MINOR < 6)
35 # include "itkTransformToDisplacementFieldSource.h"
37 # include "itkTransformToDisplacementFieldFilter.h"
43 //==============================================================================
44 // Creating an observer class that allows output at each iteration
45 //==============================================================================
46 class CommandIterationUpdate : public itk::Command
49 typedef CommandIterationUpdate Self;
50 typedef itk::Command Superclass;
51 typedef itk::SmartPointer<Self> Pointer;
54 CommandIterationUpdate() {};
56 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> OptimizerType;
57 typedef const OptimizerType * OptimizerPointer;
59 // Set the generic optimizer
60 void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
63 void Execute(itk::Object *caller, const itk::EventObject & event)
65 Execute( (const itk::Object *)caller, event);
68 void Execute(const itk::Object * object, const itk::EventObject & event)
70 if( !(itk::IterationEvent().CheckEvent( &event )) )
75 m_Optimizer->OutputIterationInfo();
78 OptimizerPointer m_Optimizer;
81 //===========================================================================//
83 //==========================================================================//
84 BLUTDIRGenericFilter::BLUTDIRGenericFilter():
85 ImageToImageGenericFilter<Self>("Register DIR")
87 InitializeImageType<2>();
88 InitializeImageType<3>();
92 //=========================================================================//
94 //==========================================================================//
95 void BLUTDIRGenericFilter::SetArgsInfo(const args_info_clitkBLUTDIR & a){
97 if (m_ArgsInfo.reference_given) AddInputFilename(m_ArgsInfo.reference_arg);
99 if (m_ArgsInfo.target_given) {
100 AddInputFilename(m_ArgsInfo.target_arg);
103 if (m_ArgsInfo.output_given) SetOutputFilename(m_ArgsInfo.output_arg);
105 if (m_ArgsInfo.verbose_given) m_Verbose=true;
108 //=========================================================================//
109 //===========================================================================//
110 template<unsigned int Dim>
111 void BLUTDIRGenericFilter::InitializeImageType()
113 ADD_DEFAULT_IMAGE_TYPES(3);
115 //--------------------------------------------------------------------
117 //==============================================================================
118 //Creating an observer class that allows us to change parameters at subsequent levels
119 //==============================================================================
120 template <typename TRegistration,class args_info_clitkBLUTDIR>
121 class RegistrationInterfaceCommand : public itk::Command
124 typedef RegistrationInterfaceCommand Self;
125 typedef itk::Command Superclass;
126 typedef itk::SmartPointer<Self> Pointer;
129 RegistrationInterfaceCommand() { };
133 typedef TRegistration RegistrationType;
134 typedef RegistrationType * RegistrationPointer;
137 typedef typename RegistrationType::FixedImageType FixedImageType;
138 typedef typename FixedImageType::RegionType RegionType;
139 itkStaticConstMacro(ImageDimension, unsigned int,FixedImageType::ImageDimension);
140 typedef clitk::MultipleBSplineDeformableTransform<double, ImageDimension, ImageDimension> TransformType;
141 typedef clitk::MultipleBSplineDeformableTransformInitializer<TransformType, FixedImageType> InitializerType;
142 typedef typename InitializerType::CoefficientImageType CoefficientImageType;
143 typedef itk::CastImageFilter<CoefficientImageType, CoefficientImageType> CastImageFilterType;
144 typedef typename TransformType::ParametersType ParametersType;
145 typedef typename InitializerType::Pointer InitializerPointer;
146 typedef CommandIterationUpdate::Pointer CommandIterationUpdatePointer;
149 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
150 typedef typename GenericOptimizerType::Pointer GenericOptimizerPointer;
153 typedef typename RegistrationType::FixedImageType InternalImageType;
154 typedef clitk::GenericMetric<args_info_clitkBLUTDIR, InternalImageType, InternalImageType> GenericMetricType;
155 typedef typename GenericMetricType::Pointer GenericMetricPointer;
157 // Two arguments are passed to the Execute() method: the first
158 // is the pointer to the object which invoked the event and the
159 // second is the event that was invoked.
160 void Execute(itk::Object * object, const itk::EventObject & event)
162 if( !(itk::IterationEvent().CheckEvent( &event )) )
168 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
169 unsigned int numberOfLevels=registration->GetNumberOfLevels();
170 unsigned int currentLevel=registration->GetCurrentLevel()+1;
173 std::cout<<std::endl;
174 std::cout<<"========================================"<<std::endl;
175 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
176 std::cout<<"========================================"<<std::endl;
177 std::cout<<std::endl;
182 // fixed image region pyramid
183 typedef clitk::MultiResolutionPyramidRegionFilter<InternalImageType> FixedImageRegionPyramidType;
184 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
185 fixedImageRegionPyramid->SetRegion(m_MetricRegion);
186 fixedImageRegionPyramid->SetSchedule(registration->GetFixedImagePyramid()->GetSchedule());
188 // Reinitialize the metric (!= number of samples)
189 m_GenericMetric= GenericMetricType::New();
190 m_GenericMetric->SetArgsInfo(m_ArgsInfo);
191 m_GenericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
192 if (m_ArgsInfo.referenceMask_given) m_GenericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
193 m_GenericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(registration->GetCurrentLevel()));
194 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
195 typename MetricType::Pointer metric=m_GenericMetric->GetMetricPointer();
196 registration->SetMetric(metric);
198 // Get the current coefficient image and make a COPY
199 typename itk::ImageDuplicator<CoefficientImageType>::Pointer caster = itk::ImageDuplicator<CoefficientImageType>::New();
200 std::vector<typename CoefficientImageType::Pointer> currentCoefficientImages = m_Initializer->GetTransform()->GetCoefficientImages();
201 for (unsigned i = 0; i < currentCoefficientImages.size(); ++i)
203 caster->SetInputImage(currentCoefficientImages[i]);
205 currentCoefficientImages[i] = caster->GetOutput();
209 // Write the intermediate result?
210 if (m_ArgsInfo.intermediate_given>=numberOfLevels)
211 writeImage<CoefficientImageType>(currentCoefficientImage, m_ArgsInfo.intermediate_arg[currentLevel-2], m_ArgsInfo.verbose_flag);
214 // Set the new transform properties
215 m_Initializer->SetImage(registration->GetFixedImagePyramid()->GetOutput(currentLevel-1));
216 if( m_Initializer->m_ControlPointSpacingIsGiven)
217 m_Initializer->SetControlPointSpacing(m_Initializer->m_ControlPointSpacingArray[registration->GetCurrentLevel()]);
218 if( m_Initializer->m_NumberOfControlPointsIsGiven)
219 m_Initializer->SetNumberOfControlPointsInsideTheImage(m_Initializer->m_NumberOfControlPointsInsideTheImageArray[registration->GetCurrentLevel()]);
221 // Reinitialize the transform
222 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
223 m_Initializer->InitializeTransform();
224 ParametersType* newParameters= new typename TransformType::ParametersType(m_Initializer->GetTransform()->GetNumberOfParameters());
226 // DS : if we want to skip the last pyramid level, force to only 1 iteration
227 DD(m_ArgsInfo.skipLastPyramidLevel_flag);
228 if ((currentLevel == numberOfLevels) && (m_ArgsInfo.skipLastPyramidLevel_flag)) {
229 DD(m_ArgsInfo.maxIt_arg);
230 std::cout << "I skip the last pyramid level : set max iteration to 0" << std::endl;
231 m_ArgsInfo.maxIt_arg = 0;
232 DD(m_ArgsInfo.maxIt_arg);
235 // Reinitialize an Optimizer (!= number of parameters)
236 m_GenericOptimizer = GenericOptimizerType::New();
237 m_GenericOptimizer->SetArgsInfo(m_ArgsInfo);
238 m_GenericOptimizer->SetMaximize(m_Maximize);
239 m_GenericOptimizer->SetNumberOfParameters(m_Initializer->GetTransform()->GetNumberOfParameters());
242 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
243 OptimizerType::Pointer optimizer = m_GenericOptimizer->GetOptimizerPointer();
244 optimizer->AddObserver( itk::IterationEvent(), m_CommandIterationUpdate);
245 registration->SetOptimizer(optimizer);
246 m_CommandIterationUpdate->SetOptimizer(m_GenericOptimizer);
248 // Set the previous transform parameters to the registration
249 // if(m_Initializer->m_Parameters!=ITK_NULLPTR )delete m_Initializer->m_Parameters;
250 m_Initializer->SetInitialParameters(currentCoefficientImages, *newParameters);
251 registration->SetInitialTransformParametersOfNextLevel(*newParameters);
255 void Execute(const itk::Object * , const itk::EventObject & )
260 void SetInitializer(InitializerPointer i){m_Initializer=i;}
261 InitializerPointer m_Initializer;
263 void SetArgsInfo(args_info_clitkBLUTDIR a){m_ArgsInfo=a;}
264 args_info_clitkBLUTDIR m_ArgsInfo;
266 void SetCommandIterationUpdate(CommandIterationUpdatePointer c){m_CommandIterationUpdate=c;};
267 CommandIterationUpdatePointer m_CommandIterationUpdate;
269 GenericOptimizerPointer m_GenericOptimizer;
270 void SetMaximize(bool b){m_Maximize=b;}
273 GenericMetricPointer m_GenericMetric;
274 void SetMetricRegion(RegionType i){m_MetricRegion=i;}
275 RegionType m_MetricRegion;
280 //==============================================================================
281 // Update with the number of dimensions and pixeltype
282 //==============================================================================
283 template<class InputImageType>
284 void BLUTDIRGenericFilter::UpdateWithInputImageType()
286 if (m_Verbose) std::cout << "BLUTDIRGenericFilter::UpdateWithInputImageType()" << std::endl;
288 //=============================================================================
290 //=============================================================================
291 bool threadsGiven=m_ArgsInfo.threads_given;
292 int threads=m_ArgsInfo.threads_arg;
293 typedef typename InputImageType::PixelType PixelType;
295 typedef double TCoordRep;
297 typename InputImageType::Pointer fixedImage = this->template GetInput<InputImageType>(0);
299 typename InputImageType::Pointer inputFixedImage = this->template GetInput<InputImageType>(0);
302 typename InputImageType::Pointer movingImage = this->template GetInput<InputImageType>(1);
304 typename InputImageType::Pointer inputMovingImage = this->template GetInput<InputImageType>(1);
306 typedef itk::Image< PixelType,InputImageType::ImageDimension > FixedImageType;
307 typedef itk::Image< PixelType, InputImageType::ImageDimension> MovingImageType;
308 const unsigned int SpaceDimension = InputImageType::ImageDimension;
309 //Whatever the pixel type, internally we work with an image represented in float
310 //Reading reference image
311 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
312 //=======================================================
314 //=======================================================
315 typename FixedImageType::Pointer croppedFixedImage=fixedImage;
316 //=======================================================
318 //=======================================================
319 // The original input region
320 typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
322 // The transform region with respect to the input region:
323 // where should the transform be DEFINED (depends on mask)
324 typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
325 typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
326 typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
327 typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
329 // The metric region with respect to the extracted transform region:
330 // where should the metric be CALCULATED (depends on transform)
331 typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
332 typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
333 typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
336 //===========================================================================
337 // If given, we connect a mask to reference or target
338 //============================================================================
339 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
340 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageLabelType;
341 typename MaskType::Pointer fixedMask = ITK_NULLPTR;
342 typename ImageLabelType::Pointer labels = ITK_NULLPTR;
343 if (m_ArgsInfo.referenceMask_given)
345 fixedMask = MaskType::New();
346 labels = ImageLabelType::New();
347 typedef itk::ImageFileReader< ImageLabelType > LabelReaderType;
348 typename LabelReaderType::Pointer labelReader = LabelReaderType::New();
349 labelReader->SetFileName(m_ArgsInfo.referenceMask_arg);
352 labelReader->Update();
354 catch( itk::ExceptionObject & err )
356 std::cerr << "ExceptionObject caught while reading mask or labels !" << std::endl;
357 std::cerr << err << std::endl;
360 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
363 typedef itk::ResampleImageFilter<ImageLabelType, ImageLabelType> ResamplerType;
364 typename ResamplerType::Pointer resampler = ResamplerType::New();
365 typedef itk::NearestNeighborInterpolateImageFunction<ImageLabelType, TCoordRep> InterpolatorType;
366 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
367 resampler->SetOutputParametersFromImage(fixedImage);
368 resampler->SetInterpolator(interpolator);
369 resampler->SetInput(labelReader->GetOutput());
371 labels = resampler->GetOutput();
373 // Set the image to the spatialObject
374 fixedMask->SetImage(labels);
376 // Find the bounding box of the "inside" label
377 typedef itk::LabelStatisticsImageFilter<ImageLabelType, ImageLabelType> StatisticsImageFilterType;
378 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
379 statisticsImageFilter->SetInput(labels);
380 statisticsImageFilter->SetLabelInput(labels);
381 statisticsImageFilter->Update();
382 typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
384 // Limit the transform region to the mask
385 for (unsigned int i=0; i<InputImageType::ImageDimension; i++)
387 transformRegionIndex[i]=boundingBox[2*i];
388 transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
390 transformRegion.SetSize(transformRegionSize);
391 transformRegion.SetIndex(transformRegionIndex);
392 fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
394 // Crop the fixedImage to the bounding box to facilitate multi-resolution
395 typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
396 typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
397 extractImageFilter->SetDirectionCollapseToSubmatrix();
398 extractImageFilter->SetInput(fixedImage);
399 extractImageFilter->SetExtractionRegion(transformRegion);
400 extractImageFilter->Update();
401 croppedFixedImage=extractImageFilter->GetOutput();
403 // Update the metric region
404 metricRegion = croppedFixedImage->GetLargestPossibleRegion();
405 metricRegionIndex=metricRegion.GetIndex();
406 croppedFixedImage->TransformIndexToPhysicalPoint(metricRegionIndex, metricRegionOrigin);
408 // Set start index to zero (with respect to croppedFixedImage/transform region)
409 metricRegionIndex.Fill(0);
410 metricRegion.SetIndex(metricRegionIndex);
411 croppedFixedImage->SetRegions(metricRegion);
412 croppedFixedImage->SetOrigin(metricRegionOrigin);
415 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
416 typename MaskType::Pointer movingMask=ITK_NULLPTR;
417 if (m_ArgsInfo.targetMask_given)
419 movingMask= MaskType::New();
420 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
421 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
422 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
423 maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
426 maskReader->Update();
428 catch( itk::ExceptionObject & err )
430 std::cerr << "ExceptionObject caught !" << std::endl;
431 std::cerr << err << std::endl;
433 if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
435 movingMask->SetImage( maskReader->GetOutput() );
439 //=======================================================
441 //=======================================================
445 // Fixed image region
446 std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
447 <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
448 <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
451 std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
452 <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
453 <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
456 std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
457 <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
458 <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
463 //=======================================================
464 // Pyramids (update them for transform initializer)
465 //=======================================================
466 typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
467 typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
468 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
469 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
470 fixedImagePyramid->SetUseShrinkImageFilter(false);
471 fixedImagePyramid->SetInput(croppedFixedImage);
472 fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
473 movingImagePyramid->SetUseShrinkImageFilter(false);
474 movingImagePyramid->SetInput(movingImage);
475 movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
476 if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
477 fixedImagePyramid->Update();
478 movingImagePyramid->Update();
479 typedef clitk::MultiResolutionPyramidRegionFilter<FixedImageType> FixedImageRegionPyramidType;
480 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
481 fixedImageRegionPyramid->SetRegion(metricRegion);
482 fixedImageRegionPyramid->SetSchedule(fixedImagePyramid->GetSchedule());
485 //=======================================================
486 // Rigid or Affine Transform
487 //=======================================================
488 typedef itk::AffineTransform <double,3> RigidTransformType;
489 RigidTransformType::Pointer rigidTransform=ITK_NULLPTR;
490 if (m_ArgsInfo.initMatrix_given)
492 if(m_Verbose) std::cout<<"Reading the prior transform matrix "<< m_ArgsInfo.initMatrix_arg<<"..."<<std::endl;
493 rigidTransform=RigidTransformType::New();
494 itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.initMatrix_arg);
497 itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
498 rigidTransform->SetMatrix(finalRotation);
500 //Set the translation
501 itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
502 rigidTransform->SetTranslation(finalTranslation);
504 else if (m_ArgsInfo.centre_flag)
506 if(m_Verbose) std::cout<<"No itinial matrix given and \"centre\" flag switched on. Centering all images..."<<std::endl;
508 rigidTransform=RigidTransformType::New();
510 typedef itk::CenteredTransformInitializer<RigidTransformType, FixedImageType, MovingImageType > TransformInitializerType;
511 typename TransformInitializerType::Pointer initializer = TransformInitializerType::New();
512 initializer->SetTransform( rigidTransform );
513 initializer->SetFixedImage( fixedImage );
514 initializer->SetMovingImage( movingImage );
515 initializer->GeometryOn();
516 initializer->InitializeTransform();
520 //=======================================================
521 // B-LUT FFD Transform
522 //=======================================================
523 typedef clitk::MultipleBSplineDeformableTransform<TCoordRep,InputImageType::ImageDimension, SpaceDimension > TransformType;
524 typename TransformType::Pointer transform = TransformType::New();
525 if (labels) transform->SetLabels(labels);
526 if (rigidTransform) transform->SetBulkTransform(rigidTransform);
528 //-------------------------------------------------------------------------
529 // The transform initializer
530 //-------------------------------------------------------------------------
531 typedef clitk::MultipleBSplineDeformableTransformInitializer< TransformType,FixedImageType> InitializerType;
532 typename InitializerType::Pointer initializer = InitializerType::New();
533 initializer->SetVerbose(m_Verbose);
534 initializer->SetImage(fixedImagePyramid->GetOutput(0));
535 initializer->SetTransform(transform);
537 //-------------------------------------------------------------------------
539 //-------------------------------------------------------------------------
540 typename FixedImageType::RegionType::SizeType splineOrders ;
541 splineOrders.Fill(3);
542 if (m_ArgsInfo.order_given)
543 for(unsigned int i=0; i<InputImageType::ImageDimension;i++)
544 splineOrders[i]=m_ArgsInfo.order_arg[i];
545 if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
546 initializer->SetSplineOrders(splineOrders);
548 //-------------------------------------------------------------------------
550 //-------------------------------------------------------------------------
553 if (m_ArgsInfo.spacing_given)
555 initializer->m_ControlPointSpacingArray.resize(m_ArgsInfo.levels_arg);
556 initializer->SetControlPointSpacing(m_ArgsInfo.spacing_arg);
557 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]=initializer->m_ControlPointSpacing;
558 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]
559 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
561 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
563 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]=initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-i]*2;
564 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]
565 <<" at level "<<m_ArgsInfo.levels_arg-i<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
571 if (m_ArgsInfo.control_given)
573 initializer->m_NumberOfControlPointsInsideTheImageArray.resize(m_ArgsInfo.levels_arg);
574 initializer->SetNumberOfControlPointsInsideTheImage(m_ArgsInfo.control_arg);
575 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]=initializer->m_NumberOfControlPointsInsideTheImage;
576 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]<<"control points inside the image"
577 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
579 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
581 for(unsigned int j=0;j<InputImageType::ImageDimension;j++)
582 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i][j]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i][j]/2.);
583 // initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i]/2.);
584 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]<<"control points inside the image"
585 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
590 // Inialize on the first level
591 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level 1 of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
592 if (m_ArgsInfo.spacing_given) initializer->SetControlPointSpacing( initializer->m_ControlPointSpacingArray[0]);
593 if (m_ArgsInfo.control_given) initializer->SetNumberOfControlPointsInsideTheImage(initializer->m_NumberOfControlPointsInsideTheImageArray[0]);
594 if (m_ArgsInfo.samplingFactor_given) initializer->SetSamplingFactors(m_ArgsInfo.samplingFactor_arg);
597 initializer->InitializeTransform();
599 //-------------------------------------------------------------------------
600 // Initial parameters (passed by reference)
601 //-------------------------------------------------------------------------
602 typedef typename TransformType::ParametersType ParametersType;
603 const unsigned int numberOfParameters = transform->GetNumberOfParameters();
604 ParametersType parameters(numberOfParameters);
605 parameters.Fill( 0.0 );
606 transform->SetParameters( parameters );
607 if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters);
609 //-------------------------------------------------------------------------
610 // DEBUG: use an itk BSpline instead of multilabel BLUTs
611 //-------------------------------------------------------------------------
612 typedef itk::Transform< TCoordRep, 3, 3 > RegistrationTransformType;
613 RegistrationTransformType::Pointer regTransform(transform);
614 typedef itk::BSplineDeformableTransform<TCoordRep,SpaceDimension, 3> SingleBSplineTransformType;
615 typename SingleBSplineTransformType::Pointer sTransform;
616 if(m_ArgsInfo.itkbspline_flag) {
617 if( transform->GetTransforms().size()>1)
618 itkExceptionMacro(<< "invalid --itkbspline option if there is more than 1 label")
619 sTransform = SingleBSplineTransformType::New();
620 sTransform->SetBulkTransform( transform->GetTransforms()[0]->GetBulkTransform() );
621 sTransform->SetGridSpacing( transform->GetTransforms()[0]->GetGridSpacing() );
622 sTransform->SetGridOrigin( transform->GetTransforms()[0]->GetGridOrigin() );
623 sTransform->SetGridRegion( transform->GetTransforms()[0]->GetGridRegion() );
624 sTransform->SetParameters( transform->GetTransforms()[0]->GetParameters() );
625 regTransform = sTransform;
626 transform = ITK_NULLPTR; // free memory
629 //=======================================================
631 //=======================================================
632 typedef clitk::GenericInterpolator<args_info_clitkBLUTDIR, FixedImageType, TCoordRep > GenericInterpolatorType;
633 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
634 genericInterpolator->SetArgsInfo(m_ArgsInfo);
635 typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
636 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
639 //=======================================================
641 //=======================================================
642 typedef clitk::GenericMetric<args_info_clitkBLUTDIR, FixedImageType,MovingImageType > GenericMetricType;
643 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
644 genericMetric->SetArgsInfo(m_ArgsInfo);
645 genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
646 if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
647 genericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(0));
648 typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
649 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
650 if (movingMask) metric->SetMovingImageMask(movingMask);
652 #if ITK_VERSION_MAJOR <= 4
653 metric->SetNumberOfThreads( threads );
655 metric->SetNumberOfWorkUnits( threads );
657 if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
660 //=======================================================
662 //=======================================================
663 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
664 GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
665 genericOptimizer->SetArgsInfo(m_ArgsInfo);
666 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
667 genericOptimizer->SetNumberOfParameters(regTransform->GetNumberOfParameters());
668 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
669 OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
672 //=======================================================
674 //=======================================================
675 typedef itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;
676 typename RegistrationType::Pointer registration = RegistrationType::New();
677 registration->SetMetric( metric );
678 registration->SetOptimizer( optimizer );
679 registration->SetInterpolator( interpolator );
680 registration->SetTransform (regTransform );
682 #if ITK_VERSION_MAJOR <= 4
683 registration->SetNumberOfThreads(threads);
685 registration->SetNumberOfWorkUnits(threads);
687 if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
689 registration->SetFixedImage( croppedFixedImage );
690 registration->SetMovingImage( movingImage );
691 registration->SetFixedImageRegion( metricRegion );
692 registration->SetFixedImagePyramid( fixedImagePyramid );
693 registration->SetMovingImagePyramid( movingImagePyramid );
694 registration->SetInitialTransformParameters( regTransform->GetParameters() );
695 registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
696 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
699 //================================================================================================
701 //================================================================================================
704 // Output iteration info
705 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
706 observer->SetOptimizer(genericOptimizer);
707 optimizer->AddObserver( itk::IterationEvent(), observer );
710 typedef RegistrationInterfaceCommand<RegistrationType,args_info_clitkBLUTDIR> CommandType;
711 typename CommandType::Pointer command = CommandType::New();
712 command->SetInitializer(initializer);
713 command->SetArgsInfo(m_ArgsInfo);
714 command->SetCommandIterationUpdate(observer);
715 command->SetMaximize(genericMetric->GetMaximize());
716 command->SetMetricRegion(metricRegion);
717 registration->AddObserver( itk::IterationEvent(), command );
719 if (m_ArgsInfo.coeff_given)
721 if(m_ArgsInfo.itkbspline_flag) {
722 itkExceptionMacro("--coeff and --itkbpline are incompatible");
725 std::cout << std::endl << "Output coefficient images every " << m_ArgsInfo.coeffEveryN_arg << " iterations." << std::endl;
726 typedef CommandIterationUpdateDVF<FixedImageType, OptimizerType, TransformType> DVFCommandType;
727 typename DVFCommandType::Pointer observerdvf = DVFCommandType::New();
728 observerdvf->SetFixedImage(fixedImage);
729 observerdvf->SetTransform(transform);
730 observerdvf->SetArgsInfo(m_ArgsInfo);
731 optimizer->AddObserver( itk::IterationEvent(), observerdvf );
736 //=======================================================
738 //=======================================================
739 if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
743 registration->Update();
745 catch( itk::ExceptionObject & err )
747 std::cerr << "ExceptionObject caught while registering!" << std::endl;
748 std::cerr << err << std::endl;
753 //=======================================================
755 //=======================================================
756 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
757 regTransform->SetParameters( finalParameters );
760 std::cout<<"Stop condition description: "
761 <<registration->GetOptimizer()->GetStopConditionDescription()<<std::endl;
765 //=======================================================
766 // Get the BSpline coefficient images and write them
767 //=======================================================
768 if (m_ArgsInfo.coeff_given)
770 typedef typename TransformType::CoefficientImageType CoefficientImageType;
771 std::vector<typename CoefficientImageType::Pointer> coefficientImages = transform->GetCoefficientImages();
772 typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
773 typename CoeffWriterType::Pointer coeffWriter = CoeffWriterType::New();
774 unsigned nLabels = transform->GetnLabels();
776 std::string fname(m_ArgsInfo.coeff_arg);
777 int dotpos = fname.length() - 1;
778 while (dotpos >= 0 && fname[dotpos] != '.')
781 for (unsigned i = 0; i < nLabels; ++i)
783 std::ostringstream osfname;
784 osfname << fname.substr(0, dotpos) << '_' << i << fname.substr(dotpos);
785 coeffWriter->SetInput(coefficientImages[i]);
786 coeffWriter->SetFileName(osfname.str());
787 coeffWriter->Update();
793 //=======================================================
794 // Compute the DVF (only deformable transform)
795 //=======================================================
796 typedef itk::Vector< float, SpaceDimension > DisplacementType;
797 typedef itk::Image< DisplacementType, InputImageType::ImageDimension > DisplacementFieldType;
798 #if (ITK_VERSION_MAJOR == 4) && (ITK_VERSION_MINOR < 6)
799 typedef itk::TransformToDisplacementFieldSource<DisplacementFieldType, double> ConvertorType;
801 typedef itk::TransformToDisplacementFieldFilter<DisplacementFieldType, double> ConvertorType;
803 typename ConvertorType::Pointer filter= ConvertorType::New();
804 #if ITK_VERSION_MAJOR <= 4
805 filter->SetNumberOfThreads(1);
807 filter->SetNumberOfWorkUnits(1);
809 if(m_ArgsInfo.itkbspline_flag)
810 sTransform->SetBulkTransform(ITK_NULLPTR);
812 transform->SetBulkTransform(ITK_NULLPTR);
813 filter->SetTransform(regTransform);
814 #if ITK_VERSION_MAJOR > 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR >= 6)
815 filter->SetReferenceImage(fixedImage);
817 filter->SetOutputParametersFromImage(fixedImage);
820 typename DisplacementFieldType::Pointer field = filter->GetOutput();
823 //=======================================================
825 //=======================================================
826 typedef itk::ImageFileWriter< DisplacementFieldType > FieldWriterType;
827 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
828 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
829 fieldWriter->SetInput( field );
832 fieldWriter->Update();
834 catch( itk::ExceptionObject & excp )
836 std::cerr << "Exception thrown writing the DVF" << std::endl;
837 std::cerr << excp << std::endl;
842 //=======================================================
843 // Resample the moving image
844 //=======================================================
845 typedef itk::ResampleImageFilter< MovingImageType, FixedImageType > ResampleFilterType;
846 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
847 if (rigidTransform) {
848 if(m_ArgsInfo.itkbspline_flag)
849 sTransform->SetBulkTransform(rigidTransform);
851 transform->SetBulkTransform(rigidTransform);
853 resampler->SetTransform( regTransform );
854 resampler->SetInput( movingImage);
855 resampler->SetOutputParametersFromImage(fixedImage);
857 typename FixedImageType::Pointer result=resampler->GetOutput();
859 // typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType;
860 // typename WarpFilterType::Pointer warp = WarpFilterType::New();
862 // warp->SetDeformationField( field );
863 // warp->SetInput( movingImageReader->GetOutput() );
864 // warp->SetOutputOrigin( fixedImage->GetOrigin() );
865 // warp->SetOutputSpacing( fixedImage->GetSpacing() );
866 // warp->SetOutputDirection( fixedImage->GetDirection() );
867 // warp->SetEdgePaddingValue( 0.0 );
871 //=======================================================
872 // Write the warped image
873 //=======================================================
874 typedef itk::ImageFileWriter< FixedImageType > WriterType;
875 typename WriterType::Pointer writer = WriterType::New();
876 writer->SetFileName( m_ArgsInfo.output_arg );
877 writer->SetInput( result );
883 catch( itk::ExceptionObject & err )
885 std::cerr << "ExceptionObject caught !" << std::endl;
886 std::cerr << err << std::endl;
891 //=======================================================
892 // Calculate the difference after the deformable transform
893 //=======================================================
894 typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType;
895 if (m_ArgsInfo.after_given)
897 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
898 difference->SetValidInput( fixedImage );
899 difference->SetTestInput( result );
903 difference->Update();
905 catch( itk::ExceptionObject & err )
907 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
908 std::cerr << err << std::endl;
912 typename WriterType::Pointer differenceWriter=WriterType::New();
913 differenceWriter->SetInput(difference->GetOutput());
914 differenceWriter->SetFileName(m_ArgsInfo.after_arg);
915 differenceWriter->Update();
920 //=======================================================
921 // Calculate the difference before the deformable transform
922 //=======================================================
923 if( m_ArgsInfo.before_given )
926 typename FixedImageType::Pointer moving=FixedImageType::New();
927 if (m_ArgsInfo.initMatrix_given)
929 typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
930 typename ResamplerType::Pointer resampler=ResamplerType::New();
931 resampler->SetInput(movingImage);
932 resampler->SetOutputOrigin(fixedImage->GetOrigin());
933 resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
934 resampler->SetOutputSpacing(fixedImage->GetSpacing());
935 resampler->SetDefaultPixelValue( 0. );
936 if (rigidTransform ) resampler->SetTransform(rigidTransform);
938 moving=resampler->GetOutput();
943 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
944 difference->SetValidInput( fixedImage );
945 difference->SetTestInput( moving );
949 difference->Update();
951 catch( itk::ExceptionObject & err )
953 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
954 std::cerr << err << std::endl;
958 typename WriterType::Pointer differenceWriter=WriterType::New();
959 writer->SetFileName( m_ArgsInfo.before_arg );
960 writer->SetInput( difference->GetOutput() );
969 #endif // #define clitkBLUTDIRGenericFilter_txx