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