otbComputeImagesStatistics.cxx 8.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
 * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
 *
 * 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.
 */
20 21 22 23 24 25

#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"

#include "otbStatisticsXMLFileWriter.h"
#include "otbStreamingStatisticsVectorImageFilter.h"
26
#include <sstream>
27 28 29 30 31 32

namespace otb
{
namespace Wrapper
{

33
class ComputeImagesStatistics: public Application
34 35 36
{
public:
  /** Standard class typedefs. */
37
  typedef ComputeImagesStatistics Self;
38 39 40 41 42 43 44
  typedef Application Superclass;
  typedef itk::SmartPointer<Self> Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

  /** Standard macro */
  itkNewMacro(Self);

45
  itkTypeMacro(ComputeImagesStatistics, otb::Application);
46 47

private:
48
  void DoInit() override
49
  {
50
    SetName("ComputeImagesStatistics");
51
    SetDocName("Compute Images second order statistics");
52 53 54 55
    SetDescription("Computes global mean and standard deviation for each band "
      "from a set of images and optionally saves the results in an XML file.");
    SetDocLongDescription("This application computes a global mean and standard deviation "
      "for each band of a set of images and optionally saves the results in an XML file."
56
      " The output XML is intended to be used as an input "
57 58 59 60 61
      "for the TrainImagesClassifier application to normalize samples before learning. "
      "You can also normalize the image with the XML file in the ImageClassifier application.");

    SetDocLimitations("Each image of the set must contain the same bands as the others"
                      " (i.e. same types, in the same order).");
62
    SetDocAuthors("OTB-Team");
63
    SetDocSeeAlso("Documentation of the TrainImagesClassifier and ImageClassifier application.");
64

65
    AddDocTag(Tags::Learning);
66
    AddDocTag(Tags::Analysis);
67

68
    AddParameter(ParameterType_InputImageList, "il", "Input images");
69
    SetParameterDescription( "il", "List of input image filenames." );
70

71
    AddParameter(ParameterType_Float, "bv", "Background Value");
72
    SetParameterDescription( "bv", "Background value to ignore in computation of statistics." );
73 74
    MandatoryOff("bv");

75
    AddParameter(ParameterType_OutputFilename, "out", "Output XML file");
76
    SetParameterDescription( "out", "XML filename where the statistics are saved for future reuse." );
77
    MandatoryOff("out");
78

79 80
    AddRAMParameter();

81
   // Doc example parameter settings
82 83
   SetDocExampleParameterValue("il", "QB_1_ortho.tif");
   SetDocExampleParameterValue("out", "EstimateImageStatisticsQB1.xml");
84

85
   SetOfficialDocLink();
86 87
  }

88
  void DoUpdateParameters() override
89 90 91 92
  {
    // Nothing to do here : all parameters are independent
  }

93
  void DoExecute() override
94 95 96 97 98 99 100
  {
    //Statistics estimator
    typedef otb::StreamingStatisticsVectorImageFilter<FloatVectorImageType> StreamingStatisticsVImageFilterType;

    // Samples
    typedef double ValueType;
    typedef itk::VariableLengthVector<ValueType> MeasurementType;
101
    typedef itk::VariableSizeMatrix<ValueType> MatrixValueType;
102 103 104

    unsigned int nbBands = 0;

105 106 107 108 109 110 111
    FloatVectorImageListType* imageList = GetParameterImageList("il");
    FloatVectorImageListType::InternalContainerSizeType nbImages = imageList->Size();

    // Initialization, all image have same size and number of band/component
    FloatVectorImageType* firstImage = imageList->GetNthElement(0);
    nbBands = firstImage->GetNumberOfComponentsPerPixel();

112
    // Build a Measurement Vector of mean
113 114
    MatrixValueType mean(nbBands, static_cast<unsigned int>(nbImages));
    mean.Fill(itk::NumericTraits<MatrixValueType::ValueType>::Zero);
115

116 117 118
    // Build a Measurement Matrix of variance
    MatrixValueType variance(nbBands, static_cast<unsigned int>(nbImages));
    variance.Fill(itk::NumericTraits<MatrixValueType::ValueType>::Zero);
119

120 121 122
    // Build a Measurement Matrix of nbSamples
    MatrixValueType nbSamples(nbBands, static_cast<unsigned int>(nbImages));
    nbSamples.Fill(itk::NumericTraits<MatrixValueType::ValueType>::Zero);
123 124

    //Iterate over all input images
125
    for (unsigned int imageId = 0; imageId < nbImages; ++imageId)
126 127
      {
      FloatVectorImageType* image = imageList->GetNthElement(imageId);
128
      if (nbBands != image->GetNumberOfComponentsPerPixel())
129
        {
130 131
        itkExceptionMacro(<< "The image #" << imageId + 1 << " has " << image->GetNumberOfComponentsPerPixel()
            << " bands, while the image #1 has " << nbBands );
132 133 134 135
        }

      // Compute Statistics of each VectorImage
      StreamingStatisticsVImageFilterType::Pointer statsEstimator = StreamingStatisticsVImageFilterType::New();
136 137
      std::ostringstream processName;
      processName << "Processing Image (" << imageId+1 << "/" << imageList->Size() << ")";
138
      AddProcess(statsEstimator->GetStreamer(), processName.str());
139
      statsEstimator->SetInput(image);
140
      statsEstimator->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
141 142

      if( HasValue( "bv" ) )
143 144 145 146
        {
        statsEstimator->SetIgnoreUserDefinedValue(true);
        statsEstimator->SetUserIgnoredValue(GetParameterFloat("bv"));
        }
147
      statsEstimator->Update();
148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175

      MeasurementType nbRelevantPixels = statsEstimator->GetNbRelevantPixels();
      MeasurementType meanPerBand = statsEstimator->GetMean();

      for(unsigned int itBand = 0; itBand < nbBands; itBand++)
        {
        mean(itBand, imageId) = meanPerBand[itBand];
        variance(itBand, imageId) = (statsEstimator->GetCovariance())( itBand, itBand );
        nbSamples(itBand, imageId) = nbRelevantPixels[itBand];
        }
      }

    // Compute total mean and pooled variation for each band of the image list
    MeasurementType totalSamplesPerBand;
    totalSamplesPerBand.SetSize(nbBands);
    totalSamplesPerBand.Fill(itk::NumericTraits<MeasurementType::ValueType>::Zero);

    MeasurementType totalMeanPerBand;
    totalMeanPerBand.SetSize(nbBands);
    totalMeanPerBand.Fill(itk::NumericTraits<MeasurementType::ValueType>::Zero);

    MeasurementType totalVariancePerBand;
    totalVariancePerBand.SetSize(nbBands);
    totalVariancePerBand.Fill(itk::NumericTraits<MeasurementType::ValueType>::Zero);

    for (unsigned int imageId = 0; imageId < nbImages; ++imageId)
      {
      for(unsigned int itBand = 0; itBand < nbBands; itBand++)
176
        {
177 178 179 180
        MeasurementType::ValueType nbSample = nbSamples(itBand, imageId);
        totalSamplesPerBand[itBand] += nbSample;
        totalMeanPerBand[itBand] += mean(itBand, imageId) * nbSample;
        totalVariancePerBand[itBand] += variance(itBand, imageId) * (nbSample  - 1);
181 182 183
        }
      }

184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206
    // Check 0 division
    for(unsigned int itBand = 0; itBand < nbBands; itBand++)
      {
      MeasurementType::ValueType nbSample = totalSamplesPerBand[itBand];

      if ( nbSample > nbImages )
        {
        totalVariancePerBand[itBand] /= (nbSample - nbImages);
        }
      else
        {
        totalVariancePerBand[itBand] = itk::NumericTraits<ValueType>::Zero;
        }

      if ( nbSample != 0 )
        {
        totalMeanPerBand[itBand] /= nbSample;
        }
      else
        {
        totalMeanPerBand[itBand] = itk::NumericTraits<ValueType>::Zero;
        }
      }
207 208 209

    MeasurementType stddev;
    stddev.SetSize(nbBands);
210 211
    stddev.Fill(itk::NumericTraits<MeasurementType::ValueType>::Zero);
    for (unsigned int i = 0; i < totalVariancePerBand.GetSize(); ++i)
212
      {
213
      stddev[i] = std::sqrt(totalVariancePerBand[i]);
214 215
      }

216
    if( HasValue( "out" ) )
217 218 219 220 221
      {
      // Write the Statistics via the statistic writer
      typedef otb::StatisticsXMLFileWriter<MeasurementType> StatisticsWriter;
      StatisticsWriter::Pointer writer = StatisticsWriter::New();
      writer->SetFileName(GetParameterString("out"));
222
      writer->AddInput("mean", totalMeanPerBand);
223 224 225 226 227
      writer->AddInput("stddev", stddev);
      writer->Update();
      }
    else
      {
228 229
      otbAppLogINFO("Mean: "<<mean<<std::endl);
      otbAppLogINFO("Standard Deviation: "<<stddev<<std::endl);
230
      }
231 232 233 234 235 236 237 238
  }

  itk::LightObject::Pointer m_FilterRef;
};

}
}

239
OTB_APPLICATION_EXPORT(otb::Wrapper::ComputeImagesStatistics)