Merge branch 'fix-nodata-statistics' into 'develop'

BUG: Fix nodata handling in StreamingStatisticsMapFromLabelImageFilter

Closes #2019

See merge request orfeotoolbox/otb!674

Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.h

@@ -67,22 +67,29 @@ public:
     : m_NoDataValue(noDataValue), m_Count(1), m_UseNoDataValue(useNoDataValue)
   {
     m_Count = 1;
-    m_Sum   = pixel;
-    m_Min   = pixel;
-    m_Max   = pixel;
-    m_SqSum = pixel;
     m_BandCount.SetSize(pixel.GetSize());
-    for (unsigned int band = 0; band < m_SqSum.GetSize(); band++)
+    m_Sum.SetSize(pixel.GetSize());
+    m_Min.SetSize(pixel.GetSize());
+    m_Max.SetSize(pixel.GetSize());
+    m_SqSum.SetSize(pixel.GetSize());
+    for (unsigned int band = 0; band < pixel.GetSize(); band++)
     {
       auto val = pixel[band];
       if (!m_UseNoDataValue || val != m_NoDataValue)
       {
         m_BandCount[band] = 1;
-        m_SqSum[band] *= m_SqSum[band];
+        m_Sum[band]       = val;
+        m_Min[band]       = val;
+        m_Max[band]       = val;
+        m_SqSum[band]     = val * val;
       }
       else
       {
         m_BandCount[band] = 0;
+        m_Sum[band]       = itk::NumericTraits<RealValueType>::ZeroValue();
+        m_Min[band]       = itk::NumericTraits<RealValueType>::max();
+        m_Max[band]       = itk::NumericTraits<RealValueType>::min();
+        m_SqSum[band]     = itk::NumericTraits<RealValueType>::ZeroValue();
       }
     }
   }
@@ -97,11 +104,10 @@ public:
       const RealValueType value   = pixel[band];
       const RealValueType sqValue = value * value;
 
-      if(!m_UseNoDataValue || value != m_NoDataValue)
-        {
-        UpdateValues(1, value, sqValue, value, value, m_BandCount[band], m_Sum[band], m_SqSum[band], m_Min[band],
-                   m_Max[band]);
-        }
+      if (!m_UseNoDataValue || value != m_NoDataValue)
+      {
+        UpdateValues(1, value, sqValue, value, value, m_BandCount[band], m_Sum[band], m_SqSum[band], m_Min[band], m_Max[band]);
+      }
     }
   }
 
@@ -172,10 +178,10 @@ class ITK_EXPORT PersistentStreamingStatisticsMapFromLabelImageFilter : public P
 {
 public:
   /** Standard Self typedef */
-  typedef PersistentStreamingStatisticsMapFromLabelImageFilter Self;
+  typedef PersistentStreamingStatisticsMapFromLabelImageFilter        Self;
   typedef PersistentImageFilter<TInputVectorImage, TInputVectorImage> Superclass;
-  typedef itk::SmartPointer<Self>       Pointer;
-  typedef itk::SmartPointer<const Self> ConstPointer;
+  typedef itk::SmartPointer<Self>                                     Pointer;
+  typedef itk::SmartPointer<const Self>                               ConstPointer;
 
   /** Method for creation through the object factory. */
   itkNewMacro(Self);
@@ -189,14 +195,14 @@ public:
   typedef TLabelImage                         LabelImageType;
   typedef typename TLabelImage::Pointer       LabelImagePointer;
 
-  typedef typename VectorImageType::RegionType           RegionType;
-  typedef typename VectorImageType::PixelType            VectorPixelType;
-  typedef typename VectorImageType::PixelType::ValueType VectorPixelValueType;
-  typedef typename LabelImageType::PixelType             LabelPixelType;
-  typedef itk::VariableLengthVector<double>              RealVectorPixelType;
-  typedef StatisticsAccumulator<RealVectorPixelType>     AccumulatorType;
-  typedef std::unordered_map<LabelPixelType, AccumulatorType> AccumulatorMapType;
-  typedef std::vector<AccumulatorMapType> AccumulatorMapCollectionType;
+  typedef typename VectorImageType::RegionType                    RegionType;
+  typedef typename VectorImageType::PixelType                     VectorPixelType;
+  typedef typename VectorImageType::PixelType::ValueType          VectorPixelValueType;
+  typedef typename LabelImageType::PixelType                      LabelPixelType;
+  typedef itk::VariableLengthVector<double>                       RealVectorPixelType;
+  typedef StatisticsAccumulator<RealVectorPixelType>              AccumulatorType;
+  typedef std::unordered_map<LabelPixelType, AccumulatorType>     AccumulatorMapType;
+  typedef std::vector<AccumulatorMapType>                         AccumulatorMapCollectionType;
   typedef std::unordered_map<LabelPixelType, RealVectorPixelType> PixelValueMapType;
   typedef std::unordered_map<LabelPixelType, double>              LabelPopulationMapType;
 
@@ -335,14 +341,14 @@ private:
 
 template <class TInputVectorImage, class TLabelImage>
 class ITK_EXPORT StreamingStatisticsMapFromLabelImageFilter
-    : public PersistentFilterStreamingDecorator<PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelImage>>
+  : public PersistentFilterStreamingDecorator<PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelImage>>
 {
 public:
   /** Standard Self typedef */
-  typedef StreamingStatisticsMapFromLabelImageFilter Self;
+  typedef StreamingStatisticsMapFromLabelImageFilter                                                                               Self;
   typedef PersistentFilterStreamingDecorator<PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelImage>> Superclass;
-  typedef itk::SmartPointer<Self>       Pointer;
-  typedef itk::SmartPointer<const Self> ConstPointer;
+  typedef itk::SmartPointer<Self>                                                                                                  Pointer;
+  typedef itk::SmartPointer<const Self>                                                                                            ConstPointer;
 
   /** Type macro */
   itkNewMacro(Self);

Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.hxx

@@ -28,6 +28,7 @@
 #include "itkProgressReporter.h"
 #include "otbMacro.h"
 #include <cmath>
+#include <utility>
 
 namespace otb
 {
@@ -168,6 +169,8 @@ void PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLa
     // Mean & stdev
     RealVectorPixelType mean(sum);
     RealVectorPixelType std(sqSum);
+    RealVectorPixelType min(it.second.GetMin());
+    RealVectorPixelType max(it.second.GetMax());
     for (unsigned int band = 0; band < mean.GetSize(); band++)
     {
       // Number of valid pixels in band
@@ -178,13 +181,20 @@ void PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLa
       // Unbiased standard deviation (not sure unbiased is useful here)
       const double variance = (sqSum[band] - (sum[band] * mean[band])) / (count - 1);
       std[band]             = std::sqrt(variance);
+
+      // Use the no data value when no valid pixels were found
+      if (this->GetUseNoDataValue() && count == 0)
+      {
+        min[band] = this->GetNoDataValue();
+        max[band] = this->GetNoDataValue();
+      }
     }
-    m_MeanRadiometricValue[label]  = mean;
-    m_StDevRadiometricValue[label] = std;
+    m_MeanRadiometricValue.emplace(label, std::move(mean));
+    m_StDevRadiometricValue.emplace(label, std::move(std));
 
     // Min & max
-    m_MinRadiometricValue[label] = it.second.GetMin();
-    m_MaxRadiometricValue[label] = it.second.GetMax();
+    m_MinRadiometricValue.emplace(label, std::move(min));
+    m_MaxRadiometricValue.emplace(label, std::move(max));
   }
 }