Add scalar into clitkImageArithm operation

[clitk.git] / tools / clitkImageArithmGenericFilter.txx

diff --git a/tools/clitkImageArithmGenericFilter.txx b/tools/clitkImageArithmGenericFilter.txx
index 24ebbbdac8d2d472ec6ddca3eebbdcfdd5f33f6f..6bf79d516d92698c09bd82c29125de065b84dae7 100644 (file)
--- a/tools/clitkImageArithmGenericFilter.txx
+++ b/tools/clitkImageArithmGenericFilter.txx
@@ -65,7 +65,7 @@ void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type
   mArgsInfo=a;
 
   // Set value
-  SetIOVerbose(mArgsInfo.verbose_flag);
+  this->SetIOVerbose(mArgsInfo.verbose_flag);
   mTypeOfOperation = mArgsInfo.operation_arg;
   mDefaultPixelValue = mArgsInfo.pixelValue_arg;
   mScalar = mArgsInfo.scalar_arg;
@@ -73,13 +73,13 @@ void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type
 
   if (mArgsInfo.imagetypes_flag) this->PrintAvailableImageTypes();
 
-  if (mArgsInfo.input1_given) AddInputFilename(mArgsInfo.input1_arg);
+  if (mArgsInfo.input1_given) this->AddInputFilename(mArgsInfo.input1_arg);
   if (mArgsInfo.input2_given) {
     mIsOperationUseASecondImage = true;
-    AddInputFilename(mArgsInfo.input2_arg);
+    this->AddInputFilename(mArgsInfo.input2_arg);
   }
 
-  if (mArgsInfo.output_given) SetOutputFilename(mArgsInfo.output_arg);
+  if (mArgsInfo.output_given) this->SetOutputFilename(mArgsInfo.output_arg);
 
   // Check type of operation (with scalar or with other image)
   if ((mArgsInfo.input2_given) && (mArgsInfo.scalar_given)) {
@@ -87,7 +87,7 @@ void ImageArithmGenericFilter<args_info_type>::SetArgsInfo(const args_info_type
     exit(-1);
   }
   if ((!mArgsInfo.input2_given) && (!mArgsInfo.scalar_given)) {
-    if (mArgsInfo.operation_arg < 5) {
+    if (mArgsInfo.operation_arg < 7 || mArgsInfo.operation_arg == 10 || mArgsInfo.operation_arg == 11 || mArgsInfo.operation_arg == 13) {
       std::cerr << "Such operation need the --scalar option." << std::endl;
       exit(-1);
     }
@@ -109,7 +109,7 @@ void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
   IteratorType it(input1, input1->GetLargestPossibleRegion());
 
   // typedef input2
-  typename ImageType::Pointer input2 = NULL;
+  typename ImageType::Pointer input2 = ITK_NULLPTR;
   IteratorType it2;
 
   // Special case for normalisation
@@ -163,6 +163,7 @@ void ImageArithmGenericFilter<args_info_type>::UpdateWithInputImageType()
       output->SetRegions(input1->GetLargestPossibleRegion());
       output->SetOrigin(input1->GetOrigin());
       output->SetSpacing(input1->GetSpacing());
+      output->SetDirection(input1->GetDirection());
       output->Allocate();
       // Set output iterator
       typedef itk::ImageRegionIterator<OutputImageType> IteratorOutputType;
@@ -218,7 +219,7 @@ void  ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it
     break;
   case 2: // Divide
     while (!ito.IsAtEnd()) {
-      if (it1.Get() != 0)
+      if (it2.Get() != 0)
         ito.Set(PixelTypeDownCast<double, PixelType>((double)it1.Get() / (double)it2.Get()));
       else ito.Set(mDefaultPixelValue);
       ++it1;
@@ -269,7 +270,7 @@ void  ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it1, Iter2 it
   case 8: // Relative Difference
     while (!ito.IsAtEnd()) {
       if (it1.Get() != 0) ito.Set(PixelTypeDownCast<double, PixelType>(((double)it1.Get()-(double)it2.Get()))/(double)it1.Get());
-      else ito.Set(0.0);
+      else ito.Set(mDefaultPixelValue);
       ++it1;
       ++it2;
       ++ito;
@@ -335,21 +336,22 @@ void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Ite
     break;
   case 5: // Absolute value
     while (!it.IsAtEnd()) {
-      if (it.Get() <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-it.Get()));
+      double value = (double)it.Get() - mScalar;
+      if (value <= 0) ito.Set(PixelTypeDownCast<double, PixelType>(-value));
       // <= zero to avoid warning for unsigned types
-      else ito.Set(PixelTypeDownCast<double, PixelType>(it.Get()));
+      else ito.Set(value);
       ++it;
       ++ito;
     }
     break;
   case 6: // Squared value
     while (!it.IsAtEnd()) {
-      ito.Set(PixelTypeDownCast<double, PixelType>((double)it.Get()*(double)it.Get()));
+      ito.Set(PixelTypeDownCast<double, PixelType>(((double)it.Get()-mScalar)*((double)it.Get()-mScalar)));
       ++it;
       ++ito;
     }
     break;
-  case 7: // Log
+  case 7: // ln
     while (!it.IsAtEnd()) {
       if (it.Get() > 0)
         ito.Set(PixelTypeDownCast<double, PixelType>(log((double)it.Get())));
@@ -377,7 +379,7 @@ void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Ite
       ++ito;
     }
     break;
-  case 10: // exp
+  case 10: // EPID
     while (!it.IsAtEnd()) {
       ito.Set(PixelTypeDownCast<double, PixelType>((0x10000 - (double)it.Get())/mScalar));
       ++it;
@@ -391,6 +393,18 @@ void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Ite
       ++ito;
     }
     break;
+  case 13: // -ln I/I0
+    while (!it.IsAtEnd()) {
+      if (it.Get() == 0) { // special case for fluence image with 0 value in a pixel -> consider 0.5
+        ito.Set(-log(0.5 / mScalar) );
+      }
+      else {
+        ito.Set(PixelTypeDownCast<double, PixelType>(-log((double)it.Get() / mScalar)) );
+      }
+      ++it;
+      ++ito;
+    }
+    break;
   default: // error ?
     std::cerr << "ERROR : the operation number (" << mTypeOfOperation << ") is not known." << std::endl;
     exit(-1);
@@ -398,6 +412,8 @@ void clitk::ImageArithmGenericFilter<args_info_type>::ComputeImage(Iter1 it, Ite
 }
 //--------------------------------------------------------------------
 
+
+
 } // end namespace
 
 #endif  //#define CLITKIMAGEARITHMGENERICFILTER_TXX