X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkAffineTransformGenericFilter.txx;h=c1fe88a40f4d06d6a94ba34983cc3f958cad1ade;hb=998ca9ca987553dc71578e5584e881f6bd537d91;hp=65d66dc741c6256637eefbeb5c222a25860caf87;hpb=784796d8d6a038d3afd29cea1bf4e73323852186;p=clitk.git

diff --git a/tools/clitkAffineTransformGenericFilter.txx b/tools/clitkAffineTransformGenericFilter.txx
index 65d66dc..c1fe88a 100644
--- a/tools/clitkAffineTransformGenericFilter.txx
+++ b/tools/clitkAffineTransformGenericFilter.txx
@@ -22,7 +22,9 @@
 #include <istream>
 #include <iterator>
 #include <itkCenteredEuler3DTransform.h>
+#include <itkRecursiveGaussianImageFilter.h>
 #include "clitkElastix.h"
+#include "clitkResampleImageWithOptionsFilter.h"
 
 namespace clitk
 {
@@ -92,6 +94,10 @@ namespace clitk
       //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl;
       //       UpdateWithDimAndPixelType<Dimension, signed char>();
       //     }
+      else if(PixelType == "double"){
+        if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and double..." << std::endl;
+        UpdateWithDimAndPixelType<Dimension, double>();
+      }
       else {
         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and float..." << std::endl;
         UpdateWithDimAndPixelType<Dimension, float>();
@@ -107,6 +113,95 @@ namespace clitk
 
   }
   //-------------------------------------------------------------------
+
+
+  //-------------------------------------------------------------------
+  // Compute updated bounding box
+  //-------------------------------------------------------------------
+  template<class args_info_type>
+  vnl_vector<double>
+  AffineTransformGenericFilter<args_info_type>::ComputeSize(vnl_vector<double> inputSize, vnl_matrix<double> transformationMatrix, bool returnMin)
+  {
+    //Compute input corners
+    int Dimension = inputSize.size();
+    vnl_matrix<double> vnlOutputSize(std::pow(2, Dimension), Dimension);
+    vnlOutputSize.fill(0);
+    if (Dimension == 2) {
+      for(unsigned int i=0; i< Dimension; i++)
+        vnlOutputSize[3][i] = inputSize[i];
+      vnlOutputSize[1][0] = inputSize[0];
+      vnlOutputSize[2][1] = inputSize[1];
+    } else if (Dimension == 3) {
+      for(unsigned int i=0; i< Dimension; i++)
+        vnlOutputSize[7][i] = inputSize[i];
+      vnlOutputSize[1][0] = inputSize[0];
+      vnlOutputSize[2][1] = inputSize[1];
+      vnlOutputSize[3][2] = inputSize[2];
+      vnlOutputSize[4][0] = inputSize[0];
+      vnlOutputSize[4][1] = inputSize[1];
+      vnlOutputSize[5][1] = inputSize[1];
+      vnlOutputSize[5][2] = inputSize[2];
+      vnlOutputSize[6][0] = inputSize[0];
+      vnlOutputSize[6][2] = inputSize[2];
+    } else { //Dimension ==4
+      for(unsigned int i=0; i< Dimension; i++)
+        vnlOutputSize[15][i] = inputSize[i];
+      vnlOutputSize[1][0] = inputSize[0];
+      vnlOutputSize[2][1] = inputSize[1];
+      vnlOutputSize[3][2] = inputSize[2];
+      vnlOutputSize[4][3] = inputSize[3];
+      vnlOutputSize[5][0] = inputSize[0];
+      vnlOutputSize[5][1] = inputSize[1];
+      vnlOutputSize[6][0] = inputSize[0];
+      vnlOutputSize[6][2] = inputSize[2];
+      vnlOutputSize[7][0] = inputSize[0];
+      vnlOutputSize[7][3] = inputSize[3];
+      vnlOutputSize[8][1] = inputSize[1];
+      vnlOutputSize[8][2] = inputSize[2];
+      vnlOutputSize[9][1] = inputSize[1];
+      vnlOutputSize[9][3] = inputSize[3];
+      vnlOutputSize[10][2] = inputSize[2];
+      vnlOutputSize[10][3] = inputSize[3];
+      vnlOutputSize[11][0] = inputSize[0];
+      vnlOutputSize[11][1] = inputSize[1];
+      vnlOutputSize[11][2] = inputSize[2];
+      vnlOutputSize[12][0] = inputSize[0];
+      vnlOutputSize[12][1] = inputSize[1];
+      vnlOutputSize[12][3] = inputSize[3];
+      vnlOutputSize[13][0] = inputSize[0];
+      vnlOutputSize[13][2] = inputSize[2];
+      vnlOutputSize[13][3] = inputSize[3];
+      vnlOutputSize[14][1] = inputSize[1];
+      vnlOutputSize[14][2] = inputSize[2];
+      vnlOutputSize[14][3] = inputSize[3];
+    }
+
+    //Compute the transformation of all corner
+    for (unsigned int i=0; i< std::pow(2, Dimension); ++i)
+      vnlOutputSize.set_row(i, transformationMatrix*vnlOutputSize.get_row(i));
+
+    //Compute the bounding box taking the max and the min
+    vnl_vector<double> minBB(vnlOutputSize.get_row(0)), maxBB(vnlOutputSize.get_row(0));
+    for (unsigned int i=0; i< std::pow(2, Dimension); ++i) {
+      for (unsigned int j=0; j< Dimension; ++j) {
+        if (vnlOutputSize[i][j] < minBB[j])
+          minBB[j] = vnlOutputSize[i][j];
+        if (vnlOutputSize[i][j] > maxBB[j])
+          maxBB[j] = vnlOutputSize[i][j];
+      }
+    }
+
+    //Compute the size
+    if (returnMin)
+      return minBB;
+    else {
+      vnl_vector<double> size;
+      size = maxBB - minBB;
+
+      return size;
+    }
+  }
+  //-------------------------------------------------------------------
  
 
   //-------------------------------------------------------------------
@@ -129,6 +224,173 @@ namespace clitk
     reader->Update();
     typename InputImageType::Pointer input= reader->GetOutput();
 
+    //Adaptative size, spacing origin (use previous clitkResampleImage)
+    if (m_ArgsInfo.adaptive_given) {
+      // Filter
+      typedef clitk::ResampleImageWithOptionsFilter<InputImageType, OutputImageType> ResampleImageFilterType;
+      typename ResampleImageFilterType::Pointer filter = ResampleImageFilterType::New();
+      filter->SetInput(input);
+
+      // Set Verbose
+      filter->SetVerboseOptions(m_ArgsInfo.verbose_flag);
+
+      // Set size / spacing
+      static const unsigned int dim = OutputImageType::ImageDimension;
+      typename OutputImageType::SpacingType spacing;
+      typename OutputImageType::SizeType size;
+      typename OutputImageType::PointType origin;
+      typename OutputImageType::DirectionType direction;
+
+      if (m_ArgsInfo.like_given) {
+        itk::ImageIOBase::Pointer header = clitk::readImageHeader(m_ArgsInfo.like_arg);
+        if (header) {
+          for(unsigned int i=0; i<dim; i++){
+            spacing[i] = header->GetSpacing(i);
+            size[i] = header->GetDimensions(i);
+            origin[i] = header->GetOrigin(i);
+          }
+          for(unsigned int i=0; i<dim; i++) {
+            for(unsigned int j=0;j<dim;j++) {
+                direction(i,j) = header->GetDirection(i)[j];
+            }
+          }
+          filter->SetOutputSpacing(spacing);
+          filter->SetOutputSize(size);
+          filter->SetOutputOrigin(origin);
+          filter->SetOutputDirection(direction);
+        }
+        else {
+          std::cerr << "*** Warning : I could not read '" << m_ArgsInfo.like_arg << "' ***" << std::endl;
+          exit(0);
+        }
+      }
+      else {
+        if (m_ArgsInfo.spacing_given == 1) {
+          filter->SetOutputIsoSpacing(m_ArgsInfo.spacing_arg[0]);
+        }
+        else if ((m_ArgsInfo.spacing_given != 0) && (m_ArgsInfo.size_given != 0)) {
+          std::cerr << "Error: use spacing or size, not both." << std::endl;
+          exit(0);
+        }
+        else if (m_ArgsInfo.spacing_given) {
+          if ((m_ArgsInfo.spacing_given != 0) && (m_ArgsInfo.spacing_given != dim)) {
+            std::cerr << "Error: spacing should have one or " << dim << " values." << std::endl;
+            exit(0);
+          }
+          for(unsigned int i=0; i<dim; i++)
+            spacing[i] = m_ArgsInfo.spacing_arg[i];
+          filter->SetOutputSpacing(spacing);
+        }
+        else if (m_ArgsInfo.size_given) {
+          if ((m_ArgsInfo.size_given != 0) && (m_ArgsInfo.size_given != dim)) {
+            std::cerr << "Error: size should have " << dim << " values." << std::endl;
+            exit(0);
+          }
+          for(unsigned int i=0; i<dim; i++)
+            size[i] = m_ArgsInfo.size_arg[i];
+          filter->SetOutputSize(size);
+        }
+        for(unsigned int i=0; i<dim; i++){
+          origin[i] = input->GetOrigin()[i];
+        }
+        for(unsigned int i=0; i<dim; i++) {
+          for(unsigned int j=0;j<dim;j++) {
+              direction(i,j) = input->GetDirection()[i][j];
+          }
+        }
+        filter->SetOutputOrigin(origin);
+        filter->SetOutputDirection(direction);
+      }
+
+      // Set temporal dimension
+      //filter->SetLastDimensionIsTime(m_ArgsInfo.time_flag);
+
+      // Set Gauss
+      filter->SetGaussianFilteringEnabled(m_ArgsInfo.autogauss_flag);
+      if (m_ArgsInfo.gauss_given != 0) {
+        typename ResampleImageFilterType::GaussianSigmaType g;
+        for(unsigned int i=0; i<dim; i++) {
+          g[i] = m_ArgsInfo.gauss_arg[i];
+        }
+        filter->SetGaussianSigma(g);
+      }
+
+      // Set Interpolation
+      int interp = m_ArgsInfo.interp_arg;
+      if (interp == 0) {
+        filter->SetInterpolationType(ResampleImageFilterType::NearestNeighbor);
+      } else {
+        if (interp == 1) {
+          filter->SetInterpolationType(ResampleImageFilterType::Linear);
+        } else {
+          if (interp == 2) {
+            filter->SetInterpolationType(ResampleImageFilterType::BSpline);
+          } else {
+            if (interp == 3) {
+              filter->SetInterpolationType(ResampleImageFilterType::B_LUT);
+            } else {
+                std::cerr << "Error. I do not know interpolation '" << m_ArgsInfo.interp_arg
+                          << "'. Choose among: nn, linear, bspline, blut, windowed sinc" << std::endl;
+                exit(0);
+            }
+          }
+        }
+      }
+
+      // Set default pixel value
+      filter->SetDefaultPixelValue(m_ArgsInfo.pad_arg);
+
+      // Set thread
+      //if (m_ArgsInfo.thread_given) {
+      //  filter->SetNumberOfThreads(m_ArgsInfo.thread_arg);
+      //}
+
+      // Go !
+      filter->Update();
+      typename OutputImageType::Pointer output = filter->GetOutput();
+      //this->template SetNextOutput<OutputImageType>(outputImage);
+
+      // Output
+      typedef itk::ImageFileWriter<OutputImageType> WriterType;
+      typename WriterType::Pointer writer = WriterType::New();
+      writer->SetFileName(m_ArgsInfo.output_arg);
+      writer->SetInput(output);
+      writer->Update();
+
+      return;
+    }
+
+    //Gaussian pre-filtering
+    typename itk::Vector<double, Dimension> gaussianSigma;
+    gaussianSigma.Fill(0);
+    bool gaussianFilteringEnabled(false);
+    bool autoGaussEnabled(false);
+    if (m_ArgsInfo.autogauss_given) { // Gaussian filter auto
+      autoGaussEnabled = m_ArgsInfo.autogauss_flag;
+    }
+    if (m_ArgsInfo.gauss_given) { // Gaussian filter set by user
+      gaussianFilteringEnabled = true;
+      if (m_ArgsInfo.gauss_given == 1)
+      {
+        for (unsigned int i=0; i<Dimension; i++)
+        {
+          gaussianSigma[i] = m_ArgsInfo.gauss_arg[0];
+        }
+      }
+      else if (m_ArgsInfo.gauss_given == Dimension)
+      {
+        for (unsigned int i=0; i<Dimension; i++)
+        {
+          gaussianSigma[i] = m_ArgsInfo.gauss_arg[i];
+        }
+      }
+      else
+      {
+        std::cerr << "Gaussian sigma dimension is incorrect" << std::endl;
+        return;
+      }
+    }
+
     //Filter
     typedef  itk::ResampleImageFilter< InputImageType,OutputImageType >  ResampleFilterType;
     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
@@ -139,7 +401,7 @@ namespace clitk
       {
         if (m_ArgsInfo.matrix_given)
           {
-            std::cerr << "You must use either rotate/translate or matrix options" << std::cout;
+            std::cerr << "You must use either rotate/translate or matrix options" << std::endl;
             return;
           }
         itk::Array<double> transformParameters(2 * Dimension);
@@ -213,6 +475,14 @@ namespace clitk
       likeReader->Update();
       resampler->SetOutputParametersFromImage(likeReader->GetOutput());
       resampler->SetOutputDirection(likeReader->GetOutput()->GetDirection());
+      if (autoGaussEnabled) { // Automated sigma when downsample
+        for(unsigned int i=0; i<Dimension; i++) {
+          if (likeReader->GetOutput()->GetSpacing()[i] > input->GetSpacing()[i]) { // downsample
+            gaussianSigma[i] = 0.5*likeReader->GetOutput()->GetSpacing()[i];// / inputSpacing[i]);
+          }
+          else gaussianSigma[i] = 0; // will be ignore after
+        }
+      }
     } else if(m_ArgsInfo.transform_grid_flag) {
       typename itk::Matrix<double, Dimension+1, Dimension+1> invMatrix( matrix.GetInverse() );
       typename itk::Matrix<double, Dimension, Dimension> invRotMatrix( clitk::GetRotationalPartMatrix(invMatrix) );
@@ -224,12 +494,6 @@ namespace clitk
       if (m_ArgsInfo.origin_given)
         std::cout << "Warning --origin ignored (because --transform_grid_flag)" << std::endl;
 
-      // Spacing is influenced by affine transform matrix and input direction
-      typename InputImageType::SpacingType outputSpacing;
-      outputSpacing = invRotMatrix *
-        input->GetDirection() *
-        input->GetSpacing();
-
       // Origin is influenced by translation but not by input direction
       typename InputImageType::PointType outputOrigin;
       outputOrigin = invRotMatrix *
@@ -238,13 +502,30 @@ namespace clitk
 
       // Size is influenced by affine transform matrix and input direction
       // Size is converted to double, transformed and converted back to size type.
-      vnl_vector<double> vnlOutputSize(Dimension);
+      // Determine the bounding box tranforming all corners
+      vnl_vector<double> vnlOutputSize(Dimension), vnlOutputmmSize(Dimension), vnlOutputOffset(Dimension);
+      typename InputImageType::SpacingType outputSpacing;
       for(unsigned int i=0; i< Dimension; i++) {
         vnlOutputSize[i] = input->GetLargestPossibleRegion().GetSize()[i];
+        vnlOutputmmSize[i] = input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i];
+        vnlOutputOffset[i] = input->GetLargestPossibleRegion().GetSize()[i]*input->GetSpacing()[i];
+      }
+      vnlOutputSize = ComputeSize(vnlOutputSize, invRotMatrix.GetVnlMatrix() * input->GetDirection().GetVnlMatrix(), 0);
+      vnlOutputmmSize = ComputeSize(vnlOutputmmSize, invRotMatrix.GetVnlMatrix() * input->GetDirection().GetVnlMatrix(), 0);
+      vnlOutputOffset = ComputeSize(vnlOutputOffset, invRotMatrix.GetVnlMatrix() * input->GetDirection().GetVnlMatrix(), 1);
+      for(unsigned int i=0; i< Dimension; i++) {
+        outputSpacing[i] = vnlOutputmmSize[i]/lrint(vnlOutputSize[i]);
+        outputOrigin[i] += vnlOutputOffset[i];
+      }
+      if (autoGaussEnabled) { // Automated sigma when downsample
+        for(unsigned int i=0; i<Dimension; i++) {
+          if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
+            gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
+          }
+          else gaussianSigma[i] = 0; // will be ignore after
+        }
       }
-      vnlOutputSize = invRotMatrix *
-        input->GetDirection().GetVnlMatrix() *
-        vnlOutputSize;
+
       typename OutputImageType::SizeType outputSize;
       for(unsigned int i=0; i< Dimension; i++) {
         // If the size is negative, we have a flip and we must modify
@@ -273,6 +554,14 @@ namespace clitk
         for(unsigned int i=0; i< Dimension; i++)
           outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
       } else outputSpacing=input->GetSpacing();
+      if (autoGaussEnabled) { // Automated sigma when downsample
+        for(unsigned int i=0; i<Dimension; i++) {
+          if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
+            gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
+          }
+          else gaussianSigma[i] = 0; // will be ignore after
+        }
+      }
 
       //Origin
       typename OutputImageType::PointType outputOrigin;
@@ -281,10 +570,19 @@ namespace clitk
           outputOrigin[i]=m_ArgsInfo.origin_arg[i];
       } else outputOrigin=input->GetOrigin();
 
+      //Direction
+      typename OutputImageType::DirectionType outputDirection;
+      if (m_ArgsInfo.direction_given) {
+        for(unsigned int j=0; j< Dimension; j++)
+            for(unsigned int i=0; i< Dimension; i++)
+                outputDirection[j][i]=m_ArgsInfo.direction_arg[i+Dimension*j];
+      } else outputDirection=input->GetDirection();
+
       // Set
       resampler->SetSize( outputSize );
       resampler->SetOutputSpacing( outputSpacing );
       resampler->SetOutputOrigin(  outputOrigin );
+      resampler->SetOutputDirection( outputDirection );
 
     }
 
@@ -315,9 +613,31 @@ namespace clitk
       std::cout << "Setting the output size to " << resampler->GetSize() << "..." << std::endl;
       std::cout << "Setting the output spacing to " << resampler->GetOutputSpacing() << "..." << std::endl;
       std::cout << "Setting the output origin to " << resampler->GetOutputOrigin() << "..." << std::endl;
+      std::cout << "Setting the output direction to " << resampler->GetOutputDirection() << "..." << std::endl;
     }
 
-    resampler->SetInput( input );
+    typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
+    std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
+    if (gaussianFilteringEnabled || autoGaussEnabled) {
+      for(unsigned int i=0; i<Dimension; i++) {
+        if (gaussianSigma[i] != 0) {
+          gaussianFilters.push_back(GaussianFilterType::New());
+          gaussianFilters[i]->SetDirection(i);
+          gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
+          gaussianFilters[i]->SetNormalizeAcrossScale(false);
+          gaussianFilters[i]->SetSigma(gaussianSigma[i]); // in millimeter !
+          if (gaussianFilters.size() == 1) { // first
+            gaussianFilters[0]->SetInput(input);
+          } else {
+            gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
+          }
+        }
+      }
+      if (gaussianFilters.size() > 0) {
+        resampler->SetInput(gaussianFilters[gaussianFilters.size()-1]->GetOutput());
+      } else resampler->SetInput(input);
+    } else resampler->SetInput(input);
+
     resampler->SetTransform( affineTransform );
     resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());
     resampler->SetDefaultPixelValue( static_cast<PixelType>(m_ArgsInfo.pad_arg) );
@@ -359,8 +679,43 @@ namespace clitk
     reader->Update();
     typename InputImageType::Pointer input= reader->GetOutput();
 
+    //Gaussian pre-filtering
+    typename itk::Vector<double, Dimension> gaussianSigma;
+    gaussianSigma.Fill(0);
+    bool gaussianFilteringEnabled(false);
+    bool autoGaussEnabled(false);
+    if (m_ArgsInfo.autogauss_given) { // Gaussian filter auto
+      autoGaussEnabled = m_ArgsInfo.autogauss_flag;
+    }
+    if (m_ArgsInfo.gauss_given) { // Gaussian filter set by user
+      gaussianFilteringEnabled = true;
+      if (m_ArgsInfo.gauss_given == 1)
+      {
+        for (unsigned int i=0; i<Dimension; i++)
+        {
+          gaussianSigma[i] = m_ArgsInfo.gauss_arg[0];
+        }
+      }
+      else if (m_ArgsInfo.gauss_given == Dimension)
+      {
+        for (unsigned int i=0; i<Dimension; i++)
+        {
+          gaussianSigma[i] = m_ArgsInfo.gauss_arg[i];
+        }
+      }
+      else
+      {
+        std::cerr << "Gaussian sigma dimension is incorrect" << std::endl;
+        return;
+      }
+    }
+
     //Filter
+#if ( ITK_VERSION_MAJOR < 5 )
     typedef  itk::VectorResampleImageFilter< InputImageType,OutputImageType, double >  ResampleFilterType;
+#else
+    typedef  itk::ResampleImageFilter< InputImageType,OutputImageType, double >  ResampleFilterType;
+#endif
     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
 
     // Matrix
@@ -369,7 +724,7 @@ namespace clitk
       {
         if (m_ArgsInfo.matrix_given)
           {
-            std::cerr << "You must use either rotate/translate or matrix options" << std::cout;
+            std::cerr << "You must use either rotate/translate or matrix options" << std::endl;
             return;
           }
         itk::Array<double> transformParameters(2 * Dimension);
@@ -438,6 +793,15 @@ namespace clitk
       resampler->SetSize( likeReader->GetOutput()->GetLargestPossibleRegion().GetSize() );
       resampler->SetOutputSpacing( likeReader->GetOutput()->GetSpacing() );
       resampler->SetOutputOrigin(  likeReader->GetOutput()->GetOrigin() );
+      resampler->SetOutputDirection( likeReader->GetOutput()->GetDirection() );
+      if (autoGaussEnabled) { // Automated sigma when downsample
+        for(unsigned int i=0; i<Dimension; i++) {
+          if (likeReader->GetOutput()->GetSpacing()[i] > input->GetSpacing()[i]) { // downsample
+            gaussianSigma[i] = 0.5*likeReader->GetOutput()->GetSpacing()[i];// / inputSpacing[i]);
+          }
+          else gaussianSigma[i] = 0; // will be ignore after
+        }
+      }
     } else {
       //Size
       typename OutputImageType::SizeType outputSize;
@@ -453,6 +817,14 @@ namespace clitk
         for(unsigned int i=0; i< Dimension; i++)
           outputSpacing[i]=m_ArgsInfo.spacing_arg[i];
       } else outputSpacing=input->GetSpacing();
+      if (autoGaussEnabled) { // Automated sigma when downsample
+        for(unsigned int i=0; i<Dimension; i++) {
+          if (outputSpacing[i] > input->GetSpacing()[i]) { // downsample
+            gaussianSigma[i] = 0.5*outputSpacing[i];// / inputSpacing[i]);
+          }
+          else gaussianSigma[i] = 0; // will be ignore after
+        }
+      }
       std::cout<<"Setting the spacing to "<<outputSpacing<<"..."<<std::endl;
 
       //Origin
@@ -463,13 +835,45 @@ namespace clitk
       } else outputOrigin=input->GetOrigin();
       std::cout<<"Setting the origin to "<<outputOrigin<<"..."<<std::endl;
 
+      //Direction
+      typename OutputImageType::DirectionType outputDirection;
+      if (m_ArgsInfo.direction_given) {
+        for(unsigned int j=0; j< Dimension; j++)
+            for(unsigned int i=0; i< Dimension; i++)
+                outputDirection[j][i]=m_ArgsInfo.direction_arg[i+Dimension*j];
+      } else outputDirection=input->GetDirection();
+      std::cout<<"Setting the direction to "<<outputDirection<<"..."<<std::endl;
+
       // Set
       resampler->SetSize( outputSize );
       resampler->SetOutputSpacing( outputSpacing );
       resampler->SetOutputOrigin(  outputOrigin );
+      resampler->SetOutputDirection( outputDirection );
 
     }
 
+    typedef itk::RecursiveGaussianImageFilter<InputImageType, InputImageType> GaussianFilterType;
+    std::vector<typename GaussianFilterType::Pointer> gaussianFilters;
+    if (gaussianFilteringEnabled || autoGaussEnabled) {
+      for(unsigned int i=0; i<Dimension; i++) {
+        if (gaussianSigma[i] != 0) {
+          gaussianFilters.push_back(GaussianFilterType::New());
+          gaussianFilters[i]->SetDirection(i);
+          gaussianFilters[i]->SetOrder(GaussianFilterType::ZeroOrder);
+          gaussianFilters[i]->SetNormalizeAcrossScale(false);
+          gaussianFilters[i]->SetSigma(gaussianSigma[i]); // in millimeter !
+          if (gaussianFilters.size() == 1) { // first
+            gaussianFilters[0]->SetInput(input);
+          } else {
+            gaussianFilters[i]->SetInput(gaussianFilters[i-1]->GetOutput());
+          }
+        }
+      }
+      if (gaussianFilters.size() > 0) {
+        resampler->SetInput(gaussianFilters[gaussianFilters.size()-1]->GetOutput());
+      } else resampler->SetInput(input);
+    } else resampler->SetInput(input);
+
     resampler->SetInput( input );
     resampler->SetTransform( affineTransform );
     resampler->SetInterpolator( genericInterpolator->GetInterpolatorPointer());