OTB  7.4.0
Orfeo Toolbox
otbStreamingStatisticsVectorImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2020 Centre National d'Etudes Spatiales (CNES)
3  *
4  * This file is part of Orfeo Toolbox
5  *
6  * https://www.orfeo-toolbox.org/
7  *
8  * Licensed under the Apache License, Version 2.0 (the "License");
9  * you may not use this file except in compliance with the License.
10  * You may obtain a copy of the License at
11  *
12  * http://www.apache.org/licenses/LICENSE-2.0
13  *
14  * Unless required by applicable law or agreed to in writing, software
15  * distributed under the License is distributed on an "AS IS" BASIS,
16  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17  * See the License for the specific language governing permissions and
18  * limitations under the License.
19  */
20 
21 #ifndef otbStreamingStatisticsVectorImageFilter_hxx
22 #define otbStreamingStatisticsVectorImageFilter_hxx
24 
25 #include "itkImageRegionIterator.h"
26 #include "itkImageRegionConstIteratorWithIndex.h"
27 #include "itkProgressReporter.h"
28 #include "otbMacro.h"
29 
30 namespace otb
31 {
32 
33 template <class TInputImage, class TPrecision>
35  : m_EnableMinMax(true),
36  m_EnableFirstOrderStats(true),
37  m_EnableSecondOrderStats(true),
38  m_UseUnbiasedEstimator(true),
39  m_IgnoreInfiniteValues(true),
40  m_IgnoreUserDefinedValue(false),
41  m_UserIgnoredValue(itk::NumericTraits<InternalPixelType>::Zero)
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 vector/matrix types
48  for (unsigned int i = 1; i < 11; ++i)
49  {
50  this->itk::ProcessObject::SetNthOutput(i, this->MakeOutput(i).GetPointer());
51  }
52  // Initiate ignored pixel counters
53  m_IgnoredInfinitePixelCount = std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
54  m_IgnoredUserPixelCount = std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
55 }
56 
57 template <class TInputImage, class TPrecision>
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  case 10:
87  // relevant pixel
88  return static_cast<itk::DataObject*>(CountObjectType::New().GetPointer());
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, class TPrecision>
99 {
100  return static_cast<CountObjectType*>(this->itk::ProcessObject::GetOutput(10));
101 }
102 
103 template <class TInputImage, class TPrecision>
106 {
107  return static_cast<const CountObjectType*>(this->itk::ProcessObject::GetOutput(10));
108 }
109 
110 
111 template <class TInputImage, class TPrecision>
114 {
115  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
116 }
117 
118 template <class TInputImage, class TPrecision>
121 {
122  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
123 }
124 
125 template <class TInputImage, class TPrecision>
128 {
129  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
130 }
131 
132 template <class TInputImage, class TPrecision>
135 {
136  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
137 }
138 
139 template <class TInputImage, class TPrecision>
142 {
143  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(7));
144 }
145 
146 template <class TInputImage, class TPrecision>
149 {
150  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(7));
151 }
152 
153 template <class TInputImage, class TPrecision>
156 {
157  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(8));
158 }
159 
160 template <class TInputImage, class TPrecision>
163 {
164  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(8));
165 }
166 
167 
168 template <class TInputImage, class TPrecision>
171 {
172  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(9));
173 }
174 
175 template <class TInputImage, class TPrecision>
178 {
179  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(9));
180 }
181 
182 template <class TInputImage, class TPrecision>
185 {
186  return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
187 }
188 
189 template <class TInputImage, class TPrecision>
192 {
193  return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
194 }
195 
196 template <class TInputImage, class TPrecision>
199 {
200  return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
201 }
202 
203 template <class TInputImage, class TPrecision>
206 {
207  return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
208 }
209 
210 template <class TInputImage, class TPrecision>
213 {
214  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
215 }
216 
217 template <class TInputImage, class TPrecision>
220 {
221  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
222 }
223 
224 template <class TInputImage, class TPrecision>
227 {
228  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(6));
229 }
230 
231 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>
240 {
241  Superclass::GenerateOutputInformation();
242  if (this->GetInput())
243  {
244  this->GetOutput()->CopyInformation(this->GetInput());
245  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
246 
247  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
248  {
249  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
250  }
251  }
252 }
253 
254 template <class TInputImage, class TPrecision>
256 {
257  // This is commented to prevent the streaming of the whole image for the first stream strip
258  // It shall not cause any problem because the output image of this filter is not intended to be used.
259  // InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
260  // this->GraftOutput( image );
261  // Nothing that needs to be allocated for the remaining outputs
262 }
263 
264 template <class TInputImage, class TPrecision>
266 {
267  TInputImage* inputPtr = const_cast<TInputImage*>(this->GetInput());
268  inputPtr->UpdateOutputInformation();
269 
270  unsigned int numberOfThreads = this->GetNumberOfThreads();
271  unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
272 
273  if (m_EnableMinMax)
274  {
275  PixelType tempPixel;
276  tempPixel.SetSize(numberOfComponent);
277 
278  tempPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
279  this->GetMinimumOutput()->Set(tempPixel);
280 
281  tempPixel.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
282  this->GetMaximumOutput()->Set(tempPixel);
283 
284  PixelType tempTemporiesPixel;
285  tempTemporiesPixel.SetSize(numberOfComponent);
286  tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
287  m_ThreadMin = std::vector<PixelType>(numberOfThreads, tempTemporiesPixel);
288 
289  tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
290  m_ThreadMax = std::vector<PixelType>(numberOfThreads, tempTemporiesPixel);
291  }
292 
294  {
296  }
297 
299  {
300  RealPixelType zeroRealPixel;
301  zeroRealPixel.SetSize(numberOfComponent);
302  zeroRealPixel.Fill(itk::NumericTraits<PrecisionType>::ZeroValue());
303  this->GetMeanOutput()->Set(zeroRealPixel);
304  this->GetSumOutput()->Set(zeroRealPixel);
305  m_ThreadFirstOrderAccumulators.resize(numberOfThreads);
306  std::fill(m_ThreadFirstOrderAccumulators.begin(), m_ThreadFirstOrderAccumulators.end(), zeroRealPixel);
307 
308  RealType zeroReal = itk::NumericTraits<RealType>::ZeroValue();
309  m_ThreadFirstOrderComponentAccumulators.resize(numberOfThreads);
311  }
312 
314  {
315  MatrixType zeroMatrix;
316  zeroMatrix.SetSize(numberOfComponent, numberOfComponent);
317  zeroMatrix.Fill(itk::NumericTraits<PrecisionType>::Zero);
318  this->GetCovarianceOutput()->Set(zeroMatrix);
319  this->GetCorrelationOutput()->Set(zeroMatrix);
320 
321  m_ThreadSecondOrderAccumulators.resize(numberOfThreads);
322  std::fill(m_ThreadSecondOrderAccumulators.begin(), m_ThreadSecondOrderAccumulators.end(), zeroMatrix);
323 
324  RealType zeroReal = itk::NumericTraits<RealType>::ZeroValue();
325  m_ThreadSecondOrderComponentAccumulators.resize(numberOfThreads);
327  }
328 
330  {
331  m_IgnoredInfinitePixelCount = std::vector<unsigned int>(numberOfThreads, 0);
332  }
333 
335  {
336  m_IgnoredUserPixelCount = std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
337  }
338 }
339 
340 template <class TInputImage, class TPrecision>
342 {
343  TInputImage* inputPtr = const_cast<TInputImage*>(this->GetInput());
344  const unsigned int nbPixels = inputPtr->GetLargestPossibleRegion().GetNumberOfPixels();
345  const unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
346 
347  PixelType minimum;
348  minimum.SetSize(numberOfComponent);
349  minimum.Fill(itk::NumericTraits<InternalPixelType>::max());
350  PixelType maximum;
351  maximum.SetSize(numberOfComponent);
352  maximum.Fill(itk::NumericTraits<InternalPixelType>::NonpositiveMin());
353 
354  RealPixelType streamFirstOrderAccumulator(numberOfComponent);
355  streamFirstOrderAccumulator.Fill(itk::NumericTraits<PrecisionType>::Zero);
356  MatrixType streamSecondOrderAccumulator(numberOfComponent, numberOfComponent);
357  streamSecondOrderAccumulator.Fill(itk::NumericTraits<PrecisionType>::Zero);
358 
359  RealType streamFirstOrderComponentAccumulator = itk::NumericTraits<RealType>::Zero;
360  RealType streamSecondOrderComponentAccumulator = itk::NumericTraits<RealType>::Zero;
361 
362  unsigned int ignoredInfinitePixelCount = 0;
363  unsigned int ignoredUserPixelCount = 0;
364 
365  // Accumulate results from all threads
366  const itk::ThreadIdType numberOfThreads = this->GetNumberOfThreads();
367  for (itk::ThreadIdType threadId = 0; threadId < numberOfThreads; ++threadId)
368  {
369  if (m_EnableMinMax)
370  {
371  const PixelType& threadMin = m_ThreadMin[threadId];
372  const PixelType& threadMax = m_ThreadMax[threadId];
373 
374  for (unsigned int j = 0; j < numberOfComponent; ++j)
375  {
376  if (threadMin[j] < minimum[j])
377  {
378  minimum[j] = threadMin[j];
379  }
380  if (threadMax[j] > maximum[j])
381  {
382  maximum[j] = threadMax[j];
383  }
384  }
385  }
386 
388  {
389  streamFirstOrderAccumulator += m_ThreadFirstOrderAccumulators[threadId];
390  streamFirstOrderComponentAccumulator += m_ThreadFirstOrderComponentAccumulators[threadId];
391  }
392 
394  {
395  streamSecondOrderAccumulator += m_ThreadSecondOrderAccumulators[threadId];
396  streamSecondOrderComponentAccumulator += m_ThreadSecondOrderComponentAccumulators[threadId];
397  }
398  // Ignored Infinite Pixels
399  ignoredInfinitePixelCount += m_IgnoredInfinitePixelCount[threadId];
400  // Ignored Pixels
401  ignoredUserPixelCount += m_IgnoredUserPixelCount[threadId];
402  }
403 
404  // There cannot be more ignored pixels than read pixels.
405  assert(nbPixels >= ignoredInfinitePixelCount + ignoredUserPixelCount);
406  if (nbPixels < ignoredInfinitePixelCount + ignoredUserPixelCount)
407  {
408  itkExceptionMacro("nbPixels < ignoredInfinitePixelCount + ignoredUserPixelCount");
409  }
410 
411  unsigned int nbRelevantPixel = nbPixels - (ignoredInfinitePixelCount + ignoredUserPixelCount);
412 
413  CountType nbRelevantPixels(numberOfComponent);
414  nbRelevantPixels.Fill(nbRelevantPixel);
415 
416  this->GetNbRelevantPixelsOutput()->Set(nbRelevantPixels);
417 
418  if (nbRelevantPixel == 0)
419  {
420  itkExceptionMacro("Statistics cannot be calculated with zero relevant pixels.");
421  }
422 
423  // Final calculations
424  if (m_EnableMinMax)
425  {
426  this->GetMinimumOutput()->Set(minimum);
427  this->GetMaximumOutput()->Set(maximum);
428  }
429 
431  {
432  this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
433 
434  this->GetMeanOutput()->Set(streamFirstOrderAccumulator / nbRelevantPixel);
435  this->GetSumOutput()->Set(streamFirstOrderAccumulator);
436  }
437 
439  {
440  MatrixType cor = streamSecondOrderAccumulator / nbRelevantPixel;
441  this->GetCorrelationOutput()->Set(cor);
442 
443  const RealPixelType& mean = this->GetMeanOutput()->Get();
444 
445  double regul = 1.0;
446  double regulComponent = 1.0;
447 
448  if (m_UseUnbiasedEstimator && nbRelevantPixel > 1)
449  {
450  regul = static_cast<double>(nbRelevantPixel) / (static_cast<double>(nbRelevantPixel) - 1.0);
451  }
452 
453  if (m_UseUnbiasedEstimator && (nbRelevantPixel * numberOfComponent) > 1)
454  {
455  regulComponent = static_cast<double>(nbRelevantPixel * numberOfComponent) / (static_cast<double>(nbRelevantPixel * numberOfComponent) - 1.0);
456  }
457 
458  MatrixType cov = cor;
459  for (unsigned int r = 0; r < numberOfComponent; ++r)
460  {
461  for (unsigned int c = 0; c < numberOfComponent; ++c)
462  {
463  cov(r, c) = regul * (cov(r, c) - mean[r] * mean[c]);
464  }
465  }
466  this->GetCovarianceOutput()->Set(cov);
467 
468  this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
469  this->GetComponentCorrelationOutput()->Set(streamSecondOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
470  this->GetComponentCovarianceOutput()->Set(regulComponent * (this->GetComponentCorrelation() - (this->GetComponentMean() * this->GetComponentMean())));
471  }
472 }
473 
474 template <class TInputImage, class TPrecision>
476  itk::ThreadIdType threadId)
477 {
478  // Support progress methods/callbacks
479  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
480 
481  // Grab the input
482  InputImagePointer inputPtr = const_cast<TInputImage*>(this->GetInput());
483  PixelType& threadMin = m_ThreadMin[threadId];
484  PixelType& threadMax = m_ThreadMax[threadId];
485 
486 
487  itk::ImageRegionConstIteratorWithIndex<TInputImage> it(inputPtr, outputRegionForThread);
488 
489  for (it.GoToBegin(); !it.IsAtEnd(); ++it, progress.CompletedPixel())
490  {
491  const PixelType& vectorValue = it.Get();
492 
493  float finiteProbe = 0.;
494  bool userProbe = m_IgnoreUserDefinedValue;
495  for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
496  {
497  finiteProbe += (float)(vectorValue[j]);
498  userProbe = userProbe && (vectorValue[j] == m_UserIgnoredValue);
499  }
500 
501  if (m_IgnoreInfiniteValues && !(vnl_math_isfinite(finiteProbe)))
502  {
503  m_IgnoredInfinitePixelCount[threadId]++;
504  }
505  else
506  {
507  if (userProbe)
508  {
509  m_IgnoredUserPixelCount[threadId]++;
510  }
511  else
512  {
513  if (m_EnableMinMax)
514  {
515  for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
516  {
517  if (vectorValue[j] < threadMin[j])
518  {
519  threadMin[j] = vectorValue[j];
520  }
521  if (vectorValue[j] > threadMax[j])
522  {
523  threadMax[j] = vectorValue[j];
524  }
525  }
526  }
527 
529  {
530  RealPixelType& threadFirstOrder = m_ThreadFirstOrderAccumulators[threadId];
531  RealType& threadFirstOrderComponent = m_ThreadFirstOrderComponentAccumulators[threadId];
532 
533  threadFirstOrder += vectorValue;
534 
535  for (unsigned int i = 0; i < vectorValue.GetSize(); ++i)
536  {
537  threadFirstOrderComponent += vectorValue[i];
538  }
539  }
540 
542  {
543  MatrixType& threadSecondOrder = m_ThreadSecondOrderAccumulators[threadId];
544  RealType& threadSecondOrderComponent = m_ThreadSecondOrderComponentAccumulators[threadId];
545 
546  for (unsigned int r = 0; r < threadSecondOrder.Rows(); ++r)
547  {
548  for (unsigned int c = 0; c < threadSecondOrder.Cols(); ++c)
549  {
550  threadSecondOrder(r, c) += static_cast<PrecisionType>(vectorValue[r]) * static_cast<PrecisionType>(vectorValue[c]);
551  }
552  }
553  threadSecondOrderComponent += vectorValue.GetSquaredNorm();
554  }
555  }
556  }
557  }
558 }
559 
560 template <class TImage, class TPrecision>
562 {
563  Superclass::PrintSelf(os, indent);
564 
565  os << indent << "Min: " << this->GetMinimumOutput()->Get() << std::endl;
566  os << indent << "Max: " << this->GetMaximumOutput()->Get() << std::endl;
567  os << indent << "Mean: " << this->GetMeanOutput()->Get() << std::endl;
568  os << indent << "Covariance: " << this->GetCovarianceOutput()->Get() << std::endl;
569  os << indent << "Correlation: " << this->GetCorrelationOutput()->Get() << std::endl;
570  os << indent << "Relevant pixel: " << this->GetNbRelevantPixelsOutput()->Get() << std::endl;
571  os << indent << "Component Mean: " << this->GetComponentMeanOutput()->Get() << std::endl;
572  os << indent << "Component Covariance: " << this->GetComponentCovarianceOutput()->Get() << std::endl;
573  os << indent << "Component Correlation: " << this->GetComponentCorrelationOutput()->Get() << std::endl;
574  os << indent << "UseUnbiasedEstimator: " << (this->m_UseUnbiasedEstimator ? "true" : "false") << std::endl;
575 }
576 
577 } // end namespace otb
578 #endif
void PrintSelf(std::ostream &os, itk::Indent indent) const override
void PrintSelf(std::ostream &os, itk::Indent indent) const override
itk::SimpleDataObjectDecorator< RealPixelType > RealPixelObjectType
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
The "otb" namespace contains all Orfeo Toolbox (OTB) classes.
itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType