An error occurred while loading the file. Please try again.
-
Pierre-Antoine Rouby authored2e158f4c
#ifndef __LSGRM_GRAPH_OPERATIONS_TXX#define __LSGRM_GRAPH_OPERATIONS_TXX//#include "lsgrmGraphOperations.h"//#include <unistd.h>#include <cstdio>namespace lsgrm{template<class TSegmenter>typename TSegmenter::ImageType::Pointer ReadImageRegion( typename TSegmenter::ImageType * inputPtr, typename TSegmenter::ImageType::RegionType region){ typedef typename otb::MultiChannelExtractROI<typename TSegmenter::ImageType::InternalPixelType, typename TSegmenter::ImageType::InternalPixelType> ExtractROIFilterType; typename ExtractROIFilterType::Pointer filter = ExtractROIFilterType::New(); filter->SetInput(inputPtr); filter->SetExtractionRegion(region); filter->SetReleaseDataFlag(true); filter->Update(); return filter->GetOutput();}template<class TSegmenter>typename TSegmenter::LabelImageType::PointerMergeAllGraphsAndAchieveSegmentation( const typename TSegmenter::ParamType& params, const float& threshold, std::vector<ProcessingTile>& tiles, const unsigned int nbTilesX, const unsigned int nbTilesY, const unsigned int imageWidth, const unsigned int imageHeight, const unsigned int imageBands, unsigned int numberOfIterations){ std::cout << "--- Graph aggregation...\n" << std::endl; TSegmenter segmenter; std::cout << "Reading graphs" << std::endl; for(unsigned int row = 0; row < nbTilesY; ++row) { for(unsigned int col = 0; col < nbTilesX; col++) { std::cout << "\tImporting graph of tile " << (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) << std::endl; InsertNodesFromTile<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], false); } } std::cout << "Removing duplicated nodes and updating neighbors" << std::endl; // TODO this might be parallelized... for(unsigned int row = 0; row < nbTilesY; ++row) { for(unsigned int col = 0; col < nbTilesX; col++) { std::cout << "Cleaning nodes of tile " << (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) << std::endl; std::unordered_map<long unsigned int, std::vector<typename TSegmenter::NodePointerType> > borderPixelMap; std::cout << "\tBuildBorderPixelMap..." << std::endl;7172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140
BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], row, col,
nbTilesX, nbTilesY, borderPixelMap, imageWidth);
std::cout << "\tRemoveDuplicatedNodes..." << std::endl;
RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);
std::cout << "\tUpdateNeighborsOfNoneDuplicatedNodes..." << std::endl;
UpdateNeighborsOfNoneDuplicatedNodes<TSegmenter>(borderPixelMap, imageWidth, imageHeight);
}
}
std::cout << "Achieve segmentation process" ;
// Segmentation of the graph
segmenter.SetImageWidth(imageWidth);
segmenter.SetImageHeight(imageHeight);
segmenter.SetNumberOfComponentsPerPixel(imageBands);
segmenter.SetParam(params);
segmenter.SetThreshold(threshold);
segmenter.SetDoFastSegmentation(false); // was true
segmenter.SetNumberOfIterations(numberOfIterations);
grm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
// // Write output graph to the output graph directory
// WriteGraph<TSegmenter>(segmenter.m_Graph, tmpDir, 0, 0);
// TODO: use std::sort to sort the graph before the polygon filling
typename TSegmenter::LabelImageType::SizeType outSize;
outSize[0] = imageWidth;
outSize[1] = imageHeight;
typedef StreamingGraphToImageFilter<typename TSegmenter::GraphType,
typename TSegmenter::LabelImageType> GraphToImageFilterType;
typename GraphToImageFilterType::Pointer graphToImageFilter = GraphToImageFilterType::New();
graphToImageFilter->SetGraph(segmenter.m_Graph);
graphToImageFilter->SetOutputSize(outSize);
typedef otb::ImageFileWriter<typename TSegmenter::LabelImageType> VWriterType;
typename VWriterType::Pointer wrt = VWriterType::New();
wrt->SetInput(graphToImageFilter->GetOutput());
wrt->SetNumberOfDivisionsStrippedStreaming(5);
wrt->SetFileName("/homeL/cresson/data/tmp/spot6_surrech_segm2.tif");
wrt->Update();
// Generate the label image
typename TSegmenter::LabelImageType::IndexType index;
index.Fill(0);
typename TSegmenter::LabelImageType::SizeType size;
size[0] = imageWidth;
size[1] = imageHeight;
typename TSegmenter::LabelImageType::RegionType region(index, size);
const typename TSegmenter::LabelImageType::InternalPixelType noDataLabel = 0;
typename TSegmenter::LabelImageType::Pointer labelImage = TSegmenter::LabelImageType::New();
labelImage->SetRegions(region);
labelImage->Allocate();
labelImage->FillBuffer(noDataLabel);
using LabelImageIterator = itk::ImageRegionIterator<typename TSegmenter::LabelImageType>;
LabelImageIterator it(labelImage, labelImage->GetLargestPossibleRegion());
// Start at 1 (value 0 can be used for invalid pixels)
unsigned int label = 1;
for(auto& node : segmenter.m_Graph.m_Nodes)
{
lp::CellLists borderPixels;
lp::ContourOperations::GenerateBorderCells(borderPixels, node->m_Contour, node->m_Id, imageWidth);
for (auto& pix: borderPixels)
141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210
{
index[0] = pix % imageWidth;
index[1] = pix / imageWidth;
labelImage->SetPixel(index, label);
}
++label;
}
// Re-order nodes labels (left->right to top->bottom)
vnl_vector<typename TSegmenter::LabelImageType::InternalPixelType> lut(label,noDataLabel);
label = 1;
for(it.GoToBegin();!it.IsAtEnd(); ++it)
{
unsigned int inputLabel = it.Get();
if (lut[ inputLabel ] == noDataLabel && inputLabel != noDataLabel)
{
lut[ inputLabel ] = label;
label++;
}
}
// Apply LUT
for(it.GoToBegin();!it.IsAtEnd(); ++it)
it.Set(lut[it.Get()]);
// Fill holes
typedef itk::GrayscaleFillholeImageFilter<typename TSegmenter::LabelImageType,
typename TSegmenter::LabelImageType> FillholeFilterType;
typename FillholeFilterType::Pointer fillFilter = FillholeFilterType::New();
fillFilter->SetInput(labelImage);
fillFilter->Update();
return fillFilter->GetOutput();
}
template<class TSegmenter>
long long unsigned int RunPartialSegmentation(const typename TSegmenter::ParamType& params,
const float& threshold,
const unsigned int niter,
std::vector<ProcessingTile>& tiles,
const unsigned int nbTilesX,
const unsigned int nbTilesY,
const unsigned int imageWidth,
const unsigned int imageHeight,
const unsigned int imageBands,
bool& isFusion)
{
long long unsigned int accumulatedMemory = 0;
isFusion = false;
const unsigned int numberOfNeighborLayers = static_cast<unsigned int>(pow(2, niter + 1) - 2);
std::cout << "--- Running partial segmentations...\nNumber of neighbor layers " << numberOfNeighborLayers << std::endl;
for(unsigned int row = 0; row < nbTilesY; ++row)
{
for(unsigned int col = 0; col < nbTilesX; col++)
{
#ifdef OTB_USE_MPI
if (MyTurn(row*nbTilesX + col))
#endif
{
// Get the current tile
std::cout << "Processing tile " << row << ", " << col << std::endl;
ProcessingTile currentTile = tiles[row*nbTilesX + col];
// Load the graph
std::cout << "\tLoad graph..." << std::endl;
TSegmenter segmenter;
ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
// Add stability margin to the graph
{
std::cout << "\tAdd stability margin..." << std::endl;
AddStabilityMargin<TSegmenter>(segmenter.m_Graph, tiles,
row, col, nbTilesX, nbTilesY);
std::cout << "\tBuild border pixel map..." << std::endl;
std::unordered_map<long unsigned int, std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, currentTile, row, col,
nbTilesX, nbTilesY, borderPixelMap, imageWidth);
std::cout << "\tRemove duplicated nodes..." << std::endl;
RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);
std::cout << "\tUpdate neighbors.." << std::endl;
UpdateNeighborsOfNoneDuplicatedNodes<TSegmenter>(borderPixelMap,
imageWidth,
imageHeight);
std::cout << "\tRemove useless.." << std::endl;
RemoveUselessNodes<TSegmenter>(currentTile, segmenter.m_Graph,
imageWidth, numberOfNeighborLayers);
}
// Segmentation of the graph
segmenter.SetImageWidth(imageWidth);
segmenter.SetImageHeight(imageHeight);
segmenter.SetNumberOfComponentsPerPixel(imageBands);
segmenter.SetParam(params);
segmenter.SetThreshold(threshold);
segmenter.SetDoFastSegmentation(false);
segmenter.SetNumberOfIterations(niter);
std::cout << "\tPartial segmentation";
auto merge = grm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
if(merge == true)
isFusion = true;
// Remove unstable segments
std::cout << "\tRemove unstable segments..." << std::endl;
RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, currentTile, imageWidth);
// Retrieve the amount of memory to store this graph
std::cout << "\tGet graph memory..." << std::endl;
accumulatedMemory += segmenter.GetGraphMemory();
// Write graph to temporay directory
std::cout << "\tWrite graph..." << std::endl;
WriteGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
}
}
}
#ifdef OTB_USE_MPI
otb::MPIConfig::Instance()->barrier();
#endif
std::cout << "Add stability margins to graph for the next round..."<< std::endl;
// During this step we extract the stability margin for the next round
for(unsigned int row = 0; row < nbTilesY; ++row)
{
for(unsigned int col = 0; col < nbTilesX; col++)
{
#ifdef OTB_USE_MPI
if (MyTurn(row*nbTilesX + col))
#endif
281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350
{
// Get current tile
ProcessingTile currentTile = tiles[row*nbTilesX + col];
// Load the graph
typename TSegmenter::GraphType graph;
ReadGraph<TSegmenter>(graph, currentTile.nodeFileName, currentTile.edgeFileName);
// Extract stability margin for all borders different from 0 imageWidth-1 and imageHeight-1
// and write them to the stability margin
{
std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;
DetectBorderNodes<TSegmenter>(graph, currentTile,
borderNodeMap, imageWidth, imageHeight);
ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);
WriteStabilityMargin<TSegmenter>(borderNodeMap, currentTile.nodeMarginFileName, currentTile.edgeMarginFileName);
}
std::cout << "Process finished on tile " << (row*nbTilesX + col) << std::endl;
}
}
}
std::cout << std::endl;
return accumulatedMemory;
}
template<class TSegmenter>
void RemoveUselessNodes(ProcessingTile& tile,
typename TSegmenter::GraphType& graph,
const unsigned int imageWidth,
const unsigned int numberOfLayers)
{
using NPtr = typename TSegmenter::NodePointerType;
using UI = unsigned int;
unsigned int rowPixel, colPixel;
std::unordered_map<NPtr, UI> marginNodes;
for(auto& node : graph.m_Nodes)
{
if(node->m_Bbox.m_UX > tile.columns[0] &&
node->m_Bbox.m_UY > tile.rows[0] &&
node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[1] &&
node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[1])
continue;
else if(node->m_Bbox.m_UX > tile.columns[1] ||
node->m_Bbox.m_UY > tile.rows[1] ||
node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[0] ||
node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[0])
continue;
else
{
lp::CellLists borderpixels;
lp::ContourOperations::GenerateBorderCells(borderpixels, node->m_Contour, node->m_Id, imageWidth);
for(auto& pix : borderpixels)
{
rowPixel = pix / imageWidth;
colPixel = pix % imageWidth;
if(rowPixel == tile.rows[0] || rowPixel == tile.rows[1])
{
if(colPixel >= tile.columns[0] && colPixel <= tile.columns[1])
{
marginNodes.insert(std::pair<NPtr, UI>(node, 0));
break;
}
}
351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
else if(colPixel == tile.columns[0] || colPixel == tile.columns[1])
{
if(rowPixel >= tile.rows[0] && rowPixel <= tile.rows[1])
{
marginNodes.insert(std::pair<NPtr, UI>(node, 0));
break;
}
}
else
continue;
}
}
}
ExtractStabilityMargin<TSegmenter>(marginNodes, numberOfLayers);
for(auto& node : graph.m_Nodes)
{
if(node->m_Bbox.m_UX > tile.columns[0] &&
node->m_Bbox.m_UY > tile.rows[0] &&
node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[1] &&
node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[1])
continue;
else if(marginNodes.find(node) == marginNodes.end())
{
RemoveEdgeToUnstableNode<TSegmenter>(node);
node->m_Expired = true;
}
else
continue;
}
grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
}
template<class TSegmenter>
void UpdateNeighborsOfNoneDuplicatedNodes(std::unordered_map<long unsigned int,
std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
const unsigned int imageWidth,
const unsigned int imageHeight)
{
using EdgeType = typename TSegmenter::EdgeType;
unsigned int boundary;
for(auto& pn : borderPixelMap)
{
long int neighborhood[4];
grm::FOURNeighborhood(neighborhood, pn.first, imageWidth, imageHeight);
for(short j = 0; j < 4; ++j)
{
if(neighborhood[j] > -1)
{
auto isNeigh = borderPixelMap.find(neighborhood[j]);
if(isNeigh != borderPixelMap.end())
{
auto currNode = pn.second.front();
auto neigh = isNeigh->second.front();
if(currNode != neigh)
{
auto toNeigh = grm::GraphOperations<TSegmenter>::FindEdge(currNode, neigh);
if(toNeigh == currNode->m_Edges.end())
{
boundary = 0;
lp::CellLists borderPixels;
lp::ContourOperations::GenerateBorderCells(borderPixels, currNode->m_Contour, currNode->m_Id, imageWidth);
421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490
for(auto pix : borderPixels)
{
if(borderPixelMap.find(pix) != borderPixelMap.end())
{
long int pixNeighborhood[4];
grm::FOURNeighborhood(pixNeighborhood, pix, imageWidth, imageHeight);
for(short k = 0; k < 4; k++)
{
if(pixNeighborhood[k] > -1)
{
auto isNeighPix = borderPixelMap.find(pixNeighborhood[k]);
if(isNeighPix != borderPixelMap.end() && isNeighPix->second.front() == neigh)
boundary++;
}
}
}
}
currNode->m_Edges.push_back(EdgeType(neigh, 0, boundary));
neigh->m_Edges.push_back(EdgeType(currNode, 0, boundary));
}
}
}
}
}
}
}
template<class TSegmenter>
void RemoveDuplicatedNodes(std::unordered_map<long unsigned int,
std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
typename TSegmenter::GraphType& graph,
const unsigned int imageWidth)
{
using EdgeType = typename TSegmenter::EdgeType;
unsigned int boundary;
// Explore the border nodes
for(auto& pn : borderPixelMap)
{
// no valid nodes are those who have already been considered.
if(pn.second.size() > 1)
{
auto refNode = pn.second.front();
// Explore duplicated nodes
for(auto nit = pn.second.begin() + 1; nit != pn.second.end(); ++nit)
{
// Explore their edges
for(auto& edg : (*nit)->m_Edges)
{
auto neighNit = edg.GetRegion();
auto toNit = grm::GraphOperations<TSegmenter>::FindEdge(neighNit, *nit);
assert(toNit != edg.GetRegion()->m_Edges.end());
boundary = edg.m_Boundary;
neighNit->m_Edges.erase(toNit);
auto toRefNode = grm::GraphOperations<TSegmenter>::FindEdge(neighNit, refNode);
if(toRefNode == neighNit->m_Edges.end())
{
// Create an edge neighNit -> refNode
neighNit->m_Edges.push_back(EdgeType(refNode, 0, boundary));
// Create an edge refNode -> neighNit
491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560
refNode->m_Edges.push_back(EdgeType(neighNit, 0, boundary));
}
}
(*nit)->m_Expired = true;
}
lp::CellLists borderPixels;
lp::ContourOperations::GenerateBorderCells(borderPixels, refNode->m_Contour, refNode->m_Id, imageWidth);
for(auto& pix : borderPixels)
{
if(borderPixelMap.find(pix) != borderPixelMap.end())
{
borderPixelMap[pix].clear();
borderPixelMap[pix].push_back(refNode);
}
}
}
}
grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
}
// Build the map assigning for each border pixel on the borders the node it belongs to.
template<class TSegmenter>
void BuildBorderPixelMap(typename TSegmenter::GraphType& graph,
ProcessingTile& tile,
const unsigned int rowTile,
const unsigned int colTile,
const unsigned int nbTilesX,
const unsigned int nbTilesY,
std::unordered_map<long unsigned int,
std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
const unsigned int imageWidth)
{
unsigned int rowPixel, colPixel;
unsigned int rowMin = (tile.rows[0] > 0) ? tile.rows[0] - 1 : tile.rows[0];
unsigned int rowMax = tile.rows[1] + 1;
unsigned int colMin = (tile.columns[0] > 0) ? tile.columns[0] - 1 : tile.columns[0];
unsigned int colMax = tile.columns[1] + 1;
for(auto& node : graph.m_Nodes)
{
if(node->m_Bbox.m_UX > tile.columns[0] &&
node->m_Bbox.m_UY > tile.rows[0] &&
node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[1] &&
node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[1])
{
continue;
}
else
{
lp::CellLists borderPixels;
lp::ContourOperations::GenerateBorderCells(borderPixels, node->m_Contour, node->m_Id, imageWidth);
for(auto& pix : borderPixels)
{
rowPixel = pix / imageWidth;
colPixel = pix % imageWidth;
if(rowTile > 0 && (rowPixel == tile.rows[0] || rowPixel == rowMin))
{
borderPixelMap[pix].push_back(node);
}
else if(colTile < nbTilesX - 1 && ( colPixel == tile.columns[1] || colPixel == colMax))
{
borderPixelMap[pix].push_back(node);
}
else if(rowTile < nbTilesY - 1 && (rowPixel == tile.rows[1] || rowPixel == rowMax))
561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630
{
borderPixelMap[pix].push_back(node);
}
else if(colTile > 0 && ( colPixel == tile.columns[0] || colPixel == colMin))
{
borderPixelMap[pix].push_back(node);
}
else
continue;
}
}
}
}
template<class TSegmenter>
void InsertNodesFromTile(typename TSegmenter::GraphType& graph,
ProcessingTile& tile, bool margin)
{
typename TSegmenter::GraphType subgraph;
if (margin)
{
ReadGraph<TSegmenter>(subgraph, tile.nodeMarginFileName, tile.edgeMarginFileName);
}
else
{
ReadGraph<TSegmenter>(subgraph, tile.nodeFileName, tile.edgeFileName);
}
graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());
}
template<class TSegmenter>
void AddStabilityMargin(typename TSegmenter::GraphType& graph,
std::vector<ProcessingTile>& tiles,
const unsigned int row,
const unsigned int col,
const unsigned int nbTilesX,
const unsigned int nbTilesY)
{
// Margin to retrieve at top
if(row > 0)
{
InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + col]);
}
// Margin to retrieve at right
if(col < nbTilesX - 1)
{
InsertNodesFromTile<TSegmenter>(graph, tiles[row * nbTilesX + (col+1)]);
}
// Margin to retrieve at bottom
if(row < nbTilesY - 1)
{
InsertNodesFromTile<TSegmenter>(graph, tiles[(row+1) * nbTilesX + col]);
}
// Margin to retrieve at left
if(col > 0)
{
InsertNodesFromTile<TSegmenter>(graph, tiles[row * nbTilesX + (col-1)]);
}
// Margin to retrieve at top right
if(row > 0 && col < nbTilesX - 1)
{
InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + (col+1)]);
}
631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700
// Margin to retrieve at bottom right
if(row < nbTilesY - 1 && col < nbTilesX - 1)
{
InsertNodesFromTile<TSegmenter>(graph, tiles[(row+1) * nbTilesX + (col+1)]);
}
// Margin to retrieve at bottom left
if(row < nbTilesY - 1 && col > 0)
{
InsertNodesFromTile<TSegmenter>(graph, tiles[(row+1) * nbTilesX + (col-1)]);
}
// Margin to retrieve at top left
if(row > 0 && col > 0)
{
InsertNodesFromTile<TSegmenter>(graph, tiles[(row-1) * nbTilesX + (col-1)]);
}
}
template<class TSegmenter>
long long unsigned int RunFirstPartialSegmentation(
typename TSegmenter::ImageType * inputPtr,
const typename TSegmenter::ParamType& params,
const float& threshold,
const unsigned int niter,
const unsigned int niter2,
std::vector<ProcessingTile>& tiles,
const unsigned int nbTilesX,
const unsigned int nbTilesY,
const unsigned int tileWidth,
const unsigned int tileHeight,
bool& isFusion)
{
using ImageType = typename TSegmenter::ImageType;
const unsigned int imageWidth = inputPtr->GetLargestPossibleRegion().GetSize()[0];
const unsigned int imageHeight = inputPtr->GetLargestPossibleRegion().GetSize()[1];
long long unsigned int accumulatedMemory = 0;
isFusion = false;
const unsigned int numberOfNeighborLayers = static_cast<unsigned int>(pow(2, niter2 + 1) - 2);
std::cout << "--- Running fist partial segmentation...\nNumber of neighbor layers " << numberOfNeighborLayers << std::endl;
for(unsigned int row = 0; row < nbTilesY; ++row)
{
for(unsigned int col = 0; col < nbTilesX; col++)
{
#ifdef OTB_USE_MPI
if (MyTurn(row*nbTilesX + col))
#endif
{
// Reading images
ProcessingTile currentTile = tiles[row*nbTilesX + col];
std::cout << "Processing tile " << (row*nbTilesX + col) << " / " << (nbTilesX*nbTilesY) <<
" (" << col << ", " << row << ")" <<
" start: [" << currentTile.region.GetIndex()[0] << ", " << currentTile.region.GetIndex()[1] <<
"] size: [" << currentTile.region.GetSize()[0] << ", " << currentTile.region.GetSize()[1] << "]" << std::endl;
typename ImageType::Pointer imageTile = ReadImageRegion<TSegmenter>(inputPtr, currentTile.region);
// Segmenting image
std::cout << "\tSegmenting";
TSegmenter segmenter;
segmenter.SetParam(params);
segmenter.SetThreshold(threshold);
segmenter.SetDoFastSegmentation(false);
701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770
segmenter.SetNumberOfIterations(niter);
segmenter.SetInput(imageTile);
segmenter.Update();
if(segmenter.GetComplete() == false)
isFusion = true;
// Rescale the graph to be in the reference of the image
std::cout << "\tRescaling graph..." << std::endl;
RescaleGraph<TSegmenter>(segmenter.m_Graph,
currentTile,
row,
col,
tileWidth,
tileHeight,
imageWidth);
// Remove unstable segments
std::cout << "\tRemoving unstable segments..." << std::endl;
RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, currentTile, imageWidth);
// Retrieve the amount of memory to store this graph
std::cout << "\tRetrieving graph memory..." << std::endl;
accumulatedMemory += segmenter.GetGraphMemory();
// Write graph to temporay directory
std::cout << "\tWriting graph..." << std::endl;
WriteGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
// Extract stability margin for all borders different from 0 imageWidth-1 and imageHeight -1
// and write them to the stability margin
std::cout << "\tComputing stability margin..." << std::endl;
{
std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;
DetectBorderNodes<TSegmenter>(segmenter.m_Graph, currentTile,
borderNodeMap, imageWidth, imageHeight);
ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);
WriteStabilityMargin<TSegmenter>(borderNodeMap, currentTile.nodeMarginFileName, currentTile.edgeMarginFileName);
}
}
} // for each col
} // for each row
return accumulatedMemory;
}
template<class TSegmenter>
void ExtractStabilityMargin(std::unordered_map<typename TSegmenter::NodePointerType, unsigned int>& nodeMap,
const unsigned int pmax)
{
std::vector<typename TSegmenter::NodePointerType> startingNodes;
startingNodes.reserve(nodeMap.size());
for(auto& kv: nodeMap)
startingNodes.push_back(kv.first);
for(auto& n : startingNodes)
{
ExploreDFS<TSegmenter>(n, 0, nodeMap, pmax);
}
}
template<class TSegmenter>
void ExploreDFS(typename TSegmenter::NodePointerType s,
const unsigned int p,
std::unordered_map<typename TSegmenter::NodePointerType, unsigned int>& Cb,
const unsigned int pmax)
{
771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840
if(p > pmax)
return;
else
{
if(Cb.find(s) != Cb.end())
{
if(p <= Cb[s])
{
Cb[s] = p;
for(auto edg : s->m_Edges)
{
ExploreDFS<TSegmenter>(edg.GetRegion(), p + 1, Cb, pmax);
}
}
else
return;
}
else
{
Cb[s] = p;
for(auto edg : s->m_Edges)
{
ExploreDFS<TSegmenter>(edg.GetRegion(), p + 1, Cb, pmax);
}
}
}
}
template<class TSegmenter>
void DetectBorderNodes(typename TSegmenter::GraphType& graph,
ProcessingTile& tile,
std::unordered_map<typename TSegmenter::NodePointerType, unsigned int>& borderNodeMap,
const unsigned int imageWidth,
const unsigned int imageHeight)
{
using NP = typename TSegmenter::NodePointerType;
using Uint = unsigned int;
unsigned int rowPixel, colPixel;
for(auto& node : graph.m_Nodes)
{
if(node->m_Bbox.m_UX > tile.columns[0] &&
node->m_Bbox.m_UY > tile.rows[0] &&
node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[1] &&
node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[1])
continue;
else
{
lp::CellLists borderpixels;
lp::ContourOperations::GenerateBorderCells(borderpixels, node->m_Contour, node->m_Id, imageWidth);
for(auto& pix : borderpixels)
{
rowPixel = pix / imageWidth;
colPixel = pix % imageWidth;
if(tile.rows[0] > 0 && rowPixel == tile.rows[0])
{
borderNodeMap.insert(std::pair<NP, Uint>(node, 0));
break;
}
else if(tile.columns[1] < imageWidth - 1 && colPixel == tile.columns[1])
{
borderNodeMap.insert(std::pair<NP, Uint>(node, 0));
break;
}
else if(tile.rows[1] < imageHeight - 1 && rowPixel == tile.rows[1])
{
borderNodeMap.insert(std::pair<NP, Uint>(node, 0));
break;
841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910
}
else if(tile.columns[0] > 0 && colPixel == tile.columns[0])
{
borderNodeMap.insert(std::pair<NP, Uint>(node, 0));
break;
}
else continue;
}
}
}
}
// Generic
template<class TSegmenter>
void ReadGraph(TSegmenter& segmenter,
const std::string& nodesPath,
const std::string& edgesPath)
{
FILE * nodeStream = fopen(nodesPath.c_str(), "rb");
assert(nodeStream != NULL);
FILE * edgeStream = fopen(edgesPath.c_str(), "rb");
assert(edgeStream != NULL);
segmenter.ReadGraph(nodeStream, edgeStream);
fclose(nodeStream);
fclose(edgeStream);
}
template<class TSegmenter>
void ReadGraph(typename TSegmenter::GraphType& graph,
const std::string& nodesPath,
const std::string& edgesPath)
{
TSegmenter segmenter;
ReadGraph<TSegmenter>(segmenter, nodesPath, edgesPath);
graph = segmenter.m_Graph;
}
// Write stability margin
template<class TSegmenter>
void WriteStabilityMargin(std::unordered_map<
typename TSegmenter::NodePointerType,
unsigned int>& stabilityMargin,
const std::string& nodesPath,
const std::string& edgesPath)
{
FILE * nodeStream = fopen(nodesPath.c_str(), "wb");
assert(nodeStream != NULL);
FILE * edgeStream = fopen(edgesPath.c_str(), "wb");
assert(edgeStream != NULL);
// Write number of nodes
std::size_t size = stabilityMargin.size();
fwrite(&size, sizeof(size), 1, nodeStream);
TSegmenter seg;
for(auto& kv : stabilityMargin)
{
auto node = kv.first;
seg.WriteNode(node, nodeStream);
// Write only edges pointing to nodes which are in the stability margin.
fwrite(&(node->m_Id), sizeof(node->m_Id), 1, edgeStream);
std::size_t edgeSize = node->m_Edges.size();
for(auto& edg : node->m_Edges)
{
if(stabilityMargin.find(edg.GetRegion()) == stabilityMargin.end())
edgeSize--;
911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961962963964965966967968969970971972973974975976977978979980
}
fwrite(&(edgeSize), sizeof(edgeSize), 1, edgeStream);
for(auto& edg : node->m_Edges)
{
if(stabilityMargin.find(edg.GetRegion()) != stabilityMargin.end())
{
seg.WriteEdge(edg, edgeStream);
}
}
}
fclose(nodeStream);
fclose(edgeStream);
}
// Write the graph
template<class TSegmenter>
void WriteGraph(typename TSegmenter::GraphType& graph,
const std::string& nodeFile,
const std::string& edgeFile)
{
FILE * nodeStream = fopen(nodeFile.c_str(), "wb");
assert(nodeStream != NULL);
FILE * edgeStream = fopen(edgeFile.c_str(), "wb");
assert(edgeStream != NULL);
// Write number of nodes
std::size_t size = graph.m_Nodes.size();
fwrite(&size, sizeof(size), 1, nodeStream);
TSegmenter seg;
for(auto& node : graph.m_Nodes)
{
seg.WriteNode(node, nodeStream);
// Write edges
fwrite(&(node->m_Id), sizeof(node->m_Id), 1, edgeStream);
std::size_t edgeSize = node->m_Edges.size();
fwrite(&(edgeSize), sizeof(edgeSize), 1, edgeStream);
for(auto& edg : node->m_Edges)
{
seg.WriteEdge(edg, edgeStream);
}
}
fclose(nodeStream);
fclose(edgeStream);
}
template<class TSegmenter>
void RemoveUnstableSegments(typename TSegmenter::GraphType& graph,
ProcessingTile& tile,
const unsigned int imageWidth)
{
unsigned int rowPixel, colPixel;
bool stable;
for(auto& node : graph.m_Nodes)
{
if(node->m_Bbox.m_UX >= tile.columns[0] &&
node->m_Bbox.m_UY >= tile.rows[0] &&
node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 <= tile.columns[1] &&
node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 <= tile.rows[1])
continue;
else if(node->m_Bbox.m_UX > tile.columns[1] ||
node->m_Bbox.m_UY > tile.rows[1] ||
node->m_Bbox.m_UX + node->m_Bbox.m_W - 1 < tile.columns[0] ||
node->m_Bbox.m_UY + node->m_Bbox.m_H - 1 < tile.rows[0])
981982983984985986987988989990991992993994995996997998999100010011002100310041005100610071008100910101011101210131014101510161017101810191020102110221023102410251026102710281029103010311032103310341035103610371038103910401041104210431044104510461047104810491050
{
node->m_Expired = true;
RemoveEdgeToUnstableNode<TSegmenter>(node);
}
else
{
lp::CellLists borderpixels;
lp::ContourOperations::GenerateBorderCells(borderpixels, node->m_Contour, node->m_Id, imageWidth);
stable = false;
for(auto& pix : borderpixels)
{
rowPixel = pix / imageWidth;
colPixel = pix % imageWidth;
if(rowPixel >= tile.rows[0] &&
rowPixel <= tile.rows[1] &&
colPixel >= tile.columns[0] &&
colPixel <= tile.columns[1])
{
stable = true;
break;
}
}
if(!stable)
{
node->m_Expired = true;
RemoveEdgeToUnstableNode<TSegmenter>(node);
}
}
}
grm::GraphOperations<TSegmenter>::RemoveExpiredNodes(graph);
}
template<class TSegmenter>
void RemoveEdgeToUnstableNode(typename TSegmenter::NodePointerType nodePtr)
{
for(auto& edg : nodePtr->m_Edges)
{
auto nodeNeighbor = edg.GetRegion();
auto EdgeToNode = grm::GraphOperations<TSegmenter>::FindEdge(nodeNeighbor, nodePtr);
assert(EdgeToNode != nodeNeighbor->m_Edges.end());
nodeNeighbor->m_Edges.erase(EdgeToNode);
}
}
template<class TSegmenter>
void RescaleGraph(typename TSegmenter::GraphType& graph,
ProcessingTile& tile,
const unsigned int rowTile,
const unsigned int colTile,
const unsigned int tileWidth,
const unsigned int tileHeight,
const unsigned int imageWidth)
{
unsigned int rowNodeTile, colNodeTile;
unsigned int rowNodeImg, colNodeImg;
for(auto& node : graph.m_Nodes)
{
// Start pixel index of the node (in the tile)
rowNodeTile = node->m_Id / tile.region.GetSize()[0];
colNodeTile = node->m_Id % tile.region.GetSize()[0];
// Start pixel index of the node (in the image)
rowNodeImg = rowTile * tileHeight + rowNodeTile - tile.margin[0];
colNodeImg = colTile * tileWidth + colNodeTile - tile.margin[3];
node->m_Id = rowNodeImg * imageWidth + colNodeImg;
1051105210531054105510561057105810591060
// Change also its bounding box
node->m_Bbox.m_UX = colTile * tileWidth + node->m_Bbox.m_UX - tile.margin[3];
node->m_Bbox.m_UY = rowTile * tileHeight + node->m_Bbox.m_UY - tile.margin[0];
}
}
}
#endif