Added support to vector field arithmetics

diff --git a/tools/clitkImageArithm.ggo b/tools/clitkImageArithm.ggo
index ca7e91607876aa9d33c608f907177a58513aa453..fd018635cc45b863786141ad0200569e84591d54 100644 (file)
--- a/tools/clitkImageArithm.ggo
+++ b/tools/clitkImageArithm.ggo
@@ -13,7 +13,7 @@ option "input2"          j    "Input second image filename"     string   no
 option "output"    o   "Output image filename"           string   yes
 
 option "scalar"           s    "Scalar value"            double   no
-option "operation" t   "Type of operation : \n With another image : 0=add, 1=multiply, 2=divide,\n 3=max, 4=min, 5=absdiff, 6=squareddiff,  7=difference, 8=relativ diff\n For 'scalar' : 0=add, 1=multiply, 2=inverse,\n 3=max, 4=min 5=absval 6=squareval\n 7=log 8=exp 9=sqrt 10=EPID 11=divide 12=normalize (divide by max)"       int default="0" no 
+option "operation" t   "Type of operation : \n With another image : 0=add*, 1=multiply, 2=divide,\n 3=max, 4=min, 5=absdiff, 6=squareddiff,  7=difference*, 8=relativ diff\n; For 'scalar' : 0=add*, 1=multiply*, 2=inverse,\n 3=max, 4=min 5=absval 6=squareval\n 7=log 8=exp 9=sqrt 10=EPID 11=divide* 12=normalize (divide by max); \n* operations supported with vector fields as inputs."         int default="0" no 
 option "pixelValue" -  "Default value for NaN/Inf"     double   default="0.0"  no
 
 option "setFloatOutput" f "Set output to float pixel type" flag off

diff --git a/tools/clitkImageArithmGenericFilter.cxx b/tools/clitkImageArithmGenericFilter.cxx
index 30b570a0bac5e1b3e162db8ed33095d57e59f3b8..964607eb52a036e361c9cad3520c3502204325cb 100644 (file)
--- a/tools/clitkImageArithmGenericFilter.cxx
+++ b/tools/clitkImageArithmGenericFilter.cxx
@@ -19,12 +19,342 @@
 #define CLITKIMAGEARITHMGENERICFILTER_CXX
 
 #include "clitkImageArithmGenericFilter.h"
-#include "clitkImageArithm_ggo.h"
 
 namespace clitk {
   // Specialisation
+//   template<>
+//   class ImageArithmGenericFilter<args_info_clitkImageArithm>;
+
+  template<>
+  template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >()
+  {
+    typedef itk::Image< itk::Vector<float, 3u>, 3u > ImageType;
+    
+    // Read input1
+    typename ImageType::Pointer input1 = this->GetInput<ImageType>(0);
+
+    // Set input image iterator
+    typedef itk::ImageRegionIterator<ImageType> IteratorType;
+    IteratorType it(input1, input1->GetLargestPossibleRegion());
+
+    // typedef input2
+    typename ImageType::Pointer input2 = NULL;
+    IteratorType it2;
+
+    /*
+    // Special case for normalisation
+    if (mTypeOfOperation == 12) {
+      typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
+      typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
+      ff->SetImage(input1);
+      ff->ComputeMaximum();
+      mScalar = ff->GetMaximum();
+      mTypeOfOperation = 11; // divide
+    }
+    */
+
+    if (mIsOperationUseASecondImage) {
+        // Read input2
+        input2 = this->GetInput<ImageType>(1);
+        // Set input image iterator
+        it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
+        // Check dimension
+        if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
+          itkExceptionMacro(<< "The images (input and input2) must have the same size");
+        }
+        if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
+          itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
+                          << "Using first input's information.");
+        }
+    }
+
+    /*
+    // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
+    if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
+      // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
+      //                     << typeid(PixelType).name()
+      //                     << std::endl;
+      mOverwriteInputImage = false;
+    }
+    */
+    
+    // ---------------- Overwrite input Image ---------------------
+    if (mOverwriteInputImage) {
+      // Set output iterator (to input1)
+      IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
+      if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
+      else ComputeImage(it, ito);
+      this->SetNextOutput<ImageType>(input1);
+    }
+    // ---------------- Create new output Image ---------------------
+    else {
+      /*if (mOutputIsFloat) {
+        // Create output image
+        typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
+        typename OutputImageType::Pointer output = OutputImageType::New();
+        output->SetRegions(input1->GetLargestPossibleRegion());
+        output->SetOrigin(input1->GetOrigin());
+        output->SetSpacing(input1->GetSpacing());
+        output->Allocate();
+        // Set output iterator
+        typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
+        IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
+        if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
+        else ComputeImage(it, ito);
+        this->template SetNextOutput<OutputImageType>(output);
+      } else*/ {
+        // Create output image
+        typedef ImageType OutputImageType;
+        typename OutputImageType::Pointer output = OutputImageType::New();
+        output->SetRegions(input1->GetLargestPossibleRegion());
+        output->SetOrigin(input1->GetOrigin());
+        output->SetSpacing(input1->GetSpacing());
+        output->Allocate();
+        // Set output iterator
+        typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
+        IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
+        if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
+        else ComputeImage(it, ito);
+        this->SetNextOutput<OutputImageType>(output);
+      }
+    }
+  }
+
   template<>
-  class ImageArithmGenericFilter<args_info_clitkImageArithm>;
+  template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage< 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > >
+    (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
+  {
+    typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
+
+    ito.GoToBegin();
+    it.GoToBegin();
+    
+    typedef typename Iter2::PixelType PixelType;
+
+    PixelType scalar_vector;
+    scalar_vector.Fill(mScalar);
+    
+    // Perform operation
+    switch (mTypeOfOperation) {
+    case 0: // Addition
+      while (!it.IsAtEnd()) {
+        ito.Set(it.Get() + scalar_vector);
+        ++it;
+        ++ito;
+      }
+      break;
+    case 1: // Multiply
+      while (!it.IsAtEnd()) {
+        ito.Set(it.Get() * mScalar);
+        ++it;
+        ++ito;
+      }
+      break;
+      /*
+    case 2: // Inverse
+      while (!it.IsAtEnd()) {
+        if (it.Get() != 0)
+          ito.Set(mScalar / it.Get()));
+        else ito.Set(mDefaultPixelValue);
+        ++it;
+        ++ito;
+      }
+      break;
+    case 3: // Max
+      while (!it.IsAtEnd()) {
+        if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
+        else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+        ++it;
+        ++ito;
+      }
+      break;
+    case 4: // Min
+      while (!it.IsAtEnd()) {
+        if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
+        else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+        ++it;
+        ++ito;
+      }
+      break;
+    case 5: // Absolute value
+      while (!it.IsAtEnd()) {
+        if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
+        // <= zero to avoid warning for unsigned types
+        else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+        ++it;
+        ++ito;
+      }
+      break;
+    case 6: // Squared value
+      while (!it.IsAtEnd()) {
+        ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
+        ++it;
+        ++ito;
+      }
+      break;
+    case 7: // Log
+      while (!it.IsAtEnd()) {
+        if (it.Get() > 0)
+          ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
+        else ito.Set(mDefaultPixelValue);
+        ++it;
+        ++ito;
+      }
+      break;
+    case 8: // exp
+      while (!it.IsAtEnd()) {
+        ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
+        ++it;
+        ++ito;
+      }
+      break;
+    case 9: // sqrt
+      while (!it.IsAtEnd()) {
+        if (it.Get() > 0)
+          ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
+        else {
+          if (it.Get() ==0) ito.Set(0);
+          else ito.Set(mDefaultPixelValue);
+        }
+        ++it;
+        ++ito;
+      }
+      break;
+    case 10: // exp
+      while (!it.IsAtEnd()) {
+        ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
+        ++it;
+        ++ito;
+      }
+      break;
+      */
+    case 11: // divide
+      while (!it.IsAtEnd()) {
+        ito.Set(it.Get() / mScalar);
+        ++it;
+        ++ito;
+      }
+      break;
+    default: // error ?
+      std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
+      exit(-1);
+    }
+    
+  }
+
+  template<>
+  template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage< 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > 
+    >
+    (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
+  {
+    typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter1;
+    typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
+    typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter3;
+    
+    it1.GoToBegin();
+    it2.GoToBegin();
+    ito.GoToBegin();
+    typedef typename Iter3::PixelType PixelType;
+
+    switch (mTypeOfOperation) {
+    case 0: // Addition
+      while (!ito.IsAtEnd()) {
+        ito.Set(it1.Get() + it2.Get());
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+      /*
+    case 1: // Multiply
+      while (!ito.IsAtEnd()) {
+        ito.Set(it1.Get() * it2.Get()) );
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+    case 2: // Divide
+      while (!ito.IsAtEnd()) {
+        if (it1.Get() != 0)
+          ito.Set(it1.Get() / it2.Get()));
+        else ito.Set(mDefaultPixelValue);
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+    case 3: // Max
+      while (!ito.IsAtEnd()) {
+        if (it1.Get() < it2.Get()) ito.Set(it2.Get());
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+    case 4: // Min
+      while (!ito.IsAtEnd()) {
+        if (it1.Get() > it2.Get()) ito.Set(it2.Get());
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+      */
+    case 5: // Absolute difference
+      while (!ito.IsAtEnd()) {
+        ito.Set(it2.Get()-it1.Get());
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+      /*
+    case 6: // Squared differences
+      while (!ito.IsAtEnd()) {
+        ito.Set(pow(it1.Get()-it2.Get(),2)));
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+      */
+    case 7: // Difference
+      while (!ito.IsAtEnd()) {
+        ito.Set(it1.Get()-it2.Get());
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+      /*
+    case 8: // Relative Difference
+      while (!ito.IsAtEnd()) {
+        if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
+        else ito.Set(0.0);
+        ++it1;
+        ++it2;
+        ++ito;
+      }
+      break;
+      */
+    default: // error ?
+      std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
+      exit(-1);
+    }
+  }
+
 }
 
 #endif //define CLITKIMAGEARITHMGENERICFILTER_CXX

diff --git a/tools/clitkImageArithmGenericFilter.h b/tools/clitkImageArithmGenericFilter.h
index 46ae08ccaf49c0b67d7448dae360f7244cf226cc..45e0959850775669c859bd7ef108fc72d29c0191 100644 (file)
--- a/tools/clitkImageArithmGenericFilter.h
+++ b/tools/clitkImageArithmGenericFilter.h
@@ -29,6 +29,7 @@
 // clitk include
 #include "clitkCommon.h"
 #include "clitkImageToImageGenericFilter.h"
+#include "clitkImageArithm_ggo.h"
 
 // itk include
 #include "itkImage.h"
@@ -95,9 +96,32 @@ namespace clitk {
     //--------------------------------------------------------------------
 
   }; // end class ImageArithmGenericFilter
+
+  template<> template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >();
+  
+  template<> template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage< 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > 
+    >
+    (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it, 
+     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito);
+
+  template<> template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage< 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > 
+    >
+    (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1, 
+     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2, 
+     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito);
+  
 } // end namespace
 //--------------------------------------------------------------------
 
+  
 #ifndef ITK_MANUAL_INSTANTIATION
 #include "clitkImageArithmGenericFilter.txx"
 #endif

diff --git a/tools/clitkImageArithmGenericFilter.txx b/tools/clitkImageArithmGenericFilter.txx
index 24ebbbdac8d2d472ec6ddca3eebbdcfdd5f33f6f..8efd6a7dfcabefe54f2e750b7f58ea91637ed8a3 100644 (file)
--- a/tools/clitkImageArithmGenericFilter.txx
+++ b/tools/clitkImageArithmGenericFilter.txx
@@ -45,6 +45,7 @@ template<unsigned int Dim>
 void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
 {
   ADD_DEFAULT_IMAGE_TYPES(Dim);
+  ADD_VEC_IMAGE_TYPE(3u,3u,float);
 }
 //--------------------------------------------------------------------
 
@@ -398,6 +399,8 @@ void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Ite
 }
 //--------------------------------------------------------------------
 
+
+
 } // end namespace
 
 #endif  //#define CLITKIMAGEARITHMGENERICFILTER_TXX