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

  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 <iostream>

21
22
23
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbWrapperStringListParameter.h"
24
25
#include "otbImageToEnvelopeVectorDataFilter.h"
#include "otbVectorDataToRandomLineGenerator.h"
26
27
28
29
#include "itkAmoebaOptimizer.h"
#include "otbStandardDSCostFunction.h"


30
31
32
33
34
35
namespace otb
{

namespace Wrapper
{

36
37
38
39
#include "itkCommand.h"
class CommandIterationUpdate : public itk::Command
{
public:
40
41
42
43
typedef  CommandIterationUpdate   Self;
typedef  itk::Command             Superclass;
typedef itk::SmartPointer<Self>   Pointer;
itkNewMacro( Self );
44
protected:
45
CommandIterationUpdate() {};
46
public:
47
48
typedef itk::AmoebaOptimizer         OptimizerType;
typedef   const OptimizerType   *    OptimizerPointer;
49

50

51
void Execute(itk::Object *caller, const itk::EventObject & event) ITK_OVERRIDE
52
53
54
{
  Execute( (const itk::Object *)caller, event);
}
55

56
void Execute(const itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE
57
58
59
60
{
  OptimizerPointer optimizer =
      dynamic_cast< OptimizerPointer >( object );
  if( ! itk::IterationEvent().CheckEvent( &event ) )
61
    {
62
    return;
63
    }
64
65
66
67
  std::ostringstream message;
  message << optimizer->GetCachedValue() << "   ";
  message << optimizer->GetCachedCurrentPosition() << std::endl;
  std::cout<<message.str()<<std::endl;
68
}
69
70


71
72
};

73

74
class DSFuzzyModelEstimation: public Application
75
{
76
77
78
79
80
81
public:
  /** Standard class typedefs. */
  typedef DSFuzzyModelEstimation Self;
  typedef Application Superclass;
  typedef itk::SmartPointer<Self> Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
82

83
84
85
86
  typedef VectorData<double>                          VectorDataType;

  typedef VectorDataType::DataTreeType                DataTreeType;
  typedef VectorDataType::DataNodeType                DataNodeType;
87
88
89

  typedef VectorDataType::ValuePrecisionType          PrecisionType;
  typedef VectorDataType::PrecisionType               CoordRepType;
90
91

  typedef otb::Wrapper::StringListParameter::StringListType         StringListType;
92
93

  typedef otb::VectorDataToDSValidatedVectorDataFilter<VectorDataType, PrecisionType>
94
                                                            ValidationFilterType;
Manuel Grizonnet's avatar
Manuel Grizonnet committed
95
  typedef otb::StandardDSCostFunction<ValidationFilterType> CostFunctionType;
96
  typedef CostFunctionType::LabelSetType                    LabelSetType;
97

98
  typedef itk::AmoebaOptimizer                            OptimizerType;
99

100
  typedef otb::FuzzyDescriptorsModelManager::DescriptorsModelType
101
                                                          DescriptorsModelType;
102
  typedef otb::FuzzyDescriptorsModelManager::DescriptorListType
103
                                                          DescriptorListType;
104

105
106
     typedef itk::PreOrderTreeIterator<VectorDataType::DataTreeType>
     TreeIteratorType;
107
108


109
  /** Standard macro */
OTB Bot's avatar
STYLE    
OTB Bot committed
110
  itkNewMacro(Self);
111

OTB Bot's avatar
STYLE    
OTB Bot committed
112
  itkTypeMacro(DSFuzzyModelEstimation, otb::Application);
113
114

private:
115
  void DoInit() ITK_OVERRIDE
116
117
118
119
  {
    SetName("DSFuzzyModelEstimation");
    SetDescription("Estimate feature fuzzy model parameters using 2 vector data (ground truth samples and wrong samples).");

120
    SetDocName("Fuzzy Model estimation");
Jonathan Guinet's avatar
Jonathan Guinet committed
121
122
    SetDocLongDescription("Estimate feature fuzzy model parameters using 2 vector data (ground truth samples and wrong samples).");
    SetDocLimitations("None.");
123
124
    SetDocAuthors("OTB-Team");
    SetDocSeeAlso(" ");
125

126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
    AddDocTag(Tags::FeatureExtraction);

    AddParameter(ParameterType_InputVectorData, "psin", "Input Positive Vector Data");
    SetParameterDescription("psin", "Ground truth vector data for positive samples");

    AddParameter(ParameterType_InputVectorData, "nsin", "Input Negative Vector Data");
    SetParameterDescription("nsin", "Ground truth vector data for negative samples");

    AddParameter(ParameterType_StringList, "belsup", "Belief Support");
    SetParameterDescription("belsup", "Dempster Shafer study hypothesis to compute belief");

    AddParameter(ParameterType_StringList, "plasup", "Plausibility Support");
    SetParameterDescription("plasup", "Dempster Shafer study hypothesis to compute plausibility");

    AddParameter(ParameterType_String, "cri", "Criterion");
    SetParameterDescription("cri", "Dempster Shafer criterion (by default (belief+plausibility)/2)");
    MandatoryOff("cri");
143
    SetParameterString("cri","((Belief + Plausibility)/2.)", false);
144
145
146
147

    AddParameter(ParameterType_Float,"wgt","Weighting");
    SetParameterDescription("wgt","Coefficient between 0 and 1 to promote undetection or false detections (default 0.5)");
    MandatoryOff("wgt");
148
    SetParameterFloat("wgt",0.5, false);
149

150
    AddParameter(ParameterType_InputFilename,"initmod","initialization model");
151
    SetParameterDescription("initmod","Initialization model (xml file) to be used. If the xml initialization model is set, the descriptor list is not used (specified using the option -desclist)");
152
153
154
    MandatoryOff("initmod");

    AddParameter(ParameterType_StringList, "desclist","Descriptor list");
155
    SetParameterDescription("desclist","List of the descriptors to be used in the model (must be specified to perform an automatic initialization)");
156
    MandatoryOff("desclist");
157
    SetParameterString("desclist","", false);
158
159
160
161

    AddParameter(ParameterType_Int,"maxnbit","Maximum number of iterations");
    MandatoryOff("maxnbit");
    SetParameterDescription("maxnbit","Maximum number of optimizer iteration (default 200)");
162
    SetParameterInt("maxnbit",200, false);
163
164
165
166
167

    AddParameter(ParameterType_Empty,"optobs","Optimizer Observer");
    SetParameterDescription("optobs","Activate the optimizer observer");
    MandatoryOff("optobs");

168
    AddParameter(ParameterType_OutputFilename,"out","Output filename");
169
    SetParameterDescription("out","Output model file name (xml file) contains the optimal model to perform information fusion.");
170

171
172
173
174
175
176
177
178
179
    // Doc example parameter settings
    SetDocExampleParameterValue("psin", "cdbTvComputePolylineFeatureFromImage_LI_NOBUIL_gt.shp");
    SetDocExampleParameterValue("nsin", "cdbTvComputePolylineFeatureFromImage_LI_NOBUIL_wr.shp");
    SetDocExampleParameterValue("belsup", "\"ROADSA\"");
    SetDocExampleParameterValue("plasup", "\"NONDVI\" \"ROADSA\" \"NOBUIL\"");
    SetDocExampleParameterValue("initmod", "Dempster-Shafer/DSFuzzyModel_Init.xml");
    SetDocExampleParameterValue("maxnbit", "4");
    SetDocExampleParameterValue("optobs", "true");
    SetDocExampleParameterValue("out", "DSFuzzyModelEstimation.xml");
180
181
  }

182
  void DoUpdateParameters() ITK_OVERRIDE
183
184
185
186
187
188
189
190
191
  {
    // Nothing to do here : all parameters are independent


    // .. //


  }

192
  void DoExecute() ITK_OVERRIDE
193
194
195
196
197
198
199
200
201
  {

    //Instantiate
    m_CostFunction = CostFunctionType::New();
    m_Optimizer = OptimizerType::New();

    //Read the vector datas
    VectorDataType::Pointer psVectorData = GetParameterVectorData("psin");
    psVectorData->Update();
202
    VectorDataType::Pointer nsVectorData = GetParameterVectorData("nsin");
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
    nsVectorData->Update();

    // Load the initial descriptor model
    DescriptorListType descList;
    DescriptorsModelType descMod;
    if (IsParameterEnabled("initmod"))
      {
      std::string descModFile = GetParameterString("initmod");
      descMod = FuzzyDescriptorsModelManager::Read(descModFile.c_str());
      descList = FuzzyDescriptorsModelManager::GetDescriptorList(descMod);
      }
    else
      {
      StringListType stringList = GetParameterStringList("desclist");
      int nbsdDesc = stringList.size();
      for (int i = 0; i < nbsdDesc; i++)
219
        {
220
        descList.push_back(stringList[i]);
221
        }
222
      }
223

224
    m_CostFunction->SetDescriptorList(descList);
225

226
    // Compute statistics of all the descriptors
227

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
    std::vector<double> accFirstOrderPS, accSecondOrderPS, minPS, maxPS;
    accFirstOrderPS.resize(descList.size());
    accSecondOrderPS.resize(descList.size());
    std::fill(accFirstOrderPS.begin(), accFirstOrderPS.end(), 0);
    std::fill(accSecondOrderPS.begin(), accSecondOrderPS.end(), 0);
    minPS.resize(descList.size());
    maxPS.resize(descList.size());
    unsigned int accNbElemPS = 0;

    TreeIteratorType itVectorPS(psVectorData->GetDataTree());
    for (itVectorPS.GoToBegin(); !itVectorPS.IsAtEnd(); ++itVectorPS)
      {
      if (!itVectorPS.Get()->IsRoot() && !itVectorPS.Get()->IsDocument() && !itVectorPS.Get()->IsFolder())
        {
        DataNodeType::Pointer currentGeometry = itVectorPS.Get();

        for (unsigned int i = 0; i < descList.size(); ++i)
245
          {
246
247
248
249
          double desc = currentGeometry->GetFieldAsDouble(descList[i]);

          accFirstOrderPS[i] += desc;
          accSecondOrderPS[i] += desc * desc;
250

251
252
253
254
255
256
257
258
259
260
261
          if (desc < minPS[i])
            {
            minPS[i] = desc;
            }
          if (desc > maxPS[i])
            {
            maxPS[i] = desc;
            }

          }
        accNbElemPS++;
262
263
        }
      }
Guillaume Pasero's avatar
Guillaume Pasero committed
264
265
266
267
    if (accNbElemPS == 0)
      {
      otbAppLogFATAL(<< "Error : no element found in positive vector data!");
      }
268

269
270
271
272
273
274
275
276
277
278
279
280
281
    TreeIteratorType itVectorNS(nsVectorData->GetDataTree());
    std::vector<double> accFirstOrderNS, accSecondOrderNS, minNS, maxNS;
    minNS.resize(descList.size());
    maxNS.resize(descList.size());
    accFirstOrderNS.resize(descList.size());
    accSecondOrderNS.resize(descList.size());
    std::fill(accFirstOrderNS.begin(), accFirstOrderNS.end(), 0);
    std::fill(accSecondOrderNS.begin(), accSecondOrderNS.end(), 0);
    std::fill(minNS.begin(), minNS.end(), 1);
    std::fill(maxNS.begin(), maxNS.end(), 0);
    unsigned int accNbElemNS = 0;

    for (itVectorNS.GoToBegin(); !itVectorNS.IsAtEnd(); ++itVectorNS)
282
      {
283
      if (!itVectorNS.Get()->IsRoot() && !itVectorNS.Get()->IsDocument() && !itVectorNS.Get()->IsFolder())
284
        {
285
        DataNodeType::Pointer currentGeometry = itVectorNS.Get();
286

287
        for (unsigned int i = 0; i < descList.size(); ++i)
288
          {
289
          double desc = currentGeometry->GetFieldAsDouble(descList[i]);
290

291
292
293
294
295
296
297
298
299
300
301
302
303
304
          accFirstOrderNS[i] += desc;
          accSecondOrderNS[i] += desc * desc;

          if (desc < minNS[i])
            {
            minNS[i] = desc;
            }
          if (desc > maxNS[i])
            {
            maxNS[i] = desc;
            }

          }
        accNbElemNS++;
305
306
        }
      }
Guillaume Pasero's avatar
Guillaume Pasero committed
307
308
309
310
    if (accNbElemNS == 0)
      {
      otbAppLogFATAL(<< "Error : no element found in negative vector data!");
      }
311
312
313
314
315
316
317
318
    otbAppLogINFO( << "Descriptors Stats : ");
    otbAppLogINFO( << "Positive Samples");
    for (unsigned int i = 0; i < descList.size(); ++i)
      {
      double mean = accFirstOrderPS[i] / accNbElemPS;
      double stddev = vcl_sqrt(accSecondOrderPS[i] / accNbElemPS - mean * mean);
      otbAppLogINFO( << descList[i] << "  :  " << mean << " +/- " << stddev << "  (min: " << minPS[i] << "  max: " << maxPS[i] << ")"<< std::endl);
      }
319

320
321
322
323
324
325
326
    otbAppLogINFO( << "Negative Samples" << std::endl);
    for (unsigned int i = 0; i < descList.size(); ++i)
      {
      double mean = accFirstOrderNS[i] / accNbElemNS;
      double stddev = vcl_sqrt(accSecondOrderNS[i] / accNbElemNS - mean * mean);
      otbAppLogINFO(<< descList[i] << "  :  " << mean << " +/- " << stddev << "  (min: " << minNS[i] << "  max: " << maxNS[i] << ")"<< std::endl);
      }
327

328
    OptimizerType::ParametersType initialPosition(4 * descList.size());
329

330
331
    if (IsParameterEnabled("initmod"))
      {
332

333
334
335
336
337
338
339
340
341
342
343
      for (unsigned int i = 0; i < 4; i++)
        {
        for (unsigned int j = 0; j < descList.size(); j++)
          {
          initialPosition.SetElement(
                                     i + 4 * j,
                                     otb::FuzzyDescriptorsModelManager::GetDescriptor(descList[j].c_str(), descMod).second[i]);
          }
        }
      }
    else
344
      {
345
      for (unsigned int j = 0; j < descList.size(); j++)
346
        {
347
348
349
350
351
352
353
        initialPosition.SetElement((j * 4), std::min(minNS[j], maxPS[j]));
        initialPosition.SetElement((j * 4) + 2, std::max(minNS[j], maxPS[j]));
        initialPosition.SetElement(
                                   (j * 4) + 1,
                                       0.5
                                       * (initialPosition.GetElement((j * 4)) + initialPosition.GetElement((j * 4) + 2)));
        initialPosition.SetElement((j * 4) + 3, 0.95);
354
355
        }
      }
356

Julien Michel's avatar
Julien Michel committed
357
358
    otbAppLogINFO(<<"Initial model: "<<initialPosition);

359
360
361
362
363
364
365
366
367
    //Cost Function
    //Format Hypothesis
    LabelSetType Bhyp, Phyp;
    int nbSet;

    StringListType stringList = GetParameterStringList("belsup");
    nbSet = stringList.size();

    for (int i = 0; i < nbSet; i++)
368
      {
369
370
      std::string str = stringList[i];
      Bhyp.insert(str);
371
      }
372
373
374
375
376
377
378
379
380
    m_CostFunction->SetBeliefHypothesis(Bhyp);
    stringList = GetParameterStringList("plasup");
    nbSet = stringList.size();
    for (int i = 0; i < nbSet; i++)
      {
      std::string str = stringList[i];
      Phyp.insert(str);
      }
    m_CostFunction->SetPlausibilityHypothesis(Phyp);
381

382
    m_CostFunction->SetWeight(GetParameterFloat("wgt"));
383

384
    m_CostFunction->SetCriterionFormula(GetParameterString("cri"));
385

386
387
388
389
    m_CostFunction->SetGTVectorData(psVectorData);
    m_CostFunction->SetNSVectorData(nsVectorData);
    //Optimizer
    m_Optimizer->SetCostFunction(m_CostFunction);
390

391
    m_Optimizer->SetMaximumNumberOfIterations(GetParameterInt("maxnbit"));
392

393
394
    OptimizerType::ParametersType simplexDelta(m_CostFunction->GetNumberOfParameters());
    simplexDelta.Fill(0.1);
395

396
397
    m_Optimizer->AutomaticInitialSimplexOff();
    m_Optimizer->SetInitialSimplexDelta(simplexDelta);
398

399
    m_Optimizer->SetInitialPosition(initialPosition);
400

401
402
403
    // Create the Command observer and register it with the optimizer.
    CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
    if (IsParameterEnabled("optobs"))
404
      {
405
      m_Optimizer->AddObserver(itk::IterationEvent(), observer);
406
      }
407

408
409
410
411
412
413
414
415
416
    try
      {
      // do the optimization
      m_Optimizer->StartOptimization();
      }
    catch (itk::ExceptionObject& err)
      {
      // An error has occurred in the optimization.
      // Update the parameters
417
      otbAppLogFATAL("ERROR: Exception Caught : "<< err.GetDescription() << std::endl
Guillaume Pasero's avatar
Guillaume Pasero committed
418
419
        << "numberOfIterations : " << m_Optimizer->GetOptimizer()->get_num_evaluations() << std::endl
        << "Results : " << m_Optimizer->GetCurrentPosition() << std::endl);
420
421
422
      }
    // get the results
    const unsigned int numberOfIterations = m_Optimizer->GetOptimizer()->get_num_evaluations();
423
424
    otbAppLogINFO("numberOfIterations : " << numberOfIterations << std::endl);
    otbAppLogINFO("Results : " << m_Optimizer->GetCurrentPosition() << std::endl);
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444

    for (unsigned int i = 0; i < descList.size(); i++)
      {
      otb::FuzzyDescriptorsModelManager::ParameterType tmpParams;
      for (unsigned int j = 0; j < 4; j++)
        {
        tmpParams.push_back(m_Optimizer->GetCurrentPosition()[(i * 4) + j]);
        }
      otb::FuzzyDescriptorsModelManager::AddDescriptor(descList[i], tmpParams, m_Model);
      }
    otb::FuzzyDescriptorsModelManager::Save(GetParameterString("out"), m_Model);

  };

  CostFunctionType::Pointer                               m_CostFunction;
  OptimizerType::Pointer                                  m_Optimizer;
  otb::FuzzyDescriptorsModelManager::DescriptorsModelType m_Model;
};

}
445
}
446
447
448

OTB_APPLICATION_EXPORT(otb::Wrapper::DSFuzzyModelEstimation)