filters/clitkImageArithmGenericFilter.txx

   1 #ifndef CLITKIMAGEARITHMGENERICFILTER_TXX
   2 #define CLITKIMAGEARITHMGENERICFILTER_TXX
   3 namespace clitk
   4 {
   5
   6   //--------------------------------------------------------------------
   7   template<class args_info_type>
   8   ImageArithmGenericFilter<args_info_type>::ImageArithmGenericFilter()
   9     :ImageToImageGenericFilter<Self>("ImageArithmGenericFilter"),mTypeOfOperation(0) {
  10     InitializeImageType<2>();
  11     InitializeImageType<3>();
  12     InitializeImageType<4>();
  13     mIsOperationUseASecondImage = false;
  14     mOverwriteInputImage = true;
  15   }
  16   //--------------------------------------------------------------------
  17
  18
  19   //--------------------------------------------------------------------
  20   template<class args_info_type>
  21   template<unsigned int Dim>
  22   void ImageArithmGenericFilter<args_info_type>::InitializeImageType() {
  23     ADD_IMAGE_TYPE(Dim, char);
  24     ADD_IMAGE_TYPE(Dim, uchar);
  25     ADD_IMAGE_TYPE(Dim, short);
  26     ADD_IMAGE_TYPE(Dim, float);
  27   }
  28   //--------------------------------------------------------------------
  29
  30
  31   //--------------------------------------------------------------------
  32   template<class args_info_type>
  33   void ImageArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b) {
  34     mOverwriteInputImage = b;
  35   }
  36   //--------------------------------------------------------------------
  37
  38
  39   //--------------------------------------------------------------------
  40   template<class args_info_type>
  41   void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a) {
  42     mArgsInfo=a;
  43
  44     // Set value
  45     SetIOVerbose(mArgsInfo.verbose_flag);
  46     mTypeOfOperation = mArgsInfo.operation_arg;
  47     mDefaultPixelValue = mArgsInfo.pixelValue_arg;
  48     mScalar = mArgsInfo.scalar_arg;
  49     mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
  50
  51     if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
  52
  53     if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
  54     if (mArgsInfo.input2_given) {
  55       mIsOperationUseASecondImage = true;
  56       AddInputFilename(mArgsInfo.input2_arg);
  57     }
  58
  59     if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
  60
  61     // Check type of operation (with scalar or with other image)
  62     if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
  63       std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
  64       exit(-1);
  65     }
  66     if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
  67       if (mArgsInfo.operation_arg < 5) {
  68         std::cerr << "Such operation need the --scalar option." << std::endl;
  69         exit(-1);
  70       }
  71     }
  72   }
  73   //--------------------------------------------------------------------
  74
  75
  76   //--------------------------------------------------------------------
  77   template<class args_info_type>
  78   template<class ImageType>
  79   void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType() {
  80     // Read input1
  81     typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
  82     typename ImageType::PixelType PixelType;
  83
  84     // Set input image iterator
  85     typedef itk::ImageRegionIterator<ImageType> IteratorType;
  86     IteratorType it(input1, input1->GetLargestPossibleRegion());
  87
  88     // typedef input2
  89     typename ImageType::Pointer input2 = this->template GetInput<ImageType>(1);
  90     IteratorType it2;
  91
  92     if (mIsOperationUseASecondImage) {
  93       // Read input2
  94       input2 = this->template GetInput<ImageType>(1);
  95       // Set input image iterator
  96       it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
  97     }
  98
  99     // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
 100     if (mOverwriteInputImage && mOutputIsFloat && (typeid(PixelType) != typeid(float))) {
 101       // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
 102       //                     << typeid(PixelType).name()
 103       //                     << std::endl;
 104       mOverwriteInputImage = false;
 105     }
 106
 107     // ---------------- Overwrite input Image ---------------------
 108     if (mOverwriteInputImage) {
 109       // Set output iterator (to input1)
 110       IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
 111       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 112       else ComputeImage(it, ito);
 113       this->template SetNextOutput<ImageType>(input1);
 114     }
 115
 116     // ---------------- Create new output Image ---------------------
 117     else {
 118       if (mOutputIsFloat) {
 119         // Create output image
 120         typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
 121         typename OutputImageType::Pointer output = OutputImageType::New();
 122         output->SetRegions(input1->GetLargestPossibleRegion());
 123         output->SetOrigin(input1->GetOrigin());
 124         output->SetSpacing(input1->GetSpacing());
 125         output->Allocate();
 126         // Set output iterator
 127         typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 128         IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 129         if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 130         else ComputeImage(it, ito);
 131         this->template SetNextOutput<OutputImageType>(output);
 132       }
 133       else {
 134         // Create output image
 135         typedef ImageType OutputImageType;
 136         typename OutputImageType::Pointer output = OutputImageType::New();
 137         output->SetRegions(input1->GetLargestPossibleRegion());
 138         output->SetOrigin(input1->GetOrigin());
 139         output->SetSpacing(input1->GetSpacing());
 140         output->Allocate();
 141         // Set output iterator
 142         typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 143         IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 144         if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 145         else ComputeImage(it, ito);
 146         this->template SetNextOutput<OutputImageType>(output);
 147       }
 148     }
 149   }
 150   //--------------------------------------------------------------------
 151
 152   //--------------------------------------------------------------------
 153   template<class args_info_type>
 154   template<class Iter1, class Iter2, class Iter3>
 155   void  ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito) {
 156     it1.GoToBegin();
 157     it2.GoToBegin();
 158     ito.GoToBegin();
 159     typedef typename Iter3::PixelType PixelType;
 160
 161     switch (mTypeOfOperation) {
 162     case 0: // Addition
 163       while (!ito.IsAtEnd()) {
 164         ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() + (double)it2.Get()) );
 165         ++it1; ++it2; ++ito;
 166       }
 167       break;
 168     case 1: // Multiply
 169       while (!ito.IsAtEnd()) {
 170         ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() * (double)it2.Get()) );
 171         ++it1; ++it2; ++ito;
 172       }
 173       break;
 174     case 2: // Divide
 175       while (!ito.IsAtEnd()) {
 176         if (it1.Get() != 0)
 177           ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
 178         else ito.Set(mDefaultPixelValue);
 179         ++it1; ++it2; ++ito;
 180       }
 181       break;
 182     case 3: // Max
 183       while (!ito.IsAtEnd()) {
 184         if (it1.Get() < it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
 185         ++it1; ++it2; ++ito;
 186       }
 187       break;
 188     case 4: // Min
 189       while (!ito.IsAtEnd()) {
 190         if (it1.Get() > it2.Get()) ito.Set(PixelTypeDownCast<double, PixelType>(it2.Get()));
 191         ++it1; ++it2; ++ito;
 192       }
 193       break;
 194     case 5: // Absolute difference
 195       while (!ito.IsAtEnd()) {
 196         ito.Set(PixelTypeDownCast<double, PixelType>(fabs(it2.Get()-it1.Get())));
 197         ++it1; ++it2; ++ito;
 198       }
 199       break;
 200     case 6: // Squared differences
 201       while (!ito.IsAtEnd()) {
 202         ito.Set(PixelTypeDownCast<double, PixelType>(pow((double)it1.Get()-(double)it2.Get(),2)));
 203         ++it1; ++it2; ++ito;
 204       }
 205       break;
 206     case 7: // Difference
 207       while (!ito.IsAtEnd()) {
 208         ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get()-(double)it2.Get()));
 209         ++it1; ++it2; ++ito;
 210       }
 211       break;
 212     case 8: // Relative Difference
 213       while (!ito.IsAtEnd()) {
 214         if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
 215         else ito.Set(0.0);
 216         ++it1; ++it2; ++ito;
 217       }
 218       break;
 219     default: // error ?
 220       std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 221       exit(-1);
 222     }
 223   }
 224   //--------------------------------------------------------------------
 225
 226
 227   //--------------------------------------------------------------------
 228   template<class args_info_type>
 229   template<class Iter1, class Iter2>
 230   void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito) {
 231     ito.GoToBegin();
 232     it.GoToBegin();
 233     typedef typename Iter2::PixelType PixelType;
 234
 235     // Perform operation
 236     switch (mTypeOfOperation) {
 237     case 0: // Addition
 238       while (!it.IsAtEnd()) {
 239         ito.Set(clitk::PixelTypeDownCast<double, PixelType>((double)it.Get() + mScalar) );
 240         ++it; ++ito;
 241       }
 242       break;
 243     case 1: // Multiply
 244       while (!it.IsAtEnd()) {
 245         ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get() * mScalar) );
 246         ++it; ++ito;
 247       }
 248       break;
 249     case 2: // Inverse
 250       while (!it.IsAtEnd()) {
 251         if (it.Get() != 0)
 252           ito.Set(PixelTypeDownCast<double, PixelType>(mScalar / (double)it.Get()));
 253         else ito.Set(mDefaultPixelValue);
 254         ++it; ++ito;
 255       }
 256       break;
 257     case 3: // Max
 258       while (!it.IsAtEnd()) {
 259         if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 260         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 261         ++it; ++ito;
 262       }
 263       break;
 264     case 4: // Min
 265       while (!it.IsAtEnd()) {
 266         if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 267         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 268         ++it; ++ito;
 269       }
 270       break;
 271     case 5: // Absolute value
 272       while (!it.IsAtEnd()) {
 273         if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
 274         // <= zero to avoid warning for unsigned types
 275         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 276         ++it; ++ito;
 277       }
 278       break;
 279     case 6: // Squared value
 280       while (!it.IsAtEnd()) {
 281         ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
 282         ++it; ++ito;
 283       }
 284       break;
 285     case 7: // Log
 286       while (!it.IsAtEnd()) {
 287         if (it.Get() > 0)
 288           ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
 289         else ito.Set(mDefaultPixelValue);
 290         ++it; ++ito;
 291       }
 292       break;
 293     case 8: // exp
 294       while (!it.IsAtEnd()) {
 295         ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
 296         ++it; ++ito;
 297       }
 298       break;
 299     case 9: // sqrt
 300       while (!it.IsAtEnd()) {
 301         if (it.Get() > 0)
 302           ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
 303         else {
 304           if (it.Get() ==0) ito.Set(0);
 305           else ito.Set(mDefaultPixelValue);
 306         }
 307         ++it; ++ito;
 308       }
 309       break;
 310     default: // error ?
 311       std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 312       exit(-1);
 313     }
 314   }
 315   //--------------------------------------------------------------------
 316
 317 } // end namespace
 318
 319 #endif  //#define CLITKIMAGEARITHMGENERICFILTER_TXX