diff --git a/Applications/lsgrmProto.cxx b/Applications/lsgrmProto.cxx
index 8a8d7e22ab84129642971df00f2ae9510ae066a9..5326dc7840f50c7180c0544fe9af111e1cdb94ff 100644
--- a/Applications/lsgrmProto.cxx
+++ b/Applications/lsgrmProto.cxx
@@ -4,58 +4,71 @@
 
 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];
-
-	/*
+
+  int myrank, nprocs;
+  MPI_Init(&argc, &argv);
+  MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+  MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+  if (myrank==0)
+    std::cout << "Number of MPI process : " << nprocs << std::endl;
+
+  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;
+
+   */
+
+  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(4096ul);
+  controller.SetTileWidth(tileWidth);
+  controller.SetTileHeight(tileHeight);
+  controller.SetNumberOfFirstIterations(niter);
+
+  // MPI configuration
+  controller.SetNumberOfProcess(nprocs);
+  controller.SetProcessRank(myrank);
+
+  // Specific parameters
+  lsrm::BaatzParam params;
+  params.m_SpectralWeight = 0.7;
+  params.m_ShapeWeight = 0.3;
+  controller.SetSpecificParameters(params);
+  controller.SetThreshold(60*60);
+
+  controller.RunSegmentation();
+
+  MPI_Finalize();
+  return 0;
 }
diff --git a/Code/lsgrmController.h b/Code/lsgrmController.h
index 2fab165a8cb439ed583ac8cfb76c317fd71f65bb..a1ad06703b8b5671524d21acf8c9156dfd80d497 100644
--- a/Code/lsgrmController.h
+++ b/Code/lsgrmController.h
@@ -33,6 +33,8 @@ namespace lsgrm
 		void SetSpecificParameters(const SegmentationParameterType& params);
 		void SetThreshold(const float& t);
 		void SetInternalMemoryAvailable(long long unsigned int v); // expecting a value in Mbytes.
+        void SetNumberOfProcess(int n) {m_NumberOfProcess = n;} ;
+        void SetProcessRank(int i) {m_ProcessRank = i;} ;
 		
 	private:
 
@@ -66,6 +68,8 @@ namespace lsgrm
 		unsigned int m_Margin;
 		unsigned int m_TileWidth;
 		unsigned int m_TileHeight;
+		int m_NumberOfProcess;
+		int m_ProcessRank;
 		std::vector<ProcessingTile> m_Tiles;
 	};
 } // end of namespace lsgrm
diff --git a/Code/lsgrmController.txx b/Code/lsgrmController.txx
index 1660abdeb7aac5991216119afd519b034330a558..ee3907791a4b8a16137bbf9d2a16f365526d7636 100644
--- a/Code/lsgrmController.txx
+++ b/Code/lsgrmController.txx
@@ -28,21 +28,29 @@ namespace lsgrm
 		}
 		
 		// Divide the input image if necessary
-		if(m_ImageDivisionActivated)
-			SplitOTBImage<ImageType>(m_InputImage, m_TileDirectory, m_TileWidth,
-									 m_TileHeight, m_Margin, m_NumberOfFirstIterations);
+		if(m_ImageDivisionActivated && (m_ProcessRank == 0))
+		  {
+		  std::cout << "Splitting tiles in " << m_TileDirectory << std::endl;
+		  boost::timer t; t.restart();
+		  SplitOTBImage<ImageType>(m_InputImage, m_TileDirectory, m_TileWidth,
+		      m_TileHeight, m_Margin, m_NumberOfFirstIterations);
+		  ShowTime(t);
+		  }
+		MPI_Barrier(MPI_COMM_WORLD);
 		
 		// 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;
+				  << ", margin " << m_Margin  << " niter " << m_NumberOfFirstIterations << std::endl;
 		
 		// Boolean indicating if there are remaining fusions
 		bool isFusion = false;
 		
 		// Run first partial segmentation
+		MPI_Barrier(MPI_COMM_WORLD);
+		boost::timer t; t.restart();
 		auto accumulatedMemory = RunFirstPartialSegmentation<TSegmenter>(m_SpecificParameters,
 																		 m_Threshold,
 																		 m_NumberOfFirstIterations,
@@ -57,12 +65,20 @@ namespace lsgrm
 																		 m_ImageWidth,
 																		 m_ImageHeight,
 																		 m_TemporaryDirectory,
-																		 isFusion);
-		
-		std::cout << "Accumulated memory " << accumulatedMemory << " bytes, there is fusion "<< isFusion << std::endl;
+																		 isFusion,
+																		 m_ProcessRank,
+																		 m_NumberOfProcess);
+
+		// Gathering useful variables
+		GatherUsefulVariables(accumulatedMemory, isFusion, m_ProcessRank, m_NumberOfProcess);
+
+		// Time monitoring
+		if (m_ProcessRank == 0)
+		  ShowTime(t);
 
 		while(accumulatedMemory > m_Memory && isFusion)
 		{
+
 		    isFusion = false;
 			accumulatedMemory = RunPartialSegmentation<TSegmenter>(m_SpecificParameters,
 																   m_Threshold,
@@ -76,12 +92,20 @@ namespace lsgrm
 																   m_ImageWidth,
 																   m_ImageHeight,
 																   m_ImageBands,
-																   isFusion);
-			std::cout << "Accumulated memory " << accumulatedMemory << " bytes, there is fusion "
-					  << isFusion << std::endl;
+																   isFusion,
+																   m_ProcessRank,
+																   m_NumberOfProcess);
+
+			// Gathering useful variables
+			GatherUsefulVariables(accumulatedMemory, isFusion, m_ProcessRank, m_NumberOfProcess);
+
 		}
 
-		if(accumulatedMemory <= m_Memory)
+        // Time monitoring
+        if (m_ProcessRank == 0)
+          ShowTime(t);
+
+		if(accumulatedMemory <= m_Memory && m_ProcessRank == 0)
 		{	
 			// Merge all the graphs
 			MergeAllGraphsAndAchieveSegmentation<TSegmenter>(m_SpecificParameters,
@@ -96,7 +120,12 @@ namespace lsgrm
 															 m_ImageHeight,
 															 m_ImageBands,
 															 isFusion,
-															 m_OutputGraphDirectory);
+															 m_OutputGraphDirectory,
+															 m_ProcessRank,
+															 m_NumberOfProcess);
+
+	        ShowTime(t);
+
 		}
 		else
 		{
diff --git a/Code/lsgrmGraphOperations.h b/Code/lsgrmGraphOperations.h
index d3f6d8c80fa819e0956e2c7a2479ab276033b2b6..4352c89afaf3a5c25bd5a624b0824a95ad604ed3 100644
--- a/Code/lsgrmGraphOperations.h
+++ b/Code/lsgrmGraphOperations.h
@@ -28,7 +28,9 @@ namespace lsgrm
 											  const unsigned int imageHeight,
 											  const unsigned int imageBands,
 											  bool& isFusion,
-											  const std::string& outputGraphDirectory);
+											  const std::string& outputGraphDirectory,
+                                              int myrank,
+                                              int nprocs);
 	
 	template<class TSegmenter>
 		long long unsigned int RunFirstPartialSegmentation(const typename TSegmenter::ParameterType& params,
@@ -44,7 +46,9 @@ namespace lsgrm
 														   const unsigned int imageWidth,
 														   const unsigned int imageHeight,
 														   const std::string& tmpDir,
-														   bool& isFusion);
+														   bool& isFusion,
+														   int myrank,
+														   int nprocs);
 
 	template<class TSegmenter>
 		long long unsigned int RunPartialSegmentation(const typename TSegmenter::ParameterType& params,
@@ -60,7 +64,9 @@ namespace lsgrm
 													  const unsigned int imageWidth,
 													  const unsigned int imageHeight,
 													  const unsigned int imageBands,
-													  bool& isFusion);
+													  bool& isFusion,
+													  int myrank,
+													  int nprocs);
 
 	template<class TSegmenter>
 	void RemoveUselessNodes(ProcessingTile& tile,
diff --git a/Code/lsgrmGraphOperations.txx b/Code/lsgrmGraphOperations.txx
index 2e8085124e4a75cd5365ebebe0f9a53576a1544b..ada54a32f63c49b9883be90ff5f5869c15a34d72 100644
--- a/Code/lsgrmGraphOperations.txx
+++ b/Code/lsgrmGraphOperations.txx
@@ -16,17 +16,26 @@ namespace lsgrm
 											  const unsigned int imageHeight,
 											  const unsigned int imageBands,
 											  bool& isFusion,
-											  const std::string& outputGraphDirectory)
+											  const std::string& outputGraphDirectory,
+											  int myrank,
+											  int nprocs)
 	{
+
+	  // TODO parallelize this ...
+
+	  if (myrank == 0)
+	    {
 		TSegmenter segmenter;
 		std::string nodesPath, edgesPath;
 		
-		std::cout << "Graph aggregation..." << std::endl;
+		if (myrank == 0)
+		  std::cout << "--- Graph aggregation...\n" << std::endl;
+
 		for(unsigned int row = 0; row < nbTilesY; ++row)
 		{
 			for(unsigned int col = 0; col < nbTilesX; col++)
 			{
-				std::cout << "*" << std::flush;
+//				std::cout << "*" << std::flush;
 				typename TSegmenter::GraphType graph;
 
 				// Load the graph
@@ -39,55 +48,63 @@ namespace lsgrm
 												 graph.m_Nodes.end());
 			}
 		}
-		std::cout << "\nRemoving duplicated nodes and updating neighbors..." << std::endl;
+		if (myrank == 0)
+		  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::cout << "*" << std::flush;
+//			  if (MyTurn(myrank,nprocs,row*nbTilesX + col))
+			    {
+			    std::cout << "Cleaning nodes of tile " << (row*nbTilesX + col) << " on " << (nbTilesX*nbTilesY) << std::endl;
 
-				std::unordered_map<long unsigned int,
-								   std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
+			    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);
+			    BuildBorderPixelMap<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], row, col,
+			        nbTilesX, nbTilesY, borderPixelMap, imageWidth);
 
-					
+			    RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);
 
-				
-				RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);
 
-					
-				UpdateNeighborsOfNoneDuplicatedNodes<TSegmenter>(borderPixelMap,
-																 imageWidth,
-																 imageHeight);
+			    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();
+		if (myrank == 0)
+		  {
+		  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>
@@ -103,124 +120,135 @@ namespace lsgrm
 												  const unsigned int imageWidth,
 												  const unsigned int imageHeight,
 												  const unsigned int imageBands,
-												  bool& isFusion)
+												  bool& isFusion,
+												  int myrank,
+												  int nprocs)
 	{
 		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;
+		if (myrank == 0)
+		  std::cout << "--- Running partial segmentations...\nNumber of neighbor layers " << numberOfNeighborLayers << std::endl;
+
 		for(unsigned int row = 0; row < nbTilesY; ++row)
 		{
 			for(unsigned int col = 0; col < nbTilesX; col++)
 			{
 				//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);
+			  if (MyTurn(myrank, nprocs, row*nbTilesX + col))
+			    {
+                  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);
+                  // 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::unordered_map<long unsigned int,
+                                         std::vector<typename TSegmenter::NodePointerType> > borderPixelMap;
 
-					std::cout << "\tRemove duplicated nodes..." << std::endl;
-					RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);
+                      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 << "\tUpdate neighbors.." << std::endl;
-					UpdateNeighborsOfNoneDuplicatedNodes<TSegmenter>(borderPixelMap,
-																	 imageWidth,
-																	 imageHeight);
+                      std::cout << "\tRemove duplicated nodes..." << std::endl;
+                      RemoveDuplicatedNodes<TSegmenter>(borderPixelMap, segmenter.m_Graph, imageWidth);
 
-					std::cout << "\tRemove useless.." << std::endl;
-					RemoveUselessNodes<TSegmenter>(tiles[row*nbTilesX + col], segmenter.m_Graph,
-												   row, col, nbTilesX, nbTilesY, imageWidth, numberOfNeighborLayers);
+                      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);
+                  // 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);
 
-				if(merge == true)
-					isFusion = true;
+                  std::cout << "\tPartial segmentation.." << std::endl;
+                  auto merge = lsrm::GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(segmenter);
 
-				std::cout << "\tRemove unstable segments..." << std::endl;
-				// Remove unstable segments
-				RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], imageWidth);
+                  if(merge == true)
+                      isFusion = true;
 
-				// Retrieve the amount of memory to store this graph
-				accumulatedMemory += GetGraphMemory<TSegmenter>(segmenter.m_Graph);
+                  std::cout << "\tRemove unstable segments..." << std::endl;
+                  // Remove unstable segments
+                  RemoveUnstableSegments<TSegmenter>(segmenter.m_Graph, tiles[row*nbTilesX + col], imageWidth);
 
-				std::cout << "\tWrite graph..." << std::endl;
-				// Write graph to temporay directory (warning specific to Baatz & Schape !!!)
-				WriteGraph<TSegmenter>(segmenter.m_Graph, tmpDir, row, col);
-				
+                  // 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;
+		MPI_Barrier(MPI_COMM_WORLD);
+
+		if (myrank == 0)
+		  std::cout << "Add stability margins to graph for the next round..."<< std::endl;
 
 		// During this step we extract the stability margin for the next round
 		for(unsigned int row = 0; row < nbTilesY; ++row)
 		{
 			for(unsigned int col = 0; col < nbTilesX; col++)
 			{
-				std::cout << "*" << std::flush;
-				typename TSegmenter::GraphType graph;
+			  if (MyTurn(myrank, nprocs, row*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);
+                  // 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;
+                  // 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);
+                      DetectBorderNodes<TSegmenter>(graph, tiles[row*nbTilesX + col],
+                                                    borderNodeMap, imageWidth, imageHeight);
 
 
-					ExtractStabilityMargin<TSegmenter>(borderNodeMap, numberOfNeighborLayers);
+                      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";
+                      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);
-				}
 
+                      WriteStabilityMargin<TSegmenter>(borderNodeMap, nodesPath, edgesPath);
+                  }
+                  std::cout << "Process " << myrank << " finished on tile " << (row*nbTilesX + col) << std::endl;
+			    }
 			}
 		}
 		std::cout << std::endl;
@@ -669,7 +697,9 @@ namespace lsgrm
 													   const unsigned int imageWidth,
 													   const unsigned int imageHeight,
 													   const std::string& tmpDir,
-													   bool& isFusion)
+													   bool& isFusion,
+													   int myrank, // MPI rank
+													   int nprocs) // MPI nb of process
 	{
 		using ImageType = typename TSegmenter::ImageType;
 		using Reader = otb::ImageFileReader<ImageType>;
@@ -678,67 +708,73 @@ namespace lsgrm
 		isFusion = false;
 
 		const unsigned int numberOfNeighborLayers = static_cast<unsigned int>(pow(2, niter2 + 1) - 2);
-		std::cout << "Number of neighbor layers " << numberOfNeighborLayers << std::endl; 
+		if (myrank==0)
+		  std::cout << "--- Running fist partial segmentation...\nNumber of neighbor layers " << numberOfNeighborLayers << std::endl;
 		
 		for(unsigned int row = 0; row < nbTilesY; ++row)
 		{
 			for(unsigned int col = 0; col < nbTilesX; col++)
 			{
-				// 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);
-				}
+			  if (MyTurn(myrank, nprocs, row*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);
+                  }
+
+                  std::cout << "Process " << myrank << " finished on tile " << (row*nbTilesX + col) << std::endl;
+			    }
 			}
 		}
 
diff --git a/Code/lsgrmHeader.h b/Code/lsgrmHeader.h
index e52c5c88b00d51f916eb596d04dbabc09867ca00..0cd91833be05b128131184432178a835e8575523 100644
--- a/Code/lsgrmHeader.h
+++ b/Code/lsgrmHeader.h
@@ -11,4 +11,79 @@
 #include <stack>
 #include <boost/algorithm/string.hpp>
 
+#include <boost/progress.hpp>
+#include "mpi/mpi.h"
+
+#define TAG_SUM     0
+#define TAG_AVERAGE 1
+#define TAG_BORDER  2
+#define TAG_PIECE   3
+
+/*
+ * This function returns TRUE if it's to the process #myrank to do the
+ * work on the yard #div in a pool of #nprocs threads
+ */
+bool MyTurn(int myrank, int nprocs, int div)
+{
+  int proc = 0;
+  if (nprocs != 0)
+    proc = div % nprocs;
+  return (proc == myrank);
+}
+
+/*
+ * This function gather the given value in other process, and update it
+ */
+template<typename T>
+void GatherMe(T& x, MPI_Datatype dataType, int myrank, int nprocs)
+{
+  if (myrank == 0)
+    {
+    // Master process
+    // Gather
+    for (unsigned int p = 1 ; p < nprocs ; p++)
+      {
+      T partial_sum;
+      MPI_Recv( &partial_sum, 1, dataType, p, TAG_PIECE, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+      x += partial_sum;
+      }
+    // Dispatch
+    for (unsigned int p = 1 ; p < nprocs ; p++)
+      MPI_Send(&x, 1, dataType, p, TAG_PIECE, MPI_COMM_WORLD);
+    }
+  else
+    {
+    // Slave process
+    MPI_Send(&x, 1, dataType, 0, TAG_PIECE, MPI_COMM_WORLD);
+    MPI_Recv(&x, 1, dataType, 0, TAG_PIECE, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
+    }
+}
+
+/*
+ * Gather accumulatedMemory and isFusion variables
+ */
+void GatherUsefulVariables(unsigned long long int& accumulatedMemory, bool& isFusion, int myrank, int nprocs)
+{
+  MPI_Barrier(MPI_COMM_WORLD);
+  int isFusionInteger = 0;
+  long long int accumulatedMemoryLLI = static_cast<long long int>(accumulatedMemory);
+  if (isFusion)
+    isFusionInteger = 1;
+  GatherMe<int>(isFusionInteger, MPI_INT, myrank, nprocs);
+  GatherMe<long long int>(accumulatedMemoryLLI, MPI_LONG_LONG_INT, myrank, nprocs);
+  accumulatedMemory = static_cast<long long unsigned int>(accumulatedMemoryLLI);
+  if (isFusionInteger>0)
+    isFusion = true;
+  if (myrank == 0)
+    std::cout << "Accumulated memory " << accumulatedMemory << " bytes, there is fusion "<< isFusion << std::endl;
+}
+
+/*
+ * Print time elapsed
+ */
+void ShowTime(boost::timer t)
+{
+  std::cout << "--- Process duration : " << t.elapsed() << std::endl;
+  t.restart();
+}
 #endif