image warping using B-Spline coefficients

- uses the resolution of the input image for the output

diff --git a/tools/clitkWarpImage.cxx b/tools/clitkWarpImage.cxx
index 3360661b1a74a0a76a0c0d33164e54d9f0b67fb6..673a5839aef8999b48f03b5924e35b06f89f9d6c 100644 (file)
--- a/tools/clitkWarpImage.cxx
+++ b/tools/clitkWarpImage.cxx
@@ -40,6 +40,13 @@ int main(int argc, char * argv[])
   GGO(clitkWarpImage, args_info);
   CLITK_INIT;
 
+  if (args_info.coeff_given) {
+    if (args_info.verbose_flag)
+      std::cout<< "Using coefficient images, so forcing sapcing = 1.\n" << std::endl;
+    
+    args_info.spacing_arg = 1;
+  }
+
   // Filter
   clitk::WarpImageGenericFilter::Pointer genericFilter=clitk::WarpImageGenericFilter::New();
 

diff --git a/tools/clitkWarpImage.ggo b/tools/clitkWarpImage.ggo
index 8a572c4c9c6167a362e657196406511fe41d56ba..12aa18e43574cb0af95f07f0c4516202c604b431 100644 (file)
--- a/tools/clitkWarpImage.ggo
+++ b/tools/clitkWarpImage.ggo
@@ -10,7 +10,10 @@ section "I/O"
 
 option "input"         i       "Input image filename"            string        yes
 option "output"        o       "Output image filename"           string        yes
-option "vf"            -       "Vector field filename"           string                yes
+
+defgroup "DVFoption" groupdesc="an option of this group is required" required
+groupoption "vf"               -       "Vector field filename"           string                yes group="DVFoption"
+groupoption "coeff"   - "B-Spline coefficients filename"     string    yes group="DVFoption"
 
 
 section "Options"

diff --git a/tools/clitkWarpImageGenericFilter.h b/tools/clitkWarpImageGenericFilter.h
index 93eda73a5c22bb62498fbb3aafc2a273f371b3f3..9394d8d5798d06ea2af9489822e37edf1716610f 100644 (file)
--- a/tools/clitkWarpImageGenericFilter.h
+++ b/tools/clitkWarpImageGenericFilter.h
@@ -98,7 +98,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);
 
     //----------------------------------------  
     // Data members

diff --git a/tools/clitkWarpImageGenericFilter.txx b/tools/clitkWarpImageGenericFilter.txx
index 8544df2ba33bf859c827513c4219ac8e6301a42b..ce979aa330690a1e08375cfca90024cfefbf208e 100644 (file)
--- a/tools/clitkWarpImageGenericFilter.txx
+++ b/tools/clitkWarpImageGenericFilter.txx
@@ -27,6 +27,13 @@
  *
  ===================================================*/
 
+#include "itkVectorResampleImageFilter.h"
+#include "clitkBSplineDeformableTransform.h"
+#if ITK_VERSION_MAJOR >= 4
+#include "itkTransformToDisplacementFieldSource.h"
+#else
+#include "itkTransformToDeformationFieldSource.h"
+#endif
 
 namespace clitk
 {
@@ -64,6 +71,63 @@ WarpImageGenericFilter::UpdateWithDim(std::string PixelType)
   }
 }
 
+//-------------------------------------------------------------------
+// Convert Coefficient image to DVF
+//-------------------------------------------------------------------
+template<class DisplacementFieldType>
+typename DisplacementFieldType::Pointer
+WarpImageGenericFilter::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
@@ -78,21 +142,26 @@ WarpImageGenericFilter::UpdateWithDimAndPixelType()
   typedef itk::Image<PixelType, Dimension> OutputImageType;
   typedef itk::Vector<float, Dimension> DisplacementType;
   typedef itk::Image<DisplacementType, Dimension> DeformationFieldType;
-
-
+  
   // Read the input
   typedef itk::ImageFileReader<InputImageType> InputReaderType;
   typename InputReaderType::Pointer reader = InputReaderType::New();
   reader->SetFileName( m_InputFileName);
   reader->Update();
   typename InputImageType::Pointer input= reader->GetOutput();
-
-  //Read the deformation field
-  typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
-  typename  DeformationFieldReaderType::Pointer deformationFieldReader= DeformationFieldReaderType::New();
-  deformationFieldReader->SetFileName(m_ArgsInfo.vf_arg);
-  deformationFieldReader->Update();
-  typename DeformationFieldType::Pointer deformationField =deformationFieldReader->GetOutput();
+  
+  typename DeformationFieldType::Pointer deformationField;
+  if (m_ArgsInfo.coeff_given) {
+    deformationField = CoeffsToDVF<DeformationFieldType>(m_ArgsInfo.coeff_arg, m_InputFileName);
+  }
+  else {
+    //Read the deformation field
+    typedef itk::ImageFileReader<DeformationFieldType> DeformationFieldReaderType;
+    typename  DeformationFieldReaderType::Pointer deformationFieldReader= DeformationFieldReaderType::New();
+    deformationFieldReader->SetFileName(m_ArgsInfo.vf_arg);
+    deformationFieldReader->Update();
+    deformationField =deformationFieldReader->GetOutput();
+  }
 
   // Intensity interpolator
   typedef clitk::GenericVectorInterpolator<args_info_clitkWarpImage, DeformationFieldType, double> GenericInterpolatorType;
@@ -139,7 +208,7 @@ WarpImageGenericFilter::UpdateWithDimAndPixelType()
   // -------------------------------------------
   // Spacing like input
   // -------------------------------------------
-  else {
+  else if (!m_ArgsInfo.coeff_given) {
     // Get size
     typename DeformationFieldType::SizeType newSize;
     for (unsigned int i=0 ; i <Dimension; i++)