]> Creatis software - clitk.git/blob - tools/clitkImageUncertainty.cxx
Merge branch 'VTK6_Qt5_Overlay4D' into VTK6_Qt5_Binarize
[clitk.git] / 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 */