X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkAffineTransformGenericFilter.txx;h=bf31033d515b9df7ee321e1a1d1bce1fc0ae5829;hb=d1dc57c5af88ff53c75904b1dd55722f9631a21e;hp=b2d7a534aa7e250b2d51a4300ebc144fe88669dd;hpb=26bcb2de1cbf158bdb6b76175a7c8db4b5f269dc;p=clitk.git

diff --git a/tools/clitkAffineTransformGenericFilter.txx b/tools/clitkAffineTransformGenericFilter.txx
index b2d7a53..bf31033 100644
--- a/tools/clitkAffineTransformGenericFilter.txx
+++ b/tools/clitkAffineTransformGenericFilter.txx
@@ -22,6 +22,8 @@
 #include <istream>
 #include <iterator>
 #include <itkCenteredEuler3DTransform.h>
+#include <itkRecursiveGaussianImageFilter.h>
+#include "clitkElastix.h"
 
 namespace clitk
 {
@@ -52,12 +54,12 @@ namespace clitk
     // Call UpdateWithDim
     if(Dimension==2) UpdateWithDim<2>(PixelType, Components);
     else 
-    if(Dimension==3) UpdateWithDim<3>(PixelType, Components);
-    else if (Dimension==4)UpdateWithDim<4>(PixelType, Components);
-    else {
-      std::cout<<"Error, Only for 2, 3 or 4  Dimensions!!!"<<std::endl ;
-      return;
-    }
+      if(Dimension==3) UpdateWithDim<3>(PixelType, Components);
+      else if (Dimension==4)UpdateWithDim<4>(PixelType, Components);
+      else {
+        std::cout<<"Error, Only for 2, 3 or 4  Dimensions!!!"<<std::endl ;
+        return;
+      }
   }
   //-------------------------------------------------------------------
  
@@ -77,10 +79,10 @@ namespace clitk
         if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed short..." << std::endl;
         UpdateWithDimAndPixelType<Dimension, signed short>();
       }
-      //    else if(PixelType == "unsigned_short"){
-      //       if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
-      //       UpdateWithDimAndPixelType<Dimension, unsigned short>();
-      //     }
+      else if(PixelType == "unsigned_short"){
+        if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl;
+        UpdateWithDimAndPixelType<Dimension, unsigned short>();
+      }
 
       else if (PixelType == "unsigned_char") {
         if (m_Verbose) std::cout  << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl;
@@ -128,6 +130,37 @@ 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
     typedef  itk::ResampleImageFilter< InputImageType,OutputImageType >  ResampleFilterType;
     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
@@ -138,7 +171,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);
@@ -181,7 +214,8 @@ namespace clitk
           }
         else {
           if (m_ArgsInfo.elastix_given) {
-            matrix = createMatrixFromElastixFile<Dimension,PixelType>(m_ArgsInfo.elastix_arg);
+            std::string filename(m_ArgsInfo.elastix_arg);
+            matrix = createMatrixFromElastixFile<Dimension>(filename, m_Verbose);
           }
           else 
             matrix.SetIdentity();
@@ -210,6 +244,15 @@ namespace clitk
       likeReader->SetFileName(m_ArgsInfo.like_arg);
       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) );
@@ -226,6 +269,14 @@ namespace clitk
       outputSpacing = invRotMatrix *
         input->GetDirection() *
         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 is influenced by translation but not by input direction
       typename InputImageType::PointType outputOrigin;
@@ -270,6 +321,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;
@@ -278,10 +337,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 );
 
     }
 
@@ -312,9 +380,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) );
@@ -356,6 +446,37 @@ 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
     typedef  itk::VectorResampleImageFilter< InputImageType,OutputImageType, double >  ResampleFilterType;
     typename ResampleFilterType::Pointer resampler = ResampleFilterType::New();
@@ -366,7 +487,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);
@@ -435,6 +556,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;
@@ -450,6 +580,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
@@ -460,13 +598,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());
@@ -489,121 +659,6 @@ namespace clitk
 
   }
   //-------------------------------------------------------------------
-  
-  
-  //-------------------------------------------------------------------
-  template<class args_info_type>
-  template<unsigned int Dimension, class PixelType>
-  typename itk::Matrix<double, Dimension+1, Dimension+1>
-   AffineTransformGenericFilter<args_info_type>::createMatrixFromElastixFile(std::string filename)
-  {
-    if (Dimension != 3) {
-      FATAL("Only 3D yet" << std::endl);
-    }
-    typename itk::Matrix<double, Dimension+1, Dimension+1> matrix;
-
-    // Open file
-    std::ifstream is;
-    clitk::openFileForReading(is, filename);
-
-    // Check Transform
-    std::string s; 
-    bool b = GetElastixValueFromTag(is, "Transform ", s);
-    if (!b) {
-      FATAL("Error must read 'Transform' in " << filename << std::endl);
-    }
-    if (s != "EulerTransform") {
-      FATAL("Sorry only 'EulerTransform'" << std::endl);
-    }
-
-    // FIXME check
-    //    (InitialTransformParametersFileName "NoInitialTransform")
-
-    // Get CenterOfRotationPoint
-    GetElastixValueFromTag(is, "CenterOfRotationPoint ", s); // space is needed
-    if (!b) {
-      FATAL("Error must read 'CenterOfRotationPoint' in " << filename << std::endl);
-    }
-    std::vector<std::string> cor; 
-    GetValuesFromValue(s, cor);
-
-    // Get Transformparameters
-    GetElastixValueFromTag(is, "TransformParameters ", s); // space is needed
-    if (!b) {
-      FATAL("Error must read 'TransformParameters' in " << filename << std::endl);
-    }
-    std::vector<std::string> results; 
-    GetValuesFromValue(s, results);
-    
-    // construct a stream from the string
-    itk::CenteredEuler3DTransform<double>::Pointer mat = itk::CenteredEuler3DTransform<double>::New();
-    itk::CenteredEuler3DTransform<double>::ParametersType p;
-    p.SetSize(9);
-    for(uint i=0; i<3; i++)
-      p[i] = atof(results[i].c_str()); // Rotation
-    for(uint i=0; i<3; i++)
-      p[i+3] = atof(cor[i].c_str()); // Centre of rotation
-    for(uint i=0; i<3; i++)
-      p[i+6] = atof(results[i+3].c_str()); // Translation
-    mat->SetParameters(p);
-    
-    if (m_Verbose) {
-      std::cout << "Rotation      (deg) : " << rad2deg(p[0]) << " " << rad2deg(p[1]) << " " << rad2deg(p[2]) << std::endl;
-      std::cout << "Translation   (phy) : " << p[3] << " " << p[4] << " " << p[5] << std::endl;
-      std::cout << "Center of rot (phy) : " << p[6] << " " << p[7] << " " << p[8] << std::endl;
-    }
-
-    for(uint i=0; i<3; i++)
-      for(uint j=0; j<3; j++)
-        matrix[i][j] = mat->GetMatrix()[i][j];
-    // Offset is -Rc + t + c
-    matrix[0][3] = mat->GetOffset()[0];
-    matrix[1][3] = mat->GetOffset()[1];
-    matrix[2][3] = mat->GetOffset()[2];
-    matrix[3][3] = 1;
-    
-    return matrix;
-  }
-
-  //-------------------------------------------------------------------
-  template<class args_info_type>
-  bool
-  AffineTransformGenericFilter<args_info_type>::GetElastixValueFromTag(std::ifstream & is, 
-                                                                       std::string tag, 
-                                                                       std::string & value)
-  {
-    std::string line;
-    is.seekg (0, is.beg);
-    while(std::getline(is, line))   {
-      unsigned pos = line.find(tag);
-      if (pos<line.size()) {
-        value=line.substr(pos+tag.size(),line.size()-2);// remove last ')'
-        value.erase (std::remove (value.begin(), value.end(), '"'), value.end());
-        value.erase (std::remove (value.begin(), value.end(), ')'), value.end());
-        return true;
-      }
-   }
-    return false;
-  }
-  //-------------------------------------------------------------------
-
-
-  //-------------------------------------------------------------------
-  template<class args_info_type>
-  void
-  AffineTransformGenericFilter<args_info_type>::GetValuesFromValue(const std::string & s, 
-                                                                   std::vector<std::string> & values)
-  {
-    std::stringstream strstr(s);
-    std::istream_iterator<std::string> it(strstr);
-    std::istream_iterator<std::string> end;
-    std::vector<std::string> results(it, end);
-    values.clear();
-    values.resize(results.size());
-    for(uint i=0; i<results.size(); i++) values[i] = results[i];
-  }
-  //-------------------------------------------------------------------
-
 
 } //end clitk