otbLSMSSegmentation.cxx 28.7 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 38 39 40 41 42 43 44 45

#include <time.h>
#include <vcl_algorithm.h>

#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
  LSMSSegmentation(): m_FinalReader(),m_ImportGeoInformationFilter(),m_FilesToRemoveAfterExecute(),m_TmpDirCleanup(false){}
92

93
  ~LSMSSegmentation() ITK_OVERRIDE{}
94

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

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

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

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

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

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

    return currentFile;
  }
126

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

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

    return currentFile;
  }
138

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

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

  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");
171

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

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

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

188
    std::string vrtfname = itksys::SystemTools::JoinPath(joins);
189
    otbAppLogINFO(<<"Creating temporary vrt file: "<<vrtfname);
190 191 192 193 194 195 196

    std::ofstream ofs(vrtfname.c_str());

    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
197
     for(unsigned int column = 0; column < nbTilesX; ++column)
198
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
199
      for(unsigned int row = 0; row < nbTilesY; ++row)
200 201
        {
        ofs<<"\t\t<SimpleSource>"<<std::endl;
202
        ofs<<"\t\t\t<SourceFilename relativeToVRT=\"1\">"<< itksys::SystemTools::GetFilenameName(CreateFileName(row,column,"FINAL")) <<"</SourceFilename>"<<std::endl;
203 204 205 206 207 208 209 210 211 212 213 214 215
        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;
  }
216

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

222
    SetDocName("Exact Large-Scale Mean-Shift segmentation, step 2");
223 224 225 226 227 228
    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 LSMSSmallRegionMerging [3] or LSMSVectorization [4] to 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.");
229
    SetDocAuthors("David Youssefi");
230 231 232 233 234 235 236
    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");
237 238 239 240
    AddDocTag(Tags::Segmentation);
    AddDocTag("LSMS");

    AddParameter(ParameterType_InputImage,  "in",    "Filtered image");
241 242 243
    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.");
244 245
    MandatoryOff("inpos");

246 247
    AddParameter(ParameterType_OutputImage, "out", "Output labeled Image");
    SetParameterDescription( "out", "This output contains the segmented image, where each pixel value is the unique label of the segment it belongs to. It is recommended to set the pixel type to uint32." );
248
    SetDefaultOutputPixelType("out",ImagePixelType_uint32);
249 250

    AddParameter(ParameterType_Float, "spatialr", "Spatial radius");
251
    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).");
252 253 254
    SetDefaultParameterFloat("spatialr", 5);
    SetMinimumParameterFloatValue("spatialr", 0);
    MandatoryOff("spatialr");
255 256
    
    AddParameter(ParameterType_Float, "ranger", "Range radius");
257
    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).");
258 259 260
    SetDefaultParameterFloat("ranger", 15);
    SetMinimumParameterFloatValue("ranger", 0);
    MandatoryOff("ranger");
261

262 263
    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.");
264 265 266
    SetDefaultParameterInt("minsize", 0);
    SetMinimumParameterIntValue("minsize", 0);
    MandatoryOff("minsize");
267

268
    AddParameter(ParameterType_Int, "tilesizex", "Size of tiles in pixel (X-axis)");
269
    SetParameterDescription("tilesizex", "Size of tiles along the X-axis for tile-wise processing.");
270 271 272 273
    SetDefaultParameterInt("tilesizex", 500);
    SetMinimumParameterIntValue("tilesizex", 1);

    AddParameter(ParameterType_Int, "tilesizey", "Size of tiles in pixel (Y-axis)");
274
    SetParameterDescription("tilesizey", "Size of tiles along the Y-axis for tile-wise processing.");
275 276
    SetDefaultParameterInt("tilesizey", 500);
    SetMinimumParameterIntValue("tilesizey", 1);
277 278

    AddParameter(ParameterType_Directory,"tmpdir","Directory where to write temporary files");
Julien Malik's avatar
Julien Malik committed
279
    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.");
280 281 282
    MandatoryOff("tmpdir");
    DisableParameter("tmpdir");

283 284
    AddParameter(ParameterType_Empty,"cleanup","Temporary files cleaning");
    EnableParameter("cleanup");
285
    SetParameterDescription("cleanup","If activated, the application will try to remove all temporary files it created.");
286 287 288 289 290 291 292
    MandatoryOff("cleanup");

    // Doc example parameter settings
    SetDocExampleParameterValue("in","smooth.tif");
    SetDocExampleParameterValue("inpos","position.tif");
    SetDocExampleParameterValue("out","segmentation.tif");
    SetDocExampleParameterValue("spatialr","5");
293
    SetDocExampleParameterValue("ranger","15");
294
    SetDocExampleParameterValue("minsize","0");
295 296
    SetDocExampleParameterValue("tilesizex","256");
    SetDocExampleParameterValue("tilesizey","256");
297

298
    SetOfficialDocLink();
299 300
  }

301
  void DoUpdateParameters() ITK_OVERRIDE
OTB Bot's avatar
STYLE  
OTB Bot committed
302
  {
303 304
  }

305
  void DoExecute() ITK_OVERRIDE
306
  {
307 308
    m_FilesToRemoveAfterExecute.clear();

309
    clock_t tic = clock();
310

311 312
    const float ranger         = GetParameterFloat("ranger");
    const float spatialr       = GetParameterFloat("spatialr");
313

314
    unsigned int minRegionSize = GetParameterInt("minsize");
315

316 317
    unsigned long sizeTilesX   = GetParameterInt("tilesizex");
    unsigned long sizeTilesY   = GetParameterInt("tilesizey");
318 319 320 321 322


    // Ensure that temporary directory exists if activated:
    if(IsParameterEnabled("tmpdir"))
      {
323 324 325 326
      if(!itksys::SystemTools::FileExists(GetParameterString("tmpdir").c_str()))
        {
        m_TmpDirCleanup = true;
        }
327 328 329 330 331 332 333 334 335 336
      otbAppLogINFO(<<"Temporary directory "<<GetParameterString("tmpdir")<<" will be used");
      itksys::SystemTools::MakeDirectory(GetParameterString("tmpdir").c_str());
      }

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

    ImageType::Pointer spatialIn;
337

338
    if(HasValue("inpos"))
339
      {
340
      spatialIn = GetParameterImage("inpos");
341 342 343 344 345
      }

    //Acquisition of the input image dimensions
    ImageType::Pointer imageIn = GetParameterImage("in");
    imageIn->UpdateOutputInformation();
346 347 348 349

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

351 352 353 354 355
    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);

356
    unsigned long regionCount = 0;
357

358 359 360
    //Segmentation by the connected component per tile and label
    //shifting
    otbAppLogINFO(<<"Tile shifting ...");
361

OTB Bot's avatar
STYLE  
OTB Bot committed
362 363 364
    for(unsigned int row = 0; row < nbTilesY; ++row)
      for(unsigned int column = 0; column < nbTilesX; ++column)
        {
365 366 367 368 369
        // Compute extraction parameters
        unsigned long startX = column*sizeTilesX;
        unsigned long startY = row*sizeTilesY;
        unsigned long sizeX = vcl_min(sizeTilesX+1,sizeImageX-startX+1);
        unsigned long sizeY = vcl_min(sizeTilesY+1,sizeImageY-startY+1);
370

371 372 373 374 375 376 377 378 379 380
        //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();

381
        //- the spatial image (final positions) if available
382 383
        MultiChannelExtractROIFilterType::Pointer extractROIFilter2 = MultiChannelExtractROIFilterType::New();
        ConcatenateType::Pointer concat = ConcatenateType::New();
384 385 386 387 388 389 390 391 392 393
        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();
394

395 396 397 398
          //Concatenation of the two input images
          concat->SetInput1(extractROIFilter->GetOutput());
          concat->SetInput2(extractROIFilter2->GetOutput());
          concat->Update();
399

400 401 402 403
          ccFilter->SetInput(concat->GetOutput());
          }
        else
          {
OTB Bot's avatar
STYLE  
OTB Bot committed
404
         ccFilter->SetInput(extractROIFilter->GetOutput());
405
          }
406

407 408 409
        //Expression 1 : radiometric distance < ranger
        std::stringstream expr;
        expr<<"sqrt((p1b1-p2b1)*(p1b1-p2b1)";
OTB Bot's avatar
STYLE  
OTB Bot committed
410
        for(unsigned int i=1; i<nbComp; i++)
411
          expr<<"+(p1b"<<i+1<<"-p2b"<<i+1<<")*(p1b"<<i+1<<"-p2b"<<i+1<<")";
412
        expr<<")"<<"<"<<ranger;
413

414 415 416
        if(HasValue("inpos"))
          {
          //Expression 2 : final positions < spatialr
417
          expr<<" and sqrt((p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")*(p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")+";
OTB Bot's avatar
STYLE  
OTB Bot committed
418
          expr<<"(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<")*(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<"))"<<"<"<<spatialr;
419
          }
420

421 422
        //Segmentation
        ccFilter->GetFunctor().SetExpression(expr.str());
OTB Bot's avatar
STYLE  
OTB Bot committed
423
        ccFilter->Update();
424

425
        //Shifting
426 427 428 429 430
        LabelShiftFilterType::Pointer labelShiftFilter = LabelShiftFilterType::New();
        labelShiftFilter->SetInput(ccFilter->GetOutput());
        labelShiftFilter->GetFunctor().SetA(1);
        labelShiftFilter->GetFunctor().SetB(regionCount);
        labelShiftFilter->Update();
431 432 433

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

438
        std::string filename = WriteTile(labelShiftFilter->GetOutput(),row,column,"SEG");
439
        }
440

441

442 443
    // Step 2: create the look-up table for all overlaps
    otbAppLogINFO(<<"LUT creation ...");
OTB Bot's avatar
STYLE  
OTB Bot committed
444
    std::vector<LabelImagePixelType> LUT;
445 446 447 448
    LUT.clear();
    LUT.resize(regionCount+1);
    for(LabelImagePixelType curLabel = 1; curLabel <= regionCount; ++curLabel)
      LUT[curLabel] = curLabel;
449

OTB Bot's avatar
STYLE  
OTB Bot committed
450
    for(unsigned int row = 0; row < nbTilesY; row++)
451
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
452
      for(unsigned int column = 0; column < nbTilesX; column++)
453 454 455 456 457
        {
        unsigned long startX = column*sizeTilesX;
        unsigned long startY = row*sizeTilesY;
        unsigned long sizeX = vcl_min(sizeTilesX+1,sizeImageX-startX+1);
        unsigned long sizeY = vcl_min(sizeTilesY+1,sizeImageY-startY+1);
458

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

461 462 463 464
        // Read current tile
        LabelImageReaderType::Pointer tileInReader = LabelImageReaderType::New();
        tileInReader->SetFileName(tileIn);
        tileInReader->Update();
465

466 467 468 469
        // 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
470 471
          LabelImageReaderType::Pointer tileUpReader = LabelImageReaderType::New();
          tileUpReader->SetFileName(tileUp);
472 473 474 475 476 477 478 479
          tileUpReader->Update();

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

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

OTB Bot's avatar
STYLE  
OTB Bot committed
480
          for(pixelIndexIn[0]=0; pixelIndexIn[0]<static_cast<long>(sizeX-1); ++pixelIndexIn[0])
481 482 483 484
            {
            pixelIndexUp[0] = pixelIndexIn[0];

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
485
           while(LUT[curCanLabel] != curCanLabel)
486 487 488
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
489
           LabelImagePixelType adjCanLabel = tileUpReader->GetOutput()->GetPixel(pixelIndexUp);
490

OTB Bot's avatar
STYLE  
OTB Bot committed
491
           while(LUT[adjCanLabel] != adjCanLabel)
492 493 494
              {
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
495
           if(curCanLabel < adjCanLabel)
496 497 498
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
499
           else
500 501 502 503 504
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
505

506 507 508 509
        // 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
510 511
          LabelImageReaderType::Pointer tileLeftReader = LabelImageReaderType::New();
          tileLeftReader->SetFileName(tileLeft);
512 513 514 515 516 517 518 519
          tileLeftReader->Update();

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

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

OTB Bot's avatar
STYLE  
OTB Bot committed
520
          for(pixelIndexIn[1]=0; pixelIndexIn[1]<static_cast<long>(sizeY-1); ++pixelIndexIn[1])
521 522 523 524
            {
            pixelIndexUp[1] = pixelIndexIn[1];

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
525
           while(LUT[curCanLabel] != curCanLabel)
526 527 528
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
529 530 531
           LabelImagePixelType adjCanLabel = tileLeftReader->GetOutput()->GetPixel(pixelIndexUp);
           while(LUT[adjCanLabel] != adjCanLabel)
              {
532 533
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
534
           if(curCanLabel < adjCanLabel)
535 536 537
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
538
           else
539 540 541 542 543
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
544 545
        }
      }
546

547
    // Reduce LUT to canonical labels
548
    for(LabelImagePixelType label = 1; label < regionCount+1; ++label)
OTB Bot's avatar
STYLE  
OTB Bot committed
549
      {
550
      LabelImagePixelType can = label;
OTB Bot's avatar
STYLE  
OTB Bot committed
551
      while(LUT[can] != can)
552 553 554
        {
        can = LUT[can];
        }
555 556
      LUT[label] = can;
      }
557
    otbAppLogINFO(<<"LUT size: "<<LUT.size()<<" segments");
558

559 560 561
    // These variables will be used to estimate the size of each
    // region on the flow
    std::vector<unsigned long> sizePerRegion(regionCount+1,0);
562

563
    otbAppLogINFO(<<"Tiles relabelisation ...");
OTB Bot's avatar
STYLE  
OTB Bot committed
564
    for(unsigned int column = 0; column < nbTilesX; ++column)
565
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
566
      for(unsigned int row = 0; row < nbTilesY; ++row)
567
        {
Julien Malik's avatar
STYLE  
Julien Malik committed
568
        unsigned long startX = column*sizeTilesX;
569
        unsigned long startY = row*sizeTilesY;
Julien Malik's avatar
STYLE  
Julien Malik committed
570
        unsigned long sizeX = vcl_min(sizeTilesX,sizeImageX-startX);
571 572 573 574 575 576 577 578 579
        unsigned long sizeY = vcl_min(sizeTilesY,sizeImageY-startY);

        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
580 581 582 583 584
        labelImage->SetInput(readerIn->GetOutput());
        labelImage->SetStartX(0);
        labelImage->SetStartY(0);
        labelImage->SetSizeX(sizeX);
        labelImage->SetSizeY(sizeY);
585

586 587
        // Relabel tile according to look-up table
        ChangeLabelImageFilterType::Pointer changeLabel = ChangeLabelImageFilterType::New();
Julien Malik's avatar
STYLE  
Julien Malik committed
588
        changeLabel->SetInput(labelImage->GetOutput());
589

590
        // Fill LUT
OTB Bot's avatar
STYLE  
OTB Bot committed
591
        for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
592
          {
Julien Malik's avatar
STYLE  
Julien Malik committed
593
          if(label!=LUT[label])
594 595 596 597 598 599 600 601 602 603 604
            {
            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());
605
        for (it.GoToBegin(); !it.IsAtEnd(); ++it)
606
          {
607
          sizePerRegion[it.Value()]+=1;
608
          }
609 610 611 612 613

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

        // Remove previous tile (not needed anymore)
614
        readerIn = ITK_NULLPTR; // release the input file
615
        RemoveFile(tileIn);
616
        }
617
      }
618

619 620
    // Clear lut, we do not need it anymore
    LUT.clear();
621

622 623
    unsigned int smallCount = 0;

624 625 626 627 628
      // 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)
629
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
630
        if(sizePerRegion[curLabel]<minRegionSize)
631
          {
632
          newLabels[curLabel]=0;
633
          ++smallCount;
634
          }
635
        else
636
          {
637 638
          newLabels[curLabel]=newLab;
          newLab+=1;
639 640 641
          }
        }

642
      otbAppLogINFO(<<smallCount<<" small regions will be removed");
643

644 645
      // Clear sizePerRegion, we do not need it anymore
      sizePerRegion.clear();
646

OTB Bot's avatar
STYLE  
OTB Bot committed
647
      for(unsigned int column = 0; column < nbTilesX; ++column)
648
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
649
        for(unsigned int row = 0; row < nbTilesY; ++row)
650 651 652 653 654 655 656 657
          {
          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
658
          for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
659 660 661 662 663 664
            {
            if(label != newLabels[label])
              {
              changeLabel->SetChange(label,newLabels[label]);
              }
            }
665

666
          // Write the relabeled tile
667
          std::string tmpfile = WriteTile(changeLabel->GetOutput(),row,column,"FINAL");
668 669
          m_FilesToRemoveAfterExecute.push_back(tmpfile);

670
          // Clean previous tiles (not needed anymore)
671
          readerIn = ITK_NULLPTR; // release the input file
672
          RemoveFile(tileIn);
673
          }
674
        }
675 676 677

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

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

683 684
      m_FilesToRemoveAfterExecute.push_back(vrtfile);

685
      clock_t toc = clock();
686

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

689 690 691
      // Final writing
      m_FinalReader = LabelImageReaderType::New();
      m_FinalReader->SetFileName(vrtfile);
692

693 694 695
      m_ImportGeoInformationFilter = ImportGeoInformationImageFilterType::New();
      m_ImportGeoInformationFilter->SetInput(m_FinalReader->GetOutput());
      m_ImportGeoInformationFilter->SetSource(imageIn);
696

697
      SetParameterOutputImage("out",m_ImportGeoInformationFilter->GetOutput());
698
  }
699

700
  void AfterExecuteAndWriteOutputs() ITK_OVERRIDE
701
  {
702
    // Release input files
703
    m_FinalReader = ITK_NULLPTR;
704

705 706 707 708 709
    if(IsParameterEnabled("cleanup"))
      {
      otbAppLogINFO(<<"Final clean-up ...");

      for(std::vector<std::string>::iterator it = m_FilesToRemoveAfterExecute.begin();
OTB Bot's avatar
STYLE  
OTB Bot committed
710
          it!=m_FilesToRemoveAfterExecute.end(); ++it)
711 712 713 714 715 716 717 718 719 720 721 722
        {
        RemoveFile(*it);
        }

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

    m_FilesToRemoveAfterExecute.clear();
723
    m_TmpDirCleanup = false;
724
  }
725 726 727 728 729
};
}
}

OTB_APPLICATION_EXPORT(otb::Wrapper::LSMSSegmentation)