OTB  6.7.0
Orfeo Toolbox
otbStreamingStatisticsVectorImageFilter.hxx
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2005-2019 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"
27 #include "itkProgressReporter.h"
28 #include "otbMacro.h"
29 
30 namespace otb
31 {
32 
33 template<class TInputImage, class TPrecision>
36  : m_EnableMinMax(true),
37  m_EnableFirstOrderStats(true),
38  m_EnableSecondOrderStats(true),
39  m_UseUnbiasedEstimator(true),
40  m_IgnoreInfiniteValues(true),
41  m_IgnoreUserDefinedValue(false),
42  m_UserIgnoredValue(itk::NumericTraits<InternalPixelType>::Zero)
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 vector/matrix types
49  for (unsigned int i = 1; i < 11; ++i)
50  {
51  this->itk::ProcessObject::SetNthOutput(i, this->MakeOutput(i).GetPointer());
52  }
53  // Initiate ignored pixel counters
54  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
55  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
56 }
57 
58 template<class TInputImage, class TPrecision>
62 {
63  switch (output)
64  {
65  case 0:
66  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
67  break;
68  case 1:
69  case 2:
70  // min/max
71  return static_cast<itk::DataObject*>(PixelObjectType::New().GetPointer());
72  break;
73  case 3:
74  case 4:
75  // mean / sum
76  return static_cast<itk::DataObject*>(RealPixelObjectType::New().GetPointer());
77  break;
78  case 5:
79  case 6:
80  // covariance / correlation
81  return static_cast<itk::DataObject*>(MatrixObjectType::New().GetPointer());
82  break;
83  case 7:
84  case 8:
85  case 9:
86  // component mean, component covariance, component correlation
87  return static_cast<itk::DataObject*>(RealObjectType::New().GetPointer());
88  break;
89  case 10:
90  // relevant pixel
91  return static_cast<itk::DataObject*>(CountObjectType::New().GetPointer());
92  default:
93  // might as well make an image
94  return static_cast<itk::DataObject*>(TInputImage::New().GetPointer());
95  break;
96  }
97 }
98 
99 template<class TInputImage, class TPrecision>
103 {
104  return static_cast<CountObjectType*>(this->itk::ProcessObject::GetOutput(10));
105 }
106 
107 template<class TInputImage, class TPrecision>
111 {
112  return static_cast<const CountObjectType*>(this->itk::ProcessObject::GetOutput(10));
113 }
114 
115 
116 template<class TInputImage, class TPrecision>
120 {
121  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
122 }
123 
124 template<class TInputImage, class TPrecision>
128 {
129  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(1));
130 }
131 
132 template<class TInputImage, class TPrecision>
136 {
137  return static_cast<PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
138 }
139 
140 template<class TInputImage, class TPrecision>
144 {
145  return static_cast<const PixelObjectType*>(this->itk::ProcessObject::GetOutput(2));
146 }
147 
148 template<class TInputImage, class TPrecision>
152 {
153  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(7));
154 }
155 
156 template<class TInputImage, class TPrecision>
160 {
161  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(7));
162 }
163 
164 template<class TInputImage, class TPrecision>
168 {
169  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(8));
170 }
171 
172 template<class TInputImage, class TPrecision>
176 {
177  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(8));
178 }
179 
180 
181 template<class TInputImage, class TPrecision>
185 {
186  return static_cast<RealObjectType*>(this->itk::ProcessObject::GetOutput(9));
187 }
188 
189 template<class TInputImage, class TPrecision>
193 {
194  return static_cast<const RealObjectType*>(this->itk::ProcessObject::GetOutput(9));
195 }
196 
197 template<class TInputImage, class TPrecision>
201 {
202  return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
203 }
204 
205 template<class TInputImage, class TPrecision>
209 {
210  return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(3));
211 }
212 
213 template<class TInputImage, class TPrecision>
217 {
218  return static_cast<RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
219 }
220 
221 template<class TInputImage, class TPrecision>
225 {
226  return static_cast<const RealPixelObjectType*>(this->itk::ProcessObject::GetOutput(4));
227 }
228 
229 template<class TInputImage, class TPrecision>
233 {
234  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
235 }
236 
237 template<class TInputImage, class TPrecision>
241 {
242  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(5));
243 }
244 
245 template<class TInputImage, class TPrecision>
249 {
250  return static_cast<MatrixObjectType*>(this->itk::ProcessObject::GetOutput(6));
251 }
252 
253 template<class TInputImage, class TPrecision>
257 {
258  return static_cast<const MatrixObjectType*>(this->itk::ProcessObject::GetOutput(6));
259 }
260 
261 template<class TInputImage, class TPrecision>
262 void
265 {
266  Superclass::GenerateOutputInformation();
267  if (this->GetInput())
268  {
269  this->GetOutput()->CopyInformation(this->GetInput());
270  this->GetOutput()->SetLargestPossibleRegion(this->GetInput()->GetLargestPossibleRegion());
271 
272  if (this->GetOutput()->GetRequestedRegion().GetNumberOfPixels() == 0)
273  {
274  this->GetOutput()->SetRequestedRegion(this->GetOutput()->GetLargestPossibleRegion());
275  }
276  }
277 }
278 
279 template<class TInputImage, class TPrecision>
280 void
283 {
284  // This is commented to prevent the streaming of the whole image for the first stream strip
285  // It shall not cause any problem because the output image of this filter is not intended to be used.
286  //InputImagePointer image = const_cast< TInputImage * >( this->GetInput() );
287  //this->GraftOutput( image );
288  // Nothing that needs to be allocated for the remaining outputs
289 }
290 
291 template<class TInputImage, class TPrecision>
292 void
295 {
296  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
297  inputPtr->UpdateOutputInformation();
298 
299  unsigned int numberOfThreads = this->GetNumberOfThreads();
300  unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
301 
302  if (m_EnableMinMax)
303  {
304  PixelType tempPixel;
305  tempPixel.SetSize(numberOfComponent);
306 
308  this->GetMinimumOutput()->Set(tempPixel);
309 
311  this->GetMaximumOutput()->Set(tempPixel);
312 
313  PixelType tempTemporiesPixel;
314  tempTemporiesPixel.SetSize(numberOfComponent);
315  tempTemporiesPixel.Fill(itk::NumericTraits<InternalPixelType>::max());
316  m_ThreadMin = std::vector<PixelType>(numberOfThreads, tempTemporiesPixel);
317 
319  m_ThreadMax = std::vector<PixelType>(numberOfThreads, tempTemporiesPixel);
320  }
321 
322  if (m_EnableSecondOrderStats)
323  {
324  m_EnableFirstOrderStats = true;
325  }
326 
327  if (m_EnableFirstOrderStats)
328  {
329  RealPixelType zeroRealPixel;
330  zeroRealPixel.SetSize(numberOfComponent);
332  this->GetMeanOutput()->Set(zeroRealPixel);
333  this->GetSumOutput()->Set(zeroRealPixel);
334  m_ThreadFirstOrderAccumulators.resize(numberOfThreads);
335  std::fill(m_ThreadFirstOrderAccumulators.begin(), m_ThreadFirstOrderAccumulators.end(), zeroRealPixel);
336 
338  m_ThreadFirstOrderComponentAccumulators.resize(numberOfThreads);
339  std::fill(m_ThreadFirstOrderComponentAccumulators.begin(), m_ThreadFirstOrderComponentAccumulators.end(), zeroReal);
340 
341  }
342 
343  if (m_EnableSecondOrderStats)
344  {
345  MatrixType zeroMatrix;
346  zeroMatrix.SetSize(numberOfComponent, numberOfComponent);
348  this->GetCovarianceOutput()->Set(zeroMatrix);
349  this->GetCorrelationOutput()->Set(zeroMatrix);
350 
351  m_ThreadSecondOrderAccumulators.resize(numberOfThreads);
352  std::fill(m_ThreadSecondOrderAccumulators.begin(), m_ThreadSecondOrderAccumulators.end(), zeroMatrix);
353 
355  m_ThreadSecondOrderComponentAccumulators.resize(numberOfThreads);
356  std::fill(m_ThreadSecondOrderComponentAccumulators.begin(), m_ThreadSecondOrderComponentAccumulators.end(), zeroReal);
357  }
358 
359  if (m_IgnoreInfiniteValues)
360  {
361  m_IgnoredInfinitePixelCount= std::vector<unsigned int>(numberOfThreads, 0);
362  }
363 
364  if (m_IgnoreUserDefinedValue)
365  {
366  m_IgnoredUserPixelCount= std::vector<unsigned int>(this->GetNumberOfThreads(), 0);
367  }
368 }
369 
370 template<class TInputImage, class TPrecision>
371 void
374 {
375  TInputImage * inputPtr = const_cast<TInputImage *>(this->GetInput());
376  const unsigned int nbPixels = inputPtr->GetLargestPossibleRegion().GetNumberOfPixels();
377  const unsigned int numberOfComponent = inputPtr->GetNumberOfComponentsPerPixel();
378 
379  PixelType minimum;
380  minimum.SetSize(numberOfComponent);
382  PixelType maximum;
383  maximum.SetSize(numberOfComponent);
385 
386  RealPixelType streamFirstOrderAccumulator(numberOfComponent);
387  streamFirstOrderAccumulator.Fill(itk::NumericTraits<PrecisionType>::Zero);
388  MatrixType streamSecondOrderAccumulator(numberOfComponent, numberOfComponent);
389  streamSecondOrderAccumulator.Fill(itk::NumericTraits<PrecisionType>::Zero);
390 
391  RealType streamFirstOrderComponentAccumulator = itk::NumericTraits<RealType>::Zero;
392  RealType streamSecondOrderComponentAccumulator = itk::NumericTraits<RealType>::Zero;
393 
394  unsigned int ignoredInfinitePixelCount = 0;
395  unsigned int ignoredUserPixelCount = 0;
396 
397  // Accumulate results from all threads
398  const itk::ThreadIdType numberOfThreads = this->GetNumberOfThreads();
399  for (itk::ThreadIdType threadId = 0; threadId < numberOfThreads; ++threadId)
400  {
401  if (m_EnableMinMax)
402  {
403  const PixelType& threadMin = m_ThreadMin [threadId];
404  const PixelType& threadMax = m_ThreadMax [threadId];
405 
406  for (unsigned int j = 0; j < numberOfComponent; ++j)
407  {
408  if (threadMin[j] < minimum[j])
409  {
410  minimum[j] = threadMin[j];
411  }
412  if (threadMax[j] > maximum[j])
413  {
414  maximum[j] = threadMax[j];
415  }
416  }
417  }
418 
419  if (m_EnableFirstOrderStats)
420  {
421  streamFirstOrderAccumulator += m_ThreadFirstOrderAccumulators[threadId];
422  streamFirstOrderComponentAccumulator += m_ThreadFirstOrderComponentAccumulators[threadId];
423  }
424 
425  if (m_EnableSecondOrderStats)
426  {
427  streamSecondOrderAccumulator += m_ThreadSecondOrderAccumulators[threadId];
428  streamSecondOrderComponentAccumulator += m_ThreadSecondOrderComponentAccumulators[threadId];
429  }
430  // Ignored Infinite Pixels
431  ignoredInfinitePixelCount += m_IgnoredInfinitePixelCount[threadId];
432  // Ignored Pixels
433  ignoredUserPixelCount += m_IgnoredUserPixelCount[threadId];
434  }
435 
436  // There cannot be more ignored pixels than read pixels.
437  assert( nbPixels >= ignoredInfinitePixelCount + ignoredUserPixelCount );
438  if( nbPixels < ignoredInfinitePixelCount + ignoredUserPixelCount )
439  {
440  itkExceptionMacro(
441  "nbPixels < ignoredInfinitePixelCount + ignoredUserPixelCount"
442  );
443  }
444 
445  unsigned int nbRelevantPixel =
446  nbPixels - (ignoredInfinitePixelCount + ignoredUserPixelCount);
447 
448  CountType nbRelevantPixels(numberOfComponent);
449  nbRelevantPixels.Fill(nbRelevantPixel);
450 
451  this->GetNbRelevantPixelsOutput()->Set(nbRelevantPixels);
452 
453  if( nbRelevantPixel==0 )
454  {
455  itkExceptionMacro(
456  "Statistics cannot be calculated with zero relevant pixels."
457  );
458  }
459 
460  // Final calculations
461  if (m_EnableMinMax)
462  {
463  this->GetMinimumOutput()->Set(minimum);
464  this->GetMaximumOutput()->Set(maximum);
465  }
466 
467  if (m_EnableFirstOrderStats)
468  {
469  this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
470 
471  this->GetMeanOutput()->Set(streamFirstOrderAccumulator / nbRelevantPixel);
472  this->GetSumOutput()->Set(streamFirstOrderAccumulator);
473  }
474 
475  if (m_EnableSecondOrderStats)
476  {
477  MatrixType cor = streamSecondOrderAccumulator / nbRelevantPixel;
478  this->GetCorrelationOutput()->Set(cor);
479 
480  const RealPixelType& mean = this->GetMeanOutput()->Get();
481 
482  double regul = 1.0;
483  double regulComponent = 1.0;
484 
485  if( m_UseUnbiasedEstimator && nbRelevantPixel>1 )
486  {
487  regul =
488  static_cast< double >( nbRelevantPixel ) /
489  ( static_cast< double >( nbRelevantPixel ) - 1.0 );
490  }
491 
492  if( m_UseUnbiasedEstimator && (nbRelevantPixel * numberOfComponent) > 1 )
493  {
494  regulComponent =
495  static_cast< double >(nbRelevantPixel * numberOfComponent) /
496  ( static_cast< double >(nbRelevantPixel * numberOfComponent) - 1.0 );
497  }
498 
499  MatrixType cov = cor;
500  for (unsigned int r = 0; r < numberOfComponent; ++r)
501  {
502  for (unsigned int c = 0; c < numberOfComponent; ++c)
503  {
504  cov(r, c) = regul * (cov(r, c) - mean[r] * mean[c]);
505  }
506  }
507  this->GetCovarianceOutput()->Set(cov);
508 
509  this->GetComponentMeanOutput()->Set(streamFirstOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
510  this->GetComponentCorrelationOutput()->Set(streamSecondOrderComponentAccumulator / (nbRelevantPixel * numberOfComponent));
511  this->GetComponentCovarianceOutput()->Set(
512  regulComponent * (this->GetComponentCorrelation()
513  - (this->GetComponentMean() * this->GetComponentMean())));
514  }
515 }
516 
517 template<class TInputImage, class TPrecision>
518 void
520 ::ThreadedGenerateData(const RegionType& outputRegionForThread, itk::ThreadIdType threadId)
521  {
522  // Support progress methods/callbacks
523  itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
524 
525  // Grab the input
526  InputImagePointer inputPtr = const_cast<TInputImage *>(this->GetInput());
527  PixelType& threadMin = m_ThreadMin [threadId];
528  PixelType& threadMax = m_ThreadMax [threadId];
529 
530 
531  itk::ImageRegionConstIteratorWithIndex<TInputImage> it(inputPtr, outputRegionForThread);
532 
533  for (it.GoToBegin(); !it.IsAtEnd(); ++it, progress.CompletedPixel())
534  {
535  const PixelType& vectorValue = it.Get();
536 
537  float finiteProbe = 0.;
538  bool userProbe = m_IgnoreUserDefinedValue;
539  for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
540  {
541  finiteProbe += (float)(vectorValue[j]);
542  userProbe = userProbe && (vectorValue[j] == m_UserIgnoredValue);
543  }
544 
545  if (m_IgnoreInfiniteValues && !(vnl_math_isfinite(finiteProbe)))
546  {
547  m_IgnoredInfinitePixelCount[threadId] ++;
548  }
549  else
550  {
551  if (userProbe)
552  {
553  m_IgnoredUserPixelCount[threadId] ++;
554  }
555  else
556  {
557  if (m_EnableMinMax)
558  {
559  for (unsigned int j = 0; j < vectorValue.GetSize(); ++j)
560  {
561  if (vectorValue[j] < threadMin[j])
562  {
563  threadMin[j] = vectorValue[j];
564  }
565  if (vectorValue[j] > threadMax[j])
566  {
567  threadMax[j] = vectorValue[j];
568  }
569  }
570  }
571 
572  if (m_EnableFirstOrderStats)
573  {
574  RealPixelType& threadFirstOrder = m_ThreadFirstOrderAccumulators [threadId];
575  RealType& threadFirstOrderComponent = m_ThreadFirstOrderComponentAccumulators [threadId];
576 
577  threadFirstOrder += vectorValue;
578 
579  for (unsigned int i = 0; i < vectorValue.GetSize(); ++i)
580  {
581  threadFirstOrderComponent += vectorValue[i];
582  }
583  }
584 
585  if (m_EnableSecondOrderStats)
586  {
587  MatrixType& threadSecondOrder = m_ThreadSecondOrderAccumulators[threadId];
588  RealType& threadSecondOrderComponent = m_ThreadSecondOrderComponentAccumulators[threadId];
589 
590  for (unsigned int r = 0; r < threadSecondOrder.Rows(); ++r)
591  {
592  for (unsigned int c = 0; c < threadSecondOrder.Cols(); ++c)
593  {
594  threadSecondOrder(r, c) += static_cast<PrecisionType>(vectorValue[r]) * static_cast<PrecisionType>(vectorValue[c]);
595  }
596  }
597  threadSecondOrderComponent += vectorValue.GetSquaredNorm();
598  }
599  }
600  }
601  }
602 
603  }
604 
605 template <class TImage, class TPrecision>
606 void
608 ::PrintSelf(std::ostream& os, itk::Indent indent) const
609 {
610  Superclass::PrintSelf(os, indent);
611 
612  os << indent << "Min: " << this->GetMinimumOutput()->Get() << std::endl;
613  os << indent << "Max: " << this->GetMaximumOutput()->Get() << std::endl;
614  os << indent << "Mean: " << this->GetMeanOutput()->Get() << std::endl;
615  os << indent << "Covariance: " << this->GetCovarianceOutput()->Get() << std::endl;
616  os << indent << "Correlation: " << this->GetCorrelationOutput()->Get() << std::endl;
617  os << indent << "Relevant pixel: " << this->GetNbRelevantPixelsOutput()->Get() << std::endl;
618  os << indent << "Component Mean: " << this->GetComponentMeanOutput()->Get() << std::endl;
619  os << indent << "Component Covariance: " << this->GetComponentCovarianceOutput()->Get() << std::endl;
620  os << indent << "Component Correlation: " << this->GetComponentCorrelationOutput()->Get() << std::endl;
621  os << indent << "UseUnbiasedEstimator: " << (this->m_UseUnbiasedEstimator ? "true" : "false") << std::endl;
622 }
623 
624 } // end namespace otb
625 #endif
unsigned int Rows() const
Superclass::DataObjectPointerArraySizeType DataObjectPointerArraySizeType
PixelType Get(void) const
static T ZeroValue()
unsigned int Cols() const
void PrintSelf(std::ostream &os, itk::Indent indent) const override
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override
void ThreadedGenerateData(const RegionType &outputRegionForThread, itk::ThreadIdType threadId) override
bool SetSize(unsigned int r, unsigned int c)
void Fill(TValue const &v) noexcept
unsigned int ThreadIdType
void Fill(const T &value)
void SetSize(unsigned int sz, TReallocatePolicy reallocatePolicy, TKeepValuesPolicy keepValues)
virtual void SetNthOutput(DataObjectPointerArraySizeType num, DataObject *output)
DataObject * GetOutput(const DataObjectIdentifierType &key)