From 09b0920034e747609a0b9deff081791b0f7a095b Mon Sep 17 00:00:00 2001 From: Commandre Benjamin <benjamin.commandre@irstea.fr> Date: Thu, 23 May 2019 13:13:42 +0200 Subject: [PATCH] Push zonalstats' modifications for histogram --- .../app/otbZonalStatistics.cxx | 318 +++++++++++++++++- ...reamingStatisticsMapFromLabelImageFilter.h | 42 ++- ...amingStatisticsMapFromLabelImageFilter.hxx | 12 + 3 files changed, 351 insertions(+), 21 deletions(-) diff --git a/Modules/Applications/AppClassification/app/otbZonalStatistics.cxx b/Modules/Applications/AppClassification/app/otbZonalStatistics.cxx index 881e8dd97a..d7e614b3f0 100644 --- a/Modules/Applications/AppClassification/app/otbZonalStatistics.cxx +++ b/Modules/Applications/AppClassification/app/otbZonalStatistics.cxx @@ -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; diff --git a/Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.h b/Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.h index 81372106be..4f5dec7a5c 100644 --- a/Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.h +++ b/Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.h @@ -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) { diff --git a/Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.hxx b/Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.hxx index 7a0ecf30e8..69a439fc2d 100644 --- a/Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.hxx +++ b/Modules/Filtering/Statistics/include/otbStreamingStatisticsMapFromLabelImageFilter.hxx @@ -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()); } -- GitLab