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"
37 //==============================================================================
38 // Creating an observer class that allows output at each iteration
39 //==============================================================================
40 class CommandIterationUpdate : public itk::Command
43 typedef CommandIterationUpdate Self;
44 typedef itk::Command Superclass;
45 typedef itk::SmartPointer<Self> Pointer;
48 CommandIterationUpdate() {};
50 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> OptimizerType;
51 typedef const OptimizerType * OptimizerPointer;
53 // Set the generic optimizer
54 void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
57 void Execute(itk::Object *caller, const itk::EventObject & event)
59 Execute( (const itk::Object *)caller, event);
62 void Execute(const itk::Object * object, const itk::EventObject & event)
64 if( !(itk::IterationEvent().CheckEvent( &event )) )
69 m_Optimizer->OutputIterationInfo();
72 OptimizerPointer m_Optimizer;
75 //===========================================================================//
77 //==========================================================================//
78 BLUTDIRGenericFilter::BLUTDIRGenericFilter():
79 ImageToImageGenericFilter<Self>("Register DIR")
81 InitializeImageType<2>();
82 InitializeImageType<3>();
86 //=========================================================================//
88 //==========================================================================//
89 void BLUTDIRGenericFilter::SetArgsInfo(const args_info_clitkBLUTDIR & a){
91 if (m_ArgsInfo.reference_given) AddInputFilename(m_ArgsInfo.reference_arg);
93 if (m_ArgsInfo.target_given) {
94 AddInputFilename(m_ArgsInfo.target_arg);
97 if (m_ArgsInfo.output_given) SetOutputFilename(m_ArgsInfo.output_arg);
99 if (m_ArgsInfo.verbose_given) m_Verbose=true;
102 //=========================================================================//
103 //===========================================================================//
104 template<unsigned int Dim>
105 void BLUTDIRGenericFilter::InitializeImageType()
107 ADD_DEFAULT_IMAGE_TYPES(3);
109 //--------------------------------------------------------------------
111 //==============================================================================
112 //Creating an observer class that allows us to change parameters at subsequent levels
113 //==============================================================================
114 template <typename TRegistration,class args_info_clitkBLUTDIR>
115 class RegistrationInterfaceCommand : public itk::Command
118 typedef RegistrationInterfaceCommand Self;
119 typedef itk::Command Superclass;
120 typedef itk::SmartPointer<Self> Pointer;
123 RegistrationInterfaceCommand() { };
127 typedef TRegistration RegistrationType;
128 typedef RegistrationType * RegistrationPointer;
131 typedef typename RegistrationType::FixedImageType FixedImageType;
132 typedef typename FixedImageType::RegionType RegionType;
133 itkStaticConstMacro(ImageDimension, unsigned int,FixedImageType::ImageDimension);
134 typedef clitk::MultipleBSplineDeformableTransform<double, ImageDimension, ImageDimension> TransformType;
135 typedef clitk::MultipleBSplineDeformableTransformInitializer<TransformType, FixedImageType> InitializerType;
136 typedef typename InitializerType::CoefficientImageType CoefficientImageType;
137 typedef itk::CastImageFilter<CoefficientImageType, CoefficientImageType> CastImageFilterType;
138 typedef typename TransformType::ParametersType ParametersType;
139 typedef typename InitializerType::Pointer InitializerPointer;
140 typedef CommandIterationUpdate::Pointer CommandIterationUpdatePointer;
143 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
144 typedef typename GenericOptimizerType::Pointer GenericOptimizerPointer;
147 typedef typename RegistrationType::FixedImageType InternalImageType;
148 typedef clitk::GenericMetric<args_info_clitkBLUTDIR, InternalImageType, InternalImageType> GenericMetricType;
149 typedef typename GenericMetricType::Pointer GenericMetricPointer;
151 // Two arguments are passed to the Execute() method: the first
152 // is the pointer to the object which invoked the event and the
153 // second is the event that was invoked.
154 void Execute(itk::Object * object, const itk::EventObject & event)
156 if( !(itk::IterationEvent().CheckEvent( &event )) )
162 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
163 unsigned int numberOfLevels=registration->GetNumberOfLevels();
164 unsigned int currentLevel=registration->GetCurrentLevel()+1;
167 std::cout<<std::endl;
168 std::cout<<"========================================"<<std::endl;
169 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
170 std::cout<<"========================================"<<std::endl;
171 std::cout<<std::endl;
176 // fixed image region pyramid
177 typedef clitk::MultiResolutionPyramidRegionFilter<InternalImageType> FixedImageRegionPyramidType;
178 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
179 fixedImageRegionPyramid->SetRegion(m_MetricRegion);
180 fixedImageRegionPyramid->SetSchedule(registration->GetFixedImagePyramid()->GetSchedule());
182 // Reinitialize the metric (!= number of samples)
183 m_GenericMetric= GenericMetricType::New();
184 m_GenericMetric->SetArgsInfo(m_ArgsInfo);
185 m_GenericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
186 if (m_ArgsInfo.referenceMask_given) m_GenericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
187 m_GenericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(registration->GetCurrentLevel()));
188 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
189 typename MetricType::Pointer metric=m_GenericMetric->GetMetricPointer();
190 registration->SetMetric(metric);
192 // Get the current coefficient image and make a COPY
193 typename itk::ImageDuplicator<CoefficientImageType>::Pointer caster = itk::ImageDuplicator<CoefficientImageType>::New();
194 std::vector<typename CoefficientImageType::Pointer> currentCoefficientImages = m_Initializer->GetTransform()->GetCoefficientImages();
195 for (unsigned i = 0; i < currentCoefficientImages.size(); ++i)
197 caster->SetInputImage(currentCoefficientImages[i]);
199 currentCoefficientImages[i] = caster->GetOutput();
203 // Write the intermediate result?
204 if (m_ArgsInfo.intermediate_given>=numberOfLevels)
205 writeImage<CoefficientImageType>(currentCoefficientImage, m_ArgsInfo.intermediate_arg[currentLevel-2], m_ArgsInfo.verbose_flag);
208 // Set the new transform properties
209 m_Initializer->SetImage(registration->GetFixedImagePyramid()->GetOutput(currentLevel-1));
210 if( m_Initializer->m_ControlPointSpacingIsGiven)
211 m_Initializer->SetControlPointSpacing(m_Initializer->m_ControlPointSpacingArray[registration->GetCurrentLevel()]);
212 if( m_Initializer->m_NumberOfControlPointsIsGiven)
213 m_Initializer->SetNumberOfControlPointsInsideTheImage(m_Initializer->m_NumberOfControlPointsInsideTheImageArray[registration->GetCurrentLevel()]);
215 // Reinitialize the transform
216 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
217 m_Initializer->InitializeTransform();
218 ParametersType* newParameters= new typename TransformType::ParametersType(m_Initializer->GetTransform()->GetNumberOfParameters());
220 // DS : if we want to skip the last pyramid level, force to only 1 iteration
221 DD(m_ArgsInfo.skipLastPyramidLevel_flag);
222 if ((currentLevel == numberOfLevels) && (m_ArgsInfo.skipLastPyramidLevel_flag)) {
223 DD(m_ArgsInfo.maxIt_arg);
224 std::cout << "I skip the last pyramid level : set max iteration to 0" << std::endl;
225 m_ArgsInfo.maxIt_arg = 0;
226 DD(m_ArgsInfo.maxIt_arg);
229 // Reinitialize an Optimizer (!= number of parameters)
230 m_GenericOptimizer = GenericOptimizerType::New();
231 m_GenericOptimizer->SetArgsInfo(m_ArgsInfo);
232 m_GenericOptimizer->SetMaximize(m_Maximize);
233 m_GenericOptimizer->SetNumberOfParameters(m_Initializer->GetTransform()->GetNumberOfParameters());
236 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
237 OptimizerType::Pointer optimizer = m_GenericOptimizer->GetOptimizerPointer();
238 optimizer->AddObserver( itk::IterationEvent(), m_CommandIterationUpdate);
239 registration->SetOptimizer(optimizer);
240 m_CommandIterationUpdate->SetOptimizer(m_GenericOptimizer);
242 // Set the previous transform parameters to the registration
243 // if(m_Initializer->m_Parameters!=NULL )delete m_Initializer->m_Parameters;
244 m_Initializer->SetInitialParameters(currentCoefficientImages, *newParameters);
245 registration->SetInitialTransformParametersOfNextLevel(*newParameters);
249 void Execute(const itk::Object * , const itk::EventObject & )
254 void SetInitializer(InitializerPointer i){m_Initializer=i;}
255 InitializerPointer m_Initializer;
257 void SetArgsInfo(args_info_clitkBLUTDIR a){m_ArgsInfo=a;}
258 args_info_clitkBLUTDIR m_ArgsInfo;
260 void SetCommandIterationUpdate(CommandIterationUpdatePointer c){m_CommandIterationUpdate=c;};
261 CommandIterationUpdatePointer m_CommandIterationUpdate;
263 GenericOptimizerPointer m_GenericOptimizer;
264 void SetMaximize(bool b){m_Maximize=b;}
267 GenericMetricPointer m_GenericMetric;
268 void SetMetricRegion(RegionType i){m_MetricRegion=i;}
269 RegionType m_MetricRegion;
274 //==============================================================================
275 // Update with the number of dimensions and pixeltype
276 //==============================================================================
277 template<class InputImageType>
278 void BLUTDIRGenericFilter::UpdateWithInputImageType()
280 if (m_Verbose) std::cout << "BLUTDIRGenericFilter::UpdateWithInputImageType()" << std::endl;
282 //=============================================================================
284 //=============================================================================
285 bool threadsGiven=m_ArgsInfo.threads_given;
286 int threads=m_ArgsInfo.threads_arg;
287 typedef typename InputImageType::PixelType PixelType;
289 typedef double TCoordRep;
291 typename InputImageType::Pointer fixedImage = this->template GetInput<InputImageType>(0);
293 typename InputImageType::Pointer inputFixedImage = this->template GetInput<InputImageType>(0);
296 typename InputImageType::Pointer movingImage = this->template GetInput<InputImageType>(1);
298 typename InputImageType::Pointer inputMovingImage = this->template GetInput<InputImageType>(1);
300 typedef itk::Image< PixelType,InputImageType::ImageDimension > FixedImageType;
301 typedef itk::Image< PixelType, InputImageType::ImageDimension> MovingImageType;
302 const unsigned int SpaceDimension = InputImageType::ImageDimension;
303 //Whatever the pixel type, internally we work with an image represented in float
304 //Reading reference image
305 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
306 //=======================================================
308 //=======================================================
309 typename FixedImageType::Pointer croppedFixedImage=fixedImage;
310 //=======================================================
312 //=======================================================
313 // The original input region
314 typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
316 // The transform region with respect to the input region:
317 // where should the transform be DEFINED (depends on mask)
318 typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
319 typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
320 typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
321 typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
323 // The metric region with respect to the extracted transform region:
324 // where should the metric be CALCULATED (depends on transform)
325 typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
326 typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
327 typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
330 //===========================================================================
331 // If given, we connect a mask to reference or target
332 //============================================================================
333 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
334 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageLabelType;
335 typename MaskType::Pointer fixedMask = NULL;
336 typename ImageLabelType::Pointer labels = NULL;
337 if (m_ArgsInfo.referenceMask_given)
339 fixedMask = MaskType::New();
340 labels = ImageLabelType::New();
341 typedef itk::ImageFileReader< ImageLabelType > LabelReaderType;
342 typename LabelReaderType::Pointer labelReader = LabelReaderType::New();
343 labelReader->SetFileName(m_ArgsInfo.referenceMask_arg);
346 labelReader->Update();
348 catch( itk::ExceptionObject & err )
350 std::cerr << "ExceptionObject caught while reading mask or labels !" << std::endl;
351 std::cerr << err << std::endl;
354 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
357 typedef itk::ResampleImageFilter<ImageLabelType, ImageLabelType> ResamplerType;
358 typename ResamplerType::Pointer resampler = ResamplerType::New();
359 typedef itk::NearestNeighborInterpolateImageFunction<ImageLabelType, TCoordRep> InterpolatorType;
360 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
361 resampler->SetOutputParametersFromImage(fixedImage);
362 resampler->SetInterpolator(interpolator);
363 resampler->SetInput(labelReader->GetOutput());
365 labels = resampler->GetOutput();
367 // Set the image to the spatialObject
368 fixedMask->SetImage(labels);
370 // Find the bounding box of the "inside" label
371 typedef itk::LabelGeometryImageFilter<ImageLabelType> GeometryImageFilterType;
372 typename GeometryImageFilterType::Pointer geometryImageFilter=GeometryImageFilterType::New();
373 geometryImageFilter->SetInput(labels);
374 geometryImageFilter->Update();
375 typename GeometryImageFilterType::BoundingBoxType boundingBox = geometryImageFilter->GetBoundingBox(1);
377 // Limit the transform region to the mask
378 for (unsigned int i=0; i<InputImageType::ImageDimension; i++)
380 transformRegionIndex[i]=boundingBox[2*i];
381 transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
383 transformRegion.SetSize(transformRegionSize);
384 transformRegion.SetIndex(transformRegionIndex);
385 fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
387 // Crop the fixedImage to the bounding box to facilitate multi-resolution
388 typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
389 typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
390 #if ITK_VERSION_MAJOR == 4
391 extractImageFilter->SetDirectionCollapseToSubmatrix();
393 extractImageFilter->SetInput(fixedImage);
394 extractImageFilter->SetExtractionRegion(transformRegion);
395 extractImageFilter->Update();
396 croppedFixedImage=extractImageFilter->GetOutput();
398 // Update the metric region
399 metricRegion = croppedFixedImage->GetLargestPossibleRegion();
400 metricRegionIndex=metricRegion.GetIndex();
401 croppedFixedImage->TransformIndexToPhysicalPoint(metricRegionIndex, metricRegionOrigin);
403 // Set start index to zero (with respect to croppedFixedImage/transform region)
404 metricRegionIndex.Fill(0);
405 metricRegion.SetIndex(metricRegionIndex);
406 croppedFixedImage->SetRegions(metricRegion);
407 croppedFixedImage->SetOrigin(metricRegionOrigin);
410 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
411 typename MaskType::Pointer movingMask=NULL;
412 if (m_ArgsInfo.targetMask_given)
414 movingMask= MaskType::New();
415 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
416 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
417 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
418 maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
421 maskReader->Update();
423 catch( itk::ExceptionObject & err )
425 std::cerr << "ExceptionObject caught !" << std::endl;
426 std::cerr << err << std::endl;
428 if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
430 movingMask->SetImage( maskReader->GetOutput() );
434 //=======================================================
436 //=======================================================
440 // Fixed image region
441 std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
442 <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
443 <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
446 std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
447 <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
448 <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
451 std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
452 <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
453 <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
458 //=======================================================
459 // Pyramids (update them for transform initializer)
460 //=======================================================
461 typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
462 typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
463 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
464 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
465 fixedImagePyramid->SetUseShrinkImageFilter(false);
466 fixedImagePyramid->SetInput(croppedFixedImage);
467 fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
468 movingImagePyramid->SetUseShrinkImageFilter(false);
469 movingImagePyramid->SetInput(movingImage);
470 movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
471 if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
472 fixedImagePyramid->Update();
473 movingImagePyramid->Update();
474 typedef clitk::MultiResolutionPyramidRegionFilter<FixedImageType> FixedImageRegionPyramidType;
475 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
476 fixedImageRegionPyramid->SetRegion(metricRegion);
477 fixedImageRegionPyramid->SetSchedule(fixedImagePyramid->GetSchedule());
480 //=======================================================
481 // Rigid or Affine Transform
482 //=======================================================
483 typedef itk::AffineTransform <double,3> RigidTransformType;
484 RigidTransformType::Pointer rigidTransform=NULL;
485 if (m_ArgsInfo.initMatrix_given)
487 if(m_Verbose) std::cout<<"Reading the prior transform matrix "<< m_ArgsInfo.initMatrix_arg<<"..."<<std::endl;
488 rigidTransform=RigidTransformType::New();
489 itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.initMatrix_arg);
492 itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
493 rigidTransform->SetMatrix(finalRotation);
495 //Set the translation
496 itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
497 rigidTransform->SetTranslation(finalTranslation);
499 else if (m_ArgsInfo.centre_flag)
501 if(m_Verbose) std::cout<<"No itinial matrix given and \"centre\" flag switched on. Centering all images..."<<std::endl;
503 rigidTransform=RigidTransformType::New();
505 typedef itk::CenteredTransformInitializer<RigidTransformType, FixedImageType, MovingImageType > TransformInitializerType;
506 typename TransformInitializerType::Pointer initializer = TransformInitializerType::New();
507 initializer->SetTransform( rigidTransform );
508 initializer->SetFixedImage( fixedImage );
509 initializer->SetMovingImage( movingImage );
510 initializer->GeometryOn();
511 initializer->InitializeTransform();
515 //=======================================================
516 // B-LUT FFD Transform
517 //=======================================================
518 typedef clitk::MultipleBSplineDeformableTransform<TCoordRep,InputImageType::ImageDimension, SpaceDimension > TransformType;
519 typename TransformType::Pointer transform = TransformType::New();
520 if (labels) transform->SetLabels(labels);
521 if (rigidTransform) transform->SetBulkTransform(rigidTransform);
523 //-------------------------------------------------------------------------
524 // The transform initializer
525 //-------------------------------------------------------------------------
526 typedef clitk::MultipleBSplineDeformableTransformInitializer< TransformType,FixedImageType> InitializerType;
527 typename InitializerType::Pointer initializer = InitializerType::New();
528 initializer->SetVerbose(m_Verbose);
529 initializer->SetImage(fixedImagePyramid->GetOutput(0));
530 initializer->SetTransform(transform);
532 //-------------------------------------------------------------------------
534 //-------------------------------------------------------------------------
535 typename FixedImageType::RegionType::SizeType splineOrders ;
536 splineOrders.Fill(3);
537 if (m_ArgsInfo.order_given)
538 for(unsigned int i=0; i<InputImageType::ImageDimension;i++)
539 splineOrders[i]=m_ArgsInfo.order_arg[i];
540 if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
541 initializer->SetSplineOrders(splineOrders);
543 //-------------------------------------------------------------------------
545 //-------------------------------------------------------------------------
548 if (m_ArgsInfo.spacing_given)
550 initializer->m_ControlPointSpacingArray.resize(m_ArgsInfo.levels_arg);
551 initializer->SetControlPointSpacing(m_ArgsInfo.spacing_arg);
552 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]=initializer->m_ControlPointSpacing;
553 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]
554 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
556 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
558 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]=initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-i]*2;
559 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]
560 <<" at level "<<m_ArgsInfo.levels_arg-i<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
566 if (m_ArgsInfo.control_given)
568 initializer->m_NumberOfControlPointsInsideTheImageArray.resize(m_ArgsInfo.levels_arg);
569 initializer->SetNumberOfControlPointsInsideTheImage(m_ArgsInfo.control_arg);
570 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]=initializer->m_NumberOfControlPointsInsideTheImage;
571 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]<<"control points inside the image"
572 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
574 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
576 for(unsigned int j=0;j<InputImageType::ImageDimension;j++)
577 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i][j]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i][j]/2.);
578 // initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i]/2.);
579 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]<<"control points inside the image"
580 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
585 // Inialize on the first level
586 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level 1 of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
587 if (m_ArgsInfo.spacing_given) initializer->SetControlPointSpacing( initializer->m_ControlPointSpacingArray[0]);
588 if (m_ArgsInfo.control_given) initializer->SetNumberOfControlPointsInsideTheImage(initializer->m_NumberOfControlPointsInsideTheImageArray[0]);
589 if (m_ArgsInfo.samplingFactor_given) initializer->SetSamplingFactors(m_ArgsInfo.samplingFactor_arg);
592 initializer->InitializeTransform();
594 //-------------------------------------------------------------------------
595 // Initial parameters (passed by reference)
596 //-------------------------------------------------------------------------
597 typedef typename TransformType::ParametersType ParametersType;
598 const unsigned int numberOfParameters = transform->GetNumberOfParameters();
599 ParametersType parameters(numberOfParameters);
600 parameters.Fill( 0.0 );
601 transform->SetParameters( parameters );
602 if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters);
604 //-------------------------------------------------------------------------
605 // DEBUG: use an itk BSpline instead of multilabel BLUTs
606 //-------------------------------------------------------------------------
607 typedef itk::Transform< TCoordRep, 3, 3 > RegistrationTransformType;
608 RegistrationTransformType::Pointer regTransform(transform);
609 typedef itk::BSplineDeformableTransform<TCoordRep,SpaceDimension, 3> SingleBSplineTransformType;
610 typename SingleBSplineTransformType::Pointer sTransform;
611 if(m_ArgsInfo.itkbspline_flag) {
612 if( transform->GetTransforms().size()>1)
613 itkExceptionMacro(<< "invalid --itkbspline option if there is more than 1 label")
614 sTransform = SingleBSplineTransformType::New();
615 sTransform->SetBulkTransform( transform->GetTransforms()[0]->GetBulkTransform() );
616 sTransform->SetGridSpacing( transform->GetTransforms()[0]->GetGridSpacing() );
617 sTransform->SetGridOrigin( transform->GetTransforms()[0]->GetGridOrigin() );
618 sTransform->SetGridRegion( transform->GetTransforms()[0]->GetGridRegion() );
619 sTransform->SetParameters( transform->GetTransforms()[0]->GetParameters() );
620 regTransform = sTransform;
621 transform = NULL; // free memory
624 //=======================================================
626 //=======================================================
627 typedef clitk::GenericInterpolator<args_info_clitkBLUTDIR, FixedImageType, TCoordRep > GenericInterpolatorType;
628 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
629 genericInterpolator->SetArgsInfo(m_ArgsInfo);
630 typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
631 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
634 //=======================================================
636 //=======================================================
637 typedef clitk::GenericMetric<args_info_clitkBLUTDIR, FixedImageType,MovingImageType > GenericMetricType;
638 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
639 genericMetric->SetArgsInfo(m_ArgsInfo);
640 genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
641 if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
642 genericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(0));
643 typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
644 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
645 if (movingMask) metric->SetMovingImageMask(movingMask);
647 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
649 metric->SetNumberOfThreads( threads );
650 if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
653 if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
657 //=======================================================
659 //=======================================================
660 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
661 GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
662 genericOptimizer->SetArgsInfo(m_ArgsInfo);
663 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
664 genericOptimizer->SetNumberOfParameters(regTransform->GetNumberOfParameters());
665 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
666 OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
669 //=======================================================
671 //=======================================================
672 typedef itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;
673 typename RegistrationType::Pointer registration = RegistrationType::New();
674 registration->SetMetric( metric );
675 registration->SetOptimizer( optimizer );
676 registration->SetInterpolator( interpolator );
677 registration->SetTransform (regTransform );
679 registration->SetNumberOfThreads(threads);
680 if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
682 registration->SetFixedImage( croppedFixedImage );
683 registration->SetMovingImage( movingImage );
684 registration->SetFixedImageRegion( metricRegion );
685 registration->SetFixedImagePyramid( fixedImagePyramid );
686 registration->SetMovingImagePyramid( movingImagePyramid );
687 registration->SetInitialTransformParameters( regTransform->GetParameters() );
688 registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
689 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
692 //================================================================================================
694 //================================================================================================
697 // Output iteration info
698 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
699 observer->SetOptimizer(genericOptimizer);
700 optimizer->AddObserver( itk::IterationEvent(), observer );
703 typedef RegistrationInterfaceCommand<RegistrationType,args_info_clitkBLUTDIR> CommandType;
704 typename CommandType::Pointer command = CommandType::New();
705 command->SetInitializer(initializer);
706 command->SetArgsInfo(m_ArgsInfo);
707 command->SetCommandIterationUpdate(observer);
708 command->SetMaximize(genericMetric->GetMaximize());
709 command->SetMetricRegion(metricRegion);
710 registration->AddObserver( itk::IterationEvent(), command );
712 if (m_ArgsInfo.coeff_given)
714 if(m_ArgsInfo.itkbspline_flag) {
715 itkExceptionMacro("--coeff and --itkbpline are incompatible");
718 std::cout << std::endl << "Output coefficient images every " << m_ArgsInfo.coeffEveryN_arg << " iterations." << std::endl;
719 typedef CommandIterationUpdateDVF<FixedImageType, OptimizerType, TransformType> DVFCommandType;
720 typename DVFCommandType::Pointer observerdvf = DVFCommandType::New();
721 observerdvf->SetFixedImage(fixedImage);
722 observerdvf->SetTransform(transform);
723 observerdvf->SetArgsInfo(m_ArgsInfo);
724 optimizer->AddObserver( itk::IterationEvent(), observerdvf );
729 //=======================================================
731 //=======================================================
732 if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
736 #if ITK_VERSION_MAJOR < 4 || (ITK_VERSION_MAJOR == 4 && ITK_VERSION_MINOR <= 2)
737 registration->StartRegistration();
739 registration->Update();
742 catch( itk::ExceptionObject & err )
744 std::cerr << "ExceptionObject caught while registering!" << std::endl;
745 std::cerr << err << std::endl;
750 //=======================================================
752 //=======================================================
753 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
754 regTransform->SetParameters( finalParameters );
757 std::cout<<"Stop condition description: "
758 <<registration->GetOptimizer()->GetStopConditionDescription()<<std::endl;
762 //=======================================================
763 // Get the BSpline coefficient images and write them
764 //=======================================================
765 if (m_ArgsInfo.coeff_given)
767 typedef typename TransformType::CoefficientImageType CoefficientImageType;
768 std::vector<typename CoefficientImageType::Pointer> coefficientImages = transform->GetCoefficientImages();
769 typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
770 typename CoeffWriterType::Pointer coeffWriter = CoeffWriterType::New();
771 unsigned nLabels = transform->GetnLabels();
773 std::string fname(m_ArgsInfo.coeff_arg);
774 int dotpos = fname.length() - 1;
775 while (dotpos >= 0 && fname[dotpos] != '.')
778 for (unsigned i = 0; i < nLabels; ++i)
780 std::ostringstream osfname;
781 osfname << fname.substr(0, dotpos) << '_' << i << fname.substr(dotpos);
782 coeffWriter->SetInput(coefficientImages[i]);
783 coeffWriter->SetFileName(osfname.str());
784 coeffWriter->Update();
790 //=======================================================
791 // Compute the DVF (only deformable transform)
792 //=======================================================
793 typedef itk::Vector< float, SpaceDimension > DisplacementType;
794 typedef itk::Image< DisplacementType, InputImageType::ImageDimension > DisplacementFieldType;
795 #if ITK_VERSION_MAJOR >= 4
796 typedef itk::TransformToDisplacementFieldSource<DisplacementFieldType, double> ConvertorType;
798 typedef itk::TransformToDeformationFieldSource<DisplacementFieldType, double> ConvertorType;
800 typename ConvertorType::Pointer filter= ConvertorType::New();
801 filter->SetNumberOfThreads(1);
802 if(m_ArgsInfo.itkbspline_flag)
803 sTransform->SetBulkTransform(NULL);
805 transform->SetBulkTransform(NULL);
806 filter->SetTransform(regTransform);
807 filter->SetOutputParametersFromImage(fixedImage);
809 typename DisplacementFieldType::Pointer field = filter->GetOutput();
812 //=======================================================
814 //=======================================================
815 typedef itk::ImageFileWriter< DisplacementFieldType > FieldWriterType;
816 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
817 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
818 fieldWriter->SetInput( field );
821 fieldWriter->Update();
823 catch( itk::ExceptionObject & excp )
825 std::cerr << "Exception thrown writing the DVF" << std::endl;
826 std::cerr << excp << std::endl;
831 //=======================================================
832 // Resample the moving image
833 //=======================================================
834 typedef itk::ResampleImageFilter< MovingImageType, FixedImageType > ResampleFilterType;
835 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
836 if (rigidTransform) {
837 if(m_ArgsInfo.itkbspline_flag)
838 sTransform->SetBulkTransform(rigidTransform);
840 transform->SetBulkTransform(rigidTransform);
842 resampler->SetTransform( regTransform );
843 resampler->SetInput( movingImage);
844 resampler->SetOutputParametersFromImage(fixedImage);
846 typename FixedImageType::Pointer result=resampler->GetOutput();
848 // typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType;
849 // typename WarpFilterType::Pointer warp = WarpFilterType::New();
851 // warp->SetDeformationField( field );
852 // warp->SetInput( movingImageReader->GetOutput() );
853 // warp->SetOutputOrigin( fixedImage->GetOrigin() );
854 // warp->SetOutputSpacing( fixedImage->GetSpacing() );
855 // warp->SetOutputDirection( fixedImage->GetDirection() );
856 // warp->SetEdgePaddingValue( 0.0 );
860 //=======================================================
861 // Write the warped image
862 //=======================================================
863 typedef itk::ImageFileWriter< FixedImageType > WriterType;
864 typename WriterType::Pointer writer = WriterType::New();
865 writer->SetFileName( m_ArgsInfo.output_arg );
866 writer->SetInput( result );
872 catch( itk::ExceptionObject & err )
874 std::cerr << "ExceptionObject caught !" << std::endl;
875 std::cerr << err << std::endl;
880 //=======================================================
881 // Calculate the difference after the deformable transform
882 //=======================================================
883 typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType;
884 if (m_ArgsInfo.after_given)
886 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
887 difference->SetValidInput( fixedImage );
888 difference->SetTestInput( result );
892 difference->Update();
894 catch( itk::ExceptionObject & err )
896 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
897 std::cerr << err << std::endl;
901 typename WriterType::Pointer differenceWriter=WriterType::New();
902 differenceWriter->SetInput(difference->GetOutput());
903 differenceWriter->SetFileName(m_ArgsInfo.after_arg);
904 differenceWriter->Update();
909 //=======================================================
910 // Calculate the difference before the deformable transform
911 //=======================================================
912 if( m_ArgsInfo.before_given )
915 typename FixedImageType::Pointer moving=FixedImageType::New();
916 if (m_ArgsInfo.initMatrix_given)
918 typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
919 typename ResamplerType::Pointer resampler=ResamplerType::New();
920 resampler->SetInput(movingImage);
921 resampler->SetOutputOrigin(fixedImage->GetOrigin());
922 resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
923 resampler->SetOutputSpacing(fixedImage->GetSpacing());
924 resampler->SetDefaultPixelValue( 0. );
925 if (rigidTransform ) resampler->SetTransform(rigidTransform);
927 moving=resampler->GetOutput();
932 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
933 difference->SetValidInput( fixedImage );
934 difference->SetTestInput( moving );
938 difference->Update();
940 catch( itk::ExceptionObject & err )
942 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
943 std::cerr << err << std::endl;
947 typename WriterType::Pointer differenceWriter=WriterType::New();
948 writer->SetFileName( m_ArgsInfo.before_arg );
949 writer->SetInput( difference->GetOutput() );
958 #endif // #define clitkBLUTDIRGenericFilter_txx