An error occurred while loading the file. Please try again.
-
Commandre Benjamin authored170c0d8d
/*========================================================================= Copyright (c) Remi Cresson (IRSTEA). All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information.=========================================================================*/#include "itkFixedArray.h"#include "itkObjectFactory.h"// Elevation handler#include "otbWrapperElevationParametersHandler.h"#include "otbWrapperApplicationFactory.h"// Application engine#include "otbStandardFilterWatcher.h"#include "itkFixedArray.h"// Filter#include "otbVectorDataToLabelImageFilter.h"#include "otbVectorDataIntoImageProjectionFilter.h"// Vnl vector#include "vnl/vnl_vector.h"// itk iterator#include "itkImageRegionConstIterator.h"#include <map>using namespace std;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; /* vector data filters */ typedef size_t LabelValueType; typedef otb::Image<LabelValueType, 2> LabelImageType; 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; typedef DataNodeType::PolygonListPointerType PolygonListPointerType; /** Statistics */ typedef vnl_matrix<double> RealMatrix; typedef vnl_vector<LabelValueType> SizeVector; /** Image iterator */7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
typedef itk::ImageRegionConstIterator<LabelImageType> LabelIteratorType;
typedef itk::ImageRegionConstIterator<FloatVectorImageType> ImageIteratorType;
/** Standard macro */
itkNewMacro(Self);
itkTypeMacro(ZonalStatistics, Application);
VectorDataType::Pointer shp;
VectorDataType::Pointer vdWithFid;
VectorDataType::Pointer vdWithStats;
VectorDataReprojFilterType::Pointer m_VectorDataReprojectionFilter;
RasterizeFilterType::Pointer m_RasterizeFilter;
void DoInit() {
SetName("ZonalStatistics");
SetDescription("Computes zonal statistics");
// Documentation
SetDocName("ZonalStatistics");
SetDocLongDescription("This application zonal statistics");
SetDocLimitations("None");
SetDocAuthors("Remi Cresson");
AddDocTag(Tags::Manip);
// Inputs
AddParameter(ParameterType_InputImage, "in", "Input Image");
AddParameter(ParameterType_InputVectorData, "vec", "Input vector data");
AddParameter(ParameterType_Int, "band", "Image band");
// AddParameter(ParameterType_Int, "ram", "ram");
AddParameter(ParameterType_StringList, "stats", "Computed statistics");
// Output
AddParameter(ParameterType_StringList, "out", "Output vector data");
AddRAMParameter();
}
void DoUpdateParameters(){
// Nothing to do here : all parameters are independent
}
void DoExecute(){
// Get inputs
FloatVectorImageType::Pointer img = GetParameterImage("in");
VectorDataType* shp = GetParameterVectorData("vec");
// Add a FID field
LabelValueType internalFID = 1;
const string internalFIDField = "__fid___";
vdWithFid = VectorDataType::New();
DataNodeType::Pointer root1 = vdWithFid->GetDataTree()->GetRoot()->Get();
DataNodeType::Pointer document1 = DataNodeType::New();
document1->SetNodeType(otb::DOCUMENT);
vdWithFid->GetDataTree()->Add(document1, root1);
DataNodeType::Pointer folder1 = DataNodeType::New();
folder1->SetNodeType(otb::FOLDER);
vdWithFid->GetDataTree()->Add(folder1, document1);
vdWithFid->SetProjectionRef(shp->GetProjectionRef());
TreeIteratorType itVector1(shp->GetDataTree());
itVector1.GoToBegin();
while (!itVector1.IsAtEnd()){
if (!itVector1.Get()->IsRoot() && !itVector1.Get()->IsDocument() && !itVector1.Get()->IsFolder()){
DataNodeType::Pointer currentGeometry = itVector1.Get();
currentGeometry->SetFieldAsInt(internalFIDField, internalFID );
vdWithFid->GetDataTree()->Add(currentGeometry, folder1);
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
internalFID++;
}
++itVector1;
} // next feature
// Rasterize vector data
m_RasterizeFilter = RasterizeFilterType::New();
m_RasterizeFilter->AddVectorData(vdWithFid);
m_RasterizeFilter->SetOutputOrigin(img->GetOrigin());
m_RasterizeFilter->SetOutputSpacing(img->GetSignedSpacing());
m_RasterizeFilter->SetOutputSize(img->GetLargestPossibleRegion().GetSize());
m_RasterizeFilter->SetOutputProjectionRef(img->GetProjectionRef());
m_RasterizeFilter->SetBurnAttribute(internalFIDField);
m_RasterizeFilter->UpdateOutputInformation();
// Computing stats
otbAppLogINFO("Computing stats");
// Explicit streaming over the input target image, based on the RAM parameter
typedef otb::RAMDrivenStrippedStreamingManager<FloatVectorImageType> StreamingManagerType;
StreamingManagerType::Pointer m_StreamingManager = StreamingManagerType::New();
m_StreamingManager->SetAvailableRAMInMB(GetParameterInt("ram"));
m_StreamingManager->PrepareStreaming(img, img->GetLargestPossibleRegion() );
// Init. stats containers
const LabelValueType N = internalFID;
const unsigned int nBands = 1;
RealMatrix sum (nBands, N, 0);
RealMatrix means (nBands, N, 0);
SizeVector counts(N,0);
if (m_RasterizeFilter->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels()
!= img->GetLargestPossibleRegion().GetNumberOfPixels() ){
otbAppLogFATAL("Rasterized image and input image don't have the same number of pixels");
}
int m_NumberOfDivisions = m_StreamingManager->GetNumberOfSplits();
std::map<int,std::map<int,int>> dico;
const unsigned int band = GetParameterInt("band");
for (int m_CurrentDivision = 0; m_CurrentDivision < m_NumberOfDivisions; m_CurrentDivision++){
FloatVectorImageType::RegionType streamRegion = m_StreamingManager->GetSplit(m_CurrentDivision);
img->SetRequestedRegion(streamRegion);
img->PropagateRequestedRegion();
img->UpdateOutputData();
m_RasterizeFilter->GetOutput()->SetRequestedRegion(streamRegion);
m_RasterizeFilter->GetOutput()->PropagateRequestedRegion();
m_RasterizeFilter->GetOutput()->UpdateOutputData();
LabelIteratorType it(m_RasterizeFilter->GetOutput(), streamRegion);
ImageIteratorType it_in(img, streamRegion);
for (it.GoToBegin(), it_in.GoToBegin(); !it.IsAtEnd(); ++it, ++it_in){
const LabelValueType fid = it.Get();
const FloatVectorImageType::PixelType pix = it_in.Get();
const FloatVectorImageType::InternalPixelType val = pix[band];
if ( dico[fid].find(val) == dico[fid].end() ) {
dico[fid][val] = 1;
} else {
dico[fid][val]++;
}
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279
sum[0][fid]+=val;
counts[fid]++;
} // next pixel
}
// Summarize stats
otbAppLogINFO("Summarizing stats");
for (LabelValueType fid = 0 ; fid < N ; fid++){
const LabelValueType count = counts[fid];
if (count>0){
means[0][fid] = sum[0][fid] / static_cast<FloatVectorImageType::InternalPixelType>(count);
}
}
std::vector<std::string> statistiques = GetParameterStringList("stats");
std::vector<std::string> resultats;
std::ostringstream ss;
bool compute_mean = false;
for (unsigned int i = 0; i < statistiques.size(); i++){
if (statistiques.at(i) == ("mean")){
compute_mean = true;
}
}
resultats.reserve(N-1);
if(compute_mean){
for (LabelValueType fid = 1 ; fid < N ; fid++){
ss << means[0][fid];
resultats.push_back(ss.str());
ss.str("");
ss.clear();
}
}else{
for (LabelValueType fid = 1 ; fid < N ; fid++){
float Maj_count = 0.0;
float Maj_count_perc = 0.0;
for(auto t : dico[fid]){
if((float(t.second)/float(counts[fid]))*100.0 > Maj_count_perc ){
Maj_count = t.first;
Maj_count_perc = (float(t.second)/float(counts[fid]))*100.0;
}
}
ss << counts[fid] << " " << Maj_count << " " << Maj_count_perc;
resultats.push_back(ss.str());
ss.str("");
ss.clear();
}
}
SetParameterStringList("out", resultats);
}
};
}
}
OTB_APPLICATION_EXPORT(otb::Wrapper::ZonalStatistics)