added initial image centralization to BLUTDIR

[clitk.git] / registration / clitkBLUTDIRGenericFilter.cxx

diff --git a/registration/clitkBLUTDIRGenericFilter.cxx b/registration/clitkBLUTDIRGenericFilter.cxx
old mode 100755 (executable)
new mode 100644 (file)
index 15afacf..9f808e7
--- a/registration/clitkBLUTDIRGenericFilter.cxx
+++ b/registration/clitkBLUTDIRGenericFilter.cxx
@@ -29,7 +29,8 @@ It is distributed under dual licence
 
 #include "clitkBLUTDIRGenericFilter.h"
 #include "clitkBLUTDIRCommandIterationUpdateDVF.h"
-
+#include "itkCenteredTransformInitializer.h"
+  
 namespace clitk
 {
 
@@ -322,7 +323,6 @@ namespace clitk
       // The metric region with respect to the extracted transform region:
       // where should the metric be CALCULATED (depends on transform)
       typename FixedImageType::RegionType metricRegion = fixedImage->GetLargestPossibleRegion();
-      typename FixedImageType::RegionType::SizeType metricRegionSize=metricRegion.GetSize();
       typename FixedImageType::RegionType::IndexType metricRegionIndex=metricRegion.GetIndex();
       typename FixedImageType::PointType metricRegionOrigin=fixedImage->GetOrigin();
 
@@ -387,6 +387,9 @@ namespace clitk
         // Crop the fixedImage to the bounding box to facilitate multi-resolution
         typedef itk::ExtractImageFilter<FixedImageType,FixedImageType> ExtractImageFilterType;
         typename ExtractImageFilterType::Pointer extractImageFilter=ExtractImageFilterType::New();
+#if ITK_VERSION_MAJOR == 4
+        extractImageFilter->SetDirectionCollapseToSubmatrix();
+#endif
         extractImageFilter->SetInput(fixedImage);
         extractImageFilter->SetExtractionRegion(transformRegion);
         extractImageFilter->Update();
@@ -493,6 +496,20 @@ namespace clitk
         itk::Vector<double,3> finalTranslation = clitk::GetTranslationPartMatrix3D(rigidTransformMatrix);
         rigidTransform->SetTranslation(finalTranslation);
       }
+      else
+      {
+        if(m_Verbose) std::cout<<"No itinial matrix given. Centering all images..."<<std::endl;
+        
+        rigidTransform=RigidTransformType::New();
+        
+        typedef itk::CenteredTransformInitializer<RigidTransformType, FixedImageType, MovingImageType > TransformInitializerType;
+        typename TransformInitializerType::Pointer initializer = TransformInitializerType::New();
+        initializer->SetTransform( rigidTransform );
+        initializer->SetFixedImage( fixedImage );
+        initializer->SetMovingImage( movingImage );        
+        initializer->GeometryOn();
+        initializer->InitializeTransform();
+      }
 
 
       //=======================================================
@@ -584,6 +601,25 @@ namespace clitk
       transform->SetParameters( parameters );
       if (m_ArgsInfo.initCoeff_given) initializer->SetInitialParameters(m_ArgsInfo.initCoeff_arg, parameters);
 
+      //-------------------------------------------------------------------------
+      // DEBUG: use an itk BSpline instead of multilabel BLUTs
+      //-------------------------------------------------------------------------
+      typedef itk::Transform< TCoordRep, 3, 3 > RegistrationTransformType;
+      RegistrationTransformType::Pointer regTransform(transform);
+      typedef itk::BSplineDeformableTransform<TCoordRep,SpaceDimension, 3> SingleBSplineTransformType;
+      typename SingleBSplineTransformType::Pointer sTransform;
+      if(m_ArgsInfo.itkbspline_flag) {
+        if( transform->GetTransforms().size()>1)
+          itkExceptionMacro(<< "invalid --itkbspline option if there is more than 1 label")
+        sTransform = SingleBSplineTransformType::New();
+        sTransform->SetBulkTransform( transform->GetTransforms()[0]->GetBulkTransform() );
+        sTransform->SetGridSpacing( transform->GetTransforms()[0]->GetGridSpacing() );
+        sTransform->SetGridOrigin( transform->GetTransforms()[0]->GetGridOrigin() );
+        sTransform->SetGridRegion( transform->GetTransforms()[0]->GetGridRegion() );
+        sTransform->SetParameters( transform->GetTransforms()[0]->GetParameters() );
+        regTransform = sTransform;
+        transform = NULL; // free memory
+      }
 
       //=======================================================
       // Interpolator
@@ -608,7 +644,7 @@ namespace clitk
       typename  MetricType::Pointer metric=genericMetric->GetMetricPointer();
       if (movingMask) metric->SetMovingImageMask(movingMask);
 
-#ifdef ITK_USE_OPTIMIZED_REGISTRATION_METHODS
+#if defined(ITK_USE_OPTIMIZED_REGISTRATION_METHODS) || ITK_VERSION_MAJOR >= 4
       if (threadsGiven) {
         metric->SetNumberOfThreads( threads );
         if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
@@ -625,7 +661,7 @@ namespace clitk
       GenericOptimizerType::Pointer genericOptimizer = GenericOptimizerType::New();
       genericOptimizer->SetArgsInfo(m_ArgsInfo);
       genericOptimizer->SetMaximize(genericMetric->GetMaximize());
-      genericOptimizer->SetNumberOfParameters(transform->GetNumberOfParameters());
+      genericOptimizer->SetNumberOfParameters(regTransform->GetNumberOfParameters());
       typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
       OptimizerType::Pointer optimizer = genericOptimizer->GetOptimizerPointer();
 
@@ -638,7 +674,7 @@ namespace clitk
       registration->SetMetric(        metric        );
       registration->SetOptimizer(     optimizer     );
       registration->SetInterpolator(  interpolator  );
-      registration->SetTransform (transform);
+      registration->SetTransform (regTransform );
       if(threadsGiven) {
         registration->SetNumberOfThreads(threads);
         if (m_Verbose) std::cout<< "Using " << threads << " threads." << std::endl;
@@ -648,7 +684,7 @@ namespace clitk
       registration->SetFixedImageRegion( metricRegion );
       registration->SetFixedImagePyramid( fixedImagePyramid );
       registration->SetMovingImagePyramid( movingImagePyramid );
-      registration->SetInitialTransformParameters( transform->GetParameters() );
+      registration->SetInitialTransformParameters( regTransform->GetParameters() );
       registration->SetNumberOfLevels( m_ArgsInfo.levels_arg );
       if (m_Verbose) std::cout<<"Setting the number of resolution levels to "<<m_ArgsInfo.levels_arg<<"..."<<std::endl;
 
@@ -675,6 +711,10 @@ namespace clitk
 
         if (m_ArgsInfo.coeff_given)
         {
+          if(m_ArgsInfo.itkbspline_flag) {
+            itkExceptionMacro("--coeff and --itkbpline are incompatible");
+          }
+
           std::cout << std::endl << "Output coefficient images every " << m_ArgsInfo.coeffEveryN_arg << " iterations." << std::endl;
           typedef CommandIterationUpdateDVF<FixedImageType, OptimizerType, TransformType> DVFCommandType;
           typename DVFCommandType::Pointer observerdvf = DVFCommandType::New();
@@ -707,7 +747,7 @@ namespace clitk
       // Get the result
       //=======================================================
       OptimizerType::ParametersType finalParameters =  registration->GetLastTransformParameters();
-      transform->SetParameters( finalParameters );
+      regTransform->SetParameters( finalParameters );
       if (m_Verbose)
       {
         std::cout<<"Stop condition description: "
@@ -747,21 +787,28 @@ namespace clitk
       // Compute the DVF (only deformable transform)
       //=======================================================
       typedef itk::Vector< float, SpaceDimension >  DisplacementType;
-      typedef itk::Image< DisplacementType, InputImageType::ImageDimension >  DeformationFieldType;
-      typedef itk::TransformToDeformationFieldSource<DeformationFieldType, double> ConvertorType;
+      typedef itk::Image< DisplacementType, InputImageType::ImageDimension >  DisplacementFieldType;
+#if ITK_VERSION_MAJOR >= 4
+      typedef itk::TransformToDisplacementFieldSource<DisplacementFieldType, double> ConvertorType;
+#else
+      typedef itk::TransformToDeformationFieldSource<DisplacementFieldType, double> ConvertorType;
+#endif
       typename ConvertorType::Pointer filter= ConvertorType::New();
       filter->SetNumberOfThreads(1);
-      transform->SetBulkTransform(NULL);
-      filter->SetTransform(transform);
+      if(m_ArgsInfo.itkbspline_flag)
+        sTransform->SetBulkTransform(NULL);
+      else
+        transform->SetBulkTransform(NULL);
+      filter->SetTransform(regTransform);
       filter->SetOutputParametersFromImage(fixedImage);
       filter->Update();
-      typename DeformationFieldType::Pointer field = filter->GetOutput();
+      typename DisplacementFieldType::Pointer field = filter->GetOutput();
 
 
       //=======================================================
       // Write the DVF
       //=======================================================
-      typedef itk::ImageFileWriter< DeformationFieldType >  FieldWriterType;
+      typedef itk::ImageFileWriter< DisplacementFieldType >  FieldWriterType;
       typename FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
       fieldWriter->SetFileName( m_ArgsInfo.vf_arg );
       fieldWriter->SetInput( field );
@@ -782,8 +829,13 @@ namespace clitk
       //=======================================================
       typedef itk::ResampleImageFilter< MovingImageType, FixedImageType >    ResampleFilterType;
       typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
-      if (rigidTransform) transform->SetBulkTransform(rigidTransform);
-      resampler->SetTransform( transform );
+      if (rigidTransform) {
+        if(m_ArgsInfo.itkbspline_flag)
+          sTransform->SetBulkTransform(rigidTransform);
+        else
+          transform->SetBulkTransform(rigidTransform);
+      }
+      resampler->SetTransform( regTransform );
       resampler->SetInput( movingImage);
       resampler->SetOutputParametersFromImage(fixedImage);
       resampler->Update();