Orfeo Toolbox  4.0
otbStreamingStatisticsVectorImageFilter.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 
13  This software is distributed WITHOUT ANY WARRANTY; without even
14  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
15  PURPOSE. See the above copyright notices for more information.
16 
17 =========================================================================*/
18 #ifndef __otbStreamingStatisticsVectorImageFilter_txx
19 #define __otbStreamingStatisticsVectorImageFilter_txx
21 
22 #include "itkImageRegionIterator.h"
24 #include "itkNumericTraits.h"
25 #include "itkProgressReporter.h"
26 #include "otbMacro.h"
27 
28 namespace otb
29 {
30 
31 template<class TInputImage, class TPrecision>
34  : m_EnableMinMax(true),
35  m_EnableFirstOrderStats(true),
36  m_EnableSecondOrderStats(true),
37  m_UseUnbiasedEstimator(true),
38  m_IgnoreInfiniteValues(true),
39  m_IgnoreUserDefinedValue(false)
40 {
41  // first output is a copy of the image, DataObject created by
42  // superclass
43 
44  // allocate the data objects for the outputs which are
45  // just decorators around vector/matrix types
46  for (unsigned int i = 1; i < 10; ++i)
47  {
48  this->itk::ProcessObject::SetNthOutput(i, this->MakeOutput(i).GetPointer());
49  }
50  // Initiate ignored pixel counters
51  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
52  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
53 }
54 
55 template<class TInputImage, class TPrecision>
58 ::MakeOutput(unsigned int output)
59 {
60  switch (output)
61  {
62  case 0:
63  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
64  break;
65  case 1:
66  case 2:
67  // min/max
68  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
69  break;
70  case 3:
71  case 4:
72  // mean / sum
73  return static_cast<itk::DataObject*>(RealPixelObjectType::New().GetPointer());
74  break;
75  case 5:
76  case 6:
77  // covariance / correlation
78  return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer());
79  break;
80  case 7:
81  case 8:
82  case 9:
83  // component mean, component covariance, component correlation
84  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
85  break;
86  default:
87  // might as well make an image
88  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
89  break;
90  }
91 }
92 
93 template<class TInputImage, class TPrecision>
97 {
98  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
99 }
100 
101 template<class TInputImage, class TPrecision>
105 {
106  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
107 }
108 
109 template<class TInputImage, class TPrecision>
113 {
114  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
115 }
116 
117 template<class TInputImage, class TPrecision>
121 {
122  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
123 }
124 
125 template<class TInputImage, class TPrecision>
129 {
130  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(7));
131 }
132 
133 template<class TInputImage, class TPrecision>
137 {
138  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(7));
139 }
140 
141 template<class TInputImage, class TPrecision>
145 {
146  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(8));
147 }
148 
149 template<class TInputImage, class TPrecision>
153 {
154  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(8));
155 }
156 
157 
158 template<class TInputImage, class TPrecision>
162 {
163  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(9));
164 }
165 
166 template<class TInputImage, class TPrecision>
170 {
171  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(9));
172 }
173 
174 template<class TInputImage, class TPrecision>
178 {
179  return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
180 }
181 
182 template<class TInputImage, class TPrecision>
186 {
187  return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
188 }
189 
190 template<class TInputImage, class TPrecision>
194 {
195  return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
196 }
197 
198 template<class TInputImage, class TPrecision>
202 {
203  return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
204 }
205 
206 template<class TInputImage, class TPrecision>
210 {
211  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
212 }
213 
214 template<class TInputImage, class TPrecision>
218 {
219  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
220 }
221 
222 template<class TInputImage, class TPrecision>
226 {
227  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(6));
228 }
229 
230 template<class TInputImage, class TPrecision>
234 {
235  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(6));
236 }
237 
238 template<class TInputImage, class TPrecision>
239 void
242 {
243  Superclass::GenerateOutputInformation();
244  if (this->GetInput())
245  {
246  this->GetOutput()->CopyInformation(this->GetInput());
247  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
248 
249  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
250  {
251  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
252  }
253  }
254 }
255 
256 template<class TInputImage, class TPrecision>
257 void
260 {
261  // This is commented to prevent the streaming of the whole image for the first stream strip
262  // It shall not cause any problem because the output image of this filter is not intended to be used.
263  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
264  //this->GraftOutput( image );
265  // Nothing that needs to be allocated for the remaining outputs
266 }
267 
268 template<class TInputImage, class TPrecision>
269 void
272 {
273  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
274  inputPtr->UpdateOutputInformation();
275 
276  unsigned int numberOfThreads = this->GetNumberOfThreads();
277  unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
278 
279  if (m_EnableMinMax)
280  {
281  PixelType tempPixel;
282  tempPixel.SetSize(numberOfComponent);
283 
285  this->GetMinimumOutput()->Set(tempPixel);
286 
288  this->GetMaximumOutput()->Set(tempPixel);
289 
290  PixelType tempTemporiesPixel;
291  tempTemporiesPixel.SetSize(numberOfComponent);
292  tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
293  m_ThreadMin = std::vector<PixelType>(numberOfThreads, tempTemporiesPixel);
294 
296  m_ThreadMax = std::vector<PixelType>(numberOfThreads, tempTemporiesPixel);
297  }
298 
299  if (m_EnableSecondOrderStats)
300  {
301  m_EnableFirstOrderStats = true;
302  }
303 
304  if (m_EnableFirstOrderStats)
305  {
306  RealPixelType zeroRealPixel;
307  zeroRealPixel.SetSize(numberOfComponent);
309  this->GetMeanOutput()->Set(zeroRealPixel);
310  this->GetSumOutput()->Set(zeroRealPixel);
311  m_ThreadFirstOrderAccumulators.resize(numberOfThreads);
312  std::fill(m_ThreadFirstOrderAccumulators.begin(), m_ThreadFirstOrderAccumulators.end(), zeroRealPixel);
313 
315  m_ThreadFirstOrderComponentAccumulators.resize(numberOfThreads);
316  std::fill(m_ThreadFirstOrderComponentAccumulators.begin(), m_ThreadFirstOrderComponentAccumulators.end(), zeroReal);
317 
318  }
319 
320  if (m_EnableSecondOrderStats)
321  {
322  MatrixType zeroMatrix;
323  zeroMatrix.SetSize(numberOfComponent, numberOfComponent);
325  this->GetCovarianceOutput()->Set(zeroMatrix);
326  this->GetCorrelationOutput()->Set(zeroMatrix);
327 
328  m_ThreadSecondOrderAccumulators.resize(numberOfThreads);
329  std::fill(m_ThreadSecondOrderAccumulators.begin(), m_ThreadSecondOrderAccumulators.end(), zeroMatrix);
330 
332  m_ThreadSecondOrderComponentAccumulators.resize(numberOfThreads);
333  std::fill(m_ThreadSecondOrderComponentAccumulators.begin(), m_ThreadSecondOrderComponentAccumulators.end(), zeroReal);
334  }
335 
336  if (m_IgnoreInfiniteValues)
337  {
338  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(numberOfThreads, 0);
339  }
340 
341  if (m_IgnoreUserDefinedValue)
342  {
343  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
344  }
345 }
346 
347 template<class TInputImage, class TPrecision>
348 void
351 {
352  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
353  const unsigned int nbPixels = inputPtr->GetLargestPossibleRegion().GetNumberOfPixels();
354  const unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
355 
356  PixelType minimum;
357  minimum.SetSize(numberOfComponent);
359  PixelType maximum;
360  maximum.SetSize(numberOfComponent);
362 
363  RealPixelType streamFirstOrderAccumulator(numberOfComponent);
364  streamFirstOrderAccumulator.Fill(itk::NumericTraits<PrecisionType>::Zero);
365  MatrixType streamSecondOrderAccumulator(numberOfComponent, numberOfComponent);
366  streamSecondOrderAccumulator.Fill(itk::NumericTraits<PrecisionType>::Zero);
367 
368  RealType streamFirstOrderComponentAccumulator = itk::NumericTraits<RealType>::Zero;
369  RealType streamSecondOrderComponentAccumulator = itk::NumericTraits<RealType>::Zero;
370 
371  unsigned int ignoredInfinitePixelCount = 0;
372  unsigned int ignoredUserPixelCount = 0;
373 
374  // Accumulate results from all threads
375  const itk::ThreadIdType numberOfThreads = this->GetNumberOfThreads();
376  for (itk::ThreadIdType threadId = 0; threadId < numberOfThreads; ++threadId)
377  {
378  if (m_EnableMinMax)
379  {
380  const PixelType& threadMin = m_ThreadMin [threadId];
381  const PixelType& threadMax = m_ThreadMax [threadId];
382 
383  for (unsigned int j = 0; j < numberOfComponent; ++j)
384  {
385  if (threadMin[j] < minimum[j])
386  {
387  minimum[j] = threadMin[j];
388  }
389  if (threadMax[j] > maximum[j])
390  {
391  maximum[j] = threadMax[j];
392  }
393  }
394  }
395 
396  if (m_EnableFirstOrderStats)
397  {
398  streamFirstOrderAccumulator += m_ThreadFirstOrderAccumulators[threadId];
399  streamFirstOrderComponentAccumulator += m_ThreadFirstOrderComponentAccumulators[threadId];
400  }
401 
402  if (m_EnableSecondOrderStats)
403  {
404  streamSecondOrderAccumulator += m_ThreadSecondOrderAccumulators[threadId];
405  streamSecondOrderComponentAccumulator += m_ThreadSecondOrderComponentAccumulators[threadId];
406  }
407  // Ignored Infinite Pixels
408  ignoredInfinitePixelCount += m_IgnoredInfinitePixelCount[threadId];
409  // Ignored Pixels
410  ignoredUserPixelCount += m_IgnoredUserPixelCount[threadId];
411  }
412 
413  // There cannot be more ignored pixels than read pixels.
414  assert( nbPixels >= ignoredInfinitePixelCount + ignoredUserPixelCount );
415  if( nbPixels < ignoredInfinitePixelCount + ignoredUserPixelCount )
416  {
417  itkExceptionMacro(
418  "nbPixels < ignoredInfinitePixelCount + ignoredUserPixelCount"
419  );
420  }
421 
422  unsigned int nbRelevantPixels =
423  nbPixels - (ignoredInfinitePixelCount + ignoredUserPixelCount);
424 
425  if( nbRelevantPixels==0 )
426  {
427  itkExceptionMacro(
428  "Statistics cannot be calculated with zero relevant pixels."
429  );
430  }
431 
432  // Final calculations
433  if (m_EnableMinMax)
434  {
435  this->GetMinimumOutput()->Set(minimum);
436  this->GetMaximumOutput()->Set(maximum);
437  }
438 
439  if (m_EnableFirstOrderStats)
440  {
441  this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixels * numberOfComponent));
442 
443  this->GetMeanOutput()->Set(streamFirstOrderAccumulator / nbRelevantPixels);
444  this->GetSumOutput()->Set(streamFirstOrderAccumulator);
445  }
446 
447  if (m_EnableSecondOrderStats)
448  {
449  MatrixType cor = streamSecondOrderAccumulator / nbRelevantPixels;
450  this->GetCorrelationOutput()->Set(cor);
451 
452  const RealPixelType& mean = this->GetMeanOutput()->Get();
453 
454  double regul = 1.0;
455 
456  if( m_UseUnbiasedEstimator && nbPixels>1 )
457  {
458  regul =
459  static_cast< double >( nbRelevantPixels ) /
460  ( static_cast< double >( nbRelevantPixels ) - 1.0 );
461  }
462 
463  MatrixType cov = cor;
464  for (unsigned int r = 0; r < numberOfComponent; ++r)
465  {
466  for (unsigned int c = 0; c < numberOfComponent; ++c)
467  {
468  cov(r, c) = regul * (cov(r, c) - mean[r] * mean[c]);
469  }
470  }
471  this->GetCovarianceOutput()->Set(cov);
472 
473  this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixels * numberOfComponent));
474  this->GetComponentCorrelationOutput()->Set(streamSecondOrderComponentAccumulator / (nbRelevantPixels * numberOfComponent));
475  this->GetComponentCovarianceOutput()->Set(
476  (nbRelevantPixels * numberOfComponent) / (nbRelevantPixels * numberOfComponent - 1)
477  * (this->GetComponentCorrelation()
478  - (this->GetComponentMean() * this->GetComponentMean())));
479  }
480 }
481 
482 template<class TInputImage, class TPrecision>
483 void
485 ::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
486  {
487  // Support progress methods/callbacks
488  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
489 
490  // Grab the input
491  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
492  PixelType& threadMin = m_ThreadMin [threadId];
493  PixelType& threadMax = m_ThreadMax [threadId];
494 
495 
496  itk::ImageRegionConstIteratorWithIndex<TInputImage> it(inputPtr, outputRegionForThread);
497 
498  for (it.GoToBegin(); !it.IsAtEnd(); ++it, progress.CompletedPixel())
499  {
500  const PixelType& vectorValue = it.Get();
501 
502  float finiteProbe = 0.;
503  bool userProbe = true;
504  for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
505  {
506  finiteProbe += (float)(vectorValue[j]);
507  userProbe = userProbe && (vectorValue[j] == m_UserIgnoredValue);
508  }
509 
510  if (m_IgnoreInfiniteValues && !(vnl_math_isfinite(finiteProbe)))
511  {
512  m_IgnoredInfinitePixelCount[threadId] ++;
513  }
514  else
515  {
516  if (m_IgnoreUserDefinedValue && (userProbe))
517  {
518  m_IgnoredUserPixelCount[threadId] ++;
519  }
520  else
521  {
522  for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
523  {
524 
525  if (m_EnableMinMax)
526  {
527  for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
528  {
529  if (vectorValue[j] < threadMin[j])
530  {
531  threadMin[j] = vectorValue[j];
532  }
533  if (vectorValue[j] > threadMax[j])
534  {
535  threadMax[j] = vectorValue[j];
536  }
537  }
538  }
539  }
540 
541  if (m_EnableFirstOrderStats)
542  {
543  RealPixelType& threadFirstOrder = m_ThreadFirstOrderAccumulators [threadId];
544  RealType& threadFirstOrderComponent = m_ThreadFirstOrderComponentAccumulators [threadId];
545 
546  threadFirstOrder += vectorValue;
547 
548  for (unsigned int i = 0; i < vectorValue.GetSize(); ++i)
549  {
550  threadFirstOrderComponent += vectorValue[i];
551  }
552  }
553 
554  if (m_EnableSecondOrderStats)
555  {
556  MatrixType& threadSecondOrder = m_ThreadSecondOrderAccumulators[threadId];
557  RealType& threadSecondOrderComponent = m_ThreadSecondOrderComponentAccumulators[threadId];
558 
559  for (unsigned int r = 0; r < threadSecondOrder.Rows(); ++r)
560  {
561  for (unsigned int c = 0; c < threadSecondOrder.Cols(); ++c)
562  {
563  threadSecondOrder(r, c) += vectorValue[r] * vectorValue[c];
564  }
565  }
566  threadSecondOrderComponent += vectorValue.GetSquaredNorm();
567  }
568  }
569  }
570  }
571 
572  }
573 
574 template <class TImage, class TPrecision>
575 void
577 ::PrintSelf(std::ostream& os, itk::Indent indent) const
578 {
579  Superclass::PrintSelf(os, indent);
580 
581  os << indent << "Min: " << this->GetMinimumOutput()->Get() << std::endl;
582  os << indent << "Max: " << this->GetMaximumOutput()->Get() << std::endl;
583  os << indent << "Mean: " << this->GetMeanOutput()->Get() << std::endl;
584  os << indent << "Covariance: " << this->GetCovarianceOutput()->Get() << std::endl;
585  os << indent << "Correlation: " << this->GetCorrelationOutput()->Get() << std::endl;
586  os << indent << "Component Mean: " << this->GetComponentMeanOutput()->Get() << std::endl;
587  os << indent << "Component Covariance: " << this->GetComponentCovarianceOutput()->Get() << std::endl;
588  os << indent << "Component Correlation: " << this->GetComponentCorrelationOutput()->Get() << std::endl;
589  os << indent << "UseUnbiasedEstimator: " << (this->m_UseUnbiasedEstimator ? "true" : "false") << std::endl;
590 }
591 
592 } // end namespace otb
593 #endif

Generated at Sat Mar 8 2014 16:20:29 for Orfeo Toolbox with doxygen 1.8.3.1