Invert VF from B-Spline coefficients

- Interpolate VF first, from coeffs, using --like image

diff --git a/tools/clitkInvertVF.ggo b/tools/clitkInvertVF.ggo
index 19c451a478bba22e6a86b204eb4a6627c17f961e..92970bb8dc460c5771f8b9d7cf1e04ff49d6f176 100644 (file)
--- a/tools/clitkInvertVF.ggo
+++ b/tools/clitkInvertVF.ggo
@@ -11,7 +11,7 @@ option "threads"              -       "Number of threads (default=min(8, #CPU))"                                                      in
 
 option "input"                 i       "Input VF filename"                                             string          yes
 option "output"                o       "Output VF filename"                                            string          yes
-option "type"                  t       "Type of filter: 0=clitk (fast forward splat using linear kernels), 1= itk (grid subsumpling with controllable precision)"      int             no      default="0"
+option "type"                  t       "Type of filter: 0=clitk (fast forward splat using linear kernels), 1=clitk (like 0, but input images are B-Spline coefficients), 2= itk (grid subsumpling with controllable precision)"        int             no      default="0"
 option "threadSafe"                    -       "Clitk: use thread safe algorithm"                              flag            off
 option "pad"                   p       "Clitk: edge padding value (1 or N number of values!, defautls to zeros)"               double          multiple        no
 option "sampling"              s       "Itk: subsampling factor"                                       int             no      default="20"

diff --git a/tools/clitkInvertVFGenericFilter.h b/tools/clitkInvertVFGenericFilter.h
index 012524c8fc266d817b9cda408cde4259def6d5d9..cdaf86fcd900e726783a8dafbd933ec3bd7197b5 100644 (file)
--- a/tools/clitkInvertVFGenericFilter.h
+++ b/tools/clitkInvertVFGenericFilter.h
@@ -99,6 +99,7 @@ namespace clitk
     //----------------------------------------  
     template <unsigned int Dimension>  void UpdateWithDim(std::string PixelType);
     template <unsigned int Dimension, class PixelType>  void UpdateWithDimAndPixelType();
+    template<class DisplacementFieldType> typename DisplacementFieldType::Pointer CoeffsToDVF(std::string fileName, std::string likeFileName);
 
 
     //----------------------------------------  

diff --git a/tools/clitkInvertVFGenericFilter.txx b/tools/clitkInvertVFGenericFilter.txx
index 558b34177d48357acd174f7c6fc14f669e20e428..5cd12c4ce202129c38ac67d59d5629fa590b58e1 100644 (file)
--- a/tools/clitkInvertVFGenericFilter.txx
+++ b/tools/clitkInvertVFGenericFilter.txx
@@ -19,6 +19,12 @@
 #define clitkInvertVFGenericFilter_txx
 
 #include "itkVectorResampleImageFilter.h"
+#include "clitkBSplineDeformableTransform.h"
+#if ITK_VERSION_MAJOR >= 4
+#include "itkTransformToDisplacementFieldSource.h"
+#else
+#include "itkTransformToDeformationFieldSource.h"
+#endif
 
 /* =================================================
  * @file   clitkInvertVFGenericFilter.txx
@@ -100,6 +106,65 @@ InvertVFGenericFilter<args_info_type>::UpdateWithDim(std::string PixelType)
   // }
 }
 
+//-------------------------------------------------------------------
+// Convert Coefficient image to DVF
+//-------------------------------------------------------------------
+template<class args_info_type>
+template<class DisplacementFieldType>
+typename DisplacementFieldType::Pointer
+InvertVFGenericFilter<args_info_type>::CoeffsToDVF(std::string fileName, std::string likeFileName)
+{
+  typedef clitk::BSplineDeformableTransform<double, DisplacementFieldType::ImageDimension, DisplacementFieldType::ImageDimension> TransformType;
+  typedef typename TransformType::CoefficientImageType CoefficientImageType;
+
+  typedef itk::ImageFileReader<CoefficientImageType> CoeffReaderType;
+  typename CoeffReaderType::Pointer reader = CoeffReaderType::New();
+  reader->SetFileName(fileName);
+  reader->Update();
+
+  typename TransformType::Pointer transform = TransformType::New();
+  transform->SetCoefficientImage(reader->GetOutput());
+  
+#if ITK_VERSION_MAJOR >= 4
+      typedef itk::TransformToDisplacementFieldSource<DisplacementFieldType, double> ConvertorType;
+#else
+      typedef itk::TransformToDeformationFieldSource<DisplacementFieldType, double> ConvertorType;
+#endif
+
+  typedef itk::ImageIOBase ImageIOType;
+  typename ImageIOType::Pointer imageIO = itk::ImageIOFactory::CreateImageIO(likeFileName.c_str(), itk::ImageIOFactory::ReadMode);
+  imageIO->SetFileName(likeFileName);
+  imageIO->ReadImageInformation();
+
+  typename ConvertorType::Pointer convertor= ConvertorType::New();
+  typename ConvertorType::SizeType output_size;
+  typename ConvertorType::SpacingType output_spacing;
+  typename ConvertorType::OriginType output_origin;
+  typename ConvertorType::DirectionType output_direction;
+  for (unsigned int i = 0; i < DisplacementFieldType::ImageDimension; i++) {
+    output_size[i] = imageIO->GetDimensions(i);
+    output_spacing[i] = imageIO->GetSpacing(i);
+    output_origin[i] = imageIO->GetOrigin(i);
+    for (unsigned int j = 0; j < DisplacementFieldType::ImageDimension; j++)
+      output_direction[i][j] = imageIO->GetDirection(i)[j];
+  }
+  
+  if (m_Verbose) {
+    std::cout << "Interpolating coefficients with grid:" << std::endl;
+    std::cout << output_size << output_spacing << std::endl;
+  }
+  
+  convertor->SetNumberOfThreads(1);
+  convertor->SetTransform(transform);
+  convertor->SetOutputOrigin(output_origin);
+  convertor->SetOutputSpacing(output_spacing);
+  convertor->SetOutputSize(output_size);
+  convertor->SetOutputDirection(output_direction);
+  convertor->Update();
+
+  return convertor->GetOutput();
+}
+
 
 //-------------------------------------------------------------------
 // Update with the number of dimensions and the pixeltype
@@ -172,6 +237,31 @@ InvertVFGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
   }
 
   case 1: {
+    // Create the InvertVFFilter
+    typedef clitk::InvertVFFilter<InputImageType,OutputImageType> FilterType;
+    typename FilterType::Pointer filter =FilterType::New();
+    if (m_ArgsInfo.like_given) {
+      filter->SetInput(CoeffsToDVF<OutputImageType>(m_InputFileName, m_ArgsInfo.like_arg));
+    }
+
+    filter->SetVerbose(m_Verbose);
+    if (m_ArgsInfo.threads_given) filter->SetNumberOfThreads(m_ArgsInfo.threads_arg);
+    if (m_ArgsInfo.pad_given) {
+      PixelType pad;
+      if (m_ArgsInfo.pad_given !=  (pad.GetNumberOfComponents()) )
+        pad.Fill(m_ArgsInfo.pad_arg[0]);
+      else
+        for(unsigned int i=0; i<Dimension; i++)
+          pad[i]=m_ArgsInfo.pad_arg[i];
+    }
+    filter->SetThreadSafe(m_ArgsInfo.threadSafe_flag);
+    filter->Update();
+    output=filter->GetOutput();
+
+    break;
+  }
+
+  case 2: {
     // Create the InverseDeformationFieldFilter
 #if ITK_VERSION_MAJOR >= 4
     typedef itk::InverseDisplacementFieldImageFilter<InputImageType,OutputImageType> FilterType;
@@ -190,6 +280,7 @@ InvertVFGenericFilter<args_info_type>::UpdateWithDimAndPixelType()
     break;
   }
 
+    
   }
 
   // Output