1 /*=========================================================================
4 Module: $RCSfile: clitkImageUncertainty.cxx,v $
6 Date: $Date: 2011/03/03 15:03:30 $
7 Version: $Revision: 1.3 $
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.
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.
17 =========================================================================*/
19 #ifndef CLITKIMAGEUNCERTAINTY_CXX
20 #define CLITKIMAGEUNCERTAINTY_CXX
23 =================================================
24 * @file clitkImageUncertainty.cxx
25 * @author David Sarrut <david.sarrut@creatis.insa-lyon.fr>
26 * @date 04 Jul 2006 14:03:57
31 =================================================*/
34 #include "clitkImageUncertainty_ggo.h"
35 #include "clitkImageCommon.h"
36 #include "clitkCommon.h"
39 #include "itkImageRegionConstIterator.h"
40 #include "itkImageRegionIterator.h"
42 #include <cmath> // for isfinite
44 //====================================================================
45 int main(int argc, char * argv[]) {
48 GGO(clitkImageUncertainty, args_info);
51 typedef float PixelType;
52 const unsigned int Dimension=3;
53 typedef itk::Image< PixelType, Dimension > ImageType;
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);
62 ImageType::Pointer output = ImageType::New();
63 output->SetRegions(input->GetLargestPossibleRegion());
64 output->CopyInformation(input);
68 typedef itk::ImageRegionConstIterator<ImageType> ConstIteratorType;
69 ConstIteratorType pi(input, input->GetLargestPossibleRegion());
70 ConstIteratorType pii(inputSquared, inputSquared->GetLargestPossibleRegion());
73 typedef itk::ImageRegionIterator<ImageType> IteratorType;
74 IteratorType po(output, output->GetLargestPossibleRegion());
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.;
88 // *po = sqrt( (NumberOfEvents*squared - mean*mean) /
89 // ((NumberOfEvents-1)*(mean*mean)) );
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);
101 clitk::openFileForWriting(os, args_info.output_arg);
102 typedef itk::ImageRegionConstIterator<ImageType> IteratorType;
103 IteratorType pi(output, output->GetLargestPossibleRegion());
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;
116 #endif /* end #define CLITKIMAGEUNCERTAINTY_CXX */