OTB  5.9.0
Orfeo Toolbox
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());
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());
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>
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;
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  if (count > 0)
271  {
272  // compute statistics
273  mean = sum / static_cast<RealType>(count);
274 
275  if (count > 1)
276  {
277  // unbiased estimate
278  variance = (sumOfSquares - (sum * sum / static_cast<RealType>(count)))
279  / static_cast<RealType>(count - 1);
280  sigma = vcl_sqrt(variance);
281  }
282  }
283  else
284  {
285  itkWarningMacro(<<"No pixel found to compute statistics!");
286  }
287 
288  // Set the outputs
289  this->GetMinimumOutput()->Set(minimum);
290  this->GetMaximumOutput()->Set(maximum);
291  this->GetMeanOutput()->Set(mean);
292  this->GetSigmaOutput()->Set(sigma);
293  this->GetVarianceOutput()->Set(variance);
294  this->GetSumOutput()->Set(sum);
295 }
296 
297 template<class TInputImage>
298 void
301 {
302  int numberOfThreads = this->GetNumberOfThreads();
303 
304  // Resize the thread temporaries
305  m_Count.SetSize(numberOfThreads);
306  m_SumOfSquares.SetSize(numberOfThreads);
307  m_ThreadSum.SetSize(numberOfThreads);
308  m_ThreadMin.SetSize(numberOfThreads);
309  m_ThreadMax.SetSize(numberOfThreads);
310  // Initialize the temporaries
311  m_Count.Fill(itk::NumericTraits<long>::Zero);
312  m_ThreadSum.Fill(itk::NumericTraits<RealType>::Zero);
313  m_SumOfSquares.Fill(itk::NumericTraits<RealType>::Zero);
314  m_ThreadMin.Fill(itk::NumericTraits<PixelType>::max());
316  if (m_IgnoreInfiniteValues)
317  {
318  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(numberOfThreads, 0);
319  }
320 
321  if (m_IgnoreUserDefinedValue)
322  {
323  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
324  }
325 }
326 
327 template<class TInputImage>
328 void
330 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
331  itk::ThreadIdType threadId)
332 {
336  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput(0));
337  // support progress methods/callbacks
338  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
340 
341  RealType realValue;
342  PixelType value;
343 
344  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
345 
346  it.GoToBegin();
347  // do the work
348  while (!it.IsAtEnd())
349  {
350  value = it.Get();
351  realValue = static_cast<RealType>(value);
352  if (m_IgnoreInfiniteValues && !(vnl_math_isfinite(realValue)))
353  {
354  m_IgnoredInfinitePixelCount[threadId] ++;
355  }
356  else
357  {
358  if (m_IgnoreUserDefinedValue && (value == m_UserIgnoredValue))
359  {
360  m_IgnoredUserPixelCount[threadId] ++;
361  }
362  else
363  {
364  if (value < m_ThreadMin[threadId])
365  {
366  m_ThreadMin[threadId] = value;
367  }
368  if (value > m_ThreadMax[threadId])
369  {
370  m_ThreadMax[threadId] = value;
371  }
372 
373  m_ThreadSum[threadId] += realValue;
374  m_SumOfSquares[threadId] += (realValue * realValue);
375  m_Count[threadId]++;
376  }
377  }
378  ++it;
379  progress.CompletedPixel();
380  }
381 }
382 template <class TImage>
383 void
385 ::PrintSelf(std::ostream& os, itk::Indent indent) const
386 {
387  Superclass::PrintSelf(os, indent);
388 
389  os << indent << "Minimum: "
390  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMinimum()) << std::endl;
391  os << indent << "Maximum: "
392  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMaximum()) << std::endl;
393  os << indent << "Sum: " << this->GetSum() << std::endl;
394  os << indent << "Mean: " << this->GetMean() << std::endl;
395  os << indent << "Sigma: " << this->GetSigma() << std::endl;
396  os << indent << "Variance: " << this->GetVariance() << std::endl;
397 }
398 } // end namespace otb
399 #endif
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) ITK_OVERRIDE
ObjectType * GetPointer() const
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) ITK_OVERRIDE
static ITK_CONSTEXPR_FUNC T max(const T &)
unsigned int ThreadIdType
static ITK_CONSTEXPR_FUNC T NonpositiveMin()
itk::NumericTraits< PixelType >::RealType RealType
bool IsAtEnd(void) const
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
PixelType Get(void) const
void PrintSelf(std::ostream &os, itk::Indent indent) const ITK_OVERRIDE
DataObject * GetOutput(const DataObjectIdentifierType &key)