tools/clitkImageArithmGenericFilter.txx

   1 /*=========================================================================
   2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
   3
   4   Authors belong to:
   5   - University of LYON              http://www.universite-lyon.fr/
   6   - Léon Bérard cancer center       http://oncora1.lyon.fnclcc.fr
   7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
   8
   9   This software is distributed WITHOUT ANY WARRANTY; without even
  10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  11   PURPOSE.  See the copyright notices for more information.
  12
  13   It is distributed under dual licence
  14
  15   - BSD        See included LICENSE.txt file
  16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
  17 ======================================================================-====*/
  18 #ifndef CLITKIMAGEARITHMGENERICFILTER_TXX
  19 #define CLITKIMAGEARITHMGENERICFILTER_TXX
  20
  21 #include "clitkImageCommon.h"
  22
  23 namespace clitk
  24 {
  25
  26 //--------------------------------------------------------------------
  27 template<class args_info_type>
  28 ImageArithmGenericFilter<args_info_type>::ImageArithmGenericFilter()
  29   :ImageToImageGenericFilter<Self>("ImageArithmGenericFilter"),mTypeOfOperation(0)
  30 {
  31   InitializeImageType<2>();
  32   InitializeImageType<3>();
  33   InitializeImageType<4>();
  34   mIsOperationUseASecondImage = false;
  35   mOverwriteInputImage = true;
  36 }
  37 //--------------------------------------------------------------------
  38
  39
  40 //--------------------------------------------------------------------
  41 template<class args_info_type>
  42 template<unsigned int Dim>
  43 void ImageArithmGenericFilter<args_info_type>::InitializeImageType()
  44 {
  45   ADD_DEFAULT_IMAGE_TYPES(Dim);
  46 }
  47 //--------------------------------------------------------------------
  48
  49
  50 //--------------------------------------------------------------------
  51 template<class args_info_type>
  52 void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
  53 {
  54   mOverwriteInputImage = b;
  55 }
  56 //--------------------------------------------------------------------
  57
  58
  59 //--------------------------------------------------------------------
  60 template<class args_info_type>
  61 void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
  62 {
  63   mArgsInfo=a;
  64
  65   // Set value
  66   SetIOVerbose(mArgsInfo.verbose_flag);
  67   mTypeOfOperation = mArgsInfo.operation_arg;
  68   mDefaultPixelValue = mArgsInfo.pixelValue_arg;
  69   mScalar = mArgsInfo.scalar_arg;
  70   mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
  71
  72   if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
  73
  74   if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
  75   if (mArgsInfo.input2_given) {
  76     mIsOperationUseASecondImage = true;
  77     AddInputFilename(mArgsInfo.input2_arg);
  78   }
  79
  80   if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
  81
  82   // Check type of operation (with scalar or with other image)
  83   if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
  84     std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
  85     exit(-1);
  86   }
  87   if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
  88     if (mArgsInfo.operation_arg < 5) {
  89       std::cerr << "Such operation need the --scalar option." << std::endl;
  90       exit(-1);
  91     }
  92   }
  93 }
  94 //--------------------------------------------------------------------
  95
  96
  97 //--------------------------------------------------------------------
  98 template<class args_info_type>
  99 template<class ImageType>
 100 void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
 101 {
 102   // Read input1
 103   typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
 104
 105   // Set input image iterator
 106   typedef itk::ImageRegionIterator<ImageType> IteratorType;
 107   IteratorType it(input1, input1->GetLargestPossibleRegion());
 108
 109   // typedef input2
 110   typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
 111   IteratorType it2;
 112
 113   // Check dimension
 114   if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input1, input2)) {
 115     std::cerr << "* ERROR * the images (input and input2) must have the same size & spacing";
 116     return;
 117   }
 118
 119   if (mIsOperationUseASecondImage) {
 120     // Read input2
 121     input2 = this->template GetInput<ImageType>(1);
 122     // Set input image iterator
 123     it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
 124   }
 125
 126   // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
 127   if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
 128     // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
 129     //                     << typeid(PixelType).name()
 130     //                     << std::endl;
 131     mOverwriteInputImage = false;
 132   }
 133
 134   // ---------------- Overwrite input Image ---------------------
 135   if (mOverwriteInputImage) {
 136     // Set output iterator (to input1)
 137     IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
 138     if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 139     else ComputeImage(it, ito);
 140     this->template SetNextOutput<ImageType>(input1);
 141   }
 142
 143   // ---------------- Create new output Image ---------------------
 144   else {
 145     if (mOutputIsFloat) {
 146       // Create output image
 147       typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
 148       typename OutputImageType::Pointer output = OutputImageType::New();
 149       output->SetRegions(input1->GetLargestPossibleRegion());
 150       output->SetOrigin(input1->GetOrigin());
 151       output->SetSpacing(input1->GetSpacing());
 152       output->Allocate();
 153       // Set output iterator
 154       typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 155       IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 156       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 157       else ComputeImage(it, ito);
 158       this->template SetNextOutput<OutputImageType>(output);
 159     } else {
 160       // Create output image
 161       typedef ImageType OutputImageType;
 162       typename OutputImageType::Pointer output = OutputImageType::New();
 163       output->SetRegions(input1->GetLargestPossibleRegion());
 164       output->SetOrigin(input1->GetOrigin());
 165       output->SetSpacing(input1->GetSpacing());
 166       output->Allocate();
 167       // Set output iterator
 168       typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 169       IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 170       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 171       else ComputeImage(it, ito);
 172       this->template SetNextOutput<OutputImageType>(output);
 173     }
 174   }
 175 }
 176 //--------------------------------------------------------------------
 177
 178 //--------------------------------------------------------------------
 179 template<class args_info_type>
 180 template<class Iter1, class Iter2, class Iter3>
 181 void  ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
 182 {
 183   it1.GoToBegin();
 184   it2.GoToBegin();
 185   ito.GoToBegin();
 186   typedef typename Iter3::PixelType PixelType;
 187
 188   switch (mTypeOfOperation) {
 189   case 0: // Addition
 190     while (!ito.IsAtEnd()) {
 191       ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
 192       ++it1;
 193       ++it2;
 194       ++ito;
 195     }
 196     break;
 197   case 1: // Multiply
 198     while (!ito.IsAtEnd()) {
 199       ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
 200       ++it1;
 201       ++it2;
 202       ++ito;
 203     }
 204     break;
 205   case 2: // Divide
 206     while (!ito.IsAtEnd()) {
 207       if (it1.Get() != 0)
 208         ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
 209       else ito.Set(mDefaultPixelValue);
 210       ++it1;
 211       ++it2;
 212       ++ito;
 213     }
 214     break;
 215   case 3: // Max
 216     while (!ito.IsAtEnd()) {
 217       if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
 218       ++it1;
 219       ++it2;
 220       ++ito;
 221     }
 222     break;
 223   case 4: // Min
 224     while (!ito.IsAtEnd()) {
 225       if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
 226       ++it1;
 227       ++it2;
 228       ++ito;
 229     }
 230     break;
 231   case 5: // Absolute difference
 232       DD("AbsoluteDifff");
 233     while (!ito.IsAtEnd()) {
 234       ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
 235       ++it1;
 236       ++it2;
 237       ++ito;
 238     }
 239     break;
 240   case 6: // Squared differences
 241     while (!ito.IsAtEnd()) {
 242       ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
 243       ++it1;
 244       ++it2;
 245       ++ito;
 246     }
 247     break;
 248   case 7: // Difference
 249     while (!ito.IsAtEnd()) {
 250       ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
 251       ++it1;
 252       ++it2;
 253       ++ito;
 254     }
 255     break;
 256   case 8: // Relative Difference
 257     while (!ito.IsAtEnd()) {
 258       if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
 259       else ito.Set(0.0);
 260       ++it1;
 261       ++it2;
 262       ++ito;
 263     }
 264     break;
 265   default: // error ?
 266     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 267     exit(-1);
 268   }
 269 }
 270 //--------------------------------------------------------------------
 271
 272
 273 //--------------------------------------------------------------------
 274 template<class args_info_type>
 275 template<class Iter1, class Iter2>
 276 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
 277 {
 278   ito.GoToBegin();
 279   it.GoToBegin();
 280   typedef typename Iter2::PixelType PixelType;
 281
 282   // Perform operation
 283   switch (mTypeOfOperation) {
 284   case 0: // Addition
 285     while (!it.IsAtEnd()) {
 286       ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
 287       ++it;
 288       ++ito;
 289     }
 290     break;
 291   case 1: // Multiply
 292     while (!it.IsAtEnd()) {
 293       ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
 294       ++it;
 295       ++ito;
 296     }
 297     break;
 298   case 2: // Inverse
 299     while (!it.IsAtEnd()) {
 300       if (it.Get() != 0)
 301         ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
 302       else ito.Set(mDefaultPixelValue);
 303       ++it;
 304       ++ito;
 305     }
 306     break;
 307   case 3: // Max
 308     while (!it.IsAtEnd()) {
 309       if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 310       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 311       ++it;
 312       ++ito;
 313     }
 314     break;
 315   case 4: // Min
 316     while (!it.IsAtEnd()) {
 317       if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 318       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 319       ++it;
 320       ++ito;
 321     }
 322     break;
 323   case 5: // Absolute value
 324     while (!it.IsAtEnd()) {
 325       if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
 326       // <= zero to avoid warning for unsigned types
 327       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 328       ++it;
 329       ++ito;
 330     }
 331     break;
 332   case 6: // Squared value
 333     while (!it.IsAtEnd()) {
 334       ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
 335       ++it;
 336       ++ito;
 337     }
 338     break;
 339   case 7: // Log
 340     while (!it.IsAtEnd()) {
 341       if (it.Get() > 0)
 342         ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
 343       else ito.Set(mDefaultPixelValue);
 344       ++it;
 345       ++ito;
 346     }
 347     break;
 348   case 8: // exp
 349     while (!it.IsAtEnd()) {
 350       ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
 351       ++it;
 352       ++ito;
 353     }
 354     break;
 355   case 9: // sqrt
 356     while (!it.IsAtEnd()) {
 357       if (it.Get() > 0)
 358         ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
 359       else {
 360         if (it.Get() ==0) ito.Set(0);
 361         else ito.Set(mDefaultPixelValue);
 362       }
 363       ++it;
 364       ++ito;
 365     }
 366     break;
 367   default: // error ?
 368     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 369     exit(-1);
 370   }
 371 }
 372 //--------------------------------------------------------------------
 373
 374 } // end namespace
 375
 376 #endif  //#define CLITKIMAGEARITHMGENERICFILTER_TXX