otbZonalStatistics.cxx 29.60 KiB
/*
 * Copyright (C) 2017 National Research Institute of Science and
 * Technology for Environment and Agriculture (IRSTEA)
 * This file is part of Orfeo Toolbox
 *     https://www.orfeo-toolbox.org/
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *     http://www.apache.org/licenses/LICENSE-2.0
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
#include "itkFixedArray.h"
#include "itkObjectFactory.h"
// Elevation handler
#include "otbWrapperElevationParametersHandler.h"
#include "otbWrapperApplicationFactory.h"
// Application engine
#include "otbStandardFilterWatcher.h"
// Process objects
#include "otbVectorDataToLabelImageFilter.h"
#include "otbVectorDataIntoImageProjectionFilter.h"
#include "otbStreamingStatisticsMapFromLabelImageFilter.h"
#include "otbStatisticsXMLFileWriter.h"
// Raster --> Vector
#include "otbLabelImageToVectorDataFilter.h"
#include "itkBinaryThresholdImageFilter.h"
#include "otbFunctorImageFilter.h"
namespace otb
namespace Wrapper
class ZonalStatistics : public Application
public:
  /** Standard class typedefs. */
  typedef ZonalStatistics               Self;
  typedef Application                   Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
  /* Typedefs */
  typedef Int32ImageType                LabelImageType;
  typedef LabelImageType::ValueType     LabelValueType;
  typedef otb::VectorData<double, 2>    VectorDataType;
  typedef otb::VectorDataIntoImageProjectionFilter<VectorDataType,
                                                   FloatVectorImageType>             VectorDataReprojFilterType;
  typedef otb::VectorDataToLabelImageFilter<VectorDataType,
                                            LabelImageType>                   RasterizeFilterType;
  typedef VectorDataType::DataTreeType  DataTreeType;
  typedef itk::PreOrderTreeIterator<DataTreeType>
  TreeIteratorType;
  typedef VectorDataType::DataNodeType  DataNodeType;
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
typedef DataNodeType::PolygonListPointerType PolygonListPointerType; typedef otb::StreamingStatisticsMapFromLabelImageFilter<FloatVectorImageType, LabelImageType> StatsFilterType; typedef otb::StatisticsXMLFileWriter<FloatVectorImageType::PixelType> StatsWriterType; typedef otb::LabelImageToVectorDataFilter<LabelImageType> LabelImageToVectorFilterType; typedef itk::BinaryThresholdImageFilter<LabelImageType, LabelImageType> ThresholdFilterType; template <class TInput, class TOutput> struct EncoderFunctorType { StatsFilterType::LabelPopulationMapType* m_CountMap; StatsFilterType::PixelValueMapType* m_MeanMap; 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; size_t m_NbStatsPerBand; size_t m_NbGlobalStats; std::vector<double> keys_list; size_t OutputSize(const std::array<size_t, 1>&) const { return m_NbInputComponents*m_NbStatsPerBand+m_NbGlobalStats; } size_t OutputSize() const { return m_NbInputComponents * m_NbStatsPerBand + m_NbGlobalStats; } TOutput operator()(TInput const& pix) { TOutput outPix(OutputSize()); outPix.Fill(m_OutBvValue); if(pix != m_InNoData && keys_list.size() == 0) { outPix[0] = (*m_CountMap)[pix]; for(size_t i=0; i<m_NbInputComponents; ++i) { outPix[i*m_NbStatsPerBand+1] = (*m_MeanMap)[pix][i]; outPix[i*m_NbStatsPerBand+2] = (*m_StdMap)[pix][i]; outPix[i*m_NbStatsPerBand+3] = (*m_MinMap)[pix][i]; 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; }
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
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|| m_MaxMap != other.m_MaxMap|| m_NbInputComponents != other.m_NbInputComponents || m_InNoData != other.m_InNoData || m_OutBvValue != other.m_OutBvValue ) return true; else return false; } }; typedef otb::FunctorImageFilter<EncoderFunctorType<LabelValueType, FloatVectorImageType::PixelType>> EncoderFilterType; /** Standard macro */ itkNewMacro(Self); itkTypeMacro(ZonalStatistics, Application); void DoInit() override { SetName("ZonalStatistics"); SetDescription("This application computes zonal statistics"); // Documentation SetDocLongDescription("This application computes zonal statistics from label image, or vector data. " "The application inputs one input multiband image, and another input for zones definition. " "Zones can be defined with a label image (inzone.labelimage.in) or a vector data layer " "(inzone.vector.in). The following statistics are computed over each zones: mean, min, max, " "and standard deviation. Statistics can be exported in a vector layer (if the input zone " "definition is a label image, it will be vectorized) or in a XML file"); SetDocLimitations("1) The inzone.vector.in must fit in memory (if \"inzone\" is \"vector\"). 2) The vectorized label " "image must also fit in memory (if \"out\" is \"vector\"): if not, consider using \"out\" to " "\"xml\""); SetDocAuthors("Remi Cresson, Jordi Inglada"); SetDocSeeAlso("ComputeImagesStatistics"); AddDocTag(Tags::Manip); AddDocTag(Tags::Analysis); // Input image AddParameter(ParameterType_InputImage, "in", "Input Image"); AddParameter(ParameterType_Float, "inbv", "Background value to ignore in statistics computation"); MandatoryOff ("inbv"); // Input zone mode AddParameter(ParameterType_Choice, "inzone", "Type of input for the zone definitions"); AddChoice("inzone.vector", "Input objects from vector data"); AddChoice("inzone.labelimage", "Input objects from label image"); // Input for vector mode AddParameter(ParameterType_InputVectorData, "inzone.vector.in", "Input vector data"); AddParameter(ParameterType_Bool, "inzone.vector.reproject", "Reproject the input vector"); // Input for label image mode AddParameter(ParameterType_InputImage, "inzone.labelimage.in", "Input label image"); AddParameter(ParameterType_Int, "inzone.labelimage.nodata", "No-data value for the input label image"); MandatoryOff ("inzone.labelimage.nodata"); // Output stats mode AddParameter(ParameterType_Choice, "out", "Format of the output stats"); AddChoice("out.vector", "Output vector data"); 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");
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
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(); // Doc example parameter settings SetDocExampleParameterValue("in", "input.tif"); SetDocExampleParameterValue("inzone.vector.in", "myvector.shp"); SetDocExampleParameterValue("out.vector.filename", "myvector_with_stats.shp"); SetOfficialDocLink(); } void DoUpdateParameters() override { // Nothing to do here : all parameters are independent } // Returns a string of the kind "prefix_i" const std::string CreateFieldName(const std::string & prefix, const unsigned int i) { std::stringstream ss; ss << prefix << "_" << i; return ss.str(); } // Returns a null pixel which has the same number of components per pixels as img FloatVectorImageType::PixelType NullPixel(FloatVectorImageType::Pointer & img) { const unsigned int nBands = img->GetNumberOfComponentsPerPixel(); FloatVectorImageType::PixelType pix(nBands); pix.Fill(0); return pix; } void GetStats() { m_CountMap = m_StatsFilter->GetLabelPopulationMap(); m_MeanMap = m_StatsFilter->GetMeanValueMap(); m_StdMap = m_StatsFilter->GetStandardDeviationValueMap(); m_MinMap = m_StatsFilter->GetMinValueMap(); m_MaxMap = m_StatsFilter->GetMaxValueMap(); m_HistoMap = m_StatsFilter->GetHistogrammeMap(); } void PrepareForLabelImageInput() { otbAppLogINFO("Zone definition: label image"); // Computing stats m_StatsFilter->SetInputLabelImage(GetParameterInt32Image("inzone.labelimage.in")); m_StatsFilter->Update(); // In this zone definition mode, the user can provide a no-data value for the labels if (HasUserValue("inzone.labelimage.nodata")) m_IntNoData = GetParameterInt("inzone.labelimage.nodata"); GetStats(); } void PrepareForVectorDataInput() { otbAppLogINFO("Zone definition: vector");
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
otbAppLogINFO("Loading vector data..."); m_VectorDataSrc = GetParameterVectorData("inzone.vector.in"); // Reproject vector data if (GetParameterInt("inzone.vector.reproject") != 0) { ReprojectVectorDataIntoInputImage(); } RasterizeInputVectorData(); // Computing stats m_StatsFilter->SetInputLabelImage(m_RasterizeFilter->GetOutput()); m_StatsFilter->Update(); GetStats(); } void ReprojectVectorDataIntoInputImage() { otbAppLogINFO("Vector data reprojection enabled"); m_VectorDataReprojectionFilter = VectorDataReprojFilterType::New(); m_VectorDataReprojectionFilter->SetInputVectorData(m_VectorDataSrc.GetPointer()); m_VectorDataReprojectionFilter->SetInputImage(m_InputImage); AddProcess(m_VectorDataReprojectionFilter, "Reproject vector data"); m_VectorDataReprojectionFilter->Update(); m_VectorDataSrc = m_VectorDataReprojectionFilter->GetOutput(); } void RasterizeInputVectorData() { // Rasterize vector data m_RasterizeFilter = RasterizeFilterType::New(); m_RasterizeFilter->AddVectorData(m_VectorDataSrc); m_RasterizeFilter->SetOutputOrigin(m_InputImage->GetOrigin()); m_RasterizeFilter->SetOutputSpacing(m_InputImage->GetSignedSpacing()); m_RasterizeFilter->SetOutputSize(m_InputImage->GetLargestPossibleRegion().GetSize()); m_RasterizeFilter->SetOutputProjectionRef(m_InputImage->GetProjectionRef()); m_RasterizeFilter->SetBurnAttribute("________"); m_RasterizeFilter->SetDefaultBurnValue(0); m_RasterizeFilter->SetGlobalWarningDisplay(false); m_RasterizeFilter->SetBackgroundValue(m_IntNoData); AddProcess(m_RasterizeFilter, "Rasterize input vector data"); } void RemoveNoDataEntry() { if (( GetParameterAsString("inzone") == "labelimage" && HasUserValue("inzone.labelimage.nodata")) || (GetParameterAsString("inzone") == "vector") ) { otbAppLogINFO("Removing entries for label value " << m_IntNoData); m_CountMap.erase(m_IntNoData); m_MeanMap.erase(m_IntNoData); m_StdMap.erase(m_IntNoData); m_MinMap.erase(m_IntNoData); m_MaxMap.erase(m_IntNoData); } } void GenerateVectorDataFromLabelImage() { // Mask for label image m_InputThresholdFilter = ThresholdFilterType::New(); m_InputThresholdFilter->SetInput(GetParameterInt32Image("inzone.labelimage.in")); m_InputThresholdFilter->SetInsideValue(0); m_InputThresholdFilter->SetOutsideValue(1); m_InputThresholdFilter->SetLowerThreshold(m_IntNoData); m_InputThresholdFilter->SetUpperThreshold(m_IntNoData); m_InputThresholdFilter->UpdateOutputInformation(); AddProcess(m_InputThresholdFilter, "Threshold label image"); // Vectorize the image m_LabelImageToVectorFilter = LabelImageToVectorFilterType::New();
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
m_LabelImageToVectorFilter->SetInput(GetParameterInt32Image("inzone.labelimage.in")); m_LabelImageToVectorFilter->SetInputMask(m_InputThresholdFilter->GetOutput()); m_LabelImageToVectorFilter->SetFieldName("polygon_id"); AddProcess(m_LabelImageToVectorFilter, "Vectorize label image"); m_LabelImageToVectorFilter->Update(); // The source vector data is the vectorized label image 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;
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
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 otbAppLogINFO("Writing output 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(); 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] );
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
for (unsigned int band = 0 ; band < m_InputImage->GetNumberOfComponentsPerPixel() ; band++) { currentGeometry->SetFieldAsDouble(CreateFieldName("mean", band), m_MeanMap[internalFID][band] ); currentGeometry->SetFieldAsDouble(CreateFieldName("stdev", band), m_StdMap [internalFID][band] ); currentGeometry->SetFieldAsDouble(CreateFieldName("min", band), m_MinMap [internalFID][band] ); currentGeometry->SetFieldAsDouble(CreateFieldName("max", band), m_MaxMap [internalFID][band] ); } m_NewVectorData->GetDataTree()->Add(currentGeometry, folder); } } ++itVector; } // next feature 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++; }
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
// 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")) { m_OutBvValue = GetParameterFloat("out.raster.bv"); } else if(HasUserValue("inbv")) { m_OutBvValue = GetParameterFloat("inbv"); } else { m_OutBvValue = m_IntNoData; } } void WriteRasterData() { otbAppLogINFO("Writing output raster data"); SetOutBvValue(); m_OutputThresholdFilter = ThresholdFilterType::New(); m_EncoderFilter = EncoderFilterType::New(); if(m_FromLabelImage) { m_EncoderFilter->SetInput(GetParameterInt32Image("inzone.labelimage.in")); } else { 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; 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 raster image will have " << (m_EncoderFilter->GetFunctor()).OutputSize() << " bands\n"); AddProcess(m_EncoderFilter, "Encode output raster image"); 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();
631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700
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(); } 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);
701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770
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 { // Get input image m_InputImage = GetParameterImage("in"); // Statistics filter m_StatsFilter = StatsFilterType::New(); m_StatsFilter->SetInput(m_InputImage); if (HasUserValue("inbv")) { m_StatsFilter->SetUseNoDataValue(true); m_StatsFilter->SetNoDataValue(GetParameterFloat("inbv")); } m_StatsFilter->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram")); AddProcess(m_StatsFilter->GetStreamer(), "Computing statistics"); // Select zone definition mode m_FromLabelImage = (GetParameterAsString("inzone") == "labelimage"); if (m_FromLabelImage) PrepareForLabelImageInput(); else if (GetParameterAsString("inzone") == "vector") PrepareForVectorDataInput(); else otbAppLogFATAL("Unknown zone definition mode"); // Remove the no-data entry RemoveNoDataEntry(); // Generate output if (GetParameterAsString("out") == "xml") { if(GetParameterInt("histogram") == 0) WriteXMLStatsFile(); else WriteXMLHistogramFile(); } else // vector or raster { if (m_FromLabelImage) GenerateVectorDataFromLabelImage(); if (GetParameterAsString("out") == "vector" && GetParameterInt("histogram") == 0 ) WriteVectorData(); else if (GetParameterAsString("out") == "vector" && GetParameterInt("histogram") != 0 ) WriteVectorHistogramData(); else if (GetParameterAsString("out") == "raster" && GetParameterInt("histogram") == 0) WriteRasterData(); else if (GetParameterAsString("out") == "raster" && GetParameterInt("histogram") != 0) WriteRasterHistogramData(); else otbAppLogFATAL("Unknown output mode"); } } VectorDataType::Pointer m_VectorDataSrc; VectorDataType::Pointer m_NewVectorData; VectorDataReprojFilterType::Pointer m_VectorDataReprojectionFilter; RasterizeFilterType::Pointer m_RasterizeFilter; StatsFilterType::Pointer m_StatsFilter; LabelImageToVectorFilterType::Pointer m_LabelImageToVectorFilter; ThresholdFilterType::Pointer m_InputThresholdFilter; ThresholdFilterType::Pointer m_OutputThresholdFilter; FloatVectorImageType::Pointer m_InputImage; LabelValueType m_IntNoData = itk::NumericTraits<LabelValueType>::max();
771772773774775776777778779780781782783784785786
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; StatsFilterType::PixelValueMapType m_MaxMap; EncoderFilterType::Pointer m_EncoderFilter; }; } } OTB_APPLICATION_EXPORT( otb::Wrapper::ZonalStatistics )