An error occurred while loading the file. Please try again.
-
ctraizet authored99d26241
/*
* Copyright (C) 2005-2019 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.
*/
// Wrappers
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbWrapperChoiceParameter.h"
// Segmentation filters includes
#include "otbMeanShiftSegmentationFilter.h"
#include "otbConnectedComponentMuParserFunctor.h"
#include "otbMaskMuParserFilter.h"
#include "otbVectorImageToAmplitudeImageFilter.h"
#include "itkGradientMagnitudeImageFilter.h"
#include "otbWatershedSegmentationFilter.h"
#include "otbMorphologicalProfilesSegmentationFilter.h"
// Large scale vectorization framework
#include "otbStreamingImageToOGRLayerSegmentationFilter.h"
// Fusion filter
#include "otbOGRLayerStreamStitchingFilter.h"
#include "otbGeoInformationConversion.h"
#include "otbClampImageFilter.h"
//Utils
#include "itksys/SystemTools.hxx"
#include "otbNoDataHelper.h"
namespace otb
{
namespace Wrapper
{
class Segmentation : public Application
{
public:
/** Standard class typedefs. */
typedef Segmentation Self;
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,
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
FloatVectorImageType> MeanShiftSegmentationFilterType;
// Simple connected components
typedef otb::Functor::ConnectedComponentMuParserFunctor
<FloatVectorImageType::PixelType> FunctorType;
typedef itk::ConnectedComponentFunctorImageFilter
<FloatVectorImageType,
LabelImageType,
FunctorType,
MaskImageType > ConnectedComponentSegmentationFilterType;
typedef itk::ScalarConnectedComponentImageFilter
<LabelImageType,
LabelImageType> LabeledConnectedComponentSegmentationFilterType;
// Watershed
typedef otb::VectorImageToAmplitudeImageFilter
<FloatVectorImageType,
FloatImageType> AmplitudeFilterType;
typedef itk::GradientMagnitudeImageFilter
<FloatImageType,FloatImageType> GradientMagnitudeFilterType;
typedef otb::WatershedSegmentationFilter
<FloatImageType,LabelImageType> WatershedSegmentationFilterType;
// Geodesic morphology multiscale segmentation
typedef otb::MorphologicalProfilesSegmentationFilter<FloatImageType,LabelImageType> MorphologicalProfilesSegmentationFilterType;
// mask filter
typedef otb::MaskMuParserFilter
<FloatVectorImageType, MaskImageType> MaskMuParserFilterType;
// Vectorize filters
// Home made mean-shift
typedef otb::StreamingImageToOGRLayerSegmentationFilter
<FloatVectorImageType,
MeanShiftSegmentationFilterType> MeanShiftVectorizedSegmentationOGRType;
// Connected components
typedef otb::StreamingImageToOGRLayerSegmentationFilter
<FloatVectorImageType,
ConnectedComponentSegmentationFilterType>
ConnectedComponentStreamingVectorizedSegmentationOGRType;
// Morphological profiles
typedef otb::StreamingImageToOGRLayerSegmentationFilter
<FloatImageType,MorphologicalProfilesSegmentationFilterType> MorphoVectorizedSegmentationOGRType;
typedef otb::OGRLayerStreamStitchingFilter
<FloatVectorImageType> FusionFilterType;
// Watershed
typedef otb::StreamingImageToOGRLayerSegmentationFilter
<FloatImageType,
WatershedSegmentationFilterType> StreamingVectorizedWatershedFilterType;
typedef otb::ClampImageFilter<FloatImageType, UInt32ImageType> ClampFilterType;
/** Standard macro */
itkNewMacro(Self);
itkTypeMacro(Segmentation, otb::Application);
private:
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
void DoInit() override
{
SetName("Segmentation");
SetDescription("Performs segmentation of an image, and output either a raster or a vector file. In vector mode, large input datasets are supported.");
// Documentation
SetDocName("Segmentation");
SetDocLongDescription(
"This application allows one to perform various segmentation algorithms on a multispectral image."
" 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\n"
"In 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\n"
"To 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"
" with the chosen algorithm, vectorized, and written into the output file or database. This piece-wise behavior ensure that memory will never get"
" overloaded, and that images of any size can be processed. There are few more options in the vector mode. The simplify option allows simplifying the "
"geometry"
" (i.e. remove nodes in polygons) according to a user-defined tolerance. The stitch option tries to stitch together the polygons corresponding"
" to segmented region that may have been split by the tiling scheme. ");
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."
" \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. \nThis application does not handle no data values. No data pixels will be treated as regular pixels,"
" This may lead to unexpected segmentation results and crashes.");
SetDocAuthors("OTB-Team");
SetDocSeeAlso("LargeScaleMeanShift");
AddDocTag(Tags::Segmentation);
AddParameter(ParameterType_InputImage, "in", "Input Image");
SetParameterDescription("in", "The input image to segment");
AddParameter(ParameterType_Choice, "filter", "Segmentation algorithm");
SetParameterDescription("filter", "Choice of segmentation algorithm (mean-shift by default)");
AddChoice("filter.meanshift", "Mean-Shift");
SetParameterDescription(
"filter.meanshift","OTB implementation of the Mean-Shift algorithm (multi-threaded).");
// MeanShift Parameters
AddParameter(ParameterType_Int, "filter.meanshift.spatialr", "Spatial radius");
SetParameterDescription("filter.meanshift.spatialr", "Spatial radius of the neighborhood.");
AddParameter(ParameterType_Float, "filter.meanshift.ranger", "Range radius");
SetParameterDescription("filter.meanshift.ranger", "Range radius defining the radius (expressed in radiometry unit) in the multispectral space.");
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.");
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."
" If set to 0 no pruning is done.");
SetDefaultParameterInt("filter.meanshift.spatialr", 5);
SetDefaultParameterFloat("filter.meanshift.ranger", 15.0);
SetDefaultParameterFloat("filter.meanshift.thres", 0.1);
SetMinimumParameterFloatValue("filter.meanshift.thres", 0.1);
SetDefaultParameterInt("filter.meanshift.minsize", 100);
SetMinimumParameterIntValue("filter.meanshift.minsize", 0);
SetDefaultParameterInt("filter.meanshift.maxiter", 100);
SetMinimumParameterIntValue("filter.meanshift.maxiter", 1);
//Connected component segmentation parameters
AddChoice("filter.cc", "Connected components");
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
SetParameterDescription("filter.cc", "Simple pixel-based connected-components algorithm with a user-defined connection condition.");
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 )");
// Watershed
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).");
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);
AddParameter(ParameterType_Choice, "mode", "Processing mode");
SetParameterDescription("mode", "Choice of processing mode, either raster or large-scale.");
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 performing segmentation of very large images.");
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.");
// GeoMorpho
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.);
SetMinimumParameterFloatValue("filter.mprofiles.sigma",0.);
//Raster mode parameters
AddParameter(ParameterType_OutputImage, "mode.raster.out", "Output labeled image");
SetParameterDescription( "mode.raster.out", "The output labeled image.");
SetDefaultOutputPixelType("mode.raster.out",ImagePixelType_uint32);
//Streaming vectorization parameters
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)");
AddParameter(ParameterType_Choice,"mode.vector.outmode","Writing mode for the output vector file");
SetParameterDescription("mode.vector.outmode","This allows one to set the writing behaviour for the output vector file. Please note that the actual behaviour depends on the file format.");
AddChoice("mode.vector.outmode.ulco","Update output vector file, only allow creating 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.");
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");
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.");
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
AddChoice("mode.vector.outmode.ulu","Update output vector file, update existing layer");
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.");
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");
AddParameter(ParameterType_Bool, "mode.vector.neighbor", "8-neighbor connectivity");
SetParameterDescription("mode.vector.neighbor",
"Activate 8-Neighborhood connectivity (default is 4).");
AddParameter(ParameterType_Bool,"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.");
SetParameterInt("mode.vector.stitch",1);
AddParameter(ParameterType_Int, "mode.vector.minsize", "Minimum object size");
SetParameterDescription("mode.vector.minsize",
"Objects whose size is below the minimum object size (area in pixels) will be ignored during vectorization.");
SetDefaultParameterInt("mode.vector.minsize", 1);
SetMinimumParameterIntValue("mode.vector.minsize", 1);
MandatoryOff("mode.vector.minsize");
DisableParameter("mode.vector.minsize");
AddParameter(ParameterType_Float, "mode.vector.simplify", "Simplify polygons");
SetParameterDescription("mode.vector.simplify",
"Simplify polygons according to a given tolerance (in pixel). This option allows reducing the size of the output file or database.");
SetDefaultParameterFloat("mode.vector.simplify",0.1);
MandatoryOff("mode.vector.simplify");
DisableParameter("mode.vector.simplify");
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");
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");
AddParameter(ParameterType_Int, "mode.vector.tilesize", "Tiles size");
SetParameterDescription("mode.vector.tilesize",
"User defined tiles size for tile-based segmentation. Optimal tile size is selected according to available RAM if null.");
SetDefaultParameterInt("mode.vector.tilesize",1024);
SetMinimumParameterIntValue("mode.vector.tilesize",0);
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);
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");
// Doc example parameter settings
SetExampleComment("Example of use with vector mode and watershed segmentation",0);
SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_PAN.tif");
SetDocExampleParameterValue("mode","vector");
SetDocExampleParameterValue("mode.vector.out", "SegmentationVector.sqlite");
SetDocExampleParameterValue("filter", "watershed");
AddExample("Example of use with raster mode and mean-shift segmentation");
SetDocExampleParameterValue("in", "QB_Toulouse_Ortho_PAN.tif",1);
SetDocExampleParameterValue("mode","raster",1);
SetDocExampleParameterValue("mode.raster.out", "SegmentationRaster.tif uint16",1);
SetDocExampleParameterValue("filter", "meanshift",1);
SetOfficialDocLink();
}
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
void DoUpdateParameters() override
{
// Nothing to do here : all parameters are independent
}
template<class TInputImage, class TSegmentationFilter>
FloatVectorImageType::SizeType
GenericApplySegmentation(
otb::StreamingImageToOGRLayerSegmentationFilter<
TInputImage,TSegmentationFilter> * streamingVectorizedFilter,
TInputImage * inputImage,
const otb::ogr::Layer& layer,
const unsigned int outputNb)
{
// Retrieve tile size parameter
const unsigned int tileSize = static_cast<unsigned int> (this->GetParameterInt("mode.vector.tilesize"));
// Retrieve the 8-connected option
bool use8connected = GetParameterInt("mode.vector.neighbor");
// Retrieve min object size parameter
const unsigned int minSize = static_cast<unsigned int> (this->GetParameterInt("mode.vector.minsize"));
// Switch on segmentation mode
const std::string segModeType = GetParameterString("mode");
streamingVectorizedFilter->SetInput(inputImage);
if (segModeType == "vector" && HasValue("mode.vector.inmask"))
{
streamingVectorizedFilter->SetInputMask(m_ClampFilter->GetOutput());
otbAppLogINFO(<<"Use a mask as input." << std::endl);
}
streamingVectorizedFilter->SetOGRLayer(layer);
if (tileSize != 0)
{
streamingVectorizedFilter->GetStreamer()->SetTileDimensionTiledStreaming(tileSize);
}
else
{
streamingVectorizedFilter->GetStreamer()->SetAutomaticTiledStreaming();
}
if (use8connected)
{
otbAppLogINFO(<<"Use 8 connected neighborhood."<<std::endl);
}
streamingVectorizedFilter->SetUse8Connected(use8connected);
if (minSize > 1)
{
otbAppLogINFO(<<"Object with size under "<< minSize <<" will be suppressed."<<std::endl);
streamingVectorizedFilter->SetFilterSmallObject(true);
streamingVectorizedFilter->SetMinimumObjectSize(minSize);
}
const std::string fieldName = this->GetParameterString("mode.vector.fieldname");
// Retrieve start label parameter
const unsigned int startLabel = this->GetParameterInt("mode.vector.startlabel");
streamingVectorizedFilter->SetFieldName(fieldName);
streamingVectorizedFilter->SetStartLabel(startLabel);
if(IsParameterEnabled("mode.vector.simplify") && GetParameterString("mode") == "vector")
{
streamingVectorizedFilter->SetSimplify(true);
streamingVectorizedFilter->SetSimplificationTolerance(GetParameterFloat("mode.vector.simplify"));
otbAppLogINFO(<<"Simplify the geometry." << std::endl);
}
else
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
{
streamingVectorizedFilter->SetSimplify(false);
}
if (segModeType == "vector")
{
otbAppLogINFO(<<"Large scale segmentation mode which output vector data" << std::endl);
DisableParameter("mode.raster.out");
EnableParameter("mode.vector.out");
AddProcess(streamingVectorizedFilter->GetStreamer(), "Computing "
+ this->GetParameterString("filter")
+ " segmentation");
streamingVectorizedFilter->Initialize(); //must be called !
streamingVectorizedFilter->Update(); //must be called !
}
else if (segModeType == "raster")
{
otbAppLogINFO(<<"Segmentation mode which output label image" << std::endl);
DisableParameter("mode.vector.out");
EnableParameter("mode.raster.out");
streamingVectorizedFilter->GetSegmentationFilter()->SetInput(inputImage);
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(),
// "Computing " + (dynamic_cast <ChoiceParameter *>
// (this->GetParameterByKey("filter")))->GetChoiceKey(GetParameterInt("filter"))
// + " segmentation");
streamingVectorizedFilter->GetSegmentationFilter()->Update();
}
return streamingVectorizedFilter->GetStreamSize();
}
void DoExecute() override
{
// Switch on segmentation mode
const std::string segModeType = GetParameterString("mode");
// Switch on segmentation filter case
const std::string segType = GetParameterString("filter");
otb::ogr::DataSource::Pointer ogrDS;
otb::ogr::Layer layer(nullptr, false);
std::string projRef = GetParameterFloatVectorImage("in")->GetProjectionRef();
std::vector<bool> noDataFlags;
std::vector<double> noDataValues;
itk::MetaDataDictionary &dict = GetParameterFloatVectorImage("in")->GetMetaDataDictionary();
bool ret = otb::ReadNoDataFlags(dict,noDataFlags,noDataValues);
if (ret)
{
otbAppLogWARNING("The input image has no data values but this application does not handle no-data. No-data pixels"
" will be treated as regular pixels.");
}
OGRSpatialReference oSRS(projRef.c_str());
if (segModeType == "vector")
{
// Retrieve output filename as well as layer names
std::string dataSourceName = GetParameterString("mode.vector.out");
//projection ref conversion to ESRI need to be tested in case of .shp
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
if ((dataSourceName.find(".shp") != std::string::npos) && (!projRef.empty()))
{
if (!(otb::GeoInformationConversion::IsESRIValidWKT(projRef)))
{
otbAppLogFATAL(<<"Image projection reference "<<std::endl<< projRef);
itkExceptionMacro(<<"Image spatial reference can't be converted to ESRI. Use another output format (kml,SQLite,...) to overcome .shp limitation ");
}
}
// Retrieve the output vector opening mode
std::string outmode = GetParameterString("mode.vector.outmode");
std::vector<std::string> options = GetParameterStringList("mode.vector.ogroptions");
// Create the DataSource in the appropriate mode
if (outmode == "ovw")
{
// Create the datasource
ogrDS = otb::ogr::DataSource::New(dataSourceName, otb::ogr::DataSource::Modes::Overwrite);
// 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 if (outmode == "ulovw")
{
// Create the datasource
ogrDS = otb::ogr::DataSource::New(dataSourceName, otb::ogr::DataSource::Modes::Update_LayerOverwrite);
// 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 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
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
{
otbAppLogFATAL(<<"outmode not handled yet: "<< outmode);
}
}
// handle mask
if (HasValue("mode.vector.inmask"))
{
m_ClampFilter = ClampFilterType::New();
m_ClampFilter->SetInput( this->GetParameterFloatImage("mode.vector.inmask"));
}
// The actual stream size used
FloatVectorImageType::SizeType streamSize;
if (segType == "cc")
{
otbAppLogINFO(<<"Use connected component segmentation."<<std::endl);
ConnectedComponentStreamingVectorizedSegmentationOGRType::Pointer
ccVectorizationFilter = ConnectedComponentStreamingVectorizedSegmentationOGRType::New();
if (HasValue("mode.vector.inmask"))
{
ccVectorizationFilter->GetSegmentationFilter()->SetMaskImage(
m_ClampFilter->GetOutput());
}
ccVectorizationFilter->GetSegmentationFilter()->GetFunctor().SetExpression(GetParameterString("filter.cc.expr"));
streamSize = GenericApplySegmentation<FloatVectorImageType,ConnectedComponentSegmentationFilterType>(
ccVectorizationFilter,
this->GetParameterFloatVectorImage("in"),
layer,
0);
}
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);
AmplitudeFilterType::Pointer amplitudeFilter = AmplitudeFilterType::New();
631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700
amplitudeFilter->SetInput(this->GetParameterFloatVectorImage("in"));
GradientMagnitudeFilterType::Pointer gradientMagnitudeFilter = GradientMagnitudeFilterType::New();
gradientMagnitudeFilter->SetInput(amplitudeFilter->GetOutput());
StreamingVectorizedWatershedFilterType::Pointer
watershedVectorizedFilter = StreamingVectorizedWatershedFilterType::New();
watershedVectorizedFilter->GetSegmentationFilter()->SetThreshold(
GetParameterFloat("filter.watershed.threshold"));
watershedVectorizedFilter->GetSegmentationFilter()->SetLevel(GetParameterFloat("filter.watershed.level"));
streamSize = this->GenericApplySegmentation<FloatImageType,WatershedSegmentationFilterType>(
watershedVectorizedFilter,
gradientMagnitudeFilter->GetOutput(),
layer,
0);
}
else if (segType == "mprofiles")
{
otbAppLogINFO(<<"Using multiscale geodesic morphology segmentation."<<std::endl);
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");
AmplitudeFilterType::Pointer amplitudeFilter = AmplitudeFilterType::New();
amplitudeFilter->SetInput(this->GetParameterFloatVectorImage("in"));
MorphoVectorizedSegmentationOGRType::Pointer morphoVectorizedSegmentation = MorphoVectorizedSegmentationOGRType::New();
morphoVectorizedSegmentation->GetSegmentationFilter()->SetProfileStart(initialValue);
morphoVectorizedSegmentation->GetSegmentationFilter()->SetProfileSize(profileSize);
morphoVectorizedSegmentation->GetSegmentationFilter()->SetProfileStep(step);
morphoVectorizedSegmentation->GetSegmentationFilter()->SetSigma(sigma);
streamSize = GenericApplySegmentation<FloatImageType,MorphologicalProfilesSegmentationFilterType>(
morphoVectorizedSegmentation,
amplitudeFilter->GetOutput(),
layer,
0);
}
else
{
otbAppLogFATAL(<<"non defined filtering method "<<GetParameterInt("filter")<<std::endl);
}
if (segModeType == "vector")
{
otbAppLogINFO(<<"Stream size: " << streamSize);
ogrDS->SyncToDisk();
// Stitching mode
if (GetParameterInt("mode.vector.stitch"))
{
otbAppLogINFO(<<"Segmentation done, stiching polygons ...");
FusionFilterType::Pointer fusionFilter = FusionFilterType::New();
fusionFilter->SetInput(GetParameterFloatVectorImage("in"));
fusionFilter->SetOGRLayer(layer);
fusionFilter->SetStreamSize(streamSize);
AddProcess(fusionFilter, "Stitching polygons");
fusionFilter->GenerateData();
//REPACK the Layer in case of Shapefile.
//This request will remove features marked as deleted in the .dbf filename,
701702703704705706707708709710711712713714715716717718719720721722723724725
//and recomputed FID for each features (without holes).
//Note : the GetDriver() Method has not been encapsulated in otb::ogr::DataSource,
//so we must access the OGR pointer by using .ogr()
std::string driverName(ogrDS->ogr().GetDriverName());
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"));
std::string repack("REPACK ");
repack = repack + shpLayerName;
ogrDS->ExecuteSQL(repack, nullptr, nullptr);
}
}
}
}
ClampFilterType::Pointer m_ClampFilter;
};
}
}
OTB_APPLICATION_EXPORT(otb::Wrapper::Segmentation)