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 #include "clitkPortability.h"
  38
  39 // itk include
  40 #include "itkImageRegionConstIterator.h"
  41 #include "itkImageRegionIterator.h"
  42
  43 //====================================================================
  44 int main(int argc, char * argv[]) {
  45
  46   // init command line
  47   GGO(clitkImageUncertainty, args_info);
  48
  49   // Declare main types
  50   typedef float                          PixelType;
  51   const unsigned int                              Dimension=3;
  52   typedef itk::Image< PixelType, Dimension >      ImageType;
  53
  54   // Read images
  55   ImageType::Pointer input =
  56     clitk::readImage<ImageType>(args_info.input_arg, args_info.verbose_flag);
  57   ImageType::Pointer inputSquared =
  58     clitk::readImage<ImageType>(args_info.inputSquared_arg, args_info.verbose_flag);
  59
  60   // Create Output
  61   ImageType::Pointer output = ImageType::New();
  62   output->SetRegions(input->GetLargestPossibleRegion());
  63   output->CopyInformation(input);
  64   output->Allocate();
  65
  66   // Loop
  67   typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
  68   ConstIteratorType pi(input, input->GetLargestPossibleRegion());
  69   ConstIteratorType pii(inputSquared, inputSquared->GetLargestPossibleRegion());
  70   pi.GoToBegin();
  71   pii.GoToBegin();
  72   typedef itk::ImageRegionIterator<ImageType> IteratorType;
  73   IteratorType po(output, output->GetLargestPossibleRegion());
  74   po.GoToBegin();
  75
  76   int NumberOfEvents = args_info.NumberOfEvents_arg;
  77   while ( !pi.IsAtEnd() ) {
  78     double squared = pii.Get();
  79     double mean = pi.Get();
  80     double uncert = sqrt((NumberOfEvents*squared - mean*mean) / (NumberOfEvents-1));
  81     if(!args_info.absolute_flag)
  82       uncert /= std::abs(mean);
  83     if (!IsNormal(uncert)) uncert = args_info.default_arg;
  84     po.Set(uncert);
  85     ++pi;
  86     ++pii;
  87     ++po;
  88   }
  89
  90   // Write output image
  91   // DD(clitk::GetExtension(args_info.output_arg));
  92   if (clitk::GetExtension(args_info.output_arg) != "txt") {
  93     clitk::writeImage<ImageType>(output, args_info.output_arg, args_info.verbose_flag);
  94   }
  95   else {
  96     std::ofstream os;
  97     clitk::openFileForWriting(os, args_info.output_arg);
  98     typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
  99     IteratorType pi(output, output->GetLargestPossibleRegion());
 100     pi.GoToBegin();
 101     os << "# Image size = " << output->GetLargestPossibleRegion().GetSize() << std::endl;
 102     os << "# Image spacing = " << output->GetSpacing() << std::endl;
 103     os << "# Number of events = " << NumberOfEvents << std::endl;
 104     while (!pi.IsAtEnd()) {
 105       os << pi.Get() << std::endl;
 106       ++pi;
 107     }
 108   }
 109 }
 110
 111
 112 #endif /* end #define CLITKIMAGEUNCERTAINTY_CXX */