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://www.centreleonberard.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 = NULL;
 111   IteratorType it2;
 112
 113
 114   if (mIsOperationUseASecondImage) {
 115       // Read input2
 116       input2 = this->template GetInput<ImageType>(1);
 117       // Set input image iterator
 118       it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
 119       // Check dimension
 120       if (!clitk::HaveSameSizeAndSpacing<ImageType, ImageType>(input1, input2)) {
 121           std::cerr << "* ERROR * the images (input and input2) must have the same size & spacing";
 122           return;
 123       }
 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     while (!ito.IsAtEnd()) {
 233       ito.Set(PixelTypeDownCast<double, PixelType>(fabs((double)it2.Get()-(double)it1.Get())));
 234       ++it1;
 235       ++it2;
 236       ++ito;
 237     }
 238     break;
 239   case 6: // Squared differences
 240     while (!ito.IsAtEnd()) {
 241       ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
 242       ++it1;
 243       ++it2;
 244       ++ito;
 245     }
 246     break;
 247   case 7: // Difference
 248     while (!ito.IsAtEnd()) {
 249       ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
 250       ++it1;
 251       ++it2;
 252       ++ito;
 253     }
 254     break;
 255   case 8: // Relative Difference
 256     while (!ito.IsAtEnd()) {
 257       if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
 258       else ito.Set(0.0);
 259       ++it1;
 260       ++it2;
 261       ++ito;
 262     }
 263     break;
 264   default: // error ?
 265     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 266     exit(-1);
 267   }
 268 }
 269 //--------------------------------------------------------------------
 270
 271
 272 //--------------------------------------------------------------------
 273 template<class args_info_type>
 274 template<class Iter1, class Iter2>
 275 void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
 276 {
 277   ito.GoToBegin();
 278   it.GoToBegin();
 279   typedef typename Iter2::PixelType PixelType;
 280
 281   // Perform operation
 282   switch (mTypeOfOperation) {
 283   case 0: // Addition
 284     while (!it.IsAtEnd()) {
 285       ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
 286       ++it;
 287       ++ito;
 288     }
 289     break;
 290   case 1: // Multiply
 291     while (!it.IsAtEnd()) {
 292       ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
 293       ++it;
 294       ++ito;
 295     }
 296     break;
 297   case 2: // Inverse
 298     while (!it.IsAtEnd()) {
 299       if (it.Get() != 0)
 300         ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
 301       else ito.Set(mDefaultPixelValue);
 302       ++it;
 303       ++ito;
 304     }
 305     break;
 306   case 3: // Max
 307     while (!it.IsAtEnd()) {
 308       if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 309       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 310       ++it;
 311       ++ito;
 312     }
 313     break;
 314   case 4: // Min
 315     while (!it.IsAtEnd()) {
 316       if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 317       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 318       ++it;
 319       ++ito;
 320     }
 321     break;
 322   case 5: // Absolute value
 323     while (!it.IsAtEnd()) {
 324       if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
 325       // <= zero to avoid warning for unsigned types
 326       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 327       ++it;
 328       ++ito;
 329     }
 330     break;
 331   case 6: // Squared value
 332     while (!it.IsAtEnd()) {
 333       ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
 334       ++it;
 335       ++ito;
 336     }
 337     break;
 338   case 7: // Log
 339     while (!it.IsAtEnd()) {
 340       if (it.Get() > 0)
 341         ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
 342       else ito.Set(mDefaultPixelValue);
 343       ++it;
 344       ++ito;
 345     }
 346     break;
 347   case 8: // exp
 348     while (!it.IsAtEnd()) {
 349       ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
 350       ++it;
 351       ++ito;
 352     }
 353     break;
 354   case 9: // sqrt
 355     while (!it.IsAtEnd()) {
 356       if (it.Get() > 0)
 357         ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
 358       else {
 359         if (it.Get() ==0) ito.Set(0);
 360         else ito.Set(mDefaultPixelValue);
 361       }
 362       ++it;
 363       ++ito;
 364     }
 365     break;
 366   case 10: // exp
 367     while (!it.IsAtEnd()) {
 368       ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
 369       ++it;
 370       ++ito;
 371     }
 372     break;
 373   default: // error ?
 374     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 375     exit(-1);
 376   }
 377 }
 378 //--------------------------------------------------------------------
 379
 380 } // end namespace
 381
 382 #endif  //#define CLITKIMAGEARITHMGENERICFILTER_TXX