tools/clitkImageUncertainty.cxx

   1 /*=========================================================================
   2
   3   Program:   clitk
   4   Module:    $RCSfile: clitkImageUncertainty.cxx,v $
   5   Language:  C++
   6   Date:      $Date: 2011/03/03 15:03:30 $
   7   Version:   $Revision: 1.3 $
   8
   9   Copyright (c) CREATIS (Centre de Recherche et d'Applications en Traitement de
  10   l'Image). All rights reserved. See Doc/License.txt or
  11   http://www.creatis.insa-lyon.fr/Public/Gdcm/License.html for details.
  12
  13      This software is distributed WITHOUT ANY WARRANTY; without even
  14      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  15      PURPOSE.  See the above copyright notices for more information.
  16
  17 =========================================================================*/
  18
  19 #ifndef CLITKIMAGEUNCERTAINTY_CXX
  20 #define CLITKIMAGEUNCERTAINTY_CXX
  21
  22 /**
  23  =================================================
  24  * @file   clitkImageUncertainty.cxx
  25  * @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
  26  * @date   04 Jul 2006 14:03:57
  27  *
  28  * @brief
  29  *
  30  *
  31  =================================================*/
  32
  33 // clitk include
  34 #include "clitkImageUncertainty_ggo.h"
  35 #include "clitkImageCommon.h"
  36 #include "clitkCommon.h"
  37
  38 // itk include
  39 #include "itkImageRegionConstIterator.h"
  40 #include "itkImageRegionIterator.h"
  41
  42 #include <cmath> // for isfinite
  43
  44 //====================================================================
  45 int main(int argc, char * argv[]) {
  46
  47   // init command line
  48   GGO(clitkImageUncertainty, args_info);
  49
  50   // Declare main types
  51   typedef float                          PixelType;
  52   const unsigned int                              Dimension=3;
  53   typedef itk::Image< PixelType, Dimension >      ImageType;
  54
  55   // Read images
  56   ImageType::Pointer input =
  57         clitk::readImage<ImageType>(args_info.input_arg, args_info.verbose_flag);
  58   ImageType::Pointer inputSquared =
  59         clitk::readImage<ImageType>(args_info.inputSquared_arg, args_info.verbose_flag);
  60
  61   // Create Output
  62   ImageType::Pointer output = ImageType::New();
  63   output->SetRegions(input->GetLargestPossibleRegion());
  64   output->CopyInformation(input);
  65   output->Allocate();
  66
  67   // Loop
  68   typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
  69   ConstIteratorType pi(input, input->GetLargestPossibleRegion());
  70   ConstIteratorType pii(inputSquared, inputSquared->GetLargestPossibleRegion());
  71   pi.Begin();
  72   pii.Begin();
  73   typedef itk::ImageRegionIterator<ImageType> IteratorType;
  74   IteratorType po(output, output->GetLargestPossibleRegion());
  75   po.Begin();
  76
  77   int NumberOfEvents = args_info.NumberOfEvents_arg;
  78   while ( !pi.IsAtEnd() ) {
  79         double squared = pii.Get();
  80         double mean = pi.Get();
  81         double uncert = sqrt((NumberOfEvents*squared - mean*mean) / ((NumberOfEvents-1)*(mean*mean)));
  82         if (!std::isnormal(uncert)) uncert = 1.;
  83         po.Set(uncert);
  84         ++pi;
  85         ++pii;
  86         ++po;
  87   }
  88 //        *po = sqrt( (NumberOfEvents*squared - mean*mean) /
  89 // ((NumberOfEvents-1)*(mean*mean)) );
  90 //        ++po;
  91
  92
  93
  94   // Write output image
  95   // DD(clitk::GetExtension(args_info.output_arg));
  96   if (clitk::GetExtension(args_info.output_arg) != "txt") {
  97         clitk::writeImage<ImageType>(output, args_info.output_arg, args_info.verbose_flag);
  98   }
  99   else {
 100         std::ofstream os;
 101         clitk::openFileForWriting(os, args_info.output_arg);
 102         typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
 103         IteratorType pi(output, output->GetLargestPossibleRegion());
 104         pi.Begin();
 105         os << "# Image size = " << output->GetLargestPossibleRegion().GetSize() << std::endl;
 106         os << "# Image spacing = " << output->GetSpacing() << std::endl;
 107         os << "# Number of events = " << NumberOfEvents << std::endl;
 108         while (!pi.IsAtEnd()) {
 109           os << pi.Get() << std::endl;
 110           ++pi;
 111         }
 112   }
 113 }
 114
 115
 116 #endif /* end #define CLITKIMAGEUNCERTAINTY_CXX */
 117