otbOGRLayerClassifier.cxx 8.66 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
/*=========================================================================
 Program:   ORFEO Toolbox
 Language:  C++


 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"
27 28

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

#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:
52
  void DoInit() ITK_OVERRIDE
53 54 55 56 57 58 59 60 61 62 63
  {
    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);
  
64
    AddParameter(ParameterType_InputVectorData, "inshp", "Name of the input shapefile");
65 66
    SetParameterDescription("inshp","Name of the input shapefile");

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

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

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

88 89
  }

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

      otb::ogr::DataSource::Pointer ogrDS;
97
      otb::ogr::Layer layer(ITK_NULLPTR, false);
98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
      
      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);
        }
118 119 120
      }
  }

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

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

    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
STYLE  
OTB Bot committed
148 149
    otb::ogr::DataSource::Pointer source = otb::ogr::DataSource::New(shapefile, otb::ogr::DataSource::Modes::Read);
    otb::ogr::Layer layer = source->GetLayer(0);
150 151 152 153 154 155
    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();
156
    input->SetMeasurementVectorSize(nbFeatures);
157 158 159

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

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

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

222 223 224 225 226
    source2->SyncToDisk();

    clock_t toc = clock();

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

};
}
}

OTB_APPLICATION_EXPORT(otb::Wrapper::OGRLayerClassifier)