Skip to content

GitLab

Commit 09b09200 authored by Commandre Benjamin's avatar Commandre Benjamin
Browse files

Push zonalstats' modifications for histogram

parent df935279
Pipeline #6172 canceled with stage
Hide whitespace changes
Inline Side-by-side
Modules/Applications/AppClassification/app/otbZonalStatistics.cxx
View file @ 09b09200
......@@ -87,11 +87,14 @@ public:
StatsFilterType::PixelValueMapType* m_StdMap;
StatsFilterType::PixelValueMapType* m_MinMap;
StatsFilterType::PixelValueMapType* m_MaxMap;
StatsFilterType::HistogrammeMap* m_HistoMap;
size_t m_NbInputComponents;
LabelValueType m_InNoData;
LabelValueType m_OutBvValue;
static constexpr size_t m_NbStatsPerBand{4};
static constexpr size_t m_NbGlobalStats{1};
size_t m_NbStatsPerBand;
size_t m_NbGlobalStats;
std::vector<double> keys_list;
size_t OutputSize(const std::array<size_t, 1>&) const
{
......@@ -107,7 +110,7 @@ public:
{
TOutput outPix(OutputSize());
outPix.Fill(m_OutBvValue);
if(pix != m_InNoData)
if(pix != m_InNoData && keys_list.size() == 0)
{
outPix[0] = (*m_CountMap)[pix];
for(size_t i=0; i<m_NbInputComponents; ++i)
......@@ -118,12 +121,27 @@ public:
outPix[i*m_NbStatsPerBand+4] = (*m_MaxMap)[pix][i];
}
}
else if(pix != m_InNoData && keys_list.size() != 0)
{
outPix[0] = (*m_CountMap)[pix];
for(size_t i=0; i<m_NbInputComponents; i++)
{
auto result = std::max_element((*m_HistoMap)[i][pix].begin(), (*m_HistoMap)[i][pix].end(), ValueMax());
outPix[1] = std::get<0>(*result);
for(size_t j = 0; j < keys_list.size(); j++)
{
auto pixval = keys_list[j];
outPix[i*m_NbStatsPerBand+(j+m_NbGlobalStats)] = (*m_HistoMap)[i][pix][pixval];
}
}
}
return outPix;
}
bool operator != (const EncoderFunctorType& other)
{
if ( m_CountMap != other.m_CountMap||
m_HistoMap != other.m_HistoMap ||
m_MeanMap != other.m_MeanMap||
m_StdMap != other.m_StdMap||
m_MinMap != other.m_MinMap||
......@@ -190,12 +208,20 @@ public:
AddParameter(ParameterType_OutputVectorData, "out.vector.filename", "Filename for the output vector data");
AddChoice("out.xml", "Output XML file");
AddParameter(ParameterType_String, "out.xml.filename", "Filename for the output xml file");
AddChoice("out.raster", "Output raster image");
AddParameter(ParameterType_OutputImage, "out.raster.filename", "File name for the raster image");
AddParameter(ParameterType_Float, "out.raster.bv", "Background value for the output raster");
MandatoryOff ("out.raster.bv");
AddChoice("out.raster", "Output raster image");
AddParameter(ParameterType_OutputImage, "out.raster.filename", "File name for the raster image");
AddParameter(ParameterType_Float, "out.raster.bv", "Background value for the output raster");
MandatoryOff ("out.raster.bv");
AddParameter(ParameterType_Bool, "histogram", "Compute histogram for each polygon.");
MandatoryOff ("histogram");
AddParameter(ParameterType_Bool, "percent", "Display histogram as percentage.");
MandatoryOff ("percent");
AddParameter(ParameterType_String, "nomenclature", "Nomenclature text file. Format = PixelValue : Label");
MandatoryOff ("nomenclature");
AddRAMParameter();
AddRAMParameter();
// Doc example parameter settings
SetDocExampleParameterValue("in", "input.tif");
......@@ -234,6 +260,7 @@ public:
m_StdMap = m_StatsFilter->GetStandardDeviationValueMap();
m_MinMap = m_StatsFilter->GetMinValueMap();
m_MaxMap = m_StatsFilter->GetMaxValueMap();
m_HistoMap = m_StatsFilter->GetHistogrammeMap();
}
void PrepareForLabelImageInput()
......@@ -330,6 +357,105 @@ public:
m_VectorDataSrc = m_LabelImageToVectorFilter->GetOutput();
}
struct ValueMax
{
template <typename Lhs, typename Rhs>
bool operator()(const Lhs& lhs, const Rhs& rhs) const
{
return lhs.second < rhs.second;
}
};
std::string GetFieldName(double pixValue, std::unordered_map<double, std::string> nomenclature)
{
auto dom = nomenclature.find(pixValue);
if (!GetParameterString("nomenclature").empty() && !(dom == nomenclature.end()))
{
return nomenclature[pixValue];
}
else
{
std::stringstream FieldName;
FieldName << std::fixed << std::setprecision(3) << pixValue;
return FieldName.str();
}
}
std::unordered_map<double, std::string> LoadNomenclature()
{
/**
Nomenclature file format :
PixelValue1 : Label1
PixelValue2 : Label2
...
**/
std::unordered_map<double, std::string> nomenclature;
std::ifstream inFile;
inFile.open(GetParameterString("nomenclature"));
if(!inFile)
{
std::cerr << "Unable to open nomenclature file.";
exit(1);
}
std::string line;
while (std::getline(inFile, line))
{
std::istringstream iss(line);
std::vector<std::string> results((std::istream_iterator<std::string>(iss)),
std::istream_iterator<std::string>());
nomenclature[std::stod(results[0])] = results[2];
}
inFile.close();
return nomenclature;
}
void WriteFeatureHistogram(DataNodeType::Pointer& geometry, const std::vector<double> keys_list, const LabelValueType FID)
{
std::stringstream FieldLabel;
std::unordered_map<double, std::string> nomenclature;
if(!GetParameterString("nomenclature").empty())
{
nomenclature = LoadNomenclature();
}
for (unsigned int band = 0 ; band < m_InputImage->GetNumberOfComponentsPerPixel() ; band++)
{
auto result = std::max_element(m_HistoMap[band][FID].begin(), m_HistoMap[band][FID].end(), ValueMax());
geometry->SetFieldAsString("dominance", GetFieldName(std::get<0>(*result), nomenclature));
for(unsigned int i = 0; i < keys_list.size(); i++)
{
double pixval = keys_list[i];
FieldLabel << GetFieldName(pixval, nomenclature) << "_" << band;
auto got = m_HistoMap[band][FID].find(pixval);
if(got == m_HistoMap[band][FID].end())
{
geometry->SetFieldAsDouble(FieldLabel.str(), 0.0);
}
else
{
if(GetParameterInt("percent") == 0)
{
geometry->SetFieldAsDouble(FieldLabel.str(), m_HistoMap[band][FID][pixval]);
}
else
{
geometry->SetFieldAsDouble(FieldLabel.str(), m_HistoMap[band][FID][pixval]/((double) m_CountMap[FID]));
}
}
std::stringstream().swap(FieldLabel);
}
}
}
void WriteVectorData()
{
// Add statistics fields
......@@ -378,6 +504,74 @@ public:
SetParameterOutputVectorData("out.vector.filename", m_NewVectorData);
}
std::vector<double> GetKeysList()
{
std::vector<double> keys_list;
for(const auto& iteratorBand : m_HistoMap)
{
for(const auto& iteratorFeature : iteratorBand.second)
{
for(const auto& iteratorPixel : iteratorFeature.second)
{
if(std::find(keys_list.begin(), keys_list.end(), iteratorPixel.first) == keys_list.end())
{
keys_list.push_back(iteratorPixel.first);
}
}
}
}
std::sort(keys_list.begin(), keys_list.end());
return keys_list;
}
void WriteVectorHistogramData()
{
// Add statistics fields
otbAppLogINFO("Writing output histogram vector data");
LabelValueType internalFID = -1;
m_NewVectorData = VectorDataType::New();
DataNodeType::Pointer root = m_NewVectorData->GetDataTree()->GetRoot()->Get();
DataNodeType::Pointer document = DataNodeType::New();
document->SetNodeType(otb::DOCUMENT);
m_NewVectorData->GetDataTree()->Add(document, root);
DataNodeType::Pointer folder = DataNodeType::New();
folder->SetNodeType(otb::FOLDER);
m_NewVectorData->GetDataTree()->Add(folder, document);
m_NewVectorData->SetProjectionRef(m_VectorDataSrc->GetProjectionRef());
TreeIteratorType itVector(m_VectorDataSrc->GetDataTree());
itVector.GoToBegin();
std::vector<double> keys_list = GetKeysList();
while(!itVector.IsAtEnd())
{
if (!itVector.Get()->IsRoot() && !itVector.Get()->IsDocument() && !itVector.Get()->IsFolder())
{
DataNodeType::Pointer currentGeometry = itVector.Get();
if (m_FromLabelImage)
{
internalFID = currentGeometry->GetFieldAsInt("polygon_id");
}
else
{
internalFID++;
}
// Add the geometry with the new fields
if (m_CountMap.count(internalFID) > 0)
{
currentGeometry->SetFieldAsDouble("count", m_CountMap[internalFID]);
WriteFeatureHistogram(currentGeometry, keys_list, internalFID);
m_NewVectorData->GetDataTree()->Add(currentGeometry, folder);
}
}
++itVector;
} // next feature
SetParameterOutputVectorData("out.vector.filename", m_NewVectorData);
}
void SetOutBvValue()
{
if (HasUserValue("out.raster.bv"))
......@@ -396,7 +590,6 @@ public:
void WriteRasterData()
{
otbAppLogINFO("Writing output raster data");
SetOutBvValue();
m_OutputThresholdFilter = ThresholdFilterType::New();
......@@ -410,6 +603,8 @@ public:
m_EncoderFilter->SetInput(m_RasterizeFilter->GetOutput());
}
m_EncoderFilter->GetModifiableFunctor().m_NbStatsPerBand = 4;
m_EncoderFilter->GetModifiableFunctor().m_NbGlobalStats = 1;
m_EncoderFilter->GetModifiableFunctor().m_CountMap = &m_CountMap;
m_EncoderFilter->GetModifiableFunctor().m_MeanMap = &m_MeanMap;
m_EncoderFilter->GetModifiableFunctor().m_StdMap = &m_StdMap;
......@@ -424,20 +619,106 @@ public:
SetParameterOutputImage("out.raster.filename", m_EncoderFilter->GetOutput());
}
void WriteRasterHistogramData()
{
otbAppLogINFO("Writing output histogram raster data");
SetOutBvValue();
m_OutputThresholdFilter = ThresholdFilterType::New();
m_EncoderFilter = EncoderFilterType::New();
std::vector<double> keys_list = GetKeysList();
if(m_FromLabelImage)
{
m_EncoderFilter->SetInput(GetParameterInt32Image("inzone.labelimage.in"));
}
else
{
m_EncoderFilter->SetInput(m_RasterizeFilter->GetOutput());
}
m_EncoderFilter->GetModifiableFunctor().m_NbGlobalStats = 2;
m_EncoderFilter->GetModifiableFunctor().m_NbStatsPerBand = keys_list.size();
m_EncoderFilter->GetModifiableFunctor().keys_list = keys_list;
m_EncoderFilter->GetModifiableFunctor().m_CountMap = &m_CountMap;
m_EncoderFilter->GetModifiableFunctor().m_HistoMap = &m_HistoMap;
m_EncoderFilter->GetModifiableFunctor().m_MeanMap = &m_MeanMap;
m_EncoderFilter->GetModifiableFunctor().m_StdMap = &m_StdMap;
m_EncoderFilter->GetModifiableFunctor().m_MinMap = &m_MinMap;
m_EncoderFilter->GetModifiableFunctor().m_MaxMap = &m_MaxMap;
m_EncoderFilter->GetModifiableFunctor().m_NbInputComponents = m_InputImage->GetNumberOfComponentsPerPixel();
m_EncoderFilter->GetModifiableFunctor().m_InNoData = m_IntNoData;
m_EncoderFilter->GetModifiableFunctor().m_OutBvValue = m_OutBvValue;
otbAppLogINFO("Output histogram raster image will have " << (m_EncoderFilter->GetFunctor()).OutputSize() << " bands\n");
AddProcess(m_EncoderFilter, "Encode output histogram raster image");
SetParameterOutputImage("out.raster.filename", m_EncoderFilter->GetOutput());
}
void WriteXMLStatsFile()
{
// Write statistics in XML file
const std::string outXMLFile = this->GetParameterString("out.xml.filename");
otbAppLogINFO("Writing " + outXMLFile)
StatsWriterType::Pointer statWriter = StatsWriterType::New();
statWriter->SetFileName(outXMLFile);
statWriter->AddInputMap<StatsFilterType::LabelPopulationMapType>("count",m_CountMap);
statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("mean",m_MeanMap);
statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("std",m_StdMap);
statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("min",m_MinMap);
statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("max",m_MaxMap);
statWriter->Update();
}
StatsWriterType::Pointer statWriter = StatsWriterType::New();
statWriter->SetFileName(outXMLFile);
statWriter->AddInputMap<StatsFilterType::LabelPopulationMapType>("count",m_CountMap);
statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("mean",m_MeanMap);
statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("std",m_StdMap);
statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("min",m_MinMap);
statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("max",m_MaxMap);
statWriter->Update();
}
void WriteXMLHistogramFile()
{
// Write statistics in XML file
const std::string outXMLFile = this->GetParameterString("out.xml.filename");
otbAppLogINFO("Writing " + outXMLFile)
StatsWriterType::Pointer statWriter = StatsWriterType::New();
statWriter->SetFileName(outXMLFile);
std::unordered_map<double, std::string> nomenclature;
if (!GetParameterString("nomenclature").empty())
{
nomenclature = LoadNomenclature();
}
std::stringstream FieldLabel;
for (size_t band = 0 ; band < m_InputImage->GetNumberOfComponentsPerPixel() ; band++)
{
std::unordered_map<double, std::string> dominance;
for(auto iterFeature : m_HistoMap[band])
{
if(iterFeature.first != m_IntNoData)
{
auto result = std::max_element(iterFeature.second.begin(), iterFeature.second.end(), ValueMax());
FieldLabel << GetFieldName(std::get<0>(*result), nomenclature);
dominance[iterFeature.first] = FieldLabel.str();
std::stringstream().swap(FieldLabel);
std::unordered_map<std::string, double> histogram;
for(auto iterValue : iterFeature.second)
{
if(GetParameterInt("percent") == 0)
{
histogram[GetFieldName(iterValue.first, nomenclature)] = iterValue.second;
}
else
{
histogram[GetFieldName(iterValue.first, nomenclature)] = iterValue.second/((double) m_CountMap[iterFeature.first]);
}
}
statWriter->AddInputMap<>(std::to_string(iterFeature.first).c_str(), histogram);
}
}
statWriter->AddInputMap<>("dominance", dominance);
}
statWriter->Update();
}
void DoExecute() override
{
......@@ -484,6 +765,7 @@ public:
LabelValueType m_OutBvValue = itk::NumericTraits<LabelValueType>::max();
bool m_FromLabelImage = false;
StatsFilterType::LabelPopulationMapType m_CountMap;
StatsFilterType::HistogrammeMap m_HistoMap;
StatsFilterType::PixelValueMapType m_MeanMap;
StatsFilterType::PixelValueMapType m_StdMap;
StatsFilterType::PixelValueMapType m_MinMap;
......
Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.h
View file @ 09b09200
......@@ -57,6 +57,7 @@ public:
typedef typename TRealVectorPixelType::ValueType RealValueType;
typedef uint64_t PixelCountType;
typedef itk::VariableLengthVector<PixelCountType> PixelCountVectorType;
typedef std::unordered_map<double, double> Histogramme;
// Constructor (default)
StatisticsAccumulator() : m_Count(), m_NoDataValue(), m_UseNoDataValue() {}
......@@ -82,6 +83,7 @@ public:
{
m_BandCount[band] = 1;
m_SqSum[band] *= m_SqSum[band];
m_Histo.insert({val, 1});
}
else
{
......@@ -99,13 +101,15 @@ public:
{
const RealValueType value = pixel[band];
const RealValueType sqValue = value * value;
const std::unordered_map<double, double> histo = {{value,1}};
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]);
m_Min[band], m_Max[band],
histo, m_Histo);
}
}
......@@ -121,7 +125,8 @@ public:
other.m_Min[band], other.m_Max[band],
m_BandCount[band],
m_Sum[band], m_SqSum[band],
m_Min[band], m_Max[band]);
m_Min[band], m_Max[band],
other.m_Histo, m_Histo);
}
}
......@@ -132,6 +137,7 @@ public:
itkGetMacro(Min, TRealVectorPixelType);
itkGetMacro(Max, TRealVectorPixelType);
itkGetMacro(Count, double);
itkGetMacro(Histo, Histogramme);
private:
void UpdateValues(PixelCountType otherCount,
......@@ -139,7 +145,8 @@ private:
RealValueType otherMin, RealValueType otherMax,
PixelCountType & count,
RealValueType & sum, RealValueType & sqSum,
RealValueType & min, RealValueType & max)
RealValueType & min, RealValueType & max,
std::unordered_map<double,double> otherHisto, std::unordered_map<double,double> & histo)
{
count += otherCount;
sum += otherSum;
......@@ -148,6 +155,20 @@ private:
min = otherMin;
if (otherMax > max)
max = otherMax;
for(const auto& it: otherHisto)
{
auto got = histo.find(it.first);
if(got == histo.end())
{
histo.insert(it);
}
else
{
histo[it.first] += it.second;
}
}
}
protected:
......@@ -159,6 +180,7 @@ protected:
RealValueType m_NoDataValue;
PixelCountType m_Count;
bool m_UseNoDataValue;
Histogramme m_Histo;
};
/** \class PersistentStreamingStatisticsMapFromLabelImageFilter
......@@ -212,6 +234,9 @@ public:
typedef std::vector<AccumulatorMapType> AccumulatorMapCollectionType;
typedef std::unordered_map<LabelPixelType, RealVectorPixelType > PixelValueMapType;
typedef std::unordered_map<LabelPixelType, double> LabelPopulationMapType;
typedef std::unordered_map<double,
std::unordered_map<LabelPixelType,
std::unordered_map<double, double>>> HistogrammeMap;
itkStaticConstMacro(InputImageDimension, unsigned int,
TInputVectorImage::ImageDimension);
......@@ -256,6 +281,9 @@ public:
/** Return the computed number of labeled pixels for each label in the input label image */
LabelPopulationMapType GetLabelPopulationMap() const;
/** Return the computed histogramme */
HistogrammeMap GetHistogrammeMap() const;
/** Make a DataObject of the correct type to be used as the specified
* output. */
DataObjectPointer MakeOutput(DataObjectPointerArraySizeType idx) override;
......@@ -298,6 +326,7 @@ private:
PixelValueMapType m_MaxRadiometricValue;
LabelPopulationMapType m_LabelPopulation;
HistogrammeMap m_Histogramme;
}; // end of class PersistentStreamingStatisticsMapFromLabelImageFilter
......@@ -374,6 +403,7 @@ public:
typedef typename Superclass::FilterType::PixelValueMapObjectType PixelValueMapObjectType;
typedef typename Superclass::FilterType::LabelPopulationMapType LabelPopulationMapType;
typedef typename Superclass::FilterType::HistogrammeMap HistogrammeMap;
/** Set input multispectral image */
using Superclass::SetInput;
......@@ -430,6 +460,12 @@ public:
return this->GetFilter()->GetLabelPopulationMap();
}
/** Return the computed number of labeled pixels for each label */
HistogrammeMap GetHistogrammeMap() const
{
return this->GetFilter()->GetHistogrammeMap();
}
/** Set the no data value */
void SetNoDataValue(VectorPixelValueType value)
{
......
Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.hxx
View file @ 09b09200
......@@ -118,6 +118,14 @@ PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelIm
return m_LabelPopulation;
}
template<class TInputVectorImage, class TLabelImage>
typename PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelImage>::HistogrammeMap
PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelImage>
::GetHistogrammeMap() const
{
return m_Histogramme;
}
template<class TInputVectorImage, class TLabelImage>
void
PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelImage>
......@@ -196,6 +204,9 @@ PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelIm
// Mean
mean[band] /= count;
// Histogramme
m_Histogramme[band][label] = it.second.GetHisto();
// Unbiased standard deviation (not sure unbiased is usefull here)
const double variance = (sqSum[band] - (sum[band] * mean[band])) / (count - 1);
std[band] = std::sqrt(variance);
......@@ -222,6 +233,7 @@ PersistentStreamingStatisticsMapFromLabelImageFilter<TInputVectorImage, TLabelIm
m_MinRadiometricValue.clear();
m_MaxRadiometricValue.clear();
m_LabelPopulation.clear();
m_Histogramme.clear();
m_AccumulatorMaps.resize(this->GetNumberOfThreads());
}
......
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