otbOGRLayerClassifier.cxx 8.71 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 23 24 25 26 27 28
/*=========================================================================
 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 "otbOGRDataSourceWrapper.h"
#include "otbOGRFeatureWrapper.h"
#include "otbStatisticsXMLFileWriter.h"

#include "itkVariableLengthVector.h"
#include "otbStatisticsXMLFileReader.h"

#include "itkListSample.h"
#include "otbShiftScaleSampleListFilter.h"
29 30

#ifdef OTB_USE_LIBSVM
31
#include "otbLibSVMMachineLearningModel.h"
32
#endif
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53

#include <time.h>

namespace otb
{
namespace Wrapper
{
class OGRLayerClassifier : public Application
{
public:
  typedef OGRLayerClassifier Self;
  typedef Application Superclass;
  typedef itk::SmartPointer<Self> Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
  itkNewMacro(Self)
;

  itkTypeMacro(OGRLayerClassifier, otb::Application)
;

private:
54
  void DoInit() ITK_OVERRIDE
55 56 57 58 59 60 61 62 63 64 65
  {
    SetName("OGRLayerClassifier");
    SetDescription("Classify an OGR layer based on a machine learning model and a list of features to consider.");

    SetDocName("OGRLayerClassifier");
    SetDocLongDescription("This application will apply a trained machine learning model on the selected feature to get a classification of each geometry contained in an OGR layer. The list of feature must match the list used for training. The predicted label is written in the user defined field for each geometry.");
    SetDocLimitations("Experimental. Only shapefiles are supported for now.");
    SetDocAuthors("David Youssefi during internship at CNES");
    SetDocSeeAlso("ComputeOGRLayersFeaturesStatistics,TrainOGRLayersClassifier");
    AddDocTag(Tags::Segmentation);
  
66
    AddParameter(ParameterType_InputVectorData, "inshp", "Name of the input shapefile");
67 68
    SetParameterDescription("inshp","Name of the input shapefile");

69 70
    AddParameter(ParameterType_InputFilename, "instats", "XML file containing mean and variance of each feature.");
    SetParameterDescription("instats", "XML file containing mean and variance of each feature.");
71 72 73 74 75 76 77 78 79 80

    AddParameter(ParameterType_OutputFilename, "insvm", "Input model filename.");
    SetParameterDescription("insvm", "Input model filename.");


    AddParameter(ParameterType_ListView,  "feat", "Features");
    SetParameterDescription("feat","Features to be calculated");

    AddParameter(ParameterType_String,"cfield","Field containing the predicted class.");
    SetParameterDescription("cfield","Field containing the predicted class");
81
    SetParameterString("cfield","predicted", false);
82

83 84 85 86 87 88 89
    // Doc example parameter settings
    SetDocExampleParameterValue("inshp", "vectorData.shp");
    SetDocExampleParameterValue("instats", "meanVar.xml");
    SetDocExampleParameterValue("insvm", "svmModel.svm");
    SetDocExampleParameterValue("feat", "perimeter");
    SetDocExampleParameterValue("cfield", "predicted");

90 91
  }

92
  void DoUpdateParameters() ITK_OVERRIDE
93 94
  {
    if ( HasValue("inshp") )
OTB Bot's avatar
OTB Bot committed
95
      {
96 97 98
      std::string shapefile = GetParameterString("inshp");

      otb::ogr::DataSource::Pointer ogrDS;
99
      otb::ogr::Layer layer(ITK_NULLPTR, false);
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
      
      OGRSpatialReference oSRS("");
      std::vector<std::string> options;
      
      ogrDS = otb::ogr::DataSource::New(shapefile, otb::ogr::DataSource::Modes::Read);
      std::string layername = itksys::SystemTools::GetFilenameName(shapefile);
      layername = layername.substr(0,layername.size()-4);
      layer = ogrDS->GetLayer(0);
      
      otb::ogr::Feature feature = layer.ogr().GetNextFeature();
      ClearChoices("feat");
      for(int iField=0; iField<feature.ogr().GetFieldCount(); iField++)
        {
        std::string key, item = feature.ogr().GetFieldDefnRef(iField)->GetNameRef();
        key = item;
        key.erase(std::remove(key.begin(), key.end(), ' '), key.end());
        std::transform(key.begin(), key.end(), key.begin(), tolower);
        key="feat."+key;
        AddChoice(key,item);
        }
120 121 122
      }
  }

123
  void DoExecute() ITK_OVERRIDE
124
  {
125 126
      
    #ifdef OTB_USE_LIBSVM 
127 128
    clock_t tic = clock();

129 130 131
    std::string shapefile = GetParameterString("inshp").c_str();
    std::string XMLfile = GetParameterString("instats").c_str();
    std::string modelfile = GetParameterString("insvm").c_str();
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149

    typedef double ValueType;
    typedef itk::VariableLengthVector<ValueType> MeasurementType;
    typedef itk::Statistics::ListSample <MeasurementType> ListSampleType;
    typedef otb::StatisticsXMLFileReader<MeasurementType> StatisticsReader;
  
    typedef unsigned int LabelPixelType;
    typedef itk::FixedArray<LabelPixelType,1> LabelSampleType;
    typedef itk::Statistics::ListSample <LabelSampleType> LabelListSampleType;
  
    typedef otb::Statistics::ShiftScaleSampleListFilter<ListSampleType, ListSampleType> ShiftScaleFilterType;
  
    StatisticsReader::Pointer statisticsReader = StatisticsReader::New();
    statisticsReader->SetFileName(XMLfile);

    MeasurementType meanMeasurementVector = statisticsReader->GetStatisticVectorByName("mean");
    MeasurementType stddevMeasurementVector = statisticsReader->GetStatisticVectorByName("stddev");

OTB Bot's avatar
OTB Bot committed
150 151
    otb::ogr::DataSource::Pointer source = otb::ogr::DataSource::New(shapefile, otb::ogr::DataSource::Modes::Read);
    otb::ogr::Layer layer = source->GetLayer(0);
152 153 154 155 156 157
    bool goesOn = true;
    otb::ogr::Feature feature = layer.ogr().GetNextFeature();

    ListSampleType::Pointer input = ListSampleType::New();
    LabelListSampleType::Pointer target = LabelListSampleType::New();
    const int nbFeatures = GetSelectedItems("feat").size();
158
    input->SetMeasurementVectorSize(nbFeatures);
159 160 161

    if(feature.addr())
      while(goesOn)
OTB Bot's avatar
OTB Bot committed
162 163 164 165 166 167 168 169 170
       {
         MeasurementType mv; mv.SetSize(nbFeatures);
         
         for(int idx=0; idx < nbFeatures; ++idx)
           mv[idx] = feature.ogr().GetFieldAsDouble(GetSelectedItems("feat")[idx]);
         
         input->PushBack(mv);
         target->PushBack(feature.ogr().GetFieldAsInteger("class"));
         feature = layer.ogr().GetNextFeature();
171
         goesOn = feature.addr() != ITK_NULLPTR;
OTB Bot's avatar
OTB Bot committed
172
       }
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195

    ShiftScaleFilterType::Pointer trainingShiftScaleFilter = ShiftScaleFilterType::New();
    trainingShiftScaleFilter->SetInput(input);
    trainingShiftScaleFilter->SetShifts(meanMeasurementVector);
    trainingShiftScaleFilter->SetScales(stddevMeasurementVector);
    trainingShiftScaleFilter->Update();
  
    ListSampleType::Pointer listSample;
    LabelListSampleType::Pointer labelListSample;
  
    listSample = trainingShiftScaleFilter->GetOutput();
    labelListSample = target;

    ListSampleType::Pointer trainingListSample = listSample;
    LabelListSampleType::Pointer trainingLabeledListSample = labelListSample;

    typedef otb::LibSVMMachineLearningModel<ValueType,LabelPixelType> LibSVMType;
    LibSVMType::Pointer libSVMClassifier = LibSVMType::New();
    libSVMClassifier->SetInputListSample(trainingListSample);
    libSVMClassifier->SetTargetListSample(trainingLabeledListSample);
    libSVMClassifier->Load(modelfile);
    libSVMClassifier->PredictAll();

OTB Bot's avatar
OTB Bot committed
196 197
    otb::ogr::DataSource::Pointer source2 = otb::ogr::DataSource::New(shapefile, otb::ogr::DataSource::Modes::Update_LayerUpdate);
    otb::ogr::Layer layer2 = source2->GetLayer(0);
198 199 200 201 202 203 204 205 206 207 208
 
    OGRFieldDefn predictedField(GetParameterString("cfield").c_str(), OFTInteger);
    layer2.CreateField(predictedField, true);

    bool goesOn2 = true;
    layer2.ogr().ResetReading();
    otb::ogr::Feature feature2 = layer2.ogr().GetNextFeature();
    unsigned int count=0;

    if(feature2.addr())
      while(goesOn2)
OTB Bot's avatar
OTB Bot committed
209
       {
210
        feature2.ogr().SetField(GetParameterString("cfield").c_str(),(int)labelListSample->GetMeasurementVector(count)[0]);
OTB Bot's avatar
OTB Bot committed
211 212
         layer2.SetFeature(feature2);
         feature2 = layer2.ogr().GetNextFeature();
213
         goesOn2 = feature2.addr() != ITK_NULLPTR;
OTB Bot's avatar
OTB Bot committed
214 215
         count++;
       }
216 217 218 219 220 221 222 223
    
    const OGRErr err = layer2.ogr().CommitTransaction();

    if (err != OGRERR_NONE)
      {
      itkExceptionMacro(<< "Unable to commit transaction for OGR layer " << layer2.ogr().GetName() << ".");
      }

224 225 226 227 228
    source2->SyncToDisk();

    clock_t toc = clock();

    otbAppLogINFO( "Elapsed: "<< ((double)(toc - tic) / CLOCKS_PER_SEC)<<" seconds.");
229
        
230
    #else
231 232
    otbAppLogFATAL("Module LIBSVM is not installed. You should consider turning OTB_USE_LIBSVM on during cmake configuration.");
    #endif
233
    
234 235 236 237 238 239 240 241 242
  }

};
}
}

OTB_APPLICATION_EXPORT(otb::Wrapper::OGRLayerClassifier)