otbOGRLayerClassifier.cxx 8.78 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 26 27 28 29 30 31 32

#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"
33 34

#ifdef OTB_USE_LIBSVM
35
#include "otbLibSVMMachineLearningModel.h"
36
#endif
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

#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:
58
  void DoInit() ITK_OVERRIDE
59 60 61 62 63 64 65 66
  {
    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");
67
    SetDocSeeAlso("ComputeOGRLayersFeaturesStatistics");
68 69
    AddDocTag(Tags::Segmentation);
  
70
    AddParameter(ParameterType_InputVectorData, "inshp", "Name of the input shapefile");
71 72
    SetParameterDescription("inshp","Name of the input shapefile");

73 74
    AddParameter(ParameterType_InputFilename, "instats", "XML file containing mean and variance of each feature.");
    SetParameterDescription("instats", "XML file containing mean and variance of each feature.");
75 76 77 78 79 80 81 82 83 84

    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");
85
    SetParameterString("cfield","predicted", false);
86

87 88 89 90 91 92 93
    // Doc example parameter settings
    SetDocExampleParameterValue("inshp", "vectorData.shp");
    SetDocExampleParameterValue("instats", "meanVar.xml");
    SetDocExampleParameterValue("insvm", "svmModel.svm");
    SetDocExampleParameterValue("feat", "perimeter");
    SetDocExampleParameterValue("cfield", "predicted");

94
    otbAppDocLink();
95 96
  }

97
  void DoUpdateParameters() ITK_OVERRIDE
98 99
  {
    if ( HasValue("inshp") )
OTB Bot's avatar
OTB Bot committed
100
      {
101 102 103
      std::string shapefile = GetParameterString("inshp");

      otb::ogr::DataSource::Pointer ogrDS;
104
      otb::ogr::Layer layer(ITK_NULLPTR, false);
105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124
      
      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);
        }
125 126 127
      }
  }

128
  void DoExecute() ITK_OVERRIDE
129
  {
130 131
      
    #ifdef OTB_USE_LIBSVM 
132 133
    clock_t tic = clock();

134 135 136
    std::string shapefile = GetParameterString("inshp").c_str();
    std::string XMLfile = GetParameterString("instats").c_str();
    std::string modelfile = GetParameterString("insvm").c_str();
137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154

    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
155 156
    otb::ogr::DataSource::Pointer source = otb::ogr::DataSource::New(shapefile, otb::ogr::DataSource::Modes::Read);
    otb::ogr::Layer layer = source->GetLayer(0);
157 158 159 160 161 162
    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();
163
    input->SetMeasurementVectorSize(nbFeatures);
164 165 166

    if(feature.addr())
      while(goesOn)
OTB Bot's avatar
OTB Bot committed
167 168 169 170 171 172 173 174 175
       {
         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();
176
         goesOn = feature.addr() != ITK_NULLPTR;
OTB Bot's avatar
OTB Bot committed
177
       }
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196

    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->Load(modelfile);
197
    trainingLabeledListSample = libSVMClassifier->PredictBatch(trainingListSample);
198

OTB Bot's avatar
OTB Bot committed
199 200
    otb::ogr::DataSource::Pointer source2 = otb::ogr::DataSource::New(shapefile, otb::ogr::DataSource::Modes::Update_LayerUpdate);
    otb::ogr::Layer layer2 = source2->GetLayer(0);
201 202 203 204 205 206 207 208 209 210 211
 
    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
212
       {
213
        feature2.ogr().SetField(GetParameterString("cfield").c_str(),(int)labelListSample->GetMeasurementVector(count)[0]);
OTB Bot's avatar
OTB Bot committed
214 215
         layer2.SetFeature(feature2);
         feature2 = layer2.ogr().GetNextFeature();
216
         goesOn2 = feature2.addr() != ITK_NULLPTR;
OTB Bot's avatar
OTB Bot committed
217 218
         count++;
       }
219 220 221 222 223 224 225 226
    
    const OGRErr err = layer2.ogr().CommitTransaction();

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

227 228 229 230 231
    source2->SyncToDisk();

    clock_t toc = clock();

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

};
}
}

OTB_APPLICATION_EXPORT(otb::Wrapper::OGRLayerClassifier)