otbZonalStatistics.cxx 19.5 KB
Newer Older
1
/*
2
 * Copyright (C) 2017 IRSTEA
3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36
 *
 * 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.
 */

#include "itkFixedArray.h"
#include "itkObjectFactory.h"

// Elevation handler
#include "otbWrapperElevationParametersHandler.h"
#include "otbWrapperApplicationFactory.h"

// Application engine
#include "otbStandardFilterWatcher.h"

// Process objects
#include "otbVectorDataToLabelImageFilter.h"
#include "otbVectorDataIntoImageProjectionFilter.h"
#include "otbStreamingStatisticsMapFromLabelImageFilter.h"
#include "otbStatisticsXMLFileWriter.h"

37 38 39 40
// Raster --> Vector
#include "otbLabelImageToVectorDataFilter.h"
#include "itkBinaryThresholdImageFilter.h"

41
#include "otbFunctorImageFilter.h"
Jordi Inglada's avatar
Jordi Inglada committed
42

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
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;

  /* Typedefs */
59
  typedef Int32ImageType                LabelImageType;
60 61 62
  typedef LabelImageType::ValueType     LabelValueType;
  typedef otb::VectorData<double, 2>    VectorDataType;
  typedef otb::VectorDataIntoImageProjectionFilter<VectorDataType,
63
                                                   FloatVectorImageType>             VectorDataReprojFilterType;
64
  typedef otb::VectorDataToLabelImageFilter<VectorDataType,
65
                                            LabelImageType>                   RasterizeFilterType;
66 67
  typedef VectorDataType::DataTreeType  DataTreeType;
  typedef itk::PreOrderTreeIterator<DataTreeType>
68
  TreeIteratorType;
69 70
  typedef VectorDataType::DataNodeType  DataNodeType;
  typedef DataNodeType::PolygonListPointerType
71
  PolygonListPointerType;
72
  typedef otb::StreamingStatisticsMapFromLabelImageFilter<FloatVectorImageType,
73
                                                          LabelImageType>                   StatsFilterType;
74
  typedef otb::StatisticsXMLFileWriter<FloatVectorImageType::PixelType>
75
  StatsWriterType;
76
  typedef otb::LabelImageToVectorDataFilter<LabelImageType>
77
  LabelImageToVectorFilterType;
78
  typedef itk::BinaryThresholdImageFilter<LabelImageType,
79
                                          LabelImageType>                   ThresholdFilterType;
80

Cédric Traizet's avatar
Cédric Traizet committed
81
  template <class TInput, class TOutput>
Jordi Inglada's avatar
Jordi Inglada committed
82 83 84 85 86 87 88 89
  struct EncoderFunctorType
  {
    StatsFilterType::LabelPopulationMapType* m_CountMap;
    StatsFilterType::PixelValueMapType* m_MeanMap;
    StatsFilterType::PixelValueMapType* m_StdMap;
    StatsFilterType::PixelValueMapType* m_MinMap;
    StatsFilterType::PixelValueMapType* m_MaxMap;
    size_t m_NbInputComponents;
Jordi Inglada's avatar
Jordi Inglada committed
90 91 92 93
    LabelValueType m_InNoData;
    LabelValueType m_OutBvValue;
    static constexpr size_t m_NbStatsPerBand{4};
    static constexpr size_t m_NbGlobalStats{1};
Jordi Inglada's avatar
Jordi Inglada committed
94

Cédric Traizet's avatar
Cédric Traizet committed
95
    size_t OutputSize(const std::array<size_t, 1>&) const
Jordi Inglada's avatar
Jordi Inglada committed
96
    {
Jordi Inglada's avatar
Jordi Inglada committed
97
      return m_NbInputComponents*m_NbStatsPerBand+m_NbGlobalStats;
Jordi Inglada's avatar
Jordi Inglada committed
98 99
    }

100
    size_t OutputSize() const
Jordi Inglada's avatar
Jordi Inglada committed
101
    {
Cédric Traizet's avatar
Cédric Traizet committed
102
      return m_NbInputComponents * m_NbStatsPerBand + m_NbGlobalStats;
103 104
    }

Cédric Traizet's avatar
Cédric Traizet committed
105
    TOutput operator()(TInput const& pix)
106 107
    {
      TOutput outPix(OutputSize());
Jordi Inglada's avatar
Jordi Inglada committed
108 109
      outPix.Fill(m_OutBvValue);
      if(pix != m_InNoData)
Jordi Inglada's avatar
Jordi Inglada committed
110
        {
Jordi Inglada's avatar
Jordi Inglada committed
111 112 113 114 115 116 117 118
        outPix[0] = (*m_CountMap)[pix];
        for(size_t i=0; i<m_NbInputComponents; ++i)
          {
          outPix[i*m_NbStatsPerBand+1] = (*m_MeanMap)[pix][i];
          outPix[i*m_NbStatsPerBand+2] = (*m_StdMap)[pix][i];
          outPix[i*m_NbStatsPerBand+3] = (*m_MinMap)[pix][i];
          outPix[i*m_NbStatsPerBand+4] = (*m_MaxMap)[pix][i];
          }
Jordi Inglada's avatar
Jordi Inglada committed
119 120 121
        }
      return outPix;
    }
Jordi Inglada's avatar
Jordi Inglada committed
122

Jordi Inglada's avatar
Jordi Inglada committed
123 124
    bool operator != (const EncoderFunctorType& other)
    {
Jordi Inglada's avatar
Jordi Inglada committed
125 126 127 128 129 130 131 132 133 134
      if ( m_CountMap != other.m_CountMap||
           m_MeanMap != other.m_MeanMap||
           m_StdMap != other.m_StdMap||
           m_MinMap != other.m_MinMap||
           m_MaxMap != other.m_MaxMap||
           m_NbInputComponents != other.m_NbInputComponents ||
           m_InNoData != other.m_InNoData ||
           m_OutBvValue != other.m_OutBvValue )
        return true;
      else return false;
Jordi Inglada's avatar
Jordi Inglada committed
135 136
    }
  };
Cédric Traizet's avatar
Cédric Traizet committed
137 138
  typedef otb::FunctorImageFilter<EncoderFunctorType<LabelValueType, FloatVectorImageType::PixelType>> EncoderFilterType;

139 140 141 142
  /** Standard macro */
  itkNewMacro(Self);
  itkTypeMacro(ZonalStatistics, Application);

Guillaume Pasero's avatar
Guillaume Pasero committed
143
  void DoInit() override
144 145 146 147 148 149 150
  {

    SetName("ZonalStatistics");
    SetDescription("This application computes zonal statistics");

    // Documentation
    SetDocLongDescription("This application computes zonal statistics from label image, or vector data. "
151 152 153 154 155 156 157 158 159 160
                          "The application inputs one input multiband image, and another input for zones definition. "
                          "Zones can be defined with a label image (inzone.labelimage.in) or a vector data layer "
                          "(inzone.vector.in). The following statistics are computed over each zones: mean, min, max, "
                          "and standard deviation. Statistics can be exported in a vector layer (if the input zone "
                          "definition is a label image, it will be vectorized) or in a XML file");
    SetDocLimitations("1) The inzone.vector.in must fit in memory (if \"inzone\" is \"vector\"). 2) The vectorized label "
                      "image must also fit in memory (if \"out\" is \"vector\"): if not, consider using \"out\" to "
                      "\"xml\"");
    SetDocAuthors("Remi Cresson, Jordi Inglada");
    SetDocSeeAlso("ComputeImagesStatistics");
161 162 163 164 165

    AddDocTag(Tags::Manip);
    AddDocTag(Tags::Analysis);

    // Input image
166 167
    AddParameter(ParameterType_InputImage, "in",   "Input Image");
    AddParameter(ParameterType_Float,      "inbv", "Background value to ignore in statistics computation");
168
    MandatoryOff                          ("inbv");
169

170 171 172 173
    // Input zone mode
    AddParameter(ParameterType_Choice, "inzone", "Type of input for the zone definitions");
    AddChoice("inzone.vector",     "Input objects from vector data");
    AddChoice("inzone.labelimage", "Input objects from label image");
174 175

    // Input for vector mode
176 177 178
    AddParameter(ParameterType_InputVectorData,  "inzone.vector.in",        "Input vector data");
    AddParameter(ParameterType_Bool,             "inzone.vector.reproject", "Reproject the input vector");

179 180

    // Input for label image mode
181 182 183 184 185 186 187 188 189 190
    AddParameter(ParameterType_InputImage, "inzone.labelimage.in",     "Input label image");
    AddParameter(ParameterType_Int,        "inzone.labelimage.nodata", "No-data value for the input label image");
    MandatoryOff                          ("inzone.labelimage.nodata");


    // Output stats mode
    AddParameter(ParameterType_Choice, "out", "Format of the output stats");
    AddChoice("out.vector", "Output vector data");
    AddParameter(ParameterType_OutputVectorData, "out.vector.filename", "Filename for the output vector data");
    AddChoice("out.xml", "Output XML file");
191
    AddParameter(ParameterType_String, "out.xml.filename", "Filename for the output xml file");
Jordi Inglada's avatar
Jordi Inglada committed
192
   AddChoice("out.raster", "Output raster image");
Jordi Inglada's avatar
Jordi Inglada committed
193 194 195
   AddParameter(ParameterType_OutputImage, "out.raster.filename", "File name for the raster image");    
   AddParameter(ParameterType_Float,      "out.raster.bv", "Background value for the output raster");
   MandatoryOff                          ("out.raster.bv");
196

Jordi Inglada's avatar
Jordi Inglada committed
197
   AddRAMParameter();
198 199 200

    // Doc example parameter settings
    SetDocExampleParameterValue("in", "input.tif");
201 202
    SetDocExampleParameterValue("inzone.vector.in", "myvector.shp");
    SetDocExampleParameterValue("out.vector.filename", "myvector_with_stats.shp");
203

204
    SetOfficialDocLink();
205 206
  }

Guillaume Pasero's avatar
Guillaume Pasero committed
207
  void DoUpdateParameters() override
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
  {
    // Nothing to do here : all parameters are independent
  }

  // Returns a string of the kind "prefix_i"
  const std::string CreateFieldName(const std::string & prefix, const unsigned int i)
  {
    std::stringstream ss;
    ss << prefix << "_" << i;
    return ss.str();
  }

  // Returns a null pixel which has the same number of components per pixels as img
  FloatVectorImageType::PixelType NullPixel(FloatVectorImageType::Pointer & img)
  {
    const unsigned int nBands = img->GetNumberOfComponentsPerPixel();
224
    FloatVectorImageType::PixelType pix(nBands);
225 226 227 228
    pix.Fill(0);
    return pix;
  }

Jordi Inglada's avatar
Jordi Inglada committed
229 230 231 232 233 234 235 236 237
  void GetStats()
  {
    m_CountMap = m_StatsFilter->GetLabelPopulationMap();
    m_MeanMap = m_StatsFilter->GetMeanValueMap();
    m_StdMap = m_StatsFilter->GetStandardDeviationValueMap();
    m_MinMap = m_StatsFilter->GetMinValueMap();
    m_MaxMap = m_StatsFilter->GetMaxValueMap();
  }

238
  void PrepareForLabelImageInput()
239
  {
240 241 242 243 244 245 246
    otbAppLogINFO("Zone definition: label image");
    // Computing stats
    m_StatsFilter->SetInputLabelImage(GetParameterInt32Image("inzone.labelimage.in"));
    m_StatsFilter->Update();
    // In this zone definition mode, the user can provide a no-data value for the labels
    if (HasUserValue("inzone.labelimage.nodata"))
      m_IntNoData = GetParameterInt("inzone.labelimage.nodata");
Jordi Inglada's avatar
Jordi Inglada committed
247
    GetStats();
248
  }
249

250 251 252 253
  void PrepareForVectorDataInput()
  {
    otbAppLogINFO("Zone definition: vector");
    otbAppLogINFO("Loading vector data...");
Jordi Inglada's avatar
Jordi Inglada committed
254
    m_VectorDataSrc = GetParameterVectorData("inzone.vector.in");
255 256 257 258 259 260
    // Reproject vector data
    if (GetParameterInt("inzone.vector.reproject") != 0)
      {
      ReprojectVectorDataIntoInputImage();
      }
    RasterizeInputVectorData();
261 262 263
    // Computing stats
    m_StatsFilter->SetInputLabelImage(m_RasterizeFilter->GetOutput());
    m_StatsFilter->Update();
Jordi Inglada's avatar
Jordi Inglada committed
264
    GetStats();
265 266 267 268 269 270
  }

  void ReprojectVectorDataIntoInputImage()
  {
    otbAppLogINFO("Vector data reprojection enabled");
    m_VectorDataReprojectionFilter = VectorDataReprojFilterType::New();
Jordi Inglada's avatar
Jordi Inglada committed
271
    m_VectorDataReprojectionFilter->SetInputVectorData(m_VectorDataSrc.GetPointer());
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290
    m_VectorDataReprojectionFilter->SetInputImage(m_InputImage);
    AddProcess(m_VectorDataReprojectionFilter, "Reproject vector data");
    m_VectorDataReprojectionFilter->Update();
    m_VectorDataSrc = m_VectorDataReprojectionFilter->GetOutput();
  }

  void RasterizeInputVectorData()
  {
    // Rasterize vector data
    m_RasterizeFilter = RasterizeFilterType::New();
    m_RasterizeFilter->AddVectorData(m_VectorDataSrc);
    m_RasterizeFilter->SetOutputOrigin(m_InputImage->GetOrigin());
    m_RasterizeFilter->SetOutputSpacing(m_InputImage->GetSignedSpacing());
    m_RasterizeFilter->SetOutputSize(m_InputImage->GetLargestPossibleRegion().GetSize());
    m_RasterizeFilter->SetOutputProjectionRef(m_InputImage->GetProjectionRef());
    m_RasterizeFilter->SetBurnAttribute("________");
    m_RasterizeFilter->SetDefaultBurnValue(0);
    m_RasterizeFilter->SetGlobalWarningDisplay(false);
    m_RasterizeFilter->SetBackgroundValue(m_IntNoData);
Jordi Inglada's avatar
Jordi Inglada committed
291
    AddProcess(m_RasterizeFilter, "Rasterize input vector data");
292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309
  }

  void RemoveNoDataEntry()
  {
    if (( GetParameterAsString("inzone") == "labelimage" && HasUserValue("inzone.labelimage.nodata"))
        || (GetParameterAsString("inzone") == "vector")      )
      {
      otbAppLogINFO("Removing entries for label value " << m_IntNoData);

      m_CountMap.erase(m_IntNoData);
      m_MeanMap.erase(m_IntNoData);
      m_StdMap.erase(m_IntNoData);
      m_MinMap.erase(m_IntNoData);
      m_MaxMap.erase(m_IntNoData);
      }
  }


Jordi Inglada's avatar
Jordi Inglada committed
310
  void GenerateVectorDataFromLabelImage()
Jordi Inglada's avatar
Jordi Inglada committed
311 312 313 314 315 316 317 318 319 320
  {    
    // Mask for label image
    m_InputThresholdFilter = ThresholdFilterType::New();
    m_InputThresholdFilter->SetInput(GetParameterInt32Image("inzone.labelimage.in"));
    m_InputThresholdFilter->SetInsideValue(0);
    m_InputThresholdFilter->SetOutsideValue(1);
    m_InputThresholdFilter->SetLowerThreshold(m_IntNoData);
    m_InputThresholdFilter->SetUpperThreshold(m_IntNoData);
    m_InputThresholdFilter->UpdateOutputInformation();
    AddProcess(m_InputThresholdFilter, "Threshold label image");
321 322 323
    // Vectorize the image
    m_LabelImageToVectorFilter = LabelImageToVectorFilterType::New();
    m_LabelImageToVectorFilter->SetInput(GetParameterInt32Image("inzone.labelimage.in"));
Jordi Inglada's avatar
Jordi Inglada committed
324
    m_LabelImageToVectorFilter->SetInputMask(m_InputThresholdFilter->GetOutput());
325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375
    m_LabelImageToVectorFilter->SetFieldName("polygon_id");
    AddProcess(m_LabelImageToVectorFilter, "Vectorize label image");
    m_LabelImageToVectorFilter->Update();
    // The source vector data is the vectorized label image
    m_VectorDataSrc = m_LabelImageToVectorFilter->GetOutput();
  }

  void WriteVectorData()
  {
    // Add statistics fields
    otbAppLogINFO("Writing output vector data");
    LabelValueType internalFID = -1;
    m_NewVectorData = VectorDataType::New();
    DataNodeType::Pointer root = m_NewVectorData->GetDataTree()->GetRoot()->Get();
    DataNodeType::Pointer document = DataNodeType::New();
    document->SetNodeType(otb::DOCUMENT);
    m_NewVectorData->GetDataTree()->Add(document, root);
    DataNodeType::Pointer folder = DataNodeType::New();
    folder->SetNodeType(otb::FOLDER);
    m_NewVectorData->GetDataTree()->Add(folder, document);
    m_NewVectorData->SetProjectionRef(m_VectorDataSrc->GetProjectionRef());
    TreeIteratorType itVector(m_VectorDataSrc->GetDataTree());
    itVector.GoToBegin();

    while (!itVector.IsAtEnd())
      {
      if (!itVector.Get()->IsRoot() && !itVector.Get()->IsDocument() && !itVector.Get()->IsFolder())
        {

        DataNodeType::Pointer currentGeometry = itVector.Get();
        if (m_FromLabelImage)
          internalFID = currentGeometry->GetFieldAsInt("polygon_id");
        else
          internalFID++;

        // Add the geometry with the new fields
        if (m_CountMap.count(internalFID) > 0)
          {
          currentGeometry->SetFieldAsDouble("count", m_CountMap[internalFID] );
          for (unsigned int band = 0 ; band < m_InputImage->GetNumberOfComponentsPerPixel() ; band++)
            {
            currentGeometry->SetFieldAsDouble(CreateFieldName("mean",  band), m_MeanMap[internalFID][band] );
            currentGeometry->SetFieldAsDouble(CreateFieldName("stdev", band), m_StdMap [internalFID][band] );
            currentGeometry->SetFieldAsDouble(CreateFieldName("min",   band), m_MinMap [internalFID][band] );
            currentGeometry->SetFieldAsDouble(CreateFieldName("max",   band), m_MaxMap [internalFID][band] );
            }
          m_NewVectorData->GetDataTree()->Add(currentGeometry, folder);
          }
        }
      ++itVector;
      } // next feature
376

377 378 379
    SetParameterOutputVectorData("out.vector.filename", m_NewVectorData);
  }

Jordi Inglada's avatar
Jordi Inglada committed
380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395
  void SetOutBvValue()
  {
    if (HasUserValue("out.raster.bv"))
      {
      m_OutBvValue = GetParameterFloat("out.raster.bv");
      }
    else if(HasUserValue("inbv"))
      {
      m_OutBvValue = GetParameterFloat("inbv");
      }
    else
      {
      m_OutBvValue = m_IntNoData;
      }
  }

Jordi Inglada's avatar
Jordi Inglada committed
396
  void WriteRasterData()
397
  {
Jordi Inglada's avatar
Jordi Inglada committed
398 399
    
    otbAppLogINFO("Writing output raster data");
Jordi Inglada's avatar
Jordi Inglada committed
400 401
    SetOutBvValue();
    m_OutputThresholdFilter = ThresholdFilterType::New();
Jordi Inglada's avatar
Jordi Inglada committed
402 403
    m_EncoderFilter = EncoderFilterType::New();
    if(m_FromLabelImage)
Jordi Inglada's avatar
Jordi Inglada committed
404 405 406
        {
        m_EncoderFilter->SetInput(GetParameterInt32Image("inzone.labelimage.in"));
        }
407 408
    else
      {
Jordi Inglada's avatar
Jordi Inglada committed
409
      m_EncoderFilter->SetInput(m_RasterizeFilter->GetOutput());
410
      }
Jordi Inglada's avatar
Jordi Inglada committed
411

Cédric Traizet's avatar
Cédric Traizet committed
412 413 414 415 416 417 418 419 420 421 422 423
      m_EncoderFilter->GetModifiableFunctor().m_CountMap          = &m_CountMap;
      m_EncoderFilter->GetModifiableFunctor().m_MeanMap           = &m_MeanMap;
      m_EncoderFilter->GetModifiableFunctor().m_StdMap            = &m_StdMap;
      m_EncoderFilter->GetModifiableFunctor().m_MinMap            = &m_MinMap;
      m_EncoderFilter->GetModifiableFunctor().m_MaxMap            = &m_MaxMap;
      m_EncoderFilter->GetModifiableFunctor().m_NbInputComponents = m_InputImage->GetNumberOfComponentsPerPixel();
      m_EncoderFilter->GetModifiableFunctor().m_InNoData          = m_IntNoData;
      m_EncoderFilter->GetModifiableFunctor().m_OutBvValue        = m_OutBvValue;

      otbAppLogINFO("Output raster image will have " << (m_EncoderFilter->GetFunctor()).OutputSize() << " bands\n");
      AddProcess(m_EncoderFilter, "Encode output raster image");
      SetParameterOutputImage("out.raster.filename", m_EncoderFilter->GetOutput());
Jordi Inglada's avatar
Jordi Inglada committed
424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440
  }

  void WriteXMLStatsFile()
  {
    // Write statistics in XML file
    const std::string outXMLFile = this->GetParameterString("out.xml.filename");
    otbAppLogINFO("Writing " + outXMLFile)
        StatsWriterType::Pointer statWriter = StatsWriterType::New();
      statWriter->SetFileName(outXMLFile);
      statWriter->AddInputMap<StatsFilterType::LabelPopulationMapType>("count",m_CountMap);
      statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("mean",m_MeanMap);
      statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("std",m_StdMap);
      statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("min",m_MinMap);
      statWriter->AddInputMap<StatsFilterType::PixelValueMapType>("max",m_MaxMap);
      statWriter->Update();
    }

Guillaume Pasero's avatar
Guillaume Pasero committed
441
    void DoExecute() override
Jordi Inglada's avatar
Jordi Inglada committed
442 443 444 445 446 447 448
    {
      // Get input image
      m_InputImage = GetParameterImage("in");
      // Statistics filter
      m_StatsFilter = StatsFilterType::New();
      m_StatsFilter->SetInput(m_InputImage);
      if (HasUserValue("inbv"))
449
        {
Jordi Inglada's avatar
Jordi Inglada committed
450 451
        m_StatsFilter->SetUseNoDataValue(true);
        m_StatsFilter->SetNoDataValue(GetParameterFloat("inbv"));
452
        }
Jordi Inglada's avatar
Jordi Inglada committed
453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471
     m_StatsFilter->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
     AddProcess(m_StatsFilter->GetStreamer(), "Computing statistics");
     // Select zone definition mode
     m_FromLabelImage = (GetParameterAsString("inzone") == "labelimage");
     if (m_FromLabelImage) PrepareForLabelImageInput();
     else if (GetParameterAsString("inzone") == "vector") PrepareForVectorDataInput();
     else otbAppLogFATAL("Unknown zone definition mode");
     // Remove the no-data entry
     RemoveNoDataEntry();
     // Generate output
     if (GetParameterAsString("out") == "xml") WriteXMLStatsFile();
     else     // vector or raster
        {
        if (m_FromLabelImage) GenerateVectorDataFromLabelImage();
        if (GetParameterAsString("out") == "vector") WriteVectorData(); 
        else if (GetParameterAsString("out") == "raster") WriteRasterData();
        else otbAppLogFATAL("Unknown output mode");
        }
    }
472 473 474 475 476 477

  VectorDataType::Pointer m_VectorDataSrc;
  VectorDataType::Pointer m_NewVectorData;
  VectorDataReprojFilterType::Pointer m_VectorDataReprojectionFilter;
  RasterizeFilterType::Pointer m_RasterizeFilter;
  StatsFilterType::Pointer m_StatsFilter;
478
  LabelImageToVectorFilterType::Pointer m_LabelImageToVectorFilter;
Jordi Inglada's avatar
Jordi Inglada committed
479 480
  ThresholdFilterType::Pointer m_InputThresholdFilter;
  ThresholdFilterType::Pointer m_OutputThresholdFilter;
481
  FloatVectorImageType::Pointer m_InputImage;
482
  LabelValueType m_IntNoData = itk::NumericTraits<LabelValueType>::max();
Jordi Inglada's avatar
Jordi Inglada committed
483
  LabelValueType m_OutBvValue = itk::NumericTraits<LabelValueType>::max(); 
484
  bool m_FromLabelImage = false;
485 486 487 488 489
  StatsFilterType::LabelPopulationMapType m_CountMap;
  StatsFilterType::PixelValueMapType m_MeanMap;
  StatsFilterType::PixelValueMapType m_StdMap;
  StatsFilterType::PixelValueMapType m_MinMap;
  StatsFilterType::PixelValueMapType m_MaxMap;
Jordi Inglada's avatar
Jordi Inglada committed
490
  EncoderFilterType::Pointer m_EncoderFilter;
491 492 493 494 495 496

};
}
}

OTB_APPLICATION_EXPORT( otb::Wrapper::ZonalStatistics )