OTB  7.4.0
Orfeo Toolbox
otbStreamingStatisticsImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 1999-2011 Insight Software Consortium
3  * Copyright (C) 2005-2020 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_hxx
23 #define otbStreamingStatisticsImageFilter_hxx
25 
26 #include "itkImageRegionIterator.h"
27 #include "itkProgressReporter.h"
28 #include "otbMacro.h"
29 
30 namespace otb
31 {
32 
33 template <class TInputImage>
35  : m_ThreadSum(1), m_SumOfSquares(1), m_Count(1), m_ThreadMin(1), m_ThreadMax(1), m_IgnoreInfiniteValues(true), m_IgnoreUserDefinedValue(false)
36 {
37  // first output is a copy of the image, DataObject created by
38  // superclass
39  //
40  // allocate the data objects for the outputs which are
41  // just decorators around pixel types
42  for (int i = 1; i < 3; ++i)
43  {
44  typename PixelObjectType::Pointer output = static_cast<PixelObjectType*>(this->MakeOutput(i).GetPointer());
45  this->itk::ProcessObject::SetNthOutput(i, output.GetPointer());
46  }
47  // allocate the data objects for the outputs which are
48  // just decorators around real types
49  for (int i = 3; i < 7; ++i)
50  {
51  typename RealObjectType::Pointer output = static_cast<RealObjectType*>(this->MakeOutput(i).GetPointer());
52  this->itk::ProcessObject::SetNthOutput(i, output.GetPointer());
53  }
54 
55  this->GetMinimumOutput()->Set(itk::NumericTraits<PixelType>::max());
56  this->GetMaximumOutput()->Set(itk::NumericTraits<PixelType>::NonpositiveMin());
57  this->GetMeanOutput()->Set(itk::NumericTraits<RealType>::max());
58  this->GetSigmaOutput()->Set(itk::NumericTraits<RealType>::max());
59  this->GetVarianceOutput()->Set(itk::NumericTraits<RealType>::max());
60  this->GetSumOutput()->Set(itk::NumericTraits<RealType>::Zero);
61 
62  // Initiate the infinite ignored pixel counters
63  m_IgnoredInfinitePixelCount = std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
64  m_IgnoredUserPixelCount = std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
65 
66  this->Reset();
67 }
68 
69 template <class TInputImage>
71 {
72  switch (output)
73  {
74  case 0:
75  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
76  break;
77  case 1:
78  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
79  break;
80  case 2:
81  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
82  break;
83  case 3:
84  case 4:
85  case 5:
86  case 6:
87  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
88  break;
89  default:
90  // might as well make an image
91  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
92  break;
93  }
94 }
95 
96 template <class TInputImage>
98 {
99  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
100 }
101 
102 template <class TInputImage>
104 {
105  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
106 }
107 
108 template <class TInputImage>
110 {
111  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
112 }
113 
114 template <class TInputImage>
116 {
117  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
118 }
119 
120 template <class TInputImage>
122 {
123  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
124 }
125 
126 template <class TInputImage>
128 {
129  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(3));
130 }
131 
132 template <class TInputImage>
134 {
135  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
136 }
137 
138 template <class TInputImage>
140 {
141  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(4));
142 }
143 
144 template <class TInputImage>
146 {
147  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(5));
148 }
149 
150 template <class TInputImage>
152 {
153  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(5));
154 }
155 
156 template <class TInputImage>
158 {
159  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(6));
160 }
161 
162 template <class TInputImage>
164 {
165  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(6));
166 }
167 template <class TInputImage>
169 {
170  Superclass::GenerateOutputInformation();
171  if (this->GetInput())
172  {
173  this->GetOutput()->CopyInformation(this->GetInput());
174  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
175 
176  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
177  {
178  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
179  }
180  }
181 }
182 template <class TInputImage>
184 {
185  // This is commented to prevent the streaming of the whole image for the first stream strip
186  // It shall not cause any problem because the output image of this filter is not intended to be used.
187  // InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
188  // this->GraftOutput( image );
189  // Nothing that needs to be allocated for the remaining outputs
190 }
191 
192 template <class TInputImage>
194 {
195  int i;
196  long count;
197  RealType sumOfSquares;
198 
199  int numberOfThreads = this->GetNumberOfThreads();
200 
201  PixelType minimum;
202  PixelType maximum;
203  RealType mean = itk::NumericTraits<RealType>::Zero;
204  RealType sigma = itk::NumericTraits<RealType>::Zero;
205  RealType variance = itk::NumericTraits<RealType>::Zero;
206  RealType sum;
207 
208  sum = sumOfSquares = itk::NumericTraits<RealType>::Zero;
209  count = 0;
210 
211  // Find the min/max over all threads and accumulate count, sum and
212  // sum of squares
213  minimum = itk::NumericTraits<PixelType>::max();
214  maximum = itk::NumericTraits<PixelType>::NonpositiveMin();
215  for (i = 0; i < numberOfThreads; ++i)
216  {
217  count += m_Count[i];
218  sum += m_ThreadSum[i];
219  sumOfSquares += m_SumOfSquares[i];
220 
221  if (m_ThreadMin[i] < minimum)
222  {
223  minimum = m_ThreadMin[i];
224  }
225  if (m_ThreadMax[i] > maximum)
226  {
227  maximum = m_ThreadMax[i];
228  }
229  }
230  if (count > 0)
231  {
232  // compute statistics
233  mean = sum / static_cast<RealType>(count);
234 
235  if (count > 1)
236  {
237  // unbiased estimate
238  variance = (sumOfSquares - (sum * sum / static_cast<RealType>(count))) / static_cast<RealType>(count - 1);
239  sigma = std::sqrt(variance);
240  }
241  }
242  else
243  {
244  itkWarningMacro(<< "No pixel found to compute statistics!");
245  }
246 
247  // Set the outputs
248  this->GetMinimumOutput()->Set(minimum);
249  this->GetMaximumOutput()->Set(maximum);
250  this->GetMeanOutput()->Set(mean);
251  this->GetSigmaOutput()->Set(sigma);
252  this->GetVarianceOutput()->Set(variance);
253  this->GetSumOutput()->Set(sum);
254 }
255 
256 template <class TInputImage>
258 {
259  int numberOfThreads = this->GetNumberOfThreads();
260 
261  // Resize the thread temporaries
262  m_Count.SetSize(numberOfThreads);
263  m_SumOfSquares.SetSize(numberOfThreads);
264  m_ThreadSum.SetSize(numberOfThreads);
265  m_ThreadMin.SetSize(numberOfThreads);
266  m_ThreadMax.SetSize(numberOfThreads);
267  // Initialize the temporaries
268  m_Count.Fill(itk::NumericTraits<long>::Zero);
269  m_ThreadSum.Fill(itk::NumericTraits<RealType>::Zero);
270  m_SumOfSquares.Fill(itk::NumericTraits<RealType>::Zero);
271  m_ThreadMin.Fill(itk::NumericTraits<PixelType>::max());
272  m_ThreadMax.Fill(itk::NumericTraits<PixelType>::NonpositiveMin());
274  {
275  m_IgnoredInfinitePixelCount = std::vector<unsigned int>(numberOfThreads, 0);
276  }
277 
279  {
280  m_IgnoredUserPixelCount = std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
281  }
282 }
283 
284 template <class TInputImage>
285 void PersistentStatisticsImageFilter<TInputImage>::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
286 {
290  InputImagePointer inputPtr = const_cast<TInputImage*>(this->GetInput(0));
291  // support progress methods/callbacks
292  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
294 
295  RealType realValue;
297 
298  itk::ImageRegionConstIterator<TInputImage> it(inputPtr, outputRegionForThread);
299 
300  it.GoToBegin();
301  // do the work
302  while (!it.IsAtEnd())
303  {
304  value = it.Get();
305  realValue = static_cast<RealType>(value);
306  if (m_IgnoreInfiniteValues && !(vnl_math_isfinite(realValue)))
307  {
308  m_IgnoredInfinitePixelCount[threadId]++;
309  }
310  else
311  {
313  {
314  m_IgnoredUserPixelCount[threadId]++;
315  }
316  else
317  {
318  if (value < m_ThreadMin[threadId])
319  {
320  m_ThreadMin[threadId] = value;
321  }
322  if (value > m_ThreadMax[threadId])
323  {
324  m_ThreadMax[threadId] = value;
325  }
326 
327  m_ThreadSum[threadId] += realValue;
328  m_SumOfSquares[threadId] += (realValue * realValue);
329  m_Count[threadId]++;
330  }
331  }
332  ++it;
333  progress.CompletedPixel();
334  }
335 }
336 template <class TImage>
337 void PersistentStatisticsImageFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const
338 {
339  Superclass::PrintSelf(os, indent);
340 
341  os << indent << "Minimum: " << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMinimum()) << std::endl;
342  os << indent << "Maximum: " << static_cast<typename itk::NumericTraits<PixelType>::PrintType>(this->GetMaximum()) << std::endl;
343  os << indent << "Sum: " << this->GetSum() << std::endl;
344  os << indent << "Mean: " << this->GetMean() << std::endl;
345  os << indent << "Sigma: " << this->GetSigma() << std::endl;
346  os << indent << "Variance: " << this->GetVariance() << std::endl;
347 }
348 } // end namespace otb
349 #endif
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
void PrintSelf(std::ostream &os, itk::Indent indent) const override
itk::SimpleDataObjectDecorator< RealType > RealObjectType
itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
itk::SimpleDataObjectDecorator< PixelType > PixelObjectType
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override
static constexpr GLenum value() noexcept
itk::NumericTraits< PixelType >::RealType RealType
void PrintSelf(std::ostream &os, itk::Indent indent) const override