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