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://oncora1.lyon.fnclcc.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_txx
19 #define clitkBLUTDIRGenericFilter_txx
21 /* =================================================
22 * @file clitkBLUTDIRGenericFilter.txx
28 ===================================================*/
30 #include "clitkBLUTDIRGenericFilter.h"
35 //==============================================================================
36 // Creating an observer class that allows output at each iteration
37 //==============================================================================
38 class CommandIterationUpdate : public itk::Command
41 typedef CommandIterationUpdate Self;
42 typedef itk::Command Superclass;
43 typedef itk::SmartPointer<Self> Pointer;
46 CommandIterationUpdate() {};
48 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> OptimizerType;
49 typedef const OptimizerType * OptimizerPointer;
51 // Set the generic optimizer
52 void SetOptimizer(OptimizerPointer o){m_Optimizer=o;}
55 void Execute(itk::Object *caller, const itk::EventObject & event)
57 Execute( (const itk::Object *)caller, event);
60 void Execute(const itk::Object * object, const itk::EventObject & event)
62 if( !(itk::IterationEvent().CheckEvent( &event )) )
67 m_Optimizer->OutputIterationInfo();
70 OptimizerPointer m_Optimizer;
74 //==============================================================================
75 //Creating an observer class that allows us to change parameters at subsequent levels
76 //==============================================================================
77 template <typename TRegistration, class args_info_clitkBLUTDIR>
78 class RegistrationInterfaceCommand : public itk::Command
81 typedef RegistrationInterfaceCommand Self;
82 typedef itk::Command Superclass;
83 typedef itk::SmartPointer<Self> Pointer;
86 RegistrationInterfaceCommand() {};
90 typedef TRegistration RegistrationType;
91 typedef RegistrationType * RegistrationPointer;
94 typedef typename RegistrationType::FixedImageType FixedImageType;
95 typedef typename FixedImageType::RegionType RegionType;
96 itkStaticConstMacro(ImageDimension, unsigned int,FixedImageType::ImageDimension);
97 typedef clitk::BSplineDeformableTransform<double, ImageDimension, ImageDimension> TransformType;
98 typedef clitk::BSplineDeformableTransformInitializer<TransformType, FixedImageType> InitializerType;
99 typedef typename InitializerType::CoefficientImageType CoefficientImageType;
100 typedef itk::CastImageFilter<CoefficientImageType, CoefficientImageType> CastImageFilterType;
101 typedef typename TransformType::ParametersType ParametersType;
102 typedef typename InitializerType::Pointer InitializerPointer;
103 typedef CommandIterationUpdate::Pointer CommandIterationUpdatePointer;
106 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
107 typedef typename GenericOptimizerType::Pointer GenericOptimizerPointer;
110 typedef typename RegistrationType::FixedImageType InternalImageType;
111 typedef clitk::GenericMetric<args_info_clitkBLUTDIR, InternalImageType, InternalImageType> GenericMetricType;
112 typedef typename GenericMetricType::Pointer GenericMetricPointer;
114 // Two arguments are passed to the Execute() method: the first
115 // is the pointer to the object which invoked the event and the
116 // second is the event that was invoked.
117 void Execute(itk::Object * object, const itk::EventObject & event)
119 if( !(itk::IterationEvent().CheckEvent( &event )) )
125 RegistrationPointer registration = dynamic_cast<RegistrationPointer>( object );
126 unsigned int numberOfLevels=registration->GetNumberOfLevels();
127 unsigned int currentLevel=registration->GetCurrentLevel()+1;
130 std::cout<<std::endl;
131 std::cout<<"========================================"<<std::endl;
132 std::cout<<"Starting resolution level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
133 std::cout<<"========================================"<<std::endl;
134 std::cout<<std::endl;
139 // fixed image region pyramid
140 typedef clitk::MultiResolutionPyramidRegionFilter<InternalImageType> FixedImageRegionPyramidType;
141 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
142 fixedImageRegionPyramid->SetRegion(m_MetricRegion);
143 fixedImageRegionPyramid->SetSchedule(registration->GetFixedImagePyramid()->GetSchedule());
145 // Reinitialize the metric (!= number of samples)
146 m_GenericMetric= GenericMetricType::New();
147 m_GenericMetric->SetArgsInfo(m_ArgsInfo);
148 m_GenericMetric->SetFixedImage(registration->GetFixedImagePyramid()->GetOutput(registration->GetCurrentLevel()));
149 if (m_ArgsInfo.referenceMask_given) m_GenericMetric->SetFixedImageMask(registration->GetMetric()->GetFixedImageMask());
150 m_GenericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(registration->GetCurrentLevel()));
151 typedef itk::ImageToImageMetric< InternalImageType, InternalImageType > MetricType;
152 typename MetricType::Pointer metric=m_GenericMetric->GetMetricPointer();
153 registration->SetMetric(metric);
155 // Get the current coefficient image and make a COPY
156 typename itk::ImageDuplicator<CoefficientImageType>::Pointer caster=itk::ImageDuplicator<CoefficientImageType>::New();
157 caster->SetInputImage(m_Initializer->GetTransform()->GetCoefficientImage());
159 typename CoefficientImageType::Pointer currentCoefficientImage=caster->GetOutput();
161 // Write the intermediate result?
162 if (m_ArgsInfo.intermediate_given>=numberOfLevels)
163 writeImage<CoefficientImageType>(currentCoefficientImage, m_ArgsInfo.intermediate_arg[currentLevel-2], m_ArgsInfo.verbose_flag);
165 // Set the new transform properties
166 m_Initializer->SetImage(registration->GetFixedImagePyramid()->GetOutput(currentLevel-1));
167 if( m_Initializer->m_ControlPointSpacingIsGiven)
168 m_Initializer->SetControlPointSpacing(m_Initializer->m_ControlPointSpacingArray[registration->GetCurrentLevel()]);
169 if( m_Initializer->m_NumberOfControlPointsIsGiven)
170 m_Initializer->SetNumberOfControlPointsInsideTheImage(m_Initializer->m_NumberOfControlPointsInsideTheImageArray[registration->GetCurrentLevel()]);
172 // Reinitialize the transform
173 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level "<<currentLevel<<" of "<<numberOfLevels<<"..."<<std::endl;
174 m_Initializer->InitializeTransform();
175 ParametersType* newParameters= new typename TransformType::ParametersType(m_Initializer->GetTransform()->GetNumberOfParameters());
177 // Reinitialize an Optimizer (!= number of parameters)
178 m_GenericOptimizer = GenericOptimizerType::New();
179 m_GenericOptimizer->SetArgsInfo(m_ArgsInfo);
180 m_GenericOptimizer->SetMaximize(m_Maximize);
181 m_GenericOptimizer->SetNumberOfParameters(m_Initializer->GetTransform()->GetNumberOfParameters());
182 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
183 OptimizerType::Pointer optimizer = m_GenericOptimizer->GetOptimizerPointer();
184 optimizer->AddObserver( itk::IterationEvent(), m_CommandIterationUpdate);
185 registration->SetOptimizer(optimizer);
186 m_CommandIterationUpdate->SetOptimizer(m_GenericOptimizer);
188 // Set the previous transform parameters to the registration
189 // if(m_Initializer->m_Parameters!=NULL )delete m_Initializer->m_Parameters;
190 m_Initializer->SetInitialParameters(currentCoefficientImage,*newParameters);
191 registration->SetInitialTransformParametersOfNextLevel(*newParameters);
195 void Execute(const itk::Object * , const itk::EventObject & )
200 void SetInitializer(InitializerPointer i){m_Initializer=i;}
201 InitializerPointer m_Initializer;
203 void SetArgsInfo(args_info_clitkBLUTDIR a){m_ArgsInfo=a;}
204 args_info_clitkBLUTDIR m_ArgsInfo;
206 void SetCommandIterationUpdate(CommandIterationUpdatePointer c){m_CommandIterationUpdate=c;};
207 CommandIterationUpdatePointer m_CommandIterationUpdate;
209 GenericOptimizerPointer m_GenericOptimizer;
210 void SetMaximize(bool b){m_Maximize=b;}
213 GenericMetricPointer m_GenericMetric;
214 void SetMetricRegion(RegionType i){m_MetricRegion=i;}
215 RegionType m_MetricRegion;
221 //==============================================================================
222 // Update with the number of dimensions
223 //==============================================================================
224 template<unsigned int Dimension>
225 void BLUTDIRGenericFilter::UpdateWithDim(std::string PixelType)
227 if (m_Verbose) std::cout << "Image was detected to be "<<Dimension<<"D and "<< PixelType<<"..."<<std::endl;
229 if(PixelType == "short"){
230 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
231 UpdateWithDimAndPixelType<Dimension, signed short>();
233 // else if(PixelType == "unsigned_short"){
234 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
235 // UpdateWithDimAndPixelType<Dimension, unsigned short>();
238 // else if (PixelType == "unsigned_char"){
239 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
240 // UpdateWithDimAndPixelType<Dimension, unsigned char>();
243 // else if (PixelType == "char"){
244 // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
245 // UpdateWithDimAndPixelType<Dimension, signed char>();
248 if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
249 UpdateWithDimAndPixelType<Dimension, float>();
254 //==============================================================================
255 // Update with the number of dimensions and pixeltype
256 //==============================================================================
257 template<unsigned int ImageDimension, class PixelType>
259 BLUTDIRGenericFilter::UpdateWithDimAndPixelType()
263 //=============================================================================
265 //=============================================================================
266 bool threadsGiven=m_ArgsInfo.threads_given;
267 int threads=m_ArgsInfo.threads_arg;
269 typedef itk::Image< PixelType, ImageDimension > FixedImageType;
270 typedef itk::Image< PixelType, ImageDimension > MovingImageType;
271 const unsigned int SpaceDimension = ImageDimension;
272 typedef double TCoordRep;
275 //=======================================================
277 //=======================================================
278 typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
279 typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
281 typename FixedImageReaderType::Pointer fixedImageReader = FixedImageReaderType::New();
282 typename MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();
284 fixedImageReader->SetFileName( m_ArgsInfo.reference_arg );
285 movingImageReader->SetFileName( m_ArgsInfo.target_arg );
286 if (m_Verbose) std::cout<<"Reading images..."<<std::endl;
287 fixedImageReader->Update();
288 movingImageReader->Update();
290 typename FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
291 typename MovingImageType::Pointer movingImage =movingImageReader->GetOutput();
292 typename FixedImageType::Pointer croppedFixedImage=fixedImage;
295 //=======================================================
297 //=======================================================
299 // The original input region
300 typename FixedImageType::RegionType fixedImageRegion = fixedImage->GetLargestPossibleRegion();
302 // The transform region with respect to the input region:
303 // where should the transform be DEFINED (depends on mask)
304 typename FixedImageType::RegionType transformRegion = fixedImage->GetLargestPossibleRegion();
305 typename FixedImageType::RegionType::SizeType transformRegionSize=transformRegion.GetSize();
306 typename FixedImageType::RegionType::IndexType transformRegionIndex=transformRegion.GetIndex();
307 typename FixedImageType::PointType transformRegionOrigin=fixedImage->GetOrigin();
309 // The metric region with respect to the extracted transform region:
310 // where should the metric be CALCULATED (depends on transform)
311 typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
312 typename FixedImageType::RegionType::SizeType metricRegionSize=metricRegion.GetSize();
313 typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
314 typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
317 //===========================================================================
318 // If given, we connect a mask to reference or target
319 //============================================================================
320 typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
321 typename MaskType::Pointer fixedMask=NULL;
322 if (m_ArgsInfo.referenceMask_given)
324 fixedMask= MaskType::New();
325 typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
326 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
327 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
328 maskReader->SetFileName(m_ArgsInfo.referenceMask_arg);
331 maskReader->Update();
333 catch( itk::ExceptionObject & err )
335 std::cerr << "ExceptionObject caught while reading mask !" << std::endl;
336 std::cerr << err << std::endl;
339 if (m_Verbose)std::cout <<"Reference image mask was read..." <<std::endl;
341 // Set the image to the spatialObject
342 fixedMask->SetImage( maskReader->GetOutput() );
344 // Find the bounding box of the "inside" label
345 typedef itk::LabelStatisticsImageFilter<ImageMaskType, ImageMaskType> StatisticsImageFilterType;
346 typename StatisticsImageFilterType::Pointer statisticsImageFilter=StatisticsImageFilterType::New();
347 statisticsImageFilter->SetInput(maskReader->GetOutput());
348 statisticsImageFilter->SetLabelInput(maskReader->GetOutput());
349 statisticsImageFilter->Update();
350 typename StatisticsImageFilterType::BoundingBoxType boundingBox = statisticsImageFilter->GetBoundingBox(1);
352 // Limit the transform region to the mask
353 for (unsigned int i=0; i<ImageDimension; i++)
355 transformRegionIndex[i]=boundingBox[2*i];
356 transformRegionSize[i]=boundingBox[2*i+1]-boundingBox[2*i]+1;
358 transformRegion.SetSize(transformRegionSize);
359 transformRegion.SetIndex(transformRegionIndex);
360 fixedImage->TransformIndexToPhysicalPoint(transformRegion.GetIndex(), transformRegionOrigin);
362 // Crop the fixedImage to the bounding box to facilitate multi-resolution
363 typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
364 typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
365 extractImageFilter->SetInput(fixedImage);
366 extractImageFilter->SetExtractionRegion(transformRegion);
367 extractImageFilter->Update();
368 croppedFixedImage=extractImageFilter->GetOutput();
370 // Update the metric region
371 metricRegion = croppedFixedImage->GetLargestPossibleRegion();
372 metricRegionIndex=metricRegion.GetIndex();
373 croppedFixedImage->TransformIndexToPhysicalPoint(metricRegionIndex, metricRegionOrigin);
375 // Set start index to zero (with respect to croppedFixedImage/transform region)
376 metricRegionIndex.Fill(0);
377 metricRegion.SetIndex(metricRegionIndex);
378 croppedFixedImage->SetRegions(metricRegion);
379 croppedFixedImage->SetOrigin(metricRegionOrigin);
382 typedef itk::ImageMaskSpatialObject< ImageDimension > MaskType;
383 typename MaskType::Pointer movingMask=NULL;
384 if (m_ArgsInfo.targetMask_given)
386 movingMask= MaskType::New();
387 typedef itk::Image< unsigned char, ImageDimension > ImageMaskType;
388 typedef itk::ImageFileReader< ImageMaskType > MaskReaderType;
389 typename MaskReaderType::Pointer maskReader = MaskReaderType::New();
390 maskReader->SetFileName(m_ArgsInfo.targetMask_arg);
393 maskReader->Update();
395 catch( itk::ExceptionObject & err )
397 std::cerr << "ExceptionObject caught !" << std::endl;
398 std::cerr << err << std::endl;
400 if (m_Verbose)std::cout <<"Target image mask was read..." <<std::endl;
402 movingMask->SetImage( maskReader->GetOutput() );
406 //=======================================================
408 //=======================================================
412 // Fixed image region
413 std::cout<<"The fixed image has its origin at "<<fixedImage->GetOrigin()<<std::endl
414 <<"The fixed image region starts at index "<<fixedImageRegion.GetIndex()<<std::endl
415 <<"The fixed image region has size "<< fixedImageRegion.GetSize()<<std::endl;
418 std::cout<<"The transform has its origin at "<<transformRegionOrigin<<std::endl
419 <<"The transform region will start at index "<<transformRegion.GetIndex()<<std::endl
420 <<"The transform region has size "<< transformRegion.GetSize()<<std::endl;
423 std::cout<<"The metric region has its origin at "<<metricRegionOrigin<<std::endl
424 <<"The metric region will start at index "<<metricRegion.GetIndex()<<std::endl
425 <<"The metric region has size "<< metricRegion.GetSize()<<std::endl;
430 //=======================================================
431 // Pyramids (update them for transform initializer)
432 //=======================================================
433 typedef itk::RecursiveMultiResolutionPyramidImageFilter< FixedImageType, FixedImageType> FixedImagePyramidType;
434 typedef itk::RecursiveMultiResolutionPyramidImageFilter< MovingImageType, MovingImageType> MovingImagePyramidType;
435 typename FixedImagePyramidType::Pointer fixedImagePyramid = FixedImagePyramidType::New();
436 typename MovingImagePyramidType::Pointer movingImagePyramid = MovingImagePyramidType::New();
437 fixedImagePyramid->SetUseShrinkImageFilter(false);
438 fixedImagePyramid->SetInput(croppedFixedImage);
439 fixedImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
440 movingImagePyramid->SetUseShrinkImageFilter(false);
441 movingImagePyramid->SetInput(movingImage);
442 movingImagePyramid->SetNumberOfLevels(m_ArgsInfo.levels_arg);
443 if (m_Verbose) std::cout<<"Creating the image pyramid..."<<std::endl;
444 fixedImagePyramid->Update();
445 movingImagePyramid->Update();
446 typedef clitk::MultiResolutionPyramidRegionFilter<FixedImageType> FixedImageRegionPyramidType;
447 typename FixedImageRegionPyramidType::Pointer fixedImageRegionPyramid=FixedImageRegionPyramidType::New();
448 fixedImageRegionPyramid->SetRegion(metricRegion);
449 fixedImageRegionPyramid->SetSchedule(fixedImagePyramid->GetSchedule());
452 //=======================================================
453 // Rigid or Affine Transform
454 //=======================================================
455 typedef itk::AffineTransform <double,3> RigidTransformType;
456 RigidTransformType::Pointer rigidTransform=NULL;
457 if (m_ArgsInfo.initMatrix_given)
459 if(m_Verbose) std::cout<<"Reading the prior transform matrix "<< m_ArgsInfo.initMatrix_arg<<"..."<<std::endl;
460 rigidTransform=RigidTransformType::New();
461 itk::Matrix<double,4,4> rigidTransformMatrix=clitk::ReadMatrix3D(m_ArgsInfo.initMatrix_arg);
464 itk::Matrix<double,3,3> finalRotation = clitk::GetRotationalPartMatrix3D(rigidTransformMatrix);
465 rigidTransform->SetMatrix(finalRotation);
467 //Set the translation
468 itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
469 rigidTransform->SetTranslation(finalTranslation);
473 //=======================================================
474 // B-LUT FFD Transform
475 //=======================================================
476 typedef clitk::BSplineDeformableTransform<TCoordRep,ImageDimension, SpaceDimension > TransformType;
477 typename TransformType::Pointer transform= TransformType::New();
478 if (fixedMask) transform->SetMask( fixedMask );
479 if (rigidTransform) transform->SetBulkTransform( rigidTransform );
481 //-------------------------------------------------------------------------
482 // The transform initializer
483 //-------------------------------------------------------------------------
484 typedef clitk::BSplineDeformableTransformInitializer< TransformType,FixedImageType> InitializerType;
485 typename InitializerType::Pointer initializer = InitializerType::New();
486 initializer->SetVerbose(m_Verbose);
487 initializer->SetImage(fixedImagePyramid->GetOutput(0));
488 initializer->SetTransform(transform);
490 //-------------------------------------------------------------------------
492 //-------------------------------------------------------------------------
493 typename FixedImageType::RegionType::SizeType splineOrders ;
494 splineOrders.Fill(3);
495 if (m_ArgsInfo.order_given)
496 for(unsigned int i=0; i<ImageDimension;i++)
497 splineOrders[i]=m_ArgsInfo.order_arg[i];
498 if (m_Verbose) std::cout<<"Setting the spline orders to "<<splineOrders<<"..."<<std::endl;
499 initializer->SetSplineOrders(splineOrders);
501 //-------------------------------------------------------------------------
503 //-------------------------------------------------------------------------
506 if (m_ArgsInfo.spacing_given)
508 initializer->m_ControlPointSpacingArray.resize(m_ArgsInfo.levels_arg);
509 initializer->SetControlPointSpacing(m_ArgsInfo.spacing_arg);
510 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]=initializer->m_ControlPointSpacing;
511 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1]
512 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
514 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
516 initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]=initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-i]*2;
517 if (m_Verbose) std::cout<<"Using a control point spacing of "<<initializer->m_ControlPointSpacingArray[m_ArgsInfo.levels_arg-1-i]
518 <<" at level "<<m_ArgsInfo.levels_arg-i<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
524 if (m_ArgsInfo.control_given)
526 initializer->m_NumberOfControlPointsInsideTheImageArray.resize(m_ArgsInfo.levels_arg);
527 initializer->SetNumberOfControlPointsInsideTheImage(m_ArgsInfo.control_arg);
528 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]=initializer->m_NumberOfControlPointsInsideTheImage;
529 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1]<<"control points inside the image"
530 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
532 for (int i=1; i<m_ArgsInfo.levels_arg; i++ )
534 for(unsigned int j=0;j<ImageDimension;j++)
535 initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i][j]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i][j]/2.);
536 // initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]=ceil ((double)initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-i]/2.);
537 if (m_Verbose) std::cout<<"Using "<< initializer->m_NumberOfControlPointsInsideTheImageArray[m_ArgsInfo.levels_arg-1-i]<<"control points inside the image"
538 <<" at level "<<m_ArgsInfo.levels_arg<<" of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
543 // Inialize on the first level
544 if (m_ArgsInfo.verbose_flag) std::cout<<"Initializing transform for level 1 of "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
545 if (m_ArgsInfo.spacing_given) initializer->SetControlPointSpacing( initializer->m_ControlPointSpacingArray[0]);
546 if (m_ArgsInfo.control_given) initializer->SetNumberOfControlPointsInsideTheImage(initializer->m_NumberOfControlPointsInsideTheImageArray[0]);
547 if (m_ArgsInfo.samplingFactor_given) initializer->SetSamplingFactors(m_ArgsInfo.samplingFactor_arg);
550 initializer->InitializeTransform();
552 //-------------------------------------------------------------------------
553 // Initial parameters (passed by reference)
554 //-------------------------------------------------------------------------
555 typedef typename TransformType::ParametersType ParametersType;
556 const unsigned int numberOfParameters = transform->GetNumberOfParameters();
557 ParametersType parameters(numberOfParameters);
558 parameters.Fill( 0.0 );
559 transform->SetParameters( parameters );
560 if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters);
563 //=======================================================
565 //=======================================================
566 typedef clitk::GenericInterpolator<args_info_clitkBLUTDIR, FixedImageType, TCoordRep > GenericInterpolatorType;
567 typename GenericInterpolatorType::Pointer genericInterpolator=GenericInterpolatorType::New();
568 genericInterpolator->SetArgsInfo(m_ArgsInfo);
569 typedef itk::InterpolateImageFunction< FixedImageType, TCoordRep > InterpolatorType;
570 typename InterpolatorType::Pointer interpolator=genericInterpolator->GetInterpolatorPointer();
573 //=======================================================
575 //=======================================================
576 typedef clitk::GenericMetric<args_info_clitkBLUTDIR, FixedImageType,MovingImageType > GenericMetricType;
577 typename GenericMetricType::Pointer genericMetric=GenericMetricType::New();
578 genericMetric->SetArgsInfo(m_ArgsInfo);
579 genericMetric->SetFixedImage(fixedImagePyramid->GetOutput(0));
580 if (fixedMask) genericMetric->SetFixedImageMask(fixedMask);
581 genericMetric->SetFixedImageRegion(fixedImageRegionPyramid->GetOutput(0));
582 typedef itk::ImageToImageMetric< FixedImageType, MovingImageType > MetricType;
583 typename MetricType::Pointer metric=genericMetric->GetMetricPointer();
584 if (movingMask) metric->SetMovingImageMask(movingMask);
586 #ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
587 if (threadsGiven) metric->SetNumberOfThreads( threads );
589 if (m_Verbose) std::cout<<"Not setting the number of threads (not compiled with USE_OPTIMIZED_REGISTRATION_METHODS)..."<<std::endl;
593 //=======================================================
595 //=======================================================
596 typedef clitk::GenericOptimizer<args_info_clitkBLUTDIR> GenericOptimizerType;
597 GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
598 genericOptimizer->SetArgsInfo(m_ArgsInfo);
599 genericOptimizer->SetMaximize(genericMetric->GetMaximize());
600 genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters());
601 typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
602 OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
605 //=======================================================
607 //=======================================================
608 typedef itk::MultiResolutionImageRegistrationMethod< FixedImageType, MovingImageType > RegistrationType;
609 typename RegistrationType::Pointer registration = RegistrationType::New();
610 registration->SetMetric( metric );
611 registration->SetOptimizer( optimizer );
612 registration->SetInterpolator( interpolator );
613 registration->SetTransform (transform);
614 if(threadsGiven) registration->SetNumberOfThreads(threads);
615 registration->SetFixedImage( croppedFixedImage );
616 registration->SetMovingImage( movingImage );
617 registration->SetFixedImageRegion( metricRegion );
618 registration->SetFixedImagePyramid( fixedImagePyramid );
619 registration->SetMovingImagePyramid( movingImagePyramid );
620 registration->SetInitialTransformParameters( transform->GetParameters() );
621 registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
622 if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
625 //================================================================================================
627 //================================================================================================
630 // Output iteration info
631 CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
632 observer->SetOptimizer(genericOptimizer);
633 optimizer->AddObserver( itk::IterationEvent(), observer );
636 typedef RegistrationInterfaceCommand<RegistrationType,args_info_clitkBLUTDIR> CommandType;
637 typename CommandType::Pointer command = CommandType::New();
638 command->SetInitializer(initializer);
639 command->SetArgsInfo(m_ArgsInfo);
640 command->SetCommandIterationUpdate(observer);
641 command->SetMaximize(genericMetric->GetMaximize());
642 command->SetMetricRegion(metricRegion);
643 registration->AddObserver( itk::IterationEvent(), command );
647 //=======================================================
649 //=======================================================
650 if (m_Verbose) std::cout << std::endl << "Starting Registration" << std::endl;
654 registration->StartRegistration();
656 catch( itk::ExceptionObject & err )
658 std::cerr << "ExceptionObject caught while registering!" << std::endl;
659 std::cerr << err << std::endl;
664 //=======================================================
666 //=======================================================
667 OptimizerType::ParametersType finalParameters = registration->GetLastTransformParameters();
668 transform->SetParameters( finalParameters );
671 std::cout<<"Stop condition description: "
672 <<registration->GetOptimizer()->GetStopConditionDescription()<<std::endl;
676 //=======================================================
677 // Get the BSpline coefficient images and write them
678 //=======================================================
679 if (m_ArgsInfo.coeff_given)
681 typedef typename TransformType::CoefficientImageType CoefficientImageType;
682 typename CoefficientImageType::Pointer coefficientImage =transform->GetCoefficientImage();
683 typedef itk::ImageFileWriter<CoefficientImageType> CoeffWriterType;
684 typename CoeffWriterType::Pointer coeffWriter=CoeffWriterType::New();
685 coeffWriter->SetInput(coefficientImage);
686 coeffWriter->SetFileName(m_ArgsInfo.coeff_arg);
687 coeffWriter->Update();
692 //=======================================================
693 // Compute the DVF (only deformable transform)
694 //=======================================================
695 typedef itk::Vector< float, SpaceDimension > DisplacementType;
696 typedef itk::Image< DisplacementType, ImageDimension > DeformationFieldType;
697 typedef itk::TransformToDeformationFieldSource<DeformationFieldType, double> ConvertorType;
698 typename ConvertorType::Pointer filter= ConvertorType::New();
699 filter->SetNumberOfThreads(1);
700 transform->SetBulkTransform(NULL);
701 filter->SetTransform(transform);
702 filter->SetOutputParametersFromImage(fixedImage);
704 typename DeformationFieldType::Pointer field = filter->GetOutput();
707 //=======================================================
709 //=======================================================
710 typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
711 typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
712 fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
713 fieldWriter->SetInput( field );
716 fieldWriter->Update();
718 catch( itk::ExceptionObject & excp )
720 std::cerr << "Exception thrown writing the DVF" << std::endl;
721 std::cerr << excp << std::endl;
726 //=======================================================
727 // Resample the moving image
728 //=======================================================
729 typedef itk::ResampleImageFilter< MovingImageType, FixedImageType > ResampleFilterType;
730 typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
731 if (rigidTransform) transform->SetBulkTransform(rigidTransform);
732 resampler->SetTransform( transform );
733 resampler->SetInput( movingImageReader->GetOutput() );
734 resampler->SetOutputParametersFromImage(fixedImage);
736 typename FixedImageType::Pointer result=resampler->GetOutput();
738 // typedef itk::WarpImageFilter< MovingImageType, FixedImageType, DeformationFieldType > WarpFilterType;
739 // typename WarpFilterType::Pointer warp = WarpFilterType::New();
741 // warp->SetDeformationField( field );
742 // warp->SetInput( movingImageReader->GetOutput() );
743 // warp->SetOutputOrigin( fixedImage->GetOrigin() );
744 // warp->SetOutputSpacing( fixedImage->GetSpacing() );
745 // warp->SetOutputDirection( fixedImage->GetDirection() );
746 // warp->SetEdgePaddingValue( 0.0 );
750 //=======================================================
751 // Write the warped image
752 //=======================================================
753 typedef itk::ImageFileWriter< FixedImageType > WriterType;
754 typename WriterType::Pointer writer = WriterType::New();
755 writer->SetFileName( m_ArgsInfo.output_arg );
756 writer->SetInput( result );
762 catch( itk::ExceptionObject & err )
764 std::cerr << "ExceptionObject caught !" << std::endl;
765 std::cerr << err << std::endl;
770 //=======================================================
771 // Calculate the difference after the deformable transform
772 //=======================================================
773 typedef clitk::DifferenceImageFilter< FixedImageType, FixedImageType> DifferenceFilterType;
774 if (m_ArgsInfo.after_given)
776 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
777 difference->SetValidInput( fixedImage );
778 difference->SetTestInput( result );
782 difference->Update();
784 catch( itk::ExceptionObject & err )
786 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
787 std::cerr << err << std::endl;
791 typename WriterType::Pointer differenceWriter=WriterType::New();
792 differenceWriter->SetInput(difference->GetOutput());
793 differenceWriter->SetFileName(m_ArgsInfo.after_arg);
794 differenceWriter->Update();
799 //=======================================================
800 // Calculate the difference before the deformable transform
801 //=======================================================
802 if( m_ArgsInfo.before_given )
805 typename FixedImageType::Pointer moving=FixedImageType::New();
806 if (m_ArgsInfo.initMatrix_given)
808 typedef itk::ResampleImageFilter<MovingImageType, FixedImageType> ResamplerType;
809 typename ResamplerType::Pointer resampler=ResamplerType::New();
810 resampler->SetInput(movingImage);
811 resampler->SetOutputOrigin(fixedImage->GetOrigin());
812 resampler->SetSize(fixedImage->GetLargestPossibleRegion().GetSize());
813 resampler->SetOutputSpacing(fixedImage->GetSpacing());
814 resampler->SetDefaultPixelValue( 0. );
815 if (rigidTransform ) resampler->SetTransform(rigidTransform);
817 moving=resampler->GetOutput();
822 typename DifferenceFilterType::Pointer difference = DifferenceFilterType::New();
823 difference->SetValidInput( fixedImage );
824 difference->SetTestInput( moving );
828 difference->Update();
830 catch( itk::ExceptionObject & err )
832 std::cerr << "ExceptionObject caught calculating the difference !" << std::endl;
833 std::cerr << err << std::endl;
837 typename WriterType::Pointer differenceWriter=WriterType::New();
838 writer->SetFileName( m_ArgsInfo.before_arg );
839 writer->SetInput( difference->GetOutput() );
848 #endif // #define clitkBLUTDIRGenericFilter_txx