#ifndef __LSGRM_GRAPH_OPERATIONS_TXX
#define __LSGRM_GRAPH_OPERATIONS_TXX
//#include "lsgrmGraphOperations.h"
//#include <unistd.h>
#include <fstream>
#include <cstdio>
#include <sys/stat.h>

namespace lsgrm
{

bool file_exists(const std::string& fileName)
{
    std::ifstream infile(fileName.c_str());
    return infile.good();
}

time_t last_mod_time(const std::string& fileName)
{
	struct stat result;
	if (!stat(fileName.c_str(), &result)) return result.st_mtime;
	else return -1;
}

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::GraphType
MergeAllGraphsAndAchieveSegmentation(
    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<std::size_t,
      std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;

      std::cout << "\tBuildBorderPixelMap..." << std::endl;

      BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], row, col,
          nbTilesX, nbTilesY, borderPixelMap, imageWidth, true);

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

  // patch raf (sort using node id directly --> seems ok)
  // Sort the nodes top-->down and left-->right
  /*
  std::sort(segmenter.m_Graph.m_Nodes.begin(), segmenter.m_Graph.m_Nodes.end(),
      [imageWidth](const auto & a, const auto & b) -> bool
  {
    lp::CellLists borderPixelsA, borderPixelsB;
    lp::ContourOperations::GenerateBorderCells(borderPixelsA, a->m_Contour, a->m_Id, imageWidth);
    lp::ContourOperations::GenerateBorderCells(borderPixelsB, b->m_Contour, b->m_Id, imageWidth);

	// patch raf (were unsigned int)
    std::size_t offA = 0;
    std::size_t offB = 0;
	// end patch
    for (auto& pix: borderPixelsA)
      if (pix>offA)
        offA=pix;

    for (auto& pix: borderPixelsB)
      if (pix>offB)
        offB=pix;

      return offA > offB;
  });
  */
  std::sort(segmenter.m_Graph.m_Nodes.begin(), segmenter.m_Graph.m_Nodes.end(),
      [](const auto & a, const auto & b) -> bool
  {
      return a->m_Id < b->m_Id;
  });
  // end patch

  return segmenter.m_Graph;
}

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,
	bool resume,
	unsigned int& nextTile)
{
  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++)
      {
        {
		// Get the current tile
        std::cout << "Processing tile " << row << ", " << col << std::endl;
        unsigned int currTile = row*nbTilesX + col;
        ProcessingTile currentTile = tiles[currTile];
		
		// If resume mode, get accumulated memory on previous tiles
		if (resume && currTile < nextTile ) {
			// Load the graph
			std::cout << "\tResuming graph..." << std::endl;
			TSegmenter segmenter;
			ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
			// Retrieve the amount of memory to store this graph
			std::cout << "\tGet graph memory..." << std::endl;
			accumulatedMemory += segmenter.GetGraphMemory();
			continue;
		}
		
        // Load the graph
        std::cout << "\tLoad graph..." << std::endl;
        TSegmenter segmenter;
        ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);

        // 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<std::size_t, 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);
		
        }
      }
    }
    
    if (resume)
		nextTile = 0;


  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++)
      {
        {
        // 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;
            }
          }
        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<std::size_t,
    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)
    {
    NeighIDType 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);

              for(auto pix : borderPixels)
                {
                if(borderPixelMap.find(pix) != borderPixelMap.end())
                  {
                  NeighIDType 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<std::size_t,
    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
            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<std::size_t,
    std::vector<typename TSegmenter::NodePointerType> >& borderPixelMap,
    const unsigned int imageWidth,
	bool merging)
{
  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;
      }
	// patch raf (exclude nodes all outside tiles - on merging stage only)
	else if(merging && (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;
      }  
	// end patch
    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))
          {
          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)]);
    }

  // 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,
    bool resume,
    unsigned int& nextTile)
{

  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;
  bool accomplished = true;

  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;

  time_t lastTime = -1;

  for(unsigned int row = 0; row < nbTilesY; ++row)
    {
    for(unsigned int col = 0; col < nbTilesX; col++)
      {
        {
        // Reading images
        ProcessingTile currentTile = tiles[row*nbTilesX + col];
        
        if (resume && file_exists(currentTile.nodeFileName) && file_exists(currentTile.edgeFileName))
		{
			int lastMod = last_mod_time(currentTile.nodeFileName);
			if (lastMod > lastTime) {
				lastTime = lastMod;
				nextTile = row*nbTilesX + col + 1;
			}
			std::cout << "\tResuming graph..." << std::endl;
			TSegmenter segmenter;
			ReadGraph<TSegmenter>(segmenter.m_Graph, currentTile.nodeFileName, currentTile.edgeFileName);
			// Retrieve the amount of memory to store this graph
			std::cout << "\tGet graph memory..." << std::endl;
			accumulatedMemory += segmenter.GetGraphMemory();
			if(segmenter.GetComplete() == false)
				isFusion = true;
			continue;
		}
		
		accomplished = false;
		
        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);
        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
    
    if (resume && !accomplished) {
		std::cout << "\tSegmentation recovered during first partial segmentation." << std::endl;
		nextTile = nbTilesX*nbTilesY;
	}
	if (resume && accomplished && nextTile == nbTilesX*nbTilesY)
		nextTile = 0;

  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)
{
  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;
          }
        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--;
      }
    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])
      {
      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 = (std::size_t)rowNodeImg * (std::size_t)imageWidth + (std::size_t)colNodeImg;

    // 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