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"
36 //==============================================================================
37 // Creating an observer class that allows output at each iteration
38 //==============================================================================
39 class CommandIterationUpdate : public itk::Command
42 typedef CommandIterationUpdate Self;
43 typedef itk::Command Superclass;
44 typedef itk::SmartPointer<Self> Pointer;
47 CommandIterationUpdate() {};
49 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> OptimizerType;
50 typedef const OptimizerType * OptimizerPointer;
52 // Set the generic optimizer
53 void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
56 void Execute(itk::Object *caller, const itk::EventObject & event)
58 Execute( (const itk::Object *)caller, event);
61 void Execute(const itk::Object * object, const itk::EventObject & event)
63 if( !(itk::IterationEvent().CheckEvent( &event )) )
68 m_Optimizer->OutputIterationInfo();
71 OptimizerPointer m_Optimizer;
74 //===========================================================================//
76 //==========================================================================//
77 BLUTDIRGenericFilter::BLUTDIRGenericFilter():
78 ImageToImageGenericFilter<Self>("Register DIR")
80 InitializeImageType<2>();
81 InitializeImageType<3>();
85 //=========================================================================//
87 //==========================================================================//
88 void BLUTDIRGenericFilter::SetArgsInfo(const args_info_clitkBLUTDIR & a){
90 if (m_ArgsInfo.reference_given) AddInputFilename(m_ArgsInfo.reference_arg);
92 if (m_ArgsInfo.target_given) {
93 AddInputFilename(m_ArgsInfo.target_arg);
96 if (m_ArgsInfo.output_given) SetOutputFilename(m_ArgsInfo.output_arg);
98 if (m_ArgsInfo.verbose_given) m_Verbose=true;
101 //=========================================================================//
102 //===========================================================================//
103 template<unsigned int Dim>
104 void BLUTDIRGenericFilter::InitializeImageType()
106 ADD_DEFAULT_IMAGE_TYPES(3);
108 //--------------------------------------------------------------------
110 //==============================================================================
111 //Creating an observer class that allows us to change parameters at subsequent levels
112 //==============================================================================
113 template <typename TRegistration,class args_info_clitkBLUTDIR>
114 class RegistrationInterfaceCommand : public itk::Command
117 typedef RegistrationInterfaceCommand Self;
118 typedef itk::Command Superclass;
119 typedef itk::SmartPointer<Self> Pointer;
122 RegistrationInterfaceCommand() { };
126 typedef TRegistration RegistrationType;
127 typedef RegistrationType * RegistrationPointer;
130 typedef typename RegistrationType::FixedImageType FixedImageType;
131 typedef typename FixedImageType::RegionType RegionType;
132 itkStaticConstMacro(ImageDimension, unsigned int,FixedImageType::ImageDimension);
133 typedef clitk::MultipleBSplineDeformableTransform<double, ImageDimension, ImageDimension> TransformType;
134 typedef clitk::MultipleBSplineDeformableTransformInitializer<TransformType, FixedImageType> InitializerType;
135 typedef typename InitializerType::CoefficientImageType CoefficientImageType;
136 typedef itk::CastImageFilter<CoefficientImageType, CoefficientImageType> CastImageFilterType;
137 typedef typename TransformType::ParametersType ParametersType;
138 typedef typename InitializerType::Pointer InitializerPointer;
139 typedef CommandIterationUpdate::Pointer CommandIterationUpdatePointer;
142 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
143 typedef typename GenericOptimizerType::Pointer GenericOptimizerPointer;
146 typedef typename RegistrationType::FixedImageType InternalImageType;
147 typedef clitk::GenericMetric<args_info_clitkBLUTDIR, InternalImageType, InternalImageType> GenericMetricType;
148 typedef typename GenericMetricType::Pointer GenericMetricPointer;
150 // Two arguments are passed to the Execute() method: the first
151 // is the pointer to the object which invoked the event and the
152 // second is the event that was invoked.
153 void Execute(itk::Object * object, const itk::EventObject & event)
155 if( !(itk::IterationEvent().CheckEvent( &event )) )
161 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
162 unsigned int numberOfLevels=registration->GetNumberOfLevels();
163 unsigned int currentLevel=registration->GetCurrentLevel()+1;
166 std::cout<<std::endl;
167 std::cout<<"========================================"<<std::endl;
168 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
169 std::cout<<"========================================"<<std::endl;
170 std::cout<<std::endl;
175 // fixed image region pyramid
176 typedef clitk::MultiResolutionPyramidRegionFilter<InternalImageType> FixedImageRegionPyramidType;
177 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
178 fixedImageRegionPyramid->SetRegion(m_MetricRegion);
179 fixedImageRegionPyramid->SetSchedule(registration->GetFixedImagePyramid()->GetSchedule());
181 // Reinitialize the metric (!= number of samples)
182 m_GenericMetric= GenericMetricType::New();
183 m_GenericMetric->SetArgsInfo(m_ArgsInfo);
184 m_GenericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
185 if (m_ArgsInfo.referenceMask_given) m_GenericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
186 m_GenericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(registration->GetCurrentLevel()));
187 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
188 typename MetricType::Pointer metric=m_GenericMetric->GetMetricPointer();
189 registration->SetMetric(metric);
191 // Get the current coefficient image and make a COPY
192 typename itk::ImageDuplicator<CoefficientImageType>::Pointer caster = itk::ImageDuplicator<CoefficientImageType>::New();
193 std::vector<typename CoefficientImageType::Pointer> currentCoefficientImages = m_Initializer->GetTransform()->GetCoefficientImages();
194 for (unsigned i = 0; i < currentCoefficientImages.size(); ++i)
196 caster->SetInputImage(currentCoefficientImages[i]);
198 currentCoefficientImages[i] = caster->GetOutput();
202 // Write the intermediate result?
203 if (m_ArgsInfo.intermediate_given>=numberOfLevels)
204 writeImage<CoefficientImageType>(currentCoefficientImage, m_ArgsInfo.intermediate_arg[currentLevel-2], m_ArgsInfo.verbose_flag);
207 // Set the new transform properties
208 m_Initializer->SetImage(registration->GetFixedImagePyramid()->GetOutput(currentLevel-1));
209 if( m_Initializer->m_ControlPointSpacingIsGiven)
210 m_Initializer->SetControlPointSpacing(m_Initializer->m_ControlPointSpacingArray[registration->GetCurrentLevel()]);
211 if( m_Initializer->m_NumberOfControlPointsIsGiven)
212 m_Initializer->SetNumberOfControlPointsInsideTheImage(m_Initializer->m_NumberOfControlPointsInsideTheImageArray[registration->GetCurrentLevel()]);
214 // Reinitialize the transform
215 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
216 m_Initializer->InitializeTransform();
217 ParametersType* newParameters= new typename TransformType::ParametersType(m_Initializer->GetTransform()->GetNumberOfParameters());
219 // DS : if we want to skip the last pyramid level, force to only 1 iteration
220 DD(m_ArgsInfo.skipLastPyramidLevel_flag);
221 if ((currentLevel == numberOfLevels) && (m_ArgsInfo.skipLastPyramidLevel_flag)) {
222 DD(m_ArgsInfo.maxIt_arg);
223 std::cout << "I skip the last pyramid level : set max iteration to 0" << std::endl;
224 m_ArgsInfo.maxIt_arg = 0;
225 DD(m_ArgsInfo.maxIt_arg);
228 // Reinitialize an Optimizer (!= number of parameters)
229 m_GenericOptimizer = GenericOptimizerType::New();
230 m_GenericOptimizer->SetArgsInfo(m_ArgsInfo);
231 m_GenericOptimizer->SetMaximize(m_Maximize);
232 m_GenericOptimizer->SetNumberOfParameters(m_Initializer->GetTransform()->GetNumberOfParameters());
235 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
236 OptimizerType::Pointer optimizer = m_GenericOptimizer->GetOptimizerPointer();
237 optimizer->AddObserver( itk::IterationEvent(), m_CommandIterationUpdate);
238 registration->SetOptimizer(optimizer);
239 m_CommandIterationUpdate->SetOptimizer(m_GenericOptimizer);
241 // Set the previous transform parameters to the registration
242 // if(m_Initializer->m_Parameters!=NULL )delete m_Initializer->m_Parameters;
243 m_Initializer->SetInitialParameters(currentCoefficientImages, *newParameters);
244 registration->SetInitialTransformParametersOfNextLevel(*newParameters);
248 void Execute(const itk::Object * , const itk::EventObject & )
253 void SetInitializer(InitializerPointer i){m_Initializer=i;}
254 InitializerPointer m_Initializer;
256 void SetArgsInfo(args_info_clitkBLUTDIR a){m_ArgsInfo=a;}
257 args_info_clitkBLUTDIR m_ArgsInfo;
259 void SetCommandIterationUpdate(CommandIterationUpdatePointer c){m_CommandIterationUpdate=c;};
260 CommandIterationUpdatePointer m_CommandIterationUpdate;
262 GenericOptimizerPointer m_GenericOptimizer;
263 void SetMaximize(bool b){m_Maximize=b;}
266 GenericMetricPointer m_GenericMetric;
267 void SetMetricRegion(RegionType i){m_MetricRegion=i;}
268 RegionType m_MetricRegion;
273 //==============================================================================
274 // Update with the number of dimensions and pixeltype
275 //==============================================================================
276 template<class InputImageType>
277 void BLUTDIRGenericFilter::UpdateWithInputImageType()
279 if (m_Verbose) std::cout << "BLUTDIRGenericFilter::UpdateWithInputImageType()" << std::endl;
281 //=============================================================================
283 //=============================================================================
284 bool threadsGiven=m_ArgsInfo.threads_given;
285 int threads=m_ArgsInfo.threads_arg;
286 typedef typename InputImageType::PixelType PixelType;
288 typedef double TCoordRep;
290 typename InputImageType::Pointer fixedImage = this->template GetInput<InputImageType>(0);
292 typename InputImageType::Pointer inputFixedImage = this->template GetInput<InputImageType>(0);
295 typename InputImageType::Pointer movingImage = this->template GetInput<InputImageType>(1);
297 typename InputImageType::Pointer inputMovingImage = this->template GetInput<InputImageType>(1);
299 typedef itk::Image< PixelType,InputImageType::ImageDimension > FixedImageType;
300 typedef itk::Image< PixelType, InputImageType::ImageDimension> MovingImageType;
301 const unsigned int SpaceDimension = InputImageType::ImageDimension;
302 //Whatever the pixel type, internally we work with an image represented in float
303 //Reading reference image
304 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
305 //=======================================================
307 //=======================================================
308 typename FixedImageType::Pointer croppedFixedImage=fixedImage;
309 //=======================================================
311 //=======================================================
312 // The original input region
313 typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
315 // The transform region with respect to the input region:
316 // where should the transform be DEFINED (depends on mask)
317 typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
318 typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
319 typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
320 typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
322 // The metric region with respect to the extracted transform region:
323 // where should the metric be CALCULATED (depends on transform)
324 typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
325 typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
326 typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
329 //===========================================================================
330 // If given, we connect a mask to reference or target
331 //============================================================================
332 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
333 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageLabelType;
334 typename MaskType::Pointer fixedMask = NULL;
335 typename ImageLabelType::Pointer labels = NULL;
336 if (m_ArgsInfo.referenceMask_given)
338 fixedMask = MaskType::New();
339 labels = ImageLabelType::New();
340 typedef itk::ImageFileReader< ImageLabelType > LabelReaderType;
341 typename LabelReaderType::Pointer labelReader = LabelReaderType::New();
342 labelReader->SetFileName(m_ArgsInfo.referenceMask_arg);
345 labelReader->Update();
347 catch( itk::ExceptionObject & err )
349 std::cerr << "ExceptionObject caught while reading mask or labels !" << std::endl;
350 std::cerr << err << std::endl;
353 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
356 typedef itk::ResampleImageFilter<ImageLabelType, ImageLabelType> ResamplerType;
357 typename ResamplerType::Pointer resampler = ResamplerType::New();
358 typedef itk::NearestNeighborInterpolateImageFunction<ImageLabelType, TCoordRep> InterpolatorType;
359 typename InterpolatorType::Pointer interpolator = InterpolatorType::New();
360 resampler->SetOutputParametersFromImage(fixedImage);
361 resampler->SetInterpolator(interpolator);
362 resampler->SetInput(labelReader->GetOutput());
364 labels = resampler->GetOutput();
366 // Set the image to the spatialObject
367 fixedMask->SetImage(labels);
369 // Find the bounding box of the "inside" label
370 typedef itk::LabelGeometryImageFilter<ImageLabelType> GeometryImageFilterType;
371 typename GeometryImageFilterType::Pointer geometryImageFilter=GeometryImageFilterType::New();
372 geometryImageFilter->SetInput(labels);
373 geometryImageFilter->Update();
374 typename GeometryImageFilterType::BoundingBoxType boundingBox = geometryImageFilter->GetBoundingBox(1);
376 // Limit the transform region to the mask
377 for (unsigned int i=0; i<InputImageType::ImageDimension; i++)
379 transformRegionIndex[i]=boundingBox[2*i];
380 transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
382 transformRegion.SetSize(transformRegionSize);
383 transformRegion.SetIndex(transformRegionIndex);
384 fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
386 // Crop the fixedImage to the bounding box to facilitate multi-resolution
387 typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
388 typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
389 #if ITK_VERSION_MAJOR == 4
390 extractImageFilter->SetDirectionCollapseToSubmatrix();
392 extractImageFilter->SetInput(fixedImage);
393 extractImageFilter->SetExtractionRegion(transformRegion);
394 extractImageFilter->Update();
395 croppedFixedImage=extractImageFilter->GetOutput();
397 // Update the metric region
398 metricRegion = croppedFixedImage->GetLargestPossibleRegion();
399 metricRegionIndex=metricRegion.GetIndex();
400 croppedFixedImage->TransformIndexToPhysicalPoint(metricRegionIndex, metricRegionOrigin);
402 // Set start index to zero (with respect to croppedFixedImage/transform region)
403 metricRegionIndex.Fill(0);
404 metricRegion.SetIndex(metricRegionIndex);
405 croppedFixedImage->SetRegions(metricRegion);
406 croppedFixedImage->SetOrigin(metricRegionOrigin);
409 typedef itk::ImageMaskSpatialObject< InputImageType::ImageDimension > MaskType;
410 typename MaskType::Pointer movingMask=NULL;
411 if (m_ArgsInfo.targetMask_given)
413 movingMask= MaskType::New();
414 typedef itk::Image< unsigned char, InputImageType::ImageDimension > ImageMaskType;
415 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
416 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
417 maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
420 maskReader->Update();
422 catch( itk::ExceptionObject & err )
424 std::cerr << "ExceptionObject caught !" << std::endl;
425 std::cerr << err << std::endl;
427 if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
429 movingMask->SetImage( maskReader->GetOutput() );
433 //=======================================================
435 //=======================================================
439 // Fixed image region
440 std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
441 <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
442 <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
445 std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
446 <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
447 <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
450 std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
451 <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
452 <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
457 //=======================================================
458 // Pyramids (update them for transform initializer)
459 //=======================================================
460 typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
461 typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
462 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
463 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
464 fixedImagePyramid->SetUseShrinkImageFilter(false);
465 fixedImagePyramid->SetInput(croppedFixedImage);
466 fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
467 movingImagePyramid->SetUseShrinkImageFilter(false);
468 movingImagePyramid->SetInput(movingImage);
469 movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
470 if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
471 fixedImagePyramid->Update();
472 movingImagePyramid->Update();
473 typedef clitk::MultiResolutionPyramidRegionFilter<FixedImageType> FixedImageRegionPyramidType;
474 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
475 fixedImageRegionPyramid->SetRegion(metricRegion);
476 fixedImageRegionPyramid->SetSchedule(fixedImagePyramid->GetSchedule());
479 //=======================================================
480 // Rigid or Affine Transform
481 //=======================================================
482 typedef itk::AffineTransform <double,3> RigidTransformType;
483 RigidTransformType::Pointer rigidTransform=NULL;
484 if (m_ArgsInfo.initMatrix_given)
486 if(m_Verbose) std::cout<<"Reading the prior transform matrix "<< m_ArgsInfo.initMatrix_arg<<"..."<<std::endl;
487 rigidTransform=RigidTransformType::New();
488 itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.initMatrix_arg);
491 itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
492 rigidTransform->SetMatrix(finalRotation);
494 //Set the translation
495 itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
496 rigidTransform->SetTranslation(finalTranslation);
500 //=======================================================
501 // B-LUT FFD Transform
502 //=======================================================
503 typedef clitk::MultipleBSplineDeformableTransform<TCoordRep,InputImageType::ImageDimension, SpaceDimension > TransformType;
504 typename TransformType::Pointer transform = TransformType::New();
505 if (labels) transform->SetLabels(labels);
506 if (rigidTransform) transform->SetBulkTransform(rigidTransform);
508 //-------------------------------------------------------------------------
509 // The transform initializer
510 //-------------------------------------------------------------------------
511 typedef clitk::MultipleBSplineDeformableTransformInitializer< TransformType,FixedImageType> InitializerType;
512 typename InitializerType::Pointer initializer = InitializerType::New();
513 initializer->SetVerbose(m_Verbose);
514 initializer->SetImage(fixedImagePyramid->GetOutput(0));
515 initializer->SetTransform(transform);
517 //-------------------------------------------------------------------------
519 //-------------------------------------------------------------------------
520 typename FixedImageType::RegionType::SizeType splineOrders ;
521 splineOrders.Fill(3);
522 if (m_ArgsInfo.order_given)
523 for(unsigned int i=0; i<InputImageType::ImageDimension;i++)
524 splineOrders[i]=m_ArgsInfo.order_arg[i];
525 if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
526 initializer->SetSplineOrders(splineOrders);
528 //-------------------------------------------------------------------------
530 //-------------------------------------------------------------------------
533 if (m_ArgsInfo.spacing_given)
535 initializer->m_ControlPointSpacingArray.resize(m_ArgsInfo.levels_arg);
536 initializer->SetControlPointSpacing(m_ArgsInfo.spacing_arg);
537 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]=initializer->m_ControlPointSpacing;
538 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]
539 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
541 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
543 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]=initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-i]*2;
544 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]
545 <<" at level "<<m_ArgsInfo.levels_arg-i<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
551 if (m_ArgsInfo.control_given)
553 initializer->m_NumberOfControlPointsInsideTheImageArray.resize(m_ArgsInfo.levels_arg);
554 initializer->SetNumberOfControlPointsInsideTheImage(m_ArgsInfo.control_arg);
555 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]=initializer->m_NumberOfControlPointsInsideTheImage;
556 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]<<"control points inside the image"
557 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
559 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
561 for(unsigned int j=0;j<InputImageType::ImageDimension;j++)
562 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i][j]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i][j]/2.);
563 // initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i]/2.);
564 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]<<"control points inside the image"
565 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
570 // Inialize on the first level
571 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level 1 of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
572 if (m_ArgsInfo.spacing_given) initializer->SetControlPointSpacing( initializer->m_ControlPointSpacingArray[0]);
573 if (m_ArgsInfo.control_given) initializer->SetNumberOfControlPointsInsideTheImage(initializer->m_NumberOfControlPointsInsideTheImageArray[0]);
574 if (m_ArgsInfo.samplingFactor_given) initializer->SetSamplingFactors(m_ArgsInfo.samplingFactor_arg);
577 initializer->InitializeTransform();
579 //-------------------------------------------------------------------------
580 // Initial parameters (passed by reference)
581 //-------------------------------------------------------------------------
582 typedef typename TransformType::ParametersType ParametersType;
583 const unsigned int numberOfParameters = transform->GetNumberOfParameters();
584 ParametersType parameters(numberOfParameters);
585 parameters.Fill( 0.0 );
586 transform->SetParameters( parameters );
587 if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters);
589 //-------------------------------------------------------------------------
590 // DEBUG: use an itk BSpline instead of multilabel BLUTs
591 //-------------------------------------------------------------------------
592 typedef itk::Transform< TCoordRep, 3, 3 > RegistrationTransformType;
593 RegistrationTransformType::Pointer regTransform(transform);
594 typedef itk::BSplineDeformableTransform<TCoordRep,SpaceDimension, 3> SingleBSplineTransformType;
595 typename SingleBSplineTransformType::Pointer sTransform;
596 if(m_ArgsInfo.itkbspline_flag) {
597 if( transform->GetTransforms().size()>1)
598 itkExceptionMacro(<< "invalid --itkbspline option if there is more than 1 label")
599 sTransform = SingleBSplineTransformType::New();
600 sTransform->SetBulkTransform( transform->GetTransforms()[0]->GetBulkTransform() );
601 sTransform->SetGridSpacing( transform->GetTransforms()[0]->GetGridSpacing() );
602 sTransform->SetGridOrigin( transform->GetTransforms()[0]->GetGridOrigin() );
603 sTransform->SetGridRegion( transform->GetTransforms()[0]->GetGridRegion() );
604 sTransform->SetParameters( transform->GetTransforms()[0]->GetParameters() );
605 regTransform = sTransform;
606 transform = NULL; // free memory
609 //=======================================================
611 //=======================================================
612 typedef clitk::GenericInterpolator<args_info_clitkBLUTDIR, FixedImageType, TCoordRep > GenericInterpolatorType;
613 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
614 genericInterpolator->SetArgsInfo(m_ArgsInfo);
615 typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
616 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
619 //=======================================================
621 //=======================================================
622 typedef clitk::GenericMetric<args_info_clitkBLUTDIR, FixedImageType,MovingImageType > GenericMetricType;
623 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
624 genericMetric->SetArgsInfo(m_ArgsInfo);
625 genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
626 if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
627 genericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(0));
628 typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
629 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
630 if (movingMask) metric->SetMovingImageMask(movingMask);
632 #if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
634 metric->SetNumberOfThreads( threads );
635 if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
638 if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
642 //=======================================================
644 //=======================================================
645 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
646 GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
647 genericOptimizer->SetArgsInfo(m_ArgsInfo);
648 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
649 genericOptimizer->SetNumberOfParameters(regTransform->GetNumberOfParameters());
650 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
651 OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
654 //=======================================================
656 //=======================================================
657 typedef itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;
658 typename RegistrationType::Pointer registration = RegistrationType::New();
659 registration->SetMetric( metric );
660 registration->SetOptimizer( optimizer );
661 registration->SetInterpolator( interpolator );
662 registration->SetTransform (regTransform );
664 registration->SetNumberOfThreads(threads);
665 if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
667 registration->SetFixedImage( croppedFixedImage );
668 registration->SetMovingImage( movingImage );
669 registration->SetFixedImageRegion( metricRegion );
670 registration->SetFixedImagePyramid( fixedImagePyramid );
671 registration->SetMovingImagePyramid( movingImagePyramid );
672 registration->SetInitialTransformParameters( regTransform->GetParameters() );
673 registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
674 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
677 //================================================================================================
679 //================================================================================================
682 // Output iteration info
683 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
684 observer->SetOptimizer(genericOptimizer);
685 optimizer->AddObserver( itk::IterationEvent(), observer );
688 typedef RegistrationInterfaceCommand<RegistrationType,args_info_clitkBLUTDIR> CommandType;
689 typename CommandType::Pointer command = CommandType::New();
690 command->SetInitializer(initializer);
691 command->SetArgsInfo(m_ArgsInfo);
692 command->SetCommandIterationUpdate(observer);
693 command->SetMaximize(genericMetric->GetMaximize());
694 command->SetMetricRegion(metricRegion);
695 registration->AddObserver( itk::IterationEvent(), command );
697 if (m_ArgsInfo.coeff_given)
699 if(m_ArgsInfo.itkbspline_flag) {
700 itkExceptionMacro("--coeff and --itkbpline are incompatible");
703 std::cout << std::endl << "Output coefficient images every " << m_ArgsInfo.coeffEveryN_arg << " iterations." << std::endl;
704 typedef CommandIterationUpdateDVF<FixedImageType, OptimizerType, TransformType> DVFCommandType;
705 typename DVFCommandType::Pointer observerdvf = DVFCommandType::New();
706 observerdvf->SetFixedImage(fixedImage);
707 observerdvf->SetTransform(transform);
708 observerdvf->SetArgsInfo(m_ArgsInfo);
709 optimizer->AddObserver( itk::IterationEvent(), observerdvf );
714 //=======================================================
716 //=======================================================
717 if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
721 registration->StartRegistration();
723 catch( itk::ExceptionObject & err )
725 std::cerr << "ExceptionObject caught while registering!" << std::endl;
726 std::cerr << err << std::endl;
731 //=======================================================
733 //=======================================================
734 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
735 regTransform->SetParameters( finalParameters );
738 std::cout<<"Stop condition description: "
739 <<registration->GetOptimizer()->GetStopConditionDescription()<<std::endl;
743 //=======================================================
744 // Get the BSpline coefficient images and write them
745 //=======================================================
746 if (m_ArgsInfo.coeff_given)
748 typedef typename TransformType::CoefficientImageType CoefficientImageType;
749 std::vector<typename CoefficientImageType::Pointer> coefficientImages = transform->GetCoefficientImages();
750 typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
751 typename CoeffWriterType::Pointer coeffWriter = CoeffWriterType::New();
752 unsigned nLabels = transform->GetnLabels();
754 std::string fname(m_ArgsInfo.coeff_arg);
755 int dotpos = fname.length() - 1;
756 while (dotpos >= 0 && fname[dotpos] != '.')
759 for (unsigned i = 0; i < nLabels; ++i)
761 std::ostringstream osfname;
762 osfname << fname.substr(0, dotpos) << '_' << i << fname.substr(dotpos);
763 coeffWriter->SetInput(coefficientImages[i]);
764 coeffWriter->SetFileName(osfname.str());
765 coeffWriter->Update();
771 //=======================================================
772 // Compute the DVF (only deformable transform)
773 //=======================================================
774 typedef itk::Vector< float, SpaceDimension > DisplacementType;
775 typedef itk::Image< DisplacementType, InputImageType::ImageDimension > DisplacementFieldType;
776 #if ITK_VERSION_MAJOR >= 4
777 typedef itk::TransformToDisplacementFieldSource<DisplacementFieldType, double> ConvertorType;
779 typedef itk::TransformToDeformationFieldSource<DisplacementFieldType, double> ConvertorType;
781 typename ConvertorType::Pointer filter= ConvertorType::New();
782 filter->SetNumberOfThreads(1);
783 if(m_ArgsInfo.itkbspline_flag)
784 sTransform->SetBulkTransform(NULL);
786 transform->SetBulkTransform(NULL);
787 filter->SetTransform(regTransform);
788 filter->SetOutputParametersFromImage(fixedImage);
790 typename DisplacementFieldType::Pointer field = filter->GetOutput();
793 //=======================================================
795 //=======================================================
796 typedef itk::ImageFileWriter< DisplacementFieldType > FieldWriterType;
797 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
798 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
799 fieldWriter->SetInput( field );
802 fieldWriter->Update();
804 catch( itk::ExceptionObject & excp )
806 std::cerr << "Exception thrown writing the DVF" << std::endl;
807 std::cerr << excp << std::endl;
812 //=======================================================
813 // Resample the moving image
814 //=======================================================
815 typedef itk::ResampleImageFilter< MovingImageType, FixedImageType > ResampleFilterType;
816 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
817 if (rigidTransform) {
818 if(m_ArgsInfo.itkbspline_flag)
819 sTransform->SetBulkTransform(rigidTransform);
821 transform->SetBulkTransform(rigidTransform);
823 resampler->SetTransform( regTransform );
824 resampler->SetInput( movingImage);
825 resampler->SetOutputParametersFromImage(fixedImage);
827 typename FixedImageType::Pointer result=resampler->GetOutput();
829 // typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType;
830 // typename WarpFilterType::Pointer warp = WarpFilterType::New();
832 // warp->SetDeformationField( field );
833 // warp->SetInput( movingImageReader->GetOutput() );
834 // warp->SetOutputOrigin( fixedImage->GetOrigin() );
835 // warp->SetOutputSpacing( fixedImage->GetSpacing() );
836 // warp->SetOutputDirection( fixedImage->GetDirection() );
837 // warp->SetEdgePaddingValue( 0.0 );
841 //=======================================================
842 // Write the warped image
843 //=======================================================
844 typedef itk::ImageFileWriter< FixedImageType > WriterType;
845 typename WriterType::Pointer writer = WriterType::New();
846 writer->SetFileName( m_ArgsInfo.output_arg );
847 writer->SetInput( result );
853 catch( itk::ExceptionObject & err )
855 std::cerr << "ExceptionObject caught !" << std::endl;
856 std::cerr << err << std::endl;
861 //=======================================================
862 // Calculate the difference after the deformable transform
863 //=======================================================
864 typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType;
865 if (m_ArgsInfo.after_given)
867 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
868 difference->SetValidInput( fixedImage );
869 difference->SetTestInput( result );
873 difference->Update();
875 catch( itk::ExceptionObject & err )
877 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
878 std::cerr << err << std::endl;
882 typename WriterType::Pointer differenceWriter=WriterType::New();
883 differenceWriter->SetInput(difference->GetOutput());
884 differenceWriter->SetFileName(m_ArgsInfo.after_arg);
885 differenceWriter->Update();
890 //=======================================================
891 // Calculate the difference before the deformable transform
892 //=======================================================
893 if( m_ArgsInfo.before_given )
896 typename FixedImageType::Pointer moving=FixedImageType::New();
897 if (m_ArgsInfo.initMatrix_given)
899 typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
900 typename ResamplerType::Pointer resampler=ResamplerType::New();
901 resampler->SetInput(movingImage);
902 resampler->SetOutputOrigin(fixedImage->GetOrigin());
903 resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
904 resampler->SetOutputSpacing(fixedImage->GetSpacing());
905 resampler->SetDefaultPixelValue( 0. );
906 if (rigidTransform ) resampler->SetTransform(rigidTransform);
908 moving=resampler->GetOutput();
913 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
914 difference->SetValidInput( fixedImage );
915 difference->SetTestInput( moving );
919 difference->Update();
921 catch( itk::ExceptionObject & err )
923 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
924 std::cerr << err << std::endl;
928 typename WriterType::Pointer differenceWriter=WriterType::New();
929 writer->SetFileName( m_ArgsInfo.before_arg );
930 writer->SetInput( difference->GetOutput() );
939 #endif // #define clitkBLUTDIRGenericFilter_txx