]> Creatis software - clitk.git/blobdiff - tools/clitkImageUncertainty.cxx
Debug RTStruct conversion with empty struc
[clitk.git] / tools / clitkImageUncertainty.cxx
index 6fc003b959710f4afd1214a27f18faa88fd8f2d6..a8e76de080c278e6d87cb095a1f2ee03160a808f 100644 (file)
@@ -1,19 +1,19 @@
 /*=========================================================================
-                                                                                
+
   Program:   clitk
   Module:    $RCSfile: clitkImageUncertainty.cxx,v $
   Language:  C++
   Date:      $Date: 2011/03/03 15:03:30 $
   Version:   $Revision: 1.3 $
-                                                                                
+
   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
   l'Image). All rights reserved. See Doc/License.txt or
   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
-                                                                                
+
      This software is distributed WITHOUT ANY WARRANTY; without even
      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
      PURPOSE.  See the above copyright notices for more information.
-                                                                             
+
 =========================================================================*/
 
 #ifndef CLITKIMAGEUNCERTAINTY_CXX
  * @file   clitkImageUncertainty.cxx
  * @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
  * @date   04 Jul 2006 14:03:57
- * 
- * @brief  
- * 
- * 
+ *
+ * @brief
+ *
+ *
  =================================================*/
 
 // clitk include
 #include "clitkImageUncertainty_ggo.h"
 #include "clitkImageCommon.h"
 #include "clitkCommon.h"
+#include "clitkPortability.h"
 
 // itk include
 #include "itkImageRegionConstIterator.h"
 #include "itkImageRegionIterator.h"
 
-#include <cmath> // for isfinite
-
 //====================================================================
 int main(int argc, char * argv[]) {
 
@@ -51,67 +50,63 @@ int main(int argc, char * argv[]) {
   typedef float                          PixelType;
   const unsigned int                              Dimension=3;
   typedef itk::Image< PixelType, Dimension >      ImageType;
-  
+
   // Read images
-  ImageType::Pointer input = 
-       clitk::readImage<ImageType>(args_info.input_arg, args_info.verbose_flag);
-  ImageType::Pointer inputSquared = 
-       clitk::readImage<ImageType>(args_info.inputSquared_arg, args_info.verbose_flag);
+  ImageType::Pointer input =
+    clitk::readImage<ImageType>(args_info.input_arg, args_info.verbose_flag);
+  ImageType::Pointer inputSquared =
+    clitk::readImage<ImageType>(args_info.inputSquared_arg, args_info.verbose_flag);
 
   // Create Output
   ImageType::Pointer output = ImageType::New();
   output->SetRegions(input->GetLargestPossibleRegion());
   output->CopyInformation(input);
-  output->Allocate();  
+  output->Allocate();
 
   // Loop
   typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
   ConstIteratorType pi(input, input->GetLargestPossibleRegion());
   ConstIteratorType pii(inputSquared, inputSquared->GetLargestPossibleRegion());
-  pi.Begin();
-  pii.Begin();
+  pi.GoToBegin();
+  pii.GoToBegin();
   typedef itk::ImageRegionIterator<ImageType> IteratorType;
   IteratorType po(output, output->GetLargestPossibleRegion());
-  po.Begin();
-  
+  po.GoToBegin();
+
   int NumberOfEvents = args_info.NumberOfEvents_arg;
-  while ( !pi.IsAtEnd() ) {  
-       double squared = pii.Get();
-       double mean = pi.Get();
-       double uncert = sqrt((NumberOfEvents*squared - mean*mean) / ((NumberOfEvents-1)*(mean*mean)));
-       if (!std::isnormal(uncert)) uncert = 1.;
-       po.Set(uncert);
-       ++pi;
-       ++pii;
-       ++po;
+  while ( !pi.IsAtEnd() ) {
+    double squared = pii.Get();
+    double mean = pi.Get();
+    double uncert = sqrt((NumberOfEvents*squared - mean*mean) / (NumberOfEvents-1));
+    if(!args_info.absolute_flag)
+      uncert /= std::abs(mean);
+    if (!IsNormal(uncert)) uncert = args_info.default_arg;
+    po.Set(uncert);
+    ++pi;
+    ++pii;
+    ++po;
   }
-//       *po = sqrt( (NumberOfEvents*squared - mean*mean) / 
-// ((NumberOfEvents-1)*(mean*mean)) );
-//       ++po;
-
-
 
   // Write output image
   // DD(clitk::GetExtension(args_info.output_arg));
   if (clitk::GetExtension(args_info.output_arg) != "txt") {
-       clitk::writeImage<ImageType>(output, args_info.output_arg, args_info.verbose_flag);
+    clitk::writeImage<ImageType>(output, args_info.output_arg, args_info.verbose_flag);
   }
   else {
-       std::ofstream os;                                                                                                       
-       clitk::openFileForWriting(os, args_info.output_arg);                            
-       typedef itk::ImageRegionConstIterator<ImageType> IteratorType; 
-       IteratorType pi(output, output->GetLargestPossibleRegion());            
-       pi.Begin();                                                                                                             
-       os << "# Image size = " << output->GetLargestPossibleRegion().GetSize() << std::endl; 
-       os << "# Image spacing = " << output->GetSpacing() << std::endl;
-       os << "# Number of events = " << NumberOfEvents << std::endl;
-       while (!pi.IsAtEnd()) {                                                                                 
-         os << pi.Get() << std::endl;                                                                  
-         ++pi;                                                                                                                 
-       }                       
+    std::ofstream os;
+    clitk::openFileForWriting(os, args_info.output_arg);
+    typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
+    IteratorType pi(output, output->GetLargestPossibleRegion());
+    pi.GoToBegin();
+    os << "# Image size = " << output->GetLargestPossibleRegion().GetSize() << std::endl;
+    os << "# Image spacing = " << output->GetSpacing() << std::endl;
+    os << "# Number of events = " << NumberOfEvents << std::endl;
+    while (!pi.IsAtEnd()) {
+      os << pi.Get() << std::endl;
+      ++pi;
+    }
   }
 }
 
 
 #endif /* end #define CLITKIMAGEUNCERTAINTY_CXX */
-