otbPolygonClassStatistics.cxx 9.44 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 "otbOGRDataToClassStatisticsFilter.h"
#include "otbStatisticsXMLFileWriter.h"
26 27
#include "otbGeometriesProjectionFilter.h"
#include "otbGeometriesSet.h"
28
#include "otbWrapperElevationParametersHandler.h"
29 30 31 32 33 34

namespace otb
{
namespace Wrapper
{

35 36 37 38 39 40
/** Utility function to negate std::isalnum */
bool IsNotAlphaNum(char c)
  {
  return !std::isalnum(c);
  }

41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59
class PolygonClassStatistics : public Application
{
public:
  /** Standard class typedefs. */
  typedef PolygonClassStatistics        Self;
  typedef Application                   Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

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

  itkTypeMacro(PolygonClassStatistics, otb::Application);

  /** Filters typedef */
  typedef otb::OGRDataToClassStatisticsFilter<FloatVectorImageType,UInt8ImageType> FilterType;
  
  typedef otb::StatisticsXMLFileWriter<FloatVectorImageType::PixelType> StatWriterType;

60 61 62 63
  typedef otb::GeometriesSet GeometriesType;

  typedef otb::GeometriesProjectionFilter ProjectionFilterType;

64 65 66 67 68 69
private:
  PolygonClassStatistics()
    {
   
    }

70
  void DoInit() override
71 72 73 74 75 76
  {
    SetName("PolygonClassStatistics");
    SetDescription("Computes statistics on a training polygon set.");

    // Documentation
    SetDocName("Polygon Class Statistics");
Victor Poughon's avatar
Victor Poughon committed
77 78 79 80 81
    SetDocLongDescription("Process a set of geometries intended for training (they should have a field giving the associated "
      "class). The geometries are analyzed against a support image to compute statistics:\n\n"
      "* Number of samples per class\n"
      "* Number of samples per geometry\n\n"

82
      "An optional raster mask can be used to discard samples. Different types"
Victor Poughon's avatar
Victor Poughon committed
83 84 85 86 87 88
      " of geometry are supported: polygons, lines, points. The behaviour is "
      "different for each type of geometry:\n\n"
      "* Polygon: select pixels whose center is inside the polygon\n"
      "* Lines: select pixels intersecting the line\n"
      "* Points: select closest pixel to the point");

89 90 91 92 93 94
    SetDocLimitations("None");
    SetDocAuthors("OTB-Team");
    SetDocSeeAlso(" ");

    AddDocTag(Tags::Learning);

95
    AddParameter(ParameterType_InputImage,  "in",   "Input image");
96 97
    SetParameterDescription("in", "Support image that will be classified");

98
    AddParameter(ParameterType_InputImage,  "mask",   "Input validity mask");
99 100 101 102
    SetParameterDescription("mask", "Validity mask (only pixels corresponding to a mask value greater than 0 will be used for statistics)");
    MandatoryOff("mask");
    
    AddParameter(ParameterType_InputFilename, "vec", "Input vectors");
103
    SetParameterDescription("vec","Input geometries to analyze");
104
    
105
    AddParameter(ParameterType_OutputFilename, "out", "Output XML statistics file");
106 107
    SetParameterDescription("out","Output file to store statistics (XML format)");

108
    AddParameter(ParameterType_ListView, "field", "Field Name");
109
    SetParameterDescription("field","Name of the field carrying the class name in the input vectors.");
110
    SetListViewSingleSelectionMode("field",true);
111 112 113 114 115
    
    AddParameter(ParameterType_Int, "layer", "Layer Index");
    SetParameterDescription("layer", "Layer index to read in the input vector file.");
    MandatoryOff("layer");
    SetDefaultParameterInt("layer",0);
116 117 118

    ElevationParametersHandler::AddElevationParameters(this, "elev");

119 120 121
    AddRAMParameter();

    // Doc example parameter settings
122 123
    SetDocExampleParameterValue("in", "support_image.tif");
    SetDocExampleParameterValue("vec", "variousVectors.sqlite");
124
    SetDocExampleParameterValue("field", "CLASS");
125
    SetDocExampleParameterValue("out","polygonStat.xml");
126

127
    SetOfficialDocLink();
128 129
  }

130
  void DoUpdateParameters() override
131
  {
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150
     if ( HasValue("vec") )
      {
      std::string vectorFile = GetParameterString("vec");
      ogr::DataSource::Pointer ogrDS =
        ogr::DataSource::New(vectorFile, ogr::DataSource::Modes::Read);
      ogr::Layer layer = ogrDS->GetLayer(this->GetParameterInt("layer"));
      ogr::Feature feature = layer.ogr().GetNextFeature();

      ClearChoices("field");
      
      for(int iField=0; iField<feature.ogr().GetFieldCount(); iField++)
        {
        std::string key, item = feature.ogr().GetFieldDefnRef(iField)->GetNameRef();
        key = item;
        std::string::iterator end = std::remove_if(key.begin(),key.end(),IsNotAlphaNum);
        std::transform(key.begin(), end, key.begin(), tolower);
        
        OGRFieldType fieldType = feature.ogr().GetFieldDefnRef(iField)->GetType();
        
151
        if(fieldType == OFTString || fieldType == OFTInteger || ogr::version_proxy::IsOFTInteger64(fieldType))
152 153 154 155 156 157
          {
          std::string tmpKey="field."+key.substr(0, end - key.begin());
          AddChoice(tmpKey,item);
          }
        }
      }
158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173

     // Check that the extension of the output parameter is XML (mandatory for
     // StatisticsXMLFileWriter)
     // Check it here to trigger the error before polygons analysis
     
     if ( HasValue("out") )
       {
       // Store filename extension
       // Check that the right extension is given : expected .xml
       const std::string extension = itksys::SystemTools::GetFilenameLastExtension(this->GetParameterString("out"));

       if (itksys::SystemTools::LowerCase(extension) != ".xml")
         {
         otbAppLogFATAL( << extension << " is a wrong extension for parameter \"out\": Expected .xml" );
         }
       }
174 175
  }

176
  void DoExecute() override
177 178 179
  {
  otb::ogr::DataSource::Pointer vectors = 
    otb::ogr::DataSource::New(this->GetParameterString("vec"));
180 181 182 183 184 185 186 187 188

  // Retrieve the field name
  std::vector<int> selectedCFieldIdx = GetSelectedItems("field");

  if(selectedCFieldIdx.empty())
    {
    otbAppLogFATAL(<<"No field has been selected for data labelling!");
    }

Julien Michel's avatar
Julien Michel committed
189
  std::vector<std::string> cFieldNames = GetChoiceNames("field");  
190
  std::string fieldName = cFieldNames[selectedCFieldIdx.front()];
191

192 193
  otb::Wrapper::ElevationParametersHandler::SetupDEMHandlerFromElevationParameters(this,"elev");

194 195 196 197 198 199 200 201 202 203 204 205
  // Reproject geometries
  FloatVectorImageType::Pointer inputImg = this->GetParameterImage("in");
  std::string imageProjectionRef = inputImg->GetProjectionRef();
  FloatVectorImageType::ImageKeywordlistType imageKwl =
    inputImg->GetImageKeywordlist();
  std::string vectorProjectionRef =
    vectors->GetLayer(GetParameterInt("layer")).GetProjectionRef();

  otb::ogr::DataSource::Pointer reprojVector = vectors;
  GeometriesType::Pointer inputGeomSet;
  ProjectionFilterType::Pointer geometriesProjFilter;
  GeometriesType::Pointer outputGeomSet;
206 207 208 209
  const OGRSpatialReference imgOGRSref = 
        OGRSpatialReference( imageProjectionRef.c_str() );
    const OGRSpatialReference vectorOGRSref = 
        OGRSpatialReference( vectorProjectionRef.c_str() );
210 211
  bool doReproj = true;
  // don't reproject for these cases
212 213 214
  if (  vectorProjectionRef.empty() 
     || ( imgOGRSref.IsSame( &vectorOGRSref ) ) 
     || ( imageProjectionRef.empty() && imageKwl.GetSize() == 0) )
215 216 217 218 219 220 221
    doReproj = false;

  if (doReproj)
    {
    inputGeomSet = GeometriesType::New(vectors);
    reprojVector = otb::ogr::DataSource::New();
    outputGeomSet = GeometriesType::New(reprojVector);
222
    // Filter instantiation
223 224 225 226 227 228 229 230 231 232 233 234
    geometriesProjFilter = ProjectionFilterType::New();
    geometriesProjFilter->SetInput(inputGeomSet);
    if (imageProjectionRef.empty())
      {
      geometriesProjFilter->SetOutputKeywordList(inputImg->GetImageKeywordlist()); // nec qd capteur
      }
    geometriesProjFilter->SetOutputProjectionRef(imageProjectionRef);
    geometriesProjFilter->SetOutput(outputGeomSet);
    otbAppLogINFO("Reprojecting input vectors...");
    geometriesProjFilter->Update();
    }

235 236 237 238 239 240
  FilterType::Pointer filter = FilterType::New();
  filter->SetInput(this->GetParameterImage("in"));
  if (IsParameterEnabled("mask") && HasValue("mask"))
    {
    filter->SetMask(this->GetParameterImage<UInt8ImageType>("mask"));
    }
241
  filter->SetOGRData(reprojVector);
242 243
  filter->SetFieldName(fieldName);
  filter->SetLayerIndex(this->GetParameterInt("layer"));
244
  filter->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
245

246
  AddProcess(filter->GetStreamer(),"Analyze polygons...");
247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264
  filter->Update();
  
  FilterType::ClassCountMapType &classCount = filter->GetClassCountOutput()->Get();
  FilterType::PolygonSizeMapType &polySize = filter->GetPolygonSizeOutput()->Get();
  
  StatWriterType::Pointer statWriter = StatWriterType::New();
  statWriter->SetFileName(this->GetParameterString("out"));
  statWriter->AddInputMap<FilterType::ClassCountMapType>("samplesPerClass",classCount);
  statWriter->AddInputMap<FilterType::PolygonSizeMapType>("samplesPerVector",polySize);
  statWriter->Update();
  }

};

} // end of namespace Wrapper
} // end of namespace otb

OTB_APPLICATION_EXPORT(otb::Wrapper::PolygonClassStatistics)