tools/clitkImageArithmGenericFilter.cxx

   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_CXX
  19 #define CLITKIMAGEARITHMGENERICFILTER_CXX
  20
  21 #include "clitkImageArithmGenericFilter.h"
  22
  23 namespace clitk {
  24   // Specialisation
  25 //   template<>
  26 //   class ImageArithmGenericFilter<args_info_clitkImageArithm>;
  27
  28   template<>
  29   template<>
  30   void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >()
  31   {
  32     typedef itk::Image< itk::Vector<float, 3u>, 3u > ImageType;
  33
  34     // Read input1
  35     ImageType::Pointer input1 = this->GetInput<ImageType>(0);
  36
  37     // Set input image iterator
  38     typedef itk::ImageRegionIterator<ImageType> IteratorType;
  39     IteratorType it(input1, input1->GetLargestPossibleRegion());
  40
  41     // typedef input2
  42     ImageType::Pointer input2 = NULL;
  43     IteratorType it2;
  44
  45     /*
  46     // Special case for normalisation
  47     if (mTypeOfOperation == 12) {
  48       typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
  49       typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
  50       ff->SetImage(input1);
  51       ff->ComputeMaximum();
  52       mScalar = ff->GetMaximum();
  53       mTypeOfOperation = 11; // divide
  54     }
  55     */
  56
  57     if (mIsOperationUseASecondImage) {
  58         // Read input2
  59         input2 = this->GetInput<ImageType>(1);
  60         // Set input image iterator
  61         it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
  62         // Check dimension
  63         if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
  64           itkExceptionMacro(<< "The images (input and input2) must have the same size");
  65         }
  66         if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
  67           itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
  68                           << "Using first input's information.");
  69         }
  70     }
  71
  72     /*
  73     // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
  74     if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
  75       // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
  76       //                     << typeid(PixelType).name()
  77       //                     << std::endl;
  78       mOverwriteInputImage = false;
  79     }
  80     */
  81
  82     // ---------------- Overwrite input Image ---------------------
  83     if (mOverwriteInputImage) {
  84       // Set output iterator (to input1)
  85       IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
  86       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
  87       else ComputeImage(it, ito);
  88       this->SetNextOutput<ImageType>(input1);
  89     }
  90     // ---------------- Create new output Image ---------------------
  91     else {
  92       /*if (mOutputIsFloat) {
  93         // Create output image
  94         typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
  95         typename OutputImageType::Pointer output = OutputImageType::New();
  96         output->SetRegions(input1->GetLargestPossibleRegion());
  97         output->SetOrigin(input1->GetOrigin());
  98         output->SetSpacing(input1->GetSpacing());
  99         output->Allocate();
 100         // Set output iterator
 101         typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 102         IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 103         if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 104         else ComputeImage(it, ito);
 105         this->template SetNextOutput<OutputImageType>(output);
 106       } else*/ {
 107         // Create output image
 108         typedef ImageType OutputImageType;
 109         OutputImageType::Pointer output = OutputImageType::New();
 110         output->SetRegions(input1->GetLargestPossibleRegion());
 111         output->SetOrigin(input1->GetOrigin());
 112         output->SetSpacing(input1->GetSpacing());
 113         output->Allocate();
 114         // Set output iterator
 115         typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 116         IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 117         if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 118         else ComputeImage(it, ito);
 119         this->SetNextOutput<OutputImageType>(output);
 120       }
 121     }
 122   }
 123
 124   template<>
 125   template<>
 126   void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
 127     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
 128     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > >
 129     (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it,
 130     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
 131   {
 132     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
 133
 134     ito.GoToBegin();
 135     it.GoToBegin();
 136
 137     typedef Iter2::PixelType PixelType;
 138
 139     PixelType scalar_vector;
 140     scalar_vector.Fill(mScalar);
 141
 142     // Perform operation
 143     switch (mTypeOfOperation) {
 144     case 0: // Addition
 145       while (!it.IsAtEnd()) {
 146         ito.Set(it.Get() + scalar_vector);
 147         ++it;
 148         ++ito;
 149       }
 150       break;
 151     case 1: // Multiply
 152       while (!it.IsAtEnd()) {
 153         ito.Set(it.Get() * mScalar);
 154         ++it;
 155         ++ito;
 156       }
 157       break;
 158       /*
 159     case 2: // Inverse
 160       while (!it.IsAtEnd()) {
 161         if (it.Get() != 0)
 162           ito.Set(mScalar / it.Get()));
 163         else ito.Set(mDefaultPixelValue);
 164         ++it;
 165         ++ito;
 166       }
 167       break;
 168     case 3: // Max
 169       while (!it.IsAtEnd()) {
 170         if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 171         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 172         ++it;
 173         ++ito;
 174       }
 175       break;
 176     case 4: // Min
 177       while (!it.IsAtEnd()) {
 178         if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 179         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 180         ++it;
 181         ++ito;
 182       }
 183       break;
 184     case 5: // Absolute value
 185       while (!it.IsAtEnd()) {
 186         if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
 187         // <= zero to avoid warning for unsigned types
 188         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 189         ++it;
 190         ++ito;
 191       }
 192       break;
 193     case 6: // Squared value
 194       while (!it.IsAtEnd()) {
 195         ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
 196         ++it;
 197         ++ito;
 198       }
 199       break;
 200     case 7: // Log
 201       while (!it.IsAtEnd()) {
 202         if (it.Get() > 0)
 203           ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
 204         else ito.Set(mDefaultPixelValue);
 205         ++it;
 206         ++ito;
 207       }
 208       break;
 209     case 8: // exp
 210       while (!it.IsAtEnd()) {
 211         ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
 212         ++it;
 213         ++ito;
 214       }
 215       break;
 216     case 9: // sqrt
 217       while (!it.IsAtEnd()) {
 218         if (it.Get() > 0)
 219           ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
 220         else {
 221           if (it.Get() ==0) ito.Set(0);
 222           else ito.Set(mDefaultPixelValue);
 223         }
 224         ++it;
 225         ++ito;
 226       }
 227       break;
 228     case 10: // exp
 229       while (!it.IsAtEnd()) {
 230         ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
 231         ++it;
 232         ++ito;
 233       }
 234       break;
 235       */
 236     case 11: // divide
 237       while (!it.IsAtEnd()) {
 238         ito.Set(it.Get() / mScalar);
 239         ++it;
 240         ++ito;
 241       }
 242       break;
 243     default: // error ?
 244       std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 245       exit(-1);
 246     }
 247
 248   }
 249
 250   template<>
 251   template<>
 252   void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
 253     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
 254     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
 255     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
 256     >
 257     (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1,
 258     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2,
 259     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
 260   {
 261     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter1;
 262     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
 263     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter3;
 264
 265     it1.GoToBegin();
 266     it2.GoToBegin();
 267     ito.GoToBegin();
 268     typedef Iter3::PixelType PixelType;
 269
 270     switch (mTypeOfOperation) {
 271     case 0: // Addition
 272       while (!ito.IsAtEnd()) {
 273         ito.Set(it1.Get() + it2.Get());
 274         ++it1;
 275         ++it2;
 276         ++ito;
 277       }
 278       break;
 279       /*
 280     case 1: // Multiply
 281       while (!ito.IsAtEnd()) {
 282         ito.Set(it1.Get() * it2.Get()) );
 283         ++it1;
 284         ++it2;
 285         ++ito;
 286       }
 287       break;
 288     case 2: // Divide
 289       while (!ito.IsAtEnd()) {
 290         if (it1.Get() != 0)
 291           ito.Set(it1.Get() / it2.Get()));
 292         else ito.Set(mDefaultPixelValue);
 293         ++it1;
 294         ++it2;
 295         ++ito;
 296       }
 297       break;
 298     case 3: // Max
 299       while (!ito.IsAtEnd()) {
 300         if (it1.Get() < it2.Get()) ito.Set(it2.Get());
 301         ++it1;
 302         ++it2;
 303         ++ito;
 304       }
 305       break;
 306     case 4: // Min
 307       while (!ito.IsAtEnd()) {
 308         if (it1.Get() > it2.Get()) ito.Set(it2.Get());
 309         ++it1;
 310         ++it2;
 311         ++ito;
 312       }
 313       break;
 314       */
 315     case 5: // Absolute difference
 316       while (!ito.IsAtEnd()) {
 317         ito.Set(it2.Get()-it1.Get());
 318         ++it1;
 319         ++it2;
 320         ++ito;
 321       }
 322       break;
 323       /*
 324     case 6: // Squared differences
 325       while (!ito.IsAtEnd()) {
 326         ito.Set(pow(it1.Get()-it2.Get(),2)));
 327         ++it1;
 328         ++it2;
 329         ++ito;
 330       }
 331       break;
 332       */
 333     case 7: // Difference
 334       while (!ito.IsAtEnd()) {
 335         ito.Set(it1.Get()-it2.Get());
 336         ++it1;
 337         ++it2;
 338         ++ito;
 339       }
 340       break;
 341       /*
 342     case 8: // Relative Difference
 343       while (!ito.IsAtEnd()) {
 344         if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
 345         else ito.Set(0.0);
 346         ++it1;
 347         ++it2;
 348         ++ito;
 349       }
 350       break;
 351       */
 352     default: // error ?
 353       std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 354       exit(-1);
 355     }
 356   }
 357
 358 }
 359
 360 #endif //define CLITKIMAGEARITHMGENERICFILTER_CXX