diff --git a/Applications/CMakeLists.txt b/Applications/CMakeLists.txt index c37914160f3ed4ce717bf7565b90df85b193212c..a882f352e2411a978d8af6d62a6d1fd457aae7b0 100644 --- a/Applications/CMakeLists.txt +++ b/Applications/CMakeLists.txt @@ -1,4 +1,5 @@ +add_subdirectory(Classification) add_subdirectory(Hyperspectral) add_subdirectory(Util) diff --git a/Applications/Classification/CMakeLists.txt b/Applications/Classification/CMakeLists.txt new file mode 100644 index 0000000000000000000000000000000000000000..a07b502344cce9f4fe37f67d33f1b67dc69c8843 --- /dev/null +++ b/Applications/Classification/CMakeLists.txt @@ -0,0 +1,6 @@ +include(WrapperMacros) + +OTB_CREATE_APPLICATION(NAME EstimateImagesStatistics + SOURCES otbEstimateImagesStatistics.cxx + LINK_LIBRARIES OTBIO;OTBCommon;OTBBasicFilters) + diff --git a/Applications/Classification/otbEstimateImagesStatistics.cxx b/Applications/Classification/otbEstimateImagesStatistics.cxx new file mode 100644 index 0000000000000000000000000000000000000000..9f2df5bdc82d79585a6c38837d2aeec06b030c74 --- /dev/null +++ b/Applications/Classification/otbEstimateImagesStatistics.cxx @@ -0,0 +1,152 @@ +/*========================================================================= + + Program: ORFEO Toolbox + Language: C++ + Date: $Date$ + Version: $Revision$ + + + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved. + See OTBCopyright.txt for details. + + + 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 "otbWrapperApplication.h" +#include "otbWrapperApplicationFactory.h" + +#include "otbStatisticsXMLFileWriter.h" +#include "otbStreamingStatisticsVectorImageFilter.h" + +namespace otb +{ +namespace Wrapper +{ + +class EstimateImageStatistics: public Application +{ +public: + /** Standard class typedefs. */ + typedef EstimateImageStatistics Self; + typedef Application Superclass; + typedef itk::SmartPointer<Self> Pointer; + typedef itk::SmartPointer<const Self> ConstPointer; + + /** Standard macro */ + itkNewMacro(Self); + + itkTypeMacro(EstimateImageStatistics, otb::Application); + +private: + EstimateImageStatistics() + { + SetName("EstimateImageStatistics"); + SetDescription("Estimate mean/standard deviation for all images in the input list"); + } + + virtual ~EstimateImageStatistics() + { + } + + void DoCreateParameters() + { + AddParameter(ParameterType_InputImageList, "in", "Input Image"); + AddParameter(ParameterType_Filename, "out", "Output xml file"); + } + + void DoUpdateParameters() + { + // Nothing to do here : all parameters are independent + } + + void DoExecute() + { + //Statistics estimator + typedef otb::StreamingStatisticsVectorImageFilter<FloatVectorImageType> StreamingStatisticsVImageFilterType; + + // Samples + typedef double ValueType; + typedef itk::VariableLengthVector<ValueType> MeasurementType; + + unsigned int nbSamples = 0; + unsigned int nbBands = 0; + + // Build a Measurement Vector of mean + MeasurementType mean; + + // Build a MeasurementVector of variance + MeasurementType variance; + + FloatVectorImageListType* imageList = GetParameterImageList("in"); + + //Iterate over all input images + for (unsigned int imageId = 0; imageId < imageList->Size(); ++imageId) + { + FloatVectorImageType* image = imageList->GetNthElement(imageId); + + if (nbBands == 0) + { + nbBands = image->GetNumberOfComponentsPerPixel(); + } + else if (nbBands != image->GetNumberOfComponentsPerPixel()) + { + itkExceptionMacro(<< "The image #" << imageId << " has " << image->GetNumberOfComponentsPerPixel() + << " bands, while the first one has " << nbBands ); + } + + FloatVectorImageType::SizeType size = image->GetLargestPossibleRegion().GetSize(); + + //Set the measurement vectors size if it's the first iteration + if (imageId == 0) + { + mean.SetSize(nbBands); + mean.Fill(0.); + variance.SetSize(nbBands); + variance.Fill(0.); + } + + // Compute Statistics of each VectorImage + StreamingStatisticsVImageFilterType::Pointer statsEstimator = StreamingStatisticsVImageFilterType::New(); + statsEstimator->SetInput(image); + statsEstimator->Update(); + mean += statsEstimator->GetMean(); + for (unsigned int itBand = 0; itBand < nbBands; itBand++) + { + variance[itBand] += (size[0] * size[1] - 1) * (statsEstimator->GetCovariance())(itBand, itBand); + } + //Increment nbSamples + nbSamples += size[0] * size[1] * nbBands; + } + + //Divide by the number of input images to get the mean over all layers + mean /= imageList->Size(); + //Apply the pooled variance formula + variance /= (nbSamples - imageList->Size()); + + MeasurementType stddev; + stddev.SetSize(nbBands); + stddev.Fill(0.); + for (unsigned int i = 0; i < variance.GetSize(); ++i) + { + stddev[i] = vcl_sqrt(variance[i]); + } + + // Write the Statistics via the statistic writer + typedef otb::StatisticsXMLFileWriter<MeasurementType> StatisticsWriter; + StatisticsWriter::Pointer writer = StatisticsWriter::New(); + writer->SetFileName(GetParameterString("out")); + writer->AddInput("mean", mean); + writer->AddInput("stddev", stddev); + writer->Update(); + } + + itk::LightObject::Pointer m_FilterRef; +}; + +} +} + +OTB_APPLICATION_EXPORT(otb::Wrapper::EstimateImageStatistics)