Orfeo Toolbox  3.16
itkStatisticsImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Insight Segmentation & Registration Toolkit
4  Module: $RCSfile: itkStatisticsImageFilter.txx,v $
5  Language: C++
6  Date: $Date: 2008-10-24 17:44:15 $
7  Version: $Revision: 1.21 $
8 
9  Copyright (c) Insight Software Consortium. All rights reserved.
10  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
11 
12  This software is distributed WITHOUT ANY WARRANTY; without even
13  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
14  PURPOSE. See the above copyright notices for more information.
15 
16 =========================================================================*/
17 #ifndef __itkStatisticsImageFilter_txx
18 #define __itkStatisticsImageFilter_txx
19 #include "itkStatisticsImageFilter.h"
20 
21 #include "itkImageRegionIterator.h"
22 #include "itkImageRegionConstIterator.h"
23 #include "itkNumericTraits.h"
24 #include "itkProgressReporter.h"
25 
26 namespace itk {
27 
28 template<class TInputImage>
29 StatisticsImageFilter<TInputImage>
30 ::StatisticsImageFilter(): m_ThreadSum(1), m_SumOfSquares(1), m_Count(1), m_ThreadMin(1), m_ThreadMax(1)
31 {
32  // first output is a copy of the image, DataObject created by
33  // superclass
34  //
35  // allocate the data objects for the outputs which are
36  // just decorators around pixel types
37  for (int i=1; i < 3; ++i)
38  {
39  typename PixelObjectType::Pointer output
40  = static_cast<PixelObjectType*>(this->MakeOutput(i).GetPointer());
41  this->ProcessObject::SetNthOutput(i, output.GetPointer());
42  }
43  // allocate the data objects for the outputs which are
44  // just decorators around real types
45  for (int i=3; i < 7; ++i)
46  {
47  typename RealObjectType::Pointer output
48  = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
49  this->ProcessObject::SetNthOutput(i, output.GetPointer());
50  }
51 
52  this->GetMinimumOutput()->Set( NumericTraits<PixelType>::max() );
53  this->GetMaximumOutput()->Set( NumericTraits<PixelType>::NonpositiveMin() );
54  this->GetMeanOutput()->Set( NumericTraits<RealType>::max() );
55  this->GetSigmaOutput()->Set( NumericTraits<RealType>::max() );
56  this->GetVarianceOutput()->Set( NumericTraits<RealType>::max() );
57  this->GetSumOutput()->Set( NumericTraits<RealType>::Zero );
58 }
59 
60 
61 template<class TInputImage>
62 DataObject::Pointer
63 StatisticsImageFilter<TInputImage>
64 ::MakeOutput(unsigned int output)
65 {
66  switch (output)
67  {
68  case 0:
69  return static_cast<DataObject*>(TInputImage::New().GetPointer());
70  break;
71  case 1:
72  return static_cast<DataObject*>(PixelObjectType::New().GetPointer());
73  break;
74  case 2:
75  return static_cast<DataObject*>(PixelObjectType::New().GetPointer());
76  break;
77  case 3:
78  case 4:
79  case 5:
80  case 6:
81  return static_cast<DataObject*>(RealObjectType::New().GetPointer());
82  break;
83  default:
84  // might as well make an image
85  return static_cast<DataObject*>(TInputImage::New().GetPointer());
86  break;
87  }
88 }
89 
90 
91 template<class TInputImage>
92 typename StatisticsImageFilter<TInputImage>::PixelObjectType*
93 StatisticsImageFilter<TInputImage>
94 ::GetMinimumOutput()
95 {
96  return static_cast<PixelObjectType*>(this->ProcessObject::GetOutput(1));
97 }
98 
99 template<class TInputImage>
100 const typename StatisticsImageFilter<TInputImage>::PixelObjectType*
101 StatisticsImageFilter<TInputImage>
102 ::GetMinimumOutput() const
103 {
104  return static_cast<const PixelObjectType*>(this->ProcessObject::GetOutput(1));
105 }
106 
107 
108 template<class TInputImage>
109 typename StatisticsImageFilter<TInputImage>::PixelObjectType*
110 StatisticsImageFilter<TInputImage>
111 ::GetMaximumOutput()
112 {
113  return static_cast<PixelObjectType*>(this->ProcessObject::GetOutput(2));
114 }
115 
116 template<class TInputImage>
117 const typename StatisticsImageFilter<TInputImage>::PixelObjectType*
118 StatisticsImageFilter<TInputImage>
119 ::GetMaximumOutput() const
120 {
121  return static_cast<const PixelObjectType*>(this->ProcessObject::GetOutput(2));
122 }
123 
124 
125 template<class TInputImage>
126 typename StatisticsImageFilter<TInputImage>::RealObjectType*
127 StatisticsImageFilter<TInputImage>
128 ::GetMeanOutput()
129 {
130  return static_cast<RealObjectType*>(this->ProcessObject::GetOutput(3));
131 }
132 
133 template<class TInputImage>
134 const typename StatisticsImageFilter<TInputImage>::RealObjectType*
135 StatisticsImageFilter<TInputImage>
136 ::GetMeanOutput() const
137 {
138  return static_cast<const RealObjectType*>(this->ProcessObject::GetOutput(3));
139 }
140 
141 
142 template<class TInputImage>
143 typename StatisticsImageFilter<TInputImage>::RealObjectType*
144 StatisticsImageFilter<TInputImage>
145 ::GetSigmaOutput()
146 {
147  return static_cast<RealObjectType*>(this->ProcessObject::GetOutput(4));
148 }
149 
150 template<class TInputImage>
151 const typename StatisticsImageFilter<TInputImage>::RealObjectType*
152 StatisticsImageFilter<TInputImage>
153 ::GetSigmaOutput() const
154 {
155  return static_cast<const RealObjectType*>(this->ProcessObject::GetOutput(4));
156 }
157 
158 
159 template<class TInputImage>
160 typename StatisticsImageFilter<TInputImage>::RealObjectType*
161 StatisticsImageFilter<TInputImage>
162 ::GetVarianceOutput()
163 {
164  return static_cast<RealObjectType*>(this->ProcessObject::GetOutput(5));
165 }
166 
167 template<class TInputImage>
168 const typename StatisticsImageFilter<TInputImage>::RealObjectType*
169 StatisticsImageFilter<TInputImage>
170 ::GetVarianceOutput() const
171 {
172  return static_cast<const RealObjectType*>(this->ProcessObject::GetOutput(5));
173 }
174 
175 
176 template<class TInputImage>
177 typename StatisticsImageFilter<TInputImage>::RealObjectType*
178 StatisticsImageFilter<TInputImage>
179 ::GetSumOutput()
180 {
181  return static_cast<RealObjectType*>(this->ProcessObject::GetOutput(6));
182 }
183 
184 template<class TInputImage>
185 const typename StatisticsImageFilter<TInputImage>::RealObjectType*
186 StatisticsImageFilter<TInputImage>
187 ::GetSumOutput() const
188 {
189  return static_cast<const RealObjectType*>(this->ProcessObject::GetOutput(6));
190 }
191 
192 template<class TInputImage>
193 void
194 StatisticsImageFilter<TInputImage>
195 ::GenerateInputRequestedRegion()
196 {
197  Superclass::GenerateInputRequestedRegion();
198  if ( this->GetInput() )
199  {
200  InputImagePointer image =
201  const_cast< typename Superclass::InputImageType * >( this->GetInput() );
202  image->SetRequestedRegionToLargestPossibleRegion();
203  }
204 }
205 
206 template<class TInputImage>
207 void
208 StatisticsImageFilter<TInputImage>
209 ::EnlargeOutputRequestedRegion(DataObject *data)
210 {
211  Superclass::EnlargeOutputRequestedRegion(data);
212  data->SetRequestedRegionToLargestPossibleRegion();
213 }
214 
215 
216 template<class TInputImage>
217 void
218 StatisticsImageFilter<TInputImage>
219 ::AllocateOutputs()
220 {
221  // Pass the input through as the output
222  InputImagePointer image =
223  const_cast< TInputImage * >( this->GetInput() );
224  this->GraftOutput( image );
225 
226  // Nothing that needs to be allocated for the remaining outputs
227 }
228 
229 template<class TInputImage>
230 void
231 StatisticsImageFilter<TInputImage>
232 ::BeforeThreadedGenerateData()
233 {
234  int numberOfThreads = this->GetNumberOfThreads();
235 
236  // Resize the thread temporaries
237  m_Count.SetSize(numberOfThreads);
238  m_SumOfSquares.SetSize(numberOfThreads);
239  m_ThreadSum.SetSize(numberOfThreads);
240  m_ThreadMin.SetSize(numberOfThreads);
241  m_ThreadMax.SetSize(numberOfThreads);
242 
243  // Initialize the temporaries
244  m_Count.Fill(NumericTraits<long>::Zero);
245  m_ThreadSum.Fill(NumericTraits<RealType>::Zero);
246  m_SumOfSquares.Fill(NumericTraits<RealType>::Zero);
247  m_ThreadMin.Fill(NumericTraits<PixelType>::max());
248  m_ThreadMax.Fill(NumericTraits<PixelType>::NonpositiveMin());
249 }
250 
251 template<class TInputImage>
252 void
253 StatisticsImageFilter<TInputImage>
254 ::AfterThreadedGenerateData()
255 {
256  int i;
257  long count;
258  RealType sumOfSquares;
259 
260  int numberOfThreads = this->GetNumberOfThreads();
261 
262  PixelType minimum;
263  PixelType maximum;
264  RealType mean;
265  RealType sigma;
266  RealType variance;
267  RealType sum;
268 
269  sum = sumOfSquares = NumericTraits<RealType>::Zero;
270  count = 0;
271 
272  // Find the min/max over all threads and accumulate count, sum and
273  // sum of squares
274  minimum = NumericTraits<PixelType>::max();
275  maximum = NumericTraits<PixelType>::NonpositiveMin();
276  for( i = 0; i < numberOfThreads; i++)
277  {
278  count += m_Count[i];
279  sum += m_ThreadSum[i];
280  sumOfSquares += m_SumOfSquares[i];
281 
282  if (m_ThreadMin[i] < minimum)
283  {
284  minimum = m_ThreadMin[i];
285  }
286  if (m_ThreadMax[i] > maximum)
287  {
288  maximum = m_ThreadMax[i];
289  }
290  }
291  // compute statistics
292  mean = sum / static_cast<RealType>(count);
293 
294  // unbiased estimate
295  variance = (sumOfSquares - (sum*sum / static_cast<RealType>(count)))
296  / (static_cast<RealType>(count) - 1);
297  sigma = vcl_sqrt(variance);
298 
299  // Set the outputs
300  this->GetMinimumOutput()->Set( minimum );
301  this->GetMaximumOutput()->Set( maximum );
302  this->GetMeanOutput()->Set( mean );
303  this->GetSigmaOutput()->Set( sigma );
304  this->GetVarianceOutput()->Set( variance );
305  this->GetSumOutput()->Set( sum );
306 }
307 
308 template<class TInputImage>
309 void
310 StatisticsImageFilter<TInputImage>
311 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
312  int threadId)
313 {
314  RealType realValue;
315  PixelType value;
316  ImageRegionConstIterator<TInputImage> it (this->GetInput(), outputRegionForThread);
317 
318  // support progress methods/callbacks
319  ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
320 
321  // do the work
322  while (!it.IsAtEnd())
323  {
324  value = it.Get();
325  realValue = static_cast<RealType>( value );
326  if (value < m_ThreadMin[threadId])
327  {
328  m_ThreadMin[threadId] = value;
329  }
330  if (value > m_ThreadMax[threadId])
331  {
332  m_ThreadMax[threadId] = value;
333  }
334 
335  m_ThreadSum[threadId] += realValue;
336  m_SumOfSquares[threadId] += (realValue * realValue);
337  m_Count[threadId]++;
338  ++it;
339  progress.CompletedPixel();
340  }
341 }
342 
343 template <class TImage>
344 void
345 StatisticsImageFilter<TImage>
346 ::PrintSelf(std::ostream& os, Indent indent) const
347 {
348  Superclass::PrintSelf(os,indent);
349 
350 // Trying to debug wrapping problem with VS6
351 #if !(defined(_MSC_VER) && (_MSC_VER <= 1300))
352  os << indent << "Minimum: "
353  << static_cast<typename NumericTraits<PixelType>::PrintType>(this->GetMinimum()) << std::endl;
354  os << indent << "Maximum: "
355  << static_cast<typename NumericTraits<PixelType>::PrintType>(this->GetMaximum()) << std::endl;
356 #endif
357  os << indent << "Sum: " << this->GetSum() << std::endl;
358  os << indent << "Mean: " << this->GetMean() << std::endl;
359  os << indent << "Sigma: " << this->GetSigma() << std::endl;
360  os << indent << "Variance: " << this->GetVariance() << std::endl;
361 }
362 
363 
364 }// end namespace itk
365 #endif
Generated at Sun Feb 3 2013 00:08:39 for Orfeo Toolbox with doxygen 1.8.1.1