diff --git a/Applications/CMakeLists.txt b/Applications/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..c6fedad9a49460874979c12a2d89c3f6c47bbf86
--- /dev/null
+++ b/Applications/CMakeLists.txt
@@ -0,0 +1,19 @@
+#=========================================================================
+
+#  Program: 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_executable(RegionMergingSegmentation RegionMergingSegmentation.cxx)
+target_link_libraries(RegionMergingSegmentation OTBGRM)
diff --git a/Applications/RegionMergingSegmentation.cxx b/Applications/RegionMergingSegmentation.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..29ed57bad6bcd8b9f80a0ede71fdc615d4223fb1
--- /dev/null
+++ b/Applications/RegionMergingSegmentation.cxx
@@ -0,0 +1,52 @@
+#include <iostream>
+#include <otbImage.h>
+#include <otbImageFileReader.h>
+#include <otbVectorImage.h>
+#include "lsrmBaatzSegmenter.h"
+
+int main(int argc, char *argv[])
+{
+	if(argc != 7)
+	{
+		std::cerr << "Usage: ./" << argv[0] << "\n"
+				  << "\t[input image path] : (.jpg, .png, .tif)\n"
+				  << "\t[output clustered image] : (.jpg, .png, .tif)\n"
+				  << "\t[output label image] : (.tif)\n"
+				  << "\t[spectral weight] : range from 0 to 1\n"
+				  << "\t[shape weight] : range from 0 to 1\n"
+				  << "\t[scale threshold] : unlimited positive value"
+				  << std::endl;
+		return 1;
+	}
+
+	lsrm::BaatzParam params;
+	const char * input_image = argv[1];
+	const char * clustered_image = argv[2];
+	const char * label_image = argv[3];
+	params.m_SpectralWeight = atof(argv[4]);
+	params.m_ShapeWeight = atof(argv[5]);
+	float sqrt_scale = atof(argv[6]);
+	const float scale = sqrt_scale * sqrt_scale;
+
+	typedef float PixelType;
+	typedef otb::VectorImage<PixelType, 2> ImageType;
+	typedef lsrm::BaatzSegmenter<ImageType> SegmenterType;
+
+	SegmenterType segmenter;
+	segmenter.SetParam(params);
+	segmenter.SetThreshold(scale);
+	segmenter.SetInputFileName(input_image);
+	segmenter.SetClusteredImageFileName(clustered_image);
+	segmenter.SetLabelImageFileName(label_image);
+	
+	segmenter.RunSegmentation();
+	
+    return 0;
+}
+
+
+
+
+
+
+
diff --git a/src/Applications/CMakeLists.txt b/CMakeLists.txt
similarity index 62%
rename from src/Applications/CMakeLists.txt
rename to CMakeLists.txt
index 445d723caab963206780ef65e8f7f3868d4a9b51..8dce6f250ae43eaf9cd460f0236651c8d92c3151 100644
--- a/src/Applications/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,14 +15,19 @@
 #     PURPOSE.  See the above copyright notices for more information.
 
 #=========================================================================
-add_executable(EuclideanDistanceSegmentation EuclideanDistanceSegmentation.cxx)
-target_link_libraries(EuclideanDistanceSegmentation 
-					GRM)
+project(GRM)
 
-add_executable(FLSASegmentation FLSASegmentation.cxx)
-target_link_libraries(FLSASegmentation 
-						GRM)
+cmake_minimum_required(VERSION 2.8)
 
-add_executable(BaatzSegmentation BaatzSegmentation.cxx)
-target_link_libraries(BaatzSegmentation 
-					GRM)
+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)
+
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/Code)
+
+add_subdirectory(Code)
+add_subdirectory(Applications)
diff --git a/src/Library/CMakeLists.txt b/Code/CMakeLists.txt
similarity index 82%
rename from src/Library/CMakeLists.txt
rename to Code/CMakeLists.txt
index dd34c538546abe82961090cf6e4ecfc54b34dc8b..67a7261aaaebadc5b0f2ff980ba8b667d6c5c1a7 100644
--- a/src/Library/CMakeLists.txt
+++ b/Code/CMakeLists.txt
@@ -15,6 +15,8 @@
 #     PURPOSE.  See the above copyright notices for more information.
 
 #=========================================================================
-file(GLOB SRC "*.cxx")
-add_library(GRM ${SRC})
-target_link_libraries(GRM OTBCommon OTBIO)
+add_library(OTBGRM 
+			lsrmContourOperations.cpp
+			lsrmNeighborhood.cpp)
+
+target_link_libraries(OTBGRM OTBCommon OTBIO)
diff --git a/Code/lsrmBaatzSegmenter.h b/Code/lsrmBaatzSegmenter.h
new file mode 100644
index 0000000000000000000000000000000000000000..55dd3a3442a5c5710e25ac75957af18a7b87d022
--- /dev/null
+++ b/Code/lsrmBaatzSegmenter.h
@@ -0,0 +1,102 @@
+#ifndef __LSRM_BAATZ_SEGMENTER_H
+#define __LSRM_BAATZ_SEGMENTER_H
+#include "lsrmSegmenter.h"
+#include "lsrmGraphToOtbImage.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 typename Superclass::GraphType GraphType;
+		typedef typename Superclass::NodePointerType NodePointerType;
+		typedef typename Superclass::GraphOperatorType GraphOperatorType;
+		typedef GraphToOtbImage<GraphType> IOType;
+
+		BaatzSegmenter();
+		
+		void RunSegmentation();
+		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..c91db1f001c50675613b3a056ffe658cc1f84a81
--- /dev/null
+++ b/Code/lsrmBaatzSegmenter.txx
@@ -0,0 +1,159 @@
+#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()
+	{
+		this->m_DoBFSegmentation = true;
+		this->m_NumberOfIterations = 75;
+	}
+	
+	template<class TImage>
+	void
+	BaatzSegmenter<TImage>::InitFromImage()
+	{
+		typedef otb::ImageFileReader<TImage> ImageReaderType;
+		typedef itk::ImageRegionIterator<TImage> ImageIterator;
+
+		typename ImageReaderType::Pointer reader = ImageReaderType::New();
+		reader->SetFileName(this->m_InputFileName);
+		reader->Update();
+
+		this->m_ImageWidth = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
+		this->m_ImageHeight = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[1];
+		this->m_NumberOfComponentsPerPixel = reader->GetOutput()->GetNumberOfComponentsPerPixel();
+
+		std::size_t idx = 0;
+		ImageIterator it(reader->GetOutput(), reader->GetOutput()->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 std::size_t bands = n1->m_Means.size();
+		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 BoundingBox merged_bbox = 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>::RunSegmentation()
+	{
+		GraphOperatorType::InitNodes(this->m_Graph, *this, this->m_InputFileName, FOUR);
+		
+		bool prev_merged =
+			GraphOperatorType::PerfomAllIterationsWithLMBFAndConstThreshold(this->m_Graph, *this,
+																			this->m_Threshold, this->m_NumberOfIterations,
+																			this->m_ImageWidth, this->m_ImageHeight);
+
+		if(prev_merged && this->m_DoBFSegmentation)
+		{
+			prev_merged =
+				GraphOperatorType::PerfomAllIterationsWithBFAndConstThreshold(this->m_Graph, *this,
+																			  this->m_Threshold, this->m_NumberOfIterations,
+																			  this->m_ImageWidth, this->m_ImageHeight);
+		}
+
+		if(!this->m_ClusteredImageFileName.empty())
+			IOType::WriteOutputRGBImage(this->m_Graph,
+										this->m_ImageWidth,
+										this->m_ImageHeight,
+										this->m_ClusteredImageFileName);
+
+		if(!this->m_LabelImageFileName.empty())
+			IOType::WriteLabelImage(this->m_Graph,
+									this->m_ImageWidth,
+									this->m_ImageHeight,
+									this->m_LabelImageFileName);
+
+		if(!this->m_ContourImageFileName.empty())
+			IOType::WriteContourImage(this->m_Graph,
+									  this->m_InputFileName,
+									  this->m_LabelImageFileName);
+	}
+} // end of namespace lsrm
+
+#endif
+
+
+
+
+
+
+
diff --git a/Code/lsrmContourOperations.cpp b/Code/lsrmContourOperations.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0605a712199e0ed3d494a320887703c02acd36c4
--- /dev/null
+++ b/Code/lsrmContourOperations.cpp
@@ -0,0 +1,689 @@
+#include "lsrmContourOperations.h"
+
+namespace lsrm
+{
+	BoundingBox
+	ContourOperations::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;
+	}
+
+	long unsigned int
+	ContourOperations::ToBoundingBoxCoordinates(const long unsigned int id,
+												const BoundingBox& bb,
+												const unsigned int width)
+	{
+		long unsigned int img_x = id % width;
+		long unsigned int img_y = id / width;
+		long unsigned int bb_x = img_x - bb.m_UX;
+		long unsigned int bb_y = img_y - bb.m_UY;
+		return (bb_y * bb.m_W + bb_x);
+	}
+
+	typename ContourOperations::PixelList
+	ContourOperations::GenerateBorderPixels(long unsigned int id,
+											const Contour& contour,
+											const unsigned int width)
+	{
+		PixelList pixels;
+		pixels.push_back(id);
+
+		short curr_move, prev_move;
+		
+		if(contour.size() > 4)
+		{
+			for(ContourConstIterator m = contour.begin()+1; m != contour.end(); m++)
+			{
+				curr_move = (*m)[0] + 2*(*m)[1];
+				prev_move = (*(m-1))[0] + 2*(*(m-1))[1];
+
+				if(curr_move == 0)
+				{
+					if(prev_move == 0)
+					{
+						id-=width;
+						pixels.push_back(id);
+					}
+					else if(prev_move == 1) // 1
+					{
+						id++;
+						pixels.push_back(id);
+						id-=width;
+						pixels.push_back(id);
+					}
+				}
+				else if(curr_move == 1)
+				{
+					if(prev_move == 1)
+					{
+						id++;
+						pixels.push_back(id);
+					}
+					else if(prev_move == 2)
+					{
+						id+=width;
+						pixels.push_back(id);
+						id++;
+						pixels.push_back(id);
+					}
+				}
+				else if(curr_move == 2)
+				{
+					if(prev_move == 3)
+					{
+						id--;
+						pixels.push_back(id);
+						id+=width;
+						pixels.push_back(id);
+					}
+					else if(prev_move == 2)
+					{
+						id+=width;
+						pixels.push_back(id);
+					}
+				}
+				else if(curr_move == 3)
+				{
+					if(prev_move==0)
+					{
+						id-=width;
+						pixels.push_back(id);
+						id--;
+						pixels.push_back(id);
+					}
+					else if(prev_move == 3)
+					{
+						id--;
+						pixels.push_back(id);
+					}
+				}
+			}
+			// Remove duplicated pixels
+			std::sort(pixels.begin(), pixels.end());
+			PixelIterator it = std::unique(pixels.begin(), pixels.end());
+  			pixels.resize(std::distance(pixels.begin(), it));
+		}
+		return pixels;
+	}
+
+	bool
+	ContourOperations::GetCollisionAtNorth(const long unsigned int idx,
+										   const short * grid,
+										   const unsigned int grid_width)
+	{
+		long int tmp_idx = idx - grid_width;
+
+		if(grid[tmp_idx] > 0)
+			return true;
+
+		while(tmp_idx > 0)
+		{
+			if(grid[tmp_idx] > 0)
+				return true;
+			tmp_idx -= grid_width;
+		}
+		return false;
+	}
+
+	bool
+	ContourOperations::GetCollisionAtNorthEast(const long unsigned int idx,
+											   const short * grid,
+											   const unsigned int grid_width)
+	{
+		long int tmp_idx = idx - grid_width + 1;
+
+		if(grid[tmp_idx] > 0)
+			return true;
+
+		while(tmp_idx > 0 && (tmp_idx % grid_width) != 0)
+		{
+			if(grid[tmp_idx] > 0)
+				return true;
+			tmp_idx = tmp_idx - grid_width + 1;
+		}
+		return false;
+	}
+
+	bool
+	ContourOperations::GetCollisionAtEast(const long unsigned int idx,
+										  const short * grid,
+										  const unsigned int grid_width)
+	{
+		long int tmp_idx = idx + 1;
+
+
+		if(grid[tmp_idx] > 0)
+			return true;
+
+		while((tmp_idx % grid_width) != 0)
+		{
+			if(grid[tmp_idx] > 0)
+				return true;
+			++tmp_idx;
+		}
+		return false;
+	}
+
+	bool
+	ContourOperations::GetCollisionAtSouthEast(const long unsigned int idx,
+											   const short * grid,
+											   const unsigned int grid_width,
+											   const unsigned int grid_height)
+	{
+		long int tmp_idx = idx + 1 + grid_width;
+
+		if(grid[tmp_idx] > 0)
+			return true;
+
+		while(tmp_idx < (grid_width * grid_height) &&
+			(tmp_idx % grid_width) != 0)
+		{
+			if(grid[tmp_idx] > 0)
+				return true;
+			tmp_idx = tmp_idx + 1 + grid_width;
+		}
+		return false;
+	}
+
+	bool
+	ContourOperations::GetCollisionAtSouth(const long unsigned int idx,
+										   const short * grid,
+										   const unsigned int grid_width,
+										   const unsigned int grid_height)
+	{
+		long int tmp_idx = idx + grid_width;
+
+		if(grid[tmp_idx] > 0)
+			return true;
+
+		while(tmp_idx < (grid_width * grid_height))
+		{
+			if(grid[tmp_idx] > 0)
+				return true;
+			tmp_idx += grid_width;
+		}
+		return false;
+	}
+
+	bool
+	ContourOperations::GetCollisionAtSouthWest(const long unsigned int idx,
+											   const short * grid,
+											   const unsigned int grid_width,
+											   const unsigned int grid_height)
+	{
+		long int tmp_idx = idx + grid_width - 1;
+
+		if(grid[tmp_idx] > 0)
+			return true;
+
+		while(tmp_idx < (grid_width * grid_height) &&
+			  (tmp_idx % grid_width) != grid_width - 1)
+		{
+			if(grid[tmp_idx] > 0)
+				return true;
+			tmp_idx = tmp_idx + grid_width - 1;
+		}
+		return false;
+	}
+
+	bool
+	ContourOperations::GetCollisionAtWest(const long unsigned int idx,
+										  const short * grid,
+										  const unsigned int grid_width)
+	{
+		long int tmp_idx = idx - 1;
+
+		if(grid[tmp_idx] > 0)
+			return true;
+
+		while((tmp_idx % grid_width) != grid_width - 1)
+		{
+			if(grid[tmp_idx] > 0)
+				return true;
+			tmp_idx -= 1;
+		}
+		return false;
+	}
+
+	bool
+	ContourOperations::GetCollisionAtNorthWest(const long unsigned int idx,
+											   const short * grid,
+											   const unsigned int grid_width)
+	{
+		long int tmp_idx = idx - grid_width - 1;
+
+		if(grid[tmp_idx] > 0)
+			return true;
+
+		while(tmp_idx > 0 && (tmp_idx % grid_width) != grid_width - 1)
+		{
+			if(grid[tmp_idx] > 0)
+				return true;
+			tmp_idx = tmp_idx - grid_width - 1;
+		}
+		return false;
+	}
+
+	bool
+	ContourOperations::IsInternalPixel(const long unsigned int idx,
+									   const short * grid,
+									   const unsigned int grid_width,
+									   const unsigned int grid_height)
+	{
+		if(!GetCollisionAtNorth(idx, grid, grid_width))
+			return false;
+		if(!GetCollisionAtNorthEast(idx, grid, grid_width))
+			return false;
+		if(!GetCollisionAtEast(idx, grid, grid_width))
+			return false;
+		if(!GetCollisionAtSouthEast(idx, grid, grid_width, grid_height))
+			return false;
+		if(!GetCollisionAtSouth(idx, grid, grid_width, grid_height))
+			return false;
+		if(!GetCollisionAtSouthWest(idx, grid, grid_width, grid_height))
+			return false;
+		if(!GetCollisionAtWest(idx, grid, grid_width))
+			return false;
+		if(!GetCollisionAtNorthWest(idx, grid, grid_width))
+			return false;
+
+		return true;
+	}
+
+	void
+	ContourOperations::UpdateSight(short * sight, const short direction)
+	{
+		sight[0] = (direction + 3) % 4; // look at the left
+		sight[1] = direction; // look in front
+		sight[2] = (direction + 1) % 4; // look at the right
+		sight[3] = (direction + 2) % 4; // look behind
+	}
+
+	void
+	ContourOperations::EncodeContourForTopPixel(long unsigned int& curr_mat_id,
+												short& direction,
+												short& pos,
+												Contour& curr_contour,
+												const unsigned int grid_width)
+	{
+		if(pos == 1)
+		{
+			if(direction==0)
+				pos = 4;
+			else if(direction == 1)
+			{
+				pos = 1;
+				curr_contour.push_back(0);
+			}
+		}
+		else if(pos == 2)
+		{
+			pos = 4;
+			curr_contour.push_back(2);
+			curr_contour.push_back(3);
+			curr_contour.push_back(0);
+		}
+		else if(pos == 3)
+		{
+			pos = 4;
+			curr_contour.push_back(3);
+			curr_contour.push_back(0);
+		}
+		else if(pos == 4)
+			curr_contour.push_back(0);
+
+
+		curr_mat_id -= grid_width;
+		direction = 0;
+	}
+
+	void
+	ContourOperations::EncodeContourForRightPixel(long unsigned int& curr_mat_id,
+												  short& direction,
+												  short& pos,
+												  Contour& curr_contour)
+	{
+		if(pos == 1)
+			curr_contour.push_back(1);
+		else if(pos == 2)
+		{
+			if(direction == 1)
+				pos = 1;
+			else if(direction == 2)
+				curr_contour.push_back(1);
+		}
+		else if(pos == 3)
+		{
+			pos = 1;
+			curr_contour.push_back(3);
+			curr_contour.push_back(0);
+			curr_contour.push_back(1);
+		}
+		else if(pos == 4)
+		{
+			pos = 1;
+			curr_contour.push_back(0);
+			curr_contour.push_back(1);
+		}
+
+		curr_mat_id++;
+		direction = 1;
+	}
+
+	void
+	ContourOperations::EncodeContourForBottomPixel(long unsigned int& curr_mat_id,
+												   short& direction,
+												   short& pos,
+												   Contour& curr_contour,
+												   const unsigned int grid_width)
+	{
+		if(pos == 1)
+		{
+			pos = 2;
+			curr_contour.push_back(1);
+			curr_contour.push_back(2);
+		}
+		else if(pos == 2)
+			curr_contour.push_back(2);
+		else if(pos == 3)
+		{
+			if(direction == 2)
+				pos = 2;
+			else if(direction == 3)
+				curr_contour.push_back(2);
+		}
+		else if(pos == 4)
+		{
+			pos = 2;
+			curr_contour.push_back(0);
+			curr_contour.push_back(1);
+			curr_contour.push_back(2);
+		}
+
+		curr_mat_id += grid_width;
+		direction = 2;	
+	}
+
+	void
+	ContourOperations::EncodeContourForLeftPixel(long unsigned int& curr_mat_id,
+												 short& direction,
+												 short& pos,
+												 Contour& curr_contour)
+	{
+		if(pos == 1)
+		{
+			pos = 3;
+			curr_contour.push_back(1);
+			curr_contour.push_back(2);
+			curr_contour.push_back(3);
+		}
+		else if(pos == 2)
+		{
+			pos = 3;
+			curr_contour.push_back(2);
+			curr_contour.push_back(3);
+		}
+		else if(pos == 3)
+			curr_contour.push_back(3);
+		else if(pos == 4)
+		{
+			if(direction == 3)
+				pos = 3;
+			else if(direction == 0)
+			{
+				pos = 4;
+				curr_contour.push_back(3);
+			}
+		}
+
+		curr_mat_id --;
+		direction = 3;
+	}
+
+	Contour
+	ContourOperations::CreateContourFromBorderPixels(const long unsigned int id,
+													 const BoundingBox& bbox,
+													 const short * grid,
+													 const unsigned int width)
+	{
+		/* Location on one of the 4 corners of a pixel */
+		short direction, pos;
+		/* Field of view */
+		short sight[4];
+
+		long unsigned int start_mat_id = ToBoundingBoxCoordinates(id, bbox, width);
+		long unsigned int curr_mat_id = start_mat_id;
+
+		/* New contour */
+		Contour curr_contour;
+
+		// Build the first move
+		long int neighbors[4];
+		FOURNeighborhood(neighbors, curr_mat_id, bbox.m_W, bbox.m_H);
+
+		// Only 2 cases: neighbor at the right or at the bottom
+		if(grid[neighbors[1]] > 0) // right
+		{
+			pos = 1;       // top left corner of the pixel
+			curr_mat_id++; // go to the pixel at the right
+			direction = 1; // direction is along the right
+			curr_contour.push_back(1); // Add the move to the right
+		}
+		else if(grid[neighbors[2]] > 0) // bottom
+		{
+			pos = 2; // top right corner of the pixel
+			curr_mat_id += bbox.m_W; // go to the pixel at the bottom
+			direction = 2; // direction is along the bottom						
+			curr_contour.push_back(1); // add move to the right
+			curr_contour.push_back(2); // add move to the bottom	
+		}
+
+		// Keep going this same reasonning until we reach the start pixel (start_mat_id)
+		while(curr_mat_id != start_mat_id)
+		{
+			UpdateSight(sight, direction);
+			FOURNeighborhood(neighbors, curr_mat_id, bbox.m_W, bbox.m_H);
+
+			// All the cases are possible: top, right, bottom, left
+			for(short d=0; d<4; d++)
+			{
+				if(neighbors[sight[d]] > -1)
+				{
+					if(grid[neighbors[sight[d]]] > 0)
+					{
+						if(sight[d] == 0)	
+						{
+							EncodeContourForTopPixel(curr_mat_id, direction, pos, curr_contour, bbox.m_W);
+							break;
+						}
+						else if(sight[d] == 1)
+						{
+							EncodeContourForRightPixel(curr_mat_id, direction, pos, curr_contour);
+							break;
+						}
+						else if(sight[d] == 2)
+						{
+							EncodeContourForBottomPixel(curr_mat_id, direction, pos, curr_contour, bbox.m_W);
+							break;
+						}
+						else if(sight[d] == 3)
+						{
+							EncodeContourForLeftPixel(curr_mat_id, direction, pos, curr_contour);
+							break;
+						}
+					}
+				}
+			}
+
+			// It is possible to reach the start pixel whithout finishing to encode the whole contour.
+			// 00111
+			// 11100
+			// In this case the direction value is always 3, the position value is always 3
+			// We have to check if there is a neighbor at left of the current direction
+			if(curr_mat_id == start_mat_id)
+			{
+				if(pos == 3 && direction == 3)
+				{
+					UpdateSight(sight, direction);
+					FOURNeighborhood(neighbors, curr_mat_id, bbox.m_W, bbox.m_H);
+
+					if(neighbors[sight[0]] > -1)
+					{
+						if(grid[neighbors[sight[0]]] == 1)
+						{
+							pos = 3;
+							direction = 2;
+							curr_contour.push_back(2);
+							curr_mat_id += bbox.m_W;
+						}
+					}
+				}
+			}		
+		}
+
+		// We reach the start pixel but maybe we did not encode the whole contour
+		if(pos == 3)
+		{
+			curr_contour.push_back(3);
+			curr_contour.push_back(0);
+		}
+		else if(pos == 4)
+			curr_contour.push_back(0);
+
+		return curr_contour;
+	}
+	
+	Contour
+	ContourOperations::MergeContours(const long unsigned int start_id1,
+									 const long unsigned int start_id2,
+									 const Contour& c1,
+									 const Contour& c2,
+									 const BoundingBox& bb,
+									 const unsigned int width,
+									 const unsigned int height)
+	{
+		/* Grid with the dimension of the bounding box initialized to 0 */
+		short matrix_bbox[bb.m_W * bb.m_H];
+		memset(matrix_bbox, 0, bb.m_W * bb.m_H * sizeof(short));
+		
+		/* Generate the border pixels wrt to c1 and c2 */
+		PixelList pixels1 = GenerateBorderPixels(start_id1, c1, width);
+		{
+			PixelList pixels2 = GenerateBorderPixels(start_id2, c2, width);
+			pixels1.insert(pixels1.end(), pixels2.begin(), pixels2.end());
+		}
+
+		/* Each case of the grid where a pixel is located is set to 1 */
+		for(PixelConstIterator pit = pixels1.begin(); pit != pixels1.end(); ++pit)
+			matrix_bbox[ToBoundingBoxCoordinates(*pit, bb, width)] = 1;
+
+		/* Remove internal pixels */
+		long unsigned int bb_id;
+		bool is_internal;
+		long int neighbors[8];
+		for(long unsigned int y = 1; y < bb.m_H - 1; ++y)
+		{
+			for(long unsigned int x = 1; x < bb.m_W - 1; ++x)
+			{
+				bb_id = y*bb.m_W + x;
+				is_internal = true;
+				if(matrix_bbox[bb_id] > 0)
+				{
+					EIGHTNeighborhood(neighbors, bb_id, bb.m_W, bb.m_H);
+					for(short j = 0; j < 8; ++j)
+					{
+						if(matrix_bbox[neighbors[j]] < 1)
+						{
+							is_internal = false;
+							break;
+						}
+					}
+
+					if(is_internal)
+						matrix_bbox[bb_id] = 2;
+				}
+			}
+		}
+
+		for(long unsigned int i = 0; i < bb.m_H * bb.m_W; ++i)
+		{
+			if(matrix_bbox[i] > 1)
+				matrix_bbox[i] = 0;
+		}
+
+		/* Create a new contour */
+		return CreateContourFromBorderPixels(start_id1, bb, matrix_bbox, width);
+		
+	}
+
+	long unsigned int
+	ContourOperations::ToImageCoordinates(const long unsigned int id,
+										  const BoundingBox& bbox,
+										  const unsigned int width)
+	{
+		long unsigned int bbox_y = id / bbox.m_W;
+		long unsigned int bbox_x = id % bbox.m_W;
+		long unsigned int img_x = bbox.m_UX + bbox_x;
+		long unsigned int img_y = bbox.m_UY + bbox_y;
+		return (img_y * width + img_x);
+	}
+
+	typename ContourOperations::PixelList
+	ContourOperations::GeneratePixels(long unsigned int id,
+									  const Contour& contour,
+									  const BoundingBox& bbox,
+									  const unsigned int width)
+	{
+		// Generate the bounding box grid
+		short matrix_bbox[bbox.m_W * bbox.m_H];
+		memset(matrix_bbox,  0, bbox.m_W * bbox.m_H * sizeof(short));
+
+		PixelList pixels;
+
+		/* Each case of the grid where a pixel is located is set to 1 */
+		{
+			// Generate the border pixels
+			PixelList border_pixels = GenerateBorderPixels(id, contour, width);
+			
+			for(PixelConstIterator pit = border_pixels.begin();
+				pit != border_pixels.end(); ++pit)
+				matrix_bbox[ToBoundingBoxCoordinates(*pit, bbox, width)] = 1;
+
+			pixels.insert(pixels.end(), border_pixels.begin(), border_pixels.end());
+		}
+
+		/* For each case of the grid we determine if it is inside or outside
+		   of the region */
+		long unsigned int bb_id;
+		for(long unsigned int y = 1; y < bbox.m_H - 1; ++y)
+		{
+			for(long unsigned int x = 1; x < bbox.m_W - 1; ++x)
+			{
+				bb_id = y*bbox.m_W + x;
+				if(matrix_bbox[bb_id] < 1)
+				{
+					if(IsInternalPixel(bb_id, matrix_bbox, bbox.m_W, bbox.m_H))
+					{
+						matrix_bbox[bb_id] = 1;
+						pixels.push_back(ToImageCoordinates(bb_id, bbox, width));
+					}
+				}
+			}
+		}
+		return pixels;
+	}
+} // end of namespace lsrm
diff --git a/Code/lsrmContourOperations.h b/Code/lsrmContourOperations.h
new file mode 100644
index 0000000000000000000000000000000000000000..d70ee54efb1b260c30b07b03e6efa3e5f7a8e8b7
--- /dev/null
+++ b/Code/lsrmContourOperations.h
@@ -0,0 +1,305 @@
+#ifndef __LSRM_CONTOUR_OPERATIONS_H
+#define __LSRM_CONTOUR_OPERATIONS_H
+#include "lsrmDataStructures.h"
+#include "lsrmNeighborhood.h"
+#include <cstring>
+#include <algorithm>
+
+namespace lsrm
+{
+	class ContourOperations
+	{
+	public:
+
+		/* Some convenient typedefs */
+		typedef std::vector<long unsigned int> PixelList;
+		typedef typename PixelList::iterator PixelIterator;
+		typedef typename PixelList::const_iterator PixelConstIterator;
+
+		/*
+		 * Given two rectangular bounding boxes,
+		 * it returns the bounding box which is
+		 * the union of both bounding boxes.
+		 *
+		 * @params
+		 * const BoundingBox& bb1 : reference to the first bounding box
+		 * const BoundingBox& bb2 : reference to the second bounding box
+		 * @return the union of bb1 and bb2.
+		 */
+		static BoundingBox MergeBoundingBoxes(const BoundingBox& bb1,
+											  const BoundingBox& bb2);
+
+		/*
+		 * Given the coordinates of the pixel in the current image, the
+		 * bounding box containing it and the width of the image, it returns
+		 * the coordinates of the pixel in the referential of the bounding box.
+		 *
+		 * @params
+		 * const long unsigned int id : coordinates of the pixel in the image
+		 *                              referential.
+		 * const BoundingBox& bb : bounding box containing the pixel
+		 * const unsigned int width : width of the image
+		 * @return the coordinates of the pixel in the bounding box referential.
+		 */
+		static long unsigned int ToBoundingBoxCoordinates(const long unsigned int id,
+												   const BoundingBox& bb,
+												   const unsigned int width);
+
+
+		/*
+		 * Given the contour, the first pixel coordinates
+		 * and the width of the current image,
+		 * it returns the list of the border pixels.
+		 *
+		 * @params
+		 * const long unsigned int id : coordinates of the first pixel.
+		 * const Contour& contour : reference to the contour.
+		 * const unsigned int width : width of the current image.
+		 */
+		static PixelList GenerateBorderPixels(long unsigned int id,
+											  const Contour& contour,
+											  const unsigned int width);
+
+		/*
+		 * Return true if it exists a pixel in the grid at the north
+		 * of the current pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 */
+		static bool GetCollisionAtNorth(const long unsigned int idx,
+										const short * grid,
+										const unsigned int grid_width);
+
+		/*
+		 * Return true if it exists a pixel in the grid at the north East
+		 * of the current pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 */
+		static bool GetCollisionAtNorthEast(const long unsigned int idx,
+											const short * grid,
+											const unsigned int grid_width);
+
+		/*
+		 * Return true if it exists a pixel in the grid at the East
+		 * of the current pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 */
+		static bool GetCollisionAtEast(const long unsigned int idx,
+									   const short * grid,
+									   const unsigned int grid_width);
+
+		/*
+		 * Return true if it exists a pixel in the grid at the South East
+		 * of the current pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 * const unsigned int grid_height : height of the grid
+		 */
+		static bool GetCollisionAtSouthEast(const long unsigned int idx,
+											const short * grid,
+											const unsigned int grid_width,
+											const unsigned int grid_height);
+
+		/*
+		 * Return true if it exists a pixel in the grid at the South
+		 * of the current pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 * const unsigned int grid_height : height of the grid
+		 */
+		static bool GetCollisionAtSouth(const long unsigned int idx,
+										const short * grid,
+										const unsigned int grid_width,
+										const unsigned int grid_height);
+
+		/*
+		 * Return true if it exists a pixel in the grid at the South West
+		 * of the current pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 * const unsigned int grid_height : height of the grid
+		 */
+		static bool GetCollisionAtSouthWest(const long unsigned int idx,
+											const short * grid,
+											const unsigned int grid_width,
+											const unsigned int grid_height);
+
+		/*
+		 * Return true if it exists a pixel in the grid at the West
+		 * of the current pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 */
+		static bool GetCollisionAtWest(const long unsigned int idx,
+									   const short * grid,
+									   const unsigned int grid_width);
+
+		/*
+		 * Return true if it exists a pixel in the grid at the North West
+		 * of the current pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 */
+		static bool GetCollisionAtNorthWest(const long unsigned int idx,
+											const short * grid,
+											const unsigned int grid_width);
+
+		/*
+		 * Return true if the pixel located at the coordinates idx
+		 * is an internal pixel, i.e, in each direction there is a pixel.
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel
+		 * const short * grid : grid
+		 * const unsigned int grid_width : width of the grid
+		 * const unsigned int grid_height : height of the grid
+		 */
+		static bool IsInternalPixel(const long unsigned int idx,
+									const short * grid,
+									const unsigned int grid_width,
+									const unsigned int grid_height);
+
+		/*
+		 * Given the direction of the front, it returns the field of view
+		 *
+		 * @params
+		 * short * sight : former field of view.
+		 * const short direction : new front direction.
+		 */
+		static void UpdateSight(short * sight, const short direction);
+
+
+		/*
+		 * Create a contour element when we move toward the top, right, bottom
+		 * or left wrt to the current pixel.
+		 *
+		 * @params
+		 * long unsigned int& curr_mat_id : pixel coordinates
+		 * short& direction : toward the direction we look
+		 * short& pos : location on the pixel
+		 * Contour& curr_contour : new contour to be created
+		 * const unsigned int grid_width : width of the grid
+		 */
+		static void EncodeContourForTopPixel(long unsigned int& curr_mat_id,
+											 short& direction,
+											 short& pos,
+											 Contour& curr_contour,
+											 const unsigned int grid_width);
+
+		static void EncodeContourForRightPixel(long unsigned int& curr_mat_id,
+											   short& direction,
+											   short& pos,
+											   Contour& curr_contour);
+
+		static void EncodeContourForBottomPixel(long unsigned int& curr_mat_id,
+												short& direction,
+												short& pos,
+												Contour& curr_contour,
+												const unsigned int grid_width);
+
+		static void EncodeContourForLeftPixel(long unsigned int& curr_mat_id,
+											  short& direction,
+											  short& pos,
+											  Contour& curr_contour);
+		
+		/*
+		 * Given the coordinates of the first pixel, the bounding box, the grid
+		 * with the location of the border pixels and the dimension of the referential
+		 * image, it returns a new contour
+		 *
+		 * @params
+		 * const long unsigned int idx : coordinates of the current pixel.
+		 * const BoundingBox& bbox : reference to the bounding box
+		 * const short * grid : grid
+		 * const unsigned int width : width of the current image.
+		 * cont unsigned int height : height of the current image.
+		 */
+		static Contour CreateContourFromBorderPixels(const long unsigned int id,
+													 const BoundingBox& bbox,
+													 const short * grid,
+													 const unsigned int width);
+		
+
+		/*
+		 * Given two contours, the coordinates of their first pixels
+		 * and the bounding box wrapping both contours, it returns the
+		 * new contour by merging c1 and c2.
+		 *
+		 * @params
+		 * const long unsigned int start_id1 : coordinates of the first
+		 *                                     pixel of the first contour
+		 * const long unsigned int start_id2 : coordinates of the first
+		 *                                     pixel of the second contour
+		 * const Contour& c1 : reference to the first contour
+		 * const Contour& c2 : reference to the second contour
+		 * const BoundingBox& bb : reference to the bounding box wrapping c1 and c2
+		 * const unsigned int width : width of the current image
+		 * const unsigned int height :  height of the current image
+		 */
+		static Contour MergeContours(const long unsigned int start_id1,
+									 const long unsigned int start_id2,
+									 const Contour& c1,
+									 const Contour& c2,
+									 const BoundingBox& bb,
+									 const unsigned int width,
+									 const unsigned int height);
+
+		/*
+		 * Given the coordinates of the pixel within its bounding box, it
+		 * returns the coordinates of this pixel in the current image.
+		 *
+		 * @params
+		 * const long unsigned int id : coordinates of the pixel within the bounding box
+		 * const BoundingBox& bbox : reference to the bounding box
+		 * const unsigned int width : width of the image
+		 * @return the coordinates of the pixel in the image. 
+		 */
+		static long unsigned int ToImageCoordinates(const long unsigned int id,
+													const BoundingBox& bbox,
+													const unsigned int width);
+
+		/*
+		 * Given the contour, the first pixel coordinates
+		 * and the width of the current image,
+		 * it returns the list of the pixels within the area
+		 * delineated by the contour.
+		 *
+		 * @params
+		 * const long unsigned int id : coordinates of the first pixel.
+		 * const Contour& contour : reference to the contour.
+		 * const unsigned int width : width of the current image.
+		 */
+		static PixelList GeneratePixels(long unsigned int id,
+										const Contour& contour,
+										const BoundingBox& bbox,
+										const unsigned int width);
+	};
+} // end of namespace lsrm
+
+#endif
diff --git a/Code/lsrmDataStructures.h b/Code/lsrmDataStructures.h
new file mode 100644
index 0000000000000000000000000000000000000000..092e0cc02bd2c3de2f51d645683c0dd85718ccc9
--- /dev/null
+++ b/Code/lsrmDataStructures.h
@@ -0,0 +1,30 @@
+#ifndef __LSRM_DATA_STRUCTURES_H
+#define __LSRM_DATA_STRUCTURES_H
+#include <bitset>
+#include <vector>
+
+namespace lsrm
+{	
+	typedef std::bitset<2> ContourElem;
+	typedef std::vector<ContourElem> Contour;
+	typedef typename Contour::iterator ContourIterator;
+	typedef typename Contour::const_iterator ContourConstIterator;
+
+	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;	
+	};
+	
+} // end of namespace lsrm
+
+#endif
diff --git a/Code/lsrmGraph.h b/Code/lsrmGraph.h
new file mode 100644
index 0000000000000000000000000000000000000000..d20486da0f19b488760a8ca2f3b110d175305268
--- /dev/null
+++ b/Code/lsrmGraph.h
@@ -0,0 +1,84 @@
+#ifndef __LSRM_GRAPH_H
+#define __LSRM_GRAPH_H
+#include <boost/shared_ptr.hpp>
+#include "lsrmDataStructures.h"
+
+namespace lsrm
+{
+	struct BaseNode
+	{
+		BaseNode(){}
+		
+		virtual ~BaseNode() {}
+
+		/* Node already merged. */
+		bool m_Valid;
+		
+		/* Node has to be removed from the graph. */
+		bool m_Expired;
+
+		/* 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.
+		 */
+		BoundingBox m_Bbox;
+
+		/*
+		  Contour of the region
+		 */
+		Contour m_Contour;
+	};
+
+	template<class DerivedNode>
+	struct Node : BaseNode
+	{
+		struct Edge
+		{
+			/* Boundary length between two adjacent regions.*/
+			unsigned int m_Boundary;
+			
+			/* Fusion cost (similarity measure) with the target node. */
+			float m_Cost;
+			
+			/* Pointer to a neighboring node. */
+			boost::shared_ptr<DerivedNode> m_Target;
+		};
+		
+		Node(){};
+		virtual ~Node() {}
+
+		std::vector<Edge> m_Edges;
+	};
+
+	template<class TNode>
+	struct Graph
+	{
+		typedef TNode NodeType;
+		typedef typename NodeType::Edge EdgeType;
+		typedef boost::shared_ptr<NodeType> NodePointerType;
+		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< boost::shared_ptr<TNode> > m_Nodes; 
+	};
+	
+} // end of namespace lsrm
+#endif
+
+
diff --git a/Code/lsrmGraphOperations.h b/Code/lsrmGraphOperations.h
new file mode 100644
index 0000000000000000000000000000000000000000..7c616ade356aef4d1689cbc421f8fb84cbfd8918
--- /dev/null
+++ b/Code/lsrmGraphOperations.h
@@ -0,0 +1,397 @@
+#ifndef __LSRM_GRAPH_OPERATIONS_H
+#define __LSRM_GRAPH_OPERATIONS_H
+#include "lsrmGraph.h"
+#include "lsrmNeighborhood.h"
+#include "lsrmContourOperations.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::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;
+
+		/* Some static constants */
+		static const long long unsigned int NodeSize = sizeof(NodeType) + sizeof(NodePointerType);
+		static const long long unsigned int EdgeSize = sizeof(EdgeType);
+
+		/* Functors (former than lambda expression...
+		   but it can accept all compiler' versions) */
+		struct EdgeEqualTo
+		{
+			explicit EdgeEqualTo(NodePointerType t_) : m_T(t_)
+				{assert(m_T != NULL);}
+			
+			bool operator()(const EdgeType& e)
+				{
+					return e.m_Target == m_T;
+				}
+
+			NodePointerType  m_T;
+		};
+
+		struct IsExpired
+		{
+			IsExpired () {}
+			bool operator()(NodePointerType  n) const
+				{
+					return n->m_Expired;
+				}
+		};
+
+		struct IsInSet
+		{
+			explicit IsInSet(std::set<long unsigned int>& ur) : m_UsedRegions(ur)
+				{
+					
+				}
+
+			bool operator()(const EdgeType& e)
+				{
+					if(m_UsedRegions.find(e.m_Target->m_Id) == m_UsedRegions.end())
+						return true;
+					else
+						return false;
+				}
+				
+			std::set<long unsigned int>& m_UsedRegions;
+		};
+
+		struct NodePtrComparator : public std::unary_function<NodePointerType, bool>
+		{
+			NodePointerType nn;
+			explicit NodePtrComparator(NodePointerType n) : nn(n){}
+			bool operator()(NodePointerType v) { return v == nn;}
+		};
+		
+
+		/*
+		 * 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(GraphType& graph,
+							  SegmenterType& seg,
+							  const std::string& inputFileName,
+							  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(GraphType& graph, 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(GraphType& graph,
+											   SegmenterType& seg,
+											   const float threshold,
+											   const unsigned int width,
+											   const unsigned int height);
+
+		/*
+		 * 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(GraphType& graph,
+											 SegmenterType& seg,
+											 const float threshold,
+											 const unsigned int width,
+											 const unsigned int height);
+
+		/*
+		 * 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(GraphType& graph,
+																 SegmenterType& seg,
+																 const float threshold,
+																 const unsigned int numberOfIterations,
+																 const unsigned int width,
+																 const unsigned int height);
+
+		/*
+		 * 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(GraphType& graph,
+															   SegmenterType& seg,
+															   const float threshold,
+															   const unsigned int numberOfIterations,
+															   const unsigned int width,
+															   const unsigned int height);
+
+		/*
+		 *
+		 */
+		static void LargeScaleInitNodes(GraphType& graph,
+										SegmenterType& seg,
+										const unsigned int tileWidth,
+										const unsigned int tileHeight,
+										const unsigned int nbTilesX,
+										const unsigned int nbTilesY,
+										const unsigned int tileId,
+										const unsigned int margin,
+										const unsigned int imageWidth,
+										const unsigned int imageHeight,
+										CONNECTIVITY mask);
+
+		/*
+		 *
+		 */
+		static bool IsOnTileBorder(const unsigned int pixelRow,
+								   const unsigned int pixelCol,
+								   const unsigned int lowerRow,
+								   const unsigned int upperRow,
+								   const unsigned int lowerCol,
+								   const unsigned int upperCol,
+								   const unsigned int imageWidth,
+								   const unsigned int imageHeight);
+
+		static bool IsOnTileBorderAndAdjacentTileBorder(const unsigned int pixelRow,
+														const unsigned int pixelCol,
+														const unsigned int lowerRow,
+														const unsigned int upperRow,
+														const unsigned int lowerCol,
+														const unsigned int upperCol,
+														const unsigned int imageWidth,
+														const unsigned int imageHeight);
+		/*
+		 *
+		 */
+		static bool IsBboxStrictlyInsideTile(NodePointerType n,
+											 const unsigned int lowerRow,
+											 const unsigned int upperRow,
+											 const unsigned int lowerCol,
+											 const unsigned int upperCol);
+
+		/*
+		 *
+		 */
+		static bool IsBboxInsideTile(NodePointerType n,
+									 const unsigned int lowerRow,
+									 const unsigned int upperRow,
+									 const unsigned int lowerCol,
+									 const unsigned int upperCol);
+		
+		/*
+		 *
+		 */
+		static void RemoveEdgesToUnstableRegion(NodePointerType r);
+		
+		/*
+		 *
+		 */
+		static void RemoveUnstableRegions(GraphType& graph,
+										  unsigned int lowerRow,
+										  unsigned int upperRow,
+										  unsigned int lowerCol,
+										  unsigned int upperCol,
+										  unsigned int imageWidth);
+
+		/*
+		 *
+		 */
+		static void AddNeighboringLayer(NodePointerType n,
+										GraphType& g,
+										std::set<long unsigned int>& bset);
+		
+		/*
+		 *
+		 */
+		static void BuildStabilityMargins(const GraphType& graph,
+										  GraphType& subgraph,
+										  const unsigned int lowerRow,
+										  const unsigned int upperRow,
+										  const unsigned int lowerCol,
+										  const unsigned int upperCol,
+										  const unsigned int imageWidth,
+										  const unsigned int imageHeight);
+
+		/*
+		 *
+		 */
+		static void MergeGraphs(GraphType& graph,
+								const unsigned int lowerRow,
+								const unsigned int upperRow,
+								const unsigned int lowerCol,
+								const unsigned int upperCol,
+								const unsigned int imageWidth,
+								const unsigned int imageHeight,
+								bool useless);
+
+		/*
+		 *
+		 */
+		static long long unsigned int GetMemorySpace(const GraphType& graph);
+	};
+} // end of namespace lsrm
+
+#include "lsrmGraphOperations.txx"
+#endif
+
+
+
+
diff --git a/Code/lsrmGraphOperations.txx b/Code/lsrmGraphOperations.txx
new file mode 100644
index 0000000000000000000000000000000000000000..d604e16c428bbf6b5b7872b45eeac56222021779
--- /dev/null
+++ b/Code/lsrmGraphOperations.txx
@@ -0,0 +1,1144 @@
+#ifndef __LSRM_GRAPH_OPERATIONS_TXX
+#define __LSRM_GRAPH_OPERATIONS_TXX
+#include <otbImageFileReader.h>
+
+namespace lsrm
+{
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::InitNodes(GraphType& graph,
+												SegmenterType& seg,
+												const std::string& inputFileName,
+												CONNECTIVITY mask)
+	{
+		unsigned int width, height;
+		
+		{
+			typedef typename TSegmenter::ImageType ImageType;
+			typedef otb::ImageFileReader<ImageType> ReaderType;
+			typename ReaderType::Pointer reader = ReaderType::New();
+			reader->SetFileName(inputFileName);
+			reader->UpdateOutputInformation();
+
+			width = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[0];
+			height = reader->GetOutput()->GetLargestPossibleRegion().GetSize()[1];
+		}
+		
+		const long unsigned int num_nodes = width * height;
+
+		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_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;
+			n->m_Contour.push_back(1);
+			n->m_Contour.push_back(2);
+			n->m_Contour.push_back(3);
+			n->m_Contour.push_back(0);
+			graph.m_Nodes.push_back(n);
+		}
+
+		if(mask == FOUR)
+		{
+			for(long unsigned int i = 0;
+				i < num_nodes;
+				++i)
+			{
+				long int neighborhood[4];
+				FOURNeighborhood(neighborhood, i, width, height);
+				for(short j = 0; j < 4; ++j)
+				{
+					if(neighborhood[j] > -1)
+					{
+						EdgeType edge;
+						edge.m_Boundary = 1;
+						edge.m_Cost = 0.0f;
+						edge.m_Target = graph.m_Nodes[neighborhood[j]];
+						graph.m_Nodes[i]->m_Edges.push_back(edge);
+					}
+				}
+			}
+		}
+		else
+		{
+			for(long unsigned int i = 0;
+				i < num_nodes;
+				++i)
+			{
+				long int neighborhood[8];
+				EIGHTNeighborhood(neighborhood, i, width, height);
+				for(short j = 0; j < 8; ++j)
+				{
+					if(neighborhood[j] > -1)
+					{
+						EdgeType edge;
+						if(j % 2 > 0)
+							edge.m_Boundary = 0;
+						else
+							edge.m_Boundary = 1;
+						edge.m_Cost = 0.0f;
+						edge.m_Target = graph.m_Nodes[neighborhood[j]];
+						graph.m_Nodes[i]->m_Edges.push_back(edge);
+					}
+				}
+			}
+		}
+		seg.InitFromImage();
+	}
+
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::UpdateMergingCosts(GraphType& graph,
+														 SegmenterType& seg)
+	{		
+		float min_cost;
+		long unsigned int min_id;
+		std::size_t idx, min_idx;
+		
+		for (NodeIterator it = graph.m_Nodes.begin();
+			 it != graph.m_Nodes.end(); ++it)
+		{
+			// Reset boolean flags
+			(*it)->m_Valid = true;
+			(*it)->m_Expired = false;
+
+			min_cost = std::numeric_limits<float>::max();
+			idx = 0;
+			min_idx = 0;
+			
+			for (EdgeIterator et = (*it)->m_Edges.begin();
+				 et != (*it)->m_Edges.end(); ++et)
+			{
+
+				et->m_Cost = seg.ComputeMergingCost(*it, et->m_Target);
+
+				if(min_cost > et->m_Cost)
+				{
+					min_cost = et->m_Cost;
+					min_id = et->m_Target->m_Id;
+					min_idx = idx;
+					++idx;
+				}
+				else if(min_cost == et->m_Cost && min_id > et->m_Target->m_Id)
+				{
+					min_id = et->m_Target->m_Id;
+					min_idx = idx;
+					++idx;
+				}
+				else
+					++idx;
+			}
+
+			if(min_idx > 0)
+				std::swap((*it)->m_Edges[0], (*it)->m_Edges[min_idx]);
+		}
+	}
+
+	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().m_Target;
+
+				if( b->m_Valid)
+				{
+					NodePointerType best_b = b->m_Edges.front().m_Target;
+
+					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().m_Target;
+
+				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(),
+							EdgeEqualTo(target));
+	}
+
+	template<class TSegmenter>
+	void
+	GraphOperations<TSegmenter>::UpdateNeighbors(NodePointerType a, NodePointerType b)
+	{
+		unsigned int boundary;
+
+		/* Explore the neighbors of b */
+		for (EdgeIterator eit = b->m_Edges.begin();
+			 eit != b->m_Edges.end(); ++eit)
+		{
+			// Retrieve the edge targeting node b.
+			NodePointerType neigh_b = eit->m_Target;
+			EdgeIterator 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 = eit->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. */
+				EdgeIterator 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. */
+					EdgeType edge;
+					edge.m_Boundary = boundary;
+					edge.m_Cost = 0.0f;
+					edge.m_Target = a;
+					neigh_b->m_Edges.push_back(edge);
+
+					/* Add an edge from node a targeting node neigh_b. */
+					edge.m_Target = neigh_b;
+					a->m_Edges.push_back(edge);
+				}
+				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. */
+					EdgeIterator 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)
+	{
+		/* Step 1: update the contour */
+		a->m_Bbox = ContourOperations::MergeBoundingBoxes(a->m_Bbox, b->m_Bbox);
+		a->m_Contour = ContourOperations::MergeContours(a->m_Id, b->m_Id, a->m_Contour, b->m_Contour,
+														a->m_Bbox, width, height);
+
+		/* 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;
+	}
+
+	template<class TSegmenter>
+	void
+	GraphOperations<TSegmenter>::RemoveExpiredNodes(GraphType& graph)
+	{
+		NodeIterator nit = std::remove_if(graph.m_Nodes.begin(), graph.m_Nodes.end(), IsExpired());
+		graph.m_Nodes.erase(nit, graph.m_Nodes.end());
+	}
+
+	template<class TSegmenter>
+	bool
+	GraphOperations<TSegmenter>::PerfomOneIterationWithLMBF(GraphType& graph,
+															SegmenterType& seg,
+															const float threshold,
+															const unsigned int width,
+															const unsigned int height)
+	{
+		bool merged = false;
+
+		/* Update the costs of merging between adjacent nodes */
+		UpdateMergingCosts(graph, seg);
+
+		for(NodeIterator nit = graph.m_Nodes.begin();
+			nit != graph.m_Nodes.end(); ++nit)
+		{
+			
+			NodePointerType res_node = CheckLMBF(*nit, threshold);
+
+			if(res_node)
+				{
+					seg.UpdateSpecificAttributes(res_node, res_node->m_Edges.front().m_Target);
+					UpdateInternalAttributes(res_node, res_node->m_Edges.front().m_Target,
+											 width, height);
+					merged = true;
+				}
+		}
+
+		RemoveExpiredNodes(graph);
+
+		if(graph.m_Nodes.size() < 2)
+			return false;
+
+		return merged;
+	}
+
+	template<class TSegmenter>
+	bool
+	GraphOperations<TSegmenter>::PerfomOneIterationWithBF(GraphType& graph,
+														  SegmenterType& seg,
+														  const float threshold,
+														  const unsigned int width,
+														  const unsigned int height)
+	{
+		bool merged = false;
+
+		/* Update the costs of merging between adjacent nodes */
+		UpdateMergingCosts(graph, seg);
+
+		for(NodeIterator nit = graph.m_Nodes.begin();
+			nit != graph.m_Nodes.end(); ++nit)
+		{
+			NodePointerType res_node = CheckBF(*nit, threshold);
+
+			if(res_node)
+				{
+					seg.UpdateSpecificAttributes(res_node, res_node->m_Edges.front().m_Target);
+					UpdateInternalAttributes(res_node, res_node->m_Edges.front().m_Target,
+											 width, height);
+					merged = true;
+				}
+		}
+
+		RemoveExpiredNodes(graph);
+
+		if(graph.m_Nodes.size() < 2)
+			return false;
+
+		return merged;
+	}
+
+	template<class TSegmenter>
+	bool
+	GraphOperations<TSegmenter>::PerfomAllIterationsWithLMBFAndConstThreshold(GraphType& graph,
+																			  SegmenterType& seg,
+																			  const float threshold,
+																			  const unsigned int numberOfIterations,
+																			  const unsigned int width,
+																			  const unsigned int height)
+	{
+		bool merged = true;
+		unsigned int iterations = 0;
+
+		while(merged &&
+			  iterations < numberOfIterations &&
+			  graph.m_Nodes.size() > 1)
+		{
+			std::cout << "." << std::flush;
+			++iterations;
+
+			merged = PerfomOneIterationWithLMBF(graph, seg, threshold,
+												width, height);
+		}
+		std::cout << std::endl;
+		
+		if(graph.m_Nodes.size() < 2)
+			return false;
+
+		return merged;
+	}
+
+	template<class TSegmenter>
+	bool
+	GraphOperations<TSegmenter>::PerfomAllIterationsWithBFAndConstThreshold(GraphType& graph,
+																			SegmenterType& seg,
+																			const float threshold,
+																			const unsigned int numberOfIterations,
+																			const unsigned int width,
+																			const unsigned int height)
+	{
+		bool merged = true;
+		unsigned int iterations = 0;
+
+		while(merged &&
+			  iterations < numberOfIterations &&
+			  graph.m_Nodes.size() > 1)
+		{
+			++iterations;
+
+			merged = PerfomOneIterationWithBF(graph, seg, threshold,
+											  width, height);
+		}
+
+		if(graph.m_Nodes.size() < 2)
+			return false;
+
+		return merged;
+	}
+	
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::LargeScaleInitNodes(GraphType& graph,
+														  SegmenterType& seg,
+														  const unsigned int tileWidth,
+														  const unsigned int tileHeight,
+														  const unsigned int nbTilesX,
+														  const unsigned int nbTilesY,
+														  const unsigned int tileId,
+														  const unsigned int margin,
+														  const unsigned int imageWidth,
+														  const unsigned int imageHeight,
+														  CONNECTIVITY mask)
+	{
+		const unsigned int tileX = tileId % nbTilesX;
+		const unsigned int tileY = tileId / nbTilesX;
+
+		unsigned int xsize = tileWidth, ysize = tileHeight;
+		unsigned int xmargin = 0, ymargin = 0;
+		
+		if(tileX > 0)
+		{
+			xmargin = margin;
+			xsize += margin;
+		}
+		if(tileY > 0)
+		{
+			ymargin = margin;
+			ysize += margin;
+		}
+		
+		if(tileX < nbTilesX - 1)
+			xsize += margin;
+
+		if(tileY < nbTilesY - 1)
+			ysize += margin;
+		
+		const long unsigned int num_nodes = xsize * ysize;
+		graph.m_Nodes.reserve(num_nodes);
+
+		unsigned int xtile, ytile, ximg, yimg;
+
+		for(long unsigned int i = 0;
+			i < num_nodes;
+			++i)
+		{
+			NodePointerType n(new NodeType);
+
+			/* Compute the coordinates of its first pixel */
+			xtile = i % xsize;
+			ytile = i / xsize;
+			ximg = tileX * tileWidth + xtile - xmargin;
+			yimg = tileY * tileHeight + ytile - ymargin;
+			n->m_Id = yimg * imageWidth + ximg;
+			
+			n->m_Valid = true;
+			n->m_Expired = false;
+			n->m_Perimeter = 4;
+			n->m_Area = 1;
+			n->m_Bbox.m_UX = ximg;
+			n->m_Bbox.m_UY = yimg;
+			n->m_Bbox.m_W = 1;
+			n->m_Bbox.m_H = 1;
+			n->m_Contour.push_back(1);
+			n->m_Contour.push_back(2);
+			n->m_Contour.push_back(3);
+			n->m_Contour.push_back(0);
+			graph.m_Nodes.push_back(n);
+		}
+
+		if(mask == FOUR)
+		{
+			for(long unsigned int i = 0;
+				i < num_nodes;
+				++i)
+			{
+				long int neighborhood[4];
+				FOURNeighborhood(neighborhood, i, xsize, ysize);
+				for(short j = 0; j < 4; ++j)
+				{
+					if(neighborhood[j] > -1)
+					{
+						EdgeType edge;
+						edge.m_Boundary = 1;
+						edge.m_Cost = 0.0f;
+						edge.m_Target = graph.m_Nodes[neighborhood[j]];
+						graph.m_Nodes[i]->m_Edges.push_back(edge);
+					}
+				}
+			}
+		}
+		else
+		{
+			for(long unsigned int i = 0;
+				i < num_nodes;
+				++i)
+			{
+				long int neighborhood[8];
+				EIGHTNeighborhood(neighborhood, i, xsize, ysize);
+				for(short j = 0; j < 8; ++j)
+				{
+					if(neighborhood[j] > -1)
+					{
+						EdgeType edge;
+						if(j % 2 > 0)
+							edge.m_Boundary = 0;
+						else
+							edge.m_Boundary = 1;
+						edge.m_Cost = 0.0f;
+						edge.m_Target = graph.m_Nodes[neighborhood[j]];
+						graph.m_Nodes[i]->m_Edges.push_back(edge);
+					}
+				}
+			}
+		}
+		seg.InitFromImage();
+	}
+
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::RemoveEdgesToUnstableRegion(NodePointerType r)
+	{
+		for (EdgeIterator et = r->m_Edges.begin();
+			 et != r->m_Edges.end(); ++et)
+		{
+			EdgeIterator toR = FindEdge(et->m_Target, r);
+			et->m_Target->m_Edges.erase(toR);
+		}
+	}
+
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::RemoveUnstableRegions(GraphType& graph,
+															unsigned int lowerRow,
+															unsigned int upperRow,
+															unsigned int lowerCol,
+															unsigned int upperCol,
+															unsigned int imageWidth)
+	{
+		for(NodeIterator nit = graph.m_Nodes.begin();
+			nit != graph.m_Nodes.end(); ++nit)
+		{
+			if((*nit)->m_Bbox.m_UX >= lowerCol &&
+			   (*nit)->m_Bbox.m_UY >= lowerRow &&
+			   ((*nit)->m_Bbox.m_W + (*nit)->m_Bbox.m_UX - 1) <= upperCol &&
+			   ((*nit)->m_Bbox.m_H + (*nit)->m_Bbox.m_UY - 1) <= upperRow)
+				continue;
+			else if((*nit)->m_Bbox.m_UX > upperCol ||
+					(*nit)->m_Bbox.m_UY > upperRow ||
+					((*nit)->m_Bbox.m_W + (*nit)->m_Bbox.m_UX - 1) < lowerCol ||
+					((*nit)->m_Bbox.m_H + (*nit)->m_Bbox.m_UY - 1) < lowerRow)
+			{
+				RemoveEdgesToUnstableRegion(*nit);
+				(*nit)->m_Expired = true;
+			}
+			else
+			{
+				bool unstable = true;
+				unsigned int pix_col, pix_row;
+				ContourOperations::PixelList pixels = ContourOperations::GeneratePixels((*nit)->m_Id,
+																						(*nit)->m_Contour,
+																						(*nit)->m_Bbox,
+																						imageWidth);
+				for (ContourOperations::PixelConstIterator pit = pixels.begin();
+					 pit != pixels.end(); ++pit)
+				{
+					pix_col = (*pit) % imageWidth;
+					pix_row = (*pit) / imageWidth;
+
+					if(pix_col >= lowerCol && pix_col <= upperCol &&
+					   pix_row >= lowerRow && pix_row <= upperRow)
+					{
+						unstable = false;
+						break;
+					}
+				}
+
+				if(unstable)
+				{
+					RemoveEdgesToUnstableRegion(*nit);
+					(*nit)->m_Expired = true;
+				}
+			}
+		}
+		// Remove unstable regions from the graph
+		RemoveExpiredNodes(graph);
+	}
+
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::AddNeighboringLayer(NodePointerType n,
+														  GraphType& g,
+														  std::set<long unsigned int>& bset)
+	{
+		typedef std::pair<long unsigned int, NodePointerType> KeyValue;
+		
+		// First we add its neighboring nodes
+		for (EdgeIterator eit = n->m_Edges.begin();
+			 eit != n->m_Edges.end(); ++eit)
+		{
+			if(bset.find(eit->m_Target->m_Id) == bset.end())
+			{
+				g.m_Nodes.push_back(eit->m_Target);
+				bset.insert(eit->m_Target->m_Id);
+			}
+
+			// We explore the neighboring nodes of this neighboring nodes
+			for (EdgeIterator neit = eit->m_Target->m_Edges.begin();
+				 neit != eit->m_Target->m_Edges.end(); ++neit)
+			{
+				if(bset.find(neit->m_Target->m_Id) == bset.end())
+				{
+					g.m_Nodes.push_back(neit->m_Target);
+					bset.insert(neit->m_Target->m_Id);
+				}	
+			}
+		}
+	}
+
+	template<class TSegmenter>
+	bool GraphOperations<TSegmenter>::IsOnTileBorder(const unsigned int pixelRow,
+													 const unsigned int pixelCol,
+													 const unsigned int lowerRow,
+													 const unsigned int upperRow,
+													 const unsigned int lowerCol,
+													 const unsigned int upperCol,
+													 const unsigned int imageWidth,
+													 const unsigned int imageHeight)
+	{
+		long int rowBounds[2];
+		long int colBounds[2];
+
+		rowBounds[0] = ((lowerRow > 0 ) ? lowerRow : -1);
+		rowBounds[1] = ((upperRow < imageHeight - 1) ? upperRow : -1);
+		colBounds[0] = ((lowerCol > 0) ? lowerCol : -1);
+		colBounds[1] = ((upperCol < imageWidth - 1) ? upperCol : -1);
+
+		for(int i = 0; i < 2; i++)
+		{
+			if(rowBounds[i] > -1 && pixelRow == rowBounds[i] && pixelCol >= lowerCol && pixelCol <= upperCol)
+				return true;
+		}
+
+		for(int i = 0; i < 2; i++)
+		{
+			if(colBounds[i] > -1 && pixelCol == colBounds[i] && pixelRow >= lowerRow && pixelRow <= upperRow)
+				return true;
+		}
+
+		return false;
+	}
+
+	template<class TSegmenter>
+	bool GraphOperations<TSegmenter>::IsOnTileBorderAndAdjacentTileBorder(const unsigned int pixelRow,
+																		  const unsigned int pixelCol,
+																		  const unsigned int lowerRow,
+																		  const unsigned int upperRow,
+																		  const unsigned int lowerCol,
+																		  const unsigned int upperCol,
+																		  const unsigned int imageWidth,
+																		  const unsigned int imageHeight)
+	{
+		long int rowBounds[4];
+		long int colBounds[4];
+
+		if(lowerRow > 0)
+		{
+			rowBounds[0] = lowerRow - 1;
+			rowBounds[1] = lowerRow;
+		}
+		else
+		{
+			rowBounds[0] = -1;
+			rowBounds[1] = -1;
+		}
+
+		if(upperRow < imageHeight - 1)
+		{
+			rowBounds[2] = upperRow;
+			rowBounds[3] = upperRow + 1;
+		}
+		else
+		{
+			rowBounds[2] = -1;
+			rowBounds[3] = -1;
+		}
+
+		if(lowerCol > 0)
+		{
+			colBounds[0] = lowerCol - 1;
+			colBounds[1] = lowerCol;
+		}
+		else
+		{
+			colBounds[0] = -1;
+			colBounds[1] = -1;
+		}
+
+		if(upperCol < imageWidth - 1)
+		{
+			colBounds[2] = upperCol;
+			colBounds[3] = upperCol + 1;
+		}
+		else
+		{
+			colBounds[2] = -1;
+			colBounds[3] = -1;
+		}
+
+		for(int i = 0; i < 4; i++)
+		{
+			if(rowBounds[i] > -1 && pixelRow == rowBounds[i])
+				return true;
+		}
+
+		for(int i = 0; i < 4; i++)
+		{
+			if(colBounds[i] > -1 && pixelCol == colBounds[i])
+				return true;
+		}
+
+		return false;
+	}
+
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::BuildStabilityMargins(const GraphType& graph,
+															GraphType& subgraph,
+															const unsigned int lowerRow,
+															const unsigned int upperRow,
+															const unsigned int lowerCol,
+															const unsigned int upperCol,
+															const unsigned int imageWidth,
+															const unsigned int imageHeight)
+	{	
+		std::set<long unsigned int> used_regions;
+		typedef typename std::set<long unsigned int>::iterator NIterator;
+
+		
+		for(NodeConstIterator nit = graph.m_Nodes.begin();
+			nit != graph.m_Nodes.end(); ++nit)
+		{
+
+			if(IsBboxStrictlyInsideTile(*nit, lowerRow, upperRow, lowerCol, upperCol))
+				continue;
+			else
+			{
+				unsigned int pix_col, pix_row;
+				ContourOperations::PixelList pixels = ContourOperations::GeneratePixels((*nit)->m_Id,
+																						(*nit)->m_Contour,
+																						(*nit)->m_Bbox,
+																						imageWidth);
+
+				bool is_border = false;
+				for (ContourOperations::PixelConstIterator pit = pixels.begin();
+					 pit != pixels.end(); ++pit)
+				{
+					pix_col = (*pit) % imageWidth;
+					pix_row = (*pit) / imageWidth;
+
+					if(IsOnTileBorder(pix_row, pix_col, lowerRow, upperRow,
+									  lowerCol, upperCol, imageWidth, imageHeight))
+					{
+						is_border = true;
+						if(used_regions.find((*nit)->m_Id) == used_regions.end())
+						{
+							subgraph.m_Nodes.push_back(*nit);
+							used_regions.insert((*nit)->m_Id);
+						}
+						
+						AddNeighboringLayer(*nit, subgraph, used_regions);
+					}
+
+					if(is_border)
+						break;
+				}
+			}
+		}
+
+		// We remove the useless targeted regions
+		for(NodeIterator nit = subgraph.m_Nodes.begin(); nit != subgraph.m_Nodes.end(); ++nit)
+		{
+			EdgeIterator ret = std::remove_if((*nit)->m_Edges.begin(), (*nit)->m_Edges.end(), IsInSet(used_regions));
+			(*nit)->m_Edges.erase(ret, (*nit)->m_Edges.end());
+		}
+	}
+
+	template<class TSegmenter>
+	bool GraphOperations<TSegmenter>::IsBboxStrictlyInsideTile(NodePointerType n,
+															   const unsigned int lowerRow,
+															   const unsigned int upperRow,
+															   const unsigned int lowerCol,
+															   const unsigned int upperCol)
+	{
+		return (n->m_Bbox.m_UX > lowerCol &&
+				n->m_Bbox.m_UY > lowerRow &&
+				(n->m_Bbox.m_W + n->m_Bbox.m_UX - 1) < upperCol &&
+				(n->m_Bbox.m_H + n->m_Bbox.m_UY - 1) < upperRow);
+	}
+
+	template<class TSegmenter>
+	bool GraphOperations<TSegmenter>::IsBboxInsideTile(NodePointerType n,
+													   const unsigned int lowerRow,
+													   const unsigned int upperRow,
+													   const unsigned int lowerCol,
+													   const unsigned int upperCol)
+	{
+		return (n->m_Bbox.m_UX >= lowerCol &&
+				n->m_Bbox.m_UY >= lowerRow &&
+				(n->m_Bbox.m_W + n->m_Bbox.m_UX - 1) <= upperCol &&
+				(n->m_Bbox.m_H + n->m_Bbox.m_UY - 1) <= upperRow);
+	}
+	
+
+	template<class TSegmenter>
+	void GraphOperations<TSegmenter>::MergeGraphs(GraphType& graph,
+												  const unsigned int lowerRow,
+												  const unsigned int upperRow,
+												  const unsigned int lowerCol,
+												  const unsigned int upperCol,
+												  const unsigned int imageWidth,
+												  const unsigned int imageHeight,
+												  bool useless)
+	{	
+		{
+			std::map<long unsigned int, std::vector<NodePointerType> > border_map;
+			typedef typename std::map<long unsigned int, std::vector<NodePointerType> >::iterator BorderMapIterator;
+		
+			// First step : detect the duplicated nodes
+		
+			for(NodeConstIterator nit = graph.m_Nodes.begin();
+				nit != graph.m_Nodes.end(); ++nit)
+			{
+				if(!IsBboxStrictlyInsideTile(*nit, lowerRow, upperRow, lowerCol, upperCol))
+				{	
+					unsigned int pix_col, pix_row;
+					ContourOperations::PixelList pixels = ContourOperations::GenerateBorderPixels((*nit)->m_Id,
+																								  (*nit)->m_Contour,
+																								  imageWidth);
+					
+					for (ContourOperations::PixelConstIterator pit = pixels.begin();
+						 pit != pixels.end(); ++pit)
+					{
+						pix_col = (*pit) % imageWidth;
+						pix_row = (*pit) / imageWidth;
+
+						if(IsOnTileBorderAndAdjacentTileBorder(pix_row, pix_col, lowerRow, upperRow,
+															   lowerCol, upperCol, imageWidth, imageHeight))
+						{	
+							if(border_map.find(*pit) != border_map.end())
+							{
+								NodeIterator is_here =
+									std::find_if(border_map[*pit].begin(), border_map[*pit].end(), NodePtrComparator(*nit));
+								if(is_here == border_map[*pit].end())
+									border_map[*pit].push_back(*nit);
+							}
+							else
+								border_map[*pit].push_back(*nit);
+						}
+					}
+				}
+			}
+
+			// Second step : remove the duplicated regions
+			for(BorderMapIterator bmit = border_map.begin(); bmit != border_map.end(); ++bmit)
+			{
+				if(bmit->second.size() > 1) // More than one node ?
+				{
+					NodePointerType ref = bmit->second.front();
+
+					// Explore the other duplicated nodes
+					for(NodeIterator nit = bmit->second.begin() + 1; nit != bmit->second.end(); ++nit)
+					{
+						// Explore their edges
+						for(EdgeIterator eit = (*nit)->m_Edges.begin(); eit != (*nit)->m_Edges.end(); ++eit)
+						{
+							// We find the edge targeting the duplicated node and we remove it
+							EdgeIterator toDoublon = FindEdge(eit->m_Target, *nit);
+							assert(toDoublon != eit->m_Target->m_Edges.end());
+							unsigned int boundary = toDoublon->m_Boundary;
+							eit->m_Target->m_Edges.erase(toDoublon);
+
+							EdgeIterator RefToTarget = FindEdge(ref, eit->m_Target);
+							if(RefToTarget == ref->m_Edges.end())
+							{	
+								// We create a new edge between the reference and the target
+								EdgeType edge;
+								edge.m_Boundary = boundary;
+								edge.m_Cost = 0.0f;
+								edge.m_Target = ref;
+								eit->m_Target->m_Edges.push_back(edge);
+								edge.m_Target = eit->m_Target;
+								ref->m_Edges.push_back(edge);
+							}
+						}
+						(*nit)->m_Expired = true;
+					}
+				
+					ContourOperations::PixelList pixels = ContourOperations::GenerateBorderPixels(ref->m_Id,
+																								  ref->m_Contour, 
+																								  imageWidth);
+					for (ContourOperations::PixelConstIterator pit = pixels.begin();
+						 pit != pixels.end(); ++pit)
+					{
+						if(border_map.find(*pit) != border_map.end())
+						{
+							border_map[*pit].clear();
+							border_map[*pit].push_back(ref);
+						}
+					}
+				}
+			}
+			RemoveExpiredNodes(graph);
+
+			// Update the border regions
+			{
+				unsigned int boundary;
+				for(BorderMapIterator bmit = border_map.begin(); bmit != border_map.end(); ++bmit)
+				{
+					long int neighborhood[4];
+					FOURNeighborhood(neighborhood, bmit->first, imageWidth, imageHeight);
+					for(short j = 0; j < 4; ++j)
+					{
+						if(neighborhood[j] > -1)
+						{
+							BorderMapIterator isInMap = border_map.find(neighborhood[j]);
+							if(isInMap != border_map.end())
+							{
+								NodePointerType curr_region = bmit->second.front();
+								NodePointerType target_region = isInMap->second.front();
+								if(curr_region != target_region)
+								{
+									EdgeIterator ToTarget = FindEdge(curr_region, target_region);
+									if(ToTarget == curr_region->m_Edges.end())
+									{
+										boundary = 0;
+										ContourOperations::PixelList pixels =
+											ContourOperations::GenerateBorderPixels(curr_region->m_Id,
+																					curr_region->m_Contour, 
+																					imageWidth);
+									
+										for (ContourOperations::PixelConstIterator pit = pixels.begin();
+											 pit != pixels.end(); ++pit)
+										{
+											if(border_map.find(*pit) != border_map.end())
+											{
+												long int curr_neighborhood[4];
+												FOURNeighborhood(curr_neighborhood, *pit, imageWidth, imageHeight);
+												for(short k = 0; k<4; ++k)
+												{
+													if(curr_neighborhood[k] > -1)
+													{
+														BorderMapIterator Detect = border_map.find(curr_neighborhood[k]);
+														if(Detect != border_map.end() &&
+														   Detect->second.front() == target_region)
+															boundary++;
+													}
+												}
+											}
+										}
+
+										EdgeType edge;
+										edge.m_Boundary = boundary;
+										edge.m_Cost = 0.0f;
+										edge.m_Target = target_region;
+										curr_region->m_Edges.push_back(edge);
+										edge.m_Target = curr_region;
+										target_region->m_Edges.push_back(edge);
+									}
+								}
+							}
+						}
+					}
+				}
+			}
+		}
+
+		/*for(NodeConstIterator nit = graph.m_Nodes.begin();
+			nit != graph.m_Nodes.end(); ++nit)
+		{
+			if((*nit)->m_Id == 577)
+			{
+				std::cout << (*nit)->m_Id << ": ";
+					for(EdgeIterator eit = (*nit)->m_Edges.begin(); eit != (*nit)->m_Edges.end(); ++eit)
+					{
+						std::cout << eit->m_Target->m_Id << " " << eit->m_Boundary << "\t";
+					}
+					std::cout << std::endl;
+			}
+			}*/
+			
+		// Remove useless nodes
+		if(useless)
+		{
+			std::set<NodePointerType> border_nodes;
+			std::set<NodePointerType> margin_nodes;
+			typedef typename std::set<NodePointerType>::iterator SetIterator;
+
+			// Retrieve the border nodes
+			for(NodeIterator nit = graph.m_Nodes.begin(); nit != graph.m_Nodes.end(); ++nit)
+			{	
+				if(!IsBboxStrictlyInsideTile(*nit, lowerRow, upperRow, lowerCol, upperCol))
+				{
+					unsigned int pix_col, pix_row;
+					ContourOperations::PixelList pixels = ContourOperations::GenerateBorderPixels((*nit)->m_Id,
+																								  (*nit)->m_Contour,
+																								  imageWidth);
+					
+					
+					for (ContourOperations::PixelConstIterator pit = pixels.begin();
+						 pit != pixels.end(); ++pit)
+					{
+						pix_col = (*pit) % imageWidth;
+						pix_row = (*pit) / imageWidth;
+
+						if(IsOnTileBorder(pix_row, pix_col, lowerRow, upperRow,
+										  lowerCol, upperCol, imageWidth, imageHeight))
+						{
+							border_nodes.insert(*nit);
+							break;
+						}
+						else
+							continue;
+					}
+				}
+			}
+
+			// Add the margin nodes of the border nodes
+			for(SetIterator sit = border_nodes.begin(); sit != border_nodes.end(); ++sit)
+			{
+				for(EdgeIterator eit = (*sit)->m_Edges.begin(); eit != (*sit)->m_Edges.end(); ++eit)
+				{
+					if(margin_nodes.find(eit->m_Target) == margin_nodes.end())
+					{
+						margin_nodes.insert(eit->m_Target);
+					}
+
+					for(EdgeIterator eeit = eit->m_Target->m_Edges.begin(); eeit != eit->m_Target->m_Edges.end(); ++eeit)
+					{
+						if(margin_nodes.find(eeit->m_Target) == margin_nodes.end())
+							margin_nodes.insert(eeit->m_Target);
+					}
+				}
+			}
+
+			for(NodeIterator nit = graph.m_Nodes.begin(); nit != graph.m_Nodes.end(); ++nit)
+			{	
+				if(!IsBboxInsideTile(*nit, lowerRow, upperRow, lowerCol, upperCol) &&
+				   margin_nodes.find(*nit) == margin_nodes.end())
+				{
+					// We remove the edges of the adjacent nodes tageting this node
+					for(EdgeIterator eit = (*nit)->m_Edges.begin(); eit != (*nit)->m_Edges.end(); ++eit)
+					{
+						EdgeIterator ret = FindEdge(eit->m_Target, *nit);
+						eit->m_Target->m_Edges.erase(ret);
+					}
+
+					// This node has to be removed
+					(*nit)->m_Expired = true;
+				}
+			}
+
+			RemoveExpiredNodes(graph);
+		}
+
+		/*for(NodeIterator nit = graph.m_Nodes.begin(); nit != graph.m_Nodes.end(); ++nit)
+		{
+			std::cout << (*nit)->m_Id << " : ";
+			for(EdgeIterator eit = (*nit)->m_Edges.begin(); eit != (*nit)->m_Edges.end(); ++eit)
+			{
+				std::cout << eit->m_Target->m_Id << " " << eit->m_Boundary << "\t";
+			}
+			std::cout <<  std::endl;
+			}*/
+
+		// Remove useless regions
+		
+	}
+
+	template<class TSegmenter>
+	long long unsigned int GraphOperations<TSegmenter>::GetMemorySpace(const GraphType& graph)
+	{
+		long long unsigned int res = graph.m_Nodes.size() * NodeSize;
+
+		for(NodeConstIterator nit = graph.m_Nodes.begin(); nit != graph.m_Nodes.end(); ++nit)
+		{
+			res += ((*nit)->m_Edges.size() * EdgeSize);
+		}
+
+		return res;
+	}
+
+} // end of namespace lsrm
+
+#endif
+
+
+
+
+
diff --git a/Code/lsrmGraphToOtbImage.h b/Code/lsrmGraphToOtbImage.h
new file mode 100644
index 0000000000000000000000000000000000000000..b6ac7bc82981ad00b994d10c0c4ec8b296987b0f
--- /dev/null
+++ b/Code/lsrmGraphToOtbImage.h
@@ -0,0 +1,68 @@
+#ifndef __LSRM_GRAPH_TO_OTBIMAGE_H
+#define __LSRM_GRAPH_TO_OTBIMAGE_H
+#include <itkRGBPixel.h>
+#include <otbImage.h>
+#include <otbVectorImage.h>
+#include <otbImageFileReader.h>
+#include <otbImageFileWriter.h>
+#include "lsrmGraph.h"
+#include "lsrmContourOperations.h"
+#include <string>
+#include <stdlib.h>
+#include <time.h>
+
+namespace lsrm
+{
+	template<class TGraph>
+	class GraphToOtbImage
+	{
+	public:
+		
+		/* Some convenient typedefs */
+		typedef TGraph GraphType;
+		typedef typename GraphType::NodeType NodeType;
+		typedef std::vector< boost::shared_ptr<NodeType> > NodeList;
+		typedef typename NodeList::const_iterator NodeConstIterator;
+
+		/*
+		 * Given a graph of nodes and the size of the image, it returns
+		 * a raster image where to every pixel of a region is assigned
+		 * the same unique label value.
+		 *
+		 * @params
+		 * const GraphType& graph : graph of nodes
+		 * const unsigned int width : width of the output labeled image.
+		 * const unsigned int height : height of the output labeled image.
+		 * const std::string& outputFileName : filename of the output image.
+		 */
+		static void WriteLabelImage(const GraphType& graph,
+									const unsigned int width,
+									const unsigned int height,
+									const std::string& outputFileName);
+
+		/*
+		 * Given a graph of nodes and the size of the image, it returns
+		 * a raster image where to every pixel of a region is assigned
+		 * the same rgb color.
+		 *
+		 * @params
+		 * const GraphType& graph : graph of nodes
+		 * const unsigned int width : width of the output rgb image.
+		 * const unsigned int height : height of the output rgb image.
+		 * const std::string& outputFileName : filename of the output image.
+		 */
+		static void WriteOutputRGBImage(const GraphType& graph,
+										const unsigned int width,
+										const unsigned int height,
+										const std::string& outputFileName);
+
+		static void WriteContourImage(const GraphType& graph,
+									  const std::string& inputFileName,
+									  const std::string& outputFileName);
+
+		
+	};
+	
+} // end of namespace lsrm
+#include "lsrmGraphToOtbImage.txx"
+#endif
diff --git a/Code/lsrmGraphToOtbImage.txx b/Code/lsrmGraphToOtbImage.txx
new file mode 100644
index 0000000000000000000000000000000000000000..1005ddeafa6811a910d5d4833bd4bfc144bb2047
--- /dev/null
+++ b/Code/lsrmGraphToOtbImage.txx
@@ -0,0 +1,172 @@
+#ifndef __LSRM_GRAPH_TO_OTBIMAGE_TXX
+#define __LSRM_GRAPH_TO_OTBIMAGE_TXX
+#include "lsrmGraphToOtbImage.h"
+
+namespace lsrm
+{
+	template<class TGraph>
+	void GraphToOtbImage<TGraph>::WriteLabelImage(const GraphType& graph,
+												  const unsigned int width,
+												  const unsigned int height,
+												  const std::string& outputFileName)
+	{
+		typedef long unsigned int LabelPixelType;
+
+		typedef otb::Image<LabelPixelType, 2> LabelImageType;
+
+		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();
+		
+		// Start at 1 (value 0 can be used for invalid pixels)
+		long unsigned int label = 1;
+		for(NodeConstIterator nit = graph.m_Nodes.begin();
+			nit != graph.m_Nodes.end(); ++nit)
+		{
+			ContourOperations::PixelList pixels = ContourOperations::GeneratePixels((*nit)->m_Id,
+																					(*nit)->m_Contour,
+																					(*nit)->m_Bbox,
+																					width);
+
+			for (ContourOperations::PixelConstIterator pit = pixels.begin();
+				 pit != pixels.end(); ++pit)
+			{
+				index[0] = (*pit) % width;
+				index[1] = (*pit) / width;
+				label_img->SetPixel(index, label);
+			}
+			++label;
+		}
+
+		typedef otb::ImageFileWriter<LabelImageType> LabelWriterType;
+		LabelWriterType::Pointer label_writer = LabelWriterType::New();
+		label_writer->SetInput(label_img);
+		label_writer->SetFileName(outputFileName);
+		label_writer->Update();
+	}
+
+	template<class TGraph>
+	void GraphToOtbImage<TGraph>::WriteOutputRGBImage(const GraphType& graph,
+													  const unsigned int width,
+													  const unsigned int height,
+													  const std::string& outputFileName)
+	{
+		typedef unsigned char RGBPixelType;
+
+		typedef otb::VectorImage<RGBPixelType, 2> RGBImageType;
+		typename RGBImageType::PixelType pixelValue;
+
+		RGBImageType::IndexType index;
+		RGBImageType::SizeType size;
+		RGBImageType::RegionType region;
+
+		index[0] = 0; index[1] = 0;
+		size[0] = width; size[1] = height;
+		region.SetIndex(index);
+		region.SetSize(size);
+
+		RGBImageType::Pointer rgb_img = RGBImageType::New();
+		rgb_img->SetRegions(region);
+		rgb_img->SetNumberOfComponentsPerPixel(3);
+		rgb_img->Allocate();
+		pixelValue.Reserve(3);
+		pixelValue[0] = 0;
+		pixelValue[1] = 0;
+		pixelValue[2] = 0;
+
+		for(unsigned int y = 0; y < height; ++y)
+		{
+			for(unsigned int x = 0; x < width; ++x)
+			{
+				index[0] = x;
+				index[1] = y;
+				rgb_img->SetPixel(index, pixelValue);
+			}
+		}
+
+		srand(time(NULL));
+		unsigned char c1, c2, c3;
+		for(NodeConstIterator nit = graph.m_Nodes.begin();
+			nit != graph.m_Nodes.end(); ++nit)
+		{
+			ContourOperations::PixelList pixels = ContourOperations::GeneratePixels((*nit)->m_Id,
+																					(*nit)->m_Contour,
+																					(*nit)->m_Bbox,
+																					width);
+			c1 = rand() % 256;
+			c2 = rand() % 256;
+			c3 = rand() % 256;
+			for (ContourOperations::PixelConstIterator pit = pixels.begin();
+				 pit != pixels.end(); ++pit)
+			{
+				index[0] = (*pit) % width;
+				index[1] = (*pit) / width;
+				pixelValue[0] = c1;
+				pixelValue[1] = c2;
+				pixelValue[2] = c3;
+				rgb_img->SetPixel(index, pixelValue);
+			}
+		}
+
+		typedef otb::ImageFileWriter<RGBImageType> RGBWriterType;
+		RGBWriterType::Pointer rgb_writer = RGBWriterType::New();
+		rgb_writer->SetInput(rgb_img);
+		rgb_writer->SetFileName(outputFileName);
+		rgb_writer->Update();
+	}
+
+	template<class TGraph>
+	void GraphToOtbImage<TGraph>::WriteContourImage(const GraphType& graph,
+													const std::string& inputFileName,
+													const std::string& outputFileName)
+	{
+		typedef unsigned char RGBPixelType;
+		typedef otb::VectorImage<RGBPixelType, 2> RGBImageType;
+		typedef typename RGBImageType::IndexType RGBIndexType;
+		
+		typedef otb::ImageFileReader<RGBImageType> RGBImageReaderType;
+		typedef otb::ImageFileWriter<RGBImageType> RGBImageWriterType;
+
+		RGBImageReaderType::Pointer rgb_reader = RGBImageReaderType::New();
+		rgb_reader->SetFileName(inputFileName);
+		rgb_reader->Update();
+		RGBImageType::Pointer rgb_img = rgb_reader->GetOutput();
+		RGBIndexType index;
+
+		const unsigned int width = rgb_img->GetLargestPossibleRegion().GetSize()[0];
+
+		for(NodeConstIterator nit = graph.m_Nodes.begin();
+			nit != graph.m_Nodes.end(); ++nit)
+		{
+			ContourOperations::PixelList pixels = ContourOperations::GenerateBorderPixels((*nit)->m_Id,
+																						  (*nit)->m_Contour,
+																						  width);
+			for (ContourOperations::PixelConstIterator pit = pixels.begin();
+				 pit != pixels.end(); ++pit)
+			{
+				index[0] = (*pit) % width;
+				index[1] = (*pit) / width;
+				typename RGBImageType::PixelType rgb_pix = rgb_img->GetPixel(index);
+				for(unsigned int b=0; b<rgb_img->GetNumberOfComponentsPerPixel(); b++)
+					rgb_pix[b] = 0;
+				rgb_img->SetPixel(index, rgb_pix);
+			}
+		}
+
+		RGBImageWriterType::Pointer rgb_writer = RGBImageWriterType::New();
+		rgb_writer->SetFileName(outputFileName);
+		rgb_writer->SetInput(rgb_img);
+		rgb_writer->Update();
+	}
+	
+} // 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..d10e4d679e9f5854e3a6b59ca280d9635e6338d7
--- /dev/null
+++ b/Code/lsrmSegmenter.h
@@ -0,0 +1,112 @@
+#ifndef __LSRM_SEGMENTER_H
+#define __LSRM_SEGMENTER_H
+#include "macro-generator.h"
+#include "lsrmGraph.h"
+#include "lsrmGraphOperations.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;
+
+		/* Default constructor and destructor */
+		
+		Segmenter(){};
+		~Segmenter(){};
+
+		/*
+		 * This method performs the region merging segmentation.
+		 */
+		virtual void RunSegmentation() = 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;
+		
+		/* Set methods */
+		SetMacro(bool, DoBFSegmentation);
+		SetMacro(unsigned int, NumberOfIterations);
+		SetMacro(float, Threshold);
+		SetMacro(ParamType, Param);
+		SetMacro(std::string, InputFileName);
+		SetMacro(std::string, ClusteredImageFileName);
+		SetMacro(std::string, LabelImageFileName);
+		SetMacro(std::string, ContourImageFileName);
+
+	protected:
+
+		/* 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;
+
+		/* Graph */
+		GraphType m_Graph;
+
+		/* 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
+		std::string m_InputFileName;
+		std::string m_ClusteredImageFileName;
+		std::string m_LabelImageFileName;
+		std::string m_ContourImageFileName;
+	};
+} // end of namespace lsrm
+
+#endif
+
+
+
+
diff --git a/src/Library/macro-generator.h b/Code/macro-generator.h
similarity index 100%
rename from src/Library/macro-generator.h
rename to Code/macro-generator.h
diff --git a/src/Applications/BaatzSegmentation.cxx b/src/Applications/BaatzSegmentation.cxx
deleted file mode 100644
index aa3f16d78b894f137bdae83b1dc7ff734191c4f9..0000000000000000000000000000000000000000
--- a/src/Applications/BaatzSegmentation.cxx
+++ /dev/null
@@ -1,65 +0,0 @@
-/*=========================================================================
-
-  Program: Example of application of the Generic Region Merging Library
-  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.
-
-=========================================================================*/
-#include <iostream>
-
-// Your algorithm header file
-#include "baatz-algorithm.h"
-
-int main(int argc, char **argv)
-{
-	if(argc < 7 || argc > 9)
-	{
-		std::cout << "###############      Manual     #############################\n\n\n"
-				<< "Usage " << argv[0] << "\n\n\n"
-				<< " -- [input image (tif, jpg, png) <mandatory>]\n"
-				<< " -- [output rgb image (tif, jpg, png) <mandatory>]\n"
-				<< " -- [output label image (tif) <mandatory>]\n"
-				<< " -- [spectral weight (value >= 0 and <= 1) <mandatory>]\n"
-				<< " -- [shape weight (value >= 0 and <= 1) <mandatory>]\n"
-				<< " -- [scale parameter (positive value) <mandatory>]\n"
-				<< " -- [number of iterations using the Local Mutual Best Fitting heuristic (70 by default) <optional>]\n"
-				<< " -- [activation of the Best Fitting heuristic (value = 0 (desactivated) or 1 (activated) \
-					- activated by default) <optional>]\n\n\n"
-				<< "################################################################\n";
-		return 1;
-	}
-
-	const std::string inputFileName = argv[1];
-	const std::string outputFileName = argv[2];
-	const std::string outputLabelFileName = argv[3];
-	const float cw = atof(argv[4]);
-	const float sw = atof(argv[5]);
-	const float sp = atof(argv[6]);
-
-	typedef otb::VectorImage<float, 2> ImageType;
-	typedef BaatzAlgorithmRM<ImageType> SegmenterType;
-
-	BaatzParams params = {cw, sw, sp};
-	SegmenterType seg_;
-	seg_.SetInput(inputFileName);
-	seg_.SetOutputRGB(outputFileName);
-	seg_.SetOutputLabel(outputLabelFileName);
-	seg_.SetParameters(params);
-	if(argc >= 8)
-		seg_.SetNumberOfIterations(static_cast<unsigned int>(atoi(argv[7])));
-	if(argc == 9)
-		seg_.SetBestFitting(atoi(argv[8]));
-
-	seg_.Run();
-
-	return 0;
-}
diff --git a/src/Applications/EuclideanDistanceSegmentation.cxx b/src/Applications/EuclideanDistanceSegmentation.cxx
deleted file mode 100644
index ff6380febfae12278fd880829b8237345119778c..0000000000000000000000000000000000000000
--- a/src/Applications/EuclideanDistanceSegmentation.cxx
+++ /dev/null
@@ -1,63 +0,0 @@
-/*=========================================================================
-
-  Program: Example of application of the Generic Region Merging Library
-  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.
-
-=========================================================================*/
-#include <iostream>
-
-// Your algorithm header file
-#include "euc-dist-algorithm.h"
-
-int main(int argc, char **argv)
-{
-	if(argc < 6 || argc > 8)
-	{
-		std::cout << "###############      Manual     #############################\n\n\n"
-				<< "Usage " << argv[0] << "\n\n\n"
-				<< " -- [input image (tif, jpg, png) <mandatory>]\n"
-				<< " -- [output rgb image (tif, jpg, png) <mandatory>]\n"
-				<< " -- [output label image (tif) <mandatory>]\n"
-				<< " -- [maximum distance <mandatory>]\n"
-				<< " -- [minimum distance <mandatory>]\n"
-				<< " -- [number of iterations using the Local Mutual Best Fitting heuristic (70 by default) <optional>]\n"
-				<< " -- [activation of the Best Fitting heuristic (value = 0 (desactivated) or 1 (activated) \
-					- activated by default) <optional>]\n\n\n"
-				<< "################################################################\n";
-		return 1;
-	}
-
-	const std::string inputFileName = argv[1];
-	const std::string outputFileName = argv[2];
-	const std::string outputLabelFileName = argv[3];
-	const float max_dist = atof(argv[4]);
-	float min_dist = atof(argv[5]);
-
-	typedef otb::VectorImage<float, 2> ImageType;
-	typedef EuclideanDistanceRM<ImageType> SegmenterType;
-
-	EucDistParams params = {max_dist, min_dist};
-	SegmenterType seg_;
-	seg_.SetInput(inputFileName);
-	seg_.SetOutputRGB(outputFileName);
-	seg_.SetOutputLabel(outputLabelFileName);
-	seg_.SetParameters(params);
-	if(argc >= 7)
-		seg_.SetNumberOfIterations(static_cast<unsigned int>(atoi(argv[6])));
-	if(argc == 8)
-		seg_.SetBestFitting(atoi(argv[7]));
-
-	seg_.Run();
-
-	return 0;
-}
diff --git a/src/Applications/FLSASegmentation.cxx b/src/Applications/FLSASegmentation.cxx
deleted file mode 100644
index b908daae808336e52b0770d72824523708923a89..0000000000000000000000000000000000000000
--- a/src/Applications/FLSASegmentation.cxx
+++ /dev/null
@@ -1,61 +0,0 @@
-/*=========================================================================
-
-  Program: Example of application of the Generic Region Merging Library
-  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.
-
-=========================================================================*/
-#include <iostream>
-// Your algorithm header file
-#include "fls-algorithm.h"
-
-int main(int argc, char **argv)
-{
-	if(argc < 5 || argc > 7)
-	{
-		std::cout << "###############      Manual     #############################\n\n\n"
-				<< "Usage " << argv[0] << "\n\n\n"
-				<< " -- [input image (tif, jpg, png) <mandatory>]\n"
-				<< " -- [output rgb image (tif, jpg, png) <mandatory>]\n"
-				<< " -- [output label image (tif) <mandatory>]\n"
-				<< " -- [lambda threshold <mandatory>]\n"
-				<< " -- [number of iterations using the Local Mutual Best Fitting heuristic (70 by default) <optional>]\n"
-				<< " -- [activation of the Best Fitting heuristic (value = 0 (desactivated) or 1 (activated) \
-					- activated by default) <optional>]\n\n\n"
-				<< "################################################################\n";
-		return 1;
-	}
-
-	const std::string inputFileName = argv[1];
-	const std::string outputFileName = argv[2];
-	const std::string outputLabelFileName = argv[3];
-	float lambda = atof(argv[4]);
-
-	typedef otb::VectorImage<float, 2> ImageType;
-	typedef FLSAlgorithmRM<ImageType> SegmenterType;
-
-	FLSParams params = lambda;
-	SegmenterType seg_;
-	seg_.SetInput(inputFileName);
-	seg_.SetOutputRGB(outputFileName);
-	seg_.SetOutputLabel(outputLabelFileName);
-	seg_.SetParameters(params);
-
-	if(argc >= 6)
-		seg_.SetNumberOfIterations(static_cast<unsigned int>(atoi(argv[5])));
-	if(argc == 7)
-		seg_.SetBestFitting(atoi(argv[6]));
-
-	seg_.Run();
-
-	return 0;
-}
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
deleted file mode 100644
index bf76e2092aab7b957ac27218b1e1f41a2ffb552a..0000000000000000000000000000000000000000
--- a/src/CMakeLists.txt
+++ /dev/null
@@ -1,50 +0,0 @@
-#=========================================================================
-
-#  Program: 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(GRM)
-
-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)
-
-#Check compiler version
-if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU")
-    # require at least gcc 4.8
-    if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 4.8)
-        message(FATAL_ERROR "GCC version must be at least 4.8!")
-    endif()
-elseif ("${CMAKE_CXX_COMPILER_ID}" STREQUAL "Clang")
-    # require at least clang 3.2
-    if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 3.2)
-        message(FATAL_ERROR "Clang version must be at least 3.2!")
-    endif()
-else()
-    message(WARNING "You are using an unsupported compiler! Compilation has only been tested with Clang and GCC.")
-endif()
-
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -std=c++11 -fpermissive")
-
-include_directories(${CMAKE_CURRENT_SOURCE_DIR}/Library)
-
-add_subdirectory(Library)
-add_subdirectory(Applications)
diff --git a/src/Library/baatz-algorithm.h b/src/Library/baatz-algorithm.h
deleted file mode 100644
index d428b5590deeedfe6ee25911e6b8e58c7a1d7187..0000000000000000000000000000000000000000
--- a/src/Library/baatz-algorithm.h
+++ /dev/null
@@ -1,64 +0,0 @@
-/*=========================================================================
-
-  Program: Implementation of the Baatz & Schäpe criterion using the GRM library
-  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 __BAATZ_ALGORITHM_H
-#define __BAATZ_ALGORITHM_H
-
-#include "grm-region.h"
-#include "grm-algorithm.h"
-
-class BaatzRegion : public grm::RegionBase<BaatzRegion>
-{
-	public:
-		float m_Area;
-		float m_Perimeter;
-		std::vector<float> m_Means;
-		std::vector<float> m_SquareMeans;
-		std::vector<float> m_SpectralSum;
-		std::vector<float> m_Std;
-};
-
-struct BaatzParams
-{
-	float m_SpectralWeight;
-	float m_ShapeWeight;
-	float m_Scale;
-};
-
-template<class TInputImage>
-class BaatzAlgorithmRM : public grm::RegionMergingAlgorithm<TInputImage, 
-														BaatzRegion,
-														BaatzParams>
-{
-	public:
-		// Mandatory typedef
-		typedef grm::RegionMergingAlgorithm<TInputImage, BaatzRegion, BaatzParams> Superclass;
-		typedef BaatzRegion RegionType;
-		typedef std::shared_ptr<RegionType> RegionPointerType;
-		typedef std::vector<RegionPointerType> RegionListType;
-		typedef BaatzParams ParameterType;
-		typedef grm::Interface<RegionType> InterfaceType;
-		typedef typename Superclass::MergingCostFunction MergingCostFunction;
-
-		BaatzAlgorithmRM(){}
-		virtual ~BaatzAlgorithmRM(){}
-
-		void InitFromImage();
-		void UpdateAttribute(RegionPointerType, RegionPointerType);
-		virtual void Segmentation();
-};
-#include "baatz-algorithm.txx"
-#endif
diff --git a/src/Library/baatz-algorithm.txx b/src/Library/baatz-algorithm.txx
deleted file mode 100644
index 6c6c25341ed56b39875f00a9225c30ed0c9b6f0e..0000000000000000000000000000000000000000
--- a/src/Library/baatz-algorithm.txx
+++ /dev/null
@@ -1,226 +0,0 @@
-/*=========================================================================
-
-  Program: Implementation of the Baatz & Schäpe criterion using the GRM library
-  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 __BAATZ_ALGORITHM_TXX
-#define __BAATZ_ALGORITHM_TXX
-
-#include "baatz-algorithm.h"
-#include <otbVectorImage.h>
-#include <otbImageFileReader.h>
-#include <itkImageRegionIterator.h>
-
-template<class TInputImage>
-void 
-BaatzAlgorithmRM<TInputImage>::InitFromImage()
-{
-	typedef otb::ImageFileReader<TInputImage> InputImageFileReaderType;
-	typedef itk::ImageRegionIterator<TInputImage> InputImageIteratorType;
-
-	// Read the input image
-	auto input_reader = InputImageFileReaderType::New();
-	input_reader->SetFileName(this->m_Input);
-	input_reader->Update();
-
-	// Retrieve useful information from the input image
-	auto regionToProcess = input_reader->GetOutput()->GetLargestPossibleRegion();
-	unsigned int rows = regionToProcess.GetSize()[1];
-	unsigned int cols = regionToProcess.GetSize()[0];
-	unsigned int bands = input_reader->GetOutput()->GetNumberOfComponentsPerPixel();
-
-	// Step 1: For efficiency the user is advised to allocate the initial number of regions
-	this->m_RegionList.reserve(rows*cols);
-
-	// Step 2: the user has to intialize the specific attributes of the initial regions
-	InputImageIteratorType it(input_reader->GetOutput(), regionToProcess);
-	for(it.GoToBegin(); !it.IsAtEnd(); ++it)
-	{
-		auto new_region = std::make_shared<BaatzRegion>();
-		new_region->m_Area = 1;
-		new_region->m_Perimeter = 4;
-		new_region->m_Means.reserve(bands);
-		new_region->m_SquareMeans.reserve(bands);
-		new_region->m_SpectralSum.reserve(bands);
-		new_region->m_Std.reserve(bands);
-		for(unsigned int b = 0; b<bands; ++b)
-		{
-			new_region->m_Means.push_back(it.Get()[b]);
-			new_region->m_SquareMeans.push_back(pow(it.Get()[b], 2));
-			new_region->m_SpectralSum.push_back(it.Get()[b]);
-			new_region->m_Std.push_back(0);
-		}
-		// The user has to push the initial regions
-		this->m_RegionList.push_back(new_region);
-	}
-
-	// Step 3: The neighborhhod and some internal attributes have to be initialized
-	// the user has to call the method InitRegions
-	this->InternalInit(rows, cols);
-}
-
-// Let assume that regions r1 and r2 have to be merged and r2 is merged into r1
-// The user has to describe how the specific attributes of r1 have to be updated
-// by defining the mandatory method UpdateAttribute(RegionPointerType, RegionPointerType)
-template<class TInputImage>
-void 
-BaatzAlgorithmRM<TInputImage>::UpdateAttribute(RegionPointerType r1, RegionPointerType r2)
-{
-	float a1 = r1->m_Area, a2 = r2->m_Area, a_tot = a1 + a2;
-	for(unsigned int b=0; b<r1->m_Means.size(); ++b)
-	{
-		r1->m_Means[b] = ((a1*r1->m_Means[b]) + (a2*r2->m_Means[b])) / a_tot;
-		r1->m_SquareMeans[b] += r2->m_SquareMeans[b];
-		r1->m_SpectralSum[b] += r2->m_SpectralSum[b];
-		r1->m_Std[b] = std::sqrt((r1->m_SquareMeans[b] - 
-								2*r1->m_Means[b]*r1->m_SpectralSum[b] + 
-								a_tot*r1->m_Means[b]*r1->m_Means[b]) / a_tot); 
-	}
-	r1->m_Area += r2->m_Area;
-	r1->m_Perimeter += r2->m_Perimeter - 2 * r1->GetConnections(r2);
-}
-
-// This method is obviously mandatory but its implementation depends on some
-// subtilities of your algorithm
-// What is common is the the iterative merging process of the region, what can be different
-// is the stopping criterion and the behavior of the threshold
-template<class TInputImage>
-void 
-BaatzAlgorithmRM<TInputImage>::Segmentation()
-{
-	bool prev_merged = true;
-	unsigned int step = 0;
-	float threshold = this->m_Parameters.m_Scale*this->m_Parameters.m_Scale;
-	float spec_w = this->m_Parameters.m_SpectralWeight;
-	float shap_w = this->m_Parameters.m_ShapeWeight;
-
-	MergingCostFunction cost_func = [&](RegionPointerType r1, RegionPointerType r2)->double
-	{
-		unsigned int bands = r1->m_Means.size();
-		float a1 = r1->m_Area, a2 = r2->m_Area, a_tot = a1 + a2;
-
-		// Compute the spectral merging cost
-		double spec_cost = 0.0;
-		float mean, square_mean, sum, std, tmp_cost;
-		for(unsigned int b = 0; b < bands; b++)
-		{
-			mean = ((a1*r1->m_Means[b]) + (a2*r2->m_Means[b])) / a_tot;
-			square_mean = r1->m_SquareMeans[b] + r2->m_SquareMeans[b];
-			sum = r1->m_SpectralSum[b] + r2->m_SpectralSum[b];
-			std = std::sqrt((square_mean - 2*mean*sum + a_tot*mean*mean)/a_tot);
-			tmp_cost = a_tot * std - a1*r1->m_Std[b] - a2*r2->m_Std[b];
-			spec_cost += tmp_cost;
-		}
-		spec_cost *= spec_w;
-
-		if(spec_cost < threshold)
-		{
-			double shape_cost, smooth_f, compact_f;
-
-			// Compute the shape merging cost
-			float p1 = r1->m_Perimeter, p2 = r2->m_Perimeter, 
-					p3 = p1 + p2 - 2*r1->GetConnections(r2);
-
-			auto mmbox = r1->GetResultingBbox(r2);
-			float bb1 = 2*r1->GetBboxSize(0) + 2*r1->GetBboxSize(1),
-				bb2 = 2*r2->GetBboxSize(0) + 2*r2->GetBboxSize(1),
-				bb3 = 2*mmbox[2] + 3*mmbox[3];
-
-			smooth_f = a_tot*p3/bb3 - a1*p1/bb1 - a2*p2/bb2;
-			compact_f = a_tot*p3/std::sqrt(a_tot) - a1*p1/std::sqrt(a1) - a2*p2/std::sqrt(a2);
-
-			shape_cost = shap_w * compact_f + (1-shap_w) * smooth_f;
-
-			return (spec_cost + (1-spec_w)*shape_cost);
-		}
-		else
-			return spec_cost;
-	};
-
-	while(prev_merged && step < this->m_NumberOfIterations)
-	{
-		std::cout << "." << std::flush;
-		prev_merged = false;
-		++step;
-
-		this->m_RMHandler.UpdateMergingCosts(this->m_RegionList, cost_func);
-
-		// Do a loop on the regions of the graph
-		for(auto& r : this->m_RegionList)
-		{
-			// For each explored region, we determine if the LMBF is met
-			auto ref_region = this->m_RMHandler.CheckLMBF(r, threshold);
-			// If it holds then the pointer to ref_region is not null
-			if(ref_region != nullptr)
-			{
-				// Process to merge the regions
-				// Step 1: Merge the specific attributes
-				UpdateAttribute(ref_region, ref_region->GetClosestNeighbor());
-				// Step 2: Internal update (mandatory)
-				this->m_RMHandler.Update(ref_region, ref_region->GetClosestNeighbor());
-
-				prev_merged = true;
-			}
-		}
-
-		// After one iteration, you have to remove the expired regions, i.e. regions
-		// who have merged into a larger one
-		this->m_RMHandler.RemoveExpiredVertices(this->m_RegionList);
-	}
-
-	if(prev_merged && this->m_BestFitting)
-	{
-		while(prev_merged)
-		{
-			std::cout << "." << std::flush;
-			prev_merged = false;
-
-			this->m_RMHandler.UpdateMergingCosts(this->m_RegionList, cost_func);
-
-			for(auto& r : this->m_RegionList)
-			{
-				// For each explored region, we determine if the LMBF is met
-				auto ref_region = this->m_RMHandler.CheckBF(r, threshold);
-				// If it holds then the pointer to ref_region is not null
-				if(ref_region != nullptr)
-				{
-					if(ref_region == r->GetClosestNeighbor())
-					{
-						// Process to merge the regions
-						// Step 1: Merge the specific attributes
-						UpdateAttribute(ref_region, r);
-						// Step 2: Internal update (mandatory)
-						this->m_RMHandler.Update(ref_region, r);
-					}
-					else
-					{
-						// Process to merge the regions
-						// Step 1: Merge the specific attributes
-						UpdateAttribute(ref_region, ref_region->GetClosestNeighbor());
-						// Step 2: Internal update (mandatory)
-						this->m_RMHandler.Update(ref_region, ref_region->GetClosestNeighbor());
-					}
-					prev_merged = true;
-				}
-			}
-
-			// After one iteration, you have to remove the expired regions, i.e. regions
-			// who have merged into a larger one
-			this->m_RMHandler.RemoveExpiredVertices(this->m_RegionList);
-		}
-	}
-	std::cout << "\n";
-}
-
-#endif
diff --git a/src/Library/contour-encoder.cxx b/src/Library/contour-encoder.cxx
deleted file mode 100644
index 2483a829dd3e039cdd6a0449de38121f9fdb2ebf..0000000000000000000000000000000000000000
--- a/src/Library/contour-encoder.cxx
+++ /dev/null
@@ -1,675 +0,0 @@
-/*=========================================================================
-
-  Program:  Shape Encoder librarie of OBIA (Object Based Image Analysis)
-  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.
-
-=========================================================================*/
-#include "contour-encoder.h"
-
-namespace shpenc
-{
-	typename ShapeEncoder::Contour
-	ShapeEncoder::EncodeContour(long unsigned int curr_id, 
-								Contour&     curr_contour,
-								long unsigned int adj_id,
-								Contour&     adj_contour,
-								BoundingBox& merged_bbox)
-	{
-		// To determine the future contour pixels
-		std::vector<short> matrix_bbox(merged_bbox[2]*merged_bbox[3], 0);
-
-		// Generate contour pixels of the merged segments
-		auto curr_v = GeneratePixels(curr_id, curr_contour);
-		auto adj_v = GeneratePixels(adj_id, adj_contour);
-		curr_v.insert(curr_v.end(), adj_v.begin(), adj_v.end());
-
-		// A case in the matrix_bbox filled by a pixel is marked true
-		for(auto& pix : curr_v)
-			matrix_bbox[ToBBoxId(pix, merged_bbox)] = 1;
-
-		// Remove internal pixels
-		RemoveInternalPixels(matrix_bbox, merged_bbox[3], merged_bbox[2]);
-
-		// Encode the new contour
-		return EncodeNewContour(curr_id, merged_bbox, matrix_bbox);
-	}
-
-	std::vector<long unsigned int> ShapeEncoder::GeneratePixels(long unsigned int id, Contour& contour)
-	{
-		std::vector<long unsigned int> pixels;
-		pixels.push_back(id);
-		
-		short curr_move, prev_move;
-		
-		if(contour.size() > 4)
-		{
-			for(auto m = contour.begin()+1; m != contour.end(); m++)
-			{
-				curr_move = (*m)[0] + 2*(*m)[1];
-				prev_move = (*(m-1))[0] + 2*(*(m-1))[1];
-
-				if(curr_move == 0)
-				{
-					if(prev_move == 0)
-					{
-						id-=m_Cols;
-						pixels.push_back(id);
-					}
-					else if(prev_move == 1) // 1
-					{
-						id++;
-						pixels.push_back(id);
-						id-=m_Cols;
-						pixels.push_back(id);
-					}
-				}
-				else if(curr_move == 1)
-				{
-					if(prev_move == 1)
-					{
-						id++;
-						pixels.push_back(id);
-					}
-					else if(prev_move == 2)
-					{
-						id+=m_Cols;
-						pixels.push_back(id);
-						id++;
-						pixels.push_back(id);
-					}
-				}
-				else if(curr_move == 2)
-				{
-					if(prev_move == 3)
-					{
-						id--;
-						pixels.push_back(id);
-						id+=m_Cols;
-						pixels.push_back(id);
-					}
-					else if(prev_move == 2)
-					{
-						id+=m_Cols;
-						pixels.push_back(id);
-					}
-				}
-				else if(curr_move == 3)
-				{
-					if(prev_move==0)
-					{
-						id-=m_Cols;
-						pixels.push_back(id);
-						id--;
-						pixels.push_back(id);
-					}
-					else if(prev_move == 3)
-					{
-						id--;
-						pixels.push_back(id);
-					}
-				}
-			}
-			// Remove duplicated pixels
-			std::sort(pixels.begin(), pixels.end());
-			auto it = std::unique(pixels.begin(), pixels.end());
-  			pixels.resize(std::distance(pixels.begin(), it));
-		}
-		return pixels;
-	}
-
-	unsigned int ShapeEncoder::ToBBoxId(long unsigned int id, BoundingBox& bounding_box)
-	{
-		unsigned int img_row = id / m_Cols;
-		unsigned int img_col = id % m_Cols;
-		unsigned int bbox_row = img_row - bounding_box[1];
-		unsigned int bbox_col = img_col - bounding_box[0];
-		return bbox_row * bounding_box[2] + bbox_col;
-	}
-
-	long unsigned int ShapeEncoder::ToImgId(unsigned int id, BoundingBox& bounding_box)
-	{
-		unsigned int bbox_row = id / bounding_box[2];
-		unsigned int bbox_col = id % bounding_box[2];
-		unsigned int img_row = bbox_row + bounding_box[1];
-		unsigned int img_col = bbox_col + bounding_box[0];
-
-		return img_row * m_Cols + img_col;
-	}
-
-	void ShapeEncoder::RemoveInternalPixels(std::vector<short>& matrix_bbox, 
-							  				  unsigned int bbox_rows, 
-							                  unsigned int bbox_cols)
-	{
-		int mat_idx;
-		bool is_internal;
-		for(int r = 1; r<bbox_rows-1; r++)
-		{
-			for(int c = 1; c<bbox_cols-1; c++)
-			{
-				mat_idx = r*bbox_cols+c;
-				is_internal = true;
-				if(matrix_bbox[mat_idx] == 1)
-				{
-					auto neighbors = mask::SquareNeighborhood(mat_idx, bbox_cols, bbox_rows);
-					for(short j=0; j<8; j++)
-					{
-						if(matrix_bbox[neighbors[j]] == 0)
-						{
-							is_internal = false;
-							break;
-						}
-					}
-					if(is_internal)
-						matrix_bbox[mat_idx] = 2;
-				}
-			}	
-		}
-		for(int i = 0; i<bbox_rows*bbox_cols; i++)
-			if(matrix_bbox[i] == 2) matrix_bbox[i] = 0;
-	}
-
-	typename ShapeEncoder::Contour 
-	ShapeEncoder::EncodeNewContour(long unsigned int curr_id,
-									BoundingBox& merged_bbox,
-									std::vector<short>& matrix_bbox)
-	{
-		short direction, pos;              // position on one of the 4 corners of the pixel
-		std::array<short, 4> sight;        // direction of the field of view
-
-		unsigned int start_mat_id = ToBBoxId(curr_id, merged_bbox);
-		unsigned int curr_mat_id = start_mat_id;
-
-		unsigned int matrix_bbox_rows = merged_bbox[3];
-		unsigned int matrix_bbox_cols = merged_bbox[2];
-		
-		Contour curr_contour; // Reset the previous contour
-
-		// Build the first move
-		auto neighbors = mask::CrossNeighborhood(curr_mat_id, matrix_bbox_cols, matrix_bbox_rows);
-
-		// Only 2 cases: neighbor at the right or at the bottom
-		if(matrix_bbox[neighbors[1]] == 1) // right
-		{
-			pos = 1;       // top left corner of the pixel
-			curr_mat_id++; // go to the pixel at the right
-			direction = 1; // direction is along the right
-			curr_contour.push_back(1); // Add the move to the right
-		}
-		else if(matrix_bbox[neighbors[2]] == 1) // bottom
-		{
-			pos = 2; // top right corner of the pixel
-			curr_mat_id += matrix_bbox_cols; // go to the pixel at the bottom
-			direction = 2; // direction is along the bottom						
-			curr_contour.push_back(1); // add move to the right
-			curr_contour.push_back(2); // add move to the bottom	
-		}
-
-		// Keep going this same reasonning until we reach the start pixel (start_mat_id)
-		while(curr_mat_id != start_mat_id)
-		{
-			UpdateSight(sight, direction);
-			auto neighbors = mask::CrossNeighborhood(curr_mat_id, matrix_bbox_cols, matrix_bbox_rows);
-
-			// All the cases are possible: top, right, bottom, left
-			for(short d=0; d<4; d++)
-			{
-				if(neighbors[sight[d]] > -1)
-				{
-					if(matrix_bbox[neighbors[sight[d]]] > 0)
-					{
-						if(sight[d] == 0)	
-						{
-							EncodeContourForTopPixel(curr_mat_id, direction, pos, curr_contour, matrix_bbox_cols);
-							break;
-						}
-						else if(sight[d] == 1)
-						{
-							EncodeContourForRightPixel(curr_mat_id, direction, pos, curr_contour, matrix_bbox_cols);
-							break;
-						}
-						else if(sight[d] == 2)
-						{
-							EncodeContourForBottomPixel(curr_mat_id, direction, pos, curr_contour, matrix_bbox_cols);
-							break;
-						}
-						else if(sight[d] == 3)
-						{
-							EncodeContourForLeftPixel(curr_mat_id, direction, pos, curr_contour, matrix_bbox_cols);
-							break;
-						}
-					}
-				}
-			}
-
-			// It is possible to reach the start pixel whithout finishing to encode the whole contour.
-			// 00111
-			// 11100
-			// In this case the direction value is always 3, the position value is always 3
-			// We have to check if there is a neighbor at left of the current direction
-			if(curr_mat_id == start_mat_id)
-			{
-				if(pos == 3 && direction == 3)
-				{
-					UpdateSight(sight, direction);
-					auto neighbors = mask::CrossNeighborhood(curr_mat_id, matrix_bbox_cols, matrix_bbox_rows);
-
-					if(neighbors[sight[0]] > -1)
-					{
-						if(matrix_bbox[neighbors[sight[0]]] == 1)
-						{
-							pos = 3;
-							direction = 2;
-							curr_contour.push_back(2);
-							curr_mat_id += matrix_bbox_cols;
-						}
-					}
-				}
-			}		
-		}
-
-		// We reach the start pixel but maybe we did not encode the whole contour
-		if(pos == 3)
-		{
-			curr_contour.push_back(3);
-			curr_contour.push_back(0);
-		}
-		else if(pos == 4)
-			curr_contour.push_back(0);
-
-		return curr_contour;
-	}
-
-	void ShapeEncoder::UpdateSight(std::array<short, 4>& sight, short direction)
-	{
-		sight[0] = (direction + 3) % 4; // Look at the left
-		sight[1] = direction; // Look in front of
-		sight[2] = (direction + 1) % 4; // Look at the right
-		sight[3] = (direction + 2) % 4; // Look behind 
-	}
-
-	void ShapeEncoder::EncodeContourForTopPixel(unsigned int& curr_mat_id,
-												  short& direction,
-												  short& pos,
-									              Contour& curr_contour,
-												  const unsigned int matrix_bbox_cols)
-	{
-		if(pos == 1)
-		{
-			if(direction==0)
-				pos = 4;
-			else if(direction == 1)
-			{
-				pos = 1;
-				curr_contour.push_back(0);
-			}
-		}
-		else if(pos == 2)
-		{
-			pos = 4;
-			curr_contour.push_back(2);
-			curr_contour.push_back(3);
-			curr_contour.push_back(0);
-		}
-		else if(pos == 3)
-		{
-			pos = 4;
-			curr_contour.push_back(3);
-			curr_contour.push_back(0);
-		}
-		else if(pos == 4)
-			curr_contour.push_back(0);
-
-
-		curr_mat_id -= matrix_bbox_cols;
-		direction = 0;
-	}
-	void ShapeEncoder::EncodeContourForRightPixel(unsigned int& curr_mat_id,
-												    short& direction,
-												    short& pos,
-									                Contour& curr_contour,
-												    const unsigned int matrix_bbox_cols)
-	{
-		if(pos == 1)
-			curr_contour.push_back(1);
-		else if(pos == 2)
-		{
-			if(direction == 1)
-				pos = 1;
-			else if(direction == 2)
-				curr_contour.push_back(1);
-		}
-		else if(pos == 3)
-		{
-			pos = 1;
-			curr_contour.push_back(3);
-			curr_contour.push_back(0);
-			curr_contour.push_back(1);
-		}
-		else if(pos == 4)
-		{
-			pos = 1;
-			curr_contour.push_back(0);
-			curr_contour.push_back(1);
-		}
-
-		curr_mat_id++;
-		direction = 1;	
-	}
-	void ShapeEncoder::EncodeContourForBottomPixel(unsigned int& curr_mat_id,
-												     short& direction,
-												     short& pos,
-									                 Contour& curr_contour,
-												     const unsigned int matrix_bbox_cols)
-	{
-		if(pos == 1)
-		{
-			pos = 2;
-			curr_contour.push_back(1);
-			curr_contour.push_back(2);
-		}
-		else if(pos == 2)
-			curr_contour.push_back(2);
-		else if(pos == 3)
-		{
-			if(direction == 2)
-				pos = 2;
-			else if(direction == 3)
-				curr_contour.push_back(2);
-		}
-		else if(pos == 4)
-		{
-			pos = 2;
-			curr_contour.push_back(0);
-			curr_contour.push_back(1);
-			curr_contour.push_back(2);
-		}
-
-		curr_mat_id += matrix_bbox_cols;
-		direction = 2;			
-	}
-	void ShapeEncoder::EncodeContourForLeftPixel(unsigned int& curr_mat_id,
-												   short& direction,
-												   short& pos,
-									               Contour& curr_contour,
-												   const unsigned int matrix_bbox_cols)
-	{
-		if(pos == 1)
-		{
-			pos = 3;
-			curr_contour.push_back(1);
-			curr_contour.push_back(2);
-			curr_contour.push_back(3);
-		}
-		else if(pos == 2)
-		{
-			pos = 3;
-			curr_contour.push_back(2);
-			curr_contour.push_back(3);
-		}
-		else if(pos == 3)
-			curr_contour.push_back(3);
-		else if(pos == 4)
-		{
-			if(direction == 3)
-				pos = 3;
-			else if(direction == 0)
-			{
-				pos = 4;
-				curr_contour.push_back(3);
-			}
-		}
-
-		curr_mat_id --;
-		direction = 3;
-	}
-
-
-	bool ShapeEncoder::GetCollisionAtNorth(unsigned int idx,
-											const std::vector<short>& bbox_array,
-											unsigned int bbox_cols,
-											unsigned int bbox_rows)
-	{
-		assert(idx >= 0 && idx < bbox_cols * bbox_rows);
-		long int tmp_idx = idx - bbox_cols;
-
-		if(bbox_array[tmp_idx] == 1)
-			return true;
-
-		while(tmp_idx > 0)
-		{
-			if(bbox_array[tmp_idx] == 1)
-				return true;
-			tmp_idx -= bbox_cols;
-		}
-		return false;
-	}
-
-	bool ShapeEncoder::GetCollisionAtNorthEast(unsigned int idx,
-											const std::vector<short>& bbox_array,
-											unsigned int bbox_cols,
-											unsigned int bbox_rows)
-	{
-		assert(idx >= 0 && idx < bbox_cols * bbox_rows);
-		long int tmp_idx = idx - bbox_cols + 1;
-
-		if(bbox_array[tmp_idx] == 1)
-			return true;
-
-		while(tmp_idx > 0 && (tmp_idx % bbox_cols) != 0)
-		{
-			if(bbox_array[tmp_idx] == 1)
-				return true;
-			tmp_idx = tmp_idx - bbox_cols + 1;
-		}
-		return false;
-	}
-
-	bool ShapeEncoder::GetCollisionAtEast(unsigned int idx,
-										const std::vector<short>& bbox_array,
-										unsigned int bbox_cols,
-										unsigned int bbox_rows)
-	{
-		assert(idx >= 0 && idx < bbox_cols * bbox_rows);
-		long int tmp_idx = idx + 1;
-
-
-		if(bbox_array[tmp_idx] == 1)
-			return true;
-
-		while((tmp_idx % bbox_cols) != 0)
-		{
-			if(bbox_array[tmp_idx] == 1)
-				return true;
-			++tmp_idx;
-		}
-		return false;
-	}
-
-	bool ShapeEncoder::GetCollisionAtSouthEast(unsigned int idx,
-										const std::vector<short>& bbox_array,
-										unsigned int bbox_cols,
-										unsigned int bbox_rows)
-	{
-		assert(idx >= 0 && idx < bbox_cols * bbox_rows);
-		long int tmp_idx = idx + 1 + bbox_cols;
-
-		if(bbox_array[tmp_idx] == 1)
-			return true;
-
-		while(tmp_idx < (bbox_cols * bbox_rows) &&
-			(tmp_idx % bbox_cols) != 0)
-		{
-			if(bbox_array[tmp_idx] == 1)
-				return true;
-			tmp_idx = tmp_idx + 1 + bbox_cols;
-		}
-		return false;
-	}
-
-	bool ShapeEncoder::GetCollisionAtSouth(unsigned int idx,
-										const std::vector<short>& bbox_array,
-										unsigned int bbox_cols,
-										unsigned int bbox_rows)
-	{
-		assert(idx >= 0 && idx < bbox_cols * bbox_rows);
-		long int tmp_idx = idx + bbox_cols;
-
-		if(bbox_array[tmp_idx] == 1)
-			return true;
-
-		while(tmp_idx < (bbox_cols * bbox_rows))
-		{
-			if(bbox_array[tmp_idx] == 1)
-				return true;
-			tmp_idx += bbox_cols;
-		}
-		return false;
-	}
-
-	bool ShapeEncoder::GetCollisionAtSouthWest(unsigned int idx,
-											const std::vector<short>& bbox_array,
-											unsigned int bbox_cols,
-											unsigned int bbox_rows)
-	{
-		assert(idx >= 0 && idx < bbox_cols * bbox_rows);
-		long int tmp_idx = idx + bbox_cols - 1;
-
-		if(bbox_array[tmp_idx] == 1)
-			return true;
-
-		while(tmp_idx < (bbox_cols * bbox_rows) &&
-				(tmp_idx % bbox_cols) != bbox_cols - 1)
-		{
-			if(bbox_array[tmp_idx] == 1)
-				return true;
-			tmp_idx = tmp_idx + bbox_cols - 1;
-		}
-		return false;
-	}
-
-	bool ShapeEncoder::GetCollisionAtWest(unsigned int idx,
-											const std::vector<short>& bbox_array,
-											unsigned int bbox_cols,
-											unsigned int bbox_rows)
-	{
-		assert(idx >= 0 && idx < bbox_cols * bbox_rows);
-		long int tmp_idx = idx - 1;
-
-		if(bbox_array[tmp_idx] == 1)
-			return true;
-
-		while((tmp_idx % bbox_cols) != bbox_cols - 1)
-		{
-			if(bbox_array[tmp_idx] == 1)
-				return true;
-			tmp_idx -= 1;
-		}
-		return false;
-	}
-
-	bool ShapeEncoder::GetCollisionAtNorthWest(unsigned int idx,
-											const std::vector<short>& bbox_array,
-											unsigned int bbox_cols,
-											unsigned int bbox_rows)
-	{
-		assert(idx >= 0 && idx < bbox_cols * bbox_rows);
-		long int tmp_idx = idx - bbox_cols - 1;
-
-		if(bbox_array[tmp_idx] == 1)
-			return true;
-
-		while(tmp_idx > 0 && (tmp_idx % bbox_cols) != bbox_cols - 1)
-		{
-			if(bbox_array[tmp_idx] == 1)
-				return true;
-			tmp_idx = tmp_idx - bbox_cols - 1;
-		}
-		return false;
-	}
-
-	bool ShapeEncoder::IsInternalPixel(unsigned int curr_idx, std::vector<short>& bbox_array, 
-										unsigned int bbox_size_x, unsigned int bbox_size_y)
-	{
-		if(!GetCollisionAtNorth(curr_idx, bbox_array, bbox_size_x, bbox_size_y))
-			return false;
-		if(!GetCollisionAtNorthEast(curr_idx, bbox_array, bbox_size_x, bbox_size_y))
-			return false;
-		if(!GetCollisionAtEast(curr_idx, bbox_array, bbox_size_x, bbox_size_y))
-			return false;
-		if(!GetCollisionAtSouthEast(curr_idx, bbox_array, bbox_size_x, bbox_size_y))
-			return false;
-		if(!GetCollisionAtSouth(curr_idx, bbox_array, bbox_size_x, bbox_size_y))
-			return false;
-		if(!GetCollisionAtSouthWest(curr_idx, bbox_array, bbox_size_x, bbox_size_y))
-			return false;
-		if(!GetCollisionAtWest(curr_idx, bbox_array, bbox_size_x, bbox_size_y))
-			return false;
-		if(!GetCollisionAtNorthWest(curr_idx, bbox_array, bbox_size_x, bbox_size_y))
-			return false;
-
-		return true;
-	}
-
-	std::vector<long unsigned int> ShapeEncoder::GenerateAllPixels(long unsigned int id, 
-																	Contour& contour, 
-																	BoundingBox& bounding_box)
-	{
-		// Generate the border pixels id in the coordinnate system of the image
-		auto border_pixels = GeneratePixels(id,contour);
-		// Build the matrix which corresponds to the bounding box
-		const unsigned int bbox_array_size = bounding_box[2]*bounding_box[3];
-		const unsigned int bbox_size_x = bounding_box[2];
-		const unsigned int bbox_size_y = bounding_box[3];
-		std::vector<short> bbox_array;
-		bbox_array.assign(bbox_array_size, 0);
-
-		// Fill the bbox_array with the border pixels
-		for(auto& bp : border_pixels)
-			bbox_array[ToBBoxId(bp, bounding_box)] = 1;
-
-		// For each case of the bbox_array we determine if it is inside or outside of the region
-		std::vector<long unsigned int> pixels;
-		pixels.insert(pixels.begin(), border_pixels.begin(), border_pixels.end());
-		bool internal = true;
-		unsigned int tmp_idx, curr_idx;
-		for(unsigned int r = 1; r < bbox_size_y - 1; r++)
-		{
-			for(unsigned int c = 1; c < bbox_size_x - 1; c++)
-			{
-				curr_idx = r*bbox_size_x + c;
-				if(bbox_array[curr_idx] != 1)
-				{
-					auto internal = IsInternalPixel(curr_idx, bbox_array, bbox_size_x, bbox_size_y);
-					bbox_array[curr_idx] = 1;
-					if(internal)
-					{
-						pixels.push_back(ToImgId(curr_idx, bounding_box));
-					}
-				}
-			}
-		}
-
-		return pixels;
-	}
-
-}
-
-
-
-
-
-
-
-
-
-
diff --git a/src/Library/contour-encoder.h b/src/Library/contour-encoder.h
deleted file mode 100644
index 6bedd89630872c57ff088d0d4a3515ce853a9ec5..0000000000000000000000000000000000000000
--- a/src/Library/contour-encoder.h
+++ /dev/null
@@ -1,132 +0,0 @@
-/*=========================================================================
-
-  Program: Shape Encoder librarie of OBIA (Object Based Image Analysis)
-  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 __SHPENC_ShapeEncoder_H
-#define __SHPENC_ShapeEncoder_H
-
-#include "macro-generator.h"
-#include "mask-neighborhood.h"
-#include <cassert>
-#include <iostream>
-#include <vector>
-#include <bitset>
-#include <algorithm>
-
-namespace shpenc
-{
-	class ShapeEncoder
-	{
-		public:
-
-			typedef std::bitset<2>             Move;
-			typedef typename std::vector<Move> Contour;
-			typedef std::array<unsigned int, 4> BoundingBox;
-
-			ShapeEncoder(){}
-			~ShapeEncoder(){}
-
-			SetMacro(unsigned int, Rows);
-			SetMacro(unsigned int, Cols);
-			GetMacro(unsigned int, Rows);
-			GetMacro(unsigned int, Cols);
-
-			Contour EncodeContour(long unsigned int, Contour&,
-								long unsigned int, Contour&,
-								BoundingBox&);
-			std::vector<long unsigned int> GeneratePixels(long unsigned int, Contour&);
-			std::vector<long unsigned int> GenerateAllPixels(long unsigned int, Contour&, BoundingBox&);
-
-		private:
-
-			// Internal methods
-			unsigned int ToBBoxId(long unsigned int, BoundingBox&);
-			long unsigned int ToImgId(unsigned int, BoundingBox&);
-			void RemoveInternalPixels(std::vector<short>&,
-									  unsigned int, unsigned int);
-			Contour EncodeNewContour(long unsigned int,
-									BoundingBox&,
-									std::vector<short>&);
-			void UpdateSight(std::array<short, 4>&, short);
-			void EncodeContourForTopPixel(unsigned int&,
-										  short&,
-										  short&,
-									      Contour&,
-										  const unsigned int);
-			void EncodeContourForRightPixel(unsigned int&,
-										    short&,
-										    short&,
-									        Contour&,
-										    const unsigned int);
-			void EncodeContourForBottomPixel(unsigned int&,
-										     short&,
-										     short&,
-									         Contour&,
-										     const unsigned int);
-			void EncodeContourForLeftPixel(unsigned int&,
-										   short&,
-										   short&,
-									       Contour&,
-										   const unsigned int);
-
-			bool GetCollisionAtNorth(unsigned int,
-								const std::vector<short>&,
-								unsigned int,
-								unsigned int);
-
-			bool GetCollisionAtNorthEast(unsigned int,
-								const std::vector<short>&,
-								unsigned int,
-								unsigned int);
-
-			bool GetCollisionAtEast(unsigned int,
-								const std::vector<short>&,
-								unsigned int,
-								unsigned int);
-
-			bool GetCollisionAtSouthEast(unsigned int,
-								const std::vector<short>&,
-								unsigned int,
-								unsigned int);
-
-			bool GetCollisionAtSouth(unsigned int,
-								const std::vector<short>&,
-								unsigned int,
-								unsigned int);
-
-			bool GetCollisionAtSouthWest(unsigned int,
-									const std::vector<short>&,
-									unsigned int,
-									unsigned int);
-
-			bool GetCollisionAtWest(unsigned int,
-								const std::vector<short>&,
-								unsigned int,
-								unsigned int);
-
-			bool GetCollisionAtNorthWest(unsigned int,
-									const std::vector<short>&,
-									unsigned int,
-									unsigned int);
-
-			bool IsInternalPixel(unsigned int cidx, std::vector<short>& bbox_array, 
-								unsigned int bbox_size_x, unsigned int bbox_size_y);
-	
-			unsigned int m_Rows;
-			unsigned int m_Cols;
-	};
-}
-
-#endif
diff --git a/src/Library/euc-dist-algorithm.h b/src/Library/euc-dist-algorithm.h
deleted file mode 100644
index ffe503fbf39b0e2675808b98fe4f876ee10c6d02..0000000000000000000000000000000000000000
--- a/src/Library/euc-dist-algorithm.h
+++ /dev/null
@@ -1,63 +0,0 @@
-/*=========================================================================
-
-  Program: Implementation of the euclidean distance criterion using the GRM library
-  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 __EUC_DIST_ALGORITHM_H
-#define __EUC_DIST_ALGORITHM_H
-
-#include "grm-region.h"
-#include "grm-algorithm.h"
-
-class EucDistRegion : public grm::RegionBase<EucDistRegion>
-{
-	public:
-		EucDistRegion(){}
-		~EucDistRegion(){}
-
-		float m_Area;
-		std::vector<float> m_Means;
-};
-
-struct EucDistParams
-{
-	float m_FinalThresh;
-	float m_InitialThresh;
-};
-
-template<class TInputImage>
-class EuclideanDistanceRM : public grm::RegionMergingAlgorithm<TInputImage, 
-														EucDistRegion,
-														EucDistParams>
-{
-	public:
-
-		// Mandatory typedef.
-		typedef grm::RegionMergingAlgorithm<TInputImage, EucDistRegion, EucDistParams> Superclass;
-		typedef EucDistRegion RegionType;
-		typedef std::shared_ptr<RegionType> RegionPointerType;
-		typedef std::vector<RegionPointerType> RegionListType;
-		typedef EucDistParams ParameterType;
-		typedef grm::Interface<RegionType> InterfaceType;
-		typedef typename Superclass::MergingCostFunction MergingCostFunction;
-
-		EuclideanDistanceRM(){}
-		virtual ~EuclideanDistanceRM(){}
-
-		void InitFromImage();
-		void UpdateAttribute(RegionPointerType, RegionPointerType);
-		virtual void Segmentation();
-};
-#include "euc-dist-algorithm.txx"
-#endif
diff --git a/src/Library/euc-dist-algorithm.txx b/src/Library/euc-dist-algorithm.txx
deleted file mode 100644
index 48d468907f9c5b0fa9b9f22b325d2775c7e240fc..0000000000000000000000000000000000000000
--- a/src/Library/euc-dist-algorithm.txx
+++ /dev/null
@@ -1,177 +0,0 @@
-/*=========================================================================
-
-  Program: Implementation of the euclidean distance criterion using the GRM library
-  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 __EUC_DIST_ALGORITHM_TXX
-#define __EUC_DIST_ALGORITHM_TXX
-
-#include "euc-dist-algorithm.h"
-#include <otbVectorImage.h>
-#include <otbImageFileReader.h>
-#include <itkImageRegionIterator.h>
-
-template<class TInputImage>
-void EuclideanDistanceRM<TInputImage>::InitFromImage()
-{
-	typedef otb::ImageFileReader<TInputImage> InputImageFileReaderType;
-	typedef itk::ImageRegionIterator<TInputImage> InputImageIteratorType;
-
-	// Read the input image
-	auto input_reader = InputImageFileReaderType::New();
-	input_reader->SetFileName(this->m_Input);
-	input_reader->Update();
-
-	// Retrieve useful information from the input image
-	auto regionToProcess = input_reader->GetOutput()->GetLargestPossibleRegion();
-	unsigned int rows = regionToProcess.GetSize()[1];
-	unsigned int cols = regionToProcess.GetSize()[0];
-	unsigned int bands = input_reader->GetOutput()->GetNumberOfComponentsPerPixel();
-
-	// Step 1: For efficiency the user is advised to allocate the initial number of regions
-	this->m_RegionList.reserve(rows*cols);
-
-	// Step 2: the user has to intialize the specific attributes of the initial regions
-	InputImageIteratorType it(input_reader->GetOutput(), regionToProcess);
-	for(it.GoToBegin(); !it.IsAtEnd(); ++it)
-	{
-		auto new_region = std::make_shared<EucDistRegion>();
-		new_region->m_Area = 1;
-		new_region->m_Means.reserve(bands);
-		for(unsigned int b = 0; b<bands; ++b)
-			new_region->m_Means.push_back(it.Get()[b]);
-		// The user has to push the initial regions
-		this->m_RegionList.push_back(new_region);
-	}
-
-	// Step 3: The neighborhhod and some internal attributes have to be initialized
-	// the user has to call the method InitRegions
-	this->InternalInit(rows, cols);
-}
-
-// Let assume that regions r1 and r2 have to be merged and r2 is merged into r1
-// The user has to describe how the specific attributes of r1 have to be updated
-// by defining the mandatory method UpdateAttribute(RegionPointerType, RegionPointerType)
-template<class TInputImage>
-void 
-EuclideanDistanceRM<TInputImage>::UpdateAttribute(RegionPointerType r1, RegionPointerType r2)
-{
-	float a1 = r1->m_Area, a2 = r2->m_Area, a_tot = a1 + a2;
-	for(unsigned int b=0; b<r1->m_Means.size(); ++b)
-		r1->m_Means[b] = ((a1*r1->m_Means[b]) + (a2*r2->m_Means[b])) / a_tot;
-	r1->m_Area += r2->m_Area;
-}
-
-// This method is obviously mandatory but its implementation depends on some
-// subtilities of your algorithm
-// What is common is the the iterative merging process of the region, what can be different
-// is the stopping criterion and the behavior of the threshold
-template<class TInputImage>
-void 
-EuclideanDistanceRM<TInputImage>::Segmentation()
-{
-	bool prev_merged = true;
-	float threshold;
-	unsigned int step = 0;
-
-	MergingCostFunction cost_func = [](RegionPointerType r1, RegionPointerType r2)->double
-	{
-		double cost = 0.0;
-		unsigned int bands = r1->m_Means.size();
-		for(unsigned int b = 0; b < bands; b++)
-			cost += pow((r1->m_Means[b] - r2->m_Means[b]), 2);
-		return cost;
-	};
-
-	while(this->m_Parameters.m_InitialThresh < this->m_Parameters.m_FinalThresh)
-	{
-		prev_merged = true;
-		threshold = static_cast<float>(this->m_Parameters.m_InitialThresh * 
-										this->m_Parameters.m_InitialThresh);
-		std::cout << "." << std::flush;
-		while(prev_merged && step < this->m_NumberOfIterations)
-		{
-			prev_merged = false;
-			++step;
-			this->m_RMHandler.UpdateMergingCosts(this->m_RegionList, cost_func);
-
-			// Do a loop on the regions of the graph
-			for(auto& r : this->m_RegionList)
-			{
-				// For each explored region, we determine if the LMBF is met
-				auto ref_region = this->m_RMHandler.CheckLMBF(r, threshold);
-				// If it holds then the pointer to ref_region is not null
-				if(ref_region != nullptr)
-				{
-					// Process to merge the regions
-					// Step 1: Merge the specific attributes
-					UpdateAttribute(ref_region, ref_region->GetClosestNeighbor());
-					// Step 2: Internal update (mandatory)
-					this->m_RMHandler.Update(ref_region, ref_region->GetClosestNeighbor());
-
-					prev_merged = true;
-				}
-			}
-
-			// After one iteration, you have to remove the expired regions, i.e. regions
-			// who have merged into a larger one
-			this->m_RMHandler.RemoveExpiredVertices(this->m_RegionList);
-		}
-
-		if(prev_merged && this->m_BestFitting)
-		{
-			while(prev_merged)
-			{
-				prev_merged = false;
-
-				this->m_RMHandler.UpdateMergingCosts(this->m_RegionList, cost_func);
-
-				for(auto& r : this->m_RegionList)
-				{
-					// For each explored region, we determine if the LMBF is met
-					auto ref_region = this->m_RMHandler.CheckBF(r, threshold);
-					// If it holds then the pointer to ref_region is not null
-					if(ref_region != nullptr)
-					{
-						if(ref_region == r->GetClosestNeighbor())
-						{
-							// Process to merge the regions
-							// Step 1: Merge the specific attributes
-							UpdateAttribute(ref_region, r);
-							// Step 2: Internal update (mandatory)
-							this->m_RMHandler.Update(ref_region, r);
-						}
-						else
-						{
-							// Process to merge the regions
-							// Step 1: Merge the specific attributes
-							UpdateAttribute(ref_region, ref_region->GetClosestNeighbor());
-							// Step 2: Internal update (mandatory)
-							this->m_RMHandler.Update(ref_region, ref_region->GetClosestNeighbor());
-						}
-						prev_merged = true;
-					}
-				}
-
-				// After one iteration, you have to remove the expired regions, i.e. regions
-				// who have merged into a larger one
-				this->m_RMHandler.RemoveExpiredVertices(this->m_RegionList);
-			}
-		}
-		++this->m_Parameters.m_InitialThresh;
-	}
-	std::cout << "\n";
-}
-
-#endif
diff --git a/src/Library/fls-algorithm.h b/src/Library/fls-algorithm.h
deleted file mode 100644
index 82ea98acb709cdbd5b8adb1510b23f23474fae90..0000000000000000000000000000000000000000
--- a/src/Library/fls-algorithm.h
+++ /dev/null
@@ -1,55 +0,0 @@
-/*=========================================================================
-
-  Program: Implementation of the Full lambda schedule criterion using the GRM library
-  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 __FLS_ALGORITHM_H
-#define __FLS_ALGORITHM_H
-
-#include "grm-region.h"
-#include "grm-algorithm.h"
-
-class FlsRegion : public grm::RegionBase<FlsRegion>
-{
-	public:
-		float m_Area;
-		std::vector<float> m_Means;
-};
-
-typedef float FLSParams;
-
-template<class TInputImage>
-class FLSAlgorithmRM : public grm::RegionMergingAlgorithm<TInputImage, 
-														FlsRegion,
-														FLSParams>
-{
-	public:
-		// Mandatory typedef
-		typedef grm::RegionMergingAlgorithm<TInputImage, FlsRegion, FLSParams> Superclass;
-		typedef FlsRegion RegionType;
-		typedef std::shared_ptr<RegionType> RegionPointerType;
-		typedef std::vector<RegionPointerType> RegionListType;
-		typedef FLSParams ParameterType;
-		typedef grm::Interface<RegionType> InterfaceType;
-		typedef typename Superclass::MergingCostFunction MergingCostFunction;
-
-		FLSAlgorithmRM(){}
-		virtual ~FLSAlgorithmRM(){}
-
-		void InitFromImage();
-		void UpdateAttribute(RegionPointerType, RegionPointerType);
-		virtual void Segmentation();
-};
-#include "fls-algorithm.txx"
-#endif
diff --git a/src/Library/fls-algorithm.txx b/src/Library/fls-algorithm.txx
deleted file mode 100644
index cf1cde58b9f2fee6766eaeda73bf514a3b2c0d84..0000000000000000000000000000000000000000
--- a/src/Library/fls-algorithm.txx
+++ /dev/null
@@ -1,167 +0,0 @@
-/*=========================================================================
-
-  Program: Implementation of the Full lambda schedule criterion using the GRM library
-  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 __FLS_ALGORITHM_TXX
-#define __FLS_ALGORITHM_TXX
-
-#include "fls-algorithm.h"
-#include <otbVectorImage.h>
-#include <otbImageFileReader.h>
-#include <itkImageRegionIterator.h>
-
-template<class TInputImage>
-void 
-FLSAlgorithmRM<TInputImage>::InitFromImage()
-{
-	typedef otb::ImageFileReader<TInputImage> InputImageFileReaderType;
-	typedef itk::ImageRegionIterator<TInputImage> InputImageIteratorType;
-
-	// Read the input image
-	auto input_reader = InputImageFileReaderType::New();
-	input_reader->SetFileName(this->m_Input);
-	input_reader->Update();
-
-	// Retrieve useful information from the input image
-	auto regionToProcess = input_reader->GetOutput()->GetLargestPossibleRegion();
-	unsigned int rows = regionToProcess.GetSize()[1];
-	unsigned int cols = regionToProcess.GetSize()[0];
-	unsigned int bands = input_reader->GetOutput()->GetNumberOfComponentsPerPixel();
-
-	// Step 1: For efficiency the user is advised to allocate the initial number of regions
-	this->m_RegionList.reserve(rows*cols);
-
-	// Step 2: the user has to intialize the specific attributes of the initial regions
-	InputImageIteratorType it(input_reader->GetOutput(), regionToProcess);
-	for(it.GoToBegin(); !it.IsAtEnd(); ++it)
-	{
-		auto new_region = std::make_shared<FlsRegion>();
-		new_region->m_Area = 1;
-		new_region->m_Means.reserve(bands);
-		for(unsigned int b = 0; b<bands; ++b)
-			new_region->m_Means.push_back(it.Get()[b]);
-		// The user has to push the initial regions
-		this->m_RegionList.push_back(new_region);
-	}
-
-	// Step 3: The neighborhhod and some internal attributes have to be initialized
-	// the user has to call the method InitRegions
-	this->InternalInit(rows, cols);
-}
-
-// Let assume that regions r1 and r2 have to be merged and r2 is merged into r1
-// The user has to describe how the specific attributes of r1 have to be updated
-// by defining the mandatory method UpdateAttribute(RegionPointerType, RegionPointerType)
-template<class TInputImage>
-void 
-FLSAlgorithmRM<TInputImage>::UpdateAttribute(RegionPointerType r1, RegionPointerType r2)
-{
-	float a1 = r1->m_Area, a2 = r2->m_Area, a_tot = a1 + a2;
-	for(unsigned int b=0; b<r1->m_Means.size(); ++b)
-		r1->m_Means[b] = ((a1*r1->m_Means[b]) + (a2*r2->m_Means[b])) / a_tot;
-	r1->m_Area += r2->m_Area;
-}
-
-// This method is obviously mandatory but its implementation depends on some
-// subtilities of your algorithm
-// What is common is the the iterative merging process of the region, what can be different
-// is the stopping criterion and the behavior of the threshold
-template<class TInputImage>
-void 
-FLSAlgorithmRM<TInputImage>::Segmentation()
-{
-	bool prev_merged = true;
-	unsigned int step = 0;
-	float threshold = this->m_Parameters*this->m_Parameters;
-
-	MergingCostFunction cost_func = [](RegionPointerType r1, RegionPointerType r2)->double
-	{
-		double cost = 0.0;
-		float a1 = r1->m_Area, a2 = r2->m_Area, a_tot = a1 + a2;
-		unsigned int c = r1->GetConnections(r2);
-		unsigned int bands = r1->m_Means.size();
-		for(unsigned int b = 0; b < bands; b++)
-			cost += pow((r1->m_Means[b] - r2->m_Means[b]), 2);
-		cost *= ((a1*a2) / (a_tot * c));
-		return cost;
-	};
-
-	while(prev_merged && step < this->m_NumberOfIterations)
-	{
-		std::cout << "." << std::flush;
-		prev_merged = false;
-		++step;
-
-		this->m_RMHandler.UpdateMergingCosts(this->m_RegionList, cost_func);
-
-		// Do a loop on the regions of the graph
-		for(auto& r : this->m_RegionList)
-		{
-			// For each explored region, we determine if the LMBF is met
-			auto ref_region = this->m_RMHandler.CheckLMBF(r, threshold);
-			// If it holds then the pointer to ref_region is not null
-			if(ref_region != nullptr)
-			{
-				// Process to merge the regions
-				// Step 1: Merge the specific attributes
-				UpdateAttribute(ref_region, ref_region->GetClosestNeighbor());
-				// Step 2: Internal update (mandatory)
-				this->m_RMHandler.Update(ref_region, ref_region->GetClosestNeighbor());
-
-				prev_merged = true;
-			}
-		}
-
-		// After one iteration, you have to remove the expired regions, i.e. regions
-		// who have merged into a larger one
-		this->m_RMHandler.RemoveExpiredVertices(this->m_RegionList);
-	}
-
-	if(prev_merged && this->m_BestFitting)
-	{
-		while(prev_merged)
-		{
-			std::cout << "." << std::flush;
-			prev_merged = false;
-
-			this->m_RMHandler.UpdateMergingCosts(this->m_RegionList, cost_func);
-
-			for(auto& r : this->m_RegionList)
-			{
-				// For each explored region, we determine if the LMBF is met
-				auto ref_region = this->m_RMHandler.CheckBF(r, threshold);
-				// If it holds then the pointer to ref_region is not null
-				if(ref_region != nullptr)
-				{
-					// Process to merge the regions
-					// Step 1: Merge the specific attributes
-					UpdateAttribute(ref_region, ref_region->GetClosestNeighbor());
-					// Step 2: Internal update (mandatory)
-					this->m_RMHandler.Update(ref_region, ref_region->GetClosestNeighbor());
-
-					prev_merged = true;
-				}
-			}
-
-			// After one iteration, you have to remove the expired regions, i.e. regions
-			// who have merged into a larger one
-			this->m_RMHandler.RemoveExpiredVertices(this->m_RegionList);
-		}
-	}
-
-	std::cout << "\n";
-}
-
-#endif
diff --git a/src/Library/grm-algorithm.h b/src/Library/grm-algorithm.h
deleted file mode 100644
index 3747a4d1070d2d02cc1876358e8184dd3402fb70..0000000000000000000000000000000000000000
--- a/src/Library/grm-algorithm.h
+++ /dev/null
@@ -1,87 +0,0 @@
-/*=========================================================================
-
-  Program: Generic Region Merging Library (GRM)
-  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 __GRM_ALGORITHM_H
-#define __GRM_ALGORITHM_H
-
-#include "macro-generator.h"
-#include "grm-interface.h"
-
-namespace grm
-{
-	template<class TInputImage, class TRegion, class TParams>
-	class RegionMergingAlgorithm
-	{
-		public:
-
-			// Mandatory typedef.
-			typedef TRegion RegionType;
-			typedef std::shared_ptr<RegionType> RegionPointerType;
-			typedef std::vector<RegionPointerType> RegionListType;
-			typedef TParams ParameterType;
-			typedef Interface<RegionType> InterfaceType;
-			typedef typename InterfaceType::MergingCostFunction MergingCostFunction;
-
-			// Empty constructor and destructor
-			RegionMergingAlgorithm() : m_BestFitting{true}, m_NumberOfIterations{70}
-			{}
-			virtual ~RegionMergingAlgorithm(){}
-
-			// Abstract method
-			virtual void InitFromImage() = 0;
-			virtual void UpdateAttribute(RegionPointerType, RegionPointerType) = 0;
-			virtual void Segmentation() = 0;
-			virtual void Run();
-
-			// No abstract method
-			virtual void WriteLabelOutputImage();
-			virtual void WriteRGBOutputImage();
-
-			// Get / Set methods
-			GetMacro(RegionListType, RegionList);
-			GetMacro(ParameterType, Parameters);
-			GetMacro(unsigned int, NumberOfIterations);
-
-			SetMacro(RegionListType, RegionList);
-			SetMacro(ParameterType, Parameters);
-			SetMacro(std::string, Input);
-			SetMacro(std::string, OutputLabel);
-			SetMacro(std::string, OutputRGB);
-			SetMacro(unsigned int, NumberOfIterations);
-			void SetDimensionForEncoder(unsigned int, unsigned int);
-			void SetBestFitting(int f);
-
-		protected:
-
-			virtual void InternalInit(unsigned int, unsigned int);
-
-			// Input / Output attributes
-			std::string m_Input;
-			std::string m_OutputLabel;
-			std::string m_OutputRGB;
-
-			// Algorithm attributes
-			bool m_BestFitting;
-			RegionListType m_RegionList;
-			ParameterType m_Parameters;
-
-			// Algorithm helper
-			unsigned int m_NumberOfIterations;
-			InterfaceType m_RMHandler;
-	};
-} // end of namespace grm
-#include "grm-algorithm.txx"
-#endif
diff --git a/src/Library/grm-algorithm.txx b/src/Library/grm-algorithm.txx
deleted file mode 100644
index b46b3f5f3cb14d0f378d250a9afe8c1bb226456d..0000000000000000000000000000000000000000
--- a/src/Library/grm-algorithm.txx
+++ /dev/null
@@ -1,158 +0,0 @@
-/*=========================================================================
-
-  Program: Generic Region Merging Library (GRM)
-  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 __GRM_ALGORITHM_TXX
-#define __GRM_ALGORITHM_TXX
-
-#include "otbImage.h"
-#include "otbVectorImage.h"
-#include "otbImageFileReader.h"
-#include "otbImageFileWriter.h"
-#include "grm-algorithm.h"
-
-namespace grm
-{
-	template<class TInputImage, class TRegion, class TParams>
-	void
-	RegionMergingAlgorithm<TInputImage, TRegion, TParams>::SetBestFitting(int f)
-	{
-		assert(f == 0 || f == 1);
-		if(f)
-			m_BestFitting = true;
-		else
-			m_BestFitting = false;
-	}
-
-	template<class TInputImage, class TRegion, class TParams>
-	void
-	RegionMergingAlgorithm<TInputImage, TRegion, TParams>::SetDimensionForEncoder(
-															unsigned int rows,
-															unsigned int columns)
-	{
-		m_RMHandler.SetDimensionForShapeEncoder(rows, columns);
-	}
-
-	template<class TInputImage, class TRegion, class TParams>
-	void 
-	RegionMergingAlgorithm<TInputImage, TRegion, TParams>::InternalInit(unsigned int rows, 
-																		unsigned int columns)
-	{
-		m_RMHandler.SetDimensionForShapeEncoder(rows, columns);
-		m_RMHandler.InitRegions(m_RegionList);
-	}
-
-	template<class TInputImage, class TRegion, class TParams>
-	void
-	RegionMergingAlgorithm<TInputImage, TRegion, TParams>::Run()
-	{
-		InitFromImage();
-		Segmentation();
-		WriteRGBOutputImage();
-		WriteLabelOutputImage();
-	}
-
-	template<class TInputImage, class TRegion, class TParams>
-	void
-	RegionMergingAlgorithm<TInputImage, TRegion, TParams>::WriteLabelOutputImage()
-	{
-		assert(m_Input != "");
-		assert(m_OutputLabel != "");
-
-		typedef long unsigned int LabelPixelType;
-		typedef otb::Image<LabelPixelType, 2> LabelImageType;
-		typedef typename LabelImageType::IndexType LabelIndexType;
-		typedef typename LabelImageType::SizeType LabelSizeType;
-		typedef typename LabelImageType::RegionType LabelRegionType;
-
-		typedef otb::ImageFileWriter<LabelImageType> LabelImageWriterType;
-
-		// Allocate the image
-		auto label_img = LabelImageType::New();
-		LabelIndexType index; index[0] = 0; index[1] = 0;
-		LabelSizeType size;
-		auto dim = m_RMHandler.GetDimensionOfInputImage();
-		size[0] = dim.first;
-		size[1] = dim.second;
-		assert(size[0] > 0);
-		assert(size[1] > 0);
-		LabelRegionType region; region.SetIndex(index); region.SetSize(size);
-		label_img->SetRegions(region);
-		label_img->Allocate();
-
-		for(unsigned int x = 0; x < size[0]; ++x)
-		{
-			for(unsigned int y = 0; y < size[1]; ++y)
-			{
-				index[0] = x;
-				index[1] = y;
-				label_img->SetPixel(index, 0);
-			}
-		}
-
-		long unsigned int l = 1; // Lower value of the label
-		std::for_each(m_RegionList.begin(), m_RegionList.end(), [&](RegionPointerType r){
-			auto pixels = m_RMHandler.GenerateAllPixels(r);
-			std::for_each(pixels.begin(), pixels.end(), [&](long unsigned int p){
-				index[0] = p % size[0];
-				index[1] = p / size[0];
-				label_img->SetPixel(index, l);
-			});
-			++l;
-		});
-
-		auto label_writer = LabelImageWriterType::New();
-		label_writer->SetFileName(m_OutputLabel);
-		label_writer->SetInput(label_img);
-		label_writer->Update();
-	}
-
-	template<class TInputImage, class TRegion, class TParams>
-	void
-	RegionMergingAlgorithm<TInputImage, TRegion, TParams>::WriteRGBOutputImage()
-	{
-		typedef unsigned char RGBPixelType;
-		typedef otb::VectorImage<RGBPixelType, 2> RGBImageType;
-		typedef typename RGBImageType::IndexType RGBIndexType;
-		
-		typedef otb::ImageFileReader<RGBImageType> RGBImageReaderType;
-		typedef otb::ImageFileWriter<RGBImageType> RGBImageWriterType;
-
-		auto rgb_reader = RGBImageReaderType::New();
-		rgb_reader->SetFileName(m_Input);
-		rgb_reader->Update();
-		auto rgb_img = rgb_reader->GetOutput();
-		RGBIndexType index;
-	
-		auto dim = m_RMHandler.GetDimensionOfInputImage();
-		std::for_each(m_RegionList.begin(), m_RegionList.end(), [&](RegionPointerType r){
-			auto pixels = m_RMHandler.GenerateFrontierPixels(r);
-			std::for_each(pixels.begin(), pixels.end(), [&](long unsigned int p){
-				index[0] = p % dim.first;
-				index[1] = p / dim.first;
-				auto rgb_pix = rgb_img->GetPixel(index);
-				for(int b=0; b<rgb_img->GetNumberOfComponentsPerPixel(); b++)
-					rgb_pix[b] = 0;
-				rgb_img->SetPixel(index, rgb_pix);
-			});
-		});
-
-		auto rgb_writer = RGBImageWriterType::New();
-		rgb_writer->SetFileName(m_OutputRGB);
-		rgb_writer->SetInput(rgb_img);
-		rgb_writer->Update();
-	}
-} // end of namespace grm
-#endif
diff --git a/src/Library/grm-interface.h b/src/Library/grm-interface.h
deleted file mode 100644
index 902e6c7e08084299725ef7610b948e71b6b617e1..0000000000000000000000000000000000000000
--- a/src/Library/grm-interface.h
+++ /dev/null
@@ -1,61 +0,0 @@
-/*=========================================================================
-
-  Program: Generic Region Merging Library (GRM)
-  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 __GRM_INTERFACE_H
-#define __GRM_INTERFACE_H
-
-#include "contour-encoder.h"
-
-#include <memory>
-#include <tuple>
-#include <functional>
-#include <limits>
-
-namespace grm
-{
-	template<class TRegion>
-	class Interface
-	{
-		public:
-
-			typedef std::shared_ptr<TRegion> RegionPointerType;
-			typedef std::tuple<RegionPointerType, double, unsigned int> EdgeType;
-			typedef std::vector<RegionPointerType> RegionListType;
-			typedef std::function<double(RegionPointerType, RegionPointerType)> MergingCostFunction;
-
-			void InitRegions(RegionListType&);
-			void UpdateMergingCosts(RegionListType&, MergingCostFunction&);
-			RegionPointerType CheckBF(RegionPointerType r, double thresh);
-			RegionPointerType CheckLMBF(RegionPointerType, double);
-			void UpdateContour(RegionPointerType, RegionPointerType);
-			void UpdateNeighborhood(RegionPointerType, RegionPointerType);
-			void Update(RegionPointerType r1, RegionPointerType r2);
-			void RemoveExpiredVertices(RegionListType& graph);
-
-			std::vector<long unsigned int> GenerateFrontierPixels(RegionPointerType);
-			std::vector<long unsigned int> GenerateAllPixels(RegionPointerType);
-
-			// Get / Set methods
-			std::pair<unsigned int, unsigned int> GetDimensionOfInputImage();
-			void SetDimensionForShapeEncoder(unsigned int, unsigned int);
-
-		private:
-			shpenc::ShapeEncoder m_ShpEncoder;
-	};
-} // end of namespace grm
-
-#include "grm-interface.txx"
-#endif
diff --git a/src/Library/grm-interface.txx b/src/Library/grm-interface.txx
deleted file mode 100644
index e5d2c3290e1a4102aac8d53aab5b13f37be308e3..0000000000000000000000000000000000000000
--- a/src/Library/grm-interface.txx
+++ /dev/null
@@ -1,249 +0,0 @@
-/*=========================================================================
-
-  Program: Generic Region Merging Library (GRM)
-  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.
-
-=========================================================================*/
-#include "grm-interface.h"
-
-namespace grm
-{
-
-	template<class TRegion>
-	void 
-	Interface<TRegion>::SetDimensionForShapeEncoder(unsigned int rows, unsigned int cols)
-	{
-		assert(rows > 0 && cols > 0);
-		m_ShpEncoder.SetRows(rows);
-		m_ShpEncoder.SetCols(cols);
-	}
-
-	template<class TRegion>
-	std::pair<unsigned int, unsigned int> 
-	Interface<TRegion>::GetDimensionOfInputImage()
-	{
-		return std::pair<unsigned int, unsigned int>(m_ShpEncoder.GetCols(), m_ShpEncoder.GetRows());
-	}
-
-	template<class TRegion>
-	std::vector<long unsigned int> 
-	Interface<TRegion>::GenerateFrontierPixels(RegionPointerType r)
-	{
-		return m_ShpEncoder.GeneratePixels(r->GetId(), r->GetContour());
-	}
-
-	template<class TRegion>
-	std::vector<long unsigned int> 
-	Interface<TRegion>::GenerateAllPixels(RegionPointerType r)
-	{
-		return m_ShpEncoder.GenerateAllPixels(r->GetId(), r->GetContour(), r->GetBbox());
-	}
-
-	template<class TRegion>
-	void
-	Interface<TRegion>::InitRegions(RegionListType& graph)
-	{
-		assert(graph.size() > 0);
-		long unsigned int id_ = 0;
-		std::for_each(graph.begin(), graph.end(), [&](RegionPointerType r){
-			r->SetId(id_);
-			r->SetBboxIndex(0, id_ % m_ShpEncoder.GetCols());
-			r->SetBboxIndex(1, id_ / m_ShpEncoder.GetCols());
-			r->SetBboxSize(0,1);
-			r->SetBboxSize(1,1);
-			r->SetValid(true);
-			r->SetExpired(false);
-			r->AddContour(1); r->AddContour(2);
-			r->AddContour(3); r->AddContour(0);
-
-			auto neighbors = mask::CrossNeighborhood(r->GetId(), m_ShpEncoder.GetCols(), 
-													m_ShpEncoder.GetRows());
-			for(short j = 0; j < 4; ++j)
-			{
-				if(neighbors[j] > -1)
-					r->AddNeighbor(graph[neighbors[j]], 0, 1);
-			}
-			++id_;
-		});
-	}
-
-	template<class TRegion>
-	void
-	Interface<TRegion>::UpdateMergingCosts(RegionListType& graph,
-											MergingCostFunction& func)
-	{
-		double min_cost;
-		long unsigned int min_id;
-		int idx_tpl, idx_min_tpl;
-
-		std::for_each(graph.begin(), graph.end(), [&](RegionPointerType r){
-
-			r->SetValid(true);
-			r->SetExpired(false);
-
-			min_cost = std::numeric_limits<double>::max();
-			idx_tpl = 0;
-			idx_min_tpl = 0;
-			std::for_each(r->EdgeListBegin(), r->EdgeListEnd(), [&](EdgeType& e){
-				std::get<1>(e) = func(r, std::get<0>(e));
-				if(min_cost > std::get<1>(e))
-				{
-					min_cost = std::get<1>(e);
-					min_id = std::get<0>(e)->GetId();
-					idx_min_tpl = idx_tpl;
-					++idx_tpl;
-				}
-				else if(min_cost == std::get<1>(e) && min_id > std::get<0>(e)->GetId())
-				{
-					min_id = std::get<0>(e)->GetId();
-					idx_min_tpl = idx_tpl;
-					++idx_tpl;
-				}
-				else
-					++idx_tpl;
-			});
-
-			if(idx_min_tpl > 0)
-				r->SwapEdges(0, idx_min_tpl);
-		});
-	}
-
-	template<class TRegion>
-	typename Interface<TRegion>::RegionPointerType
-	Interface<TRegion>::CheckLMBF(RegionPointerType r, double thresh)
-	{
-		assert(r != nullptr);
-		if(r->GetValid())
-		{
-			auto cost = std::get<1>(r->FrontEdgeList());
-			if(cost < thresh)
-			{
-				auto closest_neigh = std::get<0>(r->FrontEdgeList());
-				if(closest_neigh->GetValid())
-				{
-					if(std::get<0>(closest_neigh->FrontEdgeList()) == r)
-					{
-						if(r->GetId() < closest_neigh->GetId())
-							return r;
-						else
-							return closest_neigh;
-					}
-					else
-						return nullptr;
-				}
-				else
-					return nullptr;
-			}
-			else
-				return nullptr;
-		}
-		else
-			return nullptr;
-	}
-
-	template<class TRegion>
-	typename Interface<TRegion>::RegionPointerType
-	Interface<TRegion>::CheckBF(RegionPointerType r, double thresh)
-	{
-		assert(r != nullptr);
-		if(r->GetValid())
-		{
-			auto cost = std::get<1>(r->FrontEdgeList());
-			if(cost < thresh)
-			{
-				auto closest_neigh = std::get<0>(r->FrontEdgeList());
-				if(closest_neigh->GetValid())
-				{
-					if(r->GetId() < closest_neigh->GetId())
-						return r;
-					else
-						return closest_neigh;
-				}
-				else
-					return nullptr;
-			}
-			else
-				return nullptr;
-		}
-		else
-			return nullptr;
-	}
-
-	template<class TRegion>
-	void 
-	Interface<TRegion>::UpdateContour(RegionPointerType r1, RegionPointerType r2)
-	{
-		r1->UpdateBbox(r2);
-		r1->SetContour(m_ShpEncoder.EncodeContour(r1->GetId(), r1->GetContour(), 
-													r2->GetId(), r2->GetContour(),
-													r1->GetBbox()));
-	}
-
-	template<class TRegion>
-	void
-	Interface<TRegion>::UpdateNeighborhood(RegionPointerType r1, RegionPointerType r2)
-	{
-		unsigned int connections;
-		for(auto adj_n = r2->EdgeListBegin(); adj_n != r2->EdgeListEnd(); ++adj_n)
-		{
-			auto adj_nr = std::get<0>(*adj_n);
-			auto toR2 = adj_nr->FindEdge(r2);
-
-			if(toR2 == adj_nr->EdgeListEnd())
-				adj_nr->SetValid(false);
-
-			connections = std::get<2>(*toR2);
-			adj_nr->Erase(toR2);
-
-			if(adj_nr != r1)
-			{
-				auto toR1 = adj_nr->FindEdge(r1);
-
-				if(toR1 == adj_nr->EdgeListEnd())
-				{
-					adj_nr->AddNeighbor(r1, 0, connections);
-					r1->AddNeighbor(adj_nr, 0, connections);
-				}
-				else
-				{
-					// We increment the connections between adj_nr and r1 by connections
-					std::get<2>(*toR1) += connections;
-					auto toAdjNr = r1->FindEdge(adj_nr);
-					std::get<2>(*toAdjNr) += connections;
-				}
-			}
-		}
-	}
-
-	template<class TRegion>
-	void
-	Interface<TRegion>::Update(RegionPointerType r1, RegionPointerType r2)
-	{
-		assert(r1 != nullptr && r2 != nullptr);
-		UpdateContour(r1,r2);
-		UpdateNeighborhood(r1, r2);
-		r1->SetValid(false);
-		r2->SetValid(false);
-		r2->SetExpired(true);
-	}
-
-	template<class TRegion>
-	void
-	Interface<TRegion>::RemoveExpiredVertices(RegionListType& graph)
-	{
-		auto exp = std::remove_if(graph.begin(), graph.end(), [](RegionPointerType r)->bool{
-			return r->GetExpired();
-		});
-		graph.erase(exp, graph.end());
-	}
-}
diff --git a/src/Library/grm-region.h b/src/Library/grm-region.h
deleted file mode 100644
index ba48d841eaffdeac9d2db085427b280486fb745a..0000000000000000000000000000000000000000
--- a/src/Library/grm-region.h
+++ /dev/null
@@ -1,107 +0,0 @@
-/*=========================================================================
-
-  Program: Generic Region Merging Library (GRM)
-  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 __GRM_REGION_H
-#define __GRM_REGION_H
-
-#include "macro-generator.h"
-
-#include <cassert>
-#include <memory>
-#include <array>
-#include <vector>
-#include <tuple>
-#include <bitset>
-#include <utility>
-#include <algorithm>
-
-namespace grm
-{
-	class Region
-	{
-		public:
-
-			Region(){}
-			virtual ~Region(){}
-
-			// Get / Set methods
-			GetMacro(bool, Valid);
-			GetMacro(bool, Expired);
-			GetMacro(long unsigned int, Id);
-			GetRefMacro(std::vector<std::bitset<2>>, Contour);
-			std::array<unsigned int,4>& GetBbox(){return m_Bbox;}
-			unsigned int GetBboxIndex(unsigned int);
-			unsigned int GetBboxSize(unsigned int);
-
-			SetMacro(bool, Valid);
-			SetMacro(bool, Expired);
-			SetMacro(long unsigned int, Id);
-			SetMacro(std::vector<std::bitset<2>>, Contour);
-			void SetBboxIndex(unsigned int, unsigned int);
-			void SetBboxSize(unsigned int, unsigned int);
-
-			void AddContour(short);
-			std::array<unsigned int, 4> GetResultingBbox(std::shared_ptr<Region>);
-			void UpdateBbox(std::shared_ptr<Region>);
-			void UpdateContour(std::shared_ptr<Region>);
-
-		protected:
-
-			bool m_Valid; // True if the region has to be considered as the closest neighbor candidate
-			bool m_Expired; // True if this region has to be removed from the list
-			long unsigned int m_Id; // Identification of the region (represents also the coordinate of the start pixel)
-			std::array<unsigned int, 4> m_Bbox;	// Rectangular bounding box of the region:
-												// Idx 0 : upper left x coordinate
-												// Idx 1 : upper left y coordinate
-												// Idx 2 : size along x axis (width)
-												// Idx 3 : size along y axis (height)
-			std::vector<std::bitset<2>> m_Contour; // Freeman representation of the contour of a region
-	};
-
-	template<class DerivedRegion>
-	class RegionBase : public Region
-	{
-		public:
-
-			typedef std::tuple<
-				std::shared_ptr<DerivedRegion>, // Pointer to the neighboring region
-				double, // merging cost
-				unsigned int // Length of the common boundary with this neighboring region
-			> EdgeType;
-			typedef std::vector<EdgeType> EdgeListType;
-			typedef typename EdgeListType::iterator EdgeIteratorType;
-
-			RegionBase(){}
-			virtual ~RegionBase(){}
-
-			EdgeIteratorType EdgeListBegin();
-			EdgeIteratorType EdgeListEnd();
-			EdgeIteratorType FindEdge(std::shared_ptr<DerivedRegion>);
-			EdgeIteratorType Erase(EdgeIteratorType);
-			const EdgeType& FrontEdgeList();
-			std::shared_ptr<DerivedRegion> GetClosestNeighbor();
-			void AddNeighbor(std::shared_ptr<DerivedRegion>, double, unsigned int);
-			void SwapEdges(int idx1, int idx2);
-			void UpdateNeighborhood(std::shared_ptr<DerivedRegion>);
-			unsigned int GetConnections(std::shared_ptr<DerivedRegion>);
-
-		protected:
-			EdgeListType m_Neighbors;
-	};
-} // end of namespace grm
-
-#include "grm-region.txx"
-#endif
diff --git a/src/Library/grm-region.txx b/src/Library/grm-region.txx
deleted file mode 100644
index 88a24104091f7a7d908d3a42985394f3229c8076..0000000000000000000000000000000000000000
--- a/src/Library/grm-region.txx
+++ /dev/null
@@ -1,166 +0,0 @@
-/*=========================================================================
-
-  Program: Generic Region Merging Library (GRM)
-  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.
-
-=========================================================================*/
-#include "grm-region.h"
-
-namespace grm
-{
-	unsigned int
-	Region::GetBboxIndex(unsigned int idx)
-	{
-		assert(idx >= 0 && idx < 2);
-		return m_Bbox[idx];
-	}
-
-	unsigned int
-	Region::GetBboxSize(unsigned int idx)
-	{
-		assert(idx >= 0 && idx < 2);
-		return m_Bbox[2+idx];
-	}
-
-	void
-	Region::SetBboxIndex(unsigned int idx, unsigned int v)
-	{
-		assert(idx >= 0 && idx < 2);
-		m_Bbox[idx] = v;
-	}
-
-	void
-	Region::SetBboxSize(unsigned int idx, unsigned int v)
-	{
-		assert(idx >= 0 && idx < 2);
-		m_Bbox[idx+2] = v;
-	}
-
-	void
-	Region::AddContour(short move)
-	{
-		assert(move >=0 && move < 4);
-		m_Contour.push_back(move);
-	}
-
-	std::array<unsigned int, 4>
-	Region::GetResultingBbox(std::shared_ptr<Region> r)
-	{
-		unsigned int min_row[3], max_row[3], min_col[3], max_col[3];
-
-		min_row[0] = m_Bbox[1];
-		max_row[0] = m_Bbox[1] + m_Bbox[3];
-		min_col[0] = m_Bbox[0];
-		max_col[0] = m_Bbox[0] + m_Bbox[2];
-
-		min_row[1] = r->m_Bbox[1];
-		max_row[1] = r->m_Bbox[1] + r->m_Bbox[3];
-		min_col[1] = r->m_Bbox[0];
-		max_col[1] = r->m_Bbox[0] + r->m_Bbox[2];
-
-		if (min_row[0]<min_row[1]) min_row[2]=min_row[0]; else min_row[2]=min_row[1];
-		if (min_col[0]<min_col[1]) min_col[2]=min_col[0]; else min_col[2]=min_col[1];
-		if (max_row[0]>max_row[1]) max_row[2]=max_row[0]; else max_row[2]=max_row[1];
-		if (max_col[0]>max_col[1]) max_col[2]=max_col[0]; else max_col[2]=max_col[1];
-
-		std::array<unsigned int,4> res;
-		res[1] = min_row[2];
-		res[0] = min_col[2];
-		res[3] = max_row[2] - min_row[2];
-		res[2] = max_col[2] - min_col[2];
-		return res;
-	}
-
-	void 
-	Region::UpdateBbox(std::shared_ptr<Region> r)
-	{
-		m_Bbox = GetResultingBbox(r);
-	}
-
-	template<class DerivedRegion>
-	void
-	RegionBase<DerivedRegion>::AddNeighbor(std::shared_ptr<DerivedRegion> r,
-											double cost,
-											unsigned int len_bound)
-	{
-		assert(len_bound > 0);
-		m_Neighbors.push_back(std::make_tuple(r, cost, len_bound));
-	}
-
-	template<class DerivedRegion>
-	typename RegionBase<DerivedRegion>::EdgeIteratorType 
-	RegionBase<DerivedRegion>::EdgeListBegin()
-	{
-		return m_Neighbors.begin();
-	}
-
-	template<class DerivedRegion>
-	const typename RegionBase<DerivedRegion>::EdgeType& 
-	RegionBase<DerivedRegion>::FrontEdgeList()
-	{
-		assert(m_Neighbors.size() >= 1);
-		return m_Neighbors.front();
-	}
-
-	template<class DerivedRegion>
-	typename RegionBase<DerivedRegion>::EdgeIteratorType 
-	RegionBase<DerivedRegion>::EdgeListEnd()
-	{
-		return m_Neighbors.end();
-	}
-
-	template<class DerivedRegion>
-	typename RegionBase<DerivedRegion>::EdgeIteratorType 
-	RegionBase<DerivedRegion>::Erase(EdgeIteratorType eit)
-	{
-		return m_Neighbors.erase(eit);
-	}
-
-	template<class DerivedRegion>
-	void 
-	RegionBase<DerivedRegion>::SwapEdges(int idx1, int idx2)
-	{
-		assert(idx1 == 0);
-		assert(idx2 > idx1 && idx2 < m_Neighbors.size());
-		std::swap(m_Neighbors[idx1], m_Neighbors[idx2]);
-	}
-
-	template<class DerivedRegion>
-	typename RegionBase<DerivedRegion>::EdgeIteratorType 
-	RegionBase<DerivedRegion>::FindEdge(std::shared_ptr<DerivedRegion> r)
-	{
-		assert(r != nullptr);
-		return std::find_if(m_Neighbors.begin(), m_Neighbors.end(), 
-					[&](const EdgeType& n)->bool{
-						return std::get<0>(n) == r;
-				});
-	}
-
-	template<class DerivedRegion>
-	unsigned int 
-	RegionBase<DerivedRegion>::GetConnections(std::shared_ptr<DerivedRegion> r)
-	{
-		assert(r != nullptr);
-		auto toR = FindEdge(r);
-		assert(toR != EdgeListEnd());
-		return std::get<2>(*toR);
-	}
-
-	template<class DerivedRegion>
-	std::shared_ptr<DerivedRegion> 
-	RegionBase<DerivedRegion>::GetClosestNeighbor()
-	{
-		assert(m_Neighbors.size() > 0);
-		return std::get<0>(FrontEdgeList());
-	}
-} // end of namespace grm
diff --git a/src/Library/mask-neighborhood.cxx b/src/Library/mask-neighborhood.cxx
deleted file mode 100644
index 7b1b46bd71cb75156595eb55ee94141670dc6642..0000000000000000000000000000000000000000
--- a/src/Library/mask-neighborhood.cxx
+++ /dev/null
@@ -1,89 +0,0 @@
-/*=========================================================================
-
-  Program: Computation of the neighborhood
-  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.
-
-=========================================================================*/
-#include "mask-neighborhood.h"
-
-namespace mask
-{
-	std::array<long int, 4> CrossNeighborhood(long unsigned int coord, unsigned int sizex, unsigned int sizey)
-	{
-		std::array<long int, 4> neighbors;
-		neighbors.fill(-1);
-
-		unsigned int x = coord % sizex;
-		unsigned int y = coord / sizex;
-
-		//top
-		if(y > 0)
-			neighbors[0] = coord - sizex;
-
-		// right
-		if(x < sizex - 1)
-			neighbors[1] = coord + 1;
-
-		// bottom
-		if(y < sizey - 1)
-			neighbors[2] = coord + sizex;
-
-		// left
-		if(x > 0)
-			neighbors[3] = coord - 1;
-
-		return neighbors;
-	}
-
-	std::array<long int, 8> SquareNeighborhood(long unsigned int coord, unsigned int sizex, unsigned int sizey)
-	{
-		std::array<long int, 8> neighbors;
-		neighbors.fill(-1);
-
-		unsigned int x = coord % sizex;
-		unsigned int y = coord / sizex;
-
-		//top
-		if(y > 0)
-			neighbors[0] = coord - sizex;
-
-		// top right
-		if(y > 0 && x < sizex - 1)
-			neighbors[1] = coord - sizex + 1;
-
-		// right
-		if(x < sizex - 1)
-			neighbors[2] = coord + 1;
-
-		// bottom right
-		if(y > sizey - 1 && x < sizex - 1)
-			neighbors[3] = coord + sizex + 1;
-
-		// bottom
-		if(y < sizey - 1)
-			neighbors[4] = coord + sizex;
-
-		// bottom left
-		if(y > sizey - 1 && x > 0)
-			neighbors[5] = coord + sizex - 1;
-
-		// left
-		if(x > 0)
-			neighbors[6] = coord - 1;
-
-		// top left
-		if(y > 0 && x > 0)
-			neighbors[7] = coord - 1 - sizex;
-		return neighbors;
-	}
-} // end of namespace lss
diff --git a/src/Library/mask-neighborhood.h b/src/Library/mask-neighborhood.h
deleted file mode 100644
index 85f9173070e4e865ec474188d0ed8f15ffb368fc..0000000000000000000000000000000000000000
--- a/src/Library/mask-neighborhood.h
+++ /dev/null
@@ -1,39 +0,0 @@
-/*=========================================================================
-
-  Program: Computation of the neighborhood
-  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 __MASKNEIGHBORHOOD_H
-#define __MASKNEIGHBORHOOD_H
-
-#include <array>
-
-namespace mask
-{
-// Different masks for the neighborhood
-
-//     *
-//    *X*
-//     *
-
-std::array<long int, 4> CrossNeighborhood(long unsigned int coord, unsigned int sizex, unsigned int sizey);
-
-//    ***
-//    *X*
-//    ***
-
-std::array<long int, 8> SquareNeighborhood(long unsigned int coord, unsigned int sizex, unsigned int sizey);
-}// end of namespace lss
-
-#endif