otbLSMSSegmentation.cxx 27.2 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
/*=========================================================================

 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 "otbBandMathImageFilter.h"
24
#include "itkUnaryFunctorImageFilter.h"
25 26 27 28 29 30 31
#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
  LSMSSegmentation(): m_FinalReader(),m_ImportGeoInformationFilter(),m_FilesToRemoveAfterExecute(),m_TmpDirCleanup(false){}
82

83 84
  virtual ~LSMSSegmentation(){}

OTB Bot's avatar
STYLE  
OTB Bot committed
85
private:
86 87
  LabelImageReaderType::Pointer m_FinalReader;
  ImportGeoInformationImageFilterType::Pointer m_ImportGeoInformationFilter;
Julien Michel's avatar
Julien Michel committed
88
  std::vector<std::string> m_FilesToRemoveAfterExecute;
89
  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
    std::stringstream tileOut;
    tileOut<<tilesname<<"_"<<row<<"_"<<column<<"_"<<label<<".tif";
98

99 100 101
    std::vector<std::string> joins;
    if(IsParameterEnabled("tmpdir"))
      {
102
      std::string tmpdir = GetParameterString("tmpdir");
103

104 105 106 107 108
      if(tmpdir.size() > 1 && tmpdir[tmpdir.size()-1] != '/')
        {
        tmpdir.append("/");
        }
      joins.push_back(tmpdir);
109 110
      }
    joins.push_back(tileOut.str());
111

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
    std::string currentFile = CreateFileName(row,column,label);
120

121 122 123
    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
    std::stringstream vrtfOut;
    vrtfOut<<itksys::SystemTools::GetFilenameWithoutExtension(outfname.c_str())<<".vrt";
164

165 166 167 168
    std::vector<std::string> joins;
    if(IsParameterEnabled("tmpdir"))
      {
      std::string tmpdir = GetParameterString("tmpdir");
169

170 171 172 173 174 175 176
      if(tmpdir.size() > 1 && tmpdir[tmpdir.size()-1] != '/')
        {
        tmpdir.append("/");
        }
      joins.push_back(tmpdir);
      }
    joins.push_back(vrtfOut.str());
177

178
    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
    SetDefaultOutputPixelType("out",ImagePixelType_uint32);
229 230 231 232 233 234 235 236 237 238 239 240 241 242

    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
243
    SetParameterDescription("minsize", "Minimum Region Size. If, after the segmentation, a region is of size lower than this criterion, the region is deleted.");
244 245 246
    SetDefaultParameterInt("minsize", 0);
    SetMinimumParameterIntValue("minsize", 0);
    MandatoryOff("minsize");
247

248 249 250 251 252 253 254 255 256
    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);
257 258

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

263 264
    AddParameter(ParameterType_Empty,"cleanup","Temporary files cleaning");
    EnableParameter("cleanup");
265 266 267 268 269 270 271 272 273 274
    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");
275 276
    SetDocExampleParameterValue("tilesizex","256");
    SetDocExampleParameterValue("tilesizey","256");
277

278 279 280
  }

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

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

288
    clock_t tic = clock();
289

290 291
    const float ranger         = GetParameterFloat("ranger");
    const float spatialr       = GetParameterFloat("spatialr");
292

293
    unsigned int minRegionSize = GetParameterInt("minsize");
294

295 296
    unsigned long sizeTilesX   = GetParameterInt("tilesizex");
    unsigned long sizeTilesY   = GetParameterInt("tilesizey");
297 298 299 300 301


    // Ensure that temporary directory exists if activated:
    if(IsParameterEnabled("tmpdir"))
      {
302 303 304 305
      if(!itksys::SystemTools::FileExists(GetParameterString("tmpdir").c_str()))
        {
        m_TmpDirCleanup = true;
        }
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

317
    if(HasValue("inpos"))
318
      {
319
      spatialIn = GetParameterImage("inpos");
320 321 322 323 324
      }

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

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

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

335
    unsigned long regionCount = 0;
336

337 338 339
    //Segmentation by the connected component per tile and label
    //shifting
    otbAppLogINFO(<<"Tile shifting ...");
340

OTB Bot's avatar
STYLE  
OTB Bot committed
341 342 343
    for(unsigned int row = 0; row < nbTilesY; ++row)
      for(unsigned int column = 0; column < nbTilesX; ++column)
        {
344 345 346 347 348
        // 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);
349

350 351 352 353 354 355 356 357 358 359
        //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();

360
        //- the spatial image (final positions) if available
361 362
        MultiChannelExtractROIFilterType::Pointer extractROIFilter2 = MultiChannelExtractROIFilterType::New();
        ConcatenateType::Pointer concat = ConcatenateType::New();
363 364 365 366 367 368 369 370 371 372
        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();
373

374 375 376 377
          //Concatenation of the two input images
          concat->SetInput1(extractROIFilter->GetOutput());
          concat->SetInput2(extractROIFilter2->GetOutput());
          concat->Update();
378

379 380 381 382
          ccFilter->SetInput(concat->GetOutput());
          }
        else
          {
OTB Bot's avatar
STYLE  
OTB Bot committed
383
         ccFilter->SetInput(extractROIFilter->GetOutput());
384
          }
385

386 387 388
        //Expression 1 : radiometric distance < ranger
        std::stringstream expr;
        expr<<"sqrt((p1b1-p2b1)*(p1b1-p2b1)";
OTB Bot's avatar
STYLE  
OTB Bot committed
389
        for(unsigned int i=1; i<nbComp; i++)
390
          expr<<"+(p1b"<<i+1<<"-p2b"<<i+1<<")*(p1b"<<i+1<<"-p2b"<<i+1<<")";
391
        expr<<")"<<"<"<<ranger;
392

393 394 395
        if(HasValue("inpos"))
          {
          //Expression 2 : final positions < spatialr
396
          expr<<" and sqrt((p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")*(p1b"<<nbComp+1<<"-p2b"<<nbComp+1<<")+";
OTB Bot's avatar
STYLE  
OTB Bot committed
397
          expr<<"(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<")*(p1b"<<nbComp+2<<"-p2b"<<nbComp+2<<"))"<<"<"<<spatialr;
398
          }
399

400 401
        //Segmentation
        ccFilter->GetFunctor().SetExpression(expr.str());
OTB Bot's avatar
STYLE  
OTB Bot committed
402
        ccFilter->Update();
403

404 405
        std::stringstream ssexpr;
        ssexpr<<"label+"<<regionCount;
406

407 408 409 410 411 412 413 414
        //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
415
        stats->SetInput(ccFilter->GetOutput());
416 417 418
        stats->Update();
        regionCount+=stats->GetMaximum();

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

422

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

OTB Bot's avatar
STYLE  
OTB Bot committed
431
    for(unsigned int row = 0; row < nbTilesY; row++)
432
      {
OTB Bot's avatar
STYLE  
OTB Bot committed
433
      for(unsigned int column = 0; column < nbTilesX; column++)
434 435 436 437 438
        {
        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);
439

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

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

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

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

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

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

            LabelImagePixelType curCanLabel = tileInReader->GetOutput()->GetPixel(pixelIndexIn);
OTB Bot's avatar
STYLE  
OTB Bot committed
466
           while(LUT[curCanLabel] != curCanLabel)
467 468 469
              {
              curCanLabel = LUT[curCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
470
           LabelImagePixelType adjCanLabel = tileUpReader->GetOutput()->GetPixel(pixelIndexUp);
471

OTB Bot's avatar
STYLE  
OTB Bot committed
472
           while(LUT[adjCanLabel] != adjCanLabel)
473 474 475
              {
              adjCanLabel = LUT[adjCanLabel];
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
476
           if(curCanLabel < adjCanLabel)
477 478 479
              {
              LUT[adjCanLabel] = curCanLabel;
              }
OTB Bot's avatar
STYLE  
OTB Bot committed
480
           else
481 482 483 484 485
              {
              LUT[LUT[curCanLabel]] = adjCanLabel; LUT[curCanLabel] = adjCanLabel;
              }
            }
          }
486

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

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

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

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

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

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

540 541 542
    // These variables will be used to estimate the size of each
    // region on the flow
    std::vector<unsigned long> sizePerRegion(regionCount+1,0);
543

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

567 568
        // Relabel tile according to look-up table
        ChangeLabelImageFilterType::Pointer changeLabel = ChangeLabelImageFilterType::New();
Julien Malik's avatar
STYLE  
Julien Malik committed
569
        changeLabel->SetInput(labelImage->GetOutput());
570

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

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

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

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

603 604
    unsigned int smallCount = 0;

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

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

625 626
      // Clear sizePerRegion, we do not need it anymore
      sizePerRegion.clear();
627

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

647
          // Write the relabeled tile
648
          std::string tmpfile = WriteTile(changeLabel->GetOutput(),row,column,"FINAL");
649 650
          m_FilesToRemoveAfterExecute.push_back(tmpfile);

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

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

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

664 665
      m_FilesToRemoveAfterExecute.push_back(vrtfile);

666
      clock_t toc = clock();
667

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

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

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

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

  void AfterExecuteAndWriteOutputs()
  {
683 684
    // Release input files
    m_FinalReader = 0;
685

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

OTB_APPLICATION_EXPORT(otb::Wrapper::LSMSSegmentation)