diff --git "a/SimpleExtractionTools/app/otbZonalStatistics_r\303\251mi.cxx" "b/SimpleExtractionTools/app/otbZonalStatistics_r\303\251mi.cxx" deleted file mode 100644 index 0dfb49cfcb480002e3d045ecc264f78567baf8a6..0000000000000000000000000000000000000000 --- "a/SimpleExtractionTools/app/otbZonalStatistics_r\303\251mi.cxx" +++ /dev/null @@ -1,279 +0,0 @@ -/*========================================================================= - - 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 */ - 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); - 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]++; - } - - 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) \ No newline at end of file diff --git "a/SimpleExtractionTools/app/otbZonalStatistics_t\303\251moin.cxx" "b/SimpleExtractionTools/app/otbZonalStatistics_t\303\251moin.cxx" deleted file mode 100644 index b72e48c331f35efd943eb1850e50920f6ad5ec01..0000000000000000000000000000000000000000 --- "a/SimpleExtractionTools/app/otbZonalStatistics_t\303\251moin.cxx" +++ /dev/null @@ -1,358 +0,0 @@ -/*========================================================================= - - 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::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 */ - 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; - 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_InputImageList, "il", "Input List Images"); - AddParameter(ParameterType_InputVectorDataList, "vec", "Input List vector data"); - 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 ComputeMean(FloatVectorImageListType::Pointer listes_image, VectorDataListType::Pointer listes_vecteur){ - - } - - void DoExecute(){ - - // Get inputs - FloatVectorImageListType::Pointer listes_image = GetParameterImageList("il"); - VectorDataListType::Pointer listes_vecteur = GetParameterVectorDataList("vec"); - - ComputeMean(listes_image, listes_vecteur); - - std::map<int, int> liste_poly; - - /* - Vecteur : vec_i - Image : img_i - Bande : band_b - poly_p - */ - - std::map<int,std::map<int,RealMatrix>> dico_sum; - std::map<int,std::map<int,RealMatrix>> dico_count; - - std::map<int, LabelValueType> liste_taille; - - for (unsigned int vec = 0; vec < listes_vecteur->Size(); vec++){ - VectorDataType::Pointer shp = listes_vecteur->GetNthElement(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); - internalFID++; - } - ++itVector1; - } // next feature - - liste_poly[vec] = (internalFID - 1); - liste_taille[vec] = internalFID; - - for (unsigned int img = 0; img < listes_image->Size(); img++){ - FloatVectorImageType::Pointer image = listes_image->GetNthElement(img); - - // Rasterize vector data - m_RasterizeFilter = RasterizeFilterType::New(); - - m_RasterizeFilter->AddVectorData(vdWithFid); - m_RasterizeFilter->SetOutputOrigin(image->GetOrigin()); - m_RasterizeFilter->SetOutputSpacing(image->GetSignedSpacing()); - m_RasterizeFilter->SetOutputSize(image->GetLargestPossibleRegion().GetSize()); - m_RasterizeFilter->SetOutputProjectionRef(image->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(image, image->GetLargestPossibleRegion()); - - // Init. stats containers - const LabelValueType N = internalFID; - const unsigned int nBands = image->GetNumberOfComponentsPerPixel(); - - RealMatrix sum (nBands, N, 0); - RealMatrix counts (nBands, N, 0); - - if (m_RasterizeFilter->GetOutput()->GetLargestPossibleRegion().GetNumberOfPixels() - != image->GetLargestPossibleRegion().GetNumberOfPixels() ){ - otbAppLogFATAL("Rasterized image and input image don't have the same number of pixels"); - } - - int m_NumberOfDivisions = m_StreamingManager->GetNumberOfSplits(); - - for (int m_CurrentDivision = 0; m_CurrentDivision < m_NumberOfDivisions; m_CurrentDivision++){ - FloatVectorImageType::RegionType streamRegion = m_StreamingManager->GetSplit(m_CurrentDivision); - - image->SetRequestedRegion(streamRegion); - image->PropagateRequestedRegion(); - image->UpdateOutputData(); - - m_RasterizeFilter->GetOutput()->SetRequestedRegion(streamRegion); - m_RasterizeFilter->GetOutput()->PropagateRequestedRegion(); - m_RasterizeFilter->GetOutput()->UpdateOutputData(); - - LabelIteratorType it(m_RasterizeFilter->GetOutput(), streamRegion); - ImageIteratorType it_in(image, streamRegion); - - for (it.GoToBegin(), it_in.GoToBegin(); !it.IsAtEnd(); ++it, ++it_in){ - const LabelValueType fid = it.Get(); - const FloatVectorImageType::PixelType pix = it_in.Get(); - - for (unsigned int band = 0 ; band < nBands ; band++){ - - const FloatVectorImageType::InternalPixelType val = pix[band]; - - // if ( dico[fid].find(val) == dico[fid].end() ) { - // dico[fid][val] = 1; - // } else { - // dico[fid][val]++; - // } - - sum[band][fid] += val; - counts[band][fid]++; - } - } // next pixel - } - - dico_sum[vec][img] = sum; - dico_count[vec][img] = counts; - } - } - - // Summarize stats - otbAppLogINFO("Summarizing stats"); - - std::map<int,std::map<int,RealMatrix>> dico_mean; - - for (unsigned int vec = 0; vec < listes_vecteur->Size(); vec++){ - - LabelValueType N = liste_taille[vec]; - for (unsigned int img = 0; img < listes_image->Size(); img++){ - RealMatrix means (listes_image->GetNthElement(img)->GetNumberOfComponentsPerPixel(), N, 0); - - for (unsigned int band = 0 ; band < listes_image->GetNthElement(img)->GetNumberOfComponentsPerPixel() ; band++){ - - for (LabelValueType fid = 1 ; fid < N ; fid++){ - const LabelValueType count = dico_count[vec][img][band][fid]; - if (count>0){ - means[band][fid] = dico_sum[vec][img][band][fid] / static_cast<FloatVectorImageType::InternalPixelType>(count); - } - } - } - dico_mean[vec][img] = means; - } - } - - 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; - } - } - - unsigned int taille = 0; - - for (unsigned int t = 0; t < listes_image->Size(); t++){ - taille += listes_image->GetNthElement(t)->GetNumberOfComponentsPerPixel(); - } - - int nombre_poly = 0; - - for (unsigned int p = 0; p < liste_poly.size(); p++){ - nombre_poly += liste_poly[p]; - } - - // std::cout << nombre_poly << std::endl; - // std::cout << taille << std::endl; - // std::cout << nombre_poly*taille << std::endl; - - resultats.reserve(taille); - - std::string separator = ";" ; - - if(compute_mean){ - std::cout << "Compute mean" << std::endl; - - for (unsigned int vec = 0; vec < dico_mean.size(); vec++){ - - LabelValueType N = liste_taille[vec]; - for (unsigned int img = 0; img < dico_mean[vec].size(); img++){ - for (unsigned int band = 0; band < listes_image->GetNthElement(img)->GetNumberOfComponentsPerPixel(); band++){ - for (unsigned int poly = 1; poly < N ; poly++){ - // std::cout << dico_mean[img][vec][band][poly] << " "; - ss << dico_mean[vec][img][band][poly] << " "; - } - ss << separator ; - } - // std::cout << ";" << std::endl; - resultats.push_back(ss.str()); - ss.str(""); - ss.clear(); - } - } - // 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 percent = 0.0; - // 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 > percent ){ - // Maj_count = t.first; - // Maj_count_perc = (float(t.second)/float(counts[fid]))*100.0; - // percent = Maj_count_perc; - // } - // } - - // 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) \ No newline at end of file