Skip to content

GitLab

Commit 751e431a authored by Commandre Benjamin's avatar Commandre Benjamin
Browse files

Supression fichiers inutiles

parent 170c0d8d
Hide whitespace changes
Inline Side-by-side
SimpleExtractionTools/app/otbZonalStatistics_rémi.cxx deleted 100644 → 0
View file @ 170c0d8d
/*=========================================================================
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
SimpleExtractionTools/app/otbZonalStatistics_témoin.cxx deleted 100644 → 0
View file @ 170c0d8d
/*=========================================================================
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] << " ";