Orfeo Toolbox  4.2
otbStreamingStatisticsImageFilter.txx
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: ORFEO Toolbox
4  Language: C++
5  Date: $Date$
6  Version: $Revision$
7 
8 
9  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
10  See OTBCopyright.txt for details.
11 
12 Some parts of this code are derived from ITK. See ITKCopyright.txt
13 for details.
14 
15 
16  This software is distributed WITHOUT ANY WARRANTY; without even
17  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
18  PURPOSE. See the above copyright notices for more information.
19 
20 =========================================================================*/
21 #ifndef __otbStreamingStatisticsImageFilter_txx
22 #define __otbStreamingStatisticsImageFilter_txx
24 
25 #include "itkImageRegionIterator.h"
27 #include "itkNumericTraits.h"
28 #include "itkProgressReporter.h"
29 #include "otbMacro.h"
30 
31 namespace otb
32 {
33 
34 template<class TInputImage>
37  : m_ThreadSum(1),
38  m_SumOfSquares(1),
39  m_Count(1),
40  m_ThreadMin(1),
41  m_ThreadMax(1),
42  m_IgnoreInfiniteValues(true),
43  m_IgnoreUserDefinedValue(false)
44 {
45  // first output is a copy of the image, DataObject created by
46  // superclass
47  //
48  // allocate the data objects for the outputs which are
49  // just decorators around pixel types
50  for (int i = 1; i < 3; ++i)
51  {
52  typename PixelObjectType::Pointer output
53  = static_cast<PixelObjectType*>(this->MakeOutput(i).GetPointer());
55  }
56  // allocate the data objects for the outputs which are
57  // just decorators around real types
58  for (int i = 3; i < 7; ++i)
59  {
60  typename RealObjectType::Pointer output
61  = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
63  }
64 
65  this->GetMinimumOutput()->Set(itk::NumericTraits<PixelType>::max());
66  this->GetMaximumOutput()->Set(itk::NumericTraits<PixelType>::NonpositiveMin());
67  this->GetMeanOutput()->Set(itk::NumericTraits<RealType>::max());
68  this->GetSigmaOutput()->Set(itk::NumericTraits<RealType>::max());
69  this->GetVarianceOutput()->Set(itk::NumericTraits<RealType>::max());
70  this->GetSumOutput()->Set(itk::NumericTraits<RealType>::Zero);
71 
72  // Initiate the infinite ignored pixel counters
73  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
74  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
75 
76  this->Reset();
77 }
78 
79 template<class TInputImage>
82 ::MakeOutput(unsigned int output)
83 {
84  switch (output)
85  {
86  case 0:
87  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
88  break;
89  case 1:
90  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
91  break;
92  case 2:
93  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
94  break;
95  case 3:
96  case 4:
97  case 5:
98  case 6:
99  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
100  break;
101  default:
102  // might as well make an image
103  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
104  break;
105  }
106 }
107 
108 template<class TInputImage>
112 {
113  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
114 }
115 
116 template<class TInputImage>
120 {
121  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
122 }
123 
124 template<class TInputImage>
128 {
129  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
130 }
131 
132 template<class TInputImage>
136 {
137  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
138 }
139 
140 template<class TInputImage>
144 {
145  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
146 }
147 
148 template<class TInputImage>
152 {
153  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
154 }
155 
156 template<class TInputImage>
160 {
161  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
162 }
163 
164 template<class TInputImage>
168 {
169  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
170 }
171 
172 template<class TInputImage>
176 {
177  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(5));
178 }
179 
180 template<class TInputImage>
184 {
185  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(5));
186 }
187 
188 template<class TInputImage>
192 {
193  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(6));
194 }
195 
196 template<class TInputImage>
200 {
201  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(6));
202 }
203 template<class TInputImage>
204 void
207 {
208  Superclass::GenerateOutputInformation();
209  if (this->GetInput())
210  {
211  this->GetOutput()->CopyInformation(this->GetInput());
212  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
213 
214  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
215  {
216  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
217  }
218  }
219 }
220 template<class TInputImage>
221 void
224 {
225  // This is commented to prevent the streaming of the whole image for the first stream strip
226  // It shall not cause any problem because the output image of this filter is not intended to be used.
227  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
228  //this->GraftOutput( image );
229  // Nothing that needs to be allocated for the remaining outputs
230 }
231 
232 template<class TInputImage>
233 void
236 {
237  int i;
238  long count;
239  RealType sumOfSquares;
240 
241  int numberOfThreads = this->GetNumberOfThreads();
242 
243  PixelType minimum;
244  PixelType maximum;
245  RealType mean;
246  RealType sigma;
247  RealType variance;
248  RealType sum;
249 
250  sum = sumOfSquares = itk::NumericTraits<RealType>::Zero;
251  count = 0;
252 
253  // Find the min/max over all threads and accumulate count, sum and
254  // sum of squares
257  for (i = 0; i < numberOfThreads; ++i)
258  {
259  count += m_Count[i];
260  sum += m_ThreadSum[i];
261  sumOfSquares += m_SumOfSquares[i];
262 
263  if (m_ThreadMin[i] < minimum)
264  {
265  minimum = m_ThreadMin[i];
266  }
267  if (m_ThreadMax[i] > maximum)
268  {
269  maximum = m_ThreadMax[i];
270  }
271  }
272  // compute statistics
273  mean = sum / static_cast<RealType>(count);
274 
275  // unbiased estimate
276  variance = (sumOfSquares - (sum * sum / static_cast<RealType>(count)))
277  / (static_cast<RealType>(count) - 1);
278  sigma = vcl_sqrt(variance);
279 
280  // Set the outputs
281  this->GetMinimumOutput()->Set(minimum);
282  this->GetMaximumOutput()->Set(maximum);
283  this->GetMeanOutput()->Set(mean);
284  this->GetSigmaOutput()->Set(sigma);
285  this->GetVarianceOutput()->Set(variance);
286  this->GetSumOutput()->Set(sum);
287 }
288 
289 template<class TInputImage>
290 void
293 {
294  int numberOfThreads = this->GetNumberOfThreads();
295 
296  // Resize the thread temporaries
297  m_Count.SetSize(numberOfThreads);
298  m_SumOfSquares.SetSize(numberOfThreads);
299  m_ThreadSum.SetSize(numberOfThreads);
300  m_ThreadMin.SetSize(numberOfThreads);
301  m_ThreadMax.SetSize(numberOfThreads);
302  // Initialize the temporaries
303  m_Count.Fill(itk::NumericTraits<long>::Zero);
304  m_ThreadSum.Fill(itk::NumericTraits<RealType>::Zero);
305  m_SumOfSquares.Fill(itk::NumericTraits<RealType>::Zero);
306  m_ThreadMin.Fill(itk::NumericTraits<PixelType>::max());
308  if (m_IgnoreInfiniteValues)
309  {
310  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(numberOfThreads, 0);
311  }
312 
313  if (m_IgnoreUserDefinedValue)
314  {
315  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
316  }
317 }
318 
319 template<class TInputImage>
320 void
322 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
323  itk::ThreadIdType threadId)
324 {
328  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput(0));
329  // support progress methods/callbacks
330  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
331 
332  RealType realValue;
333  PixelType value;
334 
335  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
336 
337  it.GoToBegin();
338  // do the work
339  while (!it.IsAtEnd())
340  {
341  value = it.Get();
342  realValue = static_cast<RealType>(value);
343  if (m_IgnoreInfiniteValues && !(vnl_math_isfinite(realValue)))
344  {
345  m_IgnoredInfinitePixelCount[threadId] ++;
346  }
347  else
348  {
349  if (m_IgnoreUserDefinedValue && (value == m_UserIgnoredValue))
350  {
351  m_IgnoredUserPixelCount[threadId] ++;
352  }
353  else
354  {
355  if (value < m_ThreadMin[threadId])
356  {
357  m_ThreadMin[threadId] = value;
358  }
359  if (value > m_ThreadMax[threadId])
360  {
361  m_ThreadMax[threadId] = value;
362  }
363 
364  m_ThreadSum[threadId] += realValue;
365  m_SumOfSquares[threadId] += (realValue * realValue);
366  m_Count[threadId]++;
367  }
368  }
369  ++it;
370  progress.CompletedPixel();
371  }
372 }
373 template <class TImage>
374 void
376 ::PrintSelf(std::ostream& os, itk::Indent indent) const
377 {
378  Superclass::PrintSelf(os, indent);
379 
380  os << indent << "Minimum: "
381  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMinimum()) << std::endl;
382  os << indent << "Maximum: "
383  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMaximum()) << std::endl;
384  os << indent << "Sum: " << this->GetSum() << std::endl;
385  os << indent << "Mean: " << this->GetMean() << std::endl;
386  os << indent << "Sigma: " << this->GetSigma() << std::endl;
387  os << indent << "Variance: " << this->GetVariance() << std::endl;
388 }
389 } // end namespace otb
390 #endif

Generated at Sat Aug 30 2014 16:27:20 for Orfeo Toolbox with doxygen 1.8.3.1