Orfeo ToolBox  4.2
Orfeo ToolBox is not a black box
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"
26 #include "itkProgressReporter.h"
27 #include "otbMacro.h"
28 
29 namespace otb
30 {
31 
32 template<class TInputImage>
35  : m_ThreadSum(1),
36  m_SumOfSquares(1),
37  m_Count(1),
38  m_ThreadMin(1),
39  m_ThreadMax(1),
40  m_IgnoreInfiniteValues(true),
41  m_IgnoreUserDefinedValue(false)
42 {
43  // first output is a copy of the image, DataObject created by
44  // superclass
45  //
46  // allocate the data objects for the outputs which are
47  // just decorators around pixel types
48  for (int i = 1; i < 3; ++i)
49  {
50  typename PixelObjectType::Pointer output
51  = static_cast<PixelObjectType*>(this->MakeOutput(i).GetPointer());
52  this->itk::ProcessObject::SetNthOutput(i, output.GetPointer());
53  }
54  // allocate the data objects for the outputs which are
55  // just decorators around real types
56  for (int i = 3; i < 7; ++i)
57  {
58  typename RealObjectType::Pointer output
59  = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
60  this->itk::ProcessObject::SetNthOutput(i, output.GetPointer());
61  }
62 
63  this->GetMinimumOutput()->Set(itk::NumericTraits<PixelType>::max());
64  this->GetMaximumOutput()->Set(itk::NumericTraits<PixelType>::NonpositiveMin());
65  this->GetMeanOutput()->Set(itk::NumericTraits<RealType>::max());
66  this->GetSigmaOutput()->Set(itk::NumericTraits<RealType>::max());
67  this->GetVarianceOutput()->Set(itk::NumericTraits<RealType>::max());
68  this->GetSumOutput()->Set(itk::NumericTraits<RealType>::Zero);
69 
70  // Initiate the infinite ignored pixel counters
71  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
72  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
73 
74  this->Reset();
75 }
76 
77 template<class TInputImage>
78 typename itk::DataObject::Pointer
80 ::MakeOutput(unsigned int output)
81 {
82  switch (output)
83  {
84  case 0:
85  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
86  break;
87  case 1:
88  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
89  break;
90  case 2:
91  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
92  break;
93  case 3:
94  case 4:
95  case 5:
96  case 6:
97  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
98  break;
99  default:
100  // might as well make an image
101  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
102  break;
103  }
104 }
105 
106 template<class TInputImage>
110 {
111  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
112 }
113 
114 template<class TInputImage>
118 {
119  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
120 }
121 
122 template<class TInputImage>
126 {
127  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
128 }
129 
130 template<class TInputImage>
134 {
135  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
136 }
137 
138 template<class TInputImage>
142 {
143  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
144 }
145 
146 template<class TInputImage>
150 {
151  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
152 }
153 
154 template<class TInputImage>
158 {
159  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
160 }
161 
162 template<class TInputImage>
166 {
167  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
168 }
169 
170 template<class TInputImage>
174 {
175  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(5));
176 }
177 
178 template<class TInputImage>
182 {
183  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(5));
184 }
185 
186 template<class TInputImage>
190 {
191  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(6));
192 }
193 
194 template<class TInputImage>
198 {
199  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(6));
200 }
201 template<class TInputImage>
202 void
205 {
206  Superclass::GenerateOutputInformation();
207  if (this->GetInput())
208  {
209  this->GetOutput()->CopyInformation(this->GetInput());
210  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
211 
212  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
213  {
214  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
215  }
216  }
217 }
218 template<class TInputImage>
219 void
222 {
223  // This is commented to prevent the streaming of the whole image for the first stream strip
224  // It shall not cause any problem because the output image of this filter is not intended to be used.
225  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
226  //this->GraftOutput( image );
227  // Nothing that needs to be allocated for the remaining outputs
228 }
229 
230 template<class TInputImage>
231 void
234 {
235  int i;
236  long count;
237  RealType sumOfSquares;
238 
239  int numberOfThreads = this->GetNumberOfThreads();
240 
241  PixelType minimum;
242  PixelType maximum;
243  RealType mean;
244  RealType sigma;
245  RealType variance;
246  RealType sum;
247 
248  sum = sumOfSquares = itk::NumericTraits<RealType>::Zero;
249  count = 0;
250 
251  // Find the min/max over all threads and accumulate count, sum and
252  // sum of squares
255  for (i = 0; i < numberOfThreads; ++i)
256  {
257  count += m_Count[i];
258  sum += m_ThreadSum[i];
259  sumOfSquares += m_SumOfSquares[i];
260 
261  if (m_ThreadMin[i] < minimum)
262  {
263  minimum = m_ThreadMin[i];
264  }
265  if (m_ThreadMax[i] > maximum)
266  {
267  maximum = m_ThreadMax[i];
268  }
269  }
270  // compute statistics
271  mean = sum / static_cast<RealType>(count);
272 
273  // unbiased estimate
274  variance = (sumOfSquares - (sum * sum / static_cast<RealType>(count)))
275  / (static_cast<RealType>(count) - 1);
276  sigma = vcl_sqrt(variance);
277 
278  // Set the outputs
279  this->GetMinimumOutput()->Set(minimum);
280  this->GetMaximumOutput()->Set(maximum);
281  this->GetMeanOutput()->Set(mean);
282  this->GetSigmaOutput()->Set(sigma);
283  this->GetVarianceOutput()->Set(variance);
284  this->GetSumOutput()->Set(sum);
285 }
286 
287 template<class TInputImage>
288 void
291 {
292  int numberOfThreads = this->GetNumberOfThreads();
293 
294  // Resize the thread temporaries
295  m_Count.SetSize(numberOfThreads);
296  m_SumOfSquares.SetSize(numberOfThreads);
297  m_ThreadSum.SetSize(numberOfThreads);
298  m_ThreadMin.SetSize(numberOfThreads);
299  m_ThreadMax.SetSize(numberOfThreads);
300  // Initialize the temporaries
301  m_Count.Fill(itk::NumericTraits<long>::Zero);
302  m_ThreadSum.Fill(itk::NumericTraits<RealType>::Zero);
303  m_SumOfSquares.Fill(itk::NumericTraits<RealType>::Zero);
304  m_ThreadMin.Fill(itk::NumericTraits<PixelType>::max());
306  if (m_IgnoreInfiniteValues)
307  {
308  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(numberOfThreads, 0);
309  }
310 
311  if (m_IgnoreUserDefinedValue)
312  {
313  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
314  }
315 }
316 
317 template<class TInputImage>
318 void
320 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
321  itk::ThreadIdType threadId)
322 {
326  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput(0));
327  // support progress methods/callbacks
328  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
329 
330  RealType realValue;
331  PixelType value;
332 
333  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
334 
335  it.GoToBegin();
336  // do the work
337  while (!it.IsAtEnd())
338  {
339  value = it.Get();
340  realValue = static_cast<RealType>(value);
341  if (m_IgnoreInfiniteValues && !(vnl_math_isfinite(realValue)))
342  {
343  m_IgnoredInfinitePixelCount[threadId] ++;
344  }
345  else
346  {
347  if (m_IgnoreUserDefinedValue && (value == m_UserIgnoredValue))
348  {
349  m_IgnoredUserPixelCount[threadId] ++;
350  }
351  else
352  {
353  if (value < m_ThreadMin[threadId])
354  {
355  m_ThreadMin[threadId] = value;
356  }
357  if (value > m_ThreadMax[threadId])
358  {
359  m_ThreadMax[threadId] = value;
360  }
361 
362  m_ThreadSum[threadId] += realValue;
363  m_SumOfSquares[threadId] += (realValue * realValue);
364  m_Count[threadId]++;
365  }
366  }
367  ++it;
368  progress.CompletedPixel();
369  }
370 }
371 template <class TImage>
372 void
374 ::PrintSelf(std::ostream& os, itk::Indent indent) const
375 {
376  Superclass::PrintSelf(os, indent);
377 
378  os << indent << "Minimum: "
379  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMinimum()) << std::endl;
380  os << indent << "Maximum: "
381  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMaximum()) << std::endl;
382  os << indent << "Sum: " << this->GetSum() << std::endl;
383  os << indent << "Mean: " << this->GetMean() << std::endl;
384  os << indent << "Sigma: " << this->GetSigma() << std::endl;
385  os << indent << "Variance: " << this->GetVariance() << std::endl;
386 }
387 } // end namespace otb
388 #endif
virtual void Set(const T &val)
ObjectType * GetPointer() const
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId)
void PrintSelf(std::ostream &os, itk::Indent indent) const
static T max(const T &)
itk::NumericTraits< PixelType >::RealType RealType
static T NonpositiveMin()
bool IsAtEnd(void) const
PixelType Get(void) const
unsigned int ThreadIdType
virtual DataObjectPointer MakeOutput(unsigned int idx)