OTB  5.11.0
Orfeo Toolbox
otbStreamingStatisticsImageFilter.txx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
4  *
5  * This file is part of Orfeo Toolbox
6  *
7  * https://www.orfeo-toolbox.org/
8  *
9  * Licensed under the Apache License, Version 2.0 (the "License");
10  * you may not use this file except in compliance with the License.
11  * You may obtain a copy of the License at
12  *
13  * http://www.apache.org/licenses/LICENSE-2.0
14  *
15  * Unless required by applicable law or agreed to in writing, software
16  * distributed under the License is distributed on an "AS IS" BASIS,
17  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18  * See the License for the specific language governing permissions and
19  * limitations under the License.
20  */
21 
22 #ifndef otbStreamingStatisticsImageFilter_txx
23 #define otbStreamingStatisticsImageFilter_txx
25 
26 #include "itkImageRegionIterator.h"
27 #include "itkProgressReporter.h"
28 #include "otbMacro.h"
29 
30 namespace otb
31 {
32 
33 template<class TInputImage>
36  : m_ThreadSum(1),
37  m_SumOfSquares(1),
38  m_Count(1),
39  m_ThreadMin(1),
40  m_ThreadMax(1),
41  m_IgnoreInfiniteValues(true),
42  m_IgnoreUserDefinedValue(false)
43 {
44  // first output is a copy of the image, DataObject created by
45  // superclass
46  //
47  // allocate the data objects for the outputs which are
48  // just decorators around pixel types
49  for (int i = 1; i < 3; ++i)
50  {
51  typename PixelObjectType::Pointer output
52  = static_cast<PixelObjectType*>(this->MakeOutput(i).GetPointer());
54  }
55  // allocate the data objects for the outputs which are
56  // just decorators around real types
57  for (int i = 3; i < 7; ++i)
58  {
59  typename RealObjectType::Pointer output
60  = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
62  }
63 
64  this->GetMinimumOutput()->Set(itk::NumericTraits<PixelType>::max());
65  this->GetMaximumOutput()->Set(itk::NumericTraits<PixelType>::NonpositiveMin());
66  this->GetMeanOutput()->Set(itk::NumericTraits<RealType>::max());
67  this->GetSigmaOutput()->Set(itk::NumericTraits<RealType>::max());
68  this->GetVarianceOutput()->Set(itk::NumericTraits<RealType>::max());
69  this->GetSumOutput()->Set(itk::NumericTraits<RealType>::Zero);
70 
71  // Initiate the infinite ignored pixel counters
72  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
73  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
74 
75  this->Reset();
76 }
77 
78 template<class TInputImage>
82 {
83  switch (output)
84  {
85  case 0:
86  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
87  break;
88  case 1:
89  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
90  break;
91  case 2:
92  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
93  break;
94  case 3:
95  case 4:
96  case 5:
97  case 6:
98  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
99  break;
100  default:
101  // might as well make an image
102  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
103  break;
104  }
105 }
106 
107 template<class TInputImage>
111 {
112  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
113 }
114 
115 template<class TInputImage>
119 {
120  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
121 }
122 
123 template<class TInputImage>
127 {
128  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
129 }
130 
131 template<class TInputImage>
135 {
136  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
137 }
138 
139 template<class TInputImage>
143 {
144  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
145 }
146 
147 template<class TInputImage>
151 {
152  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
153 }
154 
155 template<class TInputImage>
159 {
160  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
161 }
162 
163 template<class TInputImage>
167 {
168  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
169 }
170 
171 template<class TInputImage>
175 {
176  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(5));
177 }
178 
179 template<class TInputImage>
183 {
184  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(5));
185 }
186 
187 template<class TInputImage>
191 {
192  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(6));
193 }
194 
195 template<class TInputImage>
199 {
200  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(6));
201 }
202 template<class TInputImage>
203 void
206 {
207  Superclass::GenerateOutputInformation();
208  if (this->GetInput())
209  {
210  this->GetOutput()->CopyInformation(this->GetInput());
211  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
212 
213  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
214  {
215  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
216  }
217  }
218 }
219 template<class TInputImage>
220 void
223 {
224  // This is commented to prevent the streaming of the whole image for the first stream strip
225  // It shall not cause any problem because the output image of this filter is not intended to be used.
226  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
227  //this->GraftOutput( image );
228  // Nothing that needs to be allocated for the remaining outputs
229 }
230 
231 template<class TInputImage>
232 void
235 {
236  int i;
237  long count;
238  RealType sumOfSquares;
239 
240  int numberOfThreads = this->GetNumberOfThreads();
241 
242  PixelType minimum;
243  PixelType maximum;
247  RealType sum;
248 
249  sum = sumOfSquares = itk::NumericTraits<RealType>::Zero;
250  count = 0;
251 
252  // Find the min/max over all threads and accumulate count, sum and
253  // sum of squares
256  for (i = 0; i < numberOfThreads; ++i)
257  {
258  count += m_Count[i];
259  sum += m_ThreadSum[i];
260  sumOfSquares += m_SumOfSquares[i];
261 
262  if (m_ThreadMin[i] < minimum)
263  {
264  minimum = m_ThreadMin[i];
265  }
266  if (m_ThreadMax[i] > maximum)
267  {
268  maximum = m_ThreadMax[i];
269  }
270  }
271  if (count > 0)
272  {
273  // compute statistics
274  mean = sum / static_cast<RealType>(count);
275 
276  if (count > 1)
277  {
278  // unbiased estimate
279  variance = (sumOfSquares - (sum * sum / static_cast<RealType>(count)))
280  / static_cast<RealType>(count - 1);
281  sigma = vcl_sqrt(variance);
282  }
283  }
284  else
285  {
286  itkWarningMacro(<<"No pixel found to compute statistics!");
287  }
288 
289  // Set the outputs
290  this->GetMinimumOutput()->Set(minimum);
291  this->GetMaximumOutput()->Set(maximum);
292  this->GetMeanOutput()->Set(mean);
293  this->GetSigmaOutput()->Set(sigma);
294  this->GetVarianceOutput()->Set(variance);
295  this->GetSumOutput()->Set(sum);
296 }
297 
298 template<class TInputImage>
299 void
302 {
303  int numberOfThreads = this->GetNumberOfThreads();
304 
305  // Resize the thread temporaries
306  m_Count.SetSize(numberOfThreads);
307  m_SumOfSquares.SetSize(numberOfThreads);
308  m_ThreadSum.SetSize(numberOfThreads);
309  m_ThreadMin.SetSize(numberOfThreads);
310  m_ThreadMax.SetSize(numberOfThreads);
311  // Initialize the temporaries
312  m_Count.Fill(itk::NumericTraits<long>::Zero);
313  m_ThreadSum.Fill(itk::NumericTraits<RealType>::Zero);
314  m_SumOfSquares.Fill(itk::NumericTraits<RealType>::Zero);
315  m_ThreadMin.Fill(itk::NumericTraits<PixelType>::max());
317  if (m_IgnoreInfiniteValues)
318  {
319  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(numberOfThreads, 0);
320  }
321 
322  if (m_IgnoreUserDefinedValue)
323  {
324  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
325  }
326 }
327 
328 template<class TInputImage>
329 void
331 ::ThreadedGenerateData(const RegionType& outputRegionForThread,
332  itk::ThreadIdType threadId)
333 {
337  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput(0));
338  // support progress methods/callbacks
339  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
341 
342  RealType realValue;
343  PixelType value;
344 
345  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
346 
347  it.GoToBegin();
348  // do the work
349  while (!it.IsAtEnd())
350  {
351  value = it.Get();
352  realValue = static_cast<RealType>(value);
353  if (m_IgnoreInfiniteValues && !(vnl_math_isfinite(realValue)))
354  {
355  m_IgnoredInfinitePixelCount[threadId] ++;
356  }
357  else
358  {
359  if (m_IgnoreUserDefinedValue && (value == m_UserIgnoredValue))
360  {
361  m_IgnoredUserPixelCount[threadId] ++;
362  }
363  else
364  {
365  if (value < m_ThreadMin[threadId])
366  {
367  m_ThreadMin[threadId] = value;
368  }
369  if (value > m_ThreadMax[threadId])
370  {
371  m_ThreadMax[threadId] = value;
372  }
373 
374  m_ThreadSum[threadId] += realValue;
375  m_SumOfSquares[threadId] += (realValue * realValue);
376  m_Count[threadId]++;
377  }
378  }
379  ++it;
380  progress.CompletedPixel();
381  }
382 }
383 template <class TImage>
384 void
386 ::PrintSelf(std::ostream& os, itk::Indent indent) const
387 {
388  Superclass::PrintSelf(os, indent);
389 
390  os << indent << "Minimum: "
391  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMinimum()) << std::endl;
392  os << indent << "Maximum: "
393  << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMaximum()) << std::endl;
394  os << indent << "Sum: " << this->GetSum() << std::endl;
395  os << indent << "Mean: " << this->GetMean() << std::endl;
396  os << indent << "Sigma: " << this->GetSigma() << std::endl;
397  os << indent << "Variance: " << this->GetVariance() << std::endl;
398 }
399 } // end namespace otb
400 #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)