diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..0443e71fc509dd4078ea4346c1106830beee6316
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,11 @@
+cmake_minimum_required(VERSION 2.8.9)
+
+project(OTBSLIC)
+
+if(NOT OTB_SOURCE_DIR)
+ find_package(OTB REQUIRED)
+ list(APPEND CMAKE_MODULE_PATH ${OTB_CMAKE_DIR})
+ include(OTBModuleExternal)
+else()
+ otb_module_impl()
+endif()
diff --git a/app/CMakeLists.txt b/app/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..73dec0b3ed22bd55817154c4386ad9d70bf14fa4
--- /dev/null
+++ b/app/CMakeLists.txt
@@ -0,0 +1,7 @@
+message(${${otb-module}_LIBRARIES})
+OTB_CREATE_APPLICATION(
+ NAME SLIC
+ SOURCES
+ otbSLIC.cxx
+ LINK_LIBRARIES ${${otb-module}_LIBRARIES}
+)
diff --git a/app/CMakeLists.txt~ b/app/CMakeLists.txt~
new file mode 100644
index 0000000000000000000000000000000000000000..646f32a74defe12cc1961060ee877707b822e450
--- /dev/null
+++ b/app/CMakeLists.txt~
@@ -0,0 +1,7 @@
+message(${${otb-module}_LIBRARIES})
+OTB_CREATE_APPLICATION(
+ NAME SLIC
+ SOURCES
+ otbSLIC.cxx
+ LINK_LIBRARIES ${${otb-module}_LIBRARIES} OTBSLIC
+)
diff --git a/app/otbSLIC.cxx b/app/otbSLIC.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..b39390b74b4e3ad772dda6a8df1acc7bf36a4e97
--- /dev/null
+++ b/app/otbSLIC.cxx
@@ -0,0 +1,243 @@
+/*
+ * Copyright (C) 2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part a Remote Module for Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program. If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+#include "otbWrapperApplication.h"
+#include "otbWrapperApplicationFactory.h"
+#include "otbImage.h"
+#include "otbMultiToMonoChannelExtractROI.h"
+#include "otbOGRDataSourceWrapper.h"
+#include "otbOGRFeatureWrapper.h"
+
+#include "otbSLICScheduler.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+// #include "otbSegmentCharacteristicsFilter.h"
+#include "otbPersistentFilterStreamingDecorator.h"
+#include "itkLabelToRGBImageFilter.h"
+
+namespace otb
+{
+ namespace Wrapper
+ {
+ class SLIC : public Application
+ {
+ public:
+ typedef SLIC Self;
+ typedef Application Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+ typedef otb::VectorImage<float> VectorImageType;
+ typedef unsigned int LabelType;
+ typedef otb::Image<LabelType> LabelImageType;
+ typedef otb::SLICScheduler<VectorImageType, LabelImageType> SLICSchedulerType;
+ typedef otb::ImageFileReader<VectorImageType> ReaderType;
+ typedef otb::ImageFileReader<LabelImageType> LabelReaderType;
+ typedef otb::ImageFileWriter<LabelImageType> WriterType;
+ typedef typename otb::VectorImage<double>::PixelType SampleType;
+
+ typedef otb::VectorImage<unsigned char> RGBImageType;
+ typedef RGBImageType::PixelType RGBPixelType;
+ typedef otb::ImageFileWriter<RGBImageType> RGBWriterType;
+
+ // Segmentation contrast maximisation LUT
+ typedef itk::LabelToRGBImageFilter <LabelImageType, RGBImageType> LabelToRGBFilterType;
+
+ itkNewMacro(Self);
+ itkTypeMacro(Merging, otb::Application);
+
+ private:
+
+
+ //Tiling modes
+ enum TilingMode
+ {
+ AUTO,
+ MANUAL
+ };
+
+ void DoInit()
+ {
+ SetName("SLIC");
+ SetDescription("");
+
+ SetDocName("");
+ SetDocLongDescription("");
+ SetDocLimitations("");
+ SetDocAuthors("");
+ SetDocSeeAlso("");
+
+ AddParameter(ParameterType_InputImage, "in", "Input image");
+ SetParameterDescription( "in", "The input image." );
+
+ AddParameter(ParameterType_OutputImage, "out", "Output raster");
+ SetParameterDescription( "out", "Output raster containing the segment label IDs" );
+ SetDefaultOutputPixelType("out",ImagePixelType_uint32);
+
+ AddParameter(ParameterType_Float,"spw","Spatial Width");
+ SetParameterDescription("spw","Initial width of segment grid, in pixels. Is approximately the average width of final segments");
+ SetDefaultParameterFloat("spw",100);
+ SetMinimumParameterFloatValue("spw",2);
+
+ AddParameter(ParameterType_Float,"dw","Distance weight");
+ SetParameterDescription("dw","Importance of spatial distance w.r.t. spectral distance");
+ SetDefaultParameterFloat("dw",1);
+ SetMinimumParameterFloatValue("dw",0);
+
+ AddParameter(ParameterType_Int,"maxit","Maximum number of iterations");
+ SetParameterDescription("maxit","Maximum number of iterations of the SLIC algorithm.");
+ SetDefaultParameterInt("maxit",100);
+ SetMinimumParameterIntValue("maxit",1);
+
+ AddParameter(ParameterType_Float,"thresh","Stopping threshold");
+ SetParameterDescription("thresh","Minimum residial for stopping the algorithm.");
+ SetDefaultParameterFloat("thresh",0.001);
+ SetMinimumParameterFloatValue("thresh",0);
+
+ AddParameter(ParameterType_Float,"margin","Margin for stabilizing");
+ SetParameterDescription("margin","Margin in superpixel widths for stabilizing the scaled algorithm");
+ SetDefaultParameterFloat("margin",3);
+ SetMinimumParameterFloatValue("margin",0);
+
+ AddParameter(ParameterType_Directory, "tmpdir", "Temporary directory");
+ SetParameterDescription( "tmpdir", "Temporary directory for storing parts of the output." );
+
+ AddParameter(ParameterType_Choice,"tiling","Tiling mode (auto or manual)");
+ AddChoice("tiling.auto","Automatic tiling mode");
+ AddChoice("tiling.manual","Manual tiling mode");
+
+ AddParameter(ParameterType_Int, "tiling.auto.ram", "Available ram (Mb)");
+ SetParameterDescription("tiling.auto.ram", "Available ram to be used per MPI process in Mb");
+ SetDefaultParameterInt("tiling.auto.ram", 128);
+ SetMinimumParameterIntValue("tiling.auto.ram", 0);
+
+ AddParameter(ParameterType_Int, "tiling.manual.nx", "Number of tiles (X-axis)");
+ SetParameterDescription("tiling.manual.nx", "Number of tiles along the X-axis.");
+ SetMinimumParameterIntValue("tiling.manual.nx", 1);
+
+ AddParameter(ParameterType_Int, "tiling.manual.ny", "Number of tiles (Y-axis)");
+ SetParameterDescription("tiling.manual.ny", "Number of tiles along the Y-axis.");
+ SetMinimumParameterIntValue("tiling.manual.ny", 1);
+
+ GDALSetCacheMax(0);
+ }
+
+ void DoUpdateParameters()
+ {
+ }
+
+ void DoExecute()
+ {
+
+
+ typename otb::MPIConfig::Pointer mpiConfig = otb::MPIConfig::Instance();
+ const unsigned int myRank = mpiConfig->GetMyRank();
+
+ VectorImageType::Pointer imageIn = GetParameterImage("in");
+ const std::string inputName = GetParameterString("in");
+ imageIn->UpdateOutputInformation();
+ const std::string outputName = GetParameterString("out");
+ const unsigned int spatialWidth = GetParameterFloat("spw");
+ const double distanceWeight = GetParameterFloat("dw");
+ const unsigned int maxIterations = GetParameterInt("maxit");
+ const double threshold = GetParameterFloat("thresh");
+ unsigned int margin = GetParameterInt("margin");
+ unsigned int nbTilesX = 1;
+ unsigned int nbTilesY = 1;
+
+ std::stringstream s_prefix;
+ s_prefix << GetParameterString("tmpdir") << "/" << "temp" ;
+ std::string prefix = s_prefix.str();
+ s_prefix << ".vrt";
+ std::string outputVrt = s_prefix.str();
+
+ switch(GetParameterInt("tiling"))
+ {
+ case AUTO :
+ {
+ const unsigned int maxMemory = GetParameterInt("tiling.auto.ram");
+ //Auto Tiling mode
+ const unsigned int X = imageIn->GetLargestPossibleRegion().GetSize()[0];
+ const unsigned int Y = imageIn->GetLargestPossibleRegion().GetSize()[1];
+ const unsigned int nbComps = imageIn->GetNumberOfComponentsPerPixel();
+ const unsigned int pin = 2*nbComps; //because reader casts to unsigned short
+ const unsigned int pout = 4; //unsigned int label
+ const unsigned int pc = (nbComps+2)*8+4; //centroid size
+
+ //Reading config
+ double alpha = 2*pin;
+ double beta = (2*pin + pout)*2*spatialWidth*margin;
+ const unsigned int readingNx = vcl_ceil(X*alpha/(vcl_sqrt(4*beta*beta + alpha*maxMemory*1024*1024) - beta));
+ const unsigned int readingNy = vcl_ceil(Y*alpha/(vcl_sqrt(4*beta*beta + alpha*maxMemory*1024*1024) - beta));
+
+ //Algo config
+ alpha = pin+pout+(double)2*pc/spatialWidth/spatialWidth;
+ beta = (pin + pout)*2*spatialWidth*margin;
+ const unsigned int algoNx = vcl_ceil(X*alpha/(vcl_sqrt(4*beta*beta + alpha*maxMemory*1024*1024) - beta));
+ const unsigned int algoNy = vcl_ceil(Y*alpha/(vcl_sqrt(4*beta*beta + alpha*maxMemory*1024*1024) - beta));
+
+ nbTilesX = readingNx > algoNx ? readingNx : algoNx;
+ nbTilesY = readingNy > algoNy ? readingNy : algoNy;
+
+ std::cout << "Auto tiling mode selected, available RAM = " << maxMemory << "MB" << std::endl;
+ std::cout << "Number of x tiles = "<< nbTilesX << std::endl;
+ std::cout << "Number of y tiles = "<< nbTilesY << std::endl;
+ break;
+ };
+ case MANUAL :
+ {
+ nbTilesX = GetParameterInt("tiling.manual.nx");
+ nbTilesY = GetParameterInt("tiling.manual.ny");
+ };
+ }
+ //SlicScheduler will write the tiled images and a .vrt in the tmpdir
+ std::cout << "Starting segmentation" << "\n";
+ SLICSchedulerType::Pointer slicScheduler = SLICSchedulerType::New();
+ slicScheduler->SetInputImage(imageIn);
+ slicScheduler->SetInputName(inputName);
+ slicScheduler->SetPrefix(prefix);
+ slicScheduler->SetSpatialWidth(spatialWidth);
+ slicScheduler->SetSpatialDistanceWeight(distanceWeight);
+ slicScheduler->SetMaxIterationNumber(maxIterations);
+ slicScheduler->SetThreshold(threshold);
+ slicScheduler->SetNbTilesX(nbTilesX);
+ slicScheduler->SetNbTilesY(nbTilesY);
+ slicScheduler->SetMargin(margin);
+ slicScheduler->Run();
+
+ std::cout << "Segmentation finished" << "\n";
+
+ //Rank 0 thread must have time to write the .vrt
+ mpiConfig->barrier();
+
+ //Convert .vrt into .tif using MPI : requires SPTW
+ LabelReaderType::Pointer labelReader = LabelReaderType::New();
+ labelReader->SetFileName(outputVrt);
+ labelReader->Update();
+ SetParameterOutputImage<UInt32ImageType> ("out",dynamic_cast<UInt32ImageType*>(labelReader->GetOutput()));
+
+ }
+ };
+ }
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::SLIC)
+
+
diff --git a/app/otbSLIC.cxx~ b/app/otbSLIC.cxx~
new file mode 100644
index 0000000000000000000000000000000000000000..45547e861caf856026d35b807e45f9c3590d34ea
--- /dev/null
+++ b/app/otbSLIC.cxx~
@@ -0,0 +1,244 @@
+/*
+ * Copyright (C) 2018 Centre National d'Etudes Spatiales (CNES)
+ *
+ * This file is part a Remote Module for Orfeo Toolbox
+ *
+ * https://www.orfeo-toolbox.org/
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Affero General Public License as
+ * published by the Free Software Foundation, either version 3 of the
+ * License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU Affero General Public License for more details.
+ *
+ * You should have received a copy of the GNU Affero General Public License
+ * along with this program. If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+#include "otbWrapperApplication.h"
+#include "otbWrapperApplicationFactory.h"
+#include "otbImage.h"
+#include "otbMultiToMonoChannelExtractROI.h"
+#include "otbOGRDataSourceWrapper.h"
+#include "otbOGRFeatureWrapper.h"
+
+#include "otbSLICScheduler.h"
+#include "otbImageFileReader.h"
+#include "otbImageFileWriter.h"
+// #include "otbSegmentCharacteristicsFilter.h"
+#include "otbPersistentFilterStreamingDecorator.h"
+#include "itkLabelToRGBImageFilter.h"
+
+namespace otb
+{
+ namespace Wrapper
+ {
+ class SLIC : public Application
+ {
+ public:
+ typedef SLIC Self;
+ typedef Application Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+ typedef otb::VectorImage<float> VectorImageType;
+ typedef unsigned int LabelType;
+ typedef otb::Image<LabelType> LabelImageType;
+ typedef otb::SLICScheduler<VectorImageType, LabelImageType> SLICSchedulerType;
+ typedef otb::ImageFileReader<VectorImageType> ReaderType;
+ typedef otb::ImageFileReader<LabelImageType> LabelReaderType;
+ typedef otb::ImageFileWriter<LabelImageType> WriterType;
+ typedef typename otb::VectorImage<double>::PixelType SampleType;
+
+ typedef otb::VectorImage<unsigned char> RGBImageType;
+ typedef RGBImageType::PixelType RGBPixelType;
+ typedef otb::ImageFileWriter<RGBImageType> RGBWriterType;
+
+ // Segmentation contrast maximisation LUT
+ typedef itk::LabelToRGBImageFilter <LabelImageType, RGBImageType> LabelToRGBFilterType;
+
+ itkNewMacro(Self);
+ itkTypeMacro(Merging, otb::Application);
+
+ private:
+
+
+ //Tiling modes
+ enum TilingMode
+ {
+ AUTO,
+ MANUAL
+ };
+
+ void DoInit()
+ {
+ SetName("SLIC");
+ SetDescription("");
+
+ SetDocName("");
+ SetDocLongDescription("");
+ SetDocLimitations("");
+ SetDocAuthors("");
+ SetDocSeeAlso("");
+
+ AddParameter(ParameterType_InputImage, "in", "Input image");
+ SetParameterDescription( "in", "The input image." );
+
+ AddParameter(ParameterType_OutputImage, "out", "Output raster");
+ SetParameterDescription( "out", "Output raster containing the segment label IDs" );
+ SetDefaultOutputPixelType("out",ImagePixelType_uint32);
+
+ AddParameter(ParameterType_Float,"spw","Spatial Width");
+ SetParameterDescription("spw","Initial width of segment grid, in pixels. Is approximately the average width of final segments");
+ SetDefaultParameterFloat("spw",100);
+ SetMinimumParameterFloatValue("spw",2);
+
+ AddParameter(ParameterType_Float,"dw","Distance weight");
+ SetParameterDescription("dw","Importance of spatial distance w.r.t. spectral distance");
+ SetDefaultParameterFloat("dw",1);
+ SetMinimumParameterFloatValue("dw",0);
+
+ AddParameter(ParameterType_Int,"maxit","Maximum number of iterations");
+ SetParameterDescription("maxit","Maximum number of iterations of the SLIC algorithm.");
+ SetDefaultParameterInt("maxit",100);
+ SetMinimumParameterIntValue("maxit",1);
+
+ AddParameter(ParameterType_Float,"thresh","Stopping threshold");
+ SetParameterDescription("thresh","Minimum residial for stopping the algorithm.");
+ SetDefaultParameterFloat("thresh",0.001);
+ SetMinimumParameterFloatValue("thresh",0);
+
+ AddParameter(ParameterType_Float,"margin","Margin for stabilizing");
+ SetParameterDescription("margin","Margin in superpixel widths for stabilizing the scaled algorithm");
+ SetDefaultParameterFloat("margin",3);
+ SetMinimumParameterFloatValue("margin",0);
+
+ AddParameter(ParameterType_Directory, "tmpdir", "Temporary directory");
+ SetParameterDescription( "tmpdir", "Temporary directory for storing parts of the output." );
+
+ AddParameter(ParameterType_Choice,"tiling","Tiling mode (auto or manual)");
+ AddChoice("tiling.auto","Automatic tiling mode");
+ AddChoice("tiling.manual","Manual tiling mode");
+
+ AddParameter(ParameterType_Int, "tiling.auto.ram", "Available ram (Mb)");
+ SetParameterDescription("tiling.auto.ram", "Available ram to be used per MPI process in Mb");
+ SetDefaultParameterInt("tiling.auto.ram", 128);
+ SetMinimumParameterIntValue("tiling.auto.ram", 0);
+
+ AddParameter(ParameterType_Int, "tiling.manual.nx", "Number of tiles (X-axis)");
+ SetParameterDescription("tiling.manual.nx", "Number of tiles along the X-axis.");
+ SetMinimumParameterIntValue("tiling.manual.nx", 1);
+
+ AddParameter(ParameterType_Int, "tiling.manual.ny", "Number of tiles (Y-axis)");
+ SetParameterDescription("tiling.manual.ny", "Number of tiles along the Y-axis.");
+ SetMinimumParameterIntValue("tiling.manual.ny", 1);
+
+ GDALSetCacheMax(0);
+ }
+
+ void DoUpdateParameters()
+ {
+ }
+
+ void DoExecute()
+ {
+
+
+ typename otb::MPIConfig::Pointer mpiConfig = otb::MPIConfig::Instance();
+ const unsigned int myRank = mpiConfig->GetMyRank();
+
+ VectorImageType::Pointer imageIn = GetParameterImage("in");
+ const std::string inputName = GetParameterString("in");
+ imageIn->UpdateOutputInformation();
+ const std::string outputName = GetParameterString("out");
+ const unsigned int spatialWidth = GetParameterFloat("spw");
+ const double distanceWeight = GetParameterFloat("dw");
+ const unsigned int maxIterations = GetParameterInt("maxit");
+ const double threshold = GetParameterFloat("thresh");
+ unsigned int margin = GetParameterInt("margin");
+ unsigned int nbTilesX = 1;
+ unsigned int nbTilesY = 1;
+
+ std::stringstream s_prefix;
+ s_prefix << GetParameterString("tmpdir") << "/" << "temp" ;
+ std::string prefix = s_prefix.str();
+ s_prefix << ".vrt";
+ std::string outputVrt = s_prefix.str();
+
+ switch(GetParameterInt("tiling"))
+ {
+ case AUTO :
+ {
+ unsigned int maxMemory = GetParameterInt("tiling.auto.ram");
+ //Auto Tiling mode
+ const unsigned int X = imageIn->GetLargestPossibleRegion().GetSize()[0];
+ const unsigned int Y = imageIn->GetLargestPossibleRegion().GetSize()[1];
+ const unsigned int nbComps = imageIn->GetNumberOfComponentsPerPixel();
+ const unsigned int pin = 2*nbComps; //because reader casts to unsigned short
+ const unsigned int pout = 4; //unsigned int label
+ const unsigned int pc = (nbComps+2)*8+4; //centroid size
+
+ //Reading config
+ double alpha = 2*pin;
+ double beta = (2*pin + pout)*2*spatialWidth*margin;
+ const unsigned int readingNx = vcl_ceil(X*alpha/(vcl_sqrt(4*beta*beta + alpha*maxMemory*1024*1024) - beta));
+ const unsigned int readingNy = vcl_ceil(Y*alpha/(vcl_sqrt(4*beta*beta + alpha*maxMemory*1024*1024) - beta));
+
+ //Algo config
+ alpha = pin+pout+(double)2*pc/spatialWidth/spatialWidth;
+ beta = (pin + pout)*2*spatialWidth*margin;
+ const unsigned int algoNx = vcl_ceil(X*alpha/(vcl_sqrt(4*beta*beta + alpha*maxMemory*1024*1024) - beta));
+ const unsigned int algoNy = vcl_ceil(Y*alpha/(vcl_sqrt(4*beta*beta + alpha*maxMemory*1024*1024) - beta));
+
+ nbTilesX = readingNx > algoNx ? readingNx : algoNx;
+ nbTilesY = readingNy > algoNy ? readingNy : algoNy;
+
+ std::cout << "Auto tiling mode selected, available RAM = " << maxMemory << "MB" << std::endl;
+ std::cout << "Number of x tiles = "<< nbTilesX << std::endl;
+ std::cout << "Number of y tiles = "<< nbTilesY << std::endl;
+ break;
+ };
+ case MANUAL :
+ {
+ nbTilesX = GetParameterInt("tiling.manual.nx");
+ nbTilesY = GetParameterInt("tiling.manual.ny");
+ };
+ }
+ std::cout << "SLIC real 3" << "\n";
+ SLICSchedulerType::Pointer slicScheduler = SLICSchedulerType::New();
+ slicScheduler->SetInputImage(imageIn);
+ slicScheduler->SetInputName(inputName);
+ slicScheduler->SetPrefix(prefix);
+ slicScheduler->SetSpatialWidth(spatialWidth);
+ slicScheduler->SetSpatialDistanceWeight(distanceWeight);
+ slicScheduler->SetMaxIterationNumber(maxIterations);
+ slicScheduler->SetThreshold(threshold);
+ slicScheduler->SetNbTilesX(nbTilesX);
+ slicScheduler->SetNbTilesY(nbTilesY);
+ slicScheduler->SetMargin(margin);
+ slicScheduler->Run();
+
+ std::cout << "App finished" << "\n";
+
+ // if (myRank == 0) {
+
+ mpiConfig->barrier();
+ LabelReaderType::Pointer labelReader = LabelReaderType::New();
+ labelReader->SetFileName(outputVrt);
+ labelReader->Update();
+
+ SetParameterOutputImage<UInt32ImageType> ("out",dynamic_cast<UInt32ImageType*>(labelReader->GetOutput()));
+
+ //TODO cleanup temp dir
+ // }
+ }
+ };
+ }
+}
+
+OTB_APPLICATION_EXPORT(otb::Wrapper::SLIC)
+
+
diff --git a/include/otbSLICFilter.h b/include/otbSLICFilter.h
new file mode 100644
index 0000000000000000000000000000000000000000..dd0907e2d897c5e910c22c2bebd3a139eb7a38f3
--- /dev/null
+++ b/include/otbSLICFilter.h
@@ -0,0 +1,154 @@
+#ifndef otbSLICFilter_h
+#define otbSLICFilter_h
+
+#include "itkMacro.h"
+//OTB, ITK includes
+#include "otbVectorImage.h"
+#include "otbImage.h"
+#include "otbStreamingTraits.h"
+
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkImageRegionExclusionConstIteratorWithIndex.h"
+#include "itkNeighborhoodIterator.h"
+
+#include "itkProcessObject.h"
+//Own includes
+#include "simplePointCalculator.h"
+
+//STD includes
+#include <algorithm>//std::find
+#include <unordered_set>
+#include <unordered_map>
+namespace otb{
+ enum RelativePosition{
+ POS_TOP,
+ POS_RIGHT,
+ POS_BOTTOM,
+ POS_LEFT
+ };
+
+ enum NeighborhoodRelativePosition{
+ NBH_TOP,
+ NBH_TOP_RIGHT,
+ NBH_RIGHT,
+ NBH_BOTTOM_RIGHT,
+ NBH_BOTTOM,
+ NBH_BOTTOM_LEFT,
+ NBH_LEFT,
+ NBH_TOP_LEFT
+ };
+ struct ProcessingTile
+ {
+ long int rows[2]; // lower and upper rows (-1 means that the row has not be considered)
+ long int columns[2]; // lower and upper columns (-1 means that the row has not be considered)
+ long int tileNeighbors[8]; // tile Neighbors at (top, top right, right, bottom right, bottom, bottom left, left, top left)
+ long int margin[4]; // Is there a margin at top, left, bottom or right
+ otb::VectorImage<double>::RegionType region; // The image region
+
+ // Temporary files
+ std::string nodeFileName;
+ std::string edgeFileName;
+ std::string nodeMarginFileName;
+ std::string edgeMarginFileName;
+ };
+ template <class TInputImage, class TOutputLabelImage>
+ class ITK_EXPORT SLICFilter : public itk::ProcessObject
+ {
+ public:
+ typedef SLICFilter Self;
+ typedef itk::ProcessObject Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Runtime information support. */
+ itkTypeMacro(SLICFilter, itk::ProcessObject);
+
+ /** ImageDimension constants */
+ itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension);
+
+ typedef double RealType;
+ typedef typename TOutputLabelImage::PixelType LabelType;
+ typedef TOutputLabelImage OutputLabelImageType;
+ typedef TInputImage VectorImageType;
+ typedef typename TInputImage::RegionType RegionType;
+ typedef typename std::vector<double> SampleType;
+ typedef typename VectorImage<double>::PixelType PixelType;
+ typedef std::unordered_map<LabelType,SampleType> ListSampleType;
+ typedef std::unordered_map<LabelType,unsigned int> CountContainerType;
+ typedef std::vector<LabelType> NeighborhoodLabelsType;
+ typedef itk::ConstantBoundaryCondition<TOutputLabelImage> BoundaryConditionType;
+ typedef itk::NeighborhoodIterator<TOutputLabelImage,BoundaryConditionType> NeighborhoodLabelIteratorType;
+ typedef typename NeighborhoodLabelIteratorType::NeighborhoodType NeighborhoodType;
+ typedef typename NeighborhoodType::ConstIterator NeighborhoodContainerIteratorType;
+
+ // Sets the spatial width, i.e. average root of size of the output superpixels
+ itkSetMacro(SpatialWidth,unsigned int);
+
+ // Sets the spatial distance weight, i.e. the compacity/homogeneity ratio of SP
+ itkSetMacro(SpatialDistanceWeight,RealType);
+
+ //Sets the max number of iterations
+ itkSetMacro(MaxIterationNumber,unsigned int);
+
+ //Sets the residual threshold for stopping the algorithm
+ itkSetMacro(Threshold,RealType);
+
+ //Desired number of superpixels along the X axis of the current tile
+ itkSetMacro(NbSPx, unsigned int);
+
+ //Margin in superpixels
+ itkSetMacro(Margin, unsigned int);
+
+ //Sets the input
+ using Superclass::SetInput;
+ void SetInput(const TInputImage* InputImagePtr);
+
+ //Add image of margin to list of margins
+ void AddInputMargin(typename TOutputLabelImage::Pointer InputSegPtr);
+ //Set information regarding the current tile (size, origin, neighbors)
+ void SetTile(ProcessingTile const& _t);
+ //Clear the margin images
+ void ClearMargins();
+ // //Returns the const image of region labels
+ const OutputLabelImageType * GetOutput() const;
+ // //Returns the image of region labels
+ OutputLabelImageType * GetOutput();
+ //Write the segmentation of a region
+ template <class TImage> static void WriteImageRegion(typename TImage::Pointer inputPtr, typename TImage::Pointer inputImage, typename TImage::RegionType const& region);
+
+ protected:
+ SLICFilter();
+ ~SLICFilter() ITK_OVERRIDE;
+ void GenerateData() ITK_OVERRIDE;
+ void GenerateInputRequestedRegion() ITK_OVERRIDE;
+ TInputImage* GetInput();
+
+ private:
+ SLICFilter(const Self &); //purposely not implemented
+ void operator =(const Self&); //purposely not implemented
+ unsigned int m_SpatialWidth;
+ RealType m_SpatialDistanceWeight;
+ unsigned int m_MaxIterationNumber;
+ RealType m_Threshold;
+ ProcessingTile m_Tile;
+ unsigned int m_NbSPx;
+ unsigned int m_Margin;
+ std::vector<typename OutputLabelImageType::Pointer> m_InputSeg;
+ //Test to check if a pixel is in the corner of the tile (with margins)
+ static bool inImageCorner(std::vector<RegionType> const& corners, typename TOutputLabelImage::IndexType index);
+ //Method to recompute the new coordinates of a centroid when removing a value
+ static void removeFromCentroid(PixelType const& currentPixel, typename TInputImage::IndexType index, LabelType label, CountContainerType & count, ListSampleType & means);
+ //Method to recompute the new coordinates of a centroid when adding a value
+ static void addToCentroid(PixelType const& currentPixel, typename TInputImage::IndexType index, LabelType label, CountContainerType & count, ListSampleType & means);
+
+ };
+}
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSLICFilter.txx"
+#endif
+
+#endif
diff --git a/include/otbSLICFilter.h~ b/include/otbSLICFilter.h~
new file mode 100644
index 0000000000000000000000000000000000000000..377e748f021a029b2f81e9870bf6f439eb53aefd
--- /dev/null
+++ b/include/otbSLICFilter.h~
@@ -0,0 +1,148 @@
+#ifndef otbSLICFilter_h
+#define otbSLICFilter_h
+
+#include "itkMacro.h"
+//OTB, ITK includes
+#include "otbVectorImage.h"
+#include "otbImage.h"
+#include "otbStreamingTraits.h"
+
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkImageRegionExclusionConstIteratorWithIndex.h"
+#include "itkNeighborhoodIterator.h"
+
+#include "itkProcessObject.h"
+//Own includes
+#include "simplePointCalculator.h"
+
+//STD includes
+#include <algorithm>//std::find
+#include <unordered_set>
+#include <unordered_map>
+namespace otb{
+ enum RelativePosition{
+ POS_TOP,
+ POS_RIGHT,
+ POS_BOTTOM,
+ POS_LEFT
+ };
+
+ enum NeighborhoodRelativePosition{
+ NBH_TOP,
+ NBH_TOP_RIGHT,
+ NBH_RIGHT,
+ NBH_BOTTOM_RIGHT,
+ NBH_BOTTOM,
+ NBH_BOTTOM_LEFT,
+ NBH_LEFT,
+ NBH_TOP_LEFT
+ };
+ struct ProcessingTile
+ {
+ long int rows[2]; // lower and upper rows (-1 means that the row has not be considered)
+ long int columns[2]; // lower and upper columns (-1 means that the row has not be considered)
+ long int tileNeighbors[8]; // tile Neighbors at (top, top right, right, bottom right, bottom, bottom left, left, top left)
+ long int margin[4]; // Is there a margin at top, left, bottom or right
+ otb::VectorImage<double>::RegionType region; // The image region
+
+ // Temporary files
+ std::string nodeFileName;
+ std::string edgeFileName;
+ std::string nodeMarginFileName;
+ std::string edgeMarginFileName;
+ };
+ template <class TInputImage, class TOutputLabelImage>
+ class ITK_EXPORT SLICFilter : public itk::ProcessObject
+ {
+ public:
+ typedef SLICFilter Self;
+ typedef itk::ProcessObject Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Runtime information support. */
+ itkTypeMacro(SLICFilter, itk::ProcessObject);
+
+ /** ImageDimension constants */
+ itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension);
+
+ typedef double RealType;
+ typedef typename TOutputLabelImage::PixelType LabelType;
+ typedef TOutputLabelImage OutputLabelImageType;
+ typedef TInputImage VectorImageType;
+ typedef typename TInputImage::RegionType RegionType;
+ typedef typename std::vector<double> SampleType;
+ typedef typename VectorImage<double>::PixelType PixelType;
+ typedef std::unordered_map<LabelType,SampleType> ListSampleType;
+ typedef std::unordered_map<LabelType,unsigned int> CountContainerType;
+ typedef std::vector<LabelType> NeighborhoodLabelsType;
+ typedef itk::ConstantBoundaryCondition<TOutputLabelImage> BoundaryConditionType;
+ typedef itk::NeighborhoodIterator<TOutputLabelImage,BoundaryConditionType> NeighborhoodLabelIteratorType;
+ typedef typename NeighborhoodLabelIteratorType::NeighborhoodType NeighborhoodType;
+ typedef typename NeighborhoodType::ConstIterator NeighborhoodContainerIteratorType;
+
+ // Sets the spatial width, i.e. average root of size of the output superpixels
+ itkSetMacro(SpatialWidth,unsigned int);
+
+ // Sets the spatial distance weight, i.e. the compacity/homogeneity ratio of SP
+ itkSetMacro(SpatialDistanceWeight,RealType);
+
+ //Sets the max number of iterations
+ itkSetMacro(MaxIterationNumber,unsigned int);
+
+ //Sets the residual threshold for stopping the algorithm
+ itkSetMacro(Threshold,RealType);
+
+ itkSetMacro(NbSPx, unsigned int);
+
+ itkSetMacro(Margin, unsigned int);
+
+ using Superclass::SetInput;
+ void SetInput(const TInputImage* InputImagePtr);
+ void AddInputMargin(typename TOutputLabelImage::Pointer InputSegPtr);
+
+ void SetTile(ProcessingTile const& _t);
+ void ClearMargins();
+ // //Returns the const image of region labels
+ const OutputLabelImageType * GetOutput() const;
+ // //Returns the image of region labels
+ OutputLabelImageType * GetOutput();
+ template <class TImage> static void WriteImageRegion(typename TImage::Pointer inputPtr, typename TImage::Pointer inputImage, typename TImage::RegionType const& region);
+
+ protected:
+ SLICFilter();
+
+ ~SLICFilter() ITK_OVERRIDE;
+
+ void GenerateData() ITK_OVERRIDE;
+
+ void GenerateInputRequestedRegion() ITK_OVERRIDE;
+ // void PrintLogLikelihood();
+ TInputImage* GetInput();
+ private:
+ SLICFilter(const Self &); //purposely not implemented
+ void operator =(const Self&); //purposely not implemented
+ unsigned int m_SpatialWidth;
+ RealType m_SpatialDistanceWeight;
+ unsigned int m_MaxIterationNumber;
+ RealType m_Threshold;
+ ProcessingTile m_Tile;
+ unsigned int m_NbSPx;
+ unsigned int m_Margin;
+ std::vector<typename OutputLabelImageType::Pointer> m_InputSeg;
+
+ static bool inImageCorner(std::vector<RegionType> const& corners, typename TOutputLabelImage::IndexType index);
+ static void removeFromCentroid(PixelType const& currentPixel, typename TInputImage::IndexType index, LabelType label, CountContainerType & count, ListSampleType & means);
+ static void addToCentroid(PixelType const& currentPixel, typename TInputImage::IndexType index, LabelType label, CountContainerType & count, ListSampleType & means);
+
+ };
+}
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSLICFilter.txx"
+#endif
+
+#endif
diff --git a/include/otbSLICFilter.txx b/include/otbSLICFilter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..a33f94f0b77fea73e95a47005568986d9af4db4e
--- /dev/null
+++ b/include/otbSLICFilter.txx
@@ -0,0 +1,383 @@
+#ifndef otbSLICFilter_txx
+#define otbSLICFilter_txx
+
+#include "otbSLICFilter.h"
+
+namespace otb{
+
+ template <class TInputImage, class TOutputLabelImage>
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::SLICFilter() :
+ m_SpatialWidth(100),
+ m_SpatialDistanceWeight(1),
+ m_MaxIterationNumber(10),
+ m_Threshold(0.1),
+ m_NbSPx(1),
+ m_Margin(0)
+ {
+ this->SetNumberOfRequiredInputs(1);
+ this->SetNumberOfRequiredOutputs(1);
+ this->SetNthOutput(0,TOutputLabelImage::New());
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::~SLICFilter()
+ {}
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::SetInput(const TInputImage* input)
+ {
+ this->ProcessObject::SetNthInput(0,const_cast<TInputImage*>(input));
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::AddInputMargin(typename TOutputLabelImage::Pointer InputSegPtr)
+ {
+ m_InputSeg.push_back(InputSegPtr);
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::SetTile(ProcessingTile const& _t)
+ {
+ m_Tile = _t;
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::ClearMargins()
+ {
+ m_InputSeg.clear();
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ TInputImage *
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GetInput()
+ {
+ return static_cast<TInputImage *>(this->itk::ProcessObject::GetInput(0));
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ TOutputLabelImage *
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GetOutput()
+ {
+ return static_cast<TOutputLabelImage *>(this->itk::ProcessObject::GetOutput(0));
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ const TOutputLabelImage *
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GetOutput() const
+ {
+ return static_cast<TOutputLabelImage *>(this->itk::ProcessObject::GetOutput(0));
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GenerateInputRequestedRegion(){
+ this->GetOutput()->SetRegions(this->GetInput()->GetLargestPossibleRegion());
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ bool
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::inImageCorner(std::vector<RegionType> const& corners, typename TOutputLabelImage::IndexType index){
+ for(RegionType cornerRegion : corners)
+ if(cornerRegion.IsInside(index))
+ return true;
+ return false;
+ }
+
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::removeFromCentroid(PixelType const& currentPixel, typename TInputImage::IndexType index, LabelType label, CountContainerType & count, ListSampleType & means){
+ const typename CountContainerType::iterator searchCount = count.find(label);
+ //Update with running formulas
+ //Update count
+ searchCount->second--;
+ //Update Means
+ const typename ListSampleType::iterator searchM = means.find(label);
+ const SampleType oldMeans = searchM->second;
+ SampleType newMeans = oldMeans;
+ const double N = searchCount->second;//Cast int to double
+ const unsigned int nbComps = newMeans.size() - 2 ;
+ for(unsigned int i = 0 ; i < nbComps ; ++i)
+ {
+ newMeans[i] += (oldMeans[i]-currentPixel[i])/N;
+ }
+ newMeans[nbComps] += (oldMeans[nbComps]-index[0])/N;
+ newMeans[nbComps+1] += (oldMeans[nbComps+1]-index[1])/N;
+ searchM->second = newMeans;
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::addToCentroid(PixelType const& currentPixel, typename TInputImage::IndexType index, LabelType label, CountContainerType & count, ListSampleType & means){
+ const typename CountContainerType::iterator searchCount = count.find(label);
+ //Update with running formulas
+ //Update count
+ searchCount->second++;
+ //Update Means
+ const typename ListSampleType::iterator searchM = means.find(label);
+ const SampleType oldMeans = searchM->second;
+ SampleType newMeans = oldMeans;
+ const double N = searchCount->second;//Cast int to double
+ const unsigned int nbComps = newMeans.size() - 2 ;
+ for(unsigned int i = 0 ; i < nbComps ; ++i)
+ {
+ newMeans[i] += (currentPixel[i]-oldMeans[i])/N;
+ }
+ newMeans[nbComps] += (index[0]-oldMeans[nbComps])/N;
+ newMeans[nbComps+1] += (index[1]-oldMeans[nbComps+1])/N;
+ searchM->second = newMeans;
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GenerateData()
+ {
+ this->GetOutput()->SetRegions(this->GetInput()->GetLargestPossibleRegion());
+ this->GetOutput()->SetOrigin(this->GetInput()->GetOrigin());
+ this->GetOutput()->SetSignedSpacing(this->GetInput()->GetSignedSpacing());
+ this->GetOutput()->Allocate();
+ this->GetOutput()->FillBuffer(0);
+
+ //Initialize region without margins
+ const typename RegionType::IndexType innerIndex({m_Tile.region.GetIndex()[0]+ m_Tile.margin[3],m_Tile.region.GetIndex()[1]+m_Tile.margin[0]});
+ const typename RegionType::SizeType innerSize({m_Tile.region.GetSize()[0]-m_Tile.margin[3]-m_Tile.margin[1],m_Tile.region.GetSize()[1]-m_Tile.margin[0]-m_Tile.margin[2]});
+ const RegionType innerRegion(innerIndex,innerSize);
+ const unsigned int nbComps = this->GetInput()->GetNumberOfComponentsPerPixel();
+ std::unordered_set<LabelType> stableLabels;
+ //Calculate stable segments, depending on whether tile is touching image edge
+ const unsigned int topMargin = m_Tile.tileNeighbors[NBH_TOP] >= 0 ? 1 : 0;
+ const unsigned int rightMargin = m_Tile.tileNeighbors[NBH_RIGHT] >= 0 ? 1 : 0;
+ const unsigned int bottomMargin = m_Tile.tileNeighbors[NBH_BOTTOM] >= 0 ? 1 : 0;
+ const unsigned int leftMargin = m_Tile.tileNeighbors[NBH_LEFT] >= 0 ? 1 : 0;
+ const bool stableSegments = (topMargin || rightMargin || bottomMargin || leftMargin) && (m_Margin > 0);
+ //Calculate margin region (contains all pixels except for ones on the edge of the image)
+ //TODO Improve this with faces calculator ?
+ const typename RegionType::IndexType marginIndex({m_Tile.region.GetIndex()[0]+leftMargin,m_Tile.region.GetIndex()[1]+topMargin});
+ const typename RegionType::SizeType marginSize({m_Tile.region.GetSize()[0]-leftMargin-rightMargin,m_Tile.region.GetSize()[1]-topMargin-bottomMargin});
+ const RegionType marginRegion(marginIndex,marginSize);
+
+ typename RegionType::IndexType cornerIndex;
+ std::vector<RegionType> corners;
+
+ if(bottomMargin && leftMargin)
+ {
+ const long unsigned int m = m_Tile.margin[2];
+ cornerIndex[0]=m_Tile.region.GetIndex()[0];
+ cornerIndex[1]=m_Tile.region.GetIndex()[1]+m_Tile.region.GetSize()[1]-m;
+ const typename RegionType::SizeType cornerSize({m,m});
+ corners.push_back(RegionType(cornerIndex,cornerSize));
+ }
+ if(bottomMargin && rightMargin)
+ {
+ const long unsigned int m = m_Tile.margin[2];
+ cornerIndex[0]=m_Tile.region.GetIndex()[0]+m_Tile.region.GetSize()[0]-m;
+ cornerIndex[1]=m_Tile.region.GetIndex()[1]+m_Tile.region.GetSize()[1]-m;
+ const typename RegionType::SizeType cornerSize({m,m});
+ corners.push_back(RegionType(cornerIndex,cornerSize));
+ }
+
+ //Write margin values into output
+ for(auto const& im : m_InputSeg)
+ {
+ itk::ImageRegionIterator<TOutputLabelImage> outputIt(this->GetOutput(),im->GetLargestPossibleRegion());
+ itk::ImageRegionConstIterator<TOutputLabelImage> inputIt(im,im->GetLargestPossibleRegion());
+ for(inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd() && !outputIt.IsAtEnd(); ++inputIt, ++outputIt)
+ {
+ const LabelType currentLabel = inputIt.Get();
+ outputIt.Set(currentLabel);
+ }
+ }
+ ListSampleType initialSeedMeans;
+ CountContainerType initialCount;
+ // //To avoid rehashes
+ initialSeedMeans.reserve(m_Tile.region.GetSize()[0]*m_Tile.region.GetSize()[1]/m_SpatialWidth/m_SpatialWidth);
+ initialCount.reserve(m_Tile.region.GetSize()[0]*m_Tile.region.GetSize()[1]/m_SpatialWidth/m_SpatialWidth);
+ //Update means and count in margin zones according to initial grid and stable elements
+ itk::ImageRegionIterator<TOutputLabelImage> outputIt(this->GetOutput(),this->GetOutput()->GetLargestPossibleRegion());
+ itk::ImageRegionConstIteratorWithIndex<TInputImage> vit(this->GetInput(), this->GetInput()->GetLargestPossibleRegion());
+ for(vit.GoToBegin(), outputIt.GoToBegin(); !outputIt.IsAtEnd() && !vit.IsAtEnd() ; ++outputIt, ++vit)
+ {
+ if(!inImageCorner(corners,vit.GetIndex())){
+ LabelType currentLabel;
+ if(innerRegion.IsInside(vit.GetIndex())){
+ //Set grid in inner region
+ currentLabel = vcl_floor(vit.GetIndex()[0]/m_SpatialWidth)+m_NbSPx*vcl_floor(vit.GetIndex()[1]/m_SpatialWidth)+1;
+ outputIt.Set(currentLabel);
+ }
+ else {
+ //Read from margin values
+ currentLabel = outputIt.Get();
+ }
+ //Calculate initial means and count
+ if(initialSeedMeans.find(currentLabel) == initialSeedMeans.end()){
+ initialCount.insert({currentLabel,1});
+
+ SampleType initialSeed(nbComps+2);
+ for (unsigned int i = 0; i < nbComps; i++) {
+ initialSeed[i] = vit.Get()[i];
+ }
+ initialSeed[nbComps] = vit.GetIndex()[0];
+ initialSeed[nbComps+1] = vit.GetIndex()[1];
+
+ initialSeedMeans.insert({currentLabel,initialSeed});
+ }
+ else {
+ addToCentroid(vit.Get(),vit.GetIndex(),currentLabel,initialCount,initialSeedMeans);
+ }
+ //Calculate stable segments
+ if(stableSegments && !marginRegion.IsInside(vit.GetIndex()) && stableLabels.find(currentLabel) == stableLabels.end()){
+ stableLabels.insert(currentLabel);
+ }
+ }
+ }
+ //Setup first iteration
+ double residual = 0;
+ bool firstLoop = (m_MaxIterationNumber > 0);
+ unsigned int iteration = 0;
+ //Copy initial values
+ ListSampleType currentSeedMeans(initialSeedMeans);
+ CountContainerType currentCount(initialCount);
+ simplePointCalculator<NeighborhoodContainerIteratorType,LabelType> spc;
+
+ //Take neighborhood iterator on labels
+ const typename NeighborhoodLabelIteratorType::RadiusType radius = {{1,1}};
+ NeighborhoodLabelIteratorType labIt(radius,this->GetOutput(),m_Tile.region);
+ itk::ImageRegionConstIteratorWithIndex<VectorImageType> vitRegion(this->GetInput(),m_Tile.region);
+
+ while(firstLoop || (residual > m_Threshold && iteration < m_MaxIterationNumber))
+ {
+ std::cout<<"Iteration "<<iteration<<", starting optimisation ..."<<std::endl;
+ for(vitRegion.GoToBegin(),labIt.GoToBegin();!vitRegion.IsAtEnd() && !labIt.IsAtEnd();++vitRegion,++labIt)
+ {
+
+ if(!inImageCorner(corners,vitRegion.GetIndex()))
+ {
+ const LabelType currentLabel = labIt.GetCenterPixel();
+ if(currentCount[currentLabel] > 1)
+ {
+ if(!stableSegments || stableLabels.find(currentLabel) == stableLabels.end())
+ {
+ //Copy 4 connected unique labels into std::vector
+ NeighborhoodLabelsType neighborLabels;
+ const NeighborhoodType n = labIt.GetNeighborhood();
+ typename NeighborhoodType::ConstIterator ni;
+ for(ni = n.Begin()+1 ;ni != n.End(); ni=ni+2)
+ {
+
+ if(*ni != 0 && (!stableSegments || stableLabels.find(*ni) == stableLabels.end()))
+ {
+ if(std::find(neighborLabels.begin(), neighborLabels.end(), *ni) == neighborLabels.end())
+ {
+
+ neighborLabels.push_back(*ni);
+ }
+ }
+ }
+ int iMin=-1;
+ double distanceMin = std::numeric_limits<double>::max();
+ if(neighborLabels.size() > 1)
+ {
+ bool topo_check = true;
+ //Point must be simple wrt all the 4 connected labeling neighbors
+ const unsigned int neighborLabelsSize = neighborLabels.size();
+ for(unsigned int i = 0; i < neighborLabelsSize; i++)
+ {
+ topo_check&=spc.isSimpleLabel(n.Begin(), n.End(), neighborLabels[i]);
+ }
+ if(topo_check)
+ {
+ //Possible label change
+
+ for(unsigned int i = 0; i < neighborLabelsSize; i++)
+ {
+ //Find closest label
+ double spectralDistance = 0;
+ double spatialDistance = 0;
+ const SampleType currentSeed = currentSeedMeans[neighborLabels[i]];
+ for(unsigned int comp = 0; comp<nbComps;++comp)
+ {
+ spectralDistance+=(vitRegion.Get()[comp]-currentSeed[comp])*(vitRegion.Get()[comp]-currentSeed[comp]);
+ }
+ spatialDistance += (vitRegion.GetIndex()[0]-currentSeed[nbComps])*(vitRegion.GetIndex()[0]-currentSeed[nbComps]);
+ spatialDistance += (vitRegion.GetIndex()[1]-currentSeed[nbComps+1])*(vitRegion.GetIndex()[1]-currentSeed[nbComps+1]);
+ const double distance = vcl_sqrt(spectralDistance)/nbComps + m_SpatialDistanceWeight * vcl_sqrt(spatialDistance)/2;
+ if (distance < distanceMin)
+ {
+ iMin=i;
+ distanceMin=distance;
+ }
+ }
+ }
+ }
+ if (iMin >= 0 && neighborLabels[iMin] != currentLabel)
+ {
+ //Change label
+ const LabelType newLabel = neighborLabels[iMin];
+ labIt.SetCenterPixel(newLabel);
+
+ //Update seeds
+ //Remove point at currentlabel from means
+ removeFromCentroid(vitRegion.Get(),vitRegion.GetIndex(),currentLabel,currentCount,currentSeedMeans);
+ //neighborLabel[iMin]
+ //Add point at newLabel in means
+ addToCentroid(vitRegion.Get(),vitRegion.GetIndex(),newLabel,currentCount,currentSeedMeans);
+ }
+ }
+ }
+ }
+ }
+ //Compute residual
+ residual = 0;
+ //Update initialSeedMeans
+ typename ListSampleType::const_iterator newSeed = currentSeedMeans.begin();
+ typename ListSampleType::iterator oldSeed = initialSeedMeans.begin();
+ for(;newSeed != currentSeedMeans.end() && oldSeed != initialSeedMeans.end();++newSeed,++oldSeed)
+ {
+ //spectral residual
+ for(unsigned int i = 0; i < nbComps ;++i)
+ {
+ residual+=vcl_abs(newSeed->second[i] - oldSeed->second[i])/nbComps;
+ }
+ //spatial residual
+ for(unsigned int i = nbComps; i < nbComps+2 ;++i)
+ {
+ residual+=vcl_abs(newSeed->second[i] - oldSeed->second[i])/2;
+ }
+
+ //Update initialSeedMeans
+ oldSeed->second=newSeed->second;
+ }
+ residual/=currentSeedMeans.size();
+ std::cout<<"Done. Residual: "<<residual<<", max remaining iterations: "<<m_MaxIterationNumber-iteration<<std::endl;
+
+ //Update count
+ typename CountContainerType::const_iterator newSeedCount = currentCount.begin();
+ typename CountContainerType::iterator oldSeedCount = initialCount.begin();
+ for(;newSeedCount != currentCount.end() && oldSeedCount != initialCount.end(); ++newSeedCount,++oldSeedCount)
+ {
+ oldSeedCount->second = newSeedCount->second;
+ }
+ firstLoop = false;
+ ++iteration;
+
+ }
+
+ return;
+ }
+
+
+} //end namespace otb
+#endif
+
diff --git a/include/otbSLICFilter.txx~ b/include/otbSLICFilter.txx~
new file mode 100644
index 0000000000000000000000000000000000000000..fe2e1d8d5e0e8cfbd165864e53d806360c2af907
--- /dev/null
+++ b/include/otbSLICFilter.txx~
@@ -0,0 +1,383 @@
+#ifndef otbSLICFilter_txx
+#define otbSLICFilter_txx
+
+#include "otbSLICFilter.h"
+
+namespace otb{
+
+ template <class TInputImage, class TOutputLabelImage>
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::SLICFilter() :
+ m_SpatialWidth(100),
+ m_SpatialDistanceWeight(1),
+ m_MaxIterationNumber(10),
+ m_Threshold(0.1),
+ m_NbSPx(1),
+ m_Margin(0)
+ {
+ this->SetNumberOfRequiredInputs(1);
+ this->SetNumberOfRequiredOutputs(1);
+ this->SetNthOutput(0,TOutputLabelImage::New());
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::~SLICFilter()
+ {}
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::SetInput(const TInputImage* input)
+ {
+ this->ProcessObject::SetNthInput(0,const_cast<TInputImage*>(input));
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::AddInputMargin(typename TOutputLabelImage::Pointer InputSegPtr)
+ {
+ m_InputSeg.push_back(InputSegPtr);
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::SetTile(ProcessingTile const& _t)
+ {
+ m_Tile = _t;
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::ClearMargins()
+ {
+ m_InputSeg.clear();
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ TInputImage *
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GetInput()
+ {
+ return static_cast<TInputImage *>(this->itk::ProcessObject::GetInput(0));
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ TOutputLabelImage *
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GetOutput()
+ {
+ return static_cast<TOutputLabelImage *>(this->itk::ProcessObject::GetOutput(0));
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ const TOutputLabelImage *
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GetOutput() const
+ {
+ return static_cast<TOutputLabelImage *>(this->itk::ProcessObject::GetOutput(0));
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GenerateInputRequestedRegion(){
+ this->GetOutput()->SetRegions(this->GetInput()->GetLargestPossibleRegion());
+ }
+ template <class TInputImage, class TOutputLabelImage>
+ bool
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::inImageCorner(std::vector<RegionType> const& corners, typename TOutputLabelImage::IndexType index){
+ for(RegionType cornerRegion : corners)
+ if(cornerRegion.IsInside(index))
+ return true;
+ return false;
+ }
+
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::removeFromCentroid(PixelType const& currentPixel, typename TInputImage::IndexType index, LabelType label, CountContainerType & count, ListSampleType & means){
+ const typename CountContainerType::iterator searchCount = count.find(label);
+ //Update with running formulas
+ //Update count
+ searchCount->second--;
+ //Update Means
+ const typename ListSampleType::iterator searchM = means.find(label);
+ const SampleType oldMeans = searchM->second;
+ SampleType newMeans = oldMeans;
+ const double N = searchCount->second;//Cast int to double
+ const unsigned int nbComps = newMeans.size() - 2 ;
+ for(unsigned int i = 0 ; i < nbComps ; ++i)
+ {
+ newMeans[i] += (oldMeans[i]-currentPixel[i])/N;
+ }
+ newMeans[nbComps] += (oldMeans[nbComps]-index[0])/N;
+ newMeans[nbComps+1] += (oldMeans[nbComps+1]-index[1])/N;
+ searchM->second = newMeans;
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::addToCentroid(PixelType const& currentPixel, typename TInputImage::IndexType index, LabelType label, CountContainerType & count, ListSampleType & means){
+ const typename CountContainerType::iterator searchCount = count.find(label);
+ //Update with running formulas
+ //Update count
+ searchCount->second++;
+ //Update Means
+ const typename ListSampleType::iterator searchM = means.find(label);
+ const SampleType oldMeans = searchM->second;
+ SampleType newMeans = oldMeans;
+ const double N = searchCount->second;//Cast int to double
+ const unsigned int nbComps = newMeans.size() - 2 ;
+ for(unsigned int i = 0 ; i < nbComps ; ++i)
+ {
+ newMeans[i] += (currentPixel[i]-oldMeans[i])/N;
+ }
+ newMeans[nbComps] += (index[0]-oldMeans[nbComps])/N;
+ newMeans[nbComps+1] += (index[1]-oldMeans[nbComps+1])/N;
+ searchM->second = newMeans;
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICFilter<TInputImage, TOutputLabelImage>
+ ::GenerateData()
+ {
+ this->GetOutput()->SetRegions(this->GetInput()->GetLargestPossibleRegion());
+ this->GetOutput()->SetOrigin(this->GetInput()->GetOrigin());
+ this->GetOutput()->SetSignedSpacing(this->GetInput()->GetSignedSpacing());
+ this->GetOutput()->Allocate();
+ this->GetOutput()->FillBuffer(0);
+
+ //Initialize region without margins
+ const typename RegionType::IndexType innerIndex({m_Tile.region.GetIndex()[0]+ m_Tile.margin[3],m_Tile.region.GetIndex()[1]+m_Tile.margin[0]});
+ const typename RegionType::SizeType innerSize({m_Tile.region.GetSize()[0]-m_Tile.margin[3]-m_Tile.margin[1],m_Tile.region.GetSize()[1]-m_Tile.margin[0]-m_Tile.margin[2]});
+ const RegionType innerRegion(innerIndex,innerSize);
+ const unsigned int nbComps = this->GetInput()->GetNumberOfComponentsPerPixel();
+ std::unordered_set<LabelType> stableLabels;
+ //Calculate stable segments, depending on whether tile is touching image edge
+ const unsigned int topMargin = m_Tile.tileNeighbors[NBH_TOP] >= 0 ? 1 : 0;
+ const unsigned int rightMargin = m_Tile.tileNeighbors[NBH_RIGHT] >= 0 ? 1 : 0;
+ const unsigned int bottomMargin = m_Tile.tileNeighbors[NBH_BOTTOM] >= 0 ? 1 : 0;
+ const unsigned int leftMargin = m_Tile.tileNeighbors[NBH_LEFT] >= 0 ? 1 : 0;
+ const bool stableSegments = (topMargin || rightMargin || bottomMargin || leftMargin) && (m_Margin > 0);
+ //Calculate margin region (contains all pixels except for ones on the edge of the image)
+ //TODO Improve this with faces calculator ?
+ const typename RegionType::IndexType marginIndex({m_Tile.region.GetIndex()[0]+leftMargin,m_Tile.region.GetIndex()[1]+topMargin});
+ const typename RegionType::SizeType marginSize({m_Tile.region.GetSize()[0]-leftMargin-rightMargin,m_Tile.region.GetSize()[1]-topMargin-bottomMargin});
+ const RegionType marginRegion(marginIndex,marginSize);
+
+ typename RegionType::IndexType cornerIndex;
+ std::vector<RegionType> corners;
+
+ if(bottomMargin && leftMargin)
+ {
+ const long unsigned int m = m_Tile.margin[2];
+ cornerIndex[0]=m_Tile.region.GetIndex()[0];
+ cornerIndex[1]=m_Tile.region.GetIndex()[1]+m_Tile.region.GetSize()[1]-m;
+ const typename RegionType::SizeType cornerSize({m,m});
+ corners.push_back(RegionType(cornerIndex,cornerSize));
+ }
+ if(bottomMargin && rightMargin)
+ {
+ const long unsigned int m = m_Tile.margin[2];
+ cornerIndex[0]=m_Tile.region.GetIndex()[0]+m_Tile.region.GetSize()[0]-m;
+ cornerIndex[1]=m_Tile.region.GetIndex()[1]+m_Tile.region.GetSize()[1]-m;
+ const typename RegionType::SizeType cornerSize({m,m});
+ corners.push_back(RegionType(cornerIndex,cornerSize));
+ }
+
+ //Write margin values into output
+ for(auto const& im : m_InputSeg)
+ {
+ itk::ImageRegionIterator<TOutputLabelImage> outputIt(this->GetOutput(),im->GetLargestPossibleRegion());
+ itk::ImageRegionConstIterator<TOutputLabelImage> inputIt(im,im->GetLargestPossibleRegion());
+ for(inputIt.GoToBegin(), outputIt.GoToBegin(); !inputIt.IsAtEnd() && !outputIt.IsAtEnd(); ++inputIt, ++outputIt)
+ {
+ const LabelType currentLabel = inputIt.Get();
+ outputIt.Set(currentLabel);
+ }
+ }
+ ListSampleType initialSeedMeans;
+ CountContainerType initialCount;
+ // //To avoid rehashes
+ initialSeedMeans.reserve(m_Tile.region.GetSize()[0]*m_Tile.region.GetSize()[1]/m_SpatialWidth/m_SpatialWidth);
+ initialCount.reserve(m_Tile.region.GetSize()[0]*m_Tile.region.GetSize()[1]/m_SpatialWidth/m_SpatialWidth);
+ //Update means and count in margin zones according to initial grid and stable elements
+ itk::ImageRegionIterator<TOutputLabelImage> outputIt(this->GetOutput(),this->GetOutput()->GetLargestPossibleRegion());
+ itk::ImageRegionConstIteratorWithIndex<TInputImage> vit(this->GetInput(), this->GetInput()->GetLargestPossibleRegion());
+ for(vit.GoToBegin(), outputIt.GoToBegin(); !outputIt.IsAtEnd() && !vit.IsAtEnd() ; ++outputIt, ++vit)
+ {
+ if(!inImageCorner(corners,vit.GetIndex())){
+ LabelType currentLabel;
+ if(innerRegion.IsInside(vit.GetIndex())){
+ //Set grid in inner region
+ currentLabel = vcl_floor(vit.GetIndex()[0]/m_SpatialWidth)+m_NbSPx*vcl_floor(vit.GetIndex()[1]/m_SpatialWidth)+1;
+ outputIt.Set(currentLabel);
+ }
+ else {
+ //Read from margin values
+ currentLabel = outputIt.Get();
+ }
+ //Calculate initial means and count
+ if(initialSeedMeans.find(currentLabel) == initialSeedMeans.end()){
+ initialCount.insert({currentLabel,1});
+
+ SampleType initialSeed(nbComps+2);
+ for (unsigned int i = 0; i < nbComps; i++) {
+ initialSeed[i] = vit.Get()[i];
+ }
+ initialSeed[nbComps] = vit.GetIndex()[0];
+ initialSeed[nbComps+1] = vit.GetIndex()[1];
+
+ initialSeedMeans.insert({currentLabel,initialSeed});
+ }
+ else {
+ addToCentroid(vit.Get(),vit.GetIndex(),currentLabel,initialCount,initialSeedMeans);
+ }
+ //Calculate stable segments
+ if(stableSegments && !marginRegion.IsInside(vit.GetIndex()) && stableLabels.find(currentLabel) == stableLabels.end()){
+ stableLabels.insert(currentLabel);
+ }
+ }
+ }
+ //Setup first iteration
+ double residual = 0;
+ bool firstLoop = (m_MaxIterationNumber > 0);
+ unsigned int iteration = 0;
+ //Copy initial values
+ ListSampleType currentSeedMeans(initialSeedMeans);
+ CountContainerType currentCount(initialCount);
+ simplePointCalculator<NeighborhoodContainerIteratorType,LabelType> spc;
+
+ //Take neighborhood iterator on labels
+ const typename NeighborhoodLabelIteratorType::RadiusType radius = {{1,1}};
+ NeighborhoodLabelIteratorType labIt(radius,this->GetOutput(),m_Tile.region);
+ itk::ImageRegionConstIteratorWithIndex<VectorImageType> vitRegion(this->GetInput(),m_Tile.region);
+
+ while(firstLoop || (residual > m_Threshold && iteration < m_MaxIterationNumber))
+ {
+ std::cout<<"Iteration "<<iteration<<", starting optimisation ..."<<std::endl;
+ for(vitRegion.GoToBegin(),labIt.GoToBegin();!vitRegion.IsAtEnd() && !labIt.IsAtEnd();++vitRegion,++labIt)
+ {
+
+ if(!inImageCorner(corners,vitRegion.GetIndex()))
+ {
+ const LabelType currentLabel = labIt.GetCenterPixel();
+ if(currentCount[currentLabel] > 1)
+ {
+ if(!stableSegments || stableLabels.find(currentLabel) == stableLabels.end())
+ {
+ //Copy 4 connected unique labels into std::vector
+ NeighborhoodLabelsType neighborLabels;
+ const NeighborhoodType n = labIt.GetNeighborhood();
+ typename NeighborhoodType::ConstIterator ni;
+ for(ni = n.Begin()+1 ;ni != n.End(); ni=ni+2)
+ {
+
+ if(*ni != 0 && (!stableSegments || stableLabels.find(*ni) == stableLabels.end()))
+ {
+ if(std::find(neighborLabels.begin(), neighborLabels.end(), *ni) == neighborLabels.end())
+ {
+
+ neighborLabels.push_back(*ni);
+ }
+ }
+ }
+ int iMin=-1;
+ double distanceMin = std::numeric_limits<double>::max();
+ if(neighborLabels.size() > 1)
+ {
+ bool topo_check = true;
+ //Point must be simple wrt all the 4 connected labeling neighbors
+ const unsigned int neighborLabelsSize = neighborLabels.size();
+ for(unsigned int i = 0; i < neighborLabelsSize; i++)
+ {
+ topo_check&=spc.isSimpleLabel(n.Begin(), n.End(), neighborLabels[i]);
+ }
+ if(topo_check)
+ {
+ //Possible label change
+
+ for(unsigned int i = 0; i < neighborLabelsSize; i++)
+ {
+ //Find closest label
+ double spectralDistance = 0;
+ double spatialDistance = 0;
+ const SampleType currentSeed = currentSeedMeans[neighborLabels[i]];
+ for(unsigned int comp = 0; comp<nbComps;++comp)
+ {
+ spectralDistance+=(vitRegion.Get()[comp]-currentSeed[comp])*(vitRegion.Get()[comp]-currentSeed[comp]);
+ }
+ spatialDistance += (vitRegion.GetIndex()[0]-currentSeed[nbComps])*(vitRegion.GetIndex()[0]-currentSeed[nbComps]);
+ spatialDistance += (vitRegion.GetIndex()[1]-currentSeed[nbComps+1])*(vitRegion.GetIndex()[1]-currentSeed[nbComps+1]);
+ const double distance = vcl_sqrt(spectralDistance)/nbComps + m_SpatialDistanceWeight * vcl_sqrt(spatialDistance)/2;
+ if (distance < distanceMin)
+ {
+ iMin=i;
+ distanceMin=distance;
+ }
+ }
+ }
+ }
+ if (iMin >= 0 && neighborLabels[iMin] != currentLabel)
+ {
+ //Change label
+ const LabelType newLabel = neighborLabels[iMin];
+ labIt.SetCenterPixel(newLabel);
+
+ //Update seeds
+ //Remove point at currentlabel from means
+ removeFromCentroid(vitRegion.Get(),vitRegion.GetIndex(),currentLabel,currentCount,currentSeedMeans);
+ //neighborLabel[iMin]
+ //Add point at newLabel in means
+ addToCentroid(vitRegion.Get(),vitRegion.GetIndex(),newLabel,currentCount,currentSeedMeans);
+ }
+ }
+ }
+ }
+ }
+
+ residual = 0;
+ //Update initialSeedMeans
+ typename ListSampleType::const_iterator newSeed = currentSeedMeans.begin();
+ typename ListSampleType::iterator oldSeed = initialSeedMeans.begin();
+ for(;newSeed != currentSeedMeans.end() && oldSeed != initialSeedMeans.end();++newSeed,++oldSeed)
+ {
+ //spectral residual
+ for(unsigned int i = 0; i < nbComps ;++i)
+ {
+ residual+=vcl_abs(newSeed->second[i] - oldSeed->second[i])/nbComps;
+ }
+ //spatial residual
+ for(unsigned int i = nbComps; i < nbComps+2 ;++i)
+ {
+ residual+=vcl_abs(newSeed->second[i] - oldSeed->second[i])/2;
+ }
+
+ //Update initialSeedMeans
+ oldSeed->second=newSeed->second;
+ }
+ residual/=currentSeedMeans.size();
+ std::cout<<"Done. Residual: "<<residual<<", max remaining iterations: "<<m_MaxIterationNumber-iteration<<std::endl;
+
+ //Update count
+ typename CountContainerType::const_iterator newSeedCount = currentCount.begin();
+ typename CountContainerType::iterator oldSeedCount = initialCount.begin();
+ for(;newSeedCount != currentCount.end() && oldSeedCount != initialCount.end(); ++newSeedCount,++oldSeedCount)
+ {
+ oldSeedCount->second = newSeedCount->second;
+ }
+ firstLoop = false;
+ ++iteration;
+
+ }
+
+ return;
+ }
+
+
+} //end namespace otb
+#endif
+
diff --git a/include/otbSLICScheduler.h b/include/otbSLICScheduler.h
new file mode 100644
index 0000000000000000000000000000000000000000..0220507982cf19d458175a594209ee2830f2c677
--- /dev/null
+++ b/include/otbSLICScheduler.h
@@ -0,0 +1,150 @@
+#ifndef otbSLICScheduler_h
+#define otbSLICScheduler_h
+
+#include "itkMacro.h"
+//OTB, ITK includes
+#include "otbVectorImage.h"
+#include "otbImage.h"
+#include "otbImageFileWriter.h"
+#include "otbImageFileReader.h"
+
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkImageRegionExclusionConstIteratorWithIndex.h"
+
+#include "itkImageToImageFilter.h"
+#include "itkNeighborhoodIterator.h"
+#include "itkExtractImageFilter.h"
+#include "otbMPIConfig.h"
+
+//Own includes
+#include "otbSLICFilter.h"
+
+#include <itksys/SystemTools.hxx>
+
+#include <gdal.h>
+#include <gdal_priv.h>
+#if defined(__GNUC__) || defined(__clang__)
+# pragma GCC diagnostic push
+# pragma GCC diagnostic ignored "-Wshadow"
+#include <vrtdataset.h>
+# pragma GCC diagnostic pop
+#else
+#include <vrtdataset.h>
+#endif
+
+#include <ogr_spatialref.h>
+
+
+
+
+namespace otb{
+
+ template <class TInputImage, class TOutputLabelImage>
+ class ITK_EXPORT SLICScheduler : public itk::Object
+ {
+ public:
+ typedef SLICScheduler Self;
+ typedef itk::Object Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Runtime information support. */
+ itkTypeMacro(SLICScheduler, ImageToImageFilter);
+
+ /** ImageDimension constants */
+ itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension);
+
+ typedef double RealType;
+ typedef typename TOutputLabelImage::PixelType LabelType;
+ typedef otb::ImageFileReader<TInputImage> ReaderType;
+ typedef TOutputLabelImage OutputLabelImageType;
+ typedef TInputImage VectorImageType;
+ typedef typename TInputImage::RegionType RegionType;
+ typedef typename std::vector<RealType> SampleType;
+ typedef std::vector<SampleType> ListSampleType;
+ typedef std::vector<LabelType> NeighborhoodLabelsType;
+ typedef itk::ConstantBoundaryCondition<TOutputLabelImage> BoundaryConditionType;
+ typedef itk::NeighborhoodIterator<TOutputLabelImage,BoundaryConditionType> NeighborhoodLabelIteratorType;
+ typedef typename NeighborhoodLabelIteratorType::NeighborhoodType NeighborhoodType;
+ typedef typename NeighborhoodType::ConstIterator NeighborhoodContainerIteratorType;
+ typedef otb::ImageFileWriter<TOutputLabelImage> WriterType;
+ typedef std::vector<otb::ProcessingTile> TileListType;
+ typedef otb::SLICFilter<TInputImage, TOutputLabelImage> SLICFilter;
+ typedef itk::ExtractImageFilter<TInputImage,TInputImage> ExtractFilterType;
+ typedef otb::ImageFileReader<TOutputLabelImage> LabelReaderType;
+ typedef itk::ExtractImageFilter<TOutputLabelImage, TOutputLabelImage> LabelExtractFilterType;
+
+ // Sets the spatial width, i.e. average root of size of the output superpixels
+ itkSetMacro(SpatialWidth,unsigned int);
+
+ // Sets the spatial distance weight, i.e. the compacity/homogeneity ratio of SP
+ itkSetMacro(SpatialDistanceWeight,RealType);
+
+ //Sets the max number of iterations
+ itkSetMacro(MaxIterationNumber,unsigned int);
+
+ //Sets the residual threshold for stopping the algorithm
+ itkSetMacro(Threshold,RealType);
+
+ //Number of tiles in X direction
+ itkSetMacro(NbTilesX,unsigned int);
+
+ //Number of tiles in Y direction
+ itkSetMacro(NbTilesY,unsigned int);
+
+ //Margin (in SP widths)
+ itkSetMacro(Margin,RealType);
+
+ //Prefix for file writing
+ itkSetMacro(Prefix,std::string);
+
+ //Name of the input for writing temp files
+ itkSetMacro(InputName,std::string);
+
+ //Calculate a list of regions of the input image with the appropriate margins
+ static std::vector<ProcessingTile> SplitOTBImage(const TInputImage * imagePtr, // input image
+ const unsigned int tileWidth, // width of the tile
+ const unsigned int tileHeight, // height of the tile
+ const unsigned int margin, // stability margin
+ const unsigned int nbTilesX,
+ const unsigned int nbTilesY,
+ const std::string temporaryFilesPrefix);
+
+ //Generate the .vrt corresponding to a list of input tiles
+ static void writeVRTRegions(TileListType const& inputRegions, GDALDataset *VRTOutput, std::string prefix, unsigned int band);
+
+ //Apply SLIC segmentation to the input image
+ void Run();
+
+ void SetInputImage(typename TInputImage::Pointer InputImagePointer);
+
+ protected:
+ SLICScheduler();
+
+ ~SLICScheduler() ITK_OVERRIDE;
+
+ private:
+ SLICScheduler(const Self &); //purposely not implemented
+ void operator =(const Self&); //purposely not implemented
+ unsigned int m_SpatialWidth;
+ RealType m_SpatialDistanceWeight;
+ unsigned int m_MaxIterationNumber;
+ RealType m_Threshold;
+ std::string m_Prefix;
+ std::string m_InputName;
+ ListSampleType m_SeedMeans;
+ unsigned int m_NbTilesX;
+ unsigned int m_NbTilesY;
+ RealType m_Margin;
+ typename TInputImage::Pointer m_InputImage;
+ };
+}
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSLICScheduler.txx"
+#endif
+
+#endif
diff --git a/include/otbSLICScheduler.h~ b/include/otbSLICScheduler.h~
new file mode 100644
index 0000000000000000000000000000000000000000..51b8eadbdf4b7abe33b92033a116679d4721ff22
--- /dev/null
+++ b/include/otbSLICScheduler.h~
@@ -0,0 +1,147 @@
+#ifndef otbSLICScheduler_h
+#define otbSLICScheduler_h
+
+#include "itkMacro.h"
+//OTB, ITK includes
+#include "otbVectorImage.h"
+#include "otbImage.h"
+#include "otbImageFileWriter.h"
+#include "otbImageFileReader.h"
+
+#include "itkImageRegionConstIteratorWithIndex.h"
+#include "itkImageRegionIteratorWithIndex.h"
+#include "itkImageRegionExclusionConstIteratorWithIndex.h"
+
+#include "itkImageToImageFilter.h"
+#include "itkNeighborhoodIterator.h"
+#include "itkExtractImageFilter.h"
+#include "otbMPIConfig.h"
+
+//Own includes
+#include "otbSLICScheduler.h"
+
+#include <itksys/SystemTools.hxx>
+
+#include <gdal.h>
+#include <gdal_priv.h>
+#if defined(__GNUC__) || defined(__clang__)
+# pragma GCC diagnostic push
+# pragma GCC diagnostic ignored "-Wshadow"
+#include <vrtdataset.h>
+# pragma GCC diagnostic pop
+#else
+#include <vrtdataset.h>
+#endif
+
+#include <ogr_spatialref.h>
+
+
+
+
+namespace otb{
+
+ template <class TInputImage, class TOutputLabelImage>
+ class ITK_EXPORT SLICScheduler : public itk::Object
+ {
+ public:
+ typedef SLICScheduler Self;
+ typedef itk::Object Superclass;
+ typedef itk::SmartPointer<Self> Pointer;
+ typedef itk::SmartPointer<const Self> ConstPointer;
+
+ /** Method for creation through the object factory. */
+ itkNewMacro(Self);
+
+ /** Runtime information support. */
+ itkTypeMacro(SLICScheduler, ImageToImageFilter);
+
+ /** ImageDimension constants */
+ itkStaticConstMacro(ImageDimension, unsigned int, TInputImage::ImageDimension);
+
+ typedef double RealType;
+ typedef typename TOutputLabelImage::PixelType LabelType;
+ typedef TOutputLabelImage OutputLabelImageType;
+ typedef TInputImage VectorImageType;
+ typedef typename TInputImage::RegionType RegionType;
+
+ typedef typename std::vector<RealType> SampleType;
+ typedef std::vector<SampleType> ListSampleType;
+
+ typedef std::vector<LabelType> NeighborhoodLabelsType;
+ typedef itk::ConstantBoundaryCondition<TOutputLabelImage> BoundaryConditionType;
+ typedef itk::NeighborhoodIterator<TOutputLabelImage,BoundaryConditionType> NeighborhoodLabelIteratorType;
+ typedef typename NeighborhoodLabelIteratorType::NeighborhoodType NeighborhoodType;
+ typedef typename NeighborhoodType::ConstIterator NeighborhoodContainerIteratorType;
+
+
+ typedef otb::ImageFileWriter<TOutputLabelImage> WriterType;
+
+ typedef std::vector<otb::ProcessingTile> TileListType;
+
+ // Sets the spatial width, i.e. average root of size of the output superpixels
+ itkSetMacro(SpatialWidth,unsigned int);
+
+ // Sets the spatial distance weight, i.e. the compacity/homogeneity ratio of SP
+ itkSetMacro(SpatialDistanceWeight,RealType);
+
+ //Sets the max number of iterations
+ itkSetMacro(MaxIterationNumber,unsigned int);
+
+ //Sets the residual threshold for stopping the algorithm
+ itkSetMacro(Threshold,RealType);
+
+ //Number of tiles in X direction
+ itkSetMacro(NbTilesX,unsigned int);
+
+ //Number of tiles in Y direction
+ itkSetMacro(NbTilesY,unsigned int);
+
+ //Margin (in SP widths)
+ itkSetMacro(Margin,RealType);
+
+ //Prefix for file writing
+ itkSetMacro(Prefix,std::string);
+
+ //Test
+ itkSetMacro(InputName,std::string);
+
+ static std::vector<ProcessingTile> SplitOTBImage(const TInputImage * imagePtr, // input image
+ const unsigned int tileWidth, // width of the tile
+ const unsigned int tileHeight, // height of the tile
+ const unsigned int margin, // stability margin
+ const unsigned int nbTilesX,
+ const unsigned int nbTilesY,
+ const std::string temporaryFilesPrefix);
+
+ static void writeVRTRegions(TileListType const& inputRegions, GDALDataset *VRTOutput, std::string prefix, unsigned int band);
+
+ void Run();
+
+ void SetInputImage(typename TInputImage::Pointer InputImagePointer);
+
+ protected:
+ SLICScheduler();
+
+ ~SLICScheduler() ITK_OVERRIDE;
+
+ private:
+ SLICScheduler(const Self &); //purposely not implemented
+ void operator =(const Self&); //purposely not implemented
+ unsigned int m_SpatialWidth;
+ RealType m_SpatialDistanceWeight;
+ unsigned int m_MaxIterationNumber;
+ RealType m_Threshold;
+ std::string m_Prefix;
+ std::string m_InputName;
+ ListSampleType m_SeedMeans;
+ unsigned int m_NbTilesX;
+ unsigned int m_NbTilesY;
+ RealType m_Margin;
+ typename TInputImage::Pointer m_InputImage;
+ };
+}
+#ifndef OTB_MANUAL_INSTANTIATION
+#include "otbSLICScheduler.txx"
+#endif
+
+#endif
diff --git a/include/otbSLICScheduler.txx b/include/otbSLICScheduler.txx
new file mode 100644
index 0000000000000000000000000000000000000000..54e395d62feebf5e2ad2df6758d29fea1183e67c
--- /dev/null
+++ b/include/otbSLICScheduler.txx
@@ -0,0 +1,553 @@
+#ifndef otbSLICScheduler_txx
+#define otbSLICScheduler_txx
+
+#include "otbSLICScheduler.h"
+
+namespace otb{
+
+ template <class TInputImage, class TOutputLabelImage>
+ SLICScheduler<TInputImage, TOutputLabelImage>
+ ::SLICScheduler() :
+ m_SpatialWidth(100),
+ m_SpatialDistanceWeight(1),
+ m_MaxIterationNumber(10),
+ m_Threshold(0.1),
+ m_NbTilesX(1),
+ m_NbTilesY(1),
+ m_Margin(0)
+ {}
+
+ template <class TInputImage, class TOutputLabelImage>
+ SLICScheduler<TInputImage, TOutputLabelImage>
+ ::~SLICScheduler()
+ {}
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICScheduler<TInputImage, TOutputLabelImage>
+ ::SetInputImage(typename TInputImage::Pointer InputImagePointer)
+ {
+ m_InputImage = InputImagePointer;
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICScheduler<TInputImage, TOutputLabelImage>
+ ::Run()
+ {
+ typename otb::MPIConfig::Pointer mpiConfig = otb::MPIConfig::Instance();
+ unsigned int myRank = mpiConfig->GetMyRank();
+ unsigned int nbProcs = mpiConfig->GetNbProcs();
+
+ const unsigned int imageSizeX = m_InputImage->GetLargestPossibleRegion().GetSize()[0];
+ const unsigned int imageSizeY = m_InputImage->GetLargestPossibleRegion().GetSize()[1];
+
+ const unsigned int tileWidth = imageSizeX/m_NbTilesX;
+ const unsigned int tileHeight = imageSizeY/m_NbTilesY;
+ const unsigned int pixelMargin = m_Margin*m_SpatialWidth;
+
+ typedef std::vector<otb::ProcessingTile> TileListType;
+
+ //Tiles for processing with no margins
+ TileListType tilesNoMargin = SplitOTBImage(m_InputImage, tileWidth, tileHeight, 0, m_NbTilesX, m_NbTilesY, m_Prefix);
+
+ //Tiles for processing with margins
+ TileListType tilesMargin = SplitOTBImage(m_InputImage, tileWidth, tileHeight, pixelMargin, m_NbTilesX, m_NbTilesY, m_Prefix);
+
+ TileListType whiteTiles;
+ TileListType blackTiles;
+ TileListType tilesToWrite;
+
+ //Generate list of black and white tiles to process
+ //Margin mode
+ if (m_Margin > 0) {
+ for(unsigned int y = 0 ; y < m_NbTilesY; ++y){
+ for(unsigned int x = 0 ; x < m_NbTilesX; ++x){
+ if(!((x+y)%2)){ //if x + y even
+ blackTiles.push_back(tilesNoMargin[x+m_NbTilesX*y]);
+ }
+ else{
+ whiteTiles.push_back(tilesMargin[x+m_NbTilesX*y]);
+ }
+ }
+ }
+ if(myRank < blackTiles.size())
+ {
+ //Dispatch according to rank of MPI process
+ unsigned int splitIdx = myRank;
+ while(splitIdx < blackTiles.size())
+ {
+ tilesToWrite.push_back(blackTiles[splitIdx]);
+ splitIdx+=nbProcs;
+ }
+ }
+ }
+ //Simple streaming mode
+ else{
+ if(myRank < tilesNoMargin.size())
+ {
+ unsigned int splitIdx = myRank;
+ while(splitIdx < tilesNoMargin.size())
+ {
+ tilesToWrite.push_back(tilesNoMargin[splitIdx]);
+ splitIdx+=nbProcs;
+ }
+ }
+ }
+ //First process the black tiles (no margins)
+ for (otb::ProcessingTile const& t : tilesToWrite)
+ {
+ typename ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(m_InputName);
+ reader->SetReleaseDataFlag(true);
+
+ typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
+ extractFilter->SetReleaseDataFlag(true);
+ extractFilter->SetInput(reader->GetOutput());
+ extractFilter->SetExtractionRegion(t.region);
+
+ typename SLICFilter::Pointer filter = SLICFilter::New();
+ filter->SetReleaseDataFlag(true);
+ filter->SetSpatialWidth(m_SpatialWidth);
+ filter->SetSpatialDistanceWeight(m_SpatialDistanceWeight);
+ filter->SetMaxIterationNumber(m_MaxIterationNumber);
+ filter->SetThreshold(m_Threshold);
+ filter->SetNbSPx(vcl_ceil((float) m_InputImage->GetLargestPossibleRegion().GetSize()[0]/m_SpatialWidth));
+ filter->SetTile(t);
+ filter->SetMargin(m_Margin);
+ filter->SetInput(extractFilter->GetOutput());
+
+ typename WriterType::Pointer writer = WriterType::New();
+ std::stringstream name;
+ name << m_Prefix << "_" << t.region.GetIndex()[0] << "_" << t.region.GetIndex()[1] << "_" << t.region.GetSize()[0] << "_" << t.region.GetSize()[1] << ".tif";
+ writer->SetInput(filter->GetOutput());
+ writer->SetFileName(name.str());
+ writer->SetNumberOfDivisionsStrippedStreaming(1);
+ writer->Update();
+ std::cout << myRank << " writing b region :" << name.str() << std::endl;
+ }
+ //Wait for all black regions to be processed
+ mpiConfig->barrier();
+
+ //Dispatch white tiles to different processes
+ if(m_Margin > 0){
+ tilesToWrite.clear();
+ for(unsigned int line = 0 ; line < m_NbTilesY ; line++){
+ unsigned int nbWhite = vcl_floor(m_NbTilesX/2);
+ if (m_NbTilesX%2){
+ nbWhite += line%2;
+ }
+ if(myRank < nbWhite)
+ {
+ unsigned int splitIdx = myRank;
+ while(splitIdx < nbWhite)
+ {
+ unsigned int ind;
+ if(m_NbTilesX%2)
+ ind = m_NbTilesX*vcl_floor(line/2) + splitIdx + vcl_floor(m_NbTilesX/2)*(line%2);
+ else
+ ind = splitIdx + nbWhite*line;
+ tilesToWrite.push_back(whiteTiles[ind]);
+ splitIdx+=nbProcs;
+ }
+ }
+
+ //Process tiles
+ for (otb::ProcessingTile t : tilesToWrite)
+ {
+ typename ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(m_InputName);
+ reader->SetReleaseDataFlag(true);
+ typename SLICFilter::Pointer filter = SLICFilter::New();
+ filter->SetSpatialWidth(m_SpatialWidth);
+ filter->SetSpatialDistanceWeight(m_SpatialDistanceWeight);
+ filter->SetMaxIterationNumber(m_MaxIterationNumber);
+ filter->SetThreshold(m_Threshold);
+ filter->SetNbSPx(vcl_ceil((float) m_InputImage->GetLargestPossibleRegion().GetSize()[0]/m_SpatialWidth));
+ filter->SetMargin(m_Margin);
+ typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
+ extractFilter->SetReleaseDataFlag(true);
+ extractFilter->SetInput(reader->GetOutput());
+ //Extract main region
+ extractFilter->SetExtractionRegion(t.region);
+ filter->SetTile(t);
+ filter->ClearMargins();
+ filter->SetInput(extractFilter->GetOutput());
+ filter->SetReleaseDataFlag(true);
+ //Read neighboring black tiles margin areas and load into filter
+ for (unsigned int i = 0; i < 8 ; i++){
+ if(pixelMargin > 0 && t.tileNeighbors[i]>=0){
+ typename TOutputLabelImage::RegionType marginRegion;
+ typename TOutputLabelImage::RegionType r;
+ bool read = false;
+ switch (i){
+ case otb::NBH_TOP:{
+ //Normal case
+ r = tilesNoMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ typename TOutputLabelImage::RegionType::SizeType size({r.GetSize()[0],pixelMargin});
+ index[1]+=r.GetSize()[1]-pixelMargin;
+ //Accomodate for corners
+ if(t.tileNeighbors[otb::NBH_LEFT]>=0){
+ //if left neighbor
+ index[0]+=pixelMargin;
+ size[0]-=pixelMargin;
+
+ }
+ if(t.tileNeighbors[otb::NBH_RIGHT]>=0){
+ //if right neighbor
+ size[0]-=pixelMargin;
+ }
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_BOTTOM:{
+ r = tilesNoMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ const typename TOutputLabelImage::RegionType::SizeType size({r.GetSize()[0],pixelMargin});
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_LEFT:{
+ r = tilesNoMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ index[0]+=r.GetSize()[0]-pixelMargin;
+ typename TOutputLabelImage::RegionType::SizeType size({pixelMargin,r.GetSize()[1]});
+ if(t.tileNeighbors[otb::NBH_TOP]>=0){
+ //if right neighbor
+ index[1]+=pixelMargin;
+ size[1]-=pixelMargin;
+ }
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_RIGHT:{
+ r = tilesNoMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ typename TOutputLabelImage::RegionType::SizeType size({pixelMargin,r.GetSize()[1]});
+ if(t.tileNeighbors[otb::NBH_TOP]>=0){
+ //if right neighbor
+ index[1]+=pixelMargin;
+ size[1]-=pixelMargin;
+ }
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_TOP_LEFT: {
+ r = tilesMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ index[0]+=r.GetSize()[0] - 2*pixelMargin;
+ index[1]+=r.GetSize()[1] - 2*pixelMargin;
+ const typename TOutputLabelImage::RegionType::SizeType size({2*pixelMargin,2*pixelMargin});
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_TOP_RIGHT: {
+ r = tilesMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ index[1]+=r.GetSize()[1] - 2*pixelMargin;
+ const typename TOutputLabelImage::RegionType::SizeType size({2*pixelMargin,2*pixelMargin});
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ }
+ if(read){
+ typename LabelExtractFilterType::Pointer marginExtract = LabelExtractFilterType::New();
+ std::stringstream readNeighborName;
+ readNeighborName << m_Prefix <<"_"<< r.GetIndex()[0] << "_" << r.GetIndex()[1] << "_"<< r.GetSize()[0] << "_"<< r.GetSize()[1] << ".tif";
+ typename LabelReaderType::Pointer labelReader = LabelReaderType::New();
+ labelReader->SetReleaseDataFlag(true);
+ labelReader->SetFileName(readNeighborName.str());
+ //Set the index of the output (reader gives 0 0)
+ marginExtract->SetInput(labelReader->GetOutput());
+ marginExtract->SetExtractionRegion(marginRegion);
+ marginExtract->Update();
+ marginRegion.SetIndex({marginRegion.GetIndex()[0]+r.GetIndex()[0],marginRegion.GetIndex()[1]+r.GetIndex()[1]});
+ marginExtract->GetOutput()->SetRegions(marginRegion);
+ filter->AddInputMargin(marginExtract->GetOutput());
+ }
+ }
+ }
+ filter->Update();
+ typename WriterType::Pointer writer = WriterType::New();
+ writer->SetInput(filter->GetOutput());
+ std::stringstream name;
+ name << m_Prefix << "_" << t.region.GetIndex()[0] << "_" << t.region.GetIndex()[1] << "_" << t.region.GetSize()[0] << "_" << t.region.GetSize()[1] << ".tif";
+ writer->SetFileName(name.str());
+ writer->SetNumberOfDivisionsStrippedStreaming(1);
+ writer->Update();
+ std::cout << myRank << " writing w region :" << name.str() << std::endl;
+ }
+ tilesToWrite.clear();
+ mpiConfig->barrier();
+ }
+ }
+ try
+ {
+ if (myRank == 0)
+ {
+ // VRT Filename
+ std::stringstream vrtfOut;
+ vrtfOut<< m_Prefix << ".vrt";
+
+ // Data type
+ GDALDataType dataType;
+ dataType = GDALGetDataTypeByName("UInt32");
+
+ // Get VRT driver
+ GDALAllRegister();
+ GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("VRT");
+ if (driver == NULL) {
+ mpiConfig->logError("Error opening VRT driver.");
+ mpiConfig->abort(EXIT_FAILURE);
+ }
+
+ // Create output raster
+ GDALDataset *VRTOutput = driver->Create(vrtfOut.str().c_str(),
+ imageSizeX,
+ imageSizeY,
+ 0,
+ dataType,
+ NULL); // No options
+ if (VRTOutput == NULL) {
+ mpiConfig->logError("driver->Create call failed.\n");
+ mpiConfig->abort(EXIT_FAILURE);
+ }
+
+ // Set GeoTransform
+ double gt[6];
+ gt[0] = m_InputImage->GetOrigin()[0] - 0.5*m_InputImage->GetSignedSpacing()[0];
+ gt[1] = m_InputImage->GetSignedSpacing()[0];
+ gt[2] = 0.0;
+ gt[3] = m_InputImage->GetOrigin()[1] - 0.5*m_InputImage->GetSignedSpacing()[1];
+ gt[4] = 0.0;
+ gt[5] = m_InputImage->GetSignedSpacing()[1];
+ VRTOutput->SetGeoTransform(gt);
+
+ // Set projection
+ OGRSpatialReference out_sr;
+ char *wkt = NULL;
+ out_sr.SetFromUserInput(m_InputImage->GetProjectionRef().c_str());
+ out_sr.exportToWkt(&wkt);
+ VRTOutput->SetProjection(wkt);
+ VRTOutput->AddBand(dataType, NULL);
+ if(m_Margin == 0){
+ writeVRTRegions(tilesNoMargin, VRTOutput, m_Prefix, 1);
+ }
+ else{
+ writeVRTRegions(blackTiles, VRTOutput, m_Prefix, 1);
+ writeVRTRegions(whiteTiles, VRTOutput, m_Prefix, 1);
+ }
+ // Close
+ OGRFree(wkt);
+ GDALClose(VRTOutput);
+ }
+ }
+ catch (itk::ExceptionObject& err)
+ {
+ std::stringstream message;
+ message << "ExceptionObject caught: " << err << std::endl;
+ mpiConfig->logError(message.str());
+ mpiConfig->abort(EXIT_FAILURE);
+ }
+ }
+ template<class TInputImage, class TOutputLabelImage>
+ void
+ SLICScheduler<TInputImage,TOutputLabelImage>
+ ::writeVRTRegions(TileListType const& inputRegions, GDALDataset *VRTOutput, std::string prefix, unsigned int band){
+ for(otb::ProcessingTile t : inputRegions)
+ {
+ int tileSizeX = t.region.GetSize()[0];
+ int tileSizeY = t.region.GetSize()[1];
+ int tileIndexX = t.region.GetIndex()[0];
+ int tileIndexY = t.region.GetIndex()[1];
+ std::stringstream tileFileName;
+ tileFileName <<prefix<<"_"<<tileIndexX<<"_"<<tileIndexY<<"_"<<tileSizeX<<"_"<<tileSizeY<<".tif";
+
+ GDALDataset *dataset = (GDALDataset*) GDALOpen(tileFileName.str().c_str(), GA_ReadOnly);
+
+ VRTSourcedRasterBand *VRTBand = dynamic_cast<VRTSourcedRasterBand*> (VRTOutput->GetRasterBand(band));
+ VRTBand->AddComplexSource(dataset->GetRasterBand(band),
+ 0,//Xoffset
+ 0,//Yoffset
+ tileSizeX,//Dimensions
+ tileSizeY,//Dimensions
+ tileIndexX,//Index
+ tileIndexY,//Index
+ tileSizeX,//Dimensions
+ tileSizeY,//Dimensions
+ 0.0,//Scale offset
+ 1,//Scale Ratio
+ 0);//NodataValue
+ }
+ }
+ template<class TInputImage, class TOutputLabelImage>
+ std::vector<ProcessingTile>
+ SLICScheduler<TInputImage,TOutputLabelImage>
+ ::SplitOTBImage(const TInputImage * imagePtr, // input image
+ const unsigned int tileWidth, // width of the tile
+ const unsigned int tileHeight, // height of the tile
+ const unsigned int margin, // stability margin
+ const unsigned int nbTilesX,
+ const unsigned int nbTilesY,
+ const std::string temporaryFilesPrefix)
+ {
+ std::vector<ProcessingTile> tiles;
+ // Image size
+ const unsigned int imageWidth = imagePtr->GetLargestPossibleRegion().GetSize()[0];
+ const unsigned int imageHeight = imagePtr->GetLargestPossibleRegion().GetSize()[1];
+
+ /* Loop over the tiles*/
+ tiles.assign(nbTilesX * nbTilesY, ProcessingTile());
+ unsigned int i = 0;
+ for(unsigned int row = 0; row < nbTilesY; ++row)
+ {
+ for(unsigned int col = 0; col < nbTilesX; ++col)
+ {
+ /* Local variables for the next loop */
+ // Compute current tile start and size
+ const unsigned int startX = col * tileWidth;
+ const unsigned int startY = row * tileHeight; // Upper left coordinates of the tile.
+ unsigned int sizeX = tileWidth;
+ unsigned int sizeY = tileHeight; // Size of the tiles.
+
+ // Current tile size might be different for right and bottom borders
+ if (col == nbTilesX -1)
+ {
+ sizeX += imageWidth % tileWidth;
+ }
+ if (row == nbTilesY -1)
+ {
+ sizeY += imageHeight % tileHeight;
+ }
+
+ /* Margin at the top ? */
+ if( row > 0 )
+ {
+ tiles[i].margin[POS_TOP] = margin;
+ tiles[i].rows[0] = row * tileHeight;
+ tiles[i].tileNeighbors[NBH_TOP] = i - nbTilesX;
+ }
+ else
+ {
+ // Tile is on the top row --> no top margin
+ tiles[i].margin[POS_TOP] = 0;
+ tiles[i].rows[0] = 0;
+ tiles[i].tileNeighbors[NBH_TOP] = -1;
+ }
+
+ /* Margin at the right */
+ if( col < nbTilesX - 1 )
+ {
+ tiles[i].margin[POS_RIGHT] = margin;
+ tiles[i].columns[1] = col * tileWidth + sizeX - 1; //sizeX
+ tiles[i].tileNeighbors[NBH_RIGHT] = i+1;
+ }
+ else
+ {
+ // Tile is on the right column --> no right margin
+ tiles[i].margin[POS_RIGHT] = 0;
+ tiles[i].columns[1] = imageWidth - 1;
+ tiles[i].tileNeighbors[NBH_RIGHT] = -1;
+ }
+
+ /* Margin at the bottom */
+ if( row < nbTilesY - 1)
+ {
+ tiles[i].margin[POS_BOTTOM] = margin;
+ tiles[i].rows[1] = row * tileHeight + sizeY - 1; // sizeY
+ tiles[i].tileNeighbors[NBH_BOTTOM] = i + nbTilesX;
+ }
+ else
+ {
+ // Tile is on the bottom --> no bottom margin
+ tiles[i].margin[POS_BOTTOM] = 0;
+ tiles[i].rows[1] = imageHeight - 1;
+ tiles[i].tileNeighbors[NBH_BOTTOM] = -1;
+ }
+
+ /* Margin at the left */
+ if( col > 0 )
+ {
+ tiles[i].margin[POS_LEFT] = margin;
+ tiles[i].columns[0] = col * tileWidth;
+ tiles[i].tileNeighbors[NBH_LEFT] = i-1;
+ }
+ else
+ {
+ // Tile is on the left --> no left margin
+ tiles[i].margin[POS_LEFT] = 0;
+ tiles[i].columns[0] = 0;
+ tiles[i].tileNeighbors[NBH_LEFT] = -1;
+ }
+
+ /* Store the tile region */
+ typename TInputImage::IndexType index;
+ index[0] = startX - tiles[i].margin[POS_LEFT];
+ index[1] = startY - tiles[i].margin[POS_TOP];
+ typename TInputImage::SizeType size;
+ size[0] = sizeX + tiles[i].margin[POS_LEFT] + tiles[i].margin[POS_RIGHT];
+ size[1] = sizeY + tiles[i].margin[POS_TOP] + tiles[i].margin[POS_BOTTOM];
+ typename TInputImage::RegionType region(index, size);
+ tiles[i].region = region;
+
+ // std::cout << "Tile " << i << ": start at " << index << " with size " << size << std::endl;
+
+ /* Is there a neighbor at the rop right */
+ if(row > 0 && col < nbTilesX - 1)
+ tiles[i].tileNeighbors[NBH_TOP_RIGHT] = i - nbTilesX + 1;
+ else
+ tiles[i].tileNeighbors[NBH_TOP_RIGHT] = -1;
+
+ /* Is there a neighbor at the bottom right */
+ if(col < nbTilesX - 1 && row < nbTilesY - 1)
+ tiles[i].tileNeighbors[NBH_BOTTOM_RIGHT] = i + nbTilesX + 1;
+ else
+ tiles[i].tileNeighbors[NBH_BOTTOM_RIGHT] = -1;
+
+ /* Is there a neighbor at the bottom left */
+ if(row < nbTilesY - 1 && col > 0)
+ tiles[i].tileNeighbors[NBH_BOTTOM_LEFT] = i + nbTilesX - 1;
+ else
+ tiles[i].tileNeighbors[NBH_BOTTOM_LEFT] = -1;
+
+ /* Is there a neighbor at the top left */
+ if(col > 0 && row > 0)
+ tiles[i].tileNeighbors[NBH_TOP_LEFT] = i - nbTilesX - 1;
+ else
+ tiles[i].tileNeighbors[NBH_TOP_LEFT] = -1;
+
+ const std::string suffix = std::to_string(row) + "_" + std::to_string(col) + ".tif";
+ tiles[i].nodeFileName = temporaryFilesPrefix + "_node_" + suffix;
+ tiles[i].edgeFileName = temporaryFilesPrefix + "_edge_" + suffix;
+ tiles[i].nodeMarginFileName = temporaryFilesPrefix + "_nodeMargin_" + suffix;
+ tiles[i].edgeMarginFileName = temporaryFilesPrefix + "_edgeMargin_" + suffix;
+
+ i++;
+ } // end for(unsigned int col = 0; col < nbTilesX; ++col)
+
+ } // for(unsigned int row = 0; row < nbTilesY; ++row)
+
+ return tiles;
+ }
+} //end namespace otb
+#endif
+
diff --git a/include/otbSLICScheduler.txx~ b/include/otbSLICScheduler.txx~
new file mode 100644
index 0000000000000000000000000000000000000000..c01946d8d79ba51b264f5adc259a320f85ce048e
--- /dev/null
+++ b/include/otbSLICScheduler.txx~
@@ -0,0 +1,552 @@
+#ifndef otbSLICScheduler_txx
+#define otbSLICScheduler_txx
+
+#include "otbSLICScheduler.h"
+
+namespace otb{
+
+ template <class TInputImage, class TOutputLabelImage>
+ SLICScheduler<TInputImage, TOutputLabelImage>
+ ::SLICScheduler() :
+ m_SpatialWidth(100),
+ m_SpatialDistanceWeight(1),
+ m_MaxIterationNumber(10),
+ m_Threshold(0.1),
+ m_NbTilesX(1),
+ m_NbTilesY(1),
+ m_Margin(0)
+ {}
+
+ template <class TInputImage, class TOutputLabelImage>
+ SLICScheduler<TInputImage, TOutputLabelImage>
+ ::~SLICScheduler()
+ {}
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICScheduler<TInputImage, TOutputLabelImage>
+ ::SetInputImage(typename TInputImage::Pointer InputImagePointer)
+ {
+ m_InputImage = InputImagePointer;
+ }
+
+ template <class TInputImage, class TOutputLabelImage>
+ void
+ SLICScheduler<TInputImage, TOutputLabelImage>
+ ::Run()
+ {
+ typename otb::MPIConfig::Pointer mpiConfig = otb::MPIConfig::Instance();
+ unsigned int myRank = mpiConfig->GetMyRank();
+ unsigned int nbProcs = mpiConfig->GetNbProcs();
+
+ const unsigned int imageSizeX = m_InputImage->GetLargestPossibleRegion().GetSize()[0];
+ const unsigned int imageSizeY = m_InputImage->GetLargestPossibleRegion().GetSize()[1];
+
+ const unsigned int tileWidth = imageSizeX/m_NbTilesX;
+ const unsigned int tileHeight = imageSizeY/m_NbTilesY;
+ const unsigned int pixelMargin = m_Margin*m_SpatialWidth;
+
+ typedef std::vector<otb::ProcessingTile> TileListType;
+
+ TileListType tilesNoMargin = SplitOTBImage(m_InputImage, tileWidth, tileHeight, 0, m_NbTilesX, m_NbTilesY, m_Prefix);
+
+ TileListType tilesMargin = SplitOTBImage(m_InputImage, tileWidth, tileHeight, pixelMargin, m_NbTilesX, m_NbTilesY, m_Prefix);
+
+ TileListType whiteTiles;
+ TileListType blackTiles;
+ TileListType tilesToWrite;
+ //Margin mode
+ if (m_Margin > 0) {
+ for(unsigned int y = 0 ; y < m_NbTilesY; ++y){
+ for(unsigned int x = 0 ; x < m_NbTilesX; ++x){
+ if(!((x+y)%2)){ //if x + y even
+ blackTiles.push_back(tilesNoMargin[x+m_NbTilesX*y]);
+ }
+ else{
+ whiteTiles.push_back(tilesMargin[x+m_NbTilesX*y]);
+ }
+ }
+ }
+ if(myRank < blackTiles.size())
+ {
+ unsigned int splitIdx = myRank;
+ while(splitIdx < blackTiles.size())
+ {
+ tilesToWrite.push_back(blackTiles[splitIdx]);
+ splitIdx+=nbProcs;
+ }
+ }
+ }
+ //Simple streaming mode
+ else{
+ if(myRank < tilesNoMargin.size())
+ {
+ unsigned int splitIdx = myRank;
+ while(splitIdx < tilesNoMargin.size())
+ {
+ tilesToWrite.push_back(tilesNoMargin[splitIdx]);
+ splitIdx+=nbProcs;
+ }
+ }
+ }
+ typedef otb::SLICFilter<TInputImage, TOutputLabelImage> SLICFilter;
+ typedef itk::ExtractImageFilter<TInputImage,TInputImage> ExtractFilterType;
+ typedef otb::ImageFileReader<TInputImage> ReaderType;
+
+ //Black squares
+ for (otb::ProcessingTile const& t : tilesToWrite)
+ {
+ typename ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(m_InputName);
+ reader->SetReleaseDataFlag(true);
+
+ typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
+ extractFilter->SetReleaseDataFlag(true);
+ extractFilter->SetInput(reader->GetOutput());
+ extractFilter->SetExtractionRegion(t.region);
+
+ typename SLICFilter::Pointer filter = SLICFilter::New();
+ filter->SetReleaseDataFlag(true);
+ filter->SetSpatialWidth(m_SpatialWidth);
+ filter->SetSpatialDistanceWeight(m_SpatialDistanceWeight);
+ filter->SetMaxIterationNumber(m_MaxIterationNumber);
+ filter->SetThreshold(m_Threshold);
+ filter->SetNbSPx(vcl_ceil((float) m_InputImage->GetLargestPossibleRegion().GetSize()[0]/m_SpatialWidth));
+ filter->SetTile(t);
+ filter->SetMargin(m_Margin);
+ filter->SetInput(extractFilter->GetOutput());
+
+ typename WriterType::Pointer writer = WriterType::New();
+ std::stringstream name;
+ name << m_Prefix << "_" << t.region.GetIndex()[0] << "_" << t.region.GetIndex()[1] << "_" << t.region.GetSize()[0] << "_" << t.region.GetSize()[1] << ".tif";
+ writer->SetInput(scheduler->GetOutput());
+ writer->SetFileName(name.str());
+ writer->SetNumberOfDivisionsStrippedStreaming(1);
+ writer->Update();
+ std::cout << myRank << " writing b region :" << name.str() << std::endl;
+ }
+
+ mpiConfig->barrier();
+
+ if(m_Margin > 0){
+ tilesToWrite.clear();
+ for(unsigned int line = 0 ; line < m_NbTilesY ; line++){
+ unsigned int nbWhite = vcl_floor(m_NbTilesX/2);
+ if (m_NbTilesX%2){
+ nbWhite += line%2;
+ }
+ if(myRank < nbWhite)
+ {
+ unsigned int splitIdx = myRank;
+ while(splitIdx < nbWhite)
+ {
+ unsigned int ind;
+ if(m_NbTilesX%2)
+ ind = m_NbTilesX*vcl_floor(line/2) + splitIdx + vcl_floor(m_NbTilesX/2)*(line%2);
+ else
+ ind = splitIdx + nbWhite*line;
+ tilesToWrite.push_back(whiteTiles[ind]);
+ splitIdx+=nbProcs;
+ }
+ }
+ typedef otb::ImageFileReader<TOutputLabelImage> LabelReaderType;
+ typedef itk::ExtractImageFilter<TOutputLabelImage, TOutputLabelImage> LabelExtractFilterType;
+ for (otb::ProcessingTile t : tilesToWrite)
+ {
+ typename ReaderType::Pointer reader = ReaderType::New();
+ reader->SetFileName(m_InputName);
+ reader->SetReleaseDataFlag(true);
+ typename SLICScheduler::Pointer scheduler = SLICScheduler::New();
+ scheduler->SetSpatialWidth(m_SpatialWidth);
+ scheduler->SetSpatialDistanceWeight(m_SpatialDistanceWeight);
+ scheduler->SetMaxIterationNumber(m_MaxIterationNumber);
+ scheduler->SetThreshold(m_Threshold);
+ scheduler->SetNbSPx(vcl_ceil((float) m_InputImage->GetLargestPossibleRegion().GetSize()[0]/m_SpatialWidth));
+ scheduler->SetMargin(m_Margin);
+ typename ExtractFilterType::Pointer extractFilter = ExtractFilterType::New();
+ extractFilter->SetReleaseDataFlag(true);
+ extractFilter->SetInput(reader->GetOutput());
+ //Extract main region
+ extractFilter->SetExtractionRegion(t.region);
+ scheduler->SetTile(t);
+ scheduler->ClearMargins();
+ scheduler->SetInput(extractFilter->GetOutput());
+ scheduler->SetReleaseDataFlag(true);
+ //Read neighboring black tiles margin areas and load into filter
+ for (unsigned int i = 0; i < 8 ; i++){
+ if(pixelMargin > 0 && t.tileNeighbors[i]>=0){
+ typename TOutputLabelImage::RegionType marginRegion;
+ typename TOutputLabelImage::RegionType r;
+ bool read = false;
+ switch (i){
+ case otb::NBH_TOP:{
+ //Normal case
+ r = tilesNoMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ typename TOutputLabelImage::RegionType::SizeType size({r.GetSize()[0],pixelMargin});
+ index[1]+=r.GetSize()[1]-pixelMargin;
+ //Accomodate for corners
+ if(t.tileNeighbors[otb::NBH_LEFT]>=0){
+ //if left neighbor
+ index[0]+=pixelMargin;
+ size[0]-=pixelMargin;
+
+ }
+ if(t.tileNeighbors[otb::NBH_RIGHT]>=0){
+ //if right neighbor
+ size[0]-=pixelMargin;
+ }
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_BOTTOM:{
+ r = tilesNoMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ const typename TOutputLabelImage::RegionType::SizeType size({r.GetSize()[0],pixelMargin});
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_LEFT:{
+ r = tilesNoMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ index[0]+=r.GetSize()[0]-pixelMargin;
+ typename TOutputLabelImage::RegionType::SizeType size({pixelMargin,r.GetSize()[1]});
+ if(t.tileNeighbors[otb::NBH_TOP]>=0){
+ //if right neighbor
+ index[1]+=pixelMargin;
+ size[1]-=pixelMargin;
+ }
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_RIGHT:{
+ r = tilesNoMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ typename TOutputLabelImage::RegionType::SizeType size({pixelMargin,r.GetSize()[1]});
+ if(t.tileNeighbors[otb::NBH_TOP]>=0){
+ //if right neighbor
+ index[1]+=pixelMargin;
+ size[1]-=pixelMargin;
+ }
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_TOP_LEFT: {
+ r = tilesMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ index[0]+=r.GetSize()[0] - 2*pixelMargin;
+ index[1]+=r.GetSize()[1] - 2*pixelMargin;
+ const typename TOutputLabelImage::RegionType::SizeType size({2*pixelMargin,2*pixelMargin});
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ case otb::NBH_TOP_RIGHT: {
+ r = tilesMargin[t.tileNeighbors[i]].region;
+ // typename TOutputLabelImage::RegionType::IndexType index = r.GetIndex();
+ typename TOutputLabelImage::RegionType::IndexType index({0,0});
+ index[1]+=r.GetSize()[1] - 2*pixelMargin;
+ const typename TOutputLabelImage::RegionType::SizeType size({2*pixelMargin,2*pixelMargin});
+ marginRegion.SetIndex(index);
+ marginRegion.SetSize(size);
+ read = true;
+ break;
+ }
+ }
+ if(read){
+ typename LabelExtractFilterType::Pointer marginExtract = LabelExtractFilterType::New();
+ std::stringstream readNeighborName;
+ readNeighborName << m_Prefix <<"_"<< r.GetIndex()[0] << "_" << r.GetIndex()[1] << "_"<< r.GetSize()[0] << "_"<< r.GetSize()[1] << ".tif";
+ typename LabelReaderType::Pointer labelReader = LabelReaderType::New();
+ labelReader->SetReleaseDataFlag(true);
+ labelReader->SetFileName(readNeighborName.str());
+ //Set the index of the output (reader gives 0 0)
+ marginExtract->SetInput(labelReader->GetOutput());
+ marginExtract->SetExtractionRegion(marginRegion);
+ marginExtract->Update();
+ marginRegion.SetIndex({marginRegion.GetIndex()[0]+r.GetIndex()[0],marginRegion.GetIndex()[1]+r.GetIndex()[1]});
+ marginExtract->GetOutput()->SetRegions(marginRegion);
+ scheduler->AddInputMargin(marginExtract->GetOutput());
+ }
+ }
+ }
+ scheduler->Update();
+ typename WriterType::Pointer writer = WriterType::New();
+ writer->SetInput(scheduler->GetOutput());
+ std::stringstream name;
+ name << m_Prefix << "_" << t.region.GetIndex()[0] << "_" << t.region.GetIndex()[1] << "_" << t.region.GetSize()[0] << "_" << t.region.GetSize()[1] << ".tif";
+ writer->SetFileName(name.str());
+ writer->SetNumberOfDivisionsStrippedStreaming(1);
+ writer->Update();
+ std::cout << myRank << " writing w region :" << name.str() << std::endl;
+ }
+ tilesToWrite.clear();
+ mpiConfig->barrier();
+ }
+ }
+ try
+ {
+ if (myRank == 0)
+ {
+ // VRT Filename
+ std::stringstream vrtfOut;
+ vrtfOut<< m_Prefix << ".vrt";
+
+ // Data type
+ GDALDataType dataType;
+ dataType = GDALGetDataTypeByName("UInt32");
+
+ // Get VRT driver
+ GDALAllRegister();
+ GDALDriver *driver = GetGDALDriverManager()->GetDriverByName("VRT");
+ if (driver == NULL) {
+ mpiConfig->logError("Error opening VRT driver.");
+ mpiConfig->abort(EXIT_FAILURE);
+ }
+
+ // Create output raster
+ GDALDataset *VRTOutput = driver->Create(vrtfOut.str().c_str(),
+ imageSizeX,
+ imageSizeY,
+ 0,
+ dataType,
+ NULL); // No options
+ if (VRTOutput == NULL) {
+ mpiConfig->logError("driver->Create call failed.\n");
+ mpiConfig->abort(EXIT_FAILURE);
+ }
+
+ // Set GeoTransform
+ double gt[6];
+ gt[0] = m_InputImage->GetOrigin()[0] - 0.5*m_InputImage->GetSignedSpacing()[0];
+ gt[1] = m_InputImage->GetSignedSpacing()[0];
+ gt[2] = 0.0;
+ gt[3] = m_InputImage->GetOrigin()[1] - 0.5*m_InputImage->GetSignedSpacing()[1];
+ gt[4] = 0.0;
+ gt[5] = m_InputImage->GetSignedSpacing()[1];
+ VRTOutput->SetGeoTransform(gt);
+
+ // Set projection
+ OGRSpatialReference out_sr;
+ char *wkt = NULL;
+ out_sr.SetFromUserInput(m_InputImage->GetProjectionRef().c_str());
+ out_sr.exportToWkt(&wkt);
+ VRTOutput->SetProjection(wkt);
+ VRTOutput->AddBand(dataType, NULL);
+ if(m_Margin == 0){
+ writeVRTRegions(tilesNoMargin, VRTOutput, m_Prefix, 1);
+ }
+ else{
+ writeVRTRegions(blackTiles, VRTOutput, m_Prefix, 1);
+ writeVRTRegions(whiteTiles, VRTOutput, m_Prefix, 1);
+ }
+ // Close
+ OGRFree(wkt);
+ GDALClose(VRTOutput);
+ }
+ }
+ catch (itk::ExceptionObject& err)
+ {
+ std::stringstream message;
+ message << "ExceptionObject caught: " << err << std::endl;
+ mpiConfig->logError(message.str());
+ mpiConfig->abort(EXIT_FAILURE);
+ }
+ }
+ template<class TInputImage, class TOutputLabelImage>
+ void
+ SLICScheduler<TInputImage,TOutputLabelImage>
+ ::writeVRTRegions(TileListType const& inputRegions, GDALDataset *VRTOutput, std::string prefix, unsigned int band){
+ for(otb::ProcessingTile t : inputRegions)
+ {
+ int tileSizeX = t.region.GetSize()[0];
+ int tileSizeY = t.region.GetSize()[1];
+ int tileIndexX = t.region.GetIndex()[0];
+ int tileIndexY = t.region.GetIndex()[1];
+ std::stringstream tileFileName;
+ tileFileName <<prefix<<"_"<<tileIndexX<<"_"<<tileIndexY<<"_"<<tileSizeX<<"_"<<tileSizeY<<".tif";
+ // std::cout<<tileFileName.str()<<std::endl;
+
+ GDALDataset *dataset = (GDALDataset*) GDALOpen(tileFileName.str().c_str(), GA_ReadOnly);
+
+ VRTSourcedRasterBand *VRTBand = dynamic_cast<VRTSourcedRasterBand*> (VRTOutput->GetRasterBand(band));
+ VRTBand->AddComplexSource(dataset->GetRasterBand(band),
+ 0,//Xoffset
+ 0,//Yoffset
+ tileSizeX,//Dimensions
+ tileSizeY,//Dimensions
+ tileIndexX,//Index
+ tileIndexY,//Index
+ tileSizeX,//Dimensions
+ tileSizeY,//Dimensions
+ 0.0,//Scale offset
+ 1,//Scale Ratio
+ 0);//NodataValue
+ }
+ }
+ template<class TInputImage, class TOutputLabelImage>
+ std::vector<ProcessingTile>
+ SLICScheduler<TInputImage,TOutputLabelImage>
+ ::SplitOTBImage(const TInputImage * imagePtr, // input image
+ const unsigned int tileWidth, // width of the tile
+ const unsigned int tileHeight, // height of the tile
+ const unsigned int margin, // stability margin
+ const unsigned int nbTilesX,
+ const unsigned int nbTilesY,
+ const std::string temporaryFilesPrefix)
+ {
+ std::vector<ProcessingTile> tiles;
+ // Image size
+ const unsigned int imageWidth = imagePtr->GetLargestPossibleRegion().GetSize()[0];
+ const unsigned int imageHeight = imagePtr->GetLargestPossibleRegion().GetSize()[1];
+
+ /* Loop over the tiles*/
+ tiles.assign(nbTilesX * nbTilesY, ProcessingTile());
+ unsigned int i = 0;
+ for(unsigned int row = 0; row < nbTilesY; ++row)
+ {
+ for(unsigned int col = 0; col < nbTilesX; ++col)
+ {
+ /* Local variables for the next loop */
+ // Compute current tile start and size
+ const unsigned int startX = col * tileWidth;
+ const unsigned int startY = row * tileHeight; // Upper left coordinates of the tile.
+ unsigned int sizeX = tileWidth;
+ unsigned int sizeY = tileHeight; // Size of the tiles.
+
+ // Current tile size might be different for right and bottom borders
+ if (col == nbTilesX -1)
+ {
+ sizeX += imageWidth % tileWidth;
+ }
+ if (row == nbTilesY -1)
+ {
+ sizeY += imageHeight % tileHeight;
+ }
+
+ /* Margin at the top ? */
+ if( row > 0 )
+ {
+ tiles[i].margin[POS_TOP] = margin;
+ tiles[i].rows[0] = row * tileHeight;
+ tiles[i].tileNeighbors[NBH_TOP] = i - nbTilesX;
+ }
+ else
+ {
+ // Tile is on the top row --> no top margin
+ tiles[i].margin[POS_TOP] = 0;
+ tiles[i].rows[0] = 0;
+ tiles[i].tileNeighbors[NBH_TOP] = -1;
+ }
+
+ /* Margin at the right */
+ if( col < nbTilesX - 1 )
+ {
+ tiles[i].margin[POS_RIGHT] = margin;
+ tiles[i].columns[1] = col * tileWidth + sizeX - 1; //sizeX
+ tiles[i].tileNeighbors[NBH_RIGHT] = i+1;
+ }
+ else
+ {
+ // Tile is on the right column --> no right margin
+ tiles[i].margin[POS_RIGHT] = 0;
+ tiles[i].columns[1] = imageWidth - 1;
+ tiles[i].tileNeighbors[NBH_RIGHT] = -1;
+ }
+
+ /* Margin at the bottom */
+ if( row < nbTilesY - 1)
+ {
+ tiles[i].margin[POS_BOTTOM] = margin;
+ tiles[i].rows[1] = row * tileHeight + sizeY - 1; // sizeY
+ tiles[i].tileNeighbors[NBH_BOTTOM] = i + nbTilesX;
+ }
+ else
+ {
+ // Tile is on the bottom --> no bottom margin
+ tiles[i].margin[POS_BOTTOM] = 0;
+ tiles[i].rows[1] = imageHeight - 1;
+ tiles[i].tileNeighbors[NBH_BOTTOM] = -1;
+ }
+
+ /* Margin at the left */
+ if( col > 0 )
+ {
+ tiles[i].margin[POS_LEFT] = margin;
+ tiles[i].columns[0] = col * tileWidth;
+ tiles[i].tileNeighbors[NBH_LEFT] = i-1;
+ }
+ else
+ {
+ // Tile is on the left --> no left margin
+ tiles[i].margin[POS_LEFT] = 0;
+ tiles[i].columns[0] = 0;
+ tiles[i].tileNeighbors[NBH_LEFT] = -1;
+ }
+
+ /* Store the tile region */
+ typename TInputImage::IndexType index;
+ index[0] = startX - tiles[i].margin[POS_LEFT];
+ index[1] = startY - tiles[i].margin[POS_TOP];
+ typename TInputImage::SizeType size;
+ size[0] = sizeX + tiles[i].margin[POS_LEFT] + tiles[i].margin[POS_RIGHT];
+ size[1] = sizeY + tiles[i].margin[POS_TOP] + tiles[i].margin[POS_BOTTOM];
+ typename TInputImage::RegionType region(index, size);
+ tiles[i].region = region;
+
+ // std::cout << "Tile " << i << ": start at " << index << " with size " << size << std::endl;
+
+ /* Is there a neighbor at the rop right */
+ if(row > 0 && col < nbTilesX - 1)
+ tiles[i].tileNeighbors[NBH_TOP_RIGHT] = i - nbTilesX + 1;
+ else
+ tiles[i].tileNeighbors[NBH_TOP_RIGHT] = -1;
+
+ /* Is there a neighbor at the bottom right */
+ if(col < nbTilesX - 1 && row < nbTilesY - 1)
+ tiles[i].tileNeighbors[NBH_BOTTOM_RIGHT] = i + nbTilesX + 1;
+ else
+ tiles[i].tileNeighbors[NBH_BOTTOM_RIGHT] = -1;
+
+ /* Is there a neighbor at the bottom left */
+ if(row < nbTilesY - 1 && col > 0)
+ tiles[i].tileNeighbors[NBH_BOTTOM_LEFT] = i + nbTilesX - 1;
+ else
+ tiles[i].tileNeighbors[NBH_BOTTOM_LEFT] = -1;
+
+ /* Is there a neighbor at the top left */
+ if(col > 0 && row > 0)
+ tiles[i].tileNeighbors[NBH_TOP_LEFT] = i - nbTilesX - 1;
+ else
+ tiles[i].tileNeighbors[NBH_TOP_LEFT] = -1;
+
+ const std::string suffix = std::to_string(row) + "_" + std::to_string(col) + ".tif";
+ tiles[i].nodeFileName = temporaryFilesPrefix + "_node_" + suffix;
+ tiles[i].edgeFileName = temporaryFilesPrefix + "_edge_" + suffix;
+ tiles[i].nodeMarginFileName = temporaryFilesPrefix + "_nodeMargin_" + suffix;
+ tiles[i].edgeMarginFileName = temporaryFilesPrefix + "_edgeMargin_" + suffix;
+
+ i++;
+ } // end for(unsigned int col = 0; col < nbTilesX; ++col)
+
+ } // for(unsigned int row = 0; row < nbTilesY; ++row)
+
+ return tiles;
+ }
+} //end namespace otb
+#endif
+
diff --git a/include/simplePointCalculator.h b/include/simplePointCalculator.h
new file mode 100644
index 0000000000000000000000000000000000000000..0e7b278a0860b4500693d1ea44356cb938424ec9
--- /dev/null
+++ b/include/simplePointCalculator.h
@@ -0,0 +1,18 @@
+#ifndef simplePointCalculator_h
+#define simplePointCalculator_h
+
+#include <array>
+#include <vector>
+#include <iostream>
+#include <bitset>
+template <class IteratorType, class LabelType>
+class simplePointCalculator{
+public:
+ int isSimpleLabel (IteratorType const& neighborhoodLabelsBegin, IteratorType const& neighborhoodLabelsEnd, LabelType const& evalLabel) const;
+ void generateInitializerList();
+private:
+ constexpr static std::array<bool,256> T4Table{0,0,1,1,0,0,1,1,1,1,0,1,0,0,0,1,1,0,0,0,1,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,1,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,1,0,1,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,1,0,1,1,0,0,0,0,0,0,0,1,1,0,0,0,0,0,0,0,1,1,0,1,0,0,0,1,1,0,0,0,1,0,1,1,1,1,0,1,1,1,1,1};
+};
+#include "simplePointCalculator.txx"
+#endif
+
diff --git a/include/simplePointCalculator.txx b/include/simplePointCalculator.txx
new file mode 100644
index 0000000000000000000000000000000000000000..c5d2021612b4e33f9cea474291b92ca4d416490d
--- /dev/null
+++ b/include/simplePointCalculator.txx
@@ -0,0 +1,88 @@
+#ifndef simplePointCalculator_txx
+#define simplePointCalculator_txx
+
+#include "simplePointCalculator.h"
+
+template <class IteratorType, class LabelType>
+constexpr
+std::array<bool,256>
+simplePointCalculator<IteratorType,LabelType>::T4Table;
+
+template <class IteratorType, class LabelType>
+int
+simplePointCalculator<IteratorType,LabelType>
+::isSimpleLabel (IteratorType const& neighborhoodLabelsBegin, IteratorType const& neighborhoodLabelsEnd, LabelType const& evalLabel) const {
+ // read in table if evalLabel is simple wrt the neighborhoodLabels
+ // simpleKey is 8 bits representing the neighborhood in order
+ // b7 b6 b5
+ // b4 XX b3
+ // b2 b1 b0
+
+ unsigned char simpleKey = 0;
+ //read neighborhood backwards
+ IteratorType n = neighborhoodLabelsEnd-1;
+ for (; n > neighborhoodLabelsBegin-1; --n){
+ //Discard central label
+ if(n != neighborhoodLabelsBegin+4)
+ {
+ if (*n == evalLabel)
+ {
+ simpleKey++;
+ }
+ if(n != neighborhoodLabelsBegin)
+ {
+ simpleKey <<= 1;
+ }
+ }
+ }
+ return simplePointCalculator<IteratorType,LabelType>::T4Table[simpleKey];
+}
+
+template <class IteratorType, class LabelType>
+void
+simplePointCalculator<IteratorType,LabelType>
+::generateInitializerList(){
+ std::array<bool,8> ib;
+ std::array<int,256> ar;
+ for(int i = 0 ; i < 255 ; ++i)
+ {
+ //Convert number to bit table
+ for (int bit = 0; bit < 8 ; ++bit)
+ {
+ ib[bit] = (i >> bit) & 1;
+ }
+ int n=0;
+ for (int bit = 1; bit < 8 ; ++bit)
+ {
+ n += (ib[bit] == 0 && ib[bit-1] == 1);
+ }
+ if(ib[7] == 1 && ib[0] == 0)
+ {
+ n++;
+ }
+ ar[i] = (n == 1)&&(ib[1] || ib[3] || ib[5] || ib[7]);
+ }
+ ar[255]=1;
+ bool correct = true;
+ std::cout<<"{";
+ for (unsigned int i = 0 ; i < 256; i++)
+ {
+ std::bitset<8> newbits(i);
+ std::bitset<8> oldbits;
+ oldbits.set(0,newbits[0]);
+ oldbits.set(1,newbits[1]);
+ oldbits.set(2,newbits[2]);
+ oldbits.set(3,newbits[4]);
+ oldbits.set(4,newbits[7]);
+ oldbits.set(5,newbits[6]);
+ oldbits.set(6,newbits[5]);
+ oldbits.set(7,newbits[3]);
+ std::cout << ar[oldbits.to_ulong()] << (i < 255 ? "," : "");
+ correct &= (ar[oldbits.to_ulong()] == simplePointCalculator<IteratorType,LabelType>::T4Table[i]);
+ }
+ std::cout << "}" << std::endl;
+
+ std::cout << (correct ? "success" : "fail") << std::endl;
+}
+#endif
+
diff --git a/otb-module.cmake b/otb-module.cmake
new file mode 100644
index 0000000000000000000000000000000000000000..359443e38383d925aef87a2adae5430474c5e9f4
--- /dev/null
+++ b/otb-module.cmake
@@ -0,0 +1,11 @@
+set(DOCUMENTATION "Large Scale OTB Application for image superpixel segmentation using SLIC algorithm.")
+
+otb_module(OTBSLIC
+ DEPENDS
+ OTBCommon
+ OTBApplicationEngine
+ OTBMPI
+ OTBGdalAdapters
+ OTBGDAL
+
+ DESCRIPTION "Superpixel segmentation - SLIC algorithm")
diff --git a/otb-module.cmake~ b/otb-module.cmake~
new file mode 100644
index 0000000000000000000000000000000000000000..fd25fd5dd313089836b15df359b2fa1a561442f2
--- /dev/null
+++ b/otb-module.cmake~
@@ -0,0 +1,6 @@
+otb_module(OTBSLIC
+ DEPENDS
+ OTBCommon
+ OTBApplicationEngine
+
+ DESCRIPTION "Superpixel segmentation - SLIC algorithm")