filters/clitkImageArithmGenericFilter.txx

   1 /*-------------------------------------------------------------------------
   2
   3   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
   4   l'Image). All rights reserved. See Doc/License.txt or
   5   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
   6
   7   This software is distributed WITHOUT ANY WARRANTY; without even
   8   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
   9   PURPOSE.  See the above copyright notices for more information.
  10
  11   -------------------------------------------------------------------------*/
  12
  13 #ifndef CLITKIMAGEARITHMGENERICFILTER_TXX
  14 #define CLITKIMAGEARITHMGENERICFILTER_TXX
  15
  16 /*-------------------------------------------------
  17  * @file   clitkImageArithmGenericFilter.txx
  18  * @author David Sarrut <David.Sarrut@creatis.insa-lyon.fr>
  19  * @date   9 August 2006
  20  *
  21  -------------------------------------------------*/
  22
  23 //--------------------------------------------------------------------
  24 template<unsigned int Dim>
  25 void ImageArithmGenericFilter::Update_WithDim() {
  26 #define TRY_TYPE(TYPE)                                                  \
  27   if (IsSameType<TYPE>(mPixelTypeName)) { Update_WithDimAndPixelType<Dim, TYPE>(); return; }
  28   TRY_TYPE(char);
  29   TRY_TYPE(uchar);
  30   TRY_TYPE(short);
  31   TRY_TYPE(ushort);
  32   TRY_TYPE(int);
  33   TRY_TYPE(unsigned int);
  34   TRY_TYPE(float);
  35   TRY_TYPE(double);
  36 #undef TRY_TYPE
  37
  38   std::string list = CreateListOfTypes<char, uchar, short, ushort, int, uint, float, double>();
  39   std::cerr << "Error, I don't know the type '" << mPixelTypeName << "' for the input image '"
  40             << mInputFilenames[0] << "'." << std::endl << "Known types are " << list << std::endl;
  41   exit(0);
  42 }
  43 //--------------------------------------------------------------------
  44
  45 //--------------------------------------------------------------------
  46 template<unsigned int Dim, class PixelType>
  47 void ImageArithmGenericFilter::Update_WithDimAndPixelType() {
  48   // Read input1
  49   typedef itk::Image<PixelType,Dim> ImageType;
  50   typename ImageType::Pointer input1 = clitk::readImage<ImageType>(mInputFilenames[0], mIOVerbose);
  51   typename ImageType::Pointer outputImage;
  52
  53   // Read input2 (float is ok altough it could take too memory ...)
  54   if (mIsOperationUseASecondImage) {
  55     typedef itk::Image<float,Dim> ImageType2;
  56     typename ImageType2::Pointer input2 = clitk::readImage<ImageType2>(mInputFilenames[1], mIOVerbose);
  57     outputImage = ComputeImage<ImageType, ImageType2>(input1, input2);
  58   }
  59   else {
  60     outputImage = ComputeImage<ImageType>(input1);
  61   }
  62
  63   // Write results
  64   this->SetNextOutput<ImageType>(outputImage);
  65   //clitk::writeImage<ImageType>(outputImage, mOutputFilename, mIOVerbose);
  66 }
  67 //--------------------------------------------------------------------
  68
  69 //--------------------------------------------------------------------
  70 template<class ImageType>
  71 typename ImageType::Pointer
  72 ImageArithmGenericFilter::ComputeImage(typename ImageType::Pointer inputImage) {
  73
  74   typedef typename ImageType::PixelType PixelType;
  75   typedef itk::ImageRegionIterator<ImageType> IteratorType;
  76   IteratorType it(inputImage, inputImage->GetLargestPossibleRegion());
  77   it.GoToBegin();
  78
  79   switch (mTypeOfOperation) {
  80   case 0: // Addition
  81     while (!it.IsAtEnd()) {
  82       it.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
  83       ++it;
  84     }
  85     break;
  86   case 1: // Multiply
  87     while (!it.IsAtEnd()) {
  88       it.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
  89       ++it;
  90     }
  91     break;
  92   case 2: // Inverse
  93     while (!it.IsAtEnd()) {
  94       if (it.Get() != 0)
  95         it.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
  96       else it.Set(mDefaultPixelValue);
  97       ++it;
  98     }
  99     break;
 100   case 3: // Max
 101     while (!it.IsAtEnd()) {
 102       if (it.Get() < mScalar) it.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 103       ++it;
 104     }
 105     break;
 106   case 4: // Min
 107     while (!it.IsAtEnd()) {
 108       if (it.Get() > mScalar) it.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 109       ++it;
 110     }
 111     break;
 112   case 5: // Absolute value
 113     while (!it.IsAtEnd()) {
 114       if (it.Get() <= 0) it.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
 115       // <= zero to avoid warning for unsigned types
 116       ++it;
 117     }
 118     break;
 119   case 6: // Squared value
 120     while (!it.IsAtEnd()) {
 121       it.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
 122       ++it;
 123     }
 124     break;
 125   case 7: // Log
 126     while (!it.IsAtEnd()) {
 127       if (it.Get() > 0)
 128         it.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
 129       else it.Set(mDefaultPixelValue);
 130       ++it;
 131     }
 132     break;
 133   case 8: // exp
 134     while (!it.IsAtEnd()) {
 135       it.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
 136       ++it;
 137     }
 138     break;
 139   case 9: // sqrt
 140     while (!it.IsAtEnd()) {
 141       if (it.Get() > 0)
 142         it.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
 143       else {
 144         if (it.Get() ==0) it.Set(0);
 145         else it.Set(mDefaultPixelValue);
 146       }
 147       ++it;
 148     }
 149     break;
 150   default: // error ?
 151     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 152     exit(-1);
 153   }
 154
 155   return inputImage;
 156 }
 157 //--------------------------------------------------------------------
 158
 159 //--------------------------------------------------------------------
 160 template<class ImageType1, class ImageType2>
 161 typename ImageType1::Pointer
 162 ImageArithmGenericFilter::ComputeImage(typename ImageType1::Pointer inputImage1,
 163                                        typename ImageType2::Pointer inputImage2) {
 164
 165   typedef typename ImageType1::PixelType PixelType;
 166   typedef itk::ImageRegionIterator<ImageType1> IteratorType;
 167   IteratorType it1(inputImage1, inputImage1->GetLargestPossibleRegion());
 168   it1.GoToBegin();
 169
 170   typedef itk::ImageRegionConstIterator<ImageType2> ConstIteratorType;
 171   ConstIteratorType it2(inputImage2, inputImage2->GetLargestPossibleRegion());
 172   it2.GoToBegin();
 173
 174   switch (mTypeOfOperation) {
 175   case 0: // Addition
 176     while (!it1.IsAtEnd()) {
 177       it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
 178       ++it1; ++it2;
 179     }
 180     break;
 181   case 1: // Multiply
 182     while (!it1.IsAtEnd()) {
 183       it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
 184       ++it1; ++it2;
 185     }
 186     break;
 187   case 2: // Divide
 188     while (!it1.IsAtEnd()) {
 189       if (it1.Get() != 0)
 190         it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
 191       else it1.Set(mDefaultPixelValue);
 192       ++it1; ++it2;
 193     }
 194     break;
 195   case 3: // Max
 196     while (!it1.IsAtEnd()) {
 197       if (it1.Get() < it2.Get()) it1.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
 198       ++it1; ++it2;
 199     }
 200     break;
 201   case 4: // Min
 202     while (!it1.IsAtEnd()) {
 203       if (it1.Get() > it2.Get()) it1.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
 204       ++it1; ++it2;
 205     }
 206     break;
 207   case 5: // Absolute difference
 208     while (!it1.IsAtEnd()) {
 209       it1.Set(PixelTypeDownCast<double, PixelType>(fabs(it2.Get()-it1.Get())));
 210       ++it1; ++it2;
 211     }
 212     break;
 213   case 6: // Squared differences
 214     while (!it1.IsAtEnd()) {
 215       it1.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
 216       ++it1; ++it2;
 217     }
 218     break;
 219   case 7: // Difference
 220     while (!it1.IsAtEnd()) {
 221       it1.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
 222       ++it1; ++it2;
 223     }
 224     break;
 225   case 8: // Relative Difference
 226     while (!it1.IsAtEnd()) {
 227       if (it1.Get() != 0) it1.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
 228       else it1.Set(0.0);
 229       ++it1; ++it2;
 230     }
 231     break;
 232   default: // error ?
 233     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 234     exit(-1);
 235   }
 236
 237   return inputImage1;
 238 }
 239 //--------------------------------------------------------------------
 240
 241 #endif  //#define CLITKIMAGEARITHMGENERICFILTER_TXX