otbPolygonClassStatistics.cxx 8.23 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
/*=========================================================================

 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 "otbOGRDataToClassStatisticsFilter.h"
#include "otbStatisticsXMLFileWriter.h"
23 24
#include "otbGeometriesProjectionFilter.h"
#include "otbGeometriesSet.h"
25 26 27 28 29 30

namespace otb
{
namespace Wrapper
{

31 32 33 34 35 36
/** Utility function to negate std::isalnum */
bool IsNotAlphaNum(char c)
  {
  return !std::isalnum(c);
  }

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
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;

56 57 58 59
  typedef otb::GeometriesSet GeometriesType;

  typedef otb::GeometriesProjectionFilter ProjectionFilterType;

60 61 62 63 64 65
private:
  PolygonClassStatistics()
    {
   
    }

66
  void DoInit() ITK_OVERRIDE
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103
  {
    SetName("PolygonClassStatistics");
    SetDescription("Computes statistics on a training polygon set.");

    // Documentation
    SetDocName("Polygon Class Statistics");
    SetDocLongDescription("The application processes a set of geometries "
      "intended for training (they should have a field giving the associated "
      "class). The geometries are analysed against a support image to compute "
      "statistics : \n"
      "  - number of samples per class\n"
      "  - number of samples per geometry\n"
      "An optional raster mask can be used to discard samples. Different types"
      " of geometry are supported : polygons, lines, points. The behaviour is "
      "different for each type of geometry :\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\n");
    SetDocLimitations("None");
    SetDocAuthors("OTB-Team");
    SetDocSeeAlso(" ");

    AddDocTag(Tags::Learning);

    AddParameter(ParameterType_InputImage,  "in",   "InputImage");
    SetParameterDescription("in", "Support image that will be classified");

    AddParameter(ParameterType_InputImage,  "mask",   "InputMask");
    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");
    SetParameterDescription("vec","Input geometries to analyse");
    
    AddParameter(ParameterType_OutputFilename, "out", "Output Statistics");
    SetParameterDescription("out","Output file to store statistics (XML format)");

104
    AddParameter(ParameterType_ListView, "field", "Field Name");
105
    SetParameterDescription("field","Name of the field carrying the class name in the input vectors.");
106
    SetListViewSingleSelectionMode("field",true);
107 108 109 110 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);
    
    AddRAMParameter();

    // Doc example parameter settings
116 117
    SetDocExampleParameterValue("in", "support_image.tif");
    SetDocExampleParameterValue("vec", "variousVectors.sqlite");
118
    SetDocExampleParameterValue("field", "label");
119
    SetDocExampleParameterValue("out","polygonStat.xml");
120 121
  }

122
  void DoUpdateParameters() ITK_OVERRIDE
123
  {
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142
     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();
        
143
        if(fieldType == OFTInteger || ogr::version_proxy::IsOFTInteger64(fieldType))
144 145 146 147 148 149
          {
          std::string tmpKey="field."+key.substr(0, end - key.begin());
          AddChoice(tmpKey,item);
          }
        }
      }
150 151
  }

152
  void DoExecute() ITK_OVERRIDE
153 154 155
  {
  otb::ogr::DataSource::Pointer vectors = 
    otb::ogr::DataSource::New(this->GetParameterString("vec"));
156 157 158 159 160 161 162 163 164 165 166

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

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

  std::vector<std::string> cFieldNames = GetChoiceNames("cfield");  
  std::string fieldName = cFieldNames[selectedCFieldIdx.front()];
167

168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
  // 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;
  bool doReproj = true;
  // don't reproject for these cases
  if (vectorProjectionRef.empty() ||
      (imageProjectionRef == vectorProjectionRef) ||
      (imageProjectionRef.empty() && imageKwl.GetSize() == 0))
    doReproj = false;

  if (doReproj)
    {
    inputGeomSet = GeometriesType::New(vectors);
    reprojVector = otb::ogr::DataSource::New();
    outputGeomSet = GeometriesType::New(reprojVector);
192
    // Filter instantiation
193 194 195 196 197 198 199 200 201 202 203 204
    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();
    }

205 206 207 208 209 210
  FilterType::Pointer filter = FilterType::New();
  filter->SetInput(this->GetParameterImage("in"));
  if (IsParameterEnabled("mask") && HasValue("mask"))
    {
    filter->SetMask(this->GetParameterImage<UInt8ImageType>("mask"));
    }
211
  filter->SetOGRData(reprojVector);
212 213
  filter->SetFieldName(fieldName);
  filter->SetLayerIndex(this->GetParameterInt("layer"));
214
  filter->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
215 216

  AddProcess(filter->GetStreamer(),"Analyse polygons...");
217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234
  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)