otbLSMSSegmentation.cxx 27.5 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
/*=========================================================================

 Program:   ORFEO Toolbox
 Language:  C++
 Date:      $Date$
 Version:   $Revision$


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


 This software is distributed WITHOUT ANY WARRANTY; without even
 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 PURPOSE.  See the above copyright notices for more information.

 =========================================================================*/
#include "otbImageFileReader.h"
#include "otbImageFileWriter.h"
#include "otbMultiChannelExtractROI.h"
#include "otbExtractROI.h"
#include "otbConnectedComponentMuParserFunctor.h"
23
#include "itkUnaryFunctorImageFilter.h"
24 25 26 27 28
#include "itkStatisticsImageFilter.h"
#include "itkChangeLabelImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkScalarConnectedComponentImageFilter.h"
#include "otbConcatenateVectorImageFilter.h"
29
#include "otbAffineFunctor.h"
30 31

#include "otbMultiToMonoChannelExtractROI.h"
32
#include "otbImportGeoInformationImageFilter.h"
33 34 35 36 37 38 39 40 41 42

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

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

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

43

44 45 46 47 48 49 50 51 52 53 54 55 56
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);

57
  itkTypeMacro(LSMSSegmentation, otb::Application);
58

59 60 61
  typedef FloatVectorImageType              ImageType;
  typedef ImageType::InternalPixelType      ImagePixelType;
  typedef UInt32ImageType                   LabelImageType;
62
  typedef LabelImageType::InternalPixelType LabelImagePixelType;
63 64 65
  typedef otb::ImageFileReader<ImageType> ImageReaderType;
  typedef otb::ImageFileWriter<ImageType> ImageWriterType;
  typedef otb::ImageFileReader<LabelImageType> LabelImageReaderType;
OTB Bot's avatar
STYLE  
OTB Bot committed
66
  typedef otb::ImageFileWriter<LabelImageType> LabelImageWriterType;
67
  typedef otb::MultiChannelExtractROI <ImagePixelType,ImagePixelType > MultiChannelExtractROIFilterType;
OTB Bot's avatar
STYLE  
OTB Bot committed
68
  typedef otb::ExtractROI<LabelImagePixelType,LabelImagePixelType> ExtractROIFilterType;
69 70 71 72 73 74
  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;
75
  typedef otb::ImportGeoInformationImageFilter<LabelImageType,ImageType> ImportGeoInformationImageFilterType;
76
  typedef itk::ImageRegionConstIterator<LabelImageType> LabelImageIterator;
77

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

88
  LSMSSegmentation(): m_FinalReader(),m_ImportGeoInformationFilter(),m_FilesToRemoveAfterExecute(),m_TmpDirCleanup(false){}
89

90
  ~LSMSSegmentation() ITK_OVERRIDE{}
91

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

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

103 104
    std::stringstream tileOut;
    tileOut<<tilesname<<"_"<<row<<"_"<<column<<"_"<<label<<".tif";
105

106 107 108
    std::vector<std::string> joins;
    if(IsParameterEnabled("tmpdir"))
      {
109
      std::string tmpdir = GetParameterString("tmpdir");
110

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

Julien Malik's avatar
STYLE  
Julien Malik committed
119
    std::string currentFile = itksys::SystemTools::JoinPath(joins);
120 121 122

    return currentFile;
  }
123

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

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

    return currentFile;
  }
135

136
  void RemoveFile(std::string tile)
137 138 139 140 141 142
  {
    // 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");
143

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

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

169 170
    std::stringstream vrtfOut;
    vrtfOut<<itksys::SystemTools::GetFilenameWithoutExtension(outfname.c_str())<<".vrt";
171

172 173 174 175
    std::vector<std::string> joins;
    if(IsParameterEnabled("tmpdir"))
      {
      std::string tmpdir = GetParameterString("tmpdir");
176

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

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

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

214
  void DoInit() ITK_OVERRIDE
215 216
  {
    SetName("LSMSSegmentation");
217
    SetDescription("Second step of the exact Large-Scale Mean-Shift segmentation workflow.");
218

219
    SetDocName("Exact Large-Scale Mean-Shift segmentation, step 2");
220
    SetDocLongDescription("This application performs the second step of the exact Large-Scale Mean-Shift segmentation workflow (LSMS). Filtered range image and spatial image should be created with the MeanShiftSmoothing application, 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. 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 nbtilesx and nbtilesy parameters for tile-wise processing, with the guarantees of identical results. 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. Please also note that the output image type should be set to uint32 to ensure that there are enough labels available.");
221 222
    SetDocLimitations("This application is part of the Large-Scale Mean-Shift segmentation workflow (LSMS) and may not be suited for any other purpose.");
    SetDocAuthors("David Youssefi");
223
    SetDocSeeAlso("MeanShiftSmoothing, LSMSSmallRegionsMerging, LSMSVectorization");
224 225 226 227 228 229 230 231 232 233
    AddDocTag(Tags::Segmentation);
    AddDocTag("LSMS");

    AddParameter(ParameterType_InputImage,  "in",    "Filtered image");
    SetParameterDescription( "in", "The filtered image (cf. Adaptive MeanShift Smoothing application)." );
    AddParameter(ParameterType_InputImage,  "inpos",    "Spatial image");
    SetParameterDescription( "inpos", " The spatial image. Spatial input is the displacement map (output of the Adaptive MeanShift Smoothing application)." );
    MandatoryOff("inpos");

    AddParameter(ParameterType_OutputImage, "out", "Output Image");
234
    SetParameterDescription( "out", "The output image. The output image is the segmentation of the filtered image. It is recommended to set the pixel type to uint32." );
235
    SetDefaultOutputPixelType("out",ImagePixelType_uint32);
236 237 238 239 240 241

    AddParameter(ParameterType_Float, "spatialr", "Spatial radius");
    SetParameterDescription("spatialr", "Spatial radius of the neighborhood.");
    SetDefaultParameterFloat("spatialr", 5);
    SetMinimumParameterFloatValue("spatialr", 0);
    MandatoryOff("spatialr");
242 243 244 245 246 247
    
    AddParameter(ParameterType_Float, "ranger", "Range radius");
    SetParameterDescription("ranger", "Range radius defining the radius (expressed in radiometry unit) in the multi-spectral space.");
    SetDefaultParameterFloat("ranger", 15);
    SetMinimumParameterFloatValue("ranger", 0);
    MandatoryOff("ranger");
248 249

    AddParameter(ParameterType_Int, "minsize", "Minimum Region Size");
OTB Bot's avatar
STYLE  
OTB Bot committed
250
    SetParameterDescription("minsize", "Minimum Region Size. If, after the segmentation, a region is of size lower than this criterion, the region is deleted.");
251 252 253
    SetDefaultParameterInt("minsize", 0);
    SetMinimumParameterIntValue("minsize", 0);
    MandatoryOff("minsize");
254

255 256 257 258 259 260 261 262 263
    AddParameter(ParameterType_Int, "tilesizex", "Size of tiles in pixel (X-axis)");
    SetParameterDescription("tilesizex", "Size of tiles along the X-axis.");
    SetDefaultParameterInt("tilesizex", 500);
    SetMinimumParameterIntValue("tilesizex", 1);

    AddParameter(ParameterType_Int, "tilesizey", "Size of tiles in pixel (Y-axis)");
    SetParameterDescription("tilesizey", "Size of tiles along the Y-axis.");
    SetDefaultParameterInt("tilesizey", 500);
    SetMinimumParameterIntValue("tilesizey", 1);
264 265

    AddParameter(ParameterType_Directory,"tmpdir","Directory where to write temporary files");
Julien Malik's avatar
Julien Malik committed
266
    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.");
267 268 269
    MandatoryOff("tmpdir");
    DisableParameter("tmpdir");

270 271
    AddParameter(ParameterType_Empty,"cleanup","Temporary files cleaning");
    EnableParameter("cleanup");
272 273 274 275 276 277 278 279
    SetParameterDescription("cleanup","If activated, the application will try to clean all temporary files it created");
    MandatoryOff("cleanup");

    // Doc example parameter settings
    SetDocExampleParameterValue("in","smooth.tif");
    SetDocExampleParameterValue("inpos","position.tif");
    SetDocExampleParameterValue("out","segmentation.tif");
    SetDocExampleParameterValue("spatialr","5");
280
    SetDocExampleParameterValue("ranger","15");
281
    SetDocExampleParameterValue("minsize","0");
282 283
    SetDocExampleParameterValue("tilesizex","256");
    SetDocExampleParameterValue("tilesizey","256");
284

285 286
  }

287
  void DoUpdateParameters() ITK_OVERRIDE
OTB Bot's avatar
STYLE  
OTB Bot committed
288
  {
289 290
  }

291
  void DoExecute() ITK_OVERRIDE
292
  {
293 294
    m_FilesToRemoveAfterExecute.clear();

295
    clock_t tic = clock();
296

297 298
    const float ranger         = GetParameterFloat("ranger");
    const float spatialr       = GetParameterFloat("spatialr");
299

300
    unsigned int minRegionSize = GetParameterInt("minsize");
301

302 303
    unsigned long sizeTilesX   = GetParameterInt("tilesizex");
    unsigned long sizeTilesY   = GetParameterInt("tilesizey");
304 305 306 307 308


    // Ensure that temporary directory exists if activated:
    if(IsParameterEnabled("tmpdir"))
      {
309 310 311 312
      if(!itksys::SystemTools::FileExists(GetParameterString("tmpdir").c_str()))
        {
        m_TmpDirCleanup = true;
        }
313 314 315 316 317 318 319 320 321 322
      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;
323

324
    if(HasValue("inpos"))
325
      {
326
      spatialIn = GetParameterImage("inpos");
327 328 329 330 331
      }

    //Acquisition of the input image dimensions
    ImageType::Pointer imageIn = GetParameterImage("in");
    imageIn->UpdateOutputInformation();
332 333 334 335

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

337 338 339 340 341
    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);

342
    unsigned long regionCount = 0;
343

344 345 346
    //Segmentation by the connected component per tile and label
    //shifting
    otbAppLogINFO(<<"Tile shifting ...");
347

OTB Bot's avatar
STYLE  
OTB Bot committed
348 349 350
    for(unsigned int row = 0; row < nbTilesY; ++row)
      for(unsigned int column = 0; column < nbTilesX; ++column)
        {
351 352 353 354 355
        // 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);
356

357 358 359 360 361 362 363 364 365 366
        //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();

367
        //- the spatial image (final positions) if available
368 369
        MultiChannelExtractROIFilterType::Pointer extractROIFilter2 = MultiChannelExtractROIFilterType::New();
        ConcatenateType::Pointer concat = ConcatenateType::New();
370 371 372 373 374 375 376 377 378 379
        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();
380

381 382 383 384
          //Concatenation of the two input images
          concat->SetInput1(extractROIFilter->GetOutput());
          concat->SetInput2(extractROIFilter2->GetOutput());
          concat->Update();
385

386 387 388 389
          ccFilter->SetInput(concat->GetOutput());
          }
        else
          {
OTB Bot's avatar
STYLE  
OTB Bot committed
390
         ccFilter->SetInput(extractROIFilter->GetOutput());
391
          }
392

393 394 395
        //Expression 1 : radiometric distance < ranger
        std::stringstream expr;
        expr<<"sqrt((p1b1-p2b1)*(p1b1-p2b1)";
OTB Bot's avatar
STYLE  
OTB Bot committed
396
        for(unsigned int i=1; i<nbComp; i++)
397
          expr<<"+(p1b"<<i+1<<"-p2b"<<i+1<<")*(p1b"<<i+1<<"-p2b"<<i+1<<")";
398
        expr<<")"<<"<"<<ranger;
399

400 401 402
        if(HasValue("inpos"))
          {
          //Expression 2 : final positions < spatialr
403
          expr<<" and sqrt((p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")*(p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")+";
OTB Bot's avatar
STYLE  
OTB Bot committed
404
          expr<<"(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<")*(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<"))"<<"<"<<spatialr;
405
          }
406

407 408
        //Segmentation
        ccFilter->GetFunctor().SetExpression(expr.str());
OTB Bot's avatar
STYLE  
OTB Bot committed
409
        ccFilter->Update();
410

411
        //Shifting
412 413 414 415 416
        LabelShiftFilterType::Pointer labelShiftFilter = LabelShiftFilterType::New();
        labelShiftFilter->SetInput(ccFilter->GetOutput());
        labelShiftFilter->GetFunctor().SetA(1);
        labelShiftFilter->GetFunctor().SetB(regionCount);
        labelShiftFilter->Update();
417 418 419

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

424
        std::string filename = WriteTile(labelShiftFilter->GetOutput(),row,column,"SEG");
425
        }
426

427

428 429
    // Step 2: create the look-up table for all overlaps
    otbAppLogINFO(<<"LUT creation ...");
OTB Bot's avatar
STYLE  
OTB Bot committed
430
    std::vector<LabelImagePixelType> LUT;
431 432 433 434
    LUT.clear();
    LUT.resize(regionCount+1);
    for(LabelImagePixelType curLabel = 1; curLabel <= regionCount; ++curLabel)
      LUT[curLabel] = curLabel;
435

OTB Bot's avatar
STYLE  
OTB Bot committed
436
    for(unsigned int row = 0; row < nbTilesY; row++)
437
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
438
      for(unsigned int column = 0; column < nbTilesX; column++)
439 440 441 442 443
        {
        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);
444

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

447 448 449 450
        // Read current tile
        LabelImageReaderType::Pointer tileInReader = LabelImageReaderType::New();
        tileInReader->SetFileName(tileIn);
        tileInReader->Update();
451

452 453 454 455
        // 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
456 457
          LabelImageReaderType::Pointer tileUpReader = LabelImageReaderType::New();
          tileUpReader->SetFileName(tileUp);
458 459 460 461 462 463 464 465
          tileUpReader->Update();

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

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

OTB Bot's avatar
STYLE  
OTB Bot committed
466
          for(pixelIndexIn[0]=0; pixelIndexIn[0]<static_cast<long>(sizeX-1); ++pixelIndexIn[0])
467 468 469 470
            {
            pixelIndexUp[0] = pixelIndexIn[0];

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
471
           while(LUT[curCanLabel] != curCanLabel)
472 473 474
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
475
           LabelImagePixelType adjCanLabel = tileUpReader->GetOutput()->GetPixel(pixelIndexUp);
476

OTB Bot's avatar
STYLE  
OTB Bot committed
477
           while(LUT[adjCanLabel] != adjCanLabel)
478 479 480
              {
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
481
           if(curCanLabel < adjCanLabel)
482 483 484
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
485
           else
486 487 488 489 490
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
491

492 493 494 495
        // 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
496 497
          LabelImageReaderType::Pointer tileLeftReader = LabelImageReaderType::New();
          tileLeftReader->SetFileName(tileLeft);
498 499 500 501 502 503 504 505
          tileLeftReader->Update();

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

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

OTB Bot's avatar
STYLE  
OTB Bot committed
506
          for(pixelIndexIn[1]=0; pixelIndexIn[1]<static_cast<long>(sizeY-1); ++pixelIndexIn[1])
507 508 509 510
            {
            pixelIndexUp[1] = pixelIndexIn[1];

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
511
           while(LUT[curCanLabel] != curCanLabel)
512 513 514
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
515 516 517
           LabelImagePixelType adjCanLabel = tileLeftReader->GetOutput()->GetPixel(pixelIndexUp);
           while(LUT[adjCanLabel] != adjCanLabel)
              {
518 519
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
520
           if(curCanLabel < adjCanLabel)
521 522 523
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
524
           else
525 526 527 528 529
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
530 531
        }
      }
532

533
    // Reduce LUT to canonical labels
534
    for(LabelImagePixelType label = 1; label < regionCount+1; ++label)
OTB Bot's avatar
STYLE  
OTB Bot committed
535
      {
536
      LabelImagePixelType can = label;
OTB Bot's avatar
STYLE  
OTB Bot committed
537
      while(LUT[can] != can)
538 539 540
        {
        can = LUT[can];
        }
541 542
      LUT[label] = can;
      }
543
    otbAppLogINFO(<<"LUT size: "<<LUT.size()<<" segments");
544

545 546 547
    // These variables will be used to estimate the size of each
    // region on the flow
    std::vector<unsigned long> sizePerRegion(regionCount+1,0);
548

549
    otbAppLogINFO(<<"Tiles relabelisation ...");
OTB Bot's avatar
STYLE  
OTB Bot committed
550
    for(unsigned int column = 0; column < nbTilesX; ++column)
551
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
552
      for(unsigned int row = 0; row < nbTilesY; ++row)
553
        {
Julien Malik's avatar
STYLE  
Julien Malik committed
554
        unsigned long startX = column*sizeTilesX;
555
        unsigned long startY = row*sizeTilesY;
Julien Malik's avatar
STYLE  
Julien Malik committed
556
        unsigned long sizeX = vcl_min(sizeTilesX,sizeImageX-startX);
557 558 559 560 561 562 563 564 565
        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
566 567 568 569 570
        labelImage->SetInput(readerIn->GetOutput());
        labelImage->SetStartX(0);
        labelImage->SetStartY(0);
        labelImage->SetSizeX(sizeX);
        labelImage->SetSizeY(sizeY);
571

572 573
        // Relabel tile according to look-up table
        ChangeLabelImageFilterType::Pointer changeLabel = ChangeLabelImageFilterType::New();
Julien Malik's avatar
STYLE  
Julien Malik committed
574
        changeLabel->SetInput(labelImage->GetOutput());
575

576
        // Fill LUT
OTB Bot's avatar
STYLE  
OTB Bot committed
577
        for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
578
          {
Julien Malik's avatar
STYLE  
Julien Malik committed
579
          if(label!=LUT[label])
580 581 582 583 584 585 586 587 588 589 590
            {
            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());
591
        for (it.GoToBegin(); !it.IsAtEnd(); ++it)
592
          {
593
          sizePerRegion[it.Value()]+=1;
594
          }
595 596 597 598 599

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

        // Remove previous tile (not needed anymore)
600
        readerIn = ITK_NULLPTR; // release the input file
601
        RemoveFile(tileIn);
602
        }
603
      }
604

605 606
    // Clear lut, we do not need it anymore
    LUT.clear();
607

608 609
    unsigned int smallCount = 0;

610 611 612 613 614
      // 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)
615
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
616
        if(sizePerRegion[curLabel]<minRegionSize)
617
          {
618
          newLabels[curLabel]=0;
619
          ++smallCount;
620
          }
621
        else
622
          {
623 624
          newLabels[curLabel]=newLab;
          newLab+=1;
625 626 627
          }
        }

628
      otbAppLogINFO(<<smallCount<<" small regions will be removed");
629

630 631
      // Clear sizePerRegion, we do not need it anymore
      sizePerRegion.clear();
632

OTB Bot's avatar
STYLE  
OTB Bot committed
633
      for(unsigned int column = 0; column < nbTilesX; ++column)
634
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
635
        for(unsigned int row = 0; row < nbTilesY; ++row)
636 637 638 639 640 641 642 643
          {
          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
644
          for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
645 646 647 648 649 650
            {
            if(label != newLabels[label])
              {
              changeLabel->SetChange(label,newLabels[label]);
              }
            }
651

652
          // Write the relabeled tile
653
          std::string tmpfile = WriteTile(changeLabel->GetOutput(),row,column,"FINAL");
654 655
          m_FilesToRemoveAfterExecute.push_back(tmpfile);

656
          // Clean previous tiles (not needed anymore)
657
          readerIn = ITK_NULLPTR; // release the input file
658
          RemoveFile(tileIn);
659
          }
660
        }
661 662 663

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

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

669 670
      m_FilesToRemoveAfterExecute.push_back(vrtfile);

671
      clock_t toc = clock();
672

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

675 676 677
      // Final writing
      m_FinalReader = LabelImageReaderType::New();
      m_FinalReader->SetFileName(vrtfile);
678

679 680 681
      m_ImportGeoInformationFilter = ImportGeoInformationImageFilterType::New();
      m_ImportGeoInformationFilter->SetInput(m_FinalReader->GetOutput());
      m_ImportGeoInformationFilter->SetSource(imageIn);
682

683
      SetParameterOutputImage("out",m_ImportGeoInformationFilter->GetOutput());
684
  }
685

686
  void AfterExecuteAndWriteOutputs() ITK_OVERRIDE
687
  {
688
    // Release input files
689
    m_FinalReader = ITK_NULLPTR;
690

691 692 693 694 695
    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
696
          it!=m_FilesToRemoveAfterExecute.end(); ++it)
697 698 699 700 701 702 703 704 705 706 707 708
        {
        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();
709
    m_TmpDirCleanup = false;
710
  }
711 712 713 714 715
};
}
}

OTB_APPLICATION_EXPORT(otb::Wrapper::LSMSSegmentation)