Commit 31c5dbb5 authored by Laurențiu Nicola's avatar Laurențiu Nicola

Merge branch 'streaming-statistics-map-nodata' into 'develop'

StreamingStatisticsMapFromLabelImageFilter no data support

See merge request orfeotoolbox/otb!223
parents fd6a1994 d914da3a
......@@ -54,36 +54,42 @@ class StatisticsAccumulator
{
public:
typedef typename TRealVectorPixelType::ValueType RealValueType;
typedef typename TRealVectorPixelType::ValueType RealValueType;
typedef uint64_t PixelCountType;
typedef itk::VariableLengthVector<PixelCountType> PixelCountVectorType;
// Constructor (default)
StatisticsAccumulator(){}
StatisticsAccumulator() : m_Count(), m_NoDataValue(), m_UseNoDataValue() {}
// Constructor (initialize the accumulator with the given pixel)
StatisticsAccumulator(const TRealVectorPixelType & pixel)
StatisticsAccumulator(RealValueType noDataValue,
bool useNoDataValue,
const TRealVectorPixelType & pixel)
: 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_SqSum[band] *= m_SqSum[band];
}
// Constructor (other)
StatisticsAccumulator(const StatisticsAccumulator & other)
{
m_Count = other.m_Count;
m_Sum = other.m_Sum;
m_Min = other.m_Min;
m_Max = other.m_Max;
m_SqSum = other.m_SqSum;
{
auto val = pixel[band];
if (!m_UseNoDataValue || val != m_NoDataValue)
{
m_BandCount[band] = 1;
m_SqSum[band] *= m_SqSum[band];
}
else
{
m_BandCount[band] = 0;
}
}
}
// Destructor
~StatisticsAccumulator(){}
// Function update (pixel)
void Update(const TRealVectorPixelType & pixel)
{
......@@ -93,8 +99,13 @@ public:
{
const RealValueType value = pixel[band];
const RealValueType sqValue = value * value;
UpdateValues(value, sqValue, value, value,
m_Sum[band], m_SqSum[band], m_Min[band], m_Max[band]);
UpdateValues(!m_UseNoDataValue || value != m_NoDataValue,
value, sqValue,
value, value,
m_BandCount[band],
m_Sum[band], m_SqSum[band],
m_Min[band], m_Max[band]);
}
}
......@@ -105,12 +116,17 @@ public:
const unsigned int nBands = other.m_Sum.GetSize();
for (unsigned int band = 0 ; band < nBands ; band ++ )
{
UpdateValues(other.m_Sum[band], other.m_SqSum[band], other.m_Min[band], other.m_Max[band],
m_Sum[band], m_SqSum[band], m_Min[band], m_Max[band]);
UpdateValues(other.m_BandCount[band],
other.m_Sum[band], other.m_SqSum[band],
other.m_Min[band], other.m_Max[band],
m_BandCount[band],
m_Sum[band], m_SqSum[band],
m_Min[band], m_Max[band]);
}
}
// Accessors
itkGetMacro(BandCount, PixelCountVectorType);
itkGetMacro(Sum, TRealVectorPixelType);
itkGetMacro(SqSum, TRealVectorPixelType);
itkGetMacro(Min, TRealVectorPixelType);
......@@ -118,11 +134,14 @@ public:
itkGetMacro(Count, double);
private:
void UpdateValues(const RealValueType & otherSum, const RealValueType & otherSqSum,
const RealValueType & otherMin, const RealValueType & otherMax,
void UpdateValues(PixelCountType otherCount,
RealValueType otherSum, RealValueType otherSqSum,
RealValueType otherMin, RealValueType otherMax,
PixelCountType & count,
RealValueType & sum, RealValueType & sqSum,
RealValueType & min, RealValueType & max)
{
count += otherCount;
sum += otherSum;
sqSum += otherSqSum;
if (otherMin < min)
......@@ -132,11 +151,14 @@ private:
}
protected:
PixelCountVectorType m_BandCount;
TRealVectorPixelType m_Sum;
TRealVectorPixelType m_SqSum;
TRealVectorPixelType m_Min;
TRealVectorPixelType m_Max;
double m_Count;
RealValueType m_NoDataValue;
PixelCountType m_Count;
bool m_UseNoDataValue;
};
/** \class PersistentStreamingStatisticsMapFromLabelImageFilter
......@@ -198,6 +220,11 @@ public:
itkStaticConstMacro(ImageDimension, unsigned int,
TInputVectorImage::ImageDimension);
itkGetMacro(NoDataValue, VectorPixelValueType);
itkSetMacro(NoDataValue, VectorPixelValueType);
itkGetMacro(UseNoDataValue, bool);
itkSetMacro(UseNoDataValue, bool);
/** Smart Pointer type to a DataObject. */
typedef typename itk::DataObject::Pointer DataObjectPointer;
typedef itk::ProcessObject::DataObjectPointerArraySizeType DataObjectPointerArraySizeType;
......@@ -260,6 +287,9 @@ private:
PersistentStreamingStatisticsMapFromLabelImageFilter(const Self &) = delete;
void operator =(const Self&) = delete;
VectorPixelValueType m_NoDataValue;
bool m_UseNoDataValue;
AccumulatorMapCollectionType m_AccumulatorMaps;
PixelValueMapType m_MeanRadiometricValue;
......@@ -337,6 +367,9 @@ public:
typedef TInputVectorImage VectorImageType;
typedef TLabelImage LabelImageType;
typedef typename VectorImageType::PixelType VectorPixelType;
typedef typename VectorImageType::PixelType::ValueType VectorPixelValueType;
typedef typename Superclass::FilterType::PixelValueMapType PixelValueMapType;
typedef typename Superclass::FilterType::PixelValueMapObjectType PixelValueMapObjectType;
......@@ -397,6 +430,30 @@ public:
return this->GetFilter()->GetLabelPopulationMap();
}
/** Set the no data value */
void SetNoDataValue(VectorPixelValueType value)
{
this->GetFilter()->SetNoDataValue(value);
}
/** Return the no data value */
VectorPixelValueType GetNoDataValue() const
{
return this->GetFilter()->GetNoDataValue();
}
/** Configure whether no data pixels ignored, treating each band independently */
void SetUseNoDataValue(bool useNoDataValue)
{
this->GetFilter()->SetUseNoDataValue(useNoDataValue);
}
/** Return whether no data pixels are ignored */
bool GetUseNoDataValue() const
{
return this->GetFilter()->GetUseNoDataValue();
}
protected:
/** Constructor */
StreamingStatisticsMapFromLabelImageFilter() {}
......
......@@ -35,6 +35,7 @@ namespace otb
template<class TInputVectorImage, class TLabelImage>
PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelImage>
::PersistentStreamingStatisticsMapFromLabelImageFilter()
: m_UseNoDataValue()
{
// first output is a copy of the image, DataObject created by
// superclass
......@@ -155,20 +156,21 @@ PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelIm
{
// Update temporary accumulator
AccumulatorMapType outputAcc;
auto endAcc = outputAcc.end();
for (auto const& threadAccMap: m_AccumulatorMaps)
{
for(auto const& it: threadAccMap)
{
const LabelPixelType label = it.first;
if (outputAcc.count(label) <= 0)
auto label = it.first;
auto itAcc = outputAcc.find(label);
if (itAcc == endAcc)
{
AccumulatorType newAcc(it.second);
outputAcc[label] = newAcc;
outputAcc.emplace(label, it.second);
}
else
{
outputAcc[label].Update(it.second);
itAcc->second.Update(it.second);
}
}
}
......@@ -177,18 +179,20 @@ PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelIm
for(auto& it: outputAcc)
{
const LabelPixelType label = it.first;
const double count = it.second.GetCount();
const RealVectorPixelType sum = it.second.GetSum();
const RealVectorPixelType sqSum = it.second.GetSqSum();
const auto &bandCount = it.second.GetBandCount();
const auto &sum = it.second.GetSum();
const auto &sqSum = it.second.GetSqSum();
// Count
m_LabelPopulation[label] = count;
m_LabelPopulation[label] = it.second.GetCount();
// Mean & stdev
RealVectorPixelType mean (sum);
RealVectorPixelType std (sqSum);
for (unsigned int band = 0 ; band < mean.GetSize() ; band++)
{
// Number of valid pixels in band
auto count = bandCount[band];
// Mean
mean[band] /= count;
......@@ -265,28 +269,31 @@ PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelIm
itk::ImageRegionConstIterator<TInputVectorImage> inIt(inputPtr, outputRegionForThread);
itk::ImageRegionConstIterator<TLabelImage> labelIt(labelInputPtr, outputRegionForThread);
itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
typename VectorImageType::PixelType value;
typename LabelImageType::PixelType label;
auto &acc = m_AccumulatorMaps[threadId];
auto endAcc = acc.end();
// do the work
for (inIt.GoToBegin(), labelIt.GoToBegin();
!inIt.IsAtEnd() && !labelIt.IsAtEnd();
++inIt, ++labelIt)
{
value = inIt.Get();
label = labelIt.Get();
const auto &value = inIt.Get();
auto label = labelIt.Get();
// Update the accumulator
if (m_AccumulatorMaps[threadId].count(label) <= 0) //add new element to the map
auto itAcc = acc.find(label);
if (itAcc == endAcc)
{
AccumulatorType newAcc(value);
m_AccumulatorMaps[threadId][label] = newAcc;
acc.emplace(label, AccumulatorType(this->GetNoDataValue(), this->GetUseNoDataValue(), value));
}
else
{
m_AccumulatorMaps[threadId][label].Update(value);
itAcc->second.Update(value);
}
progress.CompletedPixel();
}
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment