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 ////////////////////////////////////////////////////////////////////
  29   // specializations for itk::Vector<float, 3u>, 3u
  30   template<>
  31   template<>
  32   void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<float, 3u>, 3u > >()
  33   {
  34     typedef itk::Image< itk::Vector<float, 3u>, 3u > ImageType;
  35
  36     // Read input1
  37     ImageType::Pointer input1 = this->GetInput<ImageType>(0);
  38
  39     // Set input image iterator
  40     typedef itk::ImageRegionIterator<ImageType> IteratorType;
  41     IteratorType it(input1, input1->GetLargestPossibleRegion());
  42
  43     // typedef input2
  44     ImageType::Pointer input2 = NULL;
  45     IteratorType it2;
  46
  47     /*
  48     // Special case for normalisation
  49     if (mTypeOfOperation == 12) {
  50       typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
  51       typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
  52       ff->SetImage(input1);
  53       ff->ComputeMaximum();
  54       mScalar = ff->GetMaximum();
  55       mTypeOfOperation = 11; // divide
  56     }
  57     */
  58
  59     if (mIsOperationUseASecondImage) {
  60         // Read input2
  61         input2 = this->GetInput<ImageType>(1);
  62         // Set input image iterator
  63         it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
  64         // Check dimension
  65         if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
  66           itkExceptionMacro(<< "The images (input and input2) must have the same size");
  67         }
  68         if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
  69           itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
  70                           << "Using first input's information.");
  71         }
  72     }
  73
  74     /*
  75     // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
  76     if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
  77       // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
  78       //                     << typeid(PixelType).name()
  79       //                     << std::endl;
  80       mOverwriteInputImage = false;
  81     }
  82     */
  83
  84     // ---------------- Overwrite input Image ---------------------
  85     if (mOverwriteInputImage) {
  86       // Set output iterator (to input1)
  87       IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
  88       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
  89       else ComputeImage(it, ito);
  90       this->SetNextOutput<ImageType>(input1);
  91     }
  92     // ---------------- Create new output Image ---------------------
  93     else {
  94       /*if (mOutputIsFloat) {
  95         // Create output image
  96         typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
  97         typename OutputImageType::Pointer output = OutputImageType::New();
  98         output->SetRegions(input1->GetLargestPossibleRegion());
  99         output->SetOrigin(input1->GetOrigin());
 100         output->SetSpacing(input1->GetSpacing());
 101         output->Allocate();
 102         // Set output iterator
 103         typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 104         IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 105         if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 106         else ComputeImage(it, ito);
 107         this->template SetNextOutput<OutputImageType>(output);
 108       } else*/ {
 109         // Create output image
 110         typedef ImageType OutputImageType;
 111         OutputImageType::Pointer output = OutputImageType::New();
 112         output->SetRegions(input1->GetLargestPossibleRegion());
 113         output->SetOrigin(input1->GetOrigin());
 114         output->SetSpacing(input1->GetSpacing());
 115         output->Allocate();
 116         // Set output iterator
 117         typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 118         IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 119         if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 120         else ComputeImage(it, ito);
 121         this->SetNextOutput<OutputImageType>(output);
 122       }
 123     }
 124   }
 125
 126   template<>
 127   template<>
 128   void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
 129     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
 130     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > >
 131     (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it,
 132     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
 133   {
 134     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
 135
 136     ito.GoToBegin();
 137     it.GoToBegin();
 138
 139     typedef Iter2::PixelType PixelType;
 140
 141     PixelType scalar_vector;
 142     scalar_vector.Fill(mScalar);
 143
 144     // Perform operation
 145     switch (mTypeOfOperation) {
 146     case 0: // Addition
 147       while (!it.IsAtEnd()) {
 148         ito.Set(it.Get() + scalar_vector);
 149         ++it;
 150         ++ito;
 151       }
 152       break;
 153     case 1: // Multiply
 154       while (!it.IsAtEnd()) {
 155         ito.Set(it.Get() * mScalar);
 156         ++it;
 157         ++ito;
 158       }
 159       break;
 160       /*
 161     case 2: // Inverse
 162       while (!it.IsAtEnd()) {
 163         if (it.Get() != 0)
 164           ito.Set(mScalar / it.Get()));
 165         else ito.Set(mDefaultPixelValue);
 166         ++it;
 167         ++ito;
 168       }
 169       break;
 170     case 3: // Max
 171       while (!it.IsAtEnd()) {
 172         if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 173         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 174         ++it;
 175         ++ito;
 176       }
 177       break;
 178     case 4: // Min
 179       while (!it.IsAtEnd()) {
 180         if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 181         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 182         ++it;
 183         ++ito;
 184       }
 185       break;
 186     case 5: // Absolute value
 187       while (!it.IsAtEnd()) {
 188         if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
 189         // <= zero to avoid warning for unsigned types
 190         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 191         ++it;
 192         ++ito;
 193       }
 194       break;
 195     case 6: // Squared value
 196       while (!it.IsAtEnd()) {
 197         ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
 198         ++it;
 199         ++ito;
 200       }
 201       break;
 202     case 7: // Log
 203       while (!it.IsAtEnd()) {
 204         if (it.Get() > 0)
 205           ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
 206         else ito.Set(mDefaultPixelValue);
 207         ++it;
 208         ++ito;
 209       }
 210       break;
 211     case 8: // exp
 212       while (!it.IsAtEnd()) {
 213         ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
 214         ++it;
 215         ++ito;
 216       }
 217       break;
 218     case 9: // sqrt
 219       while (!it.IsAtEnd()) {
 220         if (it.Get() > 0)
 221           ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
 222         else {
 223           if (it.Get() ==0) ito.Set(0);
 224           else ito.Set(mDefaultPixelValue);
 225         }
 226         ++it;
 227         ++ito;
 228       }
 229       break;
 230     case 10: // exp
 231       while (!it.IsAtEnd()) {
 232         ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
 233         ++it;
 234         ++ito;
 235       }
 236       break;
 237       */
 238     case 11: // divide
 239       while (!it.IsAtEnd()) {
 240         ito.Set(it.Get() / mScalar);
 241         ++it;
 242         ++ito;
 243       }
 244       break;
 245     default: // error ?
 246       std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 247       exit(-1);
 248     }
 249
 250   }
 251
 252   template<>
 253   template<>
 254   void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
 255     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
 256     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >,
 257     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > >
 258     >
 259     (itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it1,
 260     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > it2,
 261     itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > ito)
 262   {
 263     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter1;
 264     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter2;
 265     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<float, 3u>, 3u > > Iter3;
 266
 267     it1.GoToBegin();
 268     it2.GoToBegin();
 269     ito.GoToBegin();
 270     typedef Iter3::PixelType PixelType;
 271
 272     switch (mTypeOfOperation) {
 273     case 0: // Addition
 274       while (!ito.IsAtEnd()) {
 275         ito.Set(it1.Get() + it2.Get());
 276         ++it1;
 277         ++it2;
 278         ++ito;
 279       }
 280       break;
 281       /*
 282     case 1: // Multiply
 283       while (!ito.IsAtEnd()) {
 284         ito.Set(it1.Get() * it2.Get()) );
 285         ++it1;
 286         ++it2;
 287         ++ito;
 288       }
 289       break;
 290     case 2: // Divide
 291       while (!ito.IsAtEnd()) {
 292         if (it1.Get() != 0)
 293           ito.Set(it1.Get() / it2.Get()));
 294         else ito.Set(mDefaultPixelValue);
 295         ++it1;
 296         ++it2;
 297         ++ito;
 298       }
 299       break;
 300     case 3: // Max
 301       while (!ito.IsAtEnd()) {
 302         if (it1.Get() < it2.Get()) ito.Set(it2.Get());
 303         ++it1;
 304         ++it2;
 305         ++ito;
 306       }
 307       break;
 308     case 4: // Min
 309       while (!ito.IsAtEnd()) {
 310         if (it1.Get() > it2.Get()) ito.Set(it2.Get());
 311         ++it1;
 312         ++it2;
 313         ++ito;
 314       }
 315       break;
 316       */
 317     case 5: // Absolute difference
 318       while (!ito.IsAtEnd()) {
 319         ito.Set(it2.Get()-it1.Get());
 320         ++it1;
 321         ++it2;
 322         ++ito;
 323       }
 324       break;
 325       /*
 326     case 6: // Squared differences
 327       while (!ito.IsAtEnd()) {
 328         ito.Set(pow(it1.Get()-it2.Get(),2)));
 329         ++it1;
 330         ++it2;
 331         ++ito;
 332       }
 333       break;
 334       */
 335     case 7: // Difference
 336       while (!ito.IsAtEnd()) {
 337         ito.Set(it1.Get()-it2.Get());
 338         ++it1;
 339         ++it2;
 340         ++ito;
 341       }
 342       break;
 343       /*
 344     case 8: // Relative Difference
 345       while (!ito.IsAtEnd()) {
 346         if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
 347         else ito.Set(0.0);
 348         ++it1;
 349         ++it2;
 350         ++ito;
 351       }
 352       break;
 353       */
 354     default: // error ?
 355       std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 356       exit(-1);
 357     }
 358   }
 359
 360 ////////////////////////////////////////////////////////////////////
 361   // specializations for itk::Vector<double, 3u>, 3u
 362   template<>
 363   template<>
 364   void ImageArithmGenericFilter<args_info_clitkImageArithm>::UpdateWithInputImageType< itk::Image< itk::Vector<double, 3u>, 3u > >()
 365   {
 366     typedef itk::Image< itk::Vector<double, 3u>, 3u > ImageType;
 367
 368     // Read input1
 369     ImageType::Pointer input1 = this->GetInput<ImageType>(0);
 370
 371     // Set input image iterator
 372     typedef itk::ImageRegionIterator<ImageType> IteratorType;
 373     IteratorType it(input1, input1->GetLargestPossibleRegion());
 374
 375     // typedef input2
 376     ImageType::Pointer input2 = NULL;
 377     IteratorType it2;
 378
 379     /*
 380     // Special case for normalisation
 381     if (mTypeOfOperation == 12) {
 382       typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
 383       typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
 384       ff->SetImage(input1);
 385       ff->ComputeMaximum();
 386       mScalar = ff->GetMaximum();
 387       mTypeOfOperation = 11; // divide
 388     }
 389     */
 390
 391     if (mIsOperationUseASecondImage) {
 392         // Read input2
 393         input2 = this->GetInput<ImageType>(1);
 394         // Set input image iterator
 395         it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
 396         // Check dimension
 397         if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
 398           itkExceptionMacro(<< "The images (input and input2) must have the same size");
 399         }
 400         if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
 401           itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
 402                           << "Using first input's information.");
 403         }
 404     }
 405
 406     /*
 407     // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
 408     if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
 409       // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
 410       //                     << typeid(PixelType).name()
 411       //                     << std::endl;
 412       mOverwriteInputImage = false;
 413     }
 414     */
 415
 416     // ---------------- Overwrite input Image ---------------------
 417     if (mOverwriteInputImage) {
 418       // Set output iterator (to input1)
 419       IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
 420       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 421       else ComputeImage(it, ito);
 422       this->SetNextOutput<ImageType>(input1);
 423     }
 424     // ---------------- Create new output Image ---------------------
 425     else {
 426       /*if (mOutputIsFloat) {
 427         // Create output image
 428         typedef itk::Image<float,ImageType::ImageDimension> OutputImageType;
 429         typename OutputImageType::Pointer output = OutputImageType::New();
 430         output->SetRegions(input1->GetLargestPossibleRegion());
 431         output->SetOrigin(input1->GetOrigin());
 432         output->SetSpacing(input1->GetSpacing());
 433         output->Allocate();
 434         // Set output iterator
 435         typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 436         IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 437         if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 438         else ComputeImage(it, ito);
 439         this->template SetNextOutput<OutputImageType>(output);
 440       } else*/ {
 441         // Create output image
 442         typedef ImageType OutputImageType;
 443         OutputImageType::Pointer output = OutputImageType::New();
 444         output->SetRegions(input1->GetLargestPossibleRegion());
 445         output->SetOrigin(input1->GetOrigin());
 446         output->SetSpacing(input1->GetSpacing());
 447         output->Allocate();
 448         // Set output iterator
 449         typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 450         IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 451         if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 452         else ComputeImage(it, ito);
 453         this->SetNextOutput<OutputImageType>(output);
 454       }
 455     }
 456   }
 457
 458   template<>
 459   template<>
 460   void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
 461     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
 462     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > >
 463     (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it,
 464     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
 465   {
 466     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter2;
 467
 468     ito.GoToBegin();
 469     it.GoToBegin();
 470
 471     typedef Iter2::PixelType PixelType;
 472
 473     PixelType scalar_vector;
 474     scalar_vector.Fill(mScalar);
 475
 476     // Perform operation
 477     switch (mTypeOfOperation) {
 478     case 0: // Addition
 479       while (!it.IsAtEnd()) {
 480         ito.Set(it.Get() + scalar_vector);
 481         ++it;
 482         ++ito;
 483       }
 484       break;
 485     case 1: // Multiply
 486       while (!it.IsAtEnd()) {
 487         ito.Set(it.Get() * mScalar);
 488         ++it;
 489         ++ito;
 490       }
 491       break;
 492       /*
 493     case 2: // Inverse
 494       while (!it.IsAtEnd()) {
 495         if (it.Get() != 0)
 496           ito.Set(mScalar / it.Get()));
 497         else ito.Set(mDefaultPixelValue);
 498         ++it;
 499         ++ito;
 500       }
 501       break;
 502     case 3: // Max
 503       while (!it.IsAtEnd()) {
 504         if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 505         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 506         ++it;
 507         ++ito;
 508       }
 509       break;
 510     case 4: // Min
 511       while (!it.IsAtEnd()) {
 512         if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 513         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 514         ++it;
 515         ++ito;
 516       }
 517       break;
 518     case 5: // Absolute value
 519       while (!it.IsAtEnd()) {
 520         if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
 521         // <= zero to avoid warning for unsigned types
 522         else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 523         ++it;
 524         ++ito;
 525       }
 526       break;
 527     case 6: // Squared value
 528       while (!it.IsAtEnd()) {
 529         ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
 530         ++it;
 531         ++ito;
 532       }
 533       break;
 534     case 7: // Log
 535       while (!it.IsAtEnd()) {
 536         if (it.Get() > 0)
 537           ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
 538         else ito.Set(mDefaultPixelValue);
 539         ++it;
 540         ++ito;
 541       }
 542       break;
 543     case 8: // exp
 544       while (!it.IsAtEnd()) {
 545         ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
 546         ++it;
 547         ++ito;
 548       }
 549       break;
 550     case 9: // sqrt
 551       while (!it.IsAtEnd()) {
 552         if (it.Get() > 0)
 553           ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
 554         else {
 555           if (it.Get() ==0) ito.Set(0);
 556           else ito.Set(mDefaultPixelValue);
 557         }
 558         ++it;
 559         ++ito;
 560       }
 561       break;
 562     case 10: // exp
 563       while (!it.IsAtEnd()) {
 564         ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
 565         ++it;
 566         ++ito;
 567       }
 568       break;
 569       */
 570     case 11: // divide
 571       while (!it.IsAtEnd()) {
 572         ito.Set(it.Get() / mScalar);
 573         ++it;
 574         ++ito;
 575       }
 576       break;
 577     default: // error ?
 578       std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 579       exit(-1);
 580     }
 581
 582   }
 583
 584   template<>
 585   template<>
 586   void ImageArithmGenericFilter<args_info_clitkImageArithm>::ComputeImage<
 587     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
 588     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >,
 589     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > >
 590     >
 591     (itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it1,
 592     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > it2,
 593     itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > ito)
 594   {
 595     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter1;
 596     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter2;
 597     typedef itk::ImageRegionIterator< itk::Image< itk::Vector<double, 3u>, 3u > > Iter3;
 598
 599     it1.GoToBegin();
 600     it2.GoToBegin();
 601     ito.GoToBegin();
 602     typedef Iter3::PixelType PixelType;
 603
 604     switch (mTypeOfOperation) {
 605     case 0: // Addition
 606       while (!ito.IsAtEnd()) {
 607         ito.Set(it1.Get() + it2.Get());
 608         ++it1;
 609         ++it2;
 610         ++ito;
 611       }
 612       break;
 613       /*
 614     case 1: // Multiply
 615       while (!ito.IsAtEnd()) {
 616         ito.Set(it1.Get() * it2.Get()) );
 617         ++it1;
 618         ++it2;
 619         ++ito;
 620       }
 621       break;
 622     case 2: // Divide
 623       while (!ito.IsAtEnd()) {
 624         if (it1.Get() != 0)
 625           ito.Set(it1.Get() / it2.Get()));
 626         else ito.Set(mDefaultPixelValue);
 627         ++it1;
 628         ++it2;
 629         ++ito;
 630       }
 631       break;
 632     case 3: // Max
 633       while (!ito.IsAtEnd()) {
 634         if (it1.Get() < it2.Get()) ito.Set(it2.Get());
 635         ++it1;
 636         ++it2;
 637         ++ito;
 638       }
 639       break;
 640     case 4: // Min
 641       while (!ito.IsAtEnd()) {
 642         if (it1.Get() > it2.Get()) ito.Set(it2.Get());
 643         ++it1;
 644         ++it2;
 645         ++ito;
 646       }
 647       break;
 648       */
 649     case 5: // Absolute difference
 650       while (!ito.IsAtEnd()) {
 651         ito.Set(it2.Get()-it1.Get());
 652         ++it1;
 653         ++it2;
 654         ++ito;
 655       }
 656       break;
 657       /*
 658     case 6: // Squared differences
 659       while (!ito.IsAtEnd()) {
 660         ito.Set(pow(it1.Get()-it2.Get(),2)));
 661         ++it1;
 662         ++it2;
 663         ++ito;
 664       }
 665       break;
 666       */
 667     case 7: // Difference
 668       while (!ito.IsAtEnd()) {
 669         ito.Set(it1.Get()-it2.Get());
 670         ++it1;
 671         ++it2;
 672         ++ito;
 673       }
 674       break;
 675       /*
 676     case 8: // Relative Difference
 677       while (!ito.IsAtEnd()) {
 678         if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
 679         else ito.Set(0.0);
 680         ++it1;
 681         ++it2;
 682         ++ito;
 683       }
 684       break;
 685       */
 686     default: // error ?
 687       std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 688       exit(-1);
 689     }
 690   }
 691
 692 }
 693
 694 #endif //define CLITKIMAGEARITHMGENERICFILTER_CXX