otbLSMSSegmentation.cxx 27.3 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 23 24 25 26 27 28 29 30 31
/*=========================================================================

 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"
#include "itkConnectedComponentFunctorImageFilter.h"
#include "otbBandMathImageFilter.h"
#include "itkStatisticsImageFilter.h"
#include "itkChangeLabelImageFilter.h"
#include "itkImageRegionConstIterator.h"
#include "itkScalarConnectedComponentImageFilter.h"
#include "otbConcatenateVectorImageFilter.h"

#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 75
  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 otb::BandMathImageFilter<LabelImageType> BandMathImageFilterType;
  typedef itk::StatisticsImageFilter<LabelImageType> StatisticsImageFilterType;
  typedef itk::ChangeLabelImageFilter<LabelImageType,LabelImageType> ChangeLabelImageFilterType;
76
  typedef otb::ImportGeoInformationImageFilter<LabelImageType,ImageType> ImportGeoInformationImageFilterType;
77
  typedef itk::ImageRegionConstIterator<LabelImageType> LabelImageIterator;
78
    
OTB Bot's avatar
STYLE  
OTB Bot committed
79
  typedef otb::ConcatenateVectorImageFilter <ImageType,ImageType,ImageType> ConcatenateType;
80

81 82 83 84
  LSMSSegmentation(): m_FinalReader(),m_ImportGeoInformationFilter(),m_FilesToRemoveAfterExecute(),m_TmpDirCleanup(false){}
  
  virtual ~LSMSSegmentation(){}

OTB Bot's avatar
STYLE  
OTB Bot committed
85
private:
86 87
  LabelImageReaderType::Pointer m_FinalReader;
  ImportGeoInformationImageFilterType::Pointer m_ImportGeoInformationFilter;
88 89
  std::vector<string> m_FilesToRemoveAfterExecute;
  bool m_TmpDirCleanup;
90 91

  std::string CreateFileName(unsigned int row, unsigned int column, std::string label)
92
  {
93
    std::string outfname = GetParameterString("out");
94 95
    std::string tilesname = itksys::SystemTools::GetFilenameWithoutExtension(outfname.c_str());

96 97 98 99 100 101
    std::stringstream tileOut;
    tileOut<<tilesname<<"_"<<row<<"_"<<column<<"_"<<label<<".tif";
    
    std::vector<std::string> joins;
    if(IsParameterEnabled("tmpdir"))
      {
102 103 104 105 106 107 108
      std::string tmpdir = GetParameterString("tmpdir");
      
      if(tmpdir.size() > 1 && tmpdir[tmpdir.size()-1] != '/')
        {
        tmpdir.append("/");
        }
      joins.push_back(tmpdir);
109 110 111
      }
    joins.push_back(tileOut.str());
    
Julien Malik's avatar
STYLE  
Julien Malik committed
112
    std::string currentFile = itksys::SystemTools::JoinPath(joins);
113 114 115

    return currentFile;
  }
116
  
117
  std::string WriteTile(LabelImageType::Pointer img, unsigned int row, unsigned int column, std::string label)
OTB Bot's avatar
STYLE  
OTB Bot committed
118
  {
119 120 121 122 123
    std::string currentFile = CreateFileName(row,column,label);
    
    LabelImageWriterType::Pointer imageWriter = LabelImageWriterType::New();
    imageWriter->SetInput(img);
    imageWriter->SetFileName(currentFile);
OTB Bot's avatar
STYLE  
OTB Bot committed
124
    imageWriter->Update();
125 126 127

    return currentFile;
  }
128

129
  void RemoveFile(std::string tile)
130 131 132 133 134 135
  {
    // 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");
136

137 138
      if(itksys::SystemTools::FileExists(geomfile.c_str()))
        {
139 140 141 142 143
        bool res = itksys::SystemTools::RemoveFile(geomfile.c_str());
        if (!res)
          {
          otbAppLogINFO(<<"Unable to remove file  "<<geomfile);
          }
144 145 146
        }
      if(itksys::SystemTools::FileExists(tile.c_str()))
        {
147 148 149 150 151
        bool res = itksys::SystemTools::RemoveFile(tile.c_str());
        if (!res)
          {
          otbAppLogINFO(<<"Unable to remove file  "<<tile);
          }
152 153 154 155 156 157 158 159 160
        }
      }
  }

  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");
161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178
    
    std::stringstream vrtfOut;
    vrtfOut<<itksys::SystemTools::GetFilenameWithoutExtension(outfname.c_str())<<".vrt";
    
    std::vector<std::string> joins;
    if(IsParameterEnabled("tmpdir"))
      {
      std::string tmpdir = GetParameterString("tmpdir");
      
      if(tmpdir.size() > 1 && tmpdir[tmpdir.size()-1] != '/')
        {
        tmpdir.append("/");
        }
      joins.push_back(tmpdir);
      }
    joins.push_back(vrtfOut.str());
    
    std::string vrtfname = itksys::SystemTools::JoinPath(joins);
179
    otbAppLogINFO(<<"Creating temporary vrt file: "<<vrtfname);
180 181 182 183 184 185 186

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

  void DoInit()
  {
    SetName("LSMSSegmentation");
210
    SetDescription("Second step of the exact Large-Scale Mean-Shift segmentation workflow.");
211

212
    SetDocName("Exact Large-Scale Mean-Shift segmentation, step 2");
213
    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 to remove 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 to define 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.");
214 215
    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");
216
    SetDocSeeAlso("MeanShiftSmoothing, LSMSSmallRegionsMerging, LSMSVectorization");
217 218 219 220 221 222 223 224 225 226
    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");
227
    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." );
228 229 230 231 232 233 234 235 236 237 238 239 240 241

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

    AddParameter(ParameterType_Float, "spatialr", "Spatial radius");
    SetParameterDescription("spatialr", "Spatial radius of the neighborhood.");
    SetDefaultParameterFloat("spatialr", 5);
    SetMinimumParameterFloatValue("spatialr", 0);
    MandatoryOff("spatialr");

    AddParameter(ParameterType_Int, "minsize", "Minimum Region Size");
OTB Bot's avatar
STYLE  
OTB Bot committed
242
    SetParameterDescription("minsize", "Minimum Region Size. If, after the segmentation, a region is of size lower than this criterion, the region is deleted.");
243 244 245 246
    SetDefaultParameterInt("minsize", 0);
    SetMinimumParameterIntValue("minsize", 0);
    MandatoryOff("minsize");
    
247 248 249 250 251 252 253 254 255
    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);
256 257

    AddParameter(ParameterType_Directory,"tmpdir","Directory where to write temporary files");
Julien Malik's avatar
Julien Malik committed
258
    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.");
259 260 261
    MandatoryOff("tmpdir");
    DisableParameter("tmpdir");

262 263
    AddParameter(ParameterType_Empty,"cleanup","Temporary files cleaning");
    EnableParameter("cleanup");
264 265 266 267 268 269 270 271 272 273
    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("ranger","15");
    SetDocExampleParameterValue("spatialr","5");
    SetDocExampleParameterValue("minsize","0");
274 275
    SetDocExampleParameterValue("tilesizex","256");
    SetDocExampleParameterValue("tilesizey","256");
276 277 278 279
    
  }

  void DoUpdateParameters()
OTB Bot's avatar
STYLE  
OTB Bot committed
280
  {
281 282 283 284
  }

  void DoExecute()
  {
285 286
    m_FilesToRemoveAfterExecute.clear();

287 288
    clock_t tic = clock();
  
289 290
    const float ranger         = GetParameterFloat("ranger");
    const float spatialr       = GetParameterFloat("spatialr");
291
    
292
    unsigned int minRegionSize = GetParameterInt("minsize");
293
    
294 295
    unsigned long sizeTilesX   = GetParameterInt("tilesizex");
    unsigned long sizeTilesY   = GetParameterInt("tilesizey");
296 297 298 299 300


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

    //Acquisition of the input image dimensions
    ImageType::Pointer imageIn = GetParameterImage("in");
    imageIn->UpdateOutputInformation();
324 325 326 327

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

329 330 331 332 333
    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);

334 335
    unsigned long regionCount = 0;
  
336 337 338
    //Segmentation by the connected component per tile and label
    //shifting
    otbAppLogINFO(<<"Tile shifting ...");
339
      
OTB Bot's avatar
STYLE  
OTB Bot committed
340 341 342
    for(unsigned int row = 0; row < nbTilesY; ++row)
      for(unsigned int column = 0; column < nbTilesX; ++column)
        {
343 344 345 346 347
        // 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);
OTB Bot's avatar
STYLE  
OTB Bot committed
348
           
349 350 351 352 353 354 355 356 357 358
        //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();

359
        //- the spatial image (final positions) if available
360 361
        MultiChannelExtractROIFilterType::Pointer extractROIFilter2 = MultiChannelExtractROIFilterType::New();
        ConcatenateType::Pointer concat = ConcatenateType::New();
362 363 364 365 366 367 368 369 370 371
        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();
OTB Bot's avatar
STYLE  
OTB Bot committed
372
         
373 374 375 376 377 378 379 380 381
          //Concatenation of the two input images
          concat->SetInput1(extractROIFilter->GetOutput());
          concat->SetInput2(extractROIFilter2->GetOutput());
          concat->Update();
          
          ccFilter->SetInput(concat->GetOutput());
          }
        else
          {
OTB Bot's avatar
STYLE  
OTB Bot committed
382
         ccFilter->SetInput(extractROIFilter->GetOutput());
383 384
          }
        
385 386 387
        //Expression 1 : radiometric distance < ranger
        std::stringstream expr;
        expr<<"sqrt((p1b1-p2b1)*(p1b1-p2b1)";
OTB Bot's avatar
STYLE  
OTB Bot committed
388
        for(unsigned int i=1; i<nbComp; i++)
389
          expr<<"+(p1b"<<i+1<<"-p2b"<<i+1<<")*(p1b"<<i+1<<"-p2b"<<i+1<<")";
390
        expr<<")"<<"<"<<ranger;
391

392 393 394 395
        if(HasValue("inpos"))
          {
          //Expression 2 : final positions < spatialr
          expr<<" and sqrt((p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")*(p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")+";
OTB Bot's avatar
STYLE  
OTB Bot committed
396
          expr<<"(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<")*(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<"))"<<"<"<<spatialr;
397 398
          }
        
399 400
        //Segmentation
        ccFilter->GetFunctor().SetExpression(expr.str());
OTB Bot's avatar
STYLE  
OTB Bot committed
401
        ccFilter->Update();
402
    
403 404
        std::stringstream ssexpr;
        ssexpr<<"label+"<<regionCount;
OTB Bot's avatar
STYLE  
OTB Bot committed
405
           
406 407 408 409 410 411 412 413
        //Shifting
        BandMathImageFilterType::Pointer labelBandMath = BandMathImageFilterType::New();
        labelBandMath->SetNthInput(0,ccFilter->GetOutput(),"label");
        labelBandMath->SetExpression(ssexpr.str());
        labelBandMath->Update();

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

418
        std::string filename = WriteTile(labelBandMath->GetOutput(),row,column,"SEG");
419 420
        }
    
421

422 423
    // Step 2: create the look-up table for all overlaps
    otbAppLogINFO(<<"LUT creation ...");
OTB Bot's avatar
STYLE  
OTB Bot committed
424
    std::vector<LabelImagePixelType> LUT;
425 426 427 428
    LUT.clear();
    LUT.resize(regionCount+1);
    for(LabelImagePixelType curLabel = 1; curLabel <= regionCount; ++curLabel)
      LUT[curLabel] = curLabel;
429
      
OTB Bot's avatar
STYLE  
OTB Bot committed
430
    for(unsigned int row = 0; row < nbTilesY; row++)
431
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
432
      for(unsigned int column = 0; column < nbTilesX; column++)
433 434 435 436 437
        {
        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);
438
      
439
        std::string tileIn = CreateFileName(row,column,"SEG");
440

441 442 443 444
        // Read current tile
        LabelImageReaderType::Pointer tileInReader = LabelImageReaderType::New();
        tileInReader->SetFileName(tileIn);
        tileInReader->Update();
445

446 447 448 449
        // 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
450 451
          LabelImageReaderType::Pointer tileUpReader = LabelImageReaderType::New();
          tileUpReader->SetFileName(tileUp);
452 453 454 455 456 457 458 459
          tileUpReader->Update();

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

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

OTB Bot's avatar
STYLE  
OTB Bot committed
460
          for(pixelIndexIn[0]=0; pixelIndexIn[0]<static_cast<long>(sizeX-1); ++pixelIndexIn[0])
461 462 463 464
            {
            pixelIndexUp[0] = pixelIndexIn[0];

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
465
           while(LUT[curCanLabel] != curCanLabel)
466 467 468
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
469
           LabelImagePixelType adjCanLabel = tileUpReader->GetOutput()->GetPixel(pixelIndexUp);
470
            
OTB Bot's avatar
STYLE  
OTB Bot committed
471
           while(LUT[adjCanLabel] != adjCanLabel)
472 473 474
              {
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
475
           if(curCanLabel < adjCanLabel)
476 477 478
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
479
           else
480 481 482 483 484
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
485

486 487 488 489
        // 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
490 491
          LabelImageReaderType::Pointer tileLeftReader = LabelImageReaderType::New();
          tileLeftReader->SetFileName(tileLeft);
492 493 494 495 496 497 498 499
          tileLeftReader->Update();

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

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

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

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
505
           while(LUT[curCanLabel] != curCanLabel)
506 507 508
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
509 510 511
           LabelImagePixelType adjCanLabel = tileLeftReader->GetOutput()->GetPixel(pixelIndexUp);
           while(LUT[adjCanLabel] != adjCanLabel)
              {
512 513
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
514
           if(curCanLabel < adjCanLabel)
515 516 517
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
518
           else
519 520 521 522 523
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
524 525
        }
      }
526
    
527
    // Reduce LUT to canonical labels
528
    for(LabelImagePixelType label = 1; label < regionCount+1; ++label)
OTB Bot's avatar
STYLE  
OTB Bot committed
529
      {
530
      LabelImagePixelType can = label;
OTB Bot's avatar
STYLE  
OTB Bot committed
531
      while(LUT[can] != can)
532 533 534
        {
        can = LUT[can];
        }
535 536
      LUT[label] = can;
      }
537
    otbAppLogINFO(<<"LUT size: "<<LUT.size()<<" segments");
538
 
539 540 541
    // These variables will be used to estimate the size of each
    // region on the flow
    std::vector<unsigned long> sizePerRegion(regionCount+1,0);
542

543
    otbAppLogINFO(<<"Tiles relabelisation ...");
OTB Bot's avatar
STYLE  
OTB Bot committed
544
    for(unsigned int column = 0; column < nbTilesX; ++column)
545
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
546
      for(unsigned int row = 0; row < nbTilesY; ++row)
547
        {
Julien Malik's avatar
STYLE  
Julien Malik committed
548
        unsigned long startX = column*sizeTilesX;
549
        unsigned long startY = row*sizeTilesY;
Julien Malik's avatar
STYLE  
Julien Malik committed
550
        unsigned long sizeX = vcl_min(sizeTilesX,sizeImageX-startX);
551 552 553 554 555 556 557 558 559
        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
560 561 562 563 564
        labelImage->SetInput(readerIn->GetOutput());
        labelImage->SetStartX(0);
        labelImage->SetStartY(0);
        labelImage->SetSizeX(sizeX);
        labelImage->SetSizeY(sizeY);
565 566 567
        
        // Relabel tile according to look-up table
        ChangeLabelImageFilterType::Pointer changeLabel = ChangeLabelImageFilterType::New();
Julien Malik's avatar
STYLE  
Julien Malik committed
568
        changeLabel->SetInput(labelImage->GetOutput());
OTB Bot's avatar
STYLE  
OTB Bot committed
569
       
570
        // Fill LUT
OTB Bot's avatar
STYLE  
OTB Bot committed
571
        for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
572
          {
Julien Malik's avatar
STYLE  
Julien Malik committed
573
          if(label!=LUT[label])
574 575 576 577 578 579 580 581 582 583 584
            {
            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());
OTB Bot's avatar
STYLE  
OTB Bot committed
585
       for (it = it.Begin(); !it.IsAtEnd(); ++it)
586
          {
587
          sizePerRegion[it.Value()]+=1;
588
          }
589 590 591 592 593

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

        // Remove previous tile (not needed anymore)
594
        readerIn = 0; // release the input file
595
        RemoveFile(tileIn);
596
        }
597
      }
598

599 600
    // Clear lut, we do not need it anymore
    LUT.clear();
601

602 603
    unsigned int smallCount = 0;

604 605 606 607 608
      // 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)
609
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
610
        if(sizePerRegion[curLabel]<minRegionSize)
611
          {
612
          newLabels[curLabel]=0;
613
          ++smallCount;
614
          }
615
        else
616
          {
617 618
          newLabels[curLabel]=newLab;
          newLab+=1;
619 620 621
          }
        }

622
      otbAppLogINFO(<<smallCount<<" small regions will be removed");
623

624 625 626
      // Clear sizePerRegion, we do not need it anymore
      sizePerRegion.clear();
     
OTB Bot's avatar
STYLE  
OTB Bot committed
627
      for(unsigned int column = 0; column < nbTilesX; ++column)
628
        {
OTB Bot's avatar
STYLE  
OTB Bot committed
629
        for(unsigned int row = 0; row < nbTilesY; ++row)
630 631 632 633 634 635 636 637
          {
          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
638
          for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
639 640 641 642 643 644 645 646
            {
            if(label != newLabels[label])
              {
              changeLabel->SetChange(label,newLabels[label]);
              }
            }
          
          // Write the relabeled tile
647
          std::string tmpfile = WriteTile(changeLabel->GetOutput(),row,column,"FINAL");
648 649
          m_FilesToRemoveAfterExecute.push_back(tmpfile);

650
          // Clean previous tiles (not needed anymore)
651
          readerIn = 0; // release the input file
652
          RemoveFile(tileIn);
653
          }
654
        }
655 656 657

      // Clear newLabels, we do not need it anymore
      newLabels.clear();
658
  
659 660
      // Here we write a temporary vrt file that will be used to
      // stitch together all the tiles
Julien Michel's avatar
Julien Michel committed
661
      std::string vrtfile = WriteVRTFile(nbTilesX,nbTilesY,sizeTilesX,sizeTilesY,sizeImageX,sizeImageY);
662

663 664
      m_FilesToRemoveAfterExecute.push_back(vrtfile);

665 666 667
      clock_t toc = clock();
      
      otbAppLogINFO(<<"Elapsed time: "<<(double)(toc - tic) / CLOCKS_PER_SEC<<" seconds");
668

669 670 671
      // Final writing
      m_FinalReader = LabelImageReaderType::New();
      m_FinalReader->SetFileName(vrtfile);
672

673 674 675
      m_ImportGeoInformationFilter = ImportGeoInformationImageFilterType::New();
      m_ImportGeoInformationFilter->SetInput(m_FinalReader->GetOutput());
      m_ImportGeoInformationFilter->SetSource(imageIn);
676

677
      SetParameterOutputImage("out",m_ImportGeoInformationFilter->GetOutput());
678
  }
679 680 681

  void AfterExecuteAndWriteOutputs()
  {
682 683 684
    // Release input files
    m_FinalReader = 0;
    
685 686 687 688 689
    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
690
          it!=m_FilesToRemoveAfterExecute.end(); ++it)
691 692 693 694 695 696 697 698 699 700 701 702
        {
        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();
703
    m_TmpDirCleanup = false;
704
  }
705 706 707 708 709 710 711
};
}
}

OTB_APPLICATION_EXPORT(otb::Wrapper::LSMSSegmentation)