From: Romulo Pinho <romulo.pinho@lyon.unicancer.fr>
Date: Wed, 22 Feb 2012 09:48:26 +0000 (+0100)
Subject: added support to new vector type arithmetics
X-Git-Tag: v1.3.0~89
X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=a0209d730acea6ba9f39a93e04632353f7f1401f;p=clitk.git

added support to new vector type arithmetics
---

diff --git a/tools/clitkImageArithmGenericFilter.cxx b/tools/clitkImageArithmGenericFilter.cxx
index 71da90e..f0d55c2 100644
--- a/tools/clitkImageArithmGenericFilter.cxx
+++ b/tools/clitkImageArithmGenericFilter.cxx
@@ -25,6 +25,8 @@ namespace clitk {
 //   template<>
 //   class ImageArithmGenericFilter<args_info_clitkImageArithm>;
 
+////////////////////////////////////////////////////////////////////
+  // specializations for itk::Vector<float, 3u>, 3u
   template<>
   template<>
   void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >()
@@ -355,6 +357,338 @@ namespace clitk {
     }
   }
 
+////////////////////////////////////////////////////////////////////
+  // specializations for itk::Vector<double, 3u>, 3u
+  template<>
+  template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<double, 3u>, 3u > >()
+  {
+    typedef itk::Image< itk::Vector<double, 3u>, 3u > ImageType;
+    
+    // Read input1
+    ImageType::Pointer input1 = this->GetInput<ImageType>(0);
+
+    // Set input image iterator
+    typedef itk::ImageRegionIterator<ImageType> IteratorType;
+    IteratorType it(input1, input1->GetLargestPossibleRegion());
+
+    // typedef input2
+    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;
+        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<>
+  template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage< 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > >
+    (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
+  {
+    typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter2;
+
+    ito.GoToBegin();
+    it.GoToBegin();
+    
+    typedef 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<double, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > 
+    >
+    (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it1, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it2, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
+  {
+    typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter1;
+    typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter2;
+    typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter3;
+    
+    it1.GoToBegin();
+    it2.GoToBegin();
+    ito.GoToBegin();
+    typedef 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 45e0959..39bfeb9 100644
--- a/tools/clitkImageArithmGenericFilter.h
+++ b/tools/clitkImageArithmGenericFilter.h
@@ -97,6 +97,7 @@ namespace clitk {
 
   }; // end class ImageArithmGenericFilter
 
+  // specializations for itk::Vector<float, 3u>, 3u
   template<> template<>
   void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >();
   
@@ -118,6 +119,27 @@ namespace clitk {
      itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2, 
      itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito);
   
+  // specializations for itk::Vector<double, 3u>, 3u
+  template<> template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<double, 3u>, 3u > >();
+  
+  template<> template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage< 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > 
+    >
+    (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it, 
+     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito);
+
+  template<> template<>
+  void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage< 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >, 
+    itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > 
+    >
+    (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it1, 
+     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it2, 
+     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito);
 } // end namespace
 //--------------------------------------------------------------------
 
diff --git a/tools/clitkImageArithmGenericFilter.txx b/tools/clitkImageArithmGenericFilter.txx
index 8efd6a7..1afb77f 100644
--- a/tools/clitkImageArithmGenericFilter.txx
+++ b/tools/clitkImageArithmGenericFilter.txx
@@ -46,6 +46,7 @@ void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
 {
   ADD_DEFAULT_IMAGE_TYPES(Dim);
   ADD_VEC_IMAGE_TYPE(3u,3u,float);
+  ADD_VEC_IMAGE_TYPE(3u,3u,double);
 }
 //--------------------------------------------------------------------