tools/clitkVectorArithmGenericFilter.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 CLITKVectorArithmGENERICFILTER_TXX
  19 #define CLITKVectorArithmGENERICFILTER_TXX
  20
  21 #include "clitkImageCommon.h"
  22
  23 #include "itkMinimumMaximumImageCalculator.h"
  24
  25 namespace clitk
  26 {
  27
  28 //--------------------------------------------------------------------
  29 template<class args_info_type>
  30 VectorArithmGenericFilter<args_info_type>::VectorArithmGenericFilter()
  31   :ImageToImageGenericFilter<Self>("VectorArithmGenericFilter"),mTypeOfOperation(0)
  32 {
  33   InitializeImageType<3>();
  34   mIsOperationUseASecondImage = false;
  35   mIsOutputScalar = false;
  36   mOverwriteInputImage = true;
  37 }
  38 //--------------------------------------------------------------------
  39
  40
  41 //--------------------------------------------------------------------
  42 template<class args_info_type>
  43 template<unsigned int Dim>
  44 void VectorArithmGenericFilter<args_info_type>::InitializeImageType()
  45 {
  46   ADD_VEC_IMAGE_TYPE(Dim,3u,float);
  47   ADD_VEC_IMAGE_TYPE(Dim,3u,double);
  48   ADD_VEC_IMAGE_TYPE(Dim,2u,float);
  49   ADD_VEC_IMAGE_TYPE(Dim,2u,double);
  50 }
  51 //--------------------------------------------------------------------
  52
  53
  54 //--------------------------------------------------------------------
  55 template<class args_info_type>
  56 void VectorArithmGenericFilter<args_info_type>::EnableOverwriteInputImage(bool b)
  57 {
  58   mOverwriteInputImage = b;
  59 }
  60 //--------------------------------------------------------------------
  61
  62
  63 //--------------------------------------------------------------------
  64 template<class args_info_type>
  65 void VectorArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type & a)
  66 {
  67   mArgsInfo=a;
  68
  69   // Set value
  70   this->SetIOVerbose(mArgsInfo.verbose_flag);
  71   mTypeOfOperation = mArgsInfo.operation_arg;
  72   mDefaultPixelValue = mArgsInfo.pixelValue_arg;
  73   mScalar = mArgsInfo.scalar_arg;
  74   mOutputIsFloat = mArgsInfo.setFloatOutput_flag;
  75
  76   if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
  77
  78   if (mArgsInfo.input1_given) this->AddInputFilename(mArgsInfo.input1_arg);
  79   if (mArgsInfo.input2_given) {
  80     mIsOperationUseASecondImage = true;
  81     this->AddInputFilename(mArgsInfo.input2_arg);
  82     if (mArgsInfo.operation_arg == 1)
  83       mIsOutputScalar = true;
  84   }
  85   else if (mArgsInfo.operation_arg == 5 || mArgsInfo.operation_arg == 6)
  86     mIsOutputScalar = true;
  87
  88   if (mArgsInfo.output_given) this->SetOutputFilename(mArgsInfo.output_arg);
  89
  90   // Check type of operation (with scalar or with other image)
  91   if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
  92     std::cerr << "ERROR : you cannot provide both --scalar and --input2 option" << std::endl;
  93     exit(-1);
  94   }
  95   if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
  96     if (mArgsInfo.operation_arg < 5) {
  97       std::cerr << "Such operation need the --scalar option." << std::endl;
  98       exit(-1);
  99     }
 100   }
 101 }
 102 //--------------------------------------------------------------------
 103
 104
 105 //--------------------------------------------------------------------
 106 template<class args_info_type>
 107 template<class ImageType>
 108 void VectorArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
 109 {
 110   // Read input1
 111   typename ImageType::Pointer input1 = this->template GetInput<ImageType>(0);
 112
 113   // Set input image iterator
 114   typedef itk::ImageRegionIterator<ImageType> IteratorType;
 115   IteratorType it(input1, input1->GetLargestPossibleRegion());
 116
 117   // typedef input2
 118   typename ImageType::Pointer input2 = NULL;
 119   IteratorType it2;
 120
 121   // Special case for normalisation
 122   /*
 123   if (mTypeOfOperation == 12) {
 124     typedef itk::MinimumMaximumImageCalculator<ImageType> MinMaxFilterType;
 125     typename MinMaxFilterType::Pointer ff = MinMaxFilterType::New();
 126     ff->SetImage(input1);
 127     ff->ComputeMaximum();
 128     mScalar = ff->GetMaximum();
 129     mTypeOfOperation = 11; // divide
 130   }
 131   */
 132
 133   if (mIsOperationUseASecondImage) {
 134       // Read input2
 135       input2 = this->template GetInput<ImageType>(1);
 136       // Set input image iterator
 137       it2 = IteratorType(input2, input2->GetLargestPossibleRegion());
 138       // Check dimension
 139       if (!clitk::HaveSameSize<ImageType, ImageType>(input1, input2)) {
 140         itkExceptionMacro(<< "The images (input and input2) must have the same size");
 141       }
 142       if(!clitk::HaveSameSpacing<ImageType, ImageType>(input1, input2)) {
 143         itkWarningMacro(<< "The images (input and input2) do not have the same spacing. "
 144                         << "Using first input's information.");
 145       }
 146   }
 147
 148   // Check if overwrite and outputisfloat and pixeltype is not float -> do not overwrite
 149   if (mOverwriteInputImage && mOutputIsFloat && (typeid(typename ImageType::PixelType) != typeid(float))) {
 150     // std::cerr << "Warning. Could not use both mOverwriteInputImage and mOutputIsFloat, because input is "
 151     //                     << typeid(PixelType).name()
 152     //                     << std::endl;
 153     mOverwriteInputImage = false;
 154   }
 155
 156   // ---------------- Overwrite input Image ---------------------
 157   if (mOverwriteInputImage && !mIsOutputScalar) {
 158     // Set output iterator (to input1)
 159     IteratorType ito = IteratorType(input1, input1->GetLargestPossibleRegion());
 160     if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 161     else ComputeImage(it, ito);
 162     this->template SetNextOutput<ImageType>(input1);
 163   }
 164
 165   // ---------------- Create new output Image ---------------------
 166   else {
 167     // Create output image
 168     if (!mIsOutputScalar) {
 169       typedef ImageType OutputImageType;
 170       typename OutputImageType::Pointer output = OutputImageType::New();
 171       output->SetRegions(input1->GetLargestPossibleRegion());
 172       output->SetOrigin(input1->GetOrigin());
 173       output->SetSpacing(input1->GetSpacing());
 174       output->Allocate();
 175       // Set output iterator
 176       typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 177       IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 178       if (mIsOperationUseASecondImage) ComputeImage(it, it2, ito);
 179       else ComputeImage(it, ito);
 180       this->template SetNextOutput<OutputImageType>(output);
 181     }
 182     else {
 183       // Create scalar output image
 184       typedef itk::Image<typename ImageType::PixelType::ValueType, ImageType::ImageDimension> OutputImageType;
 185       typename OutputImageType::Pointer output = OutputImageType::New();
 186       output->SetRegions(input1->GetLargestPossibleRegion());
 187       output->SetOrigin(input1->GetOrigin());
 188       output->SetSpacing(input1->GetSpacing());
 189       output->Allocate();
 190       // Set output iterator
 191       typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
 192       IteratorOutputType ito = IteratorOutputType(output, output->GetLargestPossibleRegion());
 193       if (mIsOperationUseASecondImage) ComputeScalarImage(it, it2, ito);
 194       else ComputeScalarImage(it, ito);
 195       this->template SetNextOutput<OutputImageType>(output);
 196     }
 197   }
 198 }
 199 //--------------------------------------------------------------------
 200
 201 //--------------------------------------------------------------------
 202 template<class args_info_type>
 203 template<class Iter1, class Iter2, class Iter3>
 204 void  VectorArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it2, Iter3 ito)
 205 {
 206   it1.GoToBegin();
 207   it2.GoToBegin();
 208   ito.GoToBegin();
 209   typedef typename Iter3::PixelType PixelType;
 210
 211   switch (mTypeOfOperation) {
 212   case 0: // Addition
 213     while (!ito.IsAtEnd()) {
 214       ito.Set(it1.Get() + it2.Get());
 215       ++it1;
 216       ++it2;
 217       ++ito;
 218     }
 219     break;
 220     /*
 221   case 1: // Multiply
 222     while (!ito.IsAtEnd()) {
 223       ito.Set(it1.Get() * it2.Get()) );
 224       ++it1;
 225       ++it2;
 226       ++ito;
 227     }
 228     break;
 229     */
 230     /*
 231   case 2: // Divide
 232     while (!ito.IsAtEnd()) {
 233       if (it1.Get() != 0)
 234         ito.Set(it1.Get() / it2.Get()));
 235       else ito.Set(mDefaultPixelValue);
 236       ++it1;
 237       ++it2;
 238       ++ito;
 239     }
 240     break;
 241   case 3: // Max
 242     while (!ito.IsAtEnd()) {
 243       if (it1.Get() < it2.Get()) ito.Set(it2.Get());
 244       ++it1;
 245       ++it2;
 246       ++ito;
 247     }
 248     break;
 249   case 4: // Min
 250     while (!ito.IsAtEnd()) {
 251       if (it1.Get() > it2.Get()) ito.Set(it2.Get());
 252       ++it1;
 253       ++it2;
 254       ++ito;
 255     }
 256     break;
 257     */
 258     /*
 259   case 5: // Absolute difference
 260     while (!ito.IsAtEnd()) {
 261       ito.Set(it2.Get()-it1.Get());
 262       ++it1;
 263       ++it2;
 264       ++ito;
 265     }
 266     break;
 267     */
 268     /*
 269   case 6: // Squared differences
 270     while (!ito.IsAtEnd()) {
 271       ito.Set(pow(it1.Get()-it2.Get(),2)));
 272       ++it1;
 273       ++it2;
 274       ++ito;
 275     }
 276     break;
 277     */
 278   case 7: // Difference
 279     while (!ito.IsAtEnd()) {
 280       ito.Set(it1.Get()-it2.Get());
 281       ++it1;
 282       ++it2;
 283       ++ito;
 284     }
 285     break;
 286     /*
 287   case 8: // Relative Difference
 288     while (!ito.IsAtEnd()) {
 289       if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
 290       else ito.Set(0.0);
 291       ++it1;
 292       ++it2;
 293       ++ito;
 294     }
 295     break;
 296     */
 297     /*
 298   case 9: // CrossProduct
 299     while (!ito.IsAtEnd()) {
 300       ito.Set(it1.Get()^it2.Get());
 301       ++it1;
 302       ++it2;
 303       ++ito;
 304     }
 305     break;
 306     */
 307   default: // error ?
 308     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 309     exit(-1);
 310   }
 311 }
 312 //--------------------------------------------------------------------
 313
 314
 315 //--------------------------------------------------------------------
 316 template<class args_info_type>
 317 template<class Iter1, class Iter2>
 318 void clitk::VectorArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Iter2 ito)
 319 {
 320   ito.GoToBegin();
 321   it.GoToBegin();
 322
 323   typedef typename Iter1::PixelType PixelType;
 324
 325   PixelType scalar_vector;
 326   scalar_vector.Fill(mScalar);
 327
 328   // Perform operation
 329   switch (mTypeOfOperation) {
 330   case 0: // Addition
 331     while (!it.IsAtEnd()) {
 332       ito.Set(it.Get() + scalar_vector);
 333       ++it;
 334       ++ito;
 335     }
 336     break;
 337   case 1: // Multiply
 338     while (!it.IsAtEnd()) {
 339       ito.Set(it.Get() * mScalar);
 340       ++it;
 341       ++ito;
 342     }
 343     break;
 344     /*
 345   case 2: // Inverse
 346     while (!it.IsAtEnd()) {
 347       if (it.Get() != 0)
 348         ito.Set(mScalar / it.Get()));
 349       else ito.Set(mDefaultPixelValue);
 350       ++it;
 351       ++ito;
 352     }
 353     break;
 354   case 3: // Max
 355     while (!it.IsAtEnd()) {
 356       if (it.Get() < mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 357       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 358       ++it;
 359       ++ito;
 360     }
 361     break;
 362   case 4: // Min
 363     while (!it.IsAtEnd()) {
 364       if (it.Get() > mScalar) ito.Set(PixelTypeDownCast<double, PixelType>(mScalar));
 365       else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
 366       ++it;
 367       ++ito;
 368     }
 369     break;
 370     */
 371     /*
 372   case 5: // Absolute value
 373     while (!it.IsAtEnd()) {
 374       ito.Set(PixelTypeDownCast<double, PixelType>(it.GetNorm()));
 375       ++it;
 376       ++ito;
 377     }
 378     break;
 379   case 6: // Squared value
 380     while (!it.IsAtEnd()) {
 381       ito.Set(PixelTypeDownCast<double, PixelType>(it.GetSquaredNorm());
 382       ++it;
 383       ++ito;
 384     }
 385     break;
 386     */
 387     /*
 388   case 7: // Log
 389     while (!it.IsAtEnd()) {
 390       if (it.Get() > 0)
 391         ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
 392       else ito.Set(mDefaultPixelValue);
 393       ++it;
 394       ++ito;
 395     }
 396     break;
 397   case 8: // exp
 398     while (!it.IsAtEnd()) {
 399       ito.Set(PixelTypeDownCast<double, PixelType>(exp((double)it.Get())));
 400       ++it;
 401       ++ito;
 402     }
 403     break;
 404   case 9: // sqrt
 405     while (!it.IsAtEnd()) {
 406       if (it.Get() > 0)
 407         ito.Set(PixelTypeDownCast<double, PixelType>(sqrt((double)it.Get())));
 408       else {
 409         if (it.Get() ==0) ito.Set(0);
 410         else ito.Set(mDefaultPixelValue);
 411       }
 412       ++it;
 413       ++ito;
 414     }
 415     break;
 416   case 10: // exp
 417     while (!it.IsAtEnd()) {
 418       ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
 419       ++it;
 420       ++ito;
 421     }
 422     break;
 423     */
 424   case 11: // divide
 425     while (!it.IsAtEnd()) {
 426       ito.Set(it.Get() / mScalar);
 427       ++it;
 428       ++ito;
 429     }
 430     break;
 431   case 12: // normalize
 432     while (!it.IsAtEnd()) {
 433       PixelType n = it.Get();
 434       if (n.GetNorm() != 0)
 435         n.Normalize();
 436
 437       ito.Set(n);
 438       ++it;
 439       ++ito;
 440     }
 441     break;
 442   default: // error ?
 443     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 444     exit(-1);
 445   }
 446 }
 447 //--------------------------------------------------------------------
 448
 449 //--------------------------------------------------------------------
 450 template<class args_info_type>
 451 template<class Iter1, class Iter2, class Iter3>
 452 void  VectorArithmGenericFilter<args_info_type>::ComputeScalarImage(Iter1 it1, Iter2 it2, Iter3 ito)
 453 {
 454   it1.GoToBegin();
 455   it2.GoToBegin();
 456   ito.GoToBegin();
 457   typedef typename Iter3::PixelType PixelType;
 458
 459   switch (mTypeOfOperation) {
 460   case 1: // Multiply
 461     while (!ito.IsAtEnd()) {
 462       ito.Set(it1.Get() * it2.Get());
 463       ++it1;
 464       ++it2;
 465       ++ito;
 466     }
 467     break;
 468   default: // error ?
 469     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 470     exit(-1);
 471   }
 472 }
 473 //--------------------------------------------------------------------
 474
 475 //--------------------------------------------------------------------
 476 template<class args_info_type>
 477 template<class Iter1, class Iter2>
 478 void clitk::VectorArithmGenericFilter<args_info_type>::ComputeScalarImage(Iter1 it, Iter2 ito)
 479 {
 480   ito.GoToBegin();
 481   it.GoToBegin();
 482
 483   typedef typename Iter2::PixelType PixelType;
 484
 485
 486   // Perform operation
 487   switch (mTypeOfOperation) {
 488   case 5: // Absolute value
 489     while (!it.IsAtEnd()) {
 490       ito.Set(PixelTypeDownCast<double, PixelType>(it.Get().GetNorm()));
 491       ++it;
 492       ++ito;
 493     }
 494     break;
 495   case 6: // Squared value
 496     while (!it.IsAtEnd()) {
 497       ito.Set(PixelTypeDownCast<double, PixelType>(it.Get().GetSquaredNorm()));
 498       ++it;
 499       ++ito;
 500     }
 501     break;
 502   default: // error ?
 503     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
 504     exit(-1);
 505   }
 506 }
 507 //--------------------------------------------------------------------
 508
 509
 510 } // end namespace
 511
 512 #endif  //#define CLITKVectorArithmGENERICFILTER_TXX