diff --git a/Applications/CMakeLists.txt b/Applications/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c04c7727f5decd2379b6f72cbe49b60cfedc935a
--- /dev/null
+++ b/Applications/CMakeLists.txt
@@ -0,0 +1,19 @@
+#=========================================================================
+
+#  Program: Large Scale Generic Region Merging Library (LSGRM)
+#  Language: C++
+#  author: Lassalle Pierre
+
+
+
+#  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
+
+#  See grmlib-copyright.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.
+
+#=========================================================================
+add_executable(Test lsgrmProto.cxx)
+target_link_libraries(Test OTBLSGRM)
diff --git a/Applications/lsgrmProto.cxx b/Applications/lsgrmProto.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..8a8d7e22ab84129642971df00f2ae9510ae066a9
--- /dev/null
+++ b/Applications/lsgrmProto.cxx
@@ -0,0 +1,61 @@
+#include <iostream>
+#include "lsrmBaatzSegmenter.h"
+#include "lsgrmController.h"
+
+int main(int argc, char *argv[])
+{
+	if(argc != 8)
+	{
+		// img/test.tif tiles/ 125 125 4 tmp/
+		std::cerr << "[input image] [tile directory] [tile width] [tile height] [number of first iterations] [temporary directory] [output directory]"
+				  << std::endl;
+		return 1;
+	}
+
+	/* Parse command line arguments */
+	const std::string imagePath = argv[1];
+	std::string tileDir = argv[2]; // May be corrected if badly written.
+	const unsigned int tileWidth = atoi(argv[3]);
+	const unsigned int tileHeight = atoi(argv[4]);
+	const unsigned int niter = atoi(argv[5]);
+	std::string tmpDir = argv[6]; // May be corrected if badly written.
+	std::string outDir = argv[7];
+
+	/*
+	  To add:
+	  internal memory available
+	  If we have to do the image division
+	  if we have to clean up the directory
+	  the output directory in case the global graph cannot fit in memory
+	  
+	 */
+	
+	using ImageType = otb::VectorImage<float, 2>;
+	using SegmenterType = lsrm::BaatzSegmenter<ImageType>;
+	using ControllerType = lsgrm::Controller<SegmenterType>;
+
+	
+	
+	ControllerType controller;
+	controller.SetInputImage(imagePath);
+	controller.SetTileDirectory(tileDir);
+	controller.SetTemporaryDirectory(tmpDir);
+	controller.SetOutputGraphDirectory(outDir);
+
+	// Memory configuration
+	controller.SetInternalMemoryAvailable(512ul);
+	controller.SetTileWidth(500);
+	controller.SetTileHeight(500);
+	controller.SetNumberOfFirstIterations(6);
+
+	// Specific parameters
+	lsrm::BaatzParam params;
+	params.m_SpectralWeight = 0.7;
+	params.m_ShapeWeight = 0.3;
+	controller.SetSpecificParameters(params);
+	controller.SetThreshold(60*60);
+	
+	controller.RunSegmentation();
+	
+    return 0;
+}
diff --git a/CMakeLists.txt b/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f03f64c8bb3f5f280cfe56477c453cc122ffd3d3
--- /dev/null
+++ b/CMakeLists.txt
@@ -0,0 +1,42 @@
+#=========================================================================
+
+#  Program: Large Scale Generic Region Merging Library (GRM)
+#  Language: C++
+#  author: Lassalle Pierre
+
+
+
+#  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
+
+#  See grmlib-copyright.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.
+
+#=========================================================================
+project(LSGRM)
+
+cmake_minimum_required(VERSION 2.8)
+
+find_package(OTB)
+IF(OTB_FOUND)
+  include(${OTB_USE_FILE})
+ELSE(OTB_FOUND)
+  message(FATAL_ERROR
+    "Cannot build OTB project without OTB. Please set OTB_DIR.")
+ENDIF(OTB_FOUND)
+
+#Activate c++11
+include(CheckCXXCompilerFlag)
+CHECK_CXX_COMPILER_FLAG("-std=c++11" COMPILER_SUPPORTS_CXX11)
+if(COMPILER_SUPPORTS_CXX11)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -fpermissive -Wall -Wmaybe-uninitialized")
+else()
+        message(STATUS "The compiler ${CMAKE_CXX_COMPILER} has no C++11 support. Please use a different C++ compiler.")
+endif()
+
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/Code)
+
+add_subdirectory(Code)
+add_subdirectory(Applications)
diff --git a/Code/CMakeLists.txt b/Code/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..00ad3180e6d3af1d6eeb07e6e72076d89d5517bb
--- /dev/null
+++ b/Code/CMakeLists.txt
@@ -0,0 +1,22 @@
+#=========================================================================
+
+#  Program: Large Scale Generic Region Merging Library (GRM)
+#  Language: C++
+#  author: Lassalle Pierre
+
+
+
+#  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
+
+#  See grmlib-copyright.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.
+
+#=========================================================================add_library(OTBGRM 
+add_library(OTBLSGRM 
+			lsrmNeighborhood.cpp
+			lpContour.cpp)
+
+target_link_libraries(OTBLSGRM ${OTB_LIBRARIES})
diff --git a/Code/lpContour.cpp b/Code/lpContour.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f8f65b30c884d450558dff7ecd77b276b48365be
--- /dev/null
+++ b/Code/lpContour.cpp
@@ -0,0 +1,377 @@
+#include "lpContour.h"
+
+namespace lp
+{
+	
+	void ContourOperations::MergeContour(Contour& mergedContour,
+										 BoundingBox& mergedBBox,
+										 Contour& contour1,
+										 Contour& contour2,
+										 BoundingBox& bbox1,
+										 BoundingBox& bbox2,
+										 const CellIndex cid1,
+										 const CellIndex cid2,
+										 const std::size_t gridSizeX)
+	{
+		// First step consists of building the bounding box resulting from the fusion
+		// of the bounding boxes bbox1 and bbox2
+		mergedBBox = MergeBoundingBoxes(bbox1, bbox2);
+
+		// Create the associate matrix indicating where the cells are located
+		// inside the bbox
+		std::vector<bool> cellMatrix(mergedBBox.m_W*mergedBBox.m_H, false);
+
+		// Fill the cell matrix with the cells from both contours
+		{
+			CellLists borderCells;
+
+			// Fill with the cells of contour 1
+			GenerateBorderCells(borderCells, contour1, cid1, gridSizeX);
+			for(auto& cell: borderCells)
+				cellMatrix[GridToBBox(cell, mergedBBox, gridSizeX)] = true;
+			borderCells.clear();
+
+
+			// Fill with the cells of contour 2
+			GenerateBorderCells(borderCells, contour2, cid2, gridSizeX);
+			for(auto& cell: borderCells)
+				cellMatrix[GridToBBox(cell, mergedBBox, gridSizeX)] = true;
+		}
+
+		// Create the new contour
+		CreateNewContour(mergedContour, GridToBBox(cid1, mergedBBox, gridSizeX), cellMatrix, mergedBBox.m_W, mergedBBox.m_H);
+	}
+
+	void ContourOperations::CreateNewContour(Contour& newContour,
+											 const CellIndex cidx,
+											 const std::vector<bool>& cellMatrix,
+											 const std::size_t bboxWidth,
+											 const std::size_t bboxHeight)
+	{
+		// The first move is always 1
+		Push1(newContour);
+
+		// Previous move
+		short prevMove = 1;
+
+		// Local pixel id
+		CellIndex currIdx = cidx;
+
+		// Table containing id neighbors
+		long int neighbors[8];
+
+		for(;;)
+		{
+		
+			// Compute neighbor'ids
+			EIGHTNeighborhood(neighbors, currIdx, bboxWidth, bboxHeight);
+			
+			if(prevMove == 1)
+			{
+				if(neighbors[1] != -1 && cellMatrix[neighbors[1]])
+				{
+					Push0(newContour);
+					currIdx = currIdx + 1 - bboxWidth;
+					prevMove = 0;
+				}
+				else if(neighbors[2] != -1 && cellMatrix[neighbors[2]])
+				{
+					Push1(newContour);
+					currIdx++;
+					prevMove = 1;
+				}
+				else
+				{
+					Push2(newContour);
+					prevMove = 2;
+				}
+			}
+			else if(prevMove == 2)
+			{
+				if(neighbors[3] != -1 && cellMatrix[neighbors[3]])
+				{
+					Push1(newContour);
+					currIdx = currIdx + bboxWidth + 1;
+					prevMove = 1;
+				}
+				else if(neighbors[4] != -1 && cellMatrix[neighbors[4]])
+				{
+					Push2(newContour);
+					currIdx += bboxWidth;
+					prevMove = 2;
+				}
+				else
+				{
+					Push3(newContour);
+					prevMove = 3;
+				}
+			}
+			else if(prevMove == 3)
+			{
+			
+				if(neighbors[5] != -1 && cellMatrix[neighbors[5]])
+				{
+					Push2(newContour);
+					currIdx = currIdx - 1 + bboxWidth;
+					prevMove = 2;
+				}
+				else if(neighbors[6] != -1 && cellMatrix[neighbors[6]])
+				{
+					Push3(newContour);
+					currIdx -= 1;
+					prevMove = 3;
+				}
+				else
+				{
+					Push0(newContour);
+					prevMove = 0;
+				}
+			}
+			else
+			{
+				assert(prevMove == 0);
+
+				if(neighbors[7] != -1 && cellMatrix[neighbors[7]])
+				{
+					Push3(newContour);
+					currIdx = currIdx - bboxWidth - 1;
+					prevMove = 3;
+				}
+				else if(neighbors[0] != -1 && cellMatrix[neighbors[0]])
+				{
+					Push0(newContour);
+					currIdx -= bboxWidth;
+					prevMove = 0;
+				}
+				else
+				{
+					if(currIdx == cidx)
+						break;
+					else
+					{
+						Push1(newContour);
+						prevMove = 1;
+					}
+				}
+			}
+		}
+	}
+	
+	void ContourOperations::GenerateBorderCells(CellLists& borderCells, 
+												const Contour& contour, 
+												const CellIndex startCellId, 
+												const std::size_t gridSizeX)
+	{
+		// Add the first pixel to the border list
+		borderCells.insert(startCellId);
+
+		if(contour.size() > 8)
+		{
+	
+			// Intialize the first move at prev
+			short curr, prev = GetMove10(GetMove2(0, contour));
+
+			// Declare the current pixel index
+			CellIndex idx = startCellId;
+
+			// Explore the contour
+			for(ContourIndex cidx = 1; cidx < contour.size() / 2; cidx++)
+			{
+				curr = GetMove10(GetMove2(cidx, contour));
+				assert(curr >= 0 && curr < 4);
+
+				if(curr == 0)
+				{
+					// Impossible case is prev = 2;
+					assert(prev != 2);
+
+					//*
+					//*
+					if(prev == 0)
+					{
+						idx -= gridSizeX;
+						borderCells.insert(idx);
+					}
+			
+					//  *
+					// **
+					if(prev == 1)
+					{
+						idx = idx + 1 - gridSizeX;
+						borderCells.insert(idx);
+					}
+				}else if(curr == 1)
+				{
+					// Impossible case
+					assert(prev != 3);
+
+					// **
+					if(prev == 1)
+					{
+						idx++;
+						borderCells.insert(idx);
+					}
+
+					//*
+					//**
+					if (prev == 2)
+					{
+						idx = idx + 1 + gridSizeX;
+						borderCells.insert(idx);
+					}
+			
+				}else if(curr == 2)
+				{
+					// Impossible case
+					assert(prev != 0);
+
+					// *
+					// *
+					if(prev == 2)
+					{
+						idx += gridSizeX;
+						borderCells.insert(idx);
+					}
+
+					// **
+					// *
+					if(prev == 3)
+					{
+						idx = idx - 1 + gridSizeX;
+						borderCells.insert(idx);
+					}
+				}else
+				{
+					// Impossible case
+					assert(prev != 1);
+
+					if(prev == 0)
+					{
+						idx = idx - 1 - gridSizeX;
+						borderCells.insert(idx);
+					}
+
+					if(prev == 3)
+					{
+						idx--;
+						borderCells.insert(idx);
+					}
+				}
+
+				prev = curr;
+			}
+		}
+	}
+
+	CellIndex ContourOperations::BBoxToGrid(const CellIndex bboxId,
+											const BoundingBox& bbox,
+											const std::size_t gridSizeX)
+	{
+		CellIndex bbX = bboxId % bbox.m_W;
+		CellIndex bbY = bboxId / bbox.m_W;
+
+		CellIndex gridX = bbox.m_UX + bbX;
+		CellIndex gridY = bbox.m_UY + bbY;
+
+		CellIndex gridId = gridY * gridSizeX + gridX;
+
+		return gridId;
+	}
+
+	CellIndex ContourOperations::GridToBBox(const CellIndex gridId,
+											const BoundingBox& bbox,
+											const std::size_t gridSizeX)
+	{
+		CellIndex gridX = gridId % gridSizeX;
+		CellIndex gridY = gridId / gridSizeX;
+
+		CellIndex bbX = gridX - bbox.m_UX;
+		CellIndex bbY = gridY - bbox.m_UY;
+
+		CellIndex bbId = bbY * bbox.m_W + bbX;
+
+		return bbId;
+	}
+	
+	void ContourOperations::EIGHTNeighborhood(long int * neighborhood,
+											  const CellIndex id,
+											  const std::size_t width,
+											  const std::size_t height)
+	{
+		const unsigned int x = id % width;
+		const unsigned int y = id / width;
+
+		/* top */
+		neighborhood[0] = ( y > 0 ? (id - width) : -1 );
+
+		/* top right */
+		neighborhood[1] = ( (y > 0 && x < (width - 1) ) ? (id - width + 1) : -1 );
+		
+		/* right */
+		neighborhood[2] = ( x < (width - 1) ? (id + 1) : -1 );
+
+		/* bottom right */
+		neighborhood[3] = ( (x < (width - 1) && y < (height - 1) ) ? (id + 1 + width) : -1);
+
+		/* bottom */
+		neighborhood[4] = ( y < (height - 1) ? (id + width) : -1 );
+
+		/* bottom left */
+		neighborhood[5] = ( (y < (height - 1) && x > 0) ? (id + width - 1) : -1 );
+		
+		/* left */
+		neighborhood[6] = ( x > 0 ? (id - 1) : -1 );
+
+		/* top left */
+		neighborhood[7] = ( (x > 0 && y > 0) ? (id -width - 1) : - 1);	
+	}
+	
+	BoundingBox ContourOperations::MergeBoundingBoxes(const BoundingBox& bb1,
+													  const BoundingBox& bb2)
+	{
+		std::size_t min_ux, min_uy, max_xw, max_yh;
+		BoundingBox bb;
+
+		min_ux = std::min(bb1.m_UX, bb2.m_UX);
+		min_uy = std::min(bb1.m_UY, bb2.m_UY);
+		max_xw = std::max(bb1.m_UX + bb1.m_W, bb2.m_UX + bb2.m_W);
+		max_yh = std::max(bb1.m_UY + bb1.m_H, bb2.m_UY + bb2.m_H);
+
+		bb.m_UX = min_ux;
+		bb.m_UY = min_uy;
+		bb.m_W = max_xw - min_ux;
+		bb.m_H = max_yh - min_uy;
+
+		return bb;
+	}
+	
+	Move ContourOperations::GetMove2(ContourIndex idx, const Contour& contour)
+	{
+		return Move(2*contour[2*idx] + contour[idx*2 +1]);
+	}
+
+	short ContourOperations::GetMove10(const Move& m)
+	{
+		return (m[0]*2 + m[1]);
+	}
+	
+	void ContourOperations::Push0(Contour& contour)
+	{
+		contour.push_back(0); contour.push_back(0);
+	}
+
+	void ContourOperations::Push1(Contour& contour)
+	{
+		contour.push_back(1); contour.push_back(0);
+	}
+
+	void ContourOperations::Push2(Contour& contour)
+	{
+		contour.push_back(0); contour.push_back(1);
+	}
+
+	void ContourOperations::Push3(Contour& contour)
+	{
+		contour.push_back(1); contour.push_back(1);
+	}
+} // end of namespace lp
diff --git a/Code/lpContour.h b/Code/lpContour.h
new file mode 100644
index 0000000000000000000000000000000000000000..d5ac0dbab9ddf4bce6cdbe5c94f3dfcf744e160b
--- /dev/null
+++ b/Code/lpContour.h
@@ -0,0 +1,97 @@
+#ifndef __LP_CONTOUR_H
+#define __LP_CONTOUR_H
+#include <cassert>
+#include <iostream>
+#include <bitset>
+#include <vector>
+#include <utility>
+#include <unordered_set>
+#include <unordered_map>
+#include <boost/dynamic_bitset.hpp>
+
+namespace lp
+{
+	/* Move along a contour: 4 possibles: top (0), right(1), bottom (2), left(3). */
+	using Move = std::bitset<2>;
+	/* A contour is a list of moves represented by a pair of bits*/
+	using Contour = boost::dynamic_bitset<>;
+ 
+	using ContourIndex = boost::dynamic_bitset<>::size_type;
+
+	/* Index of a cell in the grid */
+	using CellIndex = std::size_t;
+
+	/* List of cell indices */
+	using CellLists = std::unordered_set<CellIndex>;
+	
+	struct BoundingBox
+	{
+		std::size_t m_UX;
+		std::size_t m_UY;
+		std::size_t m_W;
+		std::size_t m_H;
+	};
+
+	class ContourOperations
+	{
+	public:
+		/* Methods to fill a contour. */
+		static void Push0(Contour& contour); // Push a move to the top.
+		static void Push1(Contour& contour); // Push a move to the right.
+		static void Push2(Contour& contour); // Push a move to the bottom.
+		static void Push3(Contour& contour); // Push a move to the left.
+		
+		/* Methods to access to elements of a contour. */
+		static Move GetMove2(ContourIndex idx, const Contour& contour); // Retrieve a move
+		static short GetMove10(const Move& m); // Transform the move into a 10 base number (0,1,2 or 3).
+
+		/* Main methods */
+		static void MergeContour(Contour& mergedContour,
+								 BoundingBox& mergedBoundingBox,
+								 Contour& contour1,
+								 Contour& contour2,
+								 BoundingBox& bbox1,
+								 BoundingBox& bbox2,
+								 const CellIndex cid1,
+								 const CellIndex cid2,
+								 const std::size_t gridSizeX);
+		
+		static void GenerateBorderCells(CellLists& borderCells, // From a contour, the smallest id of the cell
+										const Contour& contour, // belonging to the shape and the size of the grid,
+										const CellIndex startCellId, // this method returns the ids of the border cells
+										const std::size_t gridSizeX); // of the shape.
+
+		static void CreateNewContour(Contour& newContour,
+									 const CellIndex startCellId,
+									 const std::vector<bool>& cellMatrix,
+									 const std::size_t bboxWidth,
+									 const std::size_t bboxHeight);
+
+		/* Utilities methods */
+		static BoundingBox MergeBoundingBoxes(const BoundingBox& bb1, // Given 2 bounding boxes return the bounding
+											  const BoundingBox& bb2); // which is the union.
+
+		static void EIGHTNeighborhood(long int * neighborhood, // return 8-neighbors of cell id.
+									  const CellIndex id,
+									  const std::size_t width,
+									  const std::size_t height);
+
+		static CellIndex BBoxToGrid(const CellIndex bboxId,
+									const BoundingBox& bbox,
+									const std::size_t gridSizeX);
+
+		static CellIndex GridToBBox(const CellIndex gridId, // Return the coordinates of the cell
+									const BoundingBox& bbox, // in the bounding box reference.
+									const std::size_t gridSizeX);
+	};
+} // end of namespace lp
+#endif
+
+
+
+
+
+
+
+
+
diff --git a/Code/lsgrmController.h b/Code/lsgrmController.h
new file mode 100644
index 0000000000000000000000000000000000000000..2fab165a8cb439ed583ac8cfb76c317fd71f65bb
--- /dev/null
+++ b/Code/lsgrmController.h
@@ -0,0 +1,74 @@
+#ifndef __LSGRM_CONTROLLER_H
+#define __LSGRM_CONTROLLER_H
+#include "lsrmGetInternalMemory.h"
+#include "lsgrmSplitter.h"
+#include "lsgrmGraphOperations.h"
+
+namespace lsgrm
+{	
+	template<class TSegmenter>
+	class Controller
+	{
+	public:
+
+		/* Some convenient typedefs */
+		using SegmenterType = TSegmenter;
+		using ImageType = typename SegmenterType::ImageType;
+		using SegmentationParameterType = typename SegmenterType::ParameterType;
+
+		/* Default constructor and destructor. */
+		Controller();
+		~Controller();
+
+		void RunSegmentation();
+
+		void SetImageDivision(bool f);
+		void SetInputImage(const std::string& str);
+		void SetTileDirectory(const std::string& str);
+		void SetTemporaryDirectory(const std::string& str);
+		void SetOutputGraphDirectory(const std::string& str);
+		void SetTileWidth(const unsigned int v);
+		void SetTileHeight(const unsigned int v);
+		void SetNumberOfFirstIterations(const unsigned int v);
+		void SetSpecificParameters(const SegmentationParameterType& params);
+		void SetThreshold(const float& t);
+		void SetInternalMemoryAvailable(long long unsigned int v); // expecting a value in Mbytes.
+		
+	private:
+
+		void GetAutomaticConfiguration();
+		void RetrieveProblemConfiguration();
+		
+
+		/* Parameters given by the user */
+		long long unsigned int m_Memory; // RAM available for the computation.
+		std::string m_InputImage; // Input image path (if the image needs to be divided).
+		std::string m_TileDirectory; // Directory containing tiles (if the process of dividing is not activated
+		                             // then the directory must contain the tiles and the info.txt file.)
+		std::string m_TemporaryDirectory; // Directory used to store intermediate files during the process.
+		std::string m_OutputGraphDirectory;
+		std::string m_OutputLabelImage; // Output label image path.
+
+		bool m_ImageDivisionActivated; // The input image must be divided.
+		bool m_CleanTemporaryDirectory; // Clean the temporary directory.
+
+		/* Specific segmentation parameters */
+		SegmentationParameterType m_SpecificParameters;
+		float m_Threshold;
+
+		/* Internal attribute members.*/
+		unsigned int m_ImageWidth;
+		unsigned int m_ImageHeight;
+		unsigned int m_ImageBands;
+		unsigned int m_NbTilesX;
+		unsigned int m_NbTilesY;
+		unsigned int m_NumberOfFirstIterations;
+		unsigned int m_Margin;
+		unsigned int m_TileWidth;
+		unsigned int m_TileHeight;
+		std::vector<ProcessingTile> m_Tiles;
+	};
+} // end of namespace lsgrm
+
+#include "lsgrmController.txx"
+#endif
diff --git a/Code/lsgrmController.txx b/Code/lsgrmController.txx
new file mode 100644
index 0000000000000000000000000000000000000000..1660abdeb7aac5991216119afd519b034330a558
--- /dev/null
+++ b/Code/lsgrmController.txx
@@ -0,0 +1,358 @@
+#ifndef __LSGRM_CONTROLLER_TXX
+#define __LSGRM_CONTROLLER_TXX
+
+namespace lsgrm
+{
+
+	template<class TSegmenter>
+	Controller<TSegmenter>::Controller()
+	{
+		m_Memory = 0;
+		m_ImageDivisionActivated = true;
+	}
+
+	template<class TSegmenter>
+	Controller<TSegmenter>::~Controller()
+	{
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::RunSegmentation()
+	{
+		// Rajouter un if pour vérifier si l'utilisateur a enclenché la procédure automatique
+		if(m_Memory < 1)
+		{
+			this->GetAutomaticConfiguration();
+			std::cout << m_Memory  << " bytes, tile dimension " << m_TileWidth << " X " << m_TileHeight
+					  << ", margin" << m_Margin  << " niter " << m_NumberOfFirstIterations << std::endl;
+		}
+		
+		// Divide the input image if necessary
+		if(m_ImageDivisionActivated)
+			SplitOTBImage<ImageType>(m_InputImage, m_TileDirectory, m_TileWidth,
+									 m_TileHeight, m_Margin, m_NumberOfFirstIterations);
+		
+		// Retrieve the problem configuration
+		RetrieveProblemConfiguration();
+		
+		// Print values
+		std::cout << m_Memory  << " bytes, tile dimension " << m_TileWidth << " X " << m_TileHeight
+				  << ", margin" << m_Margin  << " niter " << m_NumberOfFirstIterations << std::endl;
+		
+		// Boolean indicating if there are remaining fusions
+		bool isFusion = false;
+		
+		// Run first partial segmentation
+		auto accumulatedMemory = RunFirstPartialSegmentation<TSegmenter>(m_SpecificParameters,
+																		 m_Threshold,
+																		 m_NumberOfFirstIterations,
+																		 3,
+																		 m_Tiles,
+																		 m_TileDirectory,
+																		 m_NbTilesX,
+																		 m_NbTilesY,
+																		 m_Margin,
+																		 m_TileWidth,
+																		 m_TileHeight,
+																		 m_ImageWidth,
+																		 m_ImageHeight,
+																		 m_TemporaryDirectory,
+																		 isFusion);
+		
+		std::cout << "Accumulated memory " << accumulatedMemory << " bytes, there is fusion "<< isFusion << std::endl;
+
+		while(accumulatedMemory > m_Memory && isFusion)
+		{
+		    isFusion = false;
+			accumulatedMemory = RunPartialSegmentation<TSegmenter>(m_SpecificParameters,
+																   m_Threshold,
+																   3,
+																   m_Tiles,
+																   m_TemporaryDirectory,
+																   m_NbTilesX,
+																   m_NbTilesY,
+																   m_TileWidth,
+																   m_TileHeight,
+																   m_ImageWidth,
+																   m_ImageHeight,
+																   m_ImageBands,
+																   isFusion);
+			std::cout << "Accumulated memory " << accumulatedMemory << " bytes, there is fusion "
+					  << isFusion << std::endl;
+		}
+
+		if(accumulatedMemory <= m_Memory)
+		{	
+			// Merge all the graphs
+			MergeAllGraphsAndAchieveSegmentation<TSegmenter>(m_SpecificParameters,
+															 m_Threshold,
+															 m_Tiles,
+															 m_TemporaryDirectory,
+															 m_NbTilesX,
+															 m_NbTilesY,
+															 m_TileWidth,
+															 m_TileHeight,
+															 m_ImageWidth,
+															 m_ImageHeight,
+															 m_ImageBands,
+															 isFusion,
+															 m_OutputGraphDirectory);
+		}
+		else
+		{
+			// That means there are no more possible fusions but we can not store the ouput graph
+			// Todo do not clean up temporary directory before copying resulting graph to the output directory
+			// In the output directory add an info file to give the number of tiles.
+		}
+	}
+	
+	template<class TSegmenter>
+	void Controller<TSegmenter>::RetrieveProblemConfiguration()
+	{
+		// Open the lsgrm info file
+		std::ifstream in(m_TileDirectory + "info.txt");
+		assert(in.good());
+
+	
+		std::string line;
+		std::vector<std::string> tokens;
+	
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_ImageWidth = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_ImageHeight = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_ImageBands = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_NbTilesX = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_NbTilesY = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_TileWidth = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_TileHeight = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_Margin = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		std::getline(in, line);
+		boost::split(tokens, line, boost::is_any_of(":"));
+		m_NumberOfFirstIterations = static_cast<unsigned int>(atoi(tokens[1].c_str()));
+
+		m_Tiles.assign(m_NbTilesX * m_NbTilesY, ProcessingTile());
+		std::vector<std::string> subtokens;
+		unsigned int i = 0;
+		for(unsigned int row = 0; row < m_NbTilesY; ++row)
+		{
+			for(unsigned int col = 0; col < m_NbTilesX; ++col)
+			{
+				std::getline(in, line);
+				boost::split(tokens, line, boost::is_any_of(":"));
+				boost::split(subtokens, tokens[1], boost::is_any_of(","));
+
+				/* Margin at the top ? */
+				if( row > 0 ) 
+				{
+					m_Tiles[i].rows[0] = row * m_TileHeight;
+					m_Tiles[i].tileNeighbors[0] = i - m_NbTilesX;
+					m_Tiles[i].margin[0] = true;
+				}else
+				{
+					m_Tiles[i].rows[0] = 0;
+					m_Tiles[i].tileNeighbors[0] = -1;
+					m_Tiles[i].margin[0] = false;
+				}
+
+				/* Margin at the right ? */
+				if( col < m_NbTilesX - 1 )
+				{
+					m_Tiles[i].columns[1] = col * m_TileWidth + static_cast<unsigned int>(atoi(subtokens[2].c_str())) - 1;
+					m_Tiles[i].tileNeighbors[2] = i+1;
+					m_Tiles[i].margin[1] = true;
+				}
+				else
+				{
+					m_Tiles[i].columns[1] = m_ImageWidth - 1;
+					m_Tiles[i].tileNeighbors[2] = -1;
+					m_Tiles[i].margin[1] = false;
+				}
+
+				/* Margin at the bottom */
+				if( row < m_NbTilesY - 1)
+				{
+					m_Tiles[i].rows[1] = row * m_TileHeight + static_cast<unsigned int>(atoi(subtokens[3].c_str())) - 1;
+					m_Tiles[i].tileNeighbors[4] = i + m_NbTilesX;
+					m_Tiles[i].margin[2] = true;
+				}
+				else
+				{
+					m_Tiles[i].rows[1] = m_ImageHeight - 1;
+					m_Tiles[i].tileNeighbors[4] = -1;
+					m_Tiles[i].margin[2] = false;
+				}
+
+				/* Margin at the left */
+				if( col > 0 )
+				{
+					m_Tiles[i].columns[0] = col * m_TileWidth;
+					m_Tiles[i].tileNeighbors[6] = i-1;
+					m_Tiles[i].margin[3] = true;	
+				}
+				else
+				{
+					m_Tiles[i].columns[0] = 0;
+					m_Tiles[i].tileNeighbors[6] = -1;
+					m_Tiles[i].margin[3] = false;
+				}
+
+				/* Is there a neighbor at the rop right */
+				if(row > 0 && col < m_NbTilesX - 1)
+					m_Tiles[i].tileNeighbors[1] = i - m_NbTilesX + 1;
+				else
+					m_Tiles[i].tileNeighbors[1] = -1;
+
+				/* Is there a neighbor at the bottom right */
+				if(col < m_NbTilesX - 1 && row < m_NbTilesY - 1)
+					m_Tiles[i].tileNeighbors[3] = i + m_NbTilesX + 1;
+				else
+					m_Tiles[i].tileNeighbors[3] = -1;
+
+				/* Is there a neighbor at the bottom left */
+				if(row < m_NbTilesY - 1 && col > 0)
+					m_Tiles[i].tileNeighbors[5] = i + m_NbTilesX - 1;
+				else
+					m_Tiles[i].tileNeighbors[5] = -1;
+
+				/* Is there a neighbor at the top left */
+				if(col > 0 && row > 0)
+					m_Tiles[i].tileNeighbors[7] = i - m_NbTilesX - 1;
+				else
+					m_Tiles[i].tileNeighbors[7] = -1;
+			
+				i++;
+			}
+		}
+
+		in.close();
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::GetAutomaticConfiguration()
+	{
+		m_Memory = getMemorySize();
+		assert(m_Memory > 0);
+
+		m_Memory /= 2; // For safety and can prevent out of memory troubles
+
+		// Compute the size of an initial segment
+		using NodeType = typename TSegmenter::NodeType;
+		using NodePointer = typename TSegmenter::NodePointerType;
+		using EdgeType = typename TSegmenter::EdgeType;
+		
+		long long unsigned int sizePerNode = sizeof(NodePointer) + sizeof(NodeType) + 1 + 4 *(sizeof(EdgeType) + sizeof(float)); // last term is specific to BS.
+		long unsigned int maximumNumberOfNodes = std::ceil(m_Memory / sizePerNode);
+		unsigned int tileDimension = std::sqrt(maximumNumberOfNodes);
+
+		// Compute the stability margin. The naive strategy consider a margin value and a stable size equal.
+		unsigned int niter = 1;
+		unsigned int maxMargin = tileDimension/2;
+		unsigned int currMargin = static_cast<unsigned int>(pow(2, niter + 1) - 2);
+		unsigned int prevMargin = currMargin;
+
+		while(currMargin < maxMargin)
+		{
+			prevMargin = currMargin;
+			niter++;
+			currMargin = static_cast<unsigned int>(pow(2, niter + 1) - 2);
+		}
+
+		m_TileWidth = tileDimension - prevMargin;
+		m_TileHeight = m_TileWidth;
+		m_Margin = prevMargin;
+		m_NumberOfFirstIterations = niter - 1;
+	}
+
+	template <class TSegmenter>
+	void Controller<TSegmenter>::SetInternalMemoryAvailable(long long unsigned int v) // expecting a value in Mbytes.
+	{
+		assert(v > 0);
+		m_Memory = v * 1024ul * 1024ul;
+	}
+		
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetImageDivision(bool f)
+	{
+		m_ImageDivisionActivated = f;
+	}
+		
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetInputImage(const std::string& str)
+	{
+		m_InputImage = str;
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetOutputGraphDirectory(const std::string& str)
+	{
+		m_OutputGraphDirectory = str;
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetTileDirectory(const std::string& str)
+	{
+		m_TileDirectory = str;
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetTemporaryDirectory(const std::string& str)
+	{
+		m_TemporaryDirectory = str;
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetTileWidth(const unsigned int v)
+	{
+		m_TileWidth = v;
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetTileHeight(const unsigned int v)
+	{
+		m_TileHeight = v;
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetNumberOfFirstIterations(const unsigned int v)
+	{
+		m_NumberOfFirstIterations = v;
+		m_Margin = static_cast<unsigned int>(pow(2, m_NumberOfFirstIterations + 1) - 2);// 2^{n+1}-2
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetSpecificParameters(const SegmentationParameterType& params)
+	{
+		m_SpecificParameters = params;
+	}
+
+	template<class TSegmenter>
+	void Controller<TSegmenter>::SetThreshold(const float& t)
+	{
+		m_Threshold = t;
+	}
+} // end of namespace lsgrm
+
+#endif
diff --git a/Code/lsgrmGraphOperations.h b/Code/lsgrmGraphOperations.h
new file mode 100644
index 0000000000000000000000000000000000000000..d3f6d8c80fa819e0956e2c7a2479ab276033b2b6
--- /dev/null
+++ b/Code/lsgrmGraphOperations.h
@@ -0,0 +1,168 @@
+#ifndef __LSGRM_GRAPH_OPERATIONS_H
+#define __LSGRM_GRAPH_OPERATIONS_H
+#include "lsgrmHeader.h"
+#include "lsrmGraphOperations.h"
+#include "otbImageFileReader.h"
+
+namespace lsgrm
+{
+	struct ProcessingTile
+	{
+		long int rows[2]; // lower and upper rows (-1 means that the row has not be considered)
+		long int columns[2]; // lower and upper columns (-1 means that the row has not be considered)
+		long int tileNeighbors[8]; // tile Neighbors at (top, top right, right, bottom right, bottom, bottom left, left, top left)
+		bool margin[4]; // Is there a margin at top, left, bottom or right
+	};
+
+
+	template<class TSegmenter>
+	void MergeAllGraphsAndAchieveSegmentation(const typename TSegmenter::ParameterType& params,
+											  const float& threshold,
+											  std::vector<ProcessingTile>& tiles,
+											  const std::string& tmpDir,
+											  const unsigned int nbTilesX,
+											  const unsigned int nbTilesY,
+											  const unsigned int tileWidth,
+											  const unsigned int tileHeight,
+											  const unsigned int imageWidth,
+											  const unsigned int imageHeight,
+											  const unsigned int imageBands,
+											  bool& isFusion,
+											  const std::string& outputGraphDirectory);
+	
+	template<class TSegmenter>
+		long long unsigned int RunFirstPartialSegmentation(const typename TSegmenter::ParameterType& params,
+														   const float& threshold,
+														   const unsigned int niter,
+														   std::vector<ProcessingTile>& tiles,
+														   const std::string& tileDir,
+														   const unsigned int nbTilesX,
+														   const unsigned int nbTilesY,
+														   const unsigned int margin,
+														   const unsigned int tileWidth,
+														   const unsigned int tileHeight,
+														   const unsigned int imageWidth,
+														   const unsigned int imageHeight,
+														   const std::string& tmpDir,
+														   bool& isFusion);
+
+	template<class TSegmenter>
+		long long unsigned int RunPartialSegmentation(const typename TSegmenter::ParameterType& params,
+													  const float& threshold,
+													  const unsigned int niter,
+													  const unsigned int niter2,
+													  std::vector<ProcessingTile>& tiles,
+													  const std::string& tmpDir,
+													  const unsigned int nbTilesX,
+													  const unsigned int nbTilesY,
+													  const unsigned int tileWidth,
+													  const unsigned int tileHeight,
+													  const unsigned int imageWidth,
+													  const unsigned int imageHeight,
+													  const unsigned int imageBands,
+													  bool& isFusion);
+
+	template<class TSegmenter>
+	void RemoveUselessNodes(ProcessingTile& tile,
+							typename TSegmenter::GraphType& graph,
+							const unsigned int rowTile,
+							const unsigned int colTile,
+							const unsigned int nbTilesX,
+							const unsigned int nbTilesY,
+							const unsigned int imageWidth,
+							const unsigned int numberOfLayers);
+
+
+	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);
+
+	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);
+
+	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);
+
+	template<class TSegmenter>
+	void AddStabilityMargin(typename TSegmenter::GraphType& graph,
+							ProcessingTile& tile,
+							const std::string& tmpDir,
+							const unsigned int row,
+							const unsigned int col,
+							const unsigned int nbTilesX,
+							const unsigned int nbTilesY,
+							const unsigned int tileWidth,
+							const unsigned int tileHeight);
+
+	template<class TSegmenter>
+		void WriteStabilityMargin(std::unordered_map<
+								  typename TSegmenter::NodePointerType,
+								  unsigned int>& stabilityMargin,
+								  const std::string& nodesPath,
+								  const std::string& edgesPath);
+
+	template<class TSegmenter>
+		void ExtractStabilityMargin(std::unordered_map<typename TSegmenter::NodePointerType, unsigned int>& nodeMap,
+									const unsigned int 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);
+
+	template<class TSegmenter>
+		void DetectBorderNodes(typename TSegmenter::GraphType& graph,
+							   ProcessingTile& tile,
+							   std::unordered_map<typename TSegmenter::NodePointerType, unsigned int>& borderNodeMaps,
+							   const unsigned int imageWidth,
+							   const unsigned int imageHeight);
+
+	template<class TSegmenter>
+		void ReadGraph(typename TSegmenter::GraphType& graph,
+					   const std::string& nodesPath,
+					   const std::string& edgesPath);
+
+	template<class TSegmenter>
+		void WriteGraph(typename TSegmenter::GraphType& graph,
+						const std::string& tmpDir,
+						const unsigned int row,
+						const unsigned int col);
+
+	template<class TSegmenter>
+		long long unsigned int GetGraphMemory(typename TSegmenter::GraphType& graph);
+	
+	template<class TSegmenter>
+		void RemoveUnstableSegments(typename TSegmenter::GraphType& graph,
+									ProcessingTile& tile,
+									const unsigned int imageWidth);
+
+	template<class TSegmenter>
+		void RemoveEdgeToUnstableNode(typename TSegmenter::NodePointerType nodePtr);
+
+	template<class TSegmenter>
+		void RescaleGraph(typename TSegmenter::GraphType& graph,
+						  ProcessingTile& tile,
+						  const unsigned int rowTile,
+						  const unsigned int colTile,
+						  const unsigned int margin,
+						  const unsigned int tileWidth,
+						  const unsigned int tileHeight,
+						  const unsigned int imageWidth);
+}
+
+#include "lsgrmGraphOperations.txx"
+#endif
diff --git a/Code/lsgrmGraphOperations.txx b/Code/lsgrmGraphOperations.txx
new file mode 100644
index 0000000000000000000000000000000000000000..2e8085124e4a75cd5365ebebe0f9a53576a1544b
--- /dev/null
+++ b/Code/lsgrmGraphOperations.txx
@@ -0,0 +1,1359 @@
+#ifndef __LSGRM_GRAPH_OPERATIONS_TXX
+#define __LSGRM_GRAPH_OPERATIONS_TXX
+
+namespace lsgrm
+{
+	template<class TSegmenter>
+	void MergeAllGraphsAndAchieveSegmentation(const typename TSegmenter::ParameterType& params,
+											  const float& threshold,
+											  std::vector<ProcessingTile>& tiles,
+											  const std::string& tmpDir,
+											  const unsigned int nbTilesX,
+											  const unsigned int nbTilesY,
+											  const unsigned int tileWidth,
+											  const unsigned int tileHeight,
+											  const unsigned int imageWidth,
+											  const unsigned int imageHeight,
+											  const unsigned int imageBands,
+											  bool& isFusion,
+											  const std::string& outputGraphDirectory)
+	{
+		TSegmenter segmenter;
+		std::string nodesPath, edgesPath;
+		
+		std::cout << "Graph aggregation..." << std::endl;
+		for(unsigned int row = 0; row < nbTilesY; ++row)
+		{
+			for(unsigned int col = 0; col < nbTilesX; col++)
+			{
+				std::cout << "*" << std::flush;
+				typename TSegmenter::GraphType graph;
+
+				// Load the graph
+				nodesPath = tmpDir + "tile_nodes_" + std::to_string(row) + "_" + std::to_string(col) + ".bin";
+				edgesPath = tmpDir + "tile_edges_" + std::to_string(row) + "_" + std::to_string(col) + ".bin";
+				ReadGraph<TSegmenter>(graph, nodesPath, edgesPath);
+
+				segmenter.m_Graph.m_Nodes.insert(segmenter.m_Graph.m_Nodes.end(),
+												 graph.m_Nodes.begin(),
+												 graph.m_Nodes.end());
+			}
+		}
+		std::cout << "\nRemoving duplicated nodes and updating neighbors..." << std::endl;
+
+		for(unsigned int row = 0; row < nbTilesY; ++row)
+		{
+			for(unsigned int col = 0; col < nbTilesX; col++)
+			{
+				std::cout << "*" << std::flush;
+
+				std::unordered_map<long unsigned int,
+								   std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
+
+				BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], row, col,
+												nbTilesX, nbTilesY, borderPixelMap, imageWidth);
+
+					
+
+				
+				RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);
+
+					
+				UpdateNeighborsOfNoneDuplicatedNodes<TSegmenter>(borderPixelMap,
+																 imageWidth,
+																 imageHeight);
+			}
+		}
+
+		std::cout << "\nAchieve segmentation process..." << std::endl;
+		
+		// Segmentation of the graph
+		segmenter.SetImageWidth(imageWidth);
+		segmenter.SetImageHeight(imageHeight);
+		segmenter.SetNumberOfComponentsPerPixel(imageBands);
+		segmenter.SetParam(params);
+		segmenter.SetThreshold(threshold);
+		segmenter.SetDoBFSegmentation(true);
+		segmenter.SetNumberOfIterations(75);
+
+		lsrm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
+
+		// Write output graph to the output graph directory
+		WriteGraph<TSegmenter>(segmenter.m_Graph, outputGraphDirectory, 0, 0);
+				
+		typedef unsigned char ClusterPixelType;
+		typedef otb::VectorImage<ClusterPixelType, 2> ClusterImageType;
+		typedef otb::ImageFileWriter<ClusterImageType> ClusterImageWriterType;
+		auto clusterWriter = ClusterImageWriterType::New();
+		clusterWriter->SetFileName("out/finalimage.png");
+		clusterWriter->SetInput(segmenter.GetClusteredImageOutput());
+		clusterWriter->Update();
+	}
+	
+	template<class TSegmenter>
+	long long unsigned int RunPartialSegmentation(const typename TSegmenter::ParameterType& params,
+												  const float& threshold,
+												  const unsigned int niter,
+												  std::vector<ProcessingTile>& tiles,
+												  const std::string& tmpDir,
+												  const unsigned int nbTilesX,
+												  const unsigned int nbTilesY,
+												  const unsigned int tileWidth,
+												  const unsigned int tileHeight,
+												  const unsigned int imageWidth,
+												  const unsigned int imageHeight,
+												  const unsigned int imageBands,
+												  bool& isFusion)
+	{
+		long long unsigned int accumulatedMemory = 0;
+		std::string nodesPath, edgesPath;
+		isFusion = false;
+
+		const unsigned int numberOfNeighborLayers = static_cast<unsigned int>(pow(2, niter + 1) - 2);
+		std::cout << "Number of neighbor layers " << numberOfNeighborLayers << std::endl; 
+
+		std::cout << "Partial segmentations..." << std::endl;
+		for(unsigned int row = 0; row < nbTilesY; ++row)
+		{
+			for(unsigned int col = 0; col < nbTilesX; col++)
+			{
+				//std::cout << "*" << std::flush;
+				TSegmenter segmenter;
+				std::cout << "Tile " << row << " " << col << std::endl;
+				std::cout << "\tLoad graph..." << std::endl;
+				// Load the graph
+				nodesPath = tmpDir + "tile_nodes_" + std::to_string(row) + "_" + std::to_string(col) + ".bin";
+				edgesPath = tmpDir + "tile_edges_" + std::to_string(row) + "_" + std::to_string(col) + ".bin";
+				ReadGraph<TSegmenter>(segmenter.m_Graph, nodesPath, edgesPath);
+
+
+				// Add stability margin to the graph
+				{
+					std::cout << "\tAdd stability margin..." << std::endl;
+					AddStabilityMargin<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], tmpDir,
+												   row, col, nbTilesX, nbTilesY,
+												   tileWidth, tileHeight);
+
+					
+					std::unordered_map<long unsigned int,
+									   std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
+
+					std::cout << "\tBuild border pixel map..." << std::endl;
+					BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], 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>(tiles[row*nbTilesX + col], segmenter.m_Graph,
+												   row, col, nbTilesX, nbTilesY, imageWidth, numberOfNeighborLayers);
+
+				}
+
+
+				// Segmentation of the graph
+				segmenter.SetImageWidth(imageWidth);
+				segmenter.SetImageHeight(imageHeight);
+				segmenter.SetNumberOfComponentsPerPixel(imageBands);
+				segmenter.SetParam(params);
+				segmenter.SetThreshold(threshold);
+				segmenter.SetDoBFSegmentation(false);
+				segmenter.SetNumberOfIterations(niter);
+
+				std::cout << "\tPartial segmentation.." << std::endl;
+				auto merge = lsrm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
+
+				if(merge == true)
+					isFusion = true;
+
+				std::cout << "\tRemove unstable segments..." << std::endl;
+				// Remove unstable segments
+				RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], imageWidth);
+
+				// Retrieve the amount of memory to store this graph
+				accumulatedMemory += GetGraphMemory<TSegmenter>(segmenter.m_Graph);
+
+				std::cout << "\tWrite graph..." << std::endl;
+				// Write graph to temporay directory (warning specific to Baatz & Schape !!!)
+				WriteGraph<TSegmenter>(segmenter.m_Graph, tmpDir, row, col);
+				
+			}
+		}
+
+		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++)
+			{
+				std::cout << "*" << std::flush;
+				typename TSegmenter::GraphType graph;
+
+				// Load the graph
+				nodesPath = tmpDir + "tile_nodes_" + std::to_string(row) + "_" + std::to_string(col) + ".bin";
+				edgesPath = tmpDir + "tile_edges_" + std::to_string(row) + "_" + std::to_string(col) + ".bin";
+				ReadGraph<TSegmenter>(graph, nodesPath, edgesPath);
+
+				// Extract stability margin for all borders different from 0 imageWidth-1 et imageHeight -
+				// and write them to the stability margin
+				{
+					std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;
+
+					DetectBorderNodes<TSegmenter>(graph, tiles[row*nbTilesX + col],
+												  borderNodeMap, imageWidth, imageHeight);
+
+
+					ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);
+
+					std::string nodesPath = tmpDir + "tile_nodes_margin_" +
+						std::to_string(row) + "_" + std::to_string(col) + ".bin";
+					std::string edgesPath = tmpDir + "tile_edges_margin_" +
+						std::to_string(row) + "_" + std::to_string(col) + ".bin";
+
+							
+					WriteStabilityMargin<TSegmenter>(borderNodeMap, nodesPath, edgesPath);
+				}
+
+			}
+		}
+		std::cout << std::endl;
+
+		return accumulatedMemory;
+	}
+
+	template<class TSegmenter>
+	void RemoveUselessNodes(ProcessingTile& tile,
+							typename TSegmenter::GraphType& graph,
+							const unsigned int rowTile,
+							const unsigned int colTile,
+							const unsigned int nbTilesX,
+							const unsigned int nbTilesY,
+							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;
+		}
+
+		lsrm::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];
+			lsrm::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 = lsrm::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())
+									{
+										long int pixNeighborhood[4];
+										lsrm::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 = lsrm::GraphOperations<TSegmenter>::FindEdge(neighNit, *nit);
+							
+						assert(toNit != edg.GetRegion()->m_Edges.end());
+							
+						boundary = edg.m_Boundary;
+
+						neighNit->m_Edges.erase(toNit);
+
+						auto toRefNode = lsrm::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);
+					}
+				}
+			}
+		}
+
+		lsrm::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))
+					{
+						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 AddStabilityMargin(typename TSegmenter::GraphType& graph,
+							ProcessingTile& tile,
+							const std::string& tmpDir,
+							const unsigned int row,
+							const unsigned int col,
+							const unsigned int nbTilesX,
+							const unsigned int nbTilesY,
+							const unsigned int tileWidth,
+							const unsigned int tileHeight)
+	{
+		std::string nodesPath;
+		std::string edgesPath;
+		unsigned int rowNeigh;
+		unsigned int colNeigh;
+
+		// Margin to retrieve at top
+		if(row > 0)
+		{
+			colNeigh = col;
+			rowNeigh = row-1;
+
+			nodesPath = tmpDir + "tile_nodes_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+
+			edgesPath = tmpDir + "tile_edges_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+			
+			typename TSegmenter::GraphType subgraph;
+			ReadGraph<TSegmenter>(subgraph, nodesPath, edgesPath);
+			graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());
+		}
+
+		// Margin to retrieve at right
+		if(col < nbTilesX - 1)
+		{
+			colNeigh = col + 1;
+			rowNeigh = row;
+
+			nodesPath = tmpDir + "tile_nodes_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+
+			edgesPath = tmpDir + "tile_edges_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+			
+			typename TSegmenter::GraphType subgraph;
+			ReadGraph<TSegmenter>(subgraph, nodesPath, edgesPath);
+			
+			graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());
+
+		}
+
+		// Margin to retrieve at bottom
+		if(row < nbTilesY - 1)
+		{
+			colNeigh = col;
+			rowNeigh = row + 1;
+			
+			nodesPath = tmpDir + "tile_nodes_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+
+			edgesPath = tmpDir + "tile_edges_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+			
+			typename TSegmenter::GraphType subgraph;
+			ReadGraph<TSegmenter>(subgraph, nodesPath, edgesPath);
+			graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());
+			
+		}
+
+		// Margin to retrieve at left
+		if(col > 0)
+		{
+			colNeigh = col-1;
+			rowNeigh = row;
+
+			nodesPath = tmpDir + "tile_nodes_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+
+			edgesPath = tmpDir + "tile_edges_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+			
+			typename TSegmenter::GraphType subgraph;
+			ReadGraph<TSegmenter>(subgraph, nodesPath, edgesPath);
+			graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());
+		}
+
+		// Margin to retrieve at top right
+		if(row > 0 && col < nbTilesX - 1)
+		{
+			colNeigh = col+1;
+			rowNeigh = row-1;
+
+			nodesPath = tmpDir + "tile_nodes_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+
+			edgesPath = tmpDir + "tile_edges_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+			
+			typename TSegmenter::GraphType subgraph;
+			ReadGraph<TSegmenter>(subgraph, nodesPath, edgesPath);
+			graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());	
+		}
+
+		// Margin to retrieve at bottom right
+		if(row < nbTilesY - 1 && col < nbTilesX - 1)
+		{
+			colNeigh = col+1;
+			rowNeigh = row+1;
+
+			nodesPath = tmpDir + "tile_nodes_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+
+			edgesPath = tmpDir + "tile_edges_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+			
+			typename TSegmenter::GraphType subgraph;
+			ReadGraph<TSegmenter>(subgraph, nodesPath, edgesPath);
+			graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());
+		}
+
+		// Margin to retrieve at bottom left
+		if(row < nbTilesY - 1 && col > 0)
+		{
+			colNeigh = col-1;
+			rowNeigh = row+1;
+
+			nodesPath = tmpDir + "tile_nodes_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+
+			edgesPath = tmpDir + "tile_edges_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+			
+			typename TSegmenter::GraphType subgraph;
+			ReadGraph<TSegmenter>(subgraph, nodesPath, edgesPath);
+			graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());
+		}
+
+		// Margin to retrieve at top left
+		if(row > 0 && col > 0)
+		{
+			colNeigh = col-1;
+			rowNeigh = row-1;
+
+			nodesPath = tmpDir + "tile_nodes_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+
+			edgesPath = tmpDir + "tile_edges_margin_" +
+				std::to_string(rowNeigh) + "_" + std::to_string(colNeigh) + ".bin";
+			
+			typename TSegmenter::GraphType subgraph;
+			ReadGraph<TSegmenter>(subgraph, nodesPath, edgesPath);
+			graph.m_Nodes.insert(graph.m_Nodes.end(), subgraph.m_Nodes.begin(), subgraph.m_Nodes.end());	
+		}
+	}
+	
+	template<class TSegmenter>
+	long long unsigned int RunFirstPartialSegmentation(const typename TSegmenter::ParameterType& params,
+													   const float& threshold,
+													   const unsigned int niter,
+													   const unsigned int niter2,
+													   std::vector<ProcessingTile>& tiles,
+													   const std::string& tileDir,
+													   const unsigned int nbTilesX,
+													   const unsigned int nbTilesY,
+													   const unsigned int margin,
+													   const unsigned int tileWidth,
+													   const unsigned int tileHeight,
+													   const unsigned int imageWidth,
+													   const unsigned int imageHeight,
+													   const std::string& tmpDir,
+													   bool& isFusion)
+	{
+		using ImageType = typename TSegmenter::ImageType;
+		using Reader = otb::ImageFileReader<ImageType>;
+
+		long long unsigned int accumulatedMemory = 0;
+		isFusion = false;
+
+		const unsigned int numberOfNeighborLayers = static_cast<unsigned int>(pow(2, niter2 + 1) - 2);
+		std::cout << "Number of neighbor layers " << numberOfNeighborLayers << std::endl; 
+		
+		for(unsigned int row = 0; row < nbTilesY; ++row)
+		{
+			for(unsigned int col = 0; col < nbTilesX; col++)
+			{
+				// Reading image
+				auto tileReader = Reader::New();
+				tileReader->SetFileName(tileDir + "tile_" + std::to_string(row) + "_" + std::to_string(col) + ".tif");
+				tileReader->Update();
+
+				// Segmenting image
+				TSegmenter segmenter;
+				segmenter.SetParam(params);
+				segmenter.SetThreshold(threshold);
+				segmenter.SetDoBFSegmentation(false);
+				segmenter.SetNumberOfIterations(niter);
+				segmenter.SetInput(tileReader->GetOutput());
+				segmenter.Update();
+
+				if(segmenter.GetComplete() == false)
+					isFusion = true;
+
+				// Rescale the graph to be in the reference of the image
+				RescaleGraph<TSegmenter>(segmenter.m_Graph,
+										 tiles[row*nbTilesX + col],
+										 row,
+										 col,
+										 margin,
+										 tileWidth,
+										 tileHeight,
+										 imageWidth);
+
+				// Remove unstable segments
+				RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], imageWidth);
+
+				// Retrieve the amount of memory to store this graph
+				accumulatedMemory += GetGraphMemory<TSegmenter>(segmenter.m_Graph);
+
+				// Write graph to temporay directory (warning specific to Baatz & Schape !!!)
+				WriteGraph<TSegmenter>(segmenter.m_Graph, tmpDir, row, col);
+
+				// Extract stability margin for all borders different from 0 imageWidth-1 et imageHeight -1
+				// and write them to the stability margin
+				{
+					std::unordered_map<typename TSegmenter::NodePointerType, unsigned int> borderNodeMap;
+
+					DetectBorderNodes<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col],
+												  borderNodeMap, imageWidth, imageHeight);
+
+
+					ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);
+
+					std::string nodesPath = tmpDir + "tile_nodes_margin_" +
+						std::to_string(row) + "_" + std::to_string(col) + ".bin";
+					std::string edgesPath = tmpDir + "tile_edges_margin_" +
+						std::to_string(row) + "_" + std::to_string(col) + ".bin";
+
+							
+					WriteStabilityMargin<TSegmenter>(borderNodeMap, nodesPath, edgesPath);
+				}
+			}
+		}
+
+		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;
+				}
+			}
+		}
+	}
+
+	// Specific to Baatz & Schape (sorry it has to be fast)
+	template<class TSegmenter>
+	void ReadGraph(typename TSegmenter::GraphType& graph,
+				   const std::string& nodesPath,
+				   const std::string& edgesPath)
+	{
+		FILE * nodeStream = fopen(nodesPath.c_str(), "r");
+		assert(nodeStream != NULL);
+		std::unordered_map<long unsigned int, typename TSegmenter::NodePointerType> nodeMap;
+
+		// Read the size of the graph
+		{
+			std::size_t graphSize;
+			fread(&graphSize, sizeof(graphSize), 1, nodeStream);
+			graph.m_Nodes.reserve(graphSize);
+
+			// Creation of the nodes
+			for(std::size_t i = 0; i < graphSize; i++)
+			{
+				typename TSegmenter::NodePointerType node(new typename TSegmenter::NodeType());
+
+				node->m_Valid = true;
+				node->m_Expired = false;
+				node->m_IsMerged = true; // force to compute costs for the first iteration
+				
+				fread(&(node->m_Id), sizeof(node->m_Id), 1, nodeStream);
+				fread(&(node->m_Perimeter), sizeof(node->m_Perimeter), 1, nodeStream);
+				fread(&(node->m_Area), sizeof(node->m_Area), 1, nodeStream);
+				fread(&(node->m_Bbox.m_UX), sizeof(node->m_Bbox.m_UX), 1, nodeStream);
+				fread(&(node->m_Bbox.m_UY), sizeof(node->m_Bbox.m_UY), 1, nodeStream);
+				fread(&(node->m_Bbox.m_W), sizeof(node->m_Bbox.m_W), 1, nodeStream);
+				fread(&(node->m_Bbox.m_H), sizeof(node->m_Bbox.m_H), 1, nodeStream);
+
+				std::size_t bands;
+				fread(&(bands), sizeof(bands), 1, nodeStream);
+				node->m_Means.assign(bands, 0);
+				node->m_SquareMeans.assign(bands, 0);
+				node->m_SpectralSum.assign(bands, 0);
+				node->m_Std.assign(bands, 0);
+
+				for(unsigned int b = 0; b < bands; b++)
+				{
+					fread(&(node->m_Means[b]), sizeof(node->m_Means[b]), 1, nodeStream);
+					fread(&(node->m_SquareMeans[b]), sizeof(node->m_SquareMeans[b]), 1, nodeStream);
+					fread(&(node->m_SpectralSum[b]), sizeof(node->m_SpectralSum[b]), 1, nodeStream);
+					fread(&(node->m_Std[b]), sizeof(node->m_Std[b]), 1, nodeStream);
+				}
+
+				{
+					std::size_t contourSize;
+					fread(&(contourSize), sizeof(contourSize), 1, nodeStream);
+					short moves[contourSize];
+					fread(&moves, sizeof(short), contourSize, nodeStream);
+					for(unsigned int b = 0; b < contourSize; b++)
+					{
+						node->m_Contour.push_back(moves[b]);
+					}
+				}
+				nodeMap[node->m_Id] = node;
+				graph.m_Nodes.push_back(node);
+			}
+		}
+		fclose(nodeStream);
+
+		// Build edges with nodeMap
+		FILE * edgeStream = fopen(edgesPath.c_str(), "r");
+		assert(edgeStream != NULL);
+
+		for(auto& node: graph.m_Nodes)
+		{
+			long unsigned int nodeId;
+			fread(&(nodeId), sizeof(nodeId), 1, edgeStream);
+			assert(nodeId == node->m_Id);
+			
+			std::size_t edgeSize;
+			fread(&(edgeSize), sizeof(edgeSize), 1, edgeStream);
+			node->m_Edges.reserve(edgeSize);
+			for(unsigned int b = 0; b < edgeSize; b++)
+			{
+				long unsigned int targetId;
+				unsigned int boundary;
+				fread(&(targetId), sizeof(targetId), 1, edgeStream);
+				fread(&(boundary), sizeof(boundary), 1, edgeStream);
+				node->m_Edges.push_back(typename TSegmenter::EdgeType(nodeMap[targetId], 0, boundary));
+			}
+		}
+
+		fclose(edgeStream);
+	}
+
+	// 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 size of the graph
+		{
+			// Write number of nodes
+			std::size_t size = stabilityMargin.size();
+			fwrite(&size, sizeof(size), 1, nodeStream);
+		}
+
+		for(auto& kv : stabilityMargin)
+		{
+			auto node = kv.first;
+
+			fwrite(&(node->m_Id), sizeof(node->m_Id), 1, nodeStream);
+			fwrite(&(node->m_Perimeter), sizeof(node->m_Perimeter), 1, nodeStream);
+			fwrite(&(node->m_Area), sizeof(node->m_Area), 1, nodeStream);
+			fwrite(&(node->m_Bbox.m_UX), sizeof(node->m_Bbox.m_UX), 1, nodeStream);
+			fwrite(&(node->m_Bbox.m_UY), sizeof(node->m_Bbox.m_UY), 1, nodeStream);
+			fwrite(&(node->m_Bbox.m_W), sizeof(node->m_Bbox.m_W), 1, nodeStream);
+			fwrite(&(node->m_Bbox.m_H), sizeof(node->m_Bbox.m_H), 1, nodeStream);
+
+			std::size_t bands = node->m_Means.size();
+			fwrite(&(bands), sizeof(bands), 1, nodeStream);
+
+			for(unsigned int b = 0; b < node->m_Means.size(); b++)
+			{
+				fwrite(&(node->m_Means[b]), sizeof(node->m_Means[b]), 1, nodeStream);
+				fwrite(&(node->m_SquareMeans[b]), sizeof(node->m_SquareMeans[b]), 1, nodeStream);
+				fwrite(&(node->m_SpectralSum[b]), sizeof(node->m_SpectralSum[b]), 1, nodeStream);
+				fwrite(&(node->m_Std[b]), sizeof(node->m_Std[b]), 1, nodeStream);
+			}
+
+			{
+				std::size_t contourSize = node->m_Contour.size();
+				fwrite(&(contourSize), sizeof(contourSize), 1, nodeStream);
+				short moves[contourSize];
+				for(unsigned int b = 0; b < contourSize; b++)
+				{
+					moves[b] = node->m_Contour[b];
+				}
+				fwrite(moves, sizeof(short), contourSize, 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())
+				{
+					fwrite(&(edg.GetRegion()->m_Id), sizeof(edg.GetRegion()->m_Id), 1, edgeStream);
+					fwrite(&(edg.m_Boundary), sizeof(edg.m_Boundary), 1, edgeStream);
+				}
+			}
+		}
+
+		fclose(nodeStream);
+		fclose(edgeStream);
+	}
+
+	// Specific to Baatz & Schape (sorry it has to be fast)
+	template<class TSegmenter>
+	void WriteGraph(typename TSegmenter::GraphType& graph,
+					const std::string& tmpDir,
+					const unsigned int row,
+					const unsigned int col)
+	{
+		// One file for the nodes
+		std::string nodeFile = tmpDir + "tile_nodes_" +
+			std::to_string(row) + "_" + std::to_string(col) + ".bin";
+		// One file for the edges
+		std::string edgeFile = tmpDir + "tile_edges_" +
+			std::to_string(row) + "_" + std::to_string(col) + ".bin";
+
+	    FILE * nodeStream = fopen(nodeFile.c_str(), "wb");
+		assert(nodeStream != NULL);
+		FILE * edgeStream = fopen(edgeFile.c_str(), "wb");
+		assert(edgeStream != NULL);
+
+		// Write size of the graph
+		{
+			// Write number of nodes
+			std::size_t size = graph.m_Nodes.size();
+			fwrite(&size, sizeof(size), 1, nodeStream);
+		}
+
+
+		for(auto& node : graph.m_Nodes)
+		{
+			
+			fwrite(&(node->m_Id), sizeof(node->m_Id), 1, nodeStream);
+			fwrite(&(node->m_Perimeter), sizeof(node->m_Perimeter), 1, nodeStream);
+			fwrite(&(node->m_Area), sizeof(node->m_Area), 1, nodeStream);
+			fwrite(&(node->m_Bbox.m_UX), sizeof(node->m_Bbox.m_UX), 1, nodeStream);
+			fwrite(&(node->m_Bbox.m_UY), sizeof(node->m_Bbox.m_UY), 1, nodeStream);
+			fwrite(&(node->m_Bbox.m_W), sizeof(node->m_Bbox.m_W), 1, nodeStream);
+			fwrite(&(node->m_Bbox.m_H), sizeof(node->m_Bbox.m_H), 1, nodeStream);
+
+			std::size_t bands = node->m_Means.size();
+			fwrite(&(bands), sizeof(bands), 1, nodeStream);
+
+			for(unsigned int b = 0; b < node->m_Means.size(); b++)
+			{
+				fwrite(&(node->m_Means[b]), sizeof(node->m_Means[b]), 1, nodeStream);
+				fwrite(&(node->m_SquareMeans[b]), sizeof(node->m_SquareMeans[b]), 1, nodeStream);
+				fwrite(&(node->m_SpectralSum[b]), sizeof(node->m_SpectralSum[b]), 1, nodeStream);
+				fwrite(&(node->m_Std[b]), sizeof(node->m_Std[b]), 1, nodeStream);
+			}
+
+			{
+				std::size_t contourSize = node->m_Contour.size();
+				fwrite(&(contourSize), sizeof(contourSize), 1, nodeStream);
+				short moves[contourSize];
+				for(unsigned int b = 0; b < contourSize; b++)
+				{
+					moves[b] = node->m_Contour[b];
+				}
+				fwrite(moves, sizeof(short), contourSize, 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)
+			{
+				fwrite(&(edg.GetRegion()->m_Id), sizeof(edg.GetRegion()->m_Id), 1, edgeStream);
+				fwrite(&(edg.m_Boundary), sizeof(edg.m_Boundary), 1, edgeStream);
+			}
+		}
+
+		fclose(nodeStream);
+		fclose(edgeStream);
+	}
+
+	template<class TSegmenter>
+	long long unsigned int GetGraphMemory(typename TSegmenter::GraphType& graph)
+	{
+		long long unsigned int memory = 0;
+		long unsigned int numberOfMoves = 0;
+
+		// Amount of memory needed to store the nodes
+		memory += graph.m_Nodes.size() * (sizeof(*graph.m_Nodes.begin()) + sizeof(graph.m_Nodes.begin()) + 4 * sizeof(float));
+
+		for(auto& node : graph.m_Nodes)
+		{
+			memory += node->m_Edges.size() * (sizeof(*node->m_Edges.begin()) + sizeof(node->m_Edges.begin()));
+			numberOfMoves += node->m_Contour.size();
+		}
+		
+		memory += std::ceil(numberOfMoves / 4);
+
+		return memory;
+	}
+
+	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);
+				}
+			}
+		}
+
+		lsrm::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 = lsrm::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 margin,
+					  const unsigned int tileWidth,
+					  const unsigned int tileHeight,
+					  const unsigned int imageWidth)
+	{
+		unsigned int rowNodeTile, colNodeTile;
+		unsigned int rowNodeImg, colNodeImg;
+		unsigned int diffWidthMargin = 0;
+		unsigned int diffHeightMargin = 0;
+		unsigned int currTileWidth = tile.columns[1] - tile.columns[0] + 1;
+
+		if(tile.margin[0])
+			diffHeightMargin = margin; 
+		if(tile.margin[1])
+			currTileWidth += margin;
+		if(tile.margin[3])
+		{
+			diffWidthMargin = margin;
+			currTileWidth += margin;
+		}
+
+		unsigned int count = 0;
+		for(auto& node : graph.m_Nodes)
+		{
+			
+			// Row of the start pixel of the node
+			rowNodeTile = node->m_Id / currTileWidth;
+			// Col of the start pixel of the node
+			colNodeTile = node->m_Id % currTileWidth;
+			
+			// Row of the start pixel of the node in the image
+			rowNodeImg = rowTile * tileHeight + rowNodeTile - diffHeightMargin;
+			// Col of the start pixel of the node in the image
+			colNodeImg = colTile * tileWidth + colNodeTile - diffWidthMargin;
+			
+			// New id of the node in the image			
+			node->m_Id = rowNodeImg * imageWidth + colNodeImg;
+
+
+			// Change also its bounding box
+			node->m_Bbox.m_UX = colTile * tileWidth + node->m_Bbox.m_UX - diffWidthMargin;
+			node->m_Bbox.m_UY = rowTile * tileHeight + node->m_Bbox.m_UY - diffHeightMargin;
+
+			++count;
+		}
+	}
+}
+
+#endif
+
+
+// // To be removed, read the graph
+// std::string nodesPath = tmpDir + "tile_nodes_" +
+// 	std::to_string(row) + "_" + std::to_string(col) + ".bin";
+// std::string edgesPath = tmpDir + "tile_edges_" +
+// 	std::to_string(row) + "_" + std::to_string(col) + ".bin";
+
+// typename TSegmenter::GraphType graph;
+// ReadGraph<TSegmenter>(graph, nodesPath, edgesPath);
+// std::cout << graph.m_Nodes.size() << " vs " << segmenter.m_Graph.m_Nodes.size() << std::endl;
+
+// // compare graphs
+// auto rniter = segmenter.m_Graph.m_Nodes.begin();
+// for(auto& tnode : graph.m_Nodes)
+// {
+// 	auto rnode = *rniter;
+// 	assert(tnode->m_Id == rnode->m_Id);
+// 	assert(tnode->m_Perimeter == rnode->m_Perimeter);
+// 	assert(tnode->m_Area == rnode->m_Area);
+// 	assert(tnode->m_Bbox.m_UX == rnode->m_Bbox.m_UX);
+// 	assert(tnode->m_Bbox.m_UY == rnode->m_Bbox.m_UY);
+// 	assert(tnode->m_Bbox.m_W == rnode->m_Bbox.m_W);
+// 	assert(tnode->m_Bbox.m_H == rnode->m_Bbox.m_H);
+
+// 	for(unsigned int b = 0; b < rnode->m_Means.size(); b++)
+// 	{
+// 		assert(tnode->m_Means[b] == rnode->m_Means[b]);
+// 		assert(tnode->m_SquareMeans[b] == rnode->m_SquareMeans[b]);
+// 		assert(tnode->m_SpectralSum[b] == rnode->m_SpectralSum[b]);
+// 		assert(tnode->m_Std[b] == rnode->m_Std[b]);
+// 	}
+
+// 	for(unsigned int b = 0; b < rnode->m_Contour.size(); b++)
+// 	{
+// 		assert(tnode->m_Contour[b] == rnode->m_Contour[b]);
+// 	}
+
+// 	for(unsigned int b = 0; b < tnode->m_Edges.size(); b++)
+// 	{
+// 		assert(tnode->m_Edges[b].GetRegion()->m_Id == rnode->m_Edges[b].GetRegion()->m_Id);
+// 		assert(tnode->m_Edges[b].m_Boundary == rnode->m_Edges[b].m_Boundary);
+// 	}
+					
+// 	++rniter;
+// }
+
+
+// 	std::cout << "Check validity of neighborhood" << std::endl;
+				// using NPtr = typename TSegmenter::NodePointerType;
+				// using LUI = long unsigned int;
+
+				// std::unordered_map<LUI,NPtr> nodeMap;
+				// for(auto& n : segmenter.m_Graph.m_Nodes)
+				// {
+				// 	lp::CellLists borderPixels;
+				// 	lp::ContourOperations::GenerateBorderCells(borderPixels, n->m_Contour,
+				// 											   n->m_Id, imageWidth);
+
+				// 	for(auto& pix : borderPixels)
+				// 		nodeMap[pix] = n;
+				// }
+
+				// for(auto& n : segmenter.m_Graph.m_Nodes)
+				// {
+				// 	lp::CellLists borderPixels;
+				// 	lp::ContourOperations::GenerateBorderCells(borderPixels, n->m_Contour,
+				// 											   n->m_Id, imageWidth);
+
+				// 	for(auto& pix : borderPixels)
+				// 	{
+				// 		long int neighborhood[4];
+				// 		lsrm::FOURNeighborhood(neighborhood, pix, imageWidth, imageHeight);
+
+				// 		for(short j = 0; j < 4; j++)
+				// 		{
+				// 			if(neighborhood[j] > -1)
+				// 			{
+				// 				auto isPix = nodeMap.find(neighborhood[j]);
+				// 				if(isPix != nodeMap.end() && isPix->second != n)
+				// 				{
+										
+				// 					auto toN = lsrm::GraphOperations<TSegmenter>::FindEdge(isPix->second, n);
+				// 					if(toN == isPix->second->m_Edges.end())
+				// 					{
+				// 						std::cout << n->m_Id << " et " << isPix->second->m_Id << " pixel n " << pix
+				// 								  << " pixel neigh " << neighborhood[j] << std::endl;
+
+				// 						lp::CellLists pixels;
+				// 						lp::ContourOperations::GenerateAllCells(pixels, n->m_Contour, n->m_Bbox, n->m_Id, imageWidth);
+				// 						for(auto& pix : pixels)
+				// 						{
+				// 							std::cout << pix << " ";
+				// 						}
+				// 						std::cout  << std::endl;
+				// 						pixels.clear();
+				// 						lp::ContourOperations::GenerateAllCells(pixels, isPix->second->m_Contour, isPix->second->m_Bbox, isPix->second->m_Id,
+				// 																imageWidth);
+				// 						for(auto& pix : pixels)
+				// 						{
+				// 							std::cout << pix << " ";
+				// 						}
+				// 						std::cout  << std::endl;
+				// 					}
+				// 					assert(toN != isPix->second->m_Edges.end());
+				// 					auto toNeigh = lsrm::GraphOperations<TSegmenter>::FindEdge(n,isPix->second);
+				// 					assert(toNeigh != n->m_Edges.end());
+				// 				}
+				// 			}
+				// 		}
+				// 	}
+				// }
+
+				// typedef unsigned char ClusterPixelType;
+				// typedef otb::VectorImage<ClusterPixelType, 2> ClusterImageType;
+				// typedef otb::ImageFileWriter<ClusterImageType> ClusterImageWriterType;
+				// auto clusterWriter = ClusterImageWriterType::New();
+				// clusterWriter->SetFileName("img_outs/" + std::to_string(row) + "_" + std::to_string(col) + ".png");
+				// clusterWriter->SetInput(segmenter.GetClusteredImageOutput());
+				// clusterWriter->Update();
+
+
+
diff --git a/Code/lsgrmHeader.h b/Code/lsgrmHeader.h
new file mode 100644
index 0000000000000000000000000000000000000000..e52c5c88b00d51f916eb596d04dbabc09867ca00
--- /dev/null
+++ b/Code/lsgrmHeader.h
@@ -0,0 +1,14 @@
+#ifndef __LSGRM_HEADER_H
+#define __LSGRM_HEADER_H
+#include <cassert>
+#include <cstdlib>
+#include <string>
+#include <sstream>
+#include <fstream>
+#include <algorithm>
+#include <vector>
+#include <iterator>
+#include <stack>
+#include <boost/algorithm/string.hpp>
+
+#endif
diff --git a/Code/lsgrmSplitter.h b/Code/lsgrmSplitter.h
new file mode 100644
index 0000000000000000000000000000000000000000..7a86a07de068b2c040959147c0e1b729ac30f526
--- /dev/null
+++ b/Code/lsgrmSplitter.h
@@ -0,0 +1,27 @@
+#ifndef __LSGRM_SPLITTER_H
+#define __LSGRM_SPLITTER_H
+#include "lsgrmHeader.h"
+#include <otbVectorImage.h>
+#include <otbImageFileReader.h>
+#include <otbImageFileWriter.h>
+#include <otbMultiChannelExtractROI.h>
+
+namespace lsgrm
+{
+	// Split with the OTB library
+	template <class TInputImage>
+	void SplitOTBImage(const std::string& inputImagePath, // Path to the input image
+					   const std::string& ouputDirectoryPath, // Path to the directory storing the tiles
+					   unsigned int& tileWidth, // width of the tile
+					   unsigned int& tileHeight, // height of the tile
+					   const unsigned int margin, // stability margin
+					   const unsigned int niter
+		);
+
+	// Other libraries can be used to split the
+	// image.
+	
+} // end of namespace lsgrm
+
+#include "lsgrmSplitter.txx"
+#endif
diff --git a/Code/lsgrmSplitter.txx b/Code/lsgrmSplitter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..9441286f0603b2811ee76dbfb2e5b2f6d8f5e024
--- /dev/null
+++ b/Code/lsgrmSplitter.txx
@@ -0,0 +1,147 @@
+#ifndef __LSGRM_SPLITTER_TXX
+#define __LSGRM_SPLITTER_TXX
+
+namespace lsgrm
+{
+	template <class TInputImage>
+	void SplitOTBImage(const std::string& inputImagePath, // Path to the input image
+					   const std::string& ouputDirectoryPath, // Path to the directory storing the tiles
+					   unsigned int& tileWidth, // width of the tile
+					   unsigned int& tileHeight, // height of the tile
+					   const unsigned int margin, // stability margin
+					   const unsigned int niter)
+	{
+
+		/* Some convenient typedefs to short the code. */
+		using ImageType = TInputImage;
+		using InternalPixelType = typename ImageType::InternalPixelType;
+		using Extractor = otb::MultiChannelExtractROI<InternalPixelType, InternalPixelType>;
+		using Reader = otb::ImageFileReader<ImageType>;
+		using Writer = otb::ImageFileWriter<ImageType>;
+
+		/* Read the input image and update information about its dimension */
+		auto reader = Reader::New();
+		reader->SetFileName(inputImagePath);
+		reader->UpdateOutputInformation();
+		auto imagePtr = reader->GetOutput();
+		auto imageWidth = imagePtr->GetLargestPossibleRegion().GetSize()[0];
+		auto imageHeight = imagePtr->GetLargestPossibleRegion().GetSize()[1];
+		auto imageBands = imagePtr->GetNumberOfComponentsPerPixel();
+
+		if(tileWidth + 2*margin >= imageWidth || tileHeight + 2* margin >= imageHeight)
+		{
+			std::cout << "One dimension of the tile is greater than the image" << std::endl;
+			exit(1);
+		}
+
+		// Determine final size of the tiles according to the dimension of the image
+		unsigned int k = std::ceil(imageWidth / tileWidth);
+		while(imageWidth % k)
+			k++;
+		tileWidth = imageWidth / k;
+		assert(tileWidth > 0);
+		k = std::ceil(imageHeight / tileHeight);
+		while(imageWidth % k)
+			k++;
+		tileHeight = imageHeight / k;
+		assert(tileHeight > 0);
+			
+		
+		/* Determine the number of tiles according to the tile dimension given by the user. */
+		const unsigned int nbTilesX = imageWidth / tileWidth + ( imageWidth % tileWidth > 0 ? 1 : 0);
+		const unsigned int nbTilesY = imageHeight / tileHeight + ( imageHeight % tileHeight > 0 ? 1 : 0);
+
+		/* Open text file to record general parameters for the LSGRM library. */
+		std::ofstream lsgrmLog(ouputDirectoryPath + "info.txt");
+		assert(lsgrmLog.good());
+		lsgrmLog << "Image width:"<< imageWidth << "\n" 
+				 << "Image height:"<< imageHeight << "\n"
+				 << "Number of bands:"<< imageBands << "\n"
+				 << "Number of tiles according to X:" << nbTilesX << "\n"
+				 << "Number of tiles according to Y:" << nbTilesY << "\n"
+				 << "Tile width:" << tileWidth << "\n"
+				 << "Tile height:" << tileHeight << "\n"
+				 << "Stability Margin value:" << margin << "\n"
+				 << "Number of iterations: " << niter << "\n";
+
+		/* Local variables for the next loop */
+		unsigned int startX, startY; // Upper left coordinates of the tile.
+		unsigned int sizeX, sizeY; // Size of the tiles.
+		
+		/* Loop over the tiles*/
+		for(unsigned int row = 0; row < nbTilesY; ++row)
+		{
+			for(unsigned int col = 0; col < nbTilesX; ++col)
+			{
+
+				startX = col * tileWidth;
+				startY = row * tileHeight;
+				sizeX = std::min(tileWidth, static_cast<unsigned int>(imageWidth - startX));
+				sizeY = std::min(tileHeight, static_cast<unsigned int>(imageHeight - startY));
+
+				lsgrmLog << "Tile " << row*nbTilesY + col << ":"
+						 << startX << "," << startY << "," << sizeX << "," << sizeY << ",";
+
+				/* Margin at the top ? */
+				if( row > 0 ) 
+				{
+					startY -= margin;
+					sizeY += margin;
+					lsgrmLog << "1,";
+				}else
+					lsgrmLog << "0,";
+
+				/* Margin at the right */
+				if( col < nbTilesX - 1 )
+				{
+					sizeX += margin;
+					lsgrmLog << "1,";
+				}
+				else
+					lsgrmLog << "0,";
+
+				/* Margin at the bottom */
+				if( row < nbTilesY - 1)
+				{
+					sizeY += margin;
+					lsgrmLog << "1,";
+				}
+				else
+					lsgrmLog << "0,";
+
+				/* Margin at the left */
+				if( col > 0 )
+				{
+					startX -= margin;
+					sizeX += margin;
+					lsgrmLog << "1\n";
+				}
+				else
+					lsgrmLog << "0\n";
+
+				/* Extract the tile */
+				auto extractPtr = Extractor::New();
+				extractPtr->SetStartX(startX);
+				extractPtr->SetStartY(startY);
+				extractPtr->SetSizeX(sizeX);
+				extractPtr->SetSizeY(sizeY);
+
+				/* Write the tile to the ouput directory. */
+				auto writer = Writer::New();
+				writer->SetFileName(ouputDirectoryPath + "tile_" + std::to_string(row) +
+									"_" + std::to_string(col) + ".tif");
+
+				/* OTB streaming pipeline */
+				extractPtr->SetInput(reader->GetOutput());
+				writer->SetInput(extractPtr->GetOutput());
+				writer->Update();
+				
+			} // end for(unsigned int col = 0; col < nbTilesX; ++col)
+
+		} // for(unsigned int row = 0; row < nbTilesY; ++row)
+
+		/* Close the lsgrm info file. */
+		lsgrmLog.close();
+	}
+} // end of namespace lsgrm
+#endif
diff --git a/Code/lsrmBaatzSegmenter.h b/Code/lsrmBaatzSegmenter.h
new file mode 100644
index 0000000000000000000000000000000000000000..cc5a951e7b34573874e030c094bf204bdef0d18e
--- /dev/null
+++ b/Code/lsrmBaatzSegmenter.h
@@ -0,0 +1,104 @@
+#ifndef __LSRM_BAATZ_SEGMENTER_H
+#define __LSRM_BAATZ_SEGMENTER_H
+#include "lsrmSegmenter.h"
+/*
+ * Tutorial : Implementation of the Baatz & Schape criterion
+ *
+ * Details about the criterion can be found in the following publication:
+ *
+ * Martin Baatz and Arno Schape. Multiresolution segmentation: an optimization approach for high quality multi-scale image segmentation.
+ * Angewandte Geographische Informationsverarbeitung XII, pages 12–23, 2000.
+ *
+ * The steps are ordered has to be followed in a chronological way for a full understanding.
+ * This tutorial does not aim at explaining all the details of the large scale region merging
+ * but gives all the elements to use it properly.
+ */
+
+namespace lsrm
+{
+	/*
+	 * Step 1 :
+	 * We define the specific attributes required for the Baatz & Schape criterion.
+	 * Regions are represented by nodes in a graph and the lsrm library has an internal
+	 * representation of the nodes with the class lsrm::Node.
+	 * --> A node contains a unique number to be indentified, it corresponds to the vectorized
+	 *     coordinates of its first pixel and its name is m_Id.
+	 *     To retrieve the 2D coordinates from the value of m_Id, it is necessary to know the dimension
+	 *     of the image (the number of rows and columns), then :
+	 *     x = m_Id % columns and y = m_Id / columns.
+	 *
+	 * --> A node contains the perimeter (m_Perimeter) of its corresponding region, it is the length of the external
+	 *     boundary of the region. For example a region of one pixel has a perimeter of 4.
+	 *
+	 * --> A node contains the area (m_Area) of its corresponding region, it is the number of pixels within the
+	 *     region.
+	 *
+	 * --> A node contains the minimum rectangular bounding box (parrallel to the image axis) of the region (m_Bbox).
+	 *     the boundix box is determined by the structure lsrm::BoundingBox with the coordinates of the upper left
+	 *     corner pixel and the width and the height (m_UX, m_UY, m_W, m_H).
+	 *
+	 * After Reading the paper about the Baatz & Schape criterion, we can conclude that we have to define additional
+	 * spectral attributes for a node: its mean vector, the square mean vector, the sum of all the pixels
+	 * and the standard deviation vector.
+	 *
+	 * We define a new class BaatzNode which inherits from the template structure Node. Notice that the template
+	 * is the derived node type, then the template type is BaatzNode.
+	 */
+	struct BaatzNode : Node<BaatzNode>
+	{
+		std::vector<float> m_Means;
+		std::vector<float> m_SquareMeans;
+		std::vector<float> m_SpectralSum;
+		std::vector<float> m_Std;
+	};
+
+	/*
+	 * Step 2
+	 *
+	 * The Baatz & Schape criterion has user-defined parameters, the spectral weight
+	 * which is the relative importance given to the spectral components and the shape
+	 * weight for the shape components.
+	 *
+	 * We define then a structure wrapping these two parameters we called BaatzParam.
+	 */
+	struct BaatzParam
+	{
+		float m_SpectralWeight;
+		float m_ShapeWeight;
+	};
+
+	/*
+	 * Step 3 : 
+	 *
+	 */
+	
+	template<class TImage>
+	class BaatzSegmenter : public Segmenter< TImage, BaatzNode, BaatzParam>
+	{
+	public:
+
+		/* Some convenient typedefs */
+		typedef Segmenter<TImage, BaatzNode, BaatzParam> Superclass;
+		typedef TImage ImageType;
+		typedef BaatzParam ParameterType;
+		typedef typename Superclass::GraphType GraphType;
+		typedef BaatzNode NodeType;
+		typedef typename Superclass::EdgeType EdgeType;
+		typedef typename Superclass::NodePointerType NodePointerType;
+		typedef typename Superclass::GraphOperatorType GraphOperatorType;
+		typedef GraphToOtbImage<GraphType> IOType;
+
+		BaatzSegmenter();
+		
+		void Update();
+		float ComputeMergingCost(NodePointerType n1, NodePointerType n2);
+		void UpdateSpecificAttributes(NodePointerType n1, NodePointerType n2);
+		void InitFromImage();
+	};
+	
+} // end of namespace lsrm
+#include "lsrmBaatzSegmenter.txx"
+#endif
+
+
+
diff --git a/Code/lsrmBaatzSegmenter.txx b/Code/lsrmBaatzSegmenter.txx
new file mode 100644
index 0000000000000000000000000000000000000000..f751ccd1208a6b7cf24185712bf33a08b629f3ac
--- /dev/null
+++ b/Code/lsrmBaatzSegmenter.txx
@@ -0,0 +1,141 @@
+#ifndef __LSRM_BAATZ_SEGMENTER_TXX
+#define __LSRM_BAATZ_SEGMENTER_TXX
+#include <otbImageFileReader.h>
+#include <itkImageRegionIterator.h>
+
+namespace lsrm
+{
+
+	template<class TImage>
+	BaatzSegmenter<TImage>::BaatzSegmenter() : Superclass()
+	{
+	}
+	
+	template<class TImage>
+	void
+	BaatzSegmenter<TImage>::InitFromImage()
+	{
+		typedef itk::ImageRegionIterator<TImage> ImageIterator;
+
+		this->m_ImageWidth = this->m_InputImage->GetLargestPossibleRegion().GetSize()[0];
+		this->m_ImageHeight =this->m_InputImage->GetLargestPossibleRegion().GetSize()[1];
+		this->m_NumberOfComponentsPerPixel = this->m_InputImage->GetNumberOfComponentsPerPixel();
+
+		std::size_t idx = 0;
+		ImageIterator it(this->m_InputImage, this->m_InputImage->GetLargestPossibleRegion());
+		for(it.GoToBegin(); !it.IsAtEnd(); ++it)
+		{
+			this->m_Graph.m_Nodes[idx]->m_Means.reserve(this->m_NumberOfComponentsPerPixel);
+			this->m_Graph.m_Nodes[idx]->m_SquareMeans.reserve(this->m_NumberOfComponentsPerPixel);
+			this->m_Graph.m_Nodes[idx]->m_SpectralSum.reserve(this->m_NumberOfComponentsPerPixel);
+			this->m_Graph.m_Nodes[idx]->m_Std.assign(this->m_NumberOfComponentsPerPixel, 0.0f);
+
+			for(std::size_t b = 0; b < this->m_NumberOfComponentsPerPixel; ++b)
+			{
+				this->m_Graph.m_Nodes[idx]->m_Means.push_back(it.Get()[b]);
+				this->m_Graph.m_Nodes[idx]->m_SquareMeans.push_back((it.Get()[b])*(it.Get()[b]));
+				this->m_Graph.m_Nodes[idx]->m_SpectralSum.push_back(it.Get()[b]);
+			}	
+			++idx;
+		}
+	}
+
+	template<class TImage>
+	float
+	BaatzSegmenter<TImage>::ComputeMergingCost(NodePointerType n1, NodePointerType n2)
+	{
+		const unsigned int a1 = n1->m_Area, a2 = n2->m_Area, a_sum = a1 + a2;
+
+		float spect_cost = 0.0f;
+		float mean, square_mean, sum, std;
+
+		for (unsigned int b = 0; b < this->m_NumberOfComponentsPerPixel; b++)
+		{
+			mean = ((a1 * n1->m_Means[b]) + (a2 * n2->m_Means[b])) / a_sum;
+			square_mean = n1->m_SquareMeans[b] + n2->m_SquareMeans[b];
+			sum = n1->m_SpectralSum[b] + n2->m_SpectralSum[b];
+			std = std::sqrt((square_mean - 2*mean*sum + a_sum * mean* mean) / a_sum);
+			spect_cost += (a_sum * std - a1 * n1->m_Std[b] - a2 * n2->m_Std[b]);
+		}
+		spect_cost *= this->m_Param.m_ShapeWeight;
+
+		if(spect_cost < this->m_Threshold)
+		{
+			float shape_cost, smooth_f, compact_f;
+
+ 			// Compute the shape merging cost
+			const float p1 = static_cast<float>(n1->m_Perimeter);
+			const float p2 = static_cast<float>(n2->m_Perimeter);
+			const unsigned int boundary = (GraphOperatorType::FindEdge(n1, n2))->m_Boundary;
+			const float p3 = p1 + p2 - 2 * static_cast<float>(boundary);
+			
+			const lp::BoundingBox merged_bbox = lp::ContourOperations::MergeBoundingBoxes(n1->m_Bbox, n2->m_Bbox);
+			const float bb1_perimeter = static_cast<float>(2*n1->m_Bbox.m_W + 2*n1->m_Bbox.m_H);
+			const float bb2_perimeter = static_cast<float>(2*n2->m_Bbox.m_W + 2*n2->m_Bbox.m_H);
+			const float mbb_perimeter = static_cast<float>(2 * merged_bbox.m_W + 2 * merged_bbox.m_H);
+
+			smooth_f = a_sum*p3/mbb_perimeter - a1*p1/bb1_perimeter - a2*p2/bb2_perimeter;
+			compact_f = a_sum*p3/std::sqrt(a_sum) - a1*p1/std::sqrt(a1) - a2*p2/std::sqrt(a2);
+
+			shape_cost = this->m_Param.m_ShapeWeight * compact_f + (1-this->m_Param.m_ShapeWeight) * smooth_f;
+
+			return (spect_cost + (1-this->m_Param.m_ShapeWeight)*shape_cost);
+		}
+		else
+			return spect_cost;
+	}
+
+	template<class TImage>
+	void
+	BaatzSegmenter<TImage>::UpdateSpecificAttributes(NodePointerType n1, NodePointerType n2)
+	{
+		const float a1 = static_cast<float>(n1->m_Area);
+		const float a2 = static_cast<float>(n2->m_Area);
+		const float a_sum = a1 + a2;
+
+		for(unsigned int b = 0; b < this->m_NumberOfComponentsPerPixel; ++b)
+		{
+			n1->m_Means[b] = (a1 * n1->m_Means[b] + a2 * n2->m_Means[b]) / a_sum;
+			n1->m_SquareMeans[b] += n2->m_SquareMeans[b];
+			n1->m_SpectralSum[b] += n2->m_SpectralSum[b];
+			n1->m_Std[b] = std::sqrt((n1->m_SquareMeans[b] - 2 * n1->m_Means[b] * n1->m_SpectralSum[b] +
+									  a_sum * n1->m_Means[b] * n1->m_Means[b]) / a_sum);
+		}
+	}
+
+	template<class TImage>
+	void
+	BaatzSegmenter<TImage>::Update()
+	{
+		GraphOperatorType::InitNodes(this->m_InputImage, *this, FOUR);
+		bool prev_merged = false;
+
+		if(this->m_NumberOfIterations > 0)
+		{
+			prev_merged = GraphOperatorType::PerfomAllIterationsWithLMBFAndConstThreshold(*this);
+																			
+			this->m_Complete = !prev_merged;
+			
+			if(prev_merged && this->m_DoBFSegmentation)
+			{
+				prev_merged = GraphOperatorType::PerfomAllIterationsWithBFAndConstThreshold(*this);
+				this->m_Complete = !prev_merged;  
+			}
+		}
+		else
+		{
+			assert(this->m_DoBFSegmentation == true);
+			prev_merged = GraphOperatorType::PerfomAllIterationsWithBFAndConstThreshold(*this);
+			this->m_Complete = !prev_merged;
+		}
+	}
+} // end of namespace lsrm
+
+#endif
+
+
+
+
+
+
+
diff --git a/Code/lsrmDataStructures.h b/Code/lsrmDataStructures.h
new file mode 100644
index 0000000000000000000000000000000000000000..eca299bdd941c3986ddef30b3f5883aa93ff3096
--- /dev/null
+++ b/Code/lsrmDataStructures.h
@@ -0,0 +1,46 @@
+#ifndef __LSRM_DATA_STRUCTURES_H
+#define __LSRM_DATA_STRUCTURES_H
+#include <stdexcept>
+#include <bitset>
+#include <vector>
+
+namespace lsrm
+{	
+
+	struct BoundingBox
+	{
+		/* X coordinate of the upper left. */
+		long unsigned int m_UX;
+		
+		/* Y coordinate of the upper left. */
+		long unsigned int m_UY;
+		
+		/* Width */
+		unsigned int m_W;
+
+		/* Height */
+		unsigned int m_H;	
+	};
+	
+	BoundingBox MergeBoundingBoxes(const BoundingBox& bb1,
+								   const BoundingBox& bb2)
+	{
+		long unsigned int min_ux, min_uy, max_xw, max_yh;
+		BoundingBox bb;
+
+		min_ux = std::min(bb1.m_UX, bb2.m_UX);
+		min_uy = std::min(bb1.m_UY, bb2.m_UY);
+		max_xw = std::max(bb1.m_UX + bb1.m_W, bb2.m_UX + bb2.m_W);
+		max_yh = std::max(bb1.m_UY + bb1.m_H, bb2.m_UY + bb2.m_H);
+
+		bb.m_UX = min_ux;
+		bb.m_UY = min_uy;
+		bb.m_W = max_xw - min_ux;
+		bb.m_H = max_yh - min_uy;
+
+		return bb;
+	}
+	
+} // end of namespace lsrm
+
+#endif
diff --git a/Code/lsrmGetInternalMemory.h b/Code/lsrmGetInternalMemory.h
new file mode 100644
index 0000000000000000000000000000000000000000..a5a299c87c3628802e7802948ab0bf7f394da357
--- /dev/null
+++ b/Code/lsrmGetInternalMemory.h
@@ -0,0 +1,96 @@
+/*
+ * Author:  David Robert Nadeau
+ * Site:    http://NadeauSoftware.com/
+ * License: Creative Commons Attribution 3.0 Unported License
+ *          http://creativecommons.org/licenses/by/3.0/deed.en_US
+ */
+#if defined(_WIN32)
+#include <Windows.h>
+
+#elif defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
+#include <unistd.h>
+#include <sys/types.h>
+#include <sys/param.h>
+#if defined(BSD)
+#include <sys/sysctl.h>
+#endif
+
+#else
+#error "Unable to define getMemorySize( ) for an unknown OS."
+#endif
+
+
+
+/**
+ * Returns the size of physical memory (RAM) in bytes.
+ */
+size_t getMemorySize( )
+{
+#if defined(_WIN32) && (defined(__CYGWIN__) || defined(__CYGWIN32__))
+	/* Cygwin under Windows. ------------------------------------ */
+	/* New 64-bit MEMORYSTATUSEX isn't available.  Use old 32.bit */
+	MEMORYSTATUS status;
+	status.dwLength = sizeof(status);
+	GlobalMemoryStatus( &status );
+	return (size_t)status.dwTotalPhys;
+
+#elif defined(_WIN32)
+	/* Windows. ------------------------------------------------- */
+	/* Use new 64-bit MEMORYSTATUSEX, not old 32-bit MEMORYSTATUS */
+	MEMORYSTATUSEX status;
+	status.dwLength = sizeof(status);
+	GlobalMemoryStatusEx( &status );
+	return (size_t)status.ullTotalPhys;
+
+#elif defined(__unix__) || defined(__unix) || defined(unix) || (defined(__APPLE__) && defined(__MACH__))
+	/* UNIX variants. ------------------------------------------- */
+	/* Prefer sysctl() over sysconf() except sysctl() HW_REALMEM and HW_PHYSMEM */
+
+#if defined(CTL_HW) && (defined(HW_MEMSIZE) || defined(HW_PHYSMEM64))
+	int mib[2];
+	mib[0] = CTL_HW;
+#if defined(HW_MEMSIZE)
+	mib[1] = HW_MEMSIZE;            /* OSX. --------------------- */
+#elif defined(HW_PHYSMEM64)
+	mib[1] = HW_PHYSMEM64;          /* NetBSD, OpenBSD. --------- */
+#endif
+	int64_t size = 0;               /* 64-bit */
+	size_t len = sizeof( size );
+	if ( sysctl( mib, 2, &size, &len, NULL, 0 ) == 0 )
+		return (size_t)size;
+	return 0L;			/* Failed? */
+
+#elif defined(_SC_AIX_REALMEM)
+	/* AIX. ----------------------------------------------------- */
+	return (size_t)sysconf( _SC_AIX_REALMEM ) * (size_t)1024L;
+
+#elif defined(_SC_PHYS_PAGES) && defined(_SC_PAGESIZE)
+	/* FreeBSD, Linux, OpenBSD, and Solaris. -------------------- */
+	return (size_t)sysconf( _SC_PHYS_PAGES ) *
+		(size_t)sysconf( _SC_PAGESIZE );
+
+#elif defined(_SC_PHYS_PAGES) && defined(_SC_PAGE_SIZE)
+	/* Legacy. -------------------------------------------------- */
+	return (size_t)sysconf( _SC_PHYS_PAGES ) *
+		(size_t)sysconf( _SC_PAGE_SIZE );
+
+#elif defined(CTL_HW) && (defined(HW_PHYSMEM) || defined(HW_REALMEM))
+	/* DragonFly BSD, FreeBSD, NetBSD, OpenBSD, and OSX. -------- */
+	int mib[2];
+	mib[0] = CTL_HW;
+#if defined(HW_REALMEM)
+	mib[1] = HW_REALMEM;		/* FreeBSD. ----------------- */
+#elif defined(HW_PYSMEM)
+	mib[1] = HW_PHYSMEM;		/* Others. ------------------ */
+#endif
+	unsigned int size = 0;		/* 32-bit */
+	size_t len = sizeof( size );
+	if ( sysctl( mib, 2, &size, &len, NULL, 0 ) == 0 )
+		return (size_t)size;
+	return 0L;			/* Failed? */
+#endif /* sysctl and sysconf variants */
+
+#else
+	return 0L;			/* Unknown OS. */
+#endif
+}
diff --git a/Code/lsrmGraph.h b/Code/lsrmGraph.h
new file mode 100644
index 0000000000000000000000000000000000000000..b55a384bb41552b819a10f7e62c8a533054afe17
--- /dev/null
+++ b/Code/lsrmGraph.h
@@ -0,0 +1,91 @@
+#ifndef __LSRM_GRAPH_H
+#define __LSRM_GRAPH_H
+#include "lsrmDataStructures.h"
+#include "lpContour.h"
+
+namespace lsrm
+{
+	struct BaseNode
+	{
+		/* Node already merged. */
+		bool m_Valid;
+		
+		/* Node has to be removed from the graph. */
+		bool m_Expired;
+
+		/* Does the node merge at the previous iteration */
+		bool m_IsMerged;
+		
+		/* Perimeter of the region */
+		unsigned int m_Perimeter;
+
+		/* Area (number of inner pixels) of the region */
+		unsigned int m_Area;
+
+		/*
+		  Node is identified by the location
+		  of the first pixel of the region.
+		 */
+		long unsigned int m_Id;
+
+		/*
+		  Bounding box of the region
+		  in the image.
+		 */
+		lp::BoundingBox m_Bbox;
+
+		/* Contour of the shape */
+		lp::Contour m_Contour;
+	};
+
+	template<class DerivedNode>
+		struct NeighborType
+		{
+			typedef std::weak_ptr<DerivedNode> WeakDerived;
+			typedef std::shared_ptr<DerivedNode> SharedDerived;
+
+			WeakDerived m_Target;
+			float  m_Cost;
+			unsigned int m_Boundary;
+			bool m_CostUpdated;
+
+		    NeighborType(WeakDerived ptr, double w, unsigned int c) :
+			m_Target(ptr), m_Cost(w), m_Boundary(c), m_CostUpdated(false) {}
+			
+			inline SharedDerived GetRegion()
+				{
+					SharedDerived ptr(m_Target.lock());
+					if(!ptr)
+						throw std::runtime_error("lss_GenericLMBFRegionMergingHandler.h - NeighborType::GetRegion - Region pointer is not valid");
+					return ptr;
+				}
+		};
+	
+
+	template<class DerivedNode>
+		struct Node : BaseNode
+	{
+		typedef NeighborType<DerivedNode> CRPTNeighborType;
+		std::vector<CRPTNeighborType> m_Edges;
+	};
+
+	template<class TNode>
+	struct Graph
+	{
+		typedef TNode NodeType;
+		typedef std::shared_ptr<NodeType> NodePointerType;
+		typedef typename NodeType::CRPTNeighborType EdgeType;
+		typedef std::vector<NodePointerType> NodeListType;
+		typedef typename NodeListType::iterator NodeIteratorType;
+		typedef typename NodeListType::const_iterator NodeConstIteratorType;
+		typedef std::vector<EdgeType> EdgeListType;
+		typedef typename EdgeListType::iterator EdgeIteratorType;
+		typedef typename EdgeListType::const_iterator EdgeConstIteratorType;
+		
+		std::vector< NodePointerType > m_Nodes; 
+	};
+	
+} // end of namespace lsrm
+#endif
+
+
diff --git a/Code/lsrmGraphOperations.h b/Code/lsrmGraphOperations.h
new file mode 100644
index 0000000000000000000000000000000000000000..1792357811315ece2379f0a6d4160101320067a2
--- /dev/null
+++ b/Code/lsrmGraphOperations.h
@@ -0,0 +1,223 @@
+#ifndef __LSRM_GRAPH_OPERATIONS_H
+#define __LSRM_GRAPH_OPERATIONS_H
+#include "lsrmGraph.h"
+#include "lsrmNeighborhood.h"
+#include <iostream>
+#include <cassert>
+#include <limits>
+#include <map>
+#include <utility>
+#include <set>
+
+namespace lsrm
+{
+	template<class TSegmenter>
+	class GraphOperations
+	{
+	public:
+		
+		/* Some convenient typedefs */
+		typedef TSegmenter SegmenterType;
+		typedef typename SegmenterType::ImageType ImageType;
+		typedef typename SegmenterType::GraphType GraphType;
+		typedef typename GraphType::NodeType NodeType;
+		typedef typename GraphType::EdgeType EdgeType;
+		typedef typename GraphType::NodePointerType NodePointerType;
+		typedef typename GraphType::NodeListType NodeList;
+		typedef typename GraphType::NodeIteratorType NodeIterator;
+		typedef typename GraphType::NodeConstIteratorType NodeConstIterator;
+		typedef typename GraphType::EdgeListType EdgeList;
+		typedef typename GraphType::EdgeIteratorType EdgeIterator;
+		typedef typename GraphType::EdgeConstIteratorType EdgeConstIterator;
+
+		using ContourOperator = lp::ContourOperations;
+
+
+		/*
+		 * Given the size of the input image and the mask of the
+		 * neighborhood, we initialize a new graph of nodes
+		 *
+		 * @params:
+		 * GraphType& graph: reference to a graph of nodes
+		 * const unsigned int width: width of the input image
+		 * const unsigned int height: height of the input image
+		 * CONNECTIVITY mask : mask of the neighborhood (4X4 or 8X8)
+		 */
+		static void InitNodes(ImageType * inputImg,
+							  SegmenterType& seg,
+							  CONNECTIVITY mask);
+
+		/*
+		 * Given a graph of nodes, we explore all the nodes
+		 * and for each node we compute his merging costs
+		 * with all its neighboring nodes given a function
+		 * to compute the merging cost between two nodes.
+		 *
+		 * @params:
+		 * GraphType& graph: reference to the graph of nodes
+		 * float(*fptr)(NodeType*, NodeType*): pointer to the function
+		 * to compute the merging cost between two adjacent nodes.
+		 */
+		static void UpdateMergingCosts(SegmenterType& seg);
+
+		/*
+		 * Given a node A, we analyse its best node B.
+		 * If the node A is also node B's best node
+		 * then it returns a pointer to node B if node A 's id
+		 * is smaller or a pointer to node A if node B's id is
+		 * smaller
+		 * else it returns a null pointer.
+		 * (Local Mutual Best Fitting Heuristic)
+		 * 
+		 * @params:
+		 * NodeType * a : Pointer to the node A
+		 * float t : threshold of the criterion
+		 */
+		static NodePointerType CheckLMBF(NodePointerType, float t);
+
+		/*
+		 * Given a node A, we analyse its best node B.
+		 * If the criterion is checked and the node B
+		 * is valid then it returns a pointer to the
+		 * node A if node's A id is smaller or a pointer
+		 * to node B if node B's id is smaller
+		 * else it returns a null pointer.
+		 *
+		 * @params:
+		 * NodeType * a : pointer to node A
+		 * float t : threshold of the criterion
+		 */
+		static NodePointerType CheckBF(NodePointerType a, float t);
+
+		/*
+		 * Given the current node and the target node, it returns
+		 * the edge from the current node targeting to the target
+		 * node.
+		 *
+		 * @params
+		 * const NodeType * n : pointer to the current node.
+		 * const NodeType * target : pointer to the target node.
+		 * @return an iterator pointing to the candidate edge.
+		 */
+		static EdgeIterator FindEdge(NodePointerType n, NodePointerType target);
+
+		/*
+		 * Given a node a and the node b to be merged into node a,
+		 * it updates the neighbors of node a with respect to the
+		 * neighbors of node b.
+		 *
+		 * @params
+		 * NodeType * a : pointer to node a.
+		 * NodeType * b : pointer to node b.
+		 */
+		static void UpdateNeighbors(NodePointerType a, NodePointerType b);
+		
+		/*
+		 * Given 2 nodes A and B (node B being merged into node A)
+		 * we update the internal attributes of node A with respect
+		 * to node B.
+		 *
+		 * @params:
+		 * NodeType * a: pointer to node A.
+		 * NodeType * b: pointer to node B.
+		 */
+		static void UpdateInternalAttributes(NodePointerType a,
+											 NodePointerType b,
+											 const unsigned int width,
+											 const unsigned int height);
+
+		/*
+		 * Given a graph, it removes all the expired nodes.
+		 *
+		 * @params
+		 * GraphType& graph : reference to the graph.
+		 */
+		static void RemoveExpiredNodes(GraphType& graph);
+
+
+		/*
+		 * Given a graph, a region merging algorithm, a threshold
+		 * and the dimension of the image, it performs one iteration
+		 * of the merging process using the local mutual best fitting
+		 * heuristic.
+		 *
+		 * @params
+		 * GraphType& graph : reference to the graph
+		 * SegmenterType& seg : reference to the region merging algorithm.
+		 * const float threshold : threshold for this iteration.
+		 * const unsigned int width : width of the image.
+		 * const unsigned int height : height of the image.
+		 *
+		 * @return a boolean pointing out if there was at least a fusion
+		 * of nodes.
+		 */
+		static bool PerfomOneIterationWithLMBF(SegmenterType& seg);
+
+		/*
+		 * Given a graph, a region merging algorithm, a threshold
+		 * and the dimension of the image, it performs one iteration
+		 * of the merging process using the best fitting
+		 * heuristic.
+		 *
+		 * @params
+		 * GraphType& graph : reference to the graph
+		 * SegmenterType& seg : reference to the region merging algorithm.
+		 * const float threshold : threshold for this iteration.
+		 * const unsigned int width : width of the image.
+		 * const unsigned int height : height of the image.
+		 *
+		 * @return a boolean pointing out if there was at least a fusion
+		 * of nodes.
+		 */
+		static bool PerfomOneIterationWithBF(SegmenterType& seg);
+
+		/*
+		 * Given a graph, a region merging algorithm, a threshold,
+		 * the number of iterations to apply and the dimension of the image,
+		 * it performs all the iterations of the merging process using the
+		 * local mutual best fitting heuristic.
+		 * This method can be used when the threshold is constant during
+		 * the region merging process.
+		 *
+		 * @params
+		 * GraphType& graph : reference to the graph
+		 * SegmenterType& seg : reference to the region merging algorithm.
+		 * const float threshold : threshold for this iteration.
+		 * const unsigned int numberOfIterations: number of iteration to perform.
+		 * const unsigned int width : width of the image.
+		 * const unsigned int height : height of the image.
+		 *
+		 * @return a boolean pointing out if there was at least a fusion
+		 * of nodes.
+		 */
+		static bool PerfomAllIterationsWithLMBFAndConstThreshold(SegmenterType& seg);
+
+		/*
+		 * Given a graph, a region merging algorithm, a threshold,
+		 * the number of iterations to apply and the dimension of the image,
+		 * it performs all the iterations of the merging process using the
+		 * best fitting heuristic.
+		 * This method can be used when the threshold is constant during
+		 * the region merging process.
+		 *
+		 * @params
+		 * GraphType& graph : reference to the graph
+		 * SegmenterType& seg : reference to the region merging algorithm.
+		 * const float threshold : threshold for this iteration.
+		 * const unsigned int numberOfIterations: number of iteration to perform.
+		 * const unsigned int width : width of the image.
+		 * const unsigned int height : height of the image.
+		 *
+		 * @return a boolean pointing out if there was at least a fusion
+		 * of nodes.
+		 */
+		static bool PerfomAllIterationsWithBFAndConstThreshold(SegmenterType& seg);
+	};
+} // end of namespace lsrm
+
+#include "lsrmGraphOperations.txx"
+#endif
+
+
+
+
diff --git a/Code/lsrmGraphOperations.txx b/Code/lsrmGraphOperations.txx
new file mode 100644
index 0000000000000000000000000000000000000000..0e4ac444a1c02c92899061569b3872fa5941c43e
--- /dev/null
+++ b/Code/lsrmGraphOperations.txx
@@ -0,0 +1,456 @@
+#ifndef __LSRM_GRAPH_OPERATIONS_TXX
+#define __LSRM_GRAPH_OPERATIONS_TXX
+#include <otbImageFileReader.h>
+
+namespace lsrm
+{	
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::InitNodes(ImageType * inputImg,
+												SegmenterType& seg,
+												CONNECTIVITY mask)
+	{
+		unsigned int width, height;
+		
+		{
+			width = inputImg->GetLargestPossibleRegion().GetSize()[0];
+			height = inputImg->GetLargestPossibleRegion().GetSize()[1];
+		}
+		
+		const long unsigned int num_nodes = width * height;
+
+		seg.m_Graph.m_Nodes.reserve(num_nodes);
+
+		for(long unsigned int i = 0;
+			i < num_nodes;
+			++i)
+		{
+			NodePointerType n(new NodeType);
+			n->m_Id = i;
+			n->m_Valid = true;
+			n->m_Expired = false;
+			n->m_IsMerged = true; // force to compute costs for the first iteration
+			n->m_Perimeter = 4;
+			n->m_Area = 1;
+			n->m_Bbox.m_UX = i % width;
+			n->m_Bbox.m_UY = i / width;
+			n->m_Bbox.m_W = 1;
+			n->m_Bbox.m_H = 1;
+
+			// An initial contour is the one aroun a pixel
+			ContourOperator::Push1(n->m_Contour);
+			ContourOperator::Push2(n->m_Contour);
+			ContourOperator::Push3(n->m_Contour);
+			ContourOperator::Push0(n->m_Contour);
+			
+			seg.m_Graph.m_Nodes.push_back(n);
+		}
+
+		if(mask == FOUR)
+		{
+			for(auto& r : seg.m_Graph.m_Nodes)
+			{
+				long int neighborhood[4];
+				FOURNeighborhood(neighborhood, r->m_Id, width, height);
+				for(short j = 0; j < 4; ++j)
+				{
+					if(neighborhood[j] > -1)
+						r->m_Edges.push_back(EdgeType( seg.m_Graph.m_Nodes[neighborhood[j]], 0, 1));
+				}
+			}
+		}
+		else
+		{
+			for(auto& r : seg.m_Graph.m_Nodes)
+			{
+				long int neighborhood[8];
+				EIGHTNeighborhood(neighborhood, r->m_Id, width, height);
+				for(short j = 0; j < 8; ++j)
+				{
+					if(neighborhood[j] > -1)
+					{
+						if(j % 2 > 0)
+							r->m_Edges.push_back(EdgeType( seg.m_Graph.m_Nodes[neighborhood[j]], 0, 0));
+						else
+							r->m_Edges.push_back(EdgeType( seg.m_Graph.m_Nodes[neighborhood[j]], 0, 1));
+					}
+				}
+			}
+		}
+		seg.InitFromImage();
+	}
+
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::UpdateMergingCosts(SegmenterType& seg)
+	{		
+		float min_cost;
+		long unsigned int min_id;
+		std::size_t idx, min_idx;
+
+		for(auto& r : seg.m_Graph.m_Nodes)
+		{
+			for(auto& edge : r->m_Edges)
+				edge.m_CostUpdated = false;
+		}
+
+		for(auto& r : seg.m_Graph.m_Nodes)
+		{
+			min_cost = std::numeric_limits<float>::max();
+			idx = 0;
+			min_idx = 0;
+
+			r->m_Expired = false;
+			r->m_Valid = true;
+
+			for(auto& edge : r->m_Edges)
+			{
+				auto neighborR = edge.GetRegion();
+
+				// Compute the cost if necessary
+				if(!edge.m_CostUpdated && (neighborR->m_IsMerged || r->m_IsMerged))
+				{
+					auto edgeFromNeighborToR = FindEdge(neighborR, r);
+					edge.m_Cost = seg.ComputeMergingCost(r, neighborR);
+					edgeFromNeighborToR->m_Cost = edge.m_Cost;
+					edge.m_CostUpdated = true;
+					edgeFromNeighborToR->m_CostUpdated = true;
+				}
+
+				// Check if the cost of the edge is the minimum
+				if(min_cost > edge.m_Cost)
+				{
+					min_cost = edge.m_Cost;
+					min_id = neighborR->m_Id;
+					min_idx = idx;
+				}
+				else if(min_cost == edge.m_Cost)
+				{
+					if(min_id > neighborR->m_Id)
+					{
+						min_id = neighborR->m_Id;
+						min_idx = idx;
+					}
+				}
+				++idx;	
+			}
+
+			assert(min_idx < r->m_Edges.size());
+			std::swap(r->m_Edges[0], r->m_Edges[min_idx]);
+				
+		}
+
+		// Reset the merge flag for all the regions.
+		for(auto& r : seg.m_Graph.m_Nodes)
+			r->m_IsMerged = false;
+	}
+
+	template<class TSegmenter>
+	typename GraphOperations<TSegmenter>::NodePointerType
+	GraphOperations<TSegmenter>::CheckLMBF(NodePointerType a, float t)
+	{
+		if(a->m_Valid)
+		{
+			float cost = a->m_Edges.front().m_Cost;
+			
+			if(cost < t)
+			{
+				NodePointerType b = a->m_Edges.front().GetRegion();
+
+				if( b->m_Valid)
+				{
+					NodePointerType best_b = b->m_Edges.front().GetRegion();
+
+					if(a == best_b)
+					{
+						if(a->m_Id < b->m_Id)
+							return a;
+						else
+							return b;
+					}
+					else
+						return NodePointerType();
+				}
+				else
+					return NodePointerType();
+			}
+			else
+				return NodePointerType();
+		}
+		else
+			return NodePointerType();
+	}
+
+	template<class TSegmenter>
+	typename GraphOperations<TSegmenter>::NodePointerType
+	GraphOperations<TSegmenter>::CheckBF(NodePointerType a, float t)
+	{
+		if(a->m_Valid)
+		{
+			float cost = a->m_Edges.front().m_Cost;
+
+			if( cost < t )
+			{
+				NodePointerType b = a->m_Edges.front().GetRegion();
+
+				if(b->m_Valid)
+				{
+					if( a->m_Id < b->m_Id )
+						return a;
+					else
+						return b;
+				}
+				else
+					return NodePointerType();
+			}
+			else
+				return NodePointerType();
+		}
+		else
+			return NodePointerType();
+	}
+
+	template<class TSegmenter>
+	typename GraphOperations<TSegmenter>::EdgeIterator
+	GraphOperations<TSegmenter>::FindEdge(NodePointerType n, NodePointerType target)
+	{
+		return std::find_if(n->m_Edges.begin(), n->m_Edges.end(),[&](EdgeType& e)->bool{
+				return e.GetRegion() == target;
+			});
+	}
+
+	template<class TSegmenter>
+	void
+	GraphOperations<TSegmenter>::UpdateNeighbors(NodePointerType a, NodePointerType b)
+	{
+		unsigned int boundary;
+
+		/* Explore the neighbors of b */
+		for (auto& edge : b->m_Edges)
+		{
+			// Retrieve the edge targeting node b.
+			auto neigh_b = edge.GetRegion();
+			auto toB = FindEdge(neigh_b, b);
+
+			/* If the edge tageting to node b is the first then
+			   the corresponding node is not valid anymore. */
+			if(toB == neigh_b->m_Edges.begin())
+				neigh_b->m_Valid = false;
+
+			/* Keep in memory the boundary between node b
+			   and node neigh_b */
+			boundary = edge.m_Boundary;
+
+			/*
+			  We can remove safely the edge from node neigh_b
+			  targeting to node b
+			*/
+			neigh_b->m_Edges.erase(toB);
+
+			if(neigh_b != a)
+			{
+				/* Retrieve the edge targeting to node a. */
+				auto toA = FindEdge(neigh_b, a);
+
+				if( toA == neigh_b->m_Edges.end() )
+				{
+					/* No edge exists between node a and node neigh_b. */
+
+					/* Add an edge from node neigh_b targeting node a. */
+					neigh_b->m_Edges.push_back(EdgeType(a, 0, boundary));
+
+					/* Add an edge from node a targeting node neigh_b. */
+					a->m_Edges.push_back(EdgeType(neigh_b, 0, boundary));
+				}
+				else
+				{
+					/* An edge exists between node a and node neigh_b. */
+
+					/* Increment the boundary of the edge from node neigh_b
+					   targeting to node a. */
+					toA->m_Boundary += boundary;
+
+					/* Increment the boundary of the edge from node a
+					   targeting to node neigh_b. */
+					auto toNeighB = FindEdge(a, neigh_b);
+					toNeighB->m_Boundary += boundary;
+				}
+			}
+			
+		}
+		
+	}
+
+	template<class TSegmenter>
+	void
+	GraphOperations<TSegmenter>::UpdateInternalAttributes(NodePointerType a,
+														  NodePointerType b,
+														  const unsigned int width,
+														  const unsigned int height)
+	{
+
+		lp::BoundingBox mergedBBox;
+		lp::Contour mergedContour;		
+		ContourOperator::MergeContour(mergedContour, mergedBBox, a->m_Contour,
+									  b->m_Contour, a->m_Bbox, b->m_Bbox,
+									  a->m_Id, b->m_Id, width);
+		
+		/* Step 1: update the bounding box */
+		a->m_Bbox = mergedBBox;//MergeBoundingBoxes(a->m_Bbox, b->m_Bbox);
+
+		/* Step 2: update the contour */
+		a->m_Contour = mergedContour;
+
+		/* Step 2 : update perimeter and area attributes */
+		EdgeIterator toB = FindEdge(a, b);
+		a->m_Perimeter += (b->m_Perimeter - 2 * toB->m_Boundary);
+		a->m_Area += b->m_Area;
+			
+		/* Step 2: update the neighborhood */
+		UpdateNeighbors(a,b);
+		
+		/* Step 3: update the node' states */
+		a->m_Valid = false;
+		b->m_Valid = false;
+		b->m_Expired = true;
+		a->m_IsMerged = true;
+	}
+
+	template<class TSegmenter>
+	void
+	GraphOperations<TSegmenter>::RemoveExpiredNodes(GraphType& graph)
+	{
+		NodeIterator nit = std::remove_if(graph.m_Nodes.begin(), graph.m_Nodes.end(), [](NodePointerType r)->bool{
+				return r->m_Expired;
+			});
+		graph.m_Nodes.erase(nit, graph.m_Nodes.end());
+	}
+
+	template<class TSegmenter>
+	bool
+	GraphOperations<TSegmenter>::PerfomOneIterationWithLMBF(SegmenterType& seg)
+	{
+		bool merged = false;
+
+		/* Update the costs of merging between adjacent nodes */
+		UpdateMergingCosts(seg);
+
+		for(auto& region : seg.m_Graph.m_Nodes)
+		{
+			
+			auto res_node = CheckLMBF(region, seg.GetThreshold());
+
+			if(res_node)
+				{
+					seg.UpdateSpecificAttributes(res_node, res_node->m_Edges.front().GetRegion());
+					UpdateInternalAttributes(res_node, res_node->m_Edges.front().GetRegion(),
+											 seg.GetImageWidth(), seg.GetImageHeight());
+					merged = true;
+				}
+		}
+
+		RemoveExpiredNodes(seg.m_Graph);
+
+		if(seg.m_Graph.m_Nodes.size() < 2)
+			return false;
+
+		return merged;
+	}
+
+	template<class TSegmenter>
+	bool
+	GraphOperations<TSegmenter>::PerfomOneIterationWithBF(SegmenterType& seg)
+	{
+		bool merged = false;
+
+		/* Update the costs of merging between adjacent nodes */
+		UpdateMergingCosts(seg);
+
+		for(auto& region : seg.m_Graph.m_Nodes)
+		{
+			NodePointerType res_node = CheckBF(region, seg.GetThreshold());
+
+			if(res_node)
+				{
+					seg.UpdateSpecificAttributes(res_node, res_node->m_Edges.front().GetRegion());
+					UpdateInternalAttributes(res_node, res_node->m_Edges.front().GetRegion(),
+											 seg.GetImageWidth(), seg.GetImageHeight());
+					merged = true;
+				}
+		}
+
+		RemoveExpiredNodes(seg.m_Graph);
+
+		if(seg.m_Graph.m_Nodes.size() < 2)
+			return false;
+
+		return merged;
+	}
+
+	template<class TSegmenter>
+	bool
+	GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(SegmenterType& seg)
+	{
+		bool merged = true;
+		unsigned int iterations = 0;
+		unsigned int numberOfIterations = seg.GetNumberOfIterations();
+
+		while(merged &&
+			  iterations < numberOfIterations &&
+			  seg.m_Graph.m_Nodes.size() > 1)
+		{
+			std::cout << "." << std::flush;
+			++iterations;
+
+			merged = PerfomOneIterationWithLMBF(seg);
+		}
+		std::cout << std::endl;
+		
+		if(seg.m_Graph.m_Nodes.size() < 2)
+			return false;
+
+		return merged;
+	}
+
+	template<class TSegmenter>
+	bool
+	GraphOperations<TSegmenter>::PerfomAllIterationsWithBFAndConstThreshold(SegmenterType& seg)
+	{
+		bool merged = true;
+		unsigned int maxNumberOfIterations;
+		unsigned int numberOfIterations = seg.GetNumberOfIterations();
+		
+		if(numberOfIterations < 1)
+			maxNumberOfIterations = 75;
+		else
+			maxNumberOfIterations = numberOfIterations;
+		
+		unsigned int iterations = 0;
+
+		while(merged &&
+			  iterations < maxNumberOfIterations &&
+			  seg.m_Graph.m_Nodes.size() > 1)
+		{
+			std::cout << "." << std::flush;
+			++iterations;
+
+			merged = PerfomOneIterationWithBF(seg);
+		}
+		std::cout << std::endl;
+
+		if(seg.m_Graph.m_Nodes.size() < 2)
+			return false;
+
+		return merged;
+	}
+	
+
+} // end of namespace lsrm
+
+#endif
+
+
+
+
+
+
+
+
diff --git a/Code/lsrmGraphToOtbImage.h b/Code/lsrmGraphToOtbImage.h
new file mode 100644
index 0000000000000000000000000000000000000000..89ef969ce259a05f70f6c70432d40fab054654cf
--- /dev/null
+++ b/Code/lsrmGraphToOtbImage.h
@@ -0,0 +1,45 @@
+#ifndef __LSRM_GRAPH_TO_OTBIMAGE_H
+#define __LSRM_GRAPH_TO_OTBIMAGE_H
+#include <itkRGBPixel.h>
+#include <itkImageRegion.h>
+#include <otbImage.h>
+#include <otbVectorImage.h>
+#include <otbImageFileReader.h>
+#include <otbImageFileWriter.h>
+#include "lsrmGraph.h"
+#include <string>
+#include <stdlib.h>
+#include <time.h>
+#include "lpContour.h"
+
+namespace lsrm
+{
+	template<class TGraph>
+	class GraphToOtbImage
+	{
+	public:
+		
+		/* Some convenient typedefs */
+		typedef TGraph GraphType;
+		typedef typename GraphType::NodeType NodeType;
+		typedef std::vector< std::shared_ptr<NodeType> > NodeList;
+		typedef typename NodeList::const_iterator NodeConstIterator;
+		typedef unsigned long int LabelPixelType;
+		typedef otb::Image<LabelPixelType, 2> LabelImageType;
+		typedef unsigned char ClusterPixelType;
+		typedef otb::VectorImage<ClusterPixelType, 2> ClusteredImageType;
+		using ContourOperator = lp::ContourOperations;
+		
+
+		LabelImageType::Pointer GetLabelImage(const GraphType& graph,
+											  const unsigned int width,
+											  const unsigned int height);
+
+		ClusteredImageType::Pointer GetClusteredOutput(const GraphType& graph,
+													   const unsigned int width,
+													   const unsigned int height);	
+	};
+	
+} // end of namespace lsrm
+#include "lsrmGraphToOtbImage.txx"
+#endif
diff --git a/Code/lsrmGraphToOtbImage.txx b/Code/lsrmGraphToOtbImage.txx
new file mode 100644
index 0000000000000000000000000000000000000000..66052088b818c16e71fd8396b3b1af60b6c0980e
--- /dev/null
+++ b/Code/lsrmGraphToOtbImage.txx
@@ -0,0 +1,128 @@
+#ifndef __LSRM_GRAPH_TO_OTBIMAGE_TXX
+#define __LSRM_GRAPH_TO_OTBIMAGE_TXX
+#include "lsrmGraphToOtbImage.h"
+#include "itkImageRegionIterator.h"
+
+namespace lsrm
+{
+	template<class TGraph>
+	typename GraphToOtbImage<TGraph>::LabelImageType::Pointer
+	GraphToOtbImage<TGraph>::GetLabelImage(const GraphType& graph,
+										   const unsigned int width,
+										   const unsigned int height)
+	{
+		LabelImageType::IndexType index;
+		LabelImageType::SizeType size;
+		LabelImageType::RegionType region;
+
+		index[0] = 0; index[1] = 0;
+		size[0] = width; size[1] = height;
+		region.SetIndex(index);
+		region.SetSize(size);
+
+		LabelImageType::Pointer label_img = LabelImageType::New();
+		label_img->SetRegions(region);
+		label_img->Allocate();
+
+		using LabelImageIterator = itk::ImageRegionIterator<LabelImageType>;
+		LabelImageIterator it(label_img, label_img->GetLargestPossibleRegion());
+		for(it.GoToBegin();!it.IsAtEnd(); ++it)
+			it.Set(0);
+		
+		// Start at 1 (value 0 can be used for invalid pixels)
+		long unsigned int label = 1;
+		for(auto& region : graph.m_Nodes)
+		{
+			lp::CellLists borderPixels;
+			ContourOperator::GenerateBorderCells(borderPixels, region->m_Contour, region->m_Id, width);
+			
+			for (auto& pix: borderPixels)
+			{
+				index[0] = pix % width;
+				index[1] = pix / width;
+				label_img->SetPixel(index, label);
+			}
+			++label;
+		}
+
+		long unsigned int pixelValue = 0;
+		for(it.GoToBegin(); !it.IsAtEnd(); ++it)
+		{
+			auto pixel = it.Get();
+			if(pixel == 0)
+				it.Set(pixelValue);
+			else
+				pixelValue = pixel;
+		}
+		
+		return label_img;
+	}
+
+	template<class TGraph>
+	typename GraphToOtbImage<TGraph>::ClusteredImageType::Pointer
+	GraphToOtbImage<TGraph>::GetClusteredOutput(const GraphType& graph,
+												const unsigned int width,
+												const unsigned int height)
+	{
+		ClusteredImageType::IndexType index;
+		ClusteredImageType::SizeType size;
+		ClusteredImageType::RegionType region;
+
+		index[0] = 0; index[1] = 0;
+		size[0] = width; size[1] = height;
+		region.SetIndex(index);
+		region.SetSize(size);
+
+		ClusteredImageType::Pointer clusterImg = ClusteredImageType::New();
+		clusterImg->SetRegions(region);
+		clusterImg->SetNumberOfComponentsPerPixel(3);
+		clusterImg->Allocate();
+
+		ClusteredImageType::PixelType pixelValue;
+		pixelValue.Reserve(3);
+		pixelValue[0] = 0;
+		pixelValue[1] = 0;
+		pixelValue[2] = 0;
+
+		using ClusterImageIterator = itk::ImageRegionIterator<ClusteredImageType>;
+		ClusterImageIterator it(clusterImg, clusterImg->GetLargestPossibleRegion());
+		for(it.GoToBegin();!it.IsAtEnd(); ++it)
+			it.Set(pixelValue);
+
+		srand(time(NULL));
+		unsigned char c1, c2, c3;
+		for(auto& region : graph.m_Nodes)
+		{
+			c1 = rand() % 256;
+			c2 = rand() % 256;
+			c3 = rand() % 256;
+
+			lp::CellLists borderPixels;
+			ContourOperator::GenerateBorderCells(borderPixels, region->m_Contour, region->m_Id, width);
+			
+			for (auto& pix : borderPixels)
+			{
+				index[0] = pix % width;
+				index[1] = pix / width;
+				pixelValue[0] = c1;
+				pixelValue[1] = c2;
+				pixelValue[2] = c3;
+				clusterImg->SetPixel(index, pixelValue);
+			}
+		}
+
+		
+		for(it.GoToBegin(); !it.IsAtEnd(); ++it)
+		{
+			auto pixel = it.Get();
+			if(pixel[0] == 0 && pixel[1] == 0 && pixel[2] == 0)
+				it.Set(pixelValue);
+			else
+				pixelValue = pixel;
+		}
+		
+		return clusterImg;
+	}
+		
+} // end of namespace lsrm
+#endif
diff --git a/Code/lsrmNeighborhood.cpp b/Code/lsrmNeighborhood.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d2b0bda95db959ad4cb7d64a46c4cc69ea618aa8
--- /dev/null
+++ b/Code/lsrmNeighborhood.cpp
@@ -0,0 +1,58 @@
+#include "lsrmNeighborhood.h"
+
+namespace lsrm
+{
+	void FOURNeighborhood(long int * neighborhood,
+						  const long unsigned int id,
+						  const unsigned int width,
+						  const unsigned int height)
+	{	
+		const unsigned int x = id % width;
+		const unsigned int y = id / width;
+
+		/* top */
+		neighborhood[0] = ( y > 0 ? (id - width) : -1 );
+
+		/* right */
+		neighborhood[1] = ( x < (width - 1) ? (id + 1) : -1 );
+
+		/* bottom */
+		neighborhood[2] = ( y < (height - 1) ? (id + width) : -1 );
+
+		/* left */
+		neighborhood[3] = ( x > 0 ? (id - 1) : -1 );
+	}
+	
+	void EIGHTNeighborhood(long int * neighborhood,
+						   const long unsigned int id,
+						  const unsigned int width,
+						  const unsigned int height)
+	{
+		const unsigned int x = id % width;
+		const unsigned int y = id / width;
+
+		/* top */
+		neighborhood[0] = ( y > 0 ? (id - width) : -1 );
+
+		/* top right */
+		neighborhood[1] = ( (y > 0 && x < (width - 1) ) ? (id - width + 1) : -1 );
+		
+		/* right */
+		neighborhood[2] = ( x < (width - 1) ? (id + 1) : -1 );
+
+		/* bottom right */
+		neighborhood[3] = ( (x < (width - 1) && y < (height - 1) ) ? (id + 1 + width) : -1);
+
+		/* bottom */
+		neighborhood[4] = ( y < (height - 1) ? (id + width) : -1 );
+
+		/* bottom left */
+		neighborhood[5] = ( (y < (height - 1) && x > 0) ? (id + width - 1) : -1 );
+		
+		/* left */
+		neighborhood[6] = ( x > 0 ? (id - 1) : -1 );
+
+		/* top left */
+		neighborhood[7] = ( (x > 0 && y > 0) ? (id -width - 1) : - 1);
+	}
+} // end of namespace lsrm
diff --git a/Code/lsrmNeighborhood.h b/Code/lsrmNeighborhood.h
new file mode 100644
index 0000000000000000000000000000000000000000..4d1ef3200e14977eac507815841aa5c1dec0f8e6
--- /dev/null
+++ b/Code/lsrmNeighborhood.h
@@ -0,0 +1,19 @@
+#ifndef __LSRM_NEIGHBORHOOD_H
+#define __LSRM_NEIGHBORHOOD_H
+
+enum CONNECTIVITY{FOUR = 0, EIGHT};
+
+namespace lsrm
+{
+	void FOURNeighborhood(long int * neighborhood,
+						  const long unsigned int id,
+						  const unsigned int width,
+						  const unsigned int height);
+	
+	void EIGHTNeighborhood(long int * neighborhood,
+						   const long unsigned int id,
+						  const unsigned int width,
+						  const unsigned int height);
+} // end of namespace lsrm
+
+#endif
diff --git a/Code/lsrmSegmenter.h b/Code/lsrmSegmenter.h
new file mode 100644
index 0000000000000000000000000000000000000000..f48f9290d0abfdfa357c1b8b9a343d04003bc775
--- /dev/null
+++ b/Code/lsrmSegmenter.h
@@ -0,0 +1,145 @@
+#ifndef __LSRM_SEGMENTER_H
+#define __LSRM_SEGMENTER_H
+#include "macro-generator.h"
+#include "lsrmGraph.h"
+#include "lsrmGraphOperations.h"
+#include "lsrmGraphToOtbImage.h"
+
+namespace lsrm
+{
+	template<class TImage, class TNode, class TParam>
+	class Segmenter
+	{
+	public:
+
+		/* Some convenient typedefs */
+		
+		typedef Segmenter<TImage, TNode, TParam> Self;
+		typedef TImage ImageType;
+		typedef TNode NodeType;
+		typedef TParam ParamType;
+		typedef lsrm::Graph<NodeType> GraphType;
+		typedef typename GraphType::EdgeType EdgeType;
+		typedef GraphOperations<Self> GraphOperatorType;
+		typedef typename GraphOperatorType::NodePointerType NodePointerType;
+		typedef GraphToOtbImage<GraphType> IOType;
+		typedef typename IOType::LabelImageType LabelImageType;
+		typedef typename IOType::ClusteredImageType ClusteredImageType;
+
+		/* Default constructor and destructor */
+		
+		Segmenter(){
+			this->m_DoBFSegmentation = true;
+			this->m_NumberOfIterations = 75;
+		};
+		~Segmenter(){};
+
+		/*
+		 * This method performs the region merging segmentation.
+		 */
+		virtual void Update() = 0;
+
+		/* methods to overload */
+
+		/*
+		 * Given 2 smart adjacent node pointers (boost::shared_ptr), this
+		 * method has to compute the merging cost which is coded as a float.
+		 *
+		 * @params
+		 * NodePointerType n1 : Smart pointer to node 1
+		 * NodePointerType n2 : Smart pointer to node 2
+		 *
+		 * @return the merging cost.
+		 */
+		virtual float ComputeMergingCost(NodePointerType n1, NodePointerType n2) = 0;
+
+		/*
+		 * Given 2 smart adjacent node pointers (boost::shared_ptr), this
+		 * method merges th node n2 into the node n1 by updating the customized
+		 * attributes of the node n1.
+		 *
+		 * @params
+		 * NodePointerType n1 : Smart pointer to node 1
+		 * NodePointerType n2 : Smart pointer to node 2
+		 *
+		 */
+		virtual void UpdateSpecificAttributes(NodePointerType n1, NodePointerType n2) = 0;
+
+		/*
+		 * Given the input image, this method initializes the
+		 * internal and specific attributes of the graph.
+		 *
+		 * @params
+		 * const std::string& inputFileName : input image path
+		 */
+		virtual void InitFromImage() = 0;
+
+		/* Return the label image */
+		inline typename LabelImageType::Pointer GetLabeledClusteredOutput()
+			{
+				IOType io;
+				auto labelImg = io.GetLabelImage(this->m_Graph, this->m_ImageWidth, this->m_ImageHeight);
+				return labelImg;
+			}
+
+		inline typename ClusteredImageType::Pointer GetClusteredImageOutput()
+			{
+				IOType io;
+				auto clusteredImg = io.GetClusteredOutput(this->m_Graph, this->m_ImageWidth, this->m_ImageHeight);
+				return clusteredImg;
+			}
+		
+		/* Set methods */
+		SetMacro(bool, DoBFSegmentation);
+		SetMacro(unsigned int, NumberOfIterations);
+		SetMacro(float, Threshold);
+		SetMacro(ParamType, Param);
+		SetMacro(unsigned int, ImageWidth);
+		SetMacro(unsigned int, ImageHeight);
+		SetMacro(unsigned int, NumberOfComponentsPerPixel);
+		inline void SetInput(TImage * in){ m_InputImage = in;}
+		inline bool GetComplete(){ return this->m_Complete;}
+
+		/* Get methods */
+		GetMacro(float, Threshold);
+		GetMacro(unsigned int, ImageWidth);
+		GetMacro(unsigned int, ImageHeight);
+		GetMacro(unsigned int, NumberOfComponentsPerPixel);
+		GetMacro(unsigned int, NumberOfIterations);
+		
+		/* Graph */
+		GraphType m_Graph;
+		
+		
+	protected:
+
+		/* Boolean indicating if the segmentation procedure is achieved */
+		bool m_Complete;
+		
+		/* Activate the Best Fitting heuristic */
+		bool m_DoBFSegmentation;
+		
+		/* Number of iterations for the Local Mutual Best Fitting segmentation */
+		unsigned int m_NumberOfIterations;
+
+		/* Limit threshold for the region merging criterion  */
+		float m_Threshold;
+
+		/* Specific parameters required for the region merging criterion */
+		ParamType m_Param;
+
+		/* Image information (has to be initialized in the method InitFromImage) */
+		unsigned int m_ImageWidth; // Number of columns
+		unsigned int m_ImageHeight; // NUmber of rows
+		unsigned int m_NumberOfComponentsPerPixel; // Number of spectral bands
+
+		/* Pointer to the input image to segment */
+		TImage * m_InputImage;
+	};
+} // end of namespace lsrm
+
+#endif
+
+
+
+
diff --git a/Code/macro-generator.h b/Code/macro-generator.h
new file mode 100644
index 0000000000000000000000000000000000000000..404645fdf0cfca0a25ecbc80a440130e75d50138
--- /dev/null
+++ b/Code/macro-generator.h
@@ -0,0 +1,45 @@
+/*=========================================================================
+
+  Program: Macro Generator (MG)
+  Language: C++
+  author: Lassalle Pierre
+
+
+
+  Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
+
+
+     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.
+
+=========================================================================*/
+#ifndef __MG_H
+#define __MG_H
+
+// Generate automatically get/set methods
+#define GetMacro(type, name) \
+	type Get##name() \
+	{\
+		return m_##name; \
+	}
+
+#define SetMacro(type, name) \
+	void Set##name(type v) \
+	{ \
+		m_##name = v; \
+	}
+
+#define GetRefMacro(type, name) \
+	type& Get##name() \
+	{\
+		return m_##name; \
+	}
+
+#define GetConstMacro(type, name) \
+	const type Get##name() \
+	{\
+		return m_##name; \
+	}
+
+#endif
diff --git a/README.md b/README.md
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..b4b084d3325a1620fa3475173a1942dc3c025eab 100644
--- a/README.md
+++ b/README.md
@@ -0,0 +1,36 @@
+author: Pierre Lassalle
+
+Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
+
+This is free software under the GPL v3 licence. See
+http://www.gnu.org/licenses/gpl-3.0.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.
+
+GRM Library: Quick Start Tutorial
+
+Requirements:
+* Orfeo Toolbox library available here: http://orfeo-toolbox.org/otb/download.html
+* CMake (version 2.8 minimum)
+
+On Linux:
+
+Step 1: copy the source code into your favorite directory. In the following we assume that
+your path is PATH/GRM_SRC/
+
+Step 2: Create a binary directory: mkdir PATH/GRM_BIN/
+
+Step 3: Go to the binary directory: cd PATH/GRM_BIN/
+
+Step 5: Generate the Makefile: ccmake ../GRM_SRC/
+
+Step 6: make
+
+Step 7: cd Applications | ./RegionMergingSegmentation
+
+For more explanation about this library, go to this site:
+	http://pierre33.github.io/tuto.html
+
+Enjoy !
diff --git a/grmlib-copyright.txt b/grmlib-copyright.txt
new file mode 100644
index 0000000000000000000000000000000000000000..12510b1f666fb5f844c62fcac74551af241145d4
--- /dev/null
+++ b/grmlib-copyright.txt
@@ -0,0 +1,10 @@
+author: Pierre Lassalle
+
+Copyright (c) Centre National d'Etudes Spatiales. All rights reserved
+
+This is free software under the GPL v3 licence. See
+http://www.gnu.org/licenses/gpl-3.0.html for details.
+
+This software is distributed WITHOUT ANY WARRANTY; without even
+the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+PURPOSE.