otbLSMSSegmentation.cxx 29.3 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
/*
 * Copyright (C) 2005-2017 Centre National d'Etudes Spatiales (CNES)
 *
 * This file is part of Orfeo Toolbox
 *
 *     https://www.orfeo-toolbox.org/
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
20 21 22 23 24 25

#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "otbMultiChannelExtractROI.h"
#include "otbExtractROI.h"
#include "otbConnectedComponentMuParserFunctor.h"
26
#include "itkUnaryFunctorImageFilter.h"
27 28 29 30 31
#include "itkStatisticsImageFilter.h"
#include "itkChangeLabelImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkScalarConnectedComponentImageFilter.h"
#include "otbConcatenateVectorImageFilter.h"
32
#include "otbAffineFunctor.h"
33 34

#include "otbMultiToMonoChannelExtractROI.h"
35
#include "otbImportGeoInformationImageFilter.h"
36 37

#include <time.h>
38
#include <algorithm>
39 40 41 42 43 44 45

#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"

#include "otbStandardWriterWatcher.h"
#include <itksys/SystemTools.hxx>

46

47 48 49 50 51 52 53 54 55 56 57 58 59
namespace otb
{
namespace Wrapper
{
class LSMSSegmentation : public Application
{
public:
  typedef LSMSSegmentation Self;
  typedef Application Superclass;
  typedef itk::SmartPointer<Self> Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;
  itkNewMacro(Self);

60
  itkTypeMacro(LSMSSegmentation, otb::Application);
61

62 63 64
  typedef FloatVectorImageType              ImageType;
  typedef ImageType::InternalPixelType      ImagePixelType;
  typedef UInt32ImageType                   LabelImageType;
65
  typedef LabelImageType::InternalPixelType LabelImagePixelType;
66 67 68
  typedef otb::ImageFileReader<ImageType> ImageReaderType;
  typedef otb::ImageFileWriter<ImageType> ImageWriterType;
  typedef otb::ImageFileReader<LabelImageType> LabelImageReaderType;
OTB Bot's avatar
STYLE  
OTB Bot committed
69
  typedef otb::ImageFileWriter<LabelImageType> LabelImageWriterType;
70
  typedef otb::MultiChannelExtractROI <ImagePixelType,ImagePixelType > MultiChannelExtractROIFilterType;
OTB Bot's avatar
STYLE  
OTB Bot committed
71
  typedef otb::ExtractROI<LabelImagePixelType,LabelImagePixelType> ExtractROIFilterType;
72 73 74 75 76 77
  typedef otb::MultiToMonoChannelExtractROI<ImagePixelType,LabelImagePixelType> MultiToMonoChannelExtractROIFilterType;
  typedef otb::Functor::ConnectedComponentMuParserFunctor<ImageType::PixelType>  CCFunctorType;
  typedef itk::ConnectedComponentFunctorImageFilter<ImageType, LabelImageType, CCFunctorType, otb::Image<unsigned int> > CCFilterType;
  typedef itk::ScalarConnectedComponentImageFilter<LabelImageType, LabelImageType> ScalarCCFilterType;
  typedef itk::StatisticsImageFilter<LabelImageType> StatisticsImageFilterType;
  typedef itk::ChangeLabelImageFilter<LabelImageType,LabelImageType> ChangeLabelImageFilterType;
78
  typedef otb::ImportGeoInformationImageFilter<LabelImageType,ImageType> ImportGeoInformationImageFilterType;
79
  typedef itk::ImageRegionConstIterator<LabelImageType> LabelImageIterator;
80

OTB Bot's avatar
STYLE  
OTB Bot committed
81
  typedef otb::ConcatenateVectorImageFilter <ImageType,ImageType,ImageType> ConcatenateType;
82 83 84 85 86 87 88 89
  typedef otb::Functor::AffineFunctor<
    LabelImagePixelType,
    LabelImagePixelType,
    LabelImagePixelType>                      AffineFunctorType;
  typedef itk::UnaryFunctorImageFilter<
    LabelImageType,
    LabelImageType,
    AffineFunctorType>                        LabelShiftFilterType;
90

91 92
  LSMSSegmentation(): //m_FinalReader(),m_ImportGeoInformationFilter(),
    m_FilesToRemoveAfterExecute(),m_TmpDirCleanup(false){}
93

94
  ~LSMSSegmentation() override{}
95

OTB Bot's avatar
STYLE  
OTB Bot committed
96
private:
97 98
  // LabelImageReaderType::Pointer m_FinalReader;
  // ImportGeoInformationImageFilterType::Pointer m_ImportGeoInformationFilter;
Julien Michel's avatar
Julien Michel committed
99
  std::vector<std::string> m_FilesToRemoveAfterExecute;
100
  bool m_TmpDirCleanup;
101 102

  std::string CreateFileName(unsigned int row, unsigned int column, std::string label)
103
  {
104
    std::string outfname = GetParameterString("out");
105
    std::string tilesname = itksys::SystemTools::GetFilenameWithoutExtension(outfname);
106

107 108
    std::stringstream tileOut;
    tileOut<<tilesname<<"_"<<row<<"_"<<column<<"_"<<label<<".tif";
109

110 111 112
    std::vector<std::string> joins;
    if(IsParameterEnabled("tmpdir"))
      {
113
      std::string tmpdir = GetParameterString("tmpdir");
114

115 116 117 118 119
      if(tmpdir.size() > 1 && tmpdir[tmpdir.size()-1] != '/')
        {
        tmpdir.append("/");
        }
      joins.push_back(tmpdir);
120 121
      }
    joins.push_back(tileOut.str());
122

Julien Malik's avatar
STYLE  
Julien Malik committed
123
    std::string currentFile = itksys::SystemTools::JoinPath(joins);
124 125 126

    return currentFile;
  }
127

128
  std::string WriteTile(LabelImageType::Pointer img, unsigned int row, unsigned int column, std::string label)
OTB Bot's avatar
STYLE  
OTB Bot committed
129
  {
130
    std::string currentFile = CreateFileName(row,column,label);
131

132 133 134
    LabelImageWriterType::Pointer imageWriter = LabelImageWriterType::New();
    imageWriter->SetInput(img);
    imageWriter->SetFileName(currentFile);
OTB Bot's avatar
STYLE  
OTB Bot committed
135
    imageWriter->Update();
136 137 138

    return currentFile;
  }
139

140
  void RemoveFile(std::string tile)
141 142
  {
    // Cleanup
143
    if(GetParameterInt("cleanup"))
144 145
      {
        // Try to remove the geom file if existing
146
      std::string geomfile = tile.substr(0,tile.size() - itksys::SystemTools::GetFilenameExtension(tile).size()).append(".geom");
147

148
      if(itksys::SystemTools::FileExists(geomfile))
149
        {
150
        bool res = itksys::SystemTools::RemoveFile(geomfile);
151 152 153 154
        if (!res)
          {
          otbAppLogINFO(<<"Unable to remove file  "<<geomfile);
          }
155
        }
156
      if(itksys::SystemTools::FileExists(tile))
157
        {
158
        bool res = itksys::SystemTools::RemoveFile(tile);
159 160 161 162
        if (!res)
          {
          otbAppLogINFO(<<"Unable to remove file  "<<tile);
          }
163 164 165 166 167 168 169 170 171
        }
      }
  }

  std::string WriteVRTFile(unsigned int nbTilesX,unsigned int nbTilesY, unsigned long tileSizeX,unsigned long tileSizeY, unsigned long imageSizeX,unsigned long imageSizeY)
  {
    ImageType::Pointer imageIn = GetParameterImage("in");

    std::string outfname = GetParameterString("out");
172

173
    std::stringstream vrtfOut;
174
    vrtfOut<<itksys::SystemTools::GetFilenameWithoutExtension(outfname)<<".vrt";
175

176 177 178 179
    std::vector<std::string> joins;
    if(IsParameterEnabled("tmpdir"))
      {
      std::string tmpdir = GetParameterString("tmpdir");
180

181 182 183 184 185 186 187
      if(tmpdir.size() > 1 && tmpdir[tmpdir.size()-1] != '/')
        {
        tmpdir.append("/");
        }
      joins.push_back(tmpdir);
      }
    joins.push_back(vrtfOut.str());
188

189
    std::string vrtfname = itksys::SystemTools::JoinPath(joins);
190
    otbAppLogINFO(<<"Creating temporary vrt file: "<<vrtfname);
191

192
    std::ofstream ofs(vrtfname);
193 194 195 196 197

    ofs<<"<VRTDataset rasterXSize=\""<<imageSizeX<<"\" rasterYSize=\""<<imageSizeY<<"\">"<<std::endl;
    ofs<<"\t<VRTRasterBand dataType=\"UInt32\" band=\"1\">"<<std::endl;
    ofs<<"\t\t<ColorInterp>Gray</ColorInterp>"<<std::endl;

OTB Bot's avatar
STYLE  
OTB Bot committed
198
     for(unsigned int column = 0; column < nbTilesX; ++column)
199
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
200
      for(unsigned int row = 0; row < nbTilesY; ++row)
201 202
        {
        ofs<<"\t\t<SimpleSource>"<<std::endl;
203
        ofs<<"\t\t\t<SourceFilename relativeToVRT=\"1\">"<< itksys::SystemTools::GetFilenameName(CreateFileName(row,column,"FINAL")) <<"</SourceFilename>"<<std::endl;
204 205 206 207 208 209 210 211 212 213 214 215 216
        ofs<<"\t\t\t<SourceBand>1</SourceBand>"<<std::endl;
        ofs<<"\t\t\t<SrcRect xOff=\""<<0<<"\" yOff=\""<<0<<"\" xSize=\""<<tileSizeX<<"\" ySize=\""<<tileSizeY<<"\"/>"<<std::endl;
        ofs<<"\t\t\t<DstRect xOff=\""<<column*tileSizeX<<"\" yOff=\""<<row*tileSizeY<<"\" xSize=\""<<tileSizeX<<"\" ySize=\""<<tileSizeY<<"\"/>"<<std::endl;
        ofs<<"\t\t</SimpleSource>"<<std::endl;
        }
      }
     ofs<<"\t</VRTRasterBand>"<<std::endl;
     ofs<<"</VRTDataset>"<<std::endl;

    ofs.close();

    return vrtfname;
  }
217

218
  void DoInit() override
219 220
  {
    SetName("LSMSSegmentation");
221
    SetDescription("This application performs the second step of the exact Large-Scale Mean-Shift segmentation workflow (LSMS) [1].");
222

223
    SetDocName("Exact Large-Scale Mean-Shift segmentation, step 2");
Julien Michel's avatar
Julien Michel committed
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245
    SetDocLongDescription("This application will produce a labeled image where neighbor pixels"
                          " whose range distance is below range radius (and optionally spatial"
                          " distance below spatial radius) will be grouped together into the same"
                          " cluster. For large images one can use the tilesizex and tilesizey"
                          " parameters for tile-wise processing, with the guarantees of identical"
                          " results.\n\n"
                          "Filtered range image and spatial image should be created with the"
                          " MeanShiftSmoothing application outputs (fout and foutpos) [2], with"
                          " modesearch parameter disabled. If spatial image is not set, the"
                          " application will only process the range image and spatial radius"
                          " parameter will not be taken into account.\n\n"
                          "Please note that this application will generate a lot of temporary"
                          " files (as many as the number of tiles), and will therefore require"
                          " twice the size of the final result in term of disk space. The cleanup"
                          " option (activated by default) allows removing all temporary file as"
                          " soon as they are not needed anymore (if cleanup is activated, tmpdir"
                          " set and tmpdir does not exists before running the application, it will"
                          " be removed as well during cleanup). The tmpdir option allows defining"
                          " a directory where to write the temporary files.\n\n"
                          "Please also note that the output image type should be set to uint32 to"
                          " ensure that there are enough labels available.\n\n"
                          "The output of this application can be passed to the"
Victor Poughon's avatar
Victor Poughon committed
246
                          " LSMSSmallRegionsMerging [3] or LSMSVectorization [4] applications to"
Julien Michel's avatar
Julien Michel committed
247 248 249 250 251
                          " complete the LSMS workflow.");
    SetDocLimitations("This application is part of the Large-Scale Mean-Shift segmentation"
                      " workflow (LSMS) [1] and may not be suited for any other purpose. This"
                      " application is not compatible with in-memory connection since it does"
                      " its own internal streaming.");
252
    SetDocAuthors("David Youssefi");
253 254 255 256 257 258 259
    SetDocSeeAlso( "[1] Michel, J., Youssefi, D., & Grizonnet, M. (2015). Stable"
                   " mean-shift algorithm and its application to the segmentation of"
                   " arbitrarily large remote sensing images. IEEE Transactions on"
                   " Geoscience and Remote Sensing, 53(2), 952-964.\n"
                   "[2] MeanShiftSmoothing\n"
                   "[3] LSMSSmallRegionsMerging\n"
                   "[4] LSMSVectorization");
260 261 262 263
    AddDocTag(Tags::Segmentation);
    AddDocTag("LSMS");

    AddParameter(ParameterType_InputImage,  "in",    "Filtered image");
264 265 266
    SetParameterDescription( "in", "The filtered image, corresponding to the fout output parameter of the MeanShiftSmoothing application." );
    AddParameter(ParameterType_InputImage,  "inpos",    "Filtered position image");
    SetParameterDescription( "inpos", " The filtered position image, corresponding to the foutpos output parameter of the MeanShiftSmoothing application.");
267 268
    MandatoryOff("inpos");

269
    AddParameter(ParameterType_OutputImage, "out", "Output labeled Image");
Julien Michel's avatar
Julien Michel committed
270
    SetParameterDescription( "out", "This output contains the segmented image, where each pixel value is the unique integer label of the segment it belongs to. It is recommended to set the pixel type to uint32." );
271
    SetDefaultOutputPixelType("out",ImagePixelType_uint32);
272 273

    AddParameter(ParameterType_Float, "spatialr", "Spatial radius");
274
    SetParameterDescription("spatialr", "Threshold on Spatial distance to consider pixels in the same segment. A good value is half the spatial radius used in the MeanShiftSmoothing application (spatialr parameter).");
275 276 277
    SetDefaultParameterFloat("spatialr", 5);
    SetMinimumParameterFloatValue("spatialr", 0);
    MandatoryOff("spatialr");
278 279
    
    AddParameter(ParameterType_Float, "ranger", "Range radius");
280
    SetParameterDescription("ranger", "Threshold on spectral signature euclidean distance (expressed in radiometry unit) to consider pixels in the same segment. A good value is half the range radius used in the MeanShiftSmoothing application (ranger parameter).");
281 282 283
    SetDefaultParameterFloat("ranger", 15);
    SetMinimumParameterFloatValue("ranger", 0);
    MandatoryOff("ranger");
284

285 286
    AddParameter(ParameterType_Int, "minsize", "Minimum Segment Size");
    SetParameterDescription("minsize", "Minimum Segment Size. If, after the segmentation, a segment is of size lower than this criterion, the segment is discarded.");
287 288 289
    SetDefaultParameterInt("minsize", 0);
    SetMinimumParameterIntValue("minsize", 0);
    MandatoryOff("minsize");
290

291
    AddParameter(ParameterType_Int, "tilesizex", "Size of tiles in pixel (X-axis)");
292
    SetParameterDescription("tilesizex", "Size of tiles along the X-axis for tile-wise processing.");
293 294 295 296
    SetDefaultParameterInt("tilesizex", 500);
    SetMinimumParameterIntValue("tilesizex", 1);

    AddParameter(ParameterType_Int, "tilesizey", "Size of tiles in pixel (Y-axis)");
297
    SetParameterDescription("tilesizey", "Size of tiles along the Y-axis for tile-wise processing.");
298 299
    SetDefaultParameterInt("tilesizey", 500);
    SetMinimumParameterIntValue("tilesizey", 1);
300 301

    AddParameter(ParameterType_Directory,"tmpdir","Directory where to write temporary files");
Julien Malik's avatar
Julien Malik committed
302
    SetParameterDescription("tmpdir","This applications need to write temporary files for each tile. This parameter allows choosing the path where to write those files. If disabled, the current path will be used.");
303 304 305
    MandatoryOff("tmpdir");
    DisableParameter("tmpdir");

306
    AddParameter(ParameterType_Bool,"cleanup","Temporary files cleaning");
307
    SetParameterDescription("cleanup","If activated, the application will try to remove all temporary files it created.");
308
    SetParameterInt("cleanup",1);
309 310 311 312 313 314

    // Doc example parameter settings
    SetDocExampleParameterValue("in","smooth.tif");
    SetDocExampleParameterValue("inpos","position.tif");
    SetDocExampleParameterValue("out","segmentation.tif");
    SetDocExampleParameterValue("spatialr","5");
315
    SetDocExampleParameterValue("ranger","15");
316
    SetDocExampleParameterValue("minsize","0");
317 318
    SetDocExampleParameterValue("tilesizex","256");
    SetDocExampleParameterValue("tilesizey","256");
319

320
    SetOfficialDocLink();
321 322
  }

323
  void DoUpdateParameters() override
OTB Bot's avatar
STYLE  
OTB Bot committed
324
  {
325 326
  }

327
  void DoExecute() override
328
  {
329 330
    m_FilesToRemoveAfterExecute.clear();

331
    clock_t tic = clock();
332

333 334
    const float ranger         = GetParameterFloat("ranger");
    const float spatialr       = GetParameterFloat("spatialr");
335

336
    unsigned int minRegionSize = GetParameterInt("minsize");
337

338 339
    unsigned long sizeTilesX   = GetParameterInt("tilesizex");
    unsigned long sizeTilesY   = GetParameterInt("tilesizey");
340 341 342 343 344


    // Ensure that temporary directory exists if activated:
    if(IsParameterEnabled("tmpdir"))
      {
345
      if(!itksys::SystemTools::FileExists(GetParameterString("tmpdir")))
346 347 348
        {
        m_TmpDirCleanup = true;
        }
349
      otbAppLogINFO(<<"Temporary directory "<<GetParameterString("tmpdir")<<" will be used");
350
      itksys::SystemTools::MakeDirectory(GetParameterString("tmpdir"));
351 352 353 354 355 356 357 358
      }

    //Three steps :
    // 1-Tiles segmentation
    // 2-Tiles relabelling
    // 3-Minimal size region suppression

    ImageType::Pointer spatialIn;
359

360
    if(HasValue("inpos"))
361
      {
362
      spatialIn = GetParameterImage("inpos");
363 364 365 366 367
      }

    //Acquisition of the input image dimensions
    ImageType::Pointer imageIn = GetParameterImage("in");
    imageIn->UpdateOutputInformation();
368 369 370 371

    unsigned long sizeImageX = imageIn->GetLargestPossibleRegion().GetSize()[0];
    unsigned long sizeImageY = imageIn->GetLargestPossibleRegion().GetSize()[1];
    unsigned int nbComp      = imageIn->GetNumberOfComponentsPerPixel();
372

373 374 375 376 377
    unsigned int nbTilesX = sizeImageX/sizeTilesX + (sizeImageX%sizeTilesX > 0 ? 1 : 0);
    unsigned int nbTilesY = sizeImageY/sizeTilesY + (sizeImageY%sizeTilesY > 0 ? 1 : 0);

    otbAppLogINFO(<<"Number of tiles: "<<nbTilesX<<" x "<<nbTilesY);

378
    unsigned long regionCount = 0;
379

380 381 382
    //Segmentation by the connected component per tile and label
    //shifting
    otbAppLogINFO(<<"Tile shifting ...");
383

OTB Bot's avatar
STYLE  
OTB Bot committed
384 385 386
    for(unsigned int row = 0; row < nbTilesY; ++row)
      for(unsigned int column = 0; column < nbTilesX; ++column)
        {
387 388 389
        // Compute extraction parameters
        unsigned long startX = column*sizeTilesX;
        unsigned long startY = row*sizeTilesY;
390 391
        unsigned long sizeX = std::min(sizeTilesX+1,sizeImageX-startX+1);
        unsigned long sizeY = std::min(sizeTilesY+1,sizeImageY-startY+1);
392

393 394 395 396 397 398 399 400 401 402
        //Tiles extraction of :
        //- the input image (filtering image)
        MultiChannelExtractROIFilterType::Pointer extractROIFilter = MultiChannelExtractROIFilterType::New();
        extractROIFilter->SetInput(imageIn);
        extractROIFilter->SetStartX(startX);
        extractROIFilter->SetStartY(startY);
        extractROIFilter->SetSizeX(sizeX);
        extractROIFilter->SetSizeY(sizeY);
        extractROIFilter->Update();

403
        //- the spatial image (final positions) if available
404 405
        MultiChannelExtractROIFilterType::Pointer extractROIFilter2 = MultiChannelExtractROIFilterType::New();
        ConcatenateType::Pointer concat = ConcatenateType::New();
406 407 408 409 410 411 412 413 414 415
        CCFilterType::Pointer ccFilter = CCFilterType::New();

        if(HasValue("inpos"))
          {
          extractROIFilter2->SetInput(spatialIn);
          extractROIFilter2->SetStartX(startX);
          extractROIFilter2->SetStartY(startY);
          extractROIFilter2->SetSizeX(sizeX);
          extractROIFilter2->SetSizeY(sizeY);
          extractROIFilter2->Update();
416

417 418 419 420
          //Concatenation of the two input images
          concat->SetInput1(extractROIFilter->GetOutput());
          concat->SetInput2(extractROIFilter2->GetOutput());
          concat->Update();
421

422 423 424 425
          ccFilter->SetInput(concat->GetOutput());
          }
        else
          {
OTB Bot's avatar
STYLE  
OTB Bot committed
426
         ccFilter->SetInput(extractROIFilter->GetOutput());
427
          }
428

429 430 431
        //Expression 1 : radiometric distance < ranger
        std::stringstream expr;
        expr<<"sqrt((p1b1-p2b1)*(p1b1-p2b1)";
OTB Bot's avatar
STYLE  
OTB Bot committed
432
        for(unsigned int i=1; i<nbComp; i++)
433
          expr<<"+(p1b"<<i+1<<"-p2b"<<i+1<<")*(p1b"<<i+1<<"-p2b"<<i+1<<")";
434
        expr<<")"<<"<"<<ranger;
435

436 437 438
        if(HasValue("inpos"))
          {
          //Expression 2 : final positions < spatialr
439
          expr<<" and sqrt((p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")*(p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")+";
OTB Bot's avatar
STYLE  
OTB Bot committed
440
          expr<<"(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<")*(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<"))"<<"<"<<spatialr;
441
          }
442

443 444
        //Segmentation
        ccFilter->GetFunctor().SetExpression(expr.str());
OTB Bot's avatar
STYLE  
OTB Bot committed
445
        ccFilter->Update();
446

447
        //Shifting
448 449 450 451 452
        LabelShiftFilterType::Pointer labelShiftFilter = LabelShiftFilterType::New();
        labelShiftFilter->SetInput(ccFilter->GetOutput());
        labelShiftFilter->GetFunctor().SetA(1);
        labelShiftFilter->GetFunctor().SetB(regionCount);
        labelShiftFilter->Update();
453 454 455

        //Maximum label calculation for the shifting
        StatisticsImageFilterType::Pointer stats = StatisticsImageFilterType::New();
OTB Bot's avatar
STYLE  
OTB Bot committed
456
        stats->SetInput(ccFilter->GetOutput());
457 458 459
        stats->Update();
        regionCount+=stats->GetMaximum();

460
        std::string filename = WriteTile(labelShiftFilter->GetOutput(),row,column,"SEG");
461
        }
462

463

464 465
    // Step 2: create the look-up table for all overlaps
    otbAppLogINFO(<<"LUT creation ...");
OTB Bot's avatar
STYLE  
OTB Bot committed
466
    std::vector<LabelImagePixelType> LUT;
467 468 469 470
    LUT.clear();
    LUT.resize(regionCount+1);
    for(LabelImagePixelType curLabel = 1; curLabel <= regionCount; ++curLabel)
      LUT[curLabel] = curLabel;
471

OTB Bot's avatar
STYLE  
OTB Bot committed
472
    for(unsigned int row = 0; row < nbTilesY; row++)
473
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
474
      for(unsigned int column = 0; column < nbTilesX; column++)
475 476 477
        {
        unsigned long startX = column*sizeTilesX;
        unsigned long startY = row*sizeTilesY;
478 479
        unsigned long sizeX = std::min(sizeTilesX+1,sizeImageX-startX+1);
        unsigned long sizeY = std::min(sizeTilesY+1,sizeImageY-startY+1);
480

481
        std::string tileIn = CreateFileName(row,column,"SEG");
482

483 484 485 486
        // Read current tile
        LabelImageReaderType::Pointer tileInReader = LabelImageReaderType::New();
        tileInReader->SetFileName(tileIn);
        tileInReader->Update();
487

488 489 490 491
        // Analyse intersection between in and up tiles
        if(row>0)
          {
          std::string tileUp = CreateFileName(row-1,column,"SEG");
Julien Malik's avatar
STYLE  
Julien Malik committed
492 493
          LabelImageReaderType::Pointer tileUpReader = LabelImageReaderType::New();
          tileUpReader->SetFileName(tileUp);
494 495 496 497 498 499 500 501
          tileUpReader->Update();

          LabelImageType::IndexType pixelIndexIn;
          LabelImageType::IndexType pixelIndexUp;

          pixelIndexIn[1] = 0;
          pixelIndexUp[1] = sizeTilesY;

OTB Bot's avatar
STYLE  
OTB Bot committed
502
          for(pixelIndexIn[0]=0; pixelIndexIn[0]<static_cast<long>(sizeX-1); ++pixelIndexIn[0])
503 504 505 506
            {
            pixelIndexUp[0] = pixelIndexIn[0];

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
507
           while(LUT[curCanLabel] != curCanLabel)
508 509 510
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
511
           LabelImagePixelType adjCanLabel = tileUpReader->GetOutput()->GetPixel(pixelIndexUp);
512

OTB Bot's avatar
STYLE  
OTB Bot committed
513
           while(LUT[adjCanLabel] != adjCanLabel)
514 515 516
              {
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
517
           if(curCanLabel < adjCanLabel)
518 519 520
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
521
           else
522 523 524 525 526
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
527

528 529 530 531
        // Analyse intersection between in and left tiles
         if(column>0)
          {
          std::string tileLeft = CreateFileName(row,column-1,"SEG");
Julien Malik's avatar
STYLE  
Julien Malik committed
532 533
          LabelImageReaderType::Pointer tileLeftReader = LabelImageReaderType::New();
          tileLeftReader->SetFileName(tileLeft);
534 535 536 537 538 539 540 541
          tileLeftReader->Update();

          LabelImageType::IndexType pixelIndexIn;
          LabelImageType::IndexType pixelIndexUp;

          pixelIndexIn[0] = 0;
          pixelIndexUp[0] = sizeTilesX;

OTB Bot's avatar
STYLE  
OTB Bot committed
542
          for(pixelIndexIn[1]=0; pixelIndexIn[1]<static_cast<long>(sizeY-1); ++pixelIndexIn[1])
543 544 545 546
            {
            pixelIndexUp[1] = pixelIndexIn[1];

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
547
           while(LUT[curCanLabel] != curCanLabel)
548 549 550
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
551 552 553
           LabelImagePixelType adjCanLabel = tileLeftReader->GetOutput()->GetPixel(pixelIndexUp);
           while(LUT[adjCanLabel] != adjCanLabel)
              {
554 555
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
556
           if(curCanLabel < adjCanLabel)
557 558 559
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
560
           else
561 562 563 564 565
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
566 567
        }
      }
568

569
    // Reduce LUT to canonical labels
570
    for(LabelImagePixelType label = 1; label < regionCount+1; ++label)
OTB Bot's avatar
STYLE  
OTB Bot committed
571
      {
572
      LabelImagePixelType can = label;
OTB Bot's avatar
STYLE  
OTB Bot committed
573
      while(LUT[can] != can)
574 575 576
        {
        can = LUT[can];
        }
577 578
      LUT[label] = can;
      }
579
    otbAppLogINFO(<<"LUT size: "<<LUT.size()<<" segments");
580

581 582 583
    // These variables will be used to estimate the size of each
    // region on the flow
    std::vector<unsigned long> sizePerRegion(regionCount+1,0);
584

585
    otbAppLogINFO(<<"Tiles relabelisation ...");
OTB Bot's avatar
STYLE  
OTB Bot committed
586
    for(unsigned int column = 0; column < nbTilesX; ++column)
587
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
588
      for(unsigned int row = 0; row < nbTilesY; ++row)
589
        {
Julien Malik's avatar
STYLE  
Julien Malik committed
590
        unsigned long startX = column*sizeTilesX;
591
        unsigned long startY = row*sizeTilesY;
592 593
        unsigned long sizeX = std::min(sizeTilesX,sizeImageX-startX);
        unsigned long sizeY = std::min(sizeTilesY,sizeImageY-startY);
594 595 596 597 598 599 600 601

        std::string tileIn = CreateFileName(row,column,"SEG");

        LabelImageReaderType::Pointer readerIn = LabelImageReaderType::New();
        readerIn->SetFileName(tileIn);

        // Remove extra margin now that lut is built
        ExtractROIFilterType::Pointer labelImage = ExtractROIFilterType::New();
Julien Malik's avatar
STYLE  
Julien Malik committed
602 603 604 605 606
        labelImage->SetInput(readerIn->GetOutput());
        labelImage->SetStartX(0);
        labelImage->SetStartY(0);
        labelImage->SetSizeX(sizeX);
        labelImage->SetSizeY(sizeY);
607

608 609
        // Relabel tile according to look-up table
        ChangeLabelImageFilterType::Pointer changeLabel = ChangeLabelImageFilterType::New();
Julien Malik's avatar
STYLE  
Julien Malik committed
610
        changeLabel->SetInput(labelImage->GetOutput());
611

612
        // Fill LUT
OTB Bot's avatar
STYLE  
OTB Bot committed
613
        for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
614
          {
Julien Malik's avatar
STYLE  
Julien Malik committed
615
          if(label!=LUT[label])
616 617 618 619 620 621 622 623 624 625 626
            {
            changeLabel->SetChange(label,LUT[label]);
            }
          }

        // Here we need to update the filter to be able to update
        // region sizes table
        changeLabel->Update();

        // Update region sizes
        LabelImageIterator it( changeLabel->GetOutput(), changeLabel->GetOutput()->GetLargestPossibleRegion());
627
        for (it.GoToBegin(); !it.IsAtEnd(); ++it)
628
          {
629
          sizePerRegion[it.Value()]+=1;
630
          }
631 632 633 634 635

        // Write tile
        WriteTile(changeLabel->GetOutput(),row,column,"RELAB");

        // Remove previous tile (not needed anymore)
636
        readerIn = nullptr; // release the input file
637
        RemoveFile(tileIn);
638
        }
639
      }
640

641 642
    // Clear lut, we do not need it anymore
    LUT.clear();
643

644 645
    unsigned int smallCount = 0;

646 647 648 649 650
      // Create the LUT to filter small regions and assign min labels
      otbAppLogINFO(<<"Small regions pruning ...");
      LabelImagePixelType newLab=1;
      std::vector<LabelImagePixelType> newLabels(regionCount+1,0);
      for(LabelImagePixelType curLabel = 1; curLabel <= regionCount; ++curLabel)
651
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
652
        if(sizePerRegion[curLabel]<minRegionSize)
653
          {
654
          newLabels[curLabel]=0;
655
          ++smallCount;
656
          }
657
        else
658
          {
659 660
          newLabels[curLabel]=newLab;
          newLab+=1;
661 662 663
          }
        }

664
      otbAppLogINFO(<<smallCount<<" small regions will be removed");
665

666 667
      // Clear sizePerRegion, we do not need it anymore
      sizePerRegion.clear();
668

OTB Bot's avatar
STYLE  
OTB Bot committed
669
      for(unsigned int column = 0; column < nbTilesX; ++column)
670
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
671
        for(unsigned int row = 0; row < nbTilesY; ++row)
672 673 674 675 676 677 678 679
          {
          std::string tileIn = CreateFileName(row,column,"RELAB");

          LabelImageReaderType::Pointer readerIn = LabelImageReaderType::New();
          readerIn->SetFileName(tileIn);

          ChangeLabelImageFilterType::Pointer changeLabel = ChangeLabelImageFilterType::New();
          changeLabel->SetInput(readerIn->GetOutput());
OTB Bot's avatar
STYLE  
OTB Bot committed
680
          for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
681 682 683 684 685 686
            {
            if(label != newLabels[label])
              {
              changeLabel->SetChange(label,newLabels[label]);
              }
            }
687

688
          // Write the relabeled tile
689
          std::string tmpfile = WriteTile(changeLabel->GetOutput(),row,column,"FINAL");
690 691
          m_FilesToRemoveAfterExecute.push_back(tmpfile);

692
          // Clean previous tiles (not needed anymore)
693
          readerIn = nullptr; // release the input file
694
          RemoveFile(tileIn);
695
          }
696
        }
697 698 699

      // Clear newLabels, we do not need it anymore
      newLabels.clear();
700

701 702
      // Here we write a temporary vrt file that will be used to
      // stitch together all the tiles
Julien Michel's avatar
Julien Michel committed
703
      std::string vrtfile = WriteVRTFile(nbTilesX,nbTilesY,sizeTilesX,sizeTilesY,sizeImageX,sizeImageY);
704

705 706
      m_FilesToRemoveAfterExecute.push_back(vrtfile);

707
      clock_t toc = clock();
708

709
      otbAppLogINFO(<<"Elapsed time: "<<(double)(toc - tic) / CLOCKS_PER_SEC<<" seconds");
710

711
      // Final writing
712 713
      LabelImageReaderType::Pointer finalReader = LabelImageReaderType::New();
      finalReader->SetFileName(vrtfile);
714

715
      ImportGeoInformationImageFilterType::Pointer 
716
        importGeoInformationFilter = 
717
        ImportGeoInformationImageFilterType::New();
718 719
      importGeoInformationFilter->SetInput(finalReader->GetOutput());
      importGeoInformationFilter->SetSource(imageIn);
720

721
      SetParameterOutputImage("out",importGeoInformationFilter->GetOutput());
722
      RegisterPipeline();
723
  }
724

725
  void AfterExecuteAndWriteOutputs() override
726
  {
727
    // Release input files
728
    // finalReader = nullptr;
729

730
    if(GetParameterInt("cleanup"))
731 732 733 734
      {
      otbAppLogINFO(<<"Final clean-up ...");

      for(std::vector<std::string>::iterator it = m_FilesToRemoveAfterExecute.begin();
OTB Bot's avatar
STYLE  
OTB Bot committed
735
          it!=m_FilesToRemoveAfterExecute.end(); ++it)
736 737 738 739 740 741 742
        {
        RemoveFile(*it);
        }

      if(IsParameterEnabled("tmpdir") && m_TmpDirCleanup)
        {
        otbAppLogINFO(<<"Removing tmp directory "<<GetParameterString("tmpdir")<<", since it has been created by the application");
743
        itksys::SystemTools::RemoveADirectory(GetParameterString("tmpdir"));
744 745 746 747
        }
      }

    m_FilesToRemoveAfterExecute.clear();
748
    m_TmpDirCleanup = false;
749
  }
750 751 752 753 754
};
}
}

OTB_APPLICATION_EXPORT(otb::Wrapper::LSMSSegmentation)