otbSegmentation.cxx 35.6 KB
Newer Older
1
2
/*=========================================================================

3
4
5
6
  Program:   ORFEO Toolbox
  Language:  C++
  Date:      $Date$
  Version:   $Revision$
7
8


9
10
  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  See OTBCopyright.txt for details.
11
12


13
14
15
  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.
16

17
  =========================================================================*/
18
// Wrappers
19
20
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
21
#include "otbWrapperChoiceParameter.h"
22

23
// Segmentation filters includes
24
#include "otbMeanShiftSegmentationFilter.h"
25
#include "otbConnectedComponentMuParserFunctor.h"
26
#include "otbMaskMuParserFilter.h"
27
28
29
#include "otbVectorImageToAmplitudeImageFilter.h"
#include "itkGradientMagnitudeImageFilter.h"
#include "otbWatershedSegmentationFilter.h"
30
#include "otbMorphologicalProfilesSegmentationFilter.h"
31

32
// Large scale vectorization framework
33
#include "otbStreamingImageToOGRLayerSegmentationFilter.h"
34

35
// Fusion filter
36
#include "otbOGRLayerStreamStitchingFilter.h"
37

38
#include "otbGeoInformationConversion.h"
39
40
41
42

//Utils
#include "itksys/SystemTools.hxx"

43
44
namespace otb
{
45
46
namespace Wrapper
{
47
class Segmentation : public Application
48
49
50
{
public:
  /** Standard class typedefs. */
51
  typedef Segmentation        Self;
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
  typedef Application                   Superclass;
  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

  /** Images typedefs */
  typedef UInt32ImageType               LabelImageType;
  typedef UInt32ImageType               MaskImageType;

  // Segmentation filters typedefs

  // Home made mean-shift
  typedef otb::MeanShiftSegmentationFilter
  <FloatVectorImageType,
   LabelImageType,
   FloatVectorImageType>                  MeanShiftSegmentationFilterType;

  // Simple connected components
  typedef otb::Functor::ConnectedComponentMuParserFunctor
  <FloatVectorImageType::PixelType>       FunctorType;

  typedef itk::ConnectedComponentFunctorImageFilter
  <FloatVectorImageType,
   LabelImageType,
   FunctorType,
   MaskImageType >                        ConnectedComponentSegmentationFilterType;
77

78
79
80
81
82
  typedef itk::ScalarConnectedComponentImageFilter
  <LabelImageType,
   LabelImageType>                        LabeledConnectedComponentSegmentationFilterType;


83
84
85
86
  // Watershed
  typedef otb::VectorImageToAmplitudeImageFilter
  <FloatVectorImageType,
   FloatImageType>                        AmplitudeFilterType;
87

88
89
  typedef itk::GradientMagnitudeImageFilter
  <FloatImageType,FloatImageType>         GradientMagnitudeFilterType;
90

91
92
  typedef otb::WatershedSegmentationFilter
  <FloatImageType,LabelImageType>         WatershedSegmentationFilterType;
93

94
  // Geodesic morphology multiscale segmentation
95
  typedef otb::MorphologicalProfilesSegmentationFilter<FloatImageType,LabelImageType> MorphologicalProfilesSegmentationFilterType;
96

97

98
99
100
  // mask filter
  typedef otb::MaskMuParserFilter
  <FloatVectorImageType, MaskImageType>   MaskMuParserFilterType;
101

102
  // Vectorize filters
103

104
  // Home made mean-shift
105
  typedef otb::StreamingImageToOGRLayerSegmentationFilter
106
107
  <FloatVectorImageType,
   MeanShiftSegmentationFilterType>       MeanShiftVectorizedSegmentationOGRType;
108

109

110
  // Connected components
111
  typedef otb::StreamingImageToOGRLayerSegmentationFilter
112
113
114
  <FloatVectorImageType,
   ConnectedComponentSegmentationFilterType>
  ConnectedComponentStreamingVectorizedSegmentationOGRType;
115

116
117
118
  // Morphological profiles
  typedef otb::StreamingImageToOGRLayerSegmentationFilter
  <FloatImageType,MorphologicalProfilesSegmentationFilterType> MorphoVectorizedSegmentationOGRType;
119

120
  typedef otb::OGRLayerStreamStitchingFilter
121
  <FloatVectorImageType>                  FusionFilterType;
122

123

124
  // Watershed
125
  typedef otb::StreamingImageToOGRLayerSegmentationFilter
126
127
  <FloatImageType,
   WatershedSegmentationFilterType>      StreamingVectorizedWatershedFilterType;
128

129
130
  /** Standard macro */
  itkNewMacro(Self);
131
  itkTypeMacro(Segmentation, otb::Application);
132

133
134
135
private:
  void DoInit()
  {
136
    SetName("Segmentation");
137
    SetDescription("Performs segmentation of an image, and output either a raster or a vector file. In vector mode, large input datasets are supported.");
138
139

    // Documentation
140
    SetDocName("Segmentation");
141
    SetDocLongDescription("This application allows to perform various segmentation algorithms on a multispectral image."
142
143
144
145
146
147
148
                          "Available segmentation algorithms are two different versions of Mean-Shift segmentation algorithm (one being multi-threaded),"
                          " simple pixel based connected components according to a user-defined criterion, and watershed from the gradient of the intensity"
                          " (norm of spectral bands vector). The application has two different modes that affects the nature of its output.\n\nIn raster mode,"
                          " the output of the application is a classical image of unique labels identifying the segmented regions. The labeled output can be passed to the"
                          " ColorMapping application to render regions with contrasted colours. Please note that this mode loads the whole input image into memory, and as such"
                          " can not handle large images. \n\nTo segment large data, one can use the vector mode. In this case, the output of the application is a"
                          " vector file or database. The input image is split into tiles (whose size can be set using the tilesize parameter), and each tile is loaded, segmented"
149
                          " with the chosen algorithm, vectorized, and written into the output file or database. This piece-wise behavior ensure that memory will never get overloaded,"
150
151
                          " and that images of any size can be processed. There are few more options in the vector mode. The simplify option allows to simplify the geometry"
                          " (i.e. remove nodes in polygons) according to a user-defined tolerance. The stitch option allows to application to try to stitch together polygons corresponding"
152
                          " to segmented region that may have been split by the tiling scheme. ");
153

154
    SetDocLimitations("In raster mode, the application can not handle large input images. Stitching step of vector mode might become slow with very large input images."
155
156
157
                     " \nMeanShift filter results depends on the number of threads used. \nWatershed and multiscale geodesic morphology segmentation will be performed on the amplitude "
                     " of the input image.");

158
159
160
161
162
163
    SetDocAuthors("OTB-Team");
    SetDocSeeAlso("MeanShiftSegmentation");

    AddDocTag(Tags::Segmentation);

    AddParameter(ParameterType_InputImage, "in", "Input Image");
164
    SetParameterDescription("in", "The input image to segment");
165
166

    AddParameter(ParameterType_Choice, "filter", "Segmentation algorithm");
167
    SetParameterDescription("filter", "Choice of segmentation algorithm (mean-shift by default)");
168

169
    AddChoice("filter.meanshift", "Mean-Shift");
170
    SetParameterDescription(
171
      "filter.meanshift","OTB implementation of the Mean-Shift algorithm (multi-threaded).");
172
173
    // MeanShift Parameters
    AddParameter(ParameterType_Int, "filter.meanshift.spatialr", "Spatial radius");
174
    SetParameterDescription("filter.meanshift.spatialr", "Spatial radius of the neighborhood.");
175
    AddParameter(ParameterType_Float, "filter.meanshift.ranger", "Range radius");
176
    SetParameterDescription("filter.meanshift.ranger", "Range radius defining the radius (expressed in radiometry unit) in the multispectral space.");
177
178
179
180
181
    AddParameter(ParameterType_Float, "filter.meanshift.thres", "Mode convergence threshold");
    SetParameterDescription("filter.meanshift.thres", "Algorithm iterative scheme will stop if mean-shift "
                            "vector is below this threshold or if iteration number reached maximum number of iterations.");
    AddParameter(ParameterType_Int, "filter.meanshift.maxiter", "Maximum number of iterations");
    SetParameterDescription("filter.meanshift.maxiter", "Algorithm iterative scheme will stop if convergence hasn't been reached after the maximum number of iterations.");
182
183
    AddParameter(ParameterType_Int, "filter.meanshift.minsize", "Minimum region size");
    SetParameterDescription("filter.meanshift.minsize", "Minimum size of a region (in pixel unit) in segmentation. Smaller clusters will be merged to the neighboring cluster with the closest radiometry."
184
      " If set to 0 no pruning is done.");
185

186
187
188
189
190
191
192
    //AddParameter(ParameterType_Empty, "filter.meanshift.useoptim", "use optimization");
    //SetParameterDescription("filter.meanshift.useoptim", "Use mode optimization.");
    //MandatoryOff("filter.meanshift.useoptim");

    SetDefaultParameterInt("filter.meanshift.spatialr", 5);
    SetDefaultParameterFloat("filter.meanshift.ranger", 15.0);
    SetDefaultParameterFloat("filter.meanshift.thres", 0.1);
193
194
195
    SetMinimumParameterFloatValue("filter.meanshift.thres", 0.1);
    SetDefaultParameterInt("filter.meanshift.minsize", 100);
    SetMinimumParameterIntValue("filter.meanshift.minsize", 0);
196
197
198
199
    SetDefaultParameterInt("filter.meanshift.maxiter", 100);
    SetMinimumParameterIntValue("filter.meanshift.maxiter", 1);

    //Connected component segmentation parameters
200
201
    AddChoice("filter.cc", "Connected components");
    SetParameterDescription("filter.cc", "Simple pixel-based connected-components algorithm with a user-defined connection condition.");
202

203
204
    AddParameter(ParameterType_String, "filter.cc.expr", "Condition");
    SetParameterDescription("filter.cc.expr", "User defined connection condition, written as a mathematical expression. Available variables are p(i)b(i), intensity_p(i) and distance (example of expression : distance < 10 )");
205
206

    // Watershed
207
208
    AddChoice("filter.watershed","Watershed");
    SetParameterDescription("filter.watershed","The traditional watershed algorithm. The height function is the gradient magnitude of the amplitude (square root of the sum of squared bands).");
209
210
211
212
213
214
215
216
217
218
219
220
221

    AddParameter(ParameterType_Float,"filter.watershed.threshold","Depth Threshold");
    SetParameterDescription("filter.watershed.threshold","Depth threshold Units in percentage of the maximum depth in the image.");
    SetDefaultParameterFloat("filter.watershed.threshold",0.01);
    SetMinimumParameterFloatValue("filter.watershed.threshold",0);
    SetMaximumParameterFloatValue("filter.watershed.threshold",1);

    AddParameter(ParameterType_Float,"filter.watershed.level","Flood Level");
    SetParameterDescription("filter.watershed.level","flood level for generating the merge tree from the initial segmentation (between 0 and 1)");
    SetDefaultParameterFloat("filter.watershed.level",0.1);
    SetMinimumParameterFloatValue("filter.watershed.level",0);
    SetMaximumParameterFloatValue("filter.watershed.level",1);

222
    AddParameter(ParameterType_Choice, "mode", "Processing mode");
223
    SetParameterDescription("mode", "Choice of processing mode, either raster or large-scale.");
224

225
226
    AddChoice("mode.vector", "Tile-based large-scale segmentation with vector output");
    SetParameterDescription("mode.vector","In this mode, the application will output a vector file or database, and process the input image piecewise. This allows to perform segmentation of very large images.");
227

228
229
    AddChoice("mode.raster", "Standard segmentation with labeled raster output");
    SetParameterDescription("mode.raster","In this mode, the application will output a standard labeled raster. This mode can not handle large data.");
230

231
    // GeoMorpho
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
    AddChoice("filter.mprofiles","Morphological profiles based segmentation");
    SetParameterDescription("filter.mprofiles","Segmentation based on morphological profiles, as described in Martino Pesaresi and Jon Alti Benediktsson, Member, IEEE: A new approach for the morphological segmentation of high resolution satellite imagery. IEEE Transactions on geoscience and remote sensing, vol. 39, NO. 2, February 2001, p. 309-320.");
    AddParameter(ParameterType_Int,"filter.mprofiles.size","Profile Size");
    SetParameterDescription("filter.mprofiles.size","Size of the profiles");
    SetDefaultParameterInt("filter.mprofiles.size",5);
    SetMinimumParameterIntValue("filter.mprofiles.size",2);
    AddParameter(ParameterType_Int,"filter.mprofiles.start","Initial radius");
    SetParameterDescription("filter.mprofiles.start","Initial radius of the structuring element (in pixels)");
    SetDefaultParameterInt("filter.mprofiles.start",1);
    SetMinimumParameterIntValue("filter.mprofiles.start",1);
    AddParameter(ParameterType_Int,"filter.mprofiles.step","Radius step.");
    SetParameterDescription("filter.mprofiles.step","Radius step along the profile (in pixels)");
    SetDefaultParameterInt("filter.mprofiles.step",1);
    SetMinimumParameterIntValue("filter.mprofiles.step",1);
    AddParameter(ParameterType_Float,"filter.mprofiles.sigma","Threshold of the final decision rule");
    SetParameterDescription("filter.mprofiles.sigma","Profiles values under the threshold will be ignored.");
    SetDefaultParameterFloat("filter.mprofiles.sigma",1.);
249
    SetMinimumParameterFloatValue("filter.mprofiles.sigma",0.);
250

251
252
253
    //Raster mode parameters
    AddParameter(ParameterType_OutputImage,  "mode.raster.out",    "Output labeled image");
    SetParameterDescription( "mode.raster.out", "The output labeled image.");
254

OTB Bot's avatar
STYLE    
OTB Bot committed
255
    //Streaming vectorization parameters
256
257
    AddParameter(ParameterType_OutputFilename, "mode.vector.out", "Output vector file");
    SetParameterDescription("mode.vector.out", "The output vector file or database (name can be anything understood by OGR)");
258

259
260
    AddParameter(ParameterType_Choice,"mode.vector.outmode","Writing mode for the output vector file");
    SetParameterDescription("mode.vector.outmode","This allows to set the writing behaviour for the output vector file. Please note that the actual behaviour depends on the file format.");
261

262
263
264
    AddChoice("mode.vector.outmode.ulco","Update output vector file, only allow to create new layers");
    SetParameterDescription("mode.vector.outmode.ulco","The output vector file is opened in update mode if existing. If the output layer already exists, the application stops, leaving it untouched.");

265
266
267
268
    AddChoice("mode.vector.outmode.ovw","Overwrite output vector file if existing.");
    SetParameterDescription("mode.vector.outmode.ovw","If the output vector file already exists, it is completely destroyed (including all its layers) and recreated from scratch.");

    AddChoice("mode.vector.outmode.ulovw","Update output vector file, overwrite existing layer");
269
    SetParameterDescription("mode.vector.outmode.ulovw","The output vector file is opened in update mode if existing. If the output layer already exists, it si completely destroyed and recreated from scratch.");
270

271
    AddChoice("mode.vector.outmode.ulu","Update output vector file, update existing layer");
272
    SetParameterDescription("mode.vector.outmode.ulu","The output vector file is opened in update mode if existing. If the output layer already exists, the new geometries are appended to the layer.");
273

274
275
276
    AddParameter(ParameterType_InputImage, "mode.vector.inmask", "Mask Image");
    SetParameterDescription("mode.vector.inmask", "Only pixels whose mask value is strictly positive will be segmented.");
    MandatoryOff("mode.vector.inmask");
277

278
279
    AddParameter(ParameterType_Empty, "mode.vector.neighbor", "8-neighbor connectivity");
    SetParameterDescription("mode.vector.neighbor",
280
                            "Activate 8-Neighborhood connectivity (default is 4).");
281
    MandatoryOff("mode.vector.neighbor");
282
    DisableParameter("mode.vector.neighbor");
283

284

285
286
287
288
    AddParameter(ParameterType_Empty,"mode.vector.stitch","Stitch polygons");
    SetParameterDescription("mode.vector.stitch", "Scan polygons on each side of tiles and stitch polygons which connect by more than one pixel.");
    MandatoryOff("mode.vector.stitch");
    EnableParameter("mode.vector.stitch");
289

290
291
    AddParameter(ParameterType_Int, "mode.vector.minsize", "Minimum object size");
    SetParameterDescription("mode.vector.minsize",
292
                            "Objects whose size is below the minimum object size (area in pixels) will be ignored during vectorization.");
293
294
295
296
    SetDefaultParameterInt("mode.vector.minsize", 1);
    SetMinimumParameterIntValue("mode.vector.minsize", 1);
    MandatoryOff("mode.vector.minsize");
    DisableParameter("mode.vector.minsize");
297

298
299
    AddParameter(ParameterType_Float, "mode.vector.simplify", "Simplify polygons");
    SetParameterDescription("mode.vector.simplify",
300
                            "Simplify polygons according to a given tolerance (in pixel). This option allows to reduce the size of the output file or database.");
301
302
    SetDefaultParameterFloat("mode.vector.simplify",0.1);
    MandatoryOff("mode.vector.simplify");
303
    DisableParameter("mode.vector.simplify");
304

305
306
307
    AddParameter(ParameterType_String, "mode.vector.layername", "Layer name");
    SetParameterDescription("mode.vector.layername", "Name of the layer in the vector file or database (default is Layer).");
    SetParameterString("mode.vector.layername", "layer");
308

309
310
311
    AddParameter(ParameterType_String, "mode.vector.fieldname", "Geometry index field name");
    SetParameterDescription("mode.vector.fieldname", "Name of the field holding the geometry index in the output vector file or database.");
    SetParameterString("mode.vector.fieldname", "DN");
312

313
314
    AddParameter(ParameterType_Int, "mode.vector.tilesize", "Tiles size");
    SetParameterDescription("mode.vector.tilesize",
315
                            "User defined tiles size for tile-based segmentation. Optimal tile size is selected according to available RAM if null.");
316
317
    SetDefaultParameterInt("mode.vector.tilesize",1024);
    SetMinimumParameterIntValue("mode.vector.tilesize",0);
318

319
320
321
322
    AddParameter(ParameterType_Int, "mode.vector.startlabel", "Starting geometry index");
    SetParameterDescription("mode.vector.startlabel", "Starting value of the geometry index field");
    SetDefaultParameterInt("mode.vector.startlabel", 1);
    SetMinimumParameterIntValue("mode.vector.startlabel", 1);
323

324
325
326
327
    AddParameter(ParameterType_StringList,"mode.vector.ogroptions","OGR options for layer creation");
    SetParameterDescription("mode.vector.ogroptions","A list of layer creation options in the form KEY=VALUE that will be passed directly to OGR without any validity checking. Options may depend on the file format, and can be found in OGR documentation.");
    MandatoryOff("mode.vector.ogroptions");

328
    // Doc example parameter settings
329
    SetExampleComment("Example of use with vector mode and watershed segmentation",0);
330
    SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_PAN.tif");
331
332
    SetDocExampleParameterValue("mode","vector");
    SetDocExampleParameterValue("mode.vector.out", "SegmentationVector.sqlite");
333
334
335
336
    SetDocExampleParameterValue("filter", "watershed");

    AddExample("Example of use with raster mode and mean-shift segmentation");
    SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_PAN.tif",1);
337
338
    SetDocExampleParameterValue("mode","raster",1);
    SetDocExampleParameterValue("mode.raster.out", "SegmentationRaster.tif uint16",1);
339
    SetDocExampleParameterValue("filter", "meanshift",1);
340
  }
341

342
343
344
345
  void DoUpdateParameters()
  {
    // Nothing to do here : all parameters are independent
  }
346

347
348
  template<class TInputImage, class TSegmentationFilter>
  FloatVectorImageType::SizeType
349
350
  GenericApplySegmentation(otb::StreamingImageToOGRLayerSegmentationFilter<TInputImage,
                           TSegmentationFilter> * streamingVectorizedFilter, TInputImage * inputImage, const otb::ogr::Layer& layer, const unsigned int outputNb)
351
352
  {
    // Retrieve tile size parameter
353
    const unsigned int tileSize = static_cast<unsigned int> (this->GetParameterInt("mode.vector.tilesize"));
354
    // Retrieve the 8-connected option
355
    bool use8connected = IsParameterEnabled("mode.vector.neighbor");
356
    // Retrieve min object size parameter
357
    const unsigned int minSize = static_cast<unsigned int> (this->GetParameterInt("mode.vector.minsize"));
358
359

    // Switch on segmentation mode
Julien Michel's avatar
Julien Michel committed
360
    const std::string segModeType = GetParameterString("mode");
361

362
    streamingVectorizedFilter->SetInput(inputImage);
Manuel Grizonnet's avatar
Manuel Grizonnet committed
363

364
    if (segModeType == "vector" && HasValue("mode.vector.inmask"))
365
      {
366
      streamingVectorizedFilter->SetInputMask(this->GetParameterUInt32Image("mode.vector.inmask"));
367
368
      otbAppLogINFO(<<"Use a mask as input." << std::endl);
      }
369
    streamingVectorizedFilter->SetOGRLayer(layer);
Manuel Grizonnet's avatar
Manuel Grizonnet committed
370

371
372
373
374
375
376
377
378
    if (tileSize != 0)
      {
      streamingVectorizedFilter->GetStreamer()->SetTileDimensionTiledStreaming(tileSize);
      }
    else
      {
      streamingVectorizedFilter->GetStreamer()->SetAutomaticTiledStreaming();
      }
Manuel Grizonnet's avatar
Manuel Grizonnet committed
379

380
381
382
383
384
    if (use8connected)
      {
      otbAppLogINFO(<<"Use 8 connected neighborhood."<<std::endl);
      }
    streamingVectorizedFilter->SetUse8Connected(use8connected);
Manuel Grizonnet's avatar
Manuel Grizonnet committed
385

386
387
388
389
390
391
    if (minSize > 1)
      {
      otbAppLogINFO(<<"Object with size under "<< minSize <<" will be suppressed."<<std::endl);
      streamingVectorizedFilter->SetFilterSmallObject(true);
      streamingVectorizedFilter->SetMinimumObjectSize(minSize);
      }
392

393
    const std::string fieldName = this->GetParameterString("mode.vector.fieldname");
394

395
    // Retrieve start label parameter
396
    const unsigned int startLabel = this->GetParameterInt("mode.vector.startlabel");
397

398
399
    streamingVectorizedFilter->SetFieldName(fieldName);
    streamingVectorizedFilter->SetStartLabel(startLabel);
400

401
    if(IsParameterEnabled("mode.vector.simplify") && GetParameterString("mode") == "vector")
402
403
      {
      streamingVectorizedFilter->SetSimplify(true);
404
      streamingVectorizedFilter->SetSimplificationTolerance(GetParameterFloat("mode.vector.simplify"));
405
406
407
408
409
410
      otbAppLogINFO(<<"Simplify the geometry." << std::endl);
      }
    else
      {
      streamingVectorizedFilter->SetSimplify(false);
      }
Manuel Grizonnet's avatar
Manuel Grizonnet committed
411

412
    if (segModeType == "vector")
OTB Bot's avatar
STYLE    
OTB Bot committed
413
      {
414
      otbAppLogINFO(<<"Large scale segmentation mode which output vector data" << std::endl);
415
416
417
418

      DisableParameter("mode.raster.out");
      EnableParameter("mode.vector.out");

419
420
421
422
423
      AddProcess(streamingVectorizedFilter->GetStreamer(), "Computing " + (dynamic_cast <ChoiceParameter *> (this->GetParameterByKey("filter")))->GetChoiceKey(GetParameterInt("filter")) + " segmentation");

      streamingVectorizedFilter->Initialize(); //must be called !
      streamingVectorizedFilter->Update(); //must be called !
      }
424
    else if (segModeType == "raster")
425
426
427
      {
      otbAppLogINFO(<<"Segmentation mode which output label image" << std::endl);

428
429
430
      DisableParameter("mode.vector.out");
      EnableParameter("mode.raster.out");

431
      streamingVectorizedFilter->GetSegmentationFilter()->SetInput(inputImage);
432
433
434
      SetParameterOutputImage<UInt32ImageType> ("mode.raster.out", dynamic_cast<UInt32ImageType *> (streamingVectorizedFilter->GetSegmentationFilter()->GetOutputs().at(outputNb).GetPointer()));
      //TODO add progress reporting in raster mode
      // AddProcess(dynamic_cast <OutputImageParameter *> (GetParameterByKey("mode.raster.out"))->GetWriter(),
435
436
437
      //            "Computing " + (dynamic_cast <ChoiceParameter *>
      //                            (this->GetParameterByKey("filter")))->GetChoiceKey(GetParameterInt("filter"))
      //            + " segmentation");
438
439
      streamingVectorizedFilter->GetSegmentationFilter()->Update();
      }
Manuel Grizonnet's avatar
Manuel Grizonnet committed
440
    return streamingVectorizedFilter->GetStreamSize();
441
  }
442

443
444
445
  void DoExecute()
  {
    // Switch on segmentation mode
Julien Michel's avatar
Julien Michel committed
446
    const std::string segModeType = GetParameterString("mode");
447
    // Switch on segmentation filter case
Julien Michel's avatar
Julien Michel committed
448
    const std::string segType = GetParameterString("filter");
449
450

    otb::ogr::DataSource::Pointer ogrDS;
451
    otb::ogr::Layer layer(NULL, false);
452

453
454
455
    std::string projRef = GetParameterFloatVectorImage("in")->GetProjectionRef();

    OGRSpatialReference oSRS(projRef.c_str());
456

457
    if (segModeType == "vector")
458
459
      {
      // Retrieve output filename as well as layer names
460
      std::string dataSourceName = GetParameterString("mode.vector.out");
461

462
      //projection ref conversion to ESRI need to be tested in case of .shp
Julien Michel's avatar
Julien Michel committed
463
      if ((dataSourceName.find(".shp") != std::string::npos) && (!projRef.empty()))
464
        {
465

466
467
        if (!(otb::GeoInformationConversion::IsESRIValidWKT(projRef)))
          {
468
          otbAppLogFATAL(<<"Image projection reference "<<std::endl<< projRef);
469
470
471
472
          itkExceptionMacro(<<"Image spatial reference can't be converted to ESRI. Use another output format (kml,SQLite,...) to overcome .shp limitation ");
          }
        }

473
474
475
476
477
      // Retrieve the output vector opening mode
      std::string outmode = GetParameterString("mode.vector.outmode");

      std::vector<std::string> options = GetParameterStringList("mode.vector.ogroptions");

478
      // Create the DataSource in the appropriate mode
479
      if (outmode == "ovw")
480
481
        {
        // Create the datasource
482
        ogrDS = otb::ogr::DataSource::New(dataSourceName, otb::ogr::DataSource::Modes::Overwrite);
483
484
485

        // and create the layer since we are in overwrite mode, the
        // datasource is blank
486
        layer = ogrDS->CreateLayer(GetParameterString("mode.vector.layername"), &oSRS, wkbMultiPolygon, options);
487
        // And create the field
488
489
        OGRFieldDefn field(this->GetParameterString("mode.vector.fieldname").c_str(), OFTInteger);
        layer.CreateField(field, true);
490
        }
491
492
      else
        if (outmode == "ulovw")
493
          {
494
495
          // Create the datasource
          ogrDS = otb::ogr::DataSource::New(dataSourceName, otb::ogr::DataSource::Modes::Update_LayerOverwrite);
496

497
498
499
500
501
502
          // and create the layer since we are in overwrite mode, the
          // datasource is blank
          layer = ogrDS->CreateLayer(GetParameterString("mode.vector.layername"), &oSRS, wkbMultiPolygon, options);
          // And create the field
          OGRFieldDefn field(this->GetParameterString("mode.vector.fieldname").c_str(), OFTInteger);
          layer.CreateField(field, true);
503

504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
          }
        else
          if (outmode == "ulu")
            {
            // Create the datasource
            ogrDS = otb::ogr::DataSource::New(dataSourceName, otb::ogr::DataSource::Modes::Update_LayerUpdate);
            // and create the layer since we are in overwrite mode, the
            // datasource is blank
            layer = ogrDS->CreateLayer(GetParameterString("mode.vector.layername"), &oSRS, wkbMultiPolygon, options);

            // And create the field if necessary
            std::string fieldName = this->GetParameterString("mode.vector.fieldname");
            OGRFeatureDefn & ogrFeatureDfn = layer.GetLayerDefn();

            if (-1 == ogrFeatureDfn.GetFieldIndex(fieldName.c_str()))
              {
              OGRFieldDefn field(fieldName.c_str(), OFTInteger);
              layer.CreateField(field, true);
              }

            }
          else
            if (outmode == "ulco")
              {
              // Create the datasource
              ogrDS = otb::ogr::DataSource::New(dataSourceName, otb::ogr::DataSource::Modes::Update_LayerCreateOnly);

              // and create the layer since we are in overwrite mode, the
              // datasource is blank
              layer = ogrDS->CreateLayer(GetParameterString("mode.vector.layername"), &oSRS, wkbMultiPolygon, options);
              // And create the field
              OGRFieldDefn field(this->GetParameterString("mode.vector.fieldname").c_str(), OFTInteger);
              layer.CreateField(field, true);
              }
            else
              {
              otbAppLogFATAL(<<"outmode not handled yet: "<< outmode);
              }
542
      }
543

544
545
546
    // The actual stream size used
    FloatVectorImageType::SizeType streamSize;

547
    if (segType == "cc")
548
      {
549
      otbAppLogINFO(<<"Use connected component segmentation."<<std::endl);
Manuel Grizonnet's avatar
Manuel Grizonnet committed
550
      ConnectedComponentStreamingVectorizedSegmentationOGRType::Pointer
551
          ccVectorizationFilter = ConnectedComponentStreamingVectorizedSegmentationOGRType::New();
Manuel Grizonnet's avatar
Manuel Grizonnet committed
552

553
      if (HasValue("mode.vector.inmask"))
Manuel Grizonnet's avatar
Manuel Grizonnet committed
554
        {
555
556
        ccVectorizationFilter->GetSegmentationFilter()->SetMaskImage(
                                                                     this->GetParameterUInt32Image("mode.vector.inmask"));
Manuel Grizonnet's avatar
Manuel Grizonnet committed
557
558
        }

559
560
561
562
563
564
      ccVectorizationFilter->GetSegmentationFilter()->GetFunctor().SetExpression(GetParameterString("filter.cc.expr"));
      streamSize = GenericApplySegmentation<FloatVectorImageType, ConnectedComponentSegmentationFilterType> (
                                                                                                             ccVectorizationFilter,
                                                                                                             this->GetParameterFloatVectorImage(
                                                                                                                                                "in"),
                                                                                                             layer, 0);
565
      }
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
      else
        if (segType == "meanshift")
          {
          otbAppLogINFO(<<"Use threaded Mean-shift segmentation."<<std::endl);

          MeanShiftVectorizedSegmentationOGRType::Pointer
              meanShiftVectorizationFilter = MeanShiftVectorizedSegmentationOGRType::New();

          //segmentation parameters
          const unsigned int
              spatialRadius = static_cast<unsigned int> (this->GetParameterInt("filter.meanshift.spatialr"));
          const float rangeRadius = static_cast<float> (this->GetParameterFloat("filter.meanshift.ranger"));
          const unsigned int
              minimumObjectSize = static_cast<unsigned int> (this->GetParameterInt("filter.meanshift.minsize"));

          const float threshold = this->GetParameterFloat("filter.meanshift.thres");
          const unsigned int
              maxIterNumber = static_cast<unsigned int> (this->GetParameterInt("filter.meanshift.maxiter"));

          meanShiftVectorizationFilter->GetSegmentationFilter()->SetSpatialBandwidth(spatialRadius);
          meanShiftVectorizationFilter->GetSegmentationFilter()->SetRangeBandwidth(rangeRadius);
          meanShiftVectorizationFilter->GetSegmentationFilter()->SetMaxIterationNumber(maxIterNumber);
          meanShiftVectorizationFilter->GetSegmentationFilter()->SetThreshold(threshold);
          meanShiftVectorizationFilter->GetSegmentationFilter()->SetMinRegionSize(minimumObjectSize);

          streamSize = this->GenericApplySegmentation<FloatVectorImageType, MeanShiftSegmentationFilterType> (
                                                                                                              meanShiftVectorizationFilter,
                                                                                                              this->GetParameterFloatVectorImage(
                                                                                                                                                 "in"),
                                                                                                              layer, 0);
          }
        else
          if (segType == "watershed")
            {
            otbAppLogINFO(<<"Using watershed segmentation."<<std::endl);
601

602
            AmplitudeFilterType::Pointer amplitudeFilter = AmplitudeFilterType::New();
603

604
            amplitudeFilter->SetInput(this->GetParameterFloatVectorImage("in"));
605

606
607
            GradientMagnitudeFilterType::Pointer gradientMagnitudeFilter = GradientMagnitudeFilterType::New();
            gradientMagnitudeFilter->SetInput(amplitudeFilter->GetOutput());
608

609
610
            StreamingVectorizedWatershedFilterType::Pointer
                watershedVectorizedFilter = StreamingVectorizedWatershedFilterType::New();
611

612
613
614
615
            watershedVectorizedFilter->GetSegmentationFilter()->SetThreshold(
                                                                             GetParameterFloat(
                                                                                               "filter.watershed.threshold"));
            watershedVectorizedFilter->GetSegmentationFilter()->SetLevel(GetParameterFloat("filter.watershed.level"));
616

617
618
619
620
621
            streamSize = this->GenericApplySegmentation<FloatImageType, WatershedSegmentationFilterType> (
                                                                                                          watershedVectorizedFilter,
                                                                                                          gradientMagnitudeFilter->GetOutput(),
                                                                                                          layer, 0);
            }
622
623

        else
624
          if (segType == "mprofiles")
625
626
627
            {
            otbAppLogINFO(<<"Using multiscale geodesic morphology segmentation."<<std::endl);

628
629
630
631
            unsigned int profileSize = GetParameterInt("filter.mprofiles.size");
            unsigned int initialValue = GetParameterInt("filter.mprofiles.start");
            unsigned int step = GetParameterInt("filter.mprofiles.step");
            double       sigma = GetParameterFloat("filter.mprofiles.sigma");
632
633
634
635
636
637


            AmplitudeFilterType::Pointer amplitudeFilter = AmplitudeFilterType::New();

            amplitudeFilter->SetInput(this->GetParameterFloatVectorImage("in"));

638
639
640
641
642
            MorphoVectorizedSegmentationOGRType::Pointer morphoVectorizedSegmentation = MorphoVectorizedSegmentationOGRType::New();
            morphoVectorizedSegmentation->GetSegmentationFilter()->SetProfileStart(initialValue);
            morphoVectorizedSegmentation->GetSegmentationFilter()->SetProfileSize(profileSize);
            morphoVectorizedSegmentation->GetSegmentationFilter()->SetProfileStep(step);
            morphoVectorizedSegmentation->GetSegmentationFilter()->SetSigma(sigma);
643

644
645
            streamSize = GenericApplySegmentation<FloatImageType, MorphologicalProfilesSegmentationFilterType> (
        morphoVectorizedSegmentation,                                                                     amplitudeFilter->GetOutput(),
646
                                                                                                             layer, 0);
647

648
            }
649
650
651
652
          else
            {
            otbAppLogFATAL(<<"non defined filtering method "<<GetParameterInt("filter")<<std::endl);
            }
653

654
    if (segModeType == "vector")
655
      {
656
657
      otbAppLogINFO(<<"Stream size: " << streamSize);

658
      ogrDS->SyncToDisk();
659

660
      // Stitching mode
661
      if (IsParameterEnabled("mode.vector.stitch"))
662
663
        {
        otbAppLogINFO(<<"Segmentation done, stiching polygons ...");
664

665
666
        FusionFilterType::Pointer fusionFilter = FusionFilterType::New();
        fusionFilter->SetInput(GetParameterFloatVectorImage("in"));
667
        fusionFilter->SetOGRLayer(layer);
668
        fusionFilter->SetStreamSize(streamSize);
669

670
671
        AddProcess(fusionFilter, "Stitching polygons");
        fusionFilter->GenerateData();
672

OTB Bot's avatar
STYLE    
OTB Bot committed
673
674
675
       //REPACK the Layer in case of Shapefile.
       //This request will remove features marked as deleted in the .dbf filename,
       //and recomputed FID for each features (without holes).
676
677
        //Note : the GetDriver() Method has not been encapsulated in otb::ogr::DataSource,
        //so we must access the OGR pointer by using .ogr()
OTB Bot's avatar
STYLE    
OTB Bot committed
678
679
680
681
682
683
       std::string driverName(ogrDS->ogr().GetDriver()->GetName());
       if ( driverName.find("ESRI Shapefile") != std::string::npos)
         {
           otbAppLogINFO(<<"REPACK the Shapefile ..."<<std::endl);
           //In Shapefile format, the name of the DaaSource is also the name of the Layer.
           std::string shpLayerName = itksys::SystemTools::GetFilenameWithoutExtension(GetParameterString("mode.vector.out"));
684

OTB Bot's avatar
STYLE    
OTB Bot committed
685
686
687
688
689
           std::string repack("REPACK ");
           repack = repack + shpLayerName;
           ogrDS->ExecuteSQL(repack, NULL, NULL);
         }
       }
690
      }
691
  }
692
693
};
}
694
695
}

696
OTB_APPLICATION_EXPORT(otb::Wrapper::Segmentation)