An error occurred while loading the file. Please try again.
-
Victor Poughon authoredb8578d99
/*
* Copyright (C) 2005-2019 Centre National d'Etudes Spatiales (CNES)
*
* This file is part of Orfeo Toolbox
*
* https://www.orfeo-toolbox.org/
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
#include "otbMultiChannelExtractROI.h"
#include "otbExtractROI.h"
#include "otbStreamingStatisticsImageFilter.h"
#include "otbSystem.h"
#include "itkUnaryFunctorImageFilter.h"
#include "itkChangeLabelImageFilter.h"
#include "otbTileImageFilter.h"
#include <time.h>
#include <algorithm>
#include <climits>
#include "otbWrapperApplication.h"
#include "otbWrapperApplicationFactory.h"
#include "otbStandardWriterWatcher.h"
namespace otb
{
namespace Wrapper
{
class LSMSSmallRegionsMerging : public Application
{
public:
typedef LSMSSmallRegionsMerging Self;
typedef Application Superclass;
typedef itk::SmartPointer<Self> Pointer;
typedef itk::SmartPointer<const Self> ConstPointer;
typedef FloatVectorImageType ImageType;
typedef ImageType::InternalPixelType ImagePixelType;
typedef UInt32ImageType LabelImageType;
typedef LabelImageType::InternalPixelType LabelImagePixelType;
typedef otb::MultiChannelExtractROI <ImagePixelType,ImagePixelType > MultiChannelExtractROIFilterType;
typedef otb::ExtractROI<LabelImagePixelType,LabelImagePixelType> ExtractROIFilterType;
typedef otb::StreamingStatisticsImageFilter<LabelImageType> StatisticsImageFilterType;
typedef itk::ImageRegionConstIterator<LabelImageType> LabelImageIterator;
typedef itk::ImageRegionConstIterator<ImageType> ImageIterator;
typedef itk::ChangeLabelImageFilter<LabelImageType,LabelImageType> ChangeLabelImageFilterType;
typedef otb::TileImageFilter<LabelImageType> TileImageFilterType;
itkNewMacro(Self);
7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
itkTypeMacro(Merging, otb::Application);
private:
ChangeLabelImageFilterType::Pointer m_ChangeLabelFilter;
void DoInit() override
{
SetName("LSMSSmallRegionsMerging");
SetDescription("This application performs the third (optional) step of the exact Large-Scale Mean-Shift segmentation workflow [1].");
SetDocLongDescription("Given a segmentation result (can be the out output parameter of the"
" LSMSSegmentation application [2]) and the original image, it will"
" merge segments whose size in pixels is lower than minsize parameter"
" with the adjacent segments with the adjacent segment with closest"
" radiometry and acceptable size.\n\n"
"Small segments will be processed by increasing size: first all segments"
" for which area is equal to 1 pixel will be merged with adjacent"
" segments, then all segments of area equal to 2 pixels will be processed,"
" until segments of area minsize. For large images one can use the"
" tilesizex and tilesizey parameters for tile-wise processing, with the"
" guarantees of identical results.\n\n"
"The output of this application can be passed to the"
" LSMSVectorization application [3] to complete the LSMS workflow.");
SetDocLimitations("This application is part of the Large-Scale Mean-Shift segmentation"
" workflow (LSMS) and may not be suited for any other purpose. This"
" application is not compatible with in-memory connection since it does"
" its own internal streaming.");
SetDocAuthors("David Youssefi");
SetDocSeeAlso( "[1] Michel, J., Youssefi, D., & Grizonnet, M. (2015). Stable"
" mean-shift algorithm and its application to the segmentation of"
" arbitrarily large remote sensing images. IEEE Transactions on"
" Geoscience and Remote Sensing, 53(2), 952-964.\n"
"[2] LSMSegmentation\n"
"[3] LSMSVectorization");
AddDocTag(Tags::Segmentation);
AddDocTag(Tags::Deprecated);
AddDocTag("LSMS");
AddParameter(ParameterType_InputImage, "in", "Input image");
SetParameterDescription( "in", "The input image, containing initial spectral signatures corresponding to the segmented image (inseg)." );
AddParameter(ParameterType_InputImage, "inseg", "Segmented image");
SetParameterDescription( "inseg", "Segmented image where each pixel value is the unique integer label of the segment it belongs to." );
AddParameter(ParameterType_OutputImage, "out", "Output Image");
SetParameterDescription( "out", "The output image. The output image is the segmented image where the minimal segments have been merged. An ecoding of uint32 is advised." );
SetDefaultOutputPixelType("out",ImagePixelType_uint32);
AddParameter(ParameterType_Int, "minsize", "Minimum Segment Size");
SetParameterDescription("minsize", "Minimum Segment Size. If, after the segmentation, a segment is of size lower than this criterion, the segment is merged with the segment that has the closest sepctral signature.");
SetDefaultParameterInt("minsize", 50);
SetMinimumParameterIntValue("minsize", 0);
MandatoryOff("minsize");
AddParameter(ParameterType_Int, "tilesizex", "Size of tiles in pixel (X-axis)");
SetParameterDescription("tilesizex", "Size of tiles along the X-axis for tile-wise processing.");
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 for tile-wise processing.");
SetDefaultParameterInt("tilesizey", 500);
SetMinimumParameterIntValue("tilesizey", 1);
AddRAMParameter();
// Doc example parameter settings
SetDocExampleParameterValue("in","smooth.tif");
SetDocExampleParameterValue("inseg","segmentation.tif");
SetDocExampleParameterValue("out","merged.tif");
SetDocExampleParameterValue("minsize","20");
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
SetDocExampleParameterValue("tilesizex","256");
SetDocExampleParameterValue("tilesizey","256");
SetOfficialDocLink();
}
void DoUpdateParameters() override
{
}
void DoExecute() override
{
clock_t tic = clock();
unsigned int minSize = GetParameterInt("minsize");
unsigned long sizeTilesX = GetParameterInt("tilesizex");
unsigned long sizeTilesY = GetParameterInt("tilesizey");
//Acquisition of the input image dimensions
ImageType::Pointer imageIn = GetParameterImage("in");
imageIn->UpdateOutputInformation();
unsigned long sizeImageX = imageIn->GetLargestPossibleRegion().GetSize()[0],
sizeImageY = imageIn->GetLargestPossibleRegion().GetSize()[1];
unsigned int numberOfComponentsPerPixel = imageIn->GetNumberOfComponentsPerPixel();
LabelImageType::Pointer labelIn = GetParameterUInt32Image("inseg");
StatisticsImageFilterType::Pointer stats = StatisticsImageFilterType::New();
stats->SetInput(labelIn);
stats->GetStreamer()->SetAutomaticAdaptativeStreaming(GetParameterInt("ram"));
AddProcess(stats->GetStreamer(), "Retrieve region count...");
stats->Update();
unsigned int regionCount=stats->GetMaximum();
std::vector<unsigned int>nbPixels;
nbPixels.clear();
nbPixels.resize(regionCount+1);
for(LabelImagePixelType curLabel = 1; curLabel <= regionCount; ++curLabel)
nbPixels[curLabel] = 0;
ImageType::PixelType defaultValue(numberOfComponentsPerPixel);
defaultValue.Fill(0);
std::vector<ImageType::PixelType>sum(regionCount+1,defaultValue);
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);
//Sums calculation per label
otbAppLogINFO(<<"Sums calculation ...");
for(unsigned int row = 0; row < nbTilesY; row++)
for(unsigned int column = 0; column < nbTilesX; column++)
{
unsigned long startX = column*sizeTilesX;
unsigned long startY = row*sizeTilesY;
unsigned long sizeX = std::min(sizeTilesX,sizeImageX-startX);
unsigned long sizeY = std::min(sizeTilesY,sizeImageY-startY);
//Tiles extraction of the input image
MultiChannelExtractROIFilterType::Pointer imageROI = MultiChannelExtractROIFilterType::New();
imageROI->SetInput(imageIn);
imageROI->SetStartX(startX);
imageROI->SetStartY(startY);
imageROI->SetSizeX(sizeX);
imageROI->SetSizeY(sizeY);
imageROI->Update();
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
//Tiles extraction of the segmented image
ExtractROIFilterType::Pointer labelImageROI = ExtractROIFilterType::New();
labelImageROI->SetInput(labelIn);
labelImageROI->SetStartX(startX);
labelImageROI->SetStartY(startY);
labelImageROI->SetSizeX(sizeX);
labelImageROI->SetSizeY(sizeY);
labelImageROI->Update();
//Sums calculation for the mean calculation per label
LabelImageIterator itLabel( labelImageROI->GetOutput(), labelImageROI->GetOutput()->GetLargestPossibleRegion());
ImageIterator itImage( imageROI->GetOutput(), imageROI->GetOutput()->GetLargestPossibleRegion());
for (itLabel.GoToBegin(), itImage.GoToBegin(); !itLabel.IsAtEnd(); ++itLabel, ++itImage)
{
nbPixels[itLabel.Value()]++;
for(unsigned int comp = 0; comp<numberOfComponentsPerPixel; ++comp)
{
sum[itLabel.Value()][comp]+=itImage.Get()[comp];
}
}
}
//LUT creation for the final relabelling
std::vector<LabelImagePixelType> LUT;
LUT.clear();
LUT.resize(regionCount+1);
for(LabelImagePixelType curLabel = 1; curLabel <= regionCount; ++curLabel)
{
LUT[curLabel] = curLabel;
}
//Minimal size region suppression
otbAppLogINFO(<<"Building LUT for small regions merging ...");
for (unsigned int size=1; size<minSize; size++)
{
// LUTtmp creation in order to modify the LUT only at the end of the pass
std::vector<LabelImagePixelType> LUTtmp;
LUTtmp.clear();
LUTtmp.resize(regionCount+1);
for(LabelImagePixelType curLabel = 1; curLabel <= regionCount; ++curLabel)
{
LUTtmp[curLabel] = LUT[curLabel];
}
for(unsigned int row = 0; row < nbTilesY; row++)
{
for(unsigned int column = 0; column < nbTilesX; column++)
{
std::set<int> minLabel, edgeLabel, labelMerged;
std::map<int,std::set<int> > adjMap;
unsigned long startX = column*sizeTilesX, startY = row*sizeTilesY;
unsigned long sizeX = std::min(sizeTilesX+size+1,sizeImageX-startX),
sizeY = std::min(sizeTilesY+size+1,sizeImageY-startY);
ExtractROIFilterType::Pointer labelImageROI = ExtractROIFilterType::New();
labelImageROI->SetInput(labelIn);
labelImageROI->SetStartX(startX);
labelImageROI->SetStartY(startY);
labelImageROI->SetSizeX(sizeX);
labelImageROI->SetSizeY(sizeY);
labelImageROI->Update();
LabelImageType::IndexType pixelIndex;
//"Adjacency map" creation for the region with nbPixels=="size"
for(pixelIndex[0]=0; pixelIndex[0]<static_cast<long>(sizeX); ++pixelIndex[0])
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
for(pixelIndex[1]=0; pixelIndex[1]<static_cast<long>(sizeY); ++pixelIndex[1])
{
LabelImagePixelType curLabel = labelImageROI->GetOutput()->GetPixel(pixelIndex);
if(labelMerged.count(LUT[curLabel]))
{
edgeLabel.insert(LUT[curLabel]);
}
if((pixelIndex[0]==0)&&(startX!=0))
{
edgeLabel.insert(LUT[curLabel]);
}
if((pixelIndex[1]==0)&&(startY!=0))
{
edgeLabel.insert(LUT[curLabel]);
}
if(pixelIndex[0]==static_cast<long>(sizeX)-1)
{
if(startX!=(nbTilesX-1)*sizeTilesX) edgeLabel.insert(LUT[curLabel]);
}
else
{
++pixelIndex[0];
LabelImagePixelType adjLabelX = labelImageROI->GetOutput()->GetPixel(pixelIndex);
--pixelIndex[0];
if(LUT[adjLabelX]!=LUT[curLabel])
{
if((nbPixels[LUT[curLabel]]>0)&&(nbPixels[LUT[curLabel]]==size))
{
adjMap[LUT[curLabel]].insert(LUT[adjLabelX]);
minLabel.insert(LUT[curLabel]);
}
if((nbPixels[LUT[adjLabelX]]>0)&&(nbPixels[LUT[adjLabelX]]==size))
{
adjMap[LUT[adjLabelX]].insert(LUT[curLabel]);
minLabel.insert(LUT[adjLabelX]);
}
}
}
if(pixelIndex[1]==static_cast<long>(sizeY)-1)
{
if(startY!=(nbTilesY-1)*sizeTilesY) edgeLabel.insert(LUT[curLabel]);
}
else
{
++pixelIndex[1];
LabelImagePixelType adjLabelY = labelImageROI->GetOutput()->GetPixel(pixelIndex);
--pixelIndex[1];
if(LUT[adjLabelY]!=LUT[curLabel])
{
if((nbPixels[LUT[curLabel]]>0)&&(nbPixels[LUT[curLabel]]==size))
{
adjMap[LUT[curLabel]].insert(LUT[adjLabelY]);
minLabel.insert(LUT[curLabel]);
}
if((nbPixels[LUT[adjLabelY]]>0)&&(nbPixels[LUT[adjLabelY]]==size))
{
adjMap[LUT[adjLabelY]].insert(LUT[curLabel]);
minLabel.insert(LUT[adjLabelY]); }
}
}
}
//Searching the "nearest" region
for(std::set<int>::iterator itMinLabel=minLabel.begin(); itMinLabel!=minLabel.end(); ++itMinLabel)
{
LabelImagePixelType curLabel = *itMinLabel, adjLabel(0);
double err = itk::NumericTraits<double>::max();
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
if(edgeLabel.count(curLabel)==0)
{
if(nbPixels[curLabel]==size)
{
edgeLabel.insert(curLabel);
for(std::set<int>::iterator itAdjLabel=(adjMap[curLabel]).begin();
itAdjLabel!=(adjMap[curLabel]).end(); ++itAdjLabel)
{
double tmpError = 0;
LabelImagePixelType tmpLabel = *itAdjLabel;
if(tmpLabel!=curLabel)
{
for(unsigned int comp = 0; comp<numberOfComponentsPerPixel; ++comp)
{
double curComp = static_cast<double>(sum[curLabel][comp])/nbPixels[curLabel];
int tmpComp = static_cast<double>(sum[tmpLabel][comp])/nbPixels[tmpLabel];
tmpError += (curComp-tmpComp)*(curComp-tmpComp);
}
if(tmpError<err)
{
err = tmpError;
adjLabel = tmpLabel;
}
}
}
//Fusion of the two regions
if(adjLabel!=curLabel)
{
unsigned int curLabelLUT = curLabel, adjLabelLUT = adjLabel;
while(LUTtmp[curLabelLUT] != curLabelLUT)
{
curLabelLUT = LUTtmp[curLabelLUT];
}
while(LUTtmp[adjLabelLUT] != adjLabelLUT)
{
adjLabelLUT = LUTtmp[adjLabelLUT];
}
if(curLabelLUT < adjLabelLUT)
{
LUTtmp[adjLabelLUT] = curLabelLUT;
}
else
{
LUTtmp[LUTtmp[curLabelLUT]] = adjLabelLUT; LUTtmp[curLabelLUT] = adjLabelLUT;
}
}
}
}
}
for(LabelImagePixelType label = 1; label < regionCount+1; ++label)
{
LabelImagePixelType can = label;
while(LUTtmp[can] != can)
{
can = LUTtmp[can];
}
LUTtmp[label] = can;
}
}
}
for(LabelImagePixelType label = 1; label < regionCount+1; ++label)
{
LUT[label]=LUTtmp[label];
if((nbPixels[label]!=0)&&(LUT[label]!=label))
{
nbPixels[LUT[label]]+=nbPixels[label];
nbPixels[label]=0;
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452
for(unsigned int comp = 0; comp<numberOfComponentsPerPixel; ++comp)
{
sum[LUT[label]][comp]+=sum[label][comp];
}
}
}
}
//Relabelling
m_ChangeLabelFilter = ChangeLabelImageFilterType::New();
m_ChangeLabelFilter->SetInput(labelIn);
for(LabelImagePixelType label = 1; label<regionCount+1; ++label)
{
if(label!=LUT[label])
{
m_ChangeLabelFilter->SetChange(label,LUT[label]);
}
}
SetParameterOutputImage("out", m_ChangeLabelFilter->GetOutput());
clock_t toc = clock();
otbAppLogINFO(<<"Elapsed time: "<<(double)(toc - tic) / CLOCKS_PER_SEC<<" seconds");
}
};
}
}
OTB_APPLICATION_EXPORT(otb::Wrapper::LSMSSmallRegionsMerging)