diff --git a/doc/diagram.classes b/doc/diagram.classes
index 54ac0ed..2285a52 100644
--- a/doc/diagram.classes
+++ b/doc/diagram.classes
@@ -1,113 +1,125 @@
 ################################################################################
 # Data types
 ################################################################################
 class Array
   public Array();
   public Array(UInt size);
   public Array(T * data, UInt size);
   public ~Array();
   public Array & operator=(Array & v);
   public wrapMemory(T * data, UInt size);
   public assign(UInt size, const T & value);
   public T * getPtr();
   public setPtr(T * new_ptr);
   public resize(UInt size);
   public T & operator[]();
   public UInt size();
 
 class Grid
   public Grid();
   public Grid(int n, double L, int components);
   public ~Grid();
   public iterator begin(UInt n);
   public iterator end(UInt n);
   public resize();
   public UInt computeSize();
   public printself(std::ostream & str);
   public UInt * sizes();
   public UInt dataSize();
   public UInt getNbPoints();
   public Real * lengths();
   public UInt getNbComponents();
   public setNbComponents(UInt n);
   public T * getInternalData();
   public uniformSetComponents(Grid & vec);
   private UInt unpackOffset(UInt offset, UInt index);
   private UInt unpackOffset(UInt offset, UInt index, T1... rest);
   private UInt computeOffset(UInt[dim+1] tuple);
   private Array data;
   private UInt nb_components;
   private UInt[dim] n;
   private Real[dim] L;
 
 class GridHermitian(Grid)
 
 class GridView(Grid)
 
 ################################################################################
 # FFT classes
 ################################################################################
 class FFTransform
   public FFTransform(Grid & real, GridHermitian & spectral);
   public ~FFTransform();
   public forwardTransform();
   public backwardTransform();
   public normalize();
   private Grid & real;
   private Grid & spectral;
 
 class FFTransformFFTW(FFTransform)
   private fftw_plan forward_plan;
   private fftw_plan backward_plan;
 
 class FFTransformCuFFT(FFTransform)
 
 class FFTransformManager
   public createPlan(Grid input, GridHermitian output);
   public clean();
   public destroyPlan(Grid input, GridHermitian output);
   public destroyPlan(T * some);
   private std::map transforms;
   private static FFTPlanManager * singleton;
 
 ################################################################################
 # Surface generation
 ################################################################################
 class SurfaceGenerator
 
 class Filter
   public functor();
   private std::map parameters;
 
 class Isopowerlaw(Filter)
   public setHurst(double hurst);
   public setWavenumbers(std::map wavenumbers);
   public setWavelengths(std::map wavelengths);
 
 class PythonFilter(Filter)
 
 class SurfaceGeneratorFilter(SurfaceGenerator)
   private Filter filter;
   public Grid generateSurface();
 
 ################################################################################
 # Mechanics
 ################################################################################
 class Model
   private Grid displacements;
   private Grid surface;
   private Grid pressure;
   private BEM bem;
 
 class Solver
   private Model model;
   public solveEquilibrium();
 
 class Polonsky(Solver)
 class Kato(Solver)
 class Gigi(Solver)
 class Condat(Solver)
 
 class BEM
 
-class Westergaard(BEM)
\ No newline at end of file
+class Westergaard(BEM)
+
+################################################################################
+# Flood fill
+################################################################################
+class Cluster
+  public Cluster(std::pair start, Grid contact, Grid visited, bool diagonal);
+  public getArea();
+  private std::list points;
+  private UInt area;
+
+class FloodFill
+  public static std::list getClusters(Grid contact, bool diagonal);
diff --git a/python/clusters.i b/python/clusters.i
index 5a6e364..eb4f209 100644
--- a/python/clusters.i
+++ b/python/clusters.i
@@ -1,14 +1,17 @@
 /* -------------------------------------------------------------------------- */
 /* Wrapping contact clusters classes                                           */
 /* -------------------------------------------------------------------------- */
 %ignore operator <<;
 
 %include "contact_area.hh"
 %include "contact_cluster_collection.hh"
+%include "flood_fill.hh"
 
 %include "std_vector.i"
+%include "std_list.i"
 namespace tamaas {
 %template(VectorClusters) ::std::vector<::tamaas::ContactCluster>;
+%template(ListClusters) ::std::list<::tamaas::Cluster>;
 }
 %include "contact_cluster.hh"
 
diff --git a/python/surface.i b/python/surface.i
index a6a61e3..4fe7ae3 100644
--- a/python/surface.i
+++ b/python/surface.i
@@ -1,382 +1,387 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 
 %{
 #include <complex>
 #include "map_2d.hh"
 #include "surface.hh"
 #include "types.hh"
 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #include <numpy/arrayobject.h>
 #include "surface_generator_filter_fft.hh"
 #include "surface_generator_voss.hh"
 #include "smart_pointer.hh"
 #include "surface_statistics.hh"
 #include "bem_kato.hh"
 #include "bem_polonski.hh"
 #include "bem_penalty.hh"
 #include "bem_gigi.hh"
 #include "bem_gigipol.hh"
 #include "bem_uzawa.hh"
 #include "bem_fft_base.hh"
 #include "bem_grid.hh"
 #include "bem_grid_polonski.hh"
 #include "bem_grid_kato.hh"
 #include "bem_grid_teboulle.hh"
 #include "bem_grid_condat.hh"
 #include "contact_area.hh"
 #include "contact_cluster_collection.hh"
 #include "functional.hh"
 #include "meta_functional.hh"
 #include "elastic_energy_functional.hh"
 #include "complimentary_term_functional.hh"
 #include "exponential_adhesion_functional.hh"
 #include "squared_exponential_adhesion_functional.hh"
 #include "maugis_adhesion_functional.hh"
 #include "surface_timer.hh"
 #include "fftransform_fftw.hh"
+#include "flood_fill.hh"
 %}
 %include "cpointer.i"
 %include "typemaps.i"
 %include "numpy.i"
 %init %{
   import_array();
 %}
 
  /* -------------------------------------------------------------------------- */
 /* Wrapping returned pointers
 /* -------------------------------------------------------------------------- */
 //%include "typemaps.i"
 %include "std_string.i"
 %include "std_vector.i"
 %include "smart_pointer.hh"
 %template(VectorReal) std::vector<double>;
 //%pointer_functions(int, intp);
 //%pointer_functions(long, longp);
 //%pointer_functions(unsigned int, uintp);
 //%pointer_functions(double, doublep);
 //%pointer_functions(float, floatp);
 //%pointer_functions(std::complex<double>, comlpexp);
 
 %define %my_pointer_class(NAME,TYPE)
 
 %template(SmartPointer##NAME) SmartPointer< TYPE >;
 
 %typemap(out) TYPE & {
 #define TOREPLACE SWIGTYPE_p_SmartPointerT_##NAME##_t
   SmartPointer<TYPE> * ptr = new SmartPointer<TYPE>($1);
   $result = SWIG_NewPointerObj(SWIG_as_voidptr(ptr), TOREPLACE, 0 |  0 );
 #undef TOREPLACE
  }
 %enddef
 
 %my_pointer_class(int,int)
 %my_pointer_class(double,double)
 %my_pointer_class(long,long)
 %my_pointer_class(vectorReal,std::vector<tamaas::Real>)
 
 %apply (unsigned int IN_ARRAY1[ANY]) {(tamaas::UInt n[ANY])};
 
 /* -------------------------------------------------------------------------- */
 /* Wrapping Surface class                                                     */
 /* -------------------------------------------------------------------------- */
 
 %include "tamaas.hh"
 
 %ignore tamaas::Grid::iterator::operator++;
+%ignore tamaas::Grid::begin;
+%ignore tamaas::Grid::end;
 %ignore *::operator=;
 
 %include "grid.hh"
 %include "grid_hermitian.hh"
 namespace tamaas {
   %template (Grid1dReal) Grid<Real, 1>;
   %template (Grid2dReal) Grid<Real, 2>;
+  %template (Grid2dBool) Grid<bool, 2>;
   %template (Grid2dComplex) Grid<complex, 2>;
   %template (Grid2dInt) Grid<Int, 2>;
   %template (Grid2dUInt) Grid<UInt, 2>;
   %template (GridHermitian2dReal) GridHermitian<Real, 2>;
 }
 
 %include "types.hh"
 
 %include "grid_hermitian.hh"
 
 %include "map_2d.hh"
 namespace tamaas {
 %template(Map2dComplex) Map2d<complex>;
 %template(Map2dInt)     Map2d<int>;
 %template(Map2dReal)    Map2d<Real>;
 }
 
 %include "map_2d_square.hh"
 namespace tamaas {
 %template(Map2dSquareComplex) Map2dSquare<complex>;
 %template(Map2dSquareInt)     Map2dSquare<int>;
 %template(Map2dSquareReal)    Map2dSquare<Real>;
 }
 
 %include "surface.hh"
 %include "surface_complex.hh"
 
 namespace tamaas {
 %template(SurfaceReal) Surface<Real>;
 %template(SurfaceInt)  Surface<int>;
 %template(SurfaceRealComplex) SurfaceComplex<Real>;
 }
 
 %inline %{
   __BEGIN_TAMAAS__
     template <int required_type>
     void sanityCheck(PyObject * input) {
       if (PyArray_TYPE((PyArrayObject*)input) != required_type)
 	throw Exception("Sanity check: python object stores wrong type");
     }
   __END_TAMAAS__
 %}
 
 %inline %{
   __BEGIN_TAMAAS__
 
   template <typename T, UInt d>
   class GridForPython : public Grid<T, d> {
     public:
     GridForPython(T * wrapped_memory, UInt n[d], Real L[d], UInt nb_components):
     Grid<T, d>() {
       this->__init__(wrapped_memory, n, nb_components);
     }
 
     void __init__(T * w, UInt n[d], UInt nb_components) {
       memcpy(this->n, n, d * sizeof(UInt));
       memset(this->L, 1, d * sizeof(Real));
       this->nb_components = nb_components;
       this->data.wrapMemory(w, this->computeSize());
     }
 
     GridForPython(PyObject * input): Grid<T, d>() {
       if (!PyArray_Check(input)) {
 	PyObject* obj_type = PyObject_Type(input);
 	std::string type = PyString_AsString(obj_type);
 	throw Exception("GridForPython: incompatible input which is not a numpy: " + type);
       }
 
       UInt _n = PyArray_NDIM((PyArrayObject*)input); // numpy number of dimensions
       npy_intp * ndims = PyArray_DIMS((PyArrayObject*)input); // numpy dimensions
       UInt grid_dims[d] = {0}; // tamaas dims
       if (_n >= d) {
 	for (UInt i = 0 ; i < d ; ++i)
 	  grid_dims[i] = ndims[i];
       } else {
 	throw Exception("GridForPython: incompatible input size for numpy array: " + _n);
       }
 
       if (_n > d) this->nb_components = ndims[d];
       else this->nb_components = 1;
 
       PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input);
       if (iter == NULL) throw;
       this->__init__(( T *)(iter->dataptr), grid_dims, this->nb_components);
       Py_DECREF(iter);
     }
 
     virtual ~GridForPython() {}
   };
 
   __END_TAMAAS__
 %}
 
 %inline %{
   __BEGIN_TAMAAS__
 template <typename T>
   class SurfaceForPython : public Surface<T>{
 
  public:
 
  SurfaceForPython(T * wrapped_memory,
 		   UInt n,
 		   Real L)
    : Surface<T>(0,L){
     this->__init__(wrapped_memory,n,L);
   };
 
 
   void __init__(T * wrapped_memory,
 		UInt n,
 		Real L){
     this->n[0] = n;
     this->n[1] = n;
     this->data.wrapMemory(wrapped_memory,n*n);
   }
 
   SurfaceForPython(PyObject * input, Real L = 1.)
     : Surface<T>(0,L){
 
     if (!PyArray_Check(input)){
       PyObject* obj_type = PyObject_Type(input);
       std::string type = PyString_AsString(obj_type);
       throw Exception("SurfaceForPython: incompatible input which is not a numpy: " + type);
     }
 
     UInt _n = PyArray_NDIM((PyArrayObject*)input);
     if (_n != 2) SURFACE_FATAL("SurfaceForPython: incompatible numpy dimension " << _n);
     npy_intp * ndims = PyArray_DIMS((PyArrayObject*)input);
     UInt sz = ndims[0];
     UInt sz2 = ndims[1];
     if (sz != sz2) SURFACE_FATAL("SurfaceForPython: incompatible numpy dimensions " << sz << " " << sz2);
     PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input);
     if (iter == NULL) throw;
     this->__init__(( T *)(iter->dataptr),sz,L);
     Py_DECREF(iter);
   };
 
   ~SurfaceForPython(){
   };
 
   void resize(UInt new_size){
     if (this->size() != new_size)
       throw Exception("SurfaceForPython: cannot resize a temporary vector");
   }
  };
   __END_TAMAAS__
 %}
 
 %define %grid_typemaps(GRID_CLASS, DATA_TYPE, DATA_TYPECODE, DIM)
 %typemap(out, fragment="NumPy_Fragments") (GRID_CLASS &)
 {
   using namespace tamaas;
   UInt additional_dim = ($1->getNbComponents() == 1)? 0:1;
   UInt dim = DIM+additional_dim;
   npy_intp dims[dim];
   for (UInt i = 0 ; i < DIM ; i++) dims[i] = (npy_intp)$1->sizes()[i];
   if (additional_dim)
     dims[dim-1] = $1->getNbComponents();
 
   DATA_TYPE * data = const_cast<DATA_TYPE *>($1->getInternalData());
   PyObject* obj = PyArray_SimpleNewFromData(dim,
 					    dims,
 					    DATA_TYPECODE,
 					    &data[0]);
   PyArrayObject* array = (PyArrayObject*) obj;
   if (!array) SWIG_fail;
   $result = SWIG_Python_AppendOutput($result, obj);
 }
 
 %typemap(in) tamaas::Grid<DATA_TYPE, DIM> & {
   using namespace tamaas;
   if (!PyArray_Check($input)) {
     PyObject* objectsRepresentation = PyObject_Repr(PyObject_Type($input));
     const char* s = PyString_AsString(objectsRepresentation);
     std::string __type = s;
 
     if ("<class 'tamaas.GridForPython" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else if ("<class 'tamaas.Grid" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else throw tamaas::Exception("typemap: incompatible input which is not a proper type: " + __type + " (DATATYPE = " + std::string("DATA_TYPE") + ")");
   } else {
     $1 = new GridForPython< DATA_TYPE, DIM >($input);
     sanityCheck< DATA_TYPECODE >($input);
   }
 }
 
 %inline %{
   namespace tamaas {
   GRID_CLASS & convertGrid(GRID_CLASS & s){
     return s;
   }
   }
 %}
 %enddef
 
 
 %define %surface_typemaps(SURFACE_CLASS, DATA_TYPE, DATA_TYPECODE)
 
 %typemap(out, fragment="NumPy_Fragments") (SURFACE_CLASS &)
 {
   using namespace tamaas;
   UInt nz = $1->size();
   npy_intp dims[2] = {nz,nz};
   DATA_TYPE * data = const_cast<DATA_TYPE *>($1->getInternalData());
   PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, &data[0]);
   PyArrayObject* array = (PyArrayObject*) obj;
   if (!array) SWIG_fail;
   $result = SWIG_Python_AppendOutput($result, obj);
 }
 
 %typemap(in) tamaas::Surface< DATA_TYPE > & {
   using namespace tamaas;
   if (!PyArray_Check($input)){
     PyObject* objectsRepresentation = PyObject_Repr(PyObject_Type($input));
     const char* s = PyString_AsString(objectsRepresentation);
     std::string __type = s;
 
     if ("<class 'tamaas.SurfaceForPython" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else if ("<class 'tamaas.Surface" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else if ("<class 'tamaas.SurfaceComplex" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else throw tamaas::Exception("typemap: incompatible input which is not a proper type: " + __type + " (DATATYPE = " + std::string("DATA_TYPE") + ")");
   }
   else{
     $1 = new SurfaceForPython< DATA_TYPE >($input);
     sanityCheck< DATA_TYPECODE >($input);
   }
  }
 
 %inline %{
   namespace tamaas {
   SURFACE_CLASS & convertSurface(SURFACE_CLASS & s){
     return s;
   }
   }
 %}
 
 
 %enddef
 
 %surface_typemaps(tamaas::Surface<Real>, Real,       NPY_DOUBLE)
 %surface_typemaps(tamaas::SurfaceComplex<Real>, std::complex<Real>,        NPY_COMPLEX128 )
 %surface_typemaps(tamaas::Surface<unsigned long>, unsigned long, NPY_UINT  )
 %surface_typemaps(tamaas::Surface<UInt>, UInt, NPY_UINT  )
 %surface_typemaps(tamaas::Surface<int>, int, NPY_INT  )
 
 %grid_typemaps(tamaas::Grid1dReal, Real, NPY_DOUBLE, 1)
 %grid_typemaps(tamaas::Grid2dReal, Real, NPY_DOUBLE, 2)
+%grid_typemaps(tamaas::Grid2dBool, bool, NPY_BOOL, 2)
 %grid_typemaps(tamaas::GridHermitian2dReal, tamaas::complex, NPY_COMPLEX128, 2)
 %grid_typemaps(tamaas::Grid2dUInt, UInt, NPY_UINT, 2)
 %grid_typemaps(tamaas::SurfaceComplex<Real>, tamaas::complex, NPY_COMPLEX128, 2)
 
 namespace tamaas {
 %template(SurfaceForPythonReal) SurfaceForPython< Real >;
 }
 
 %include "surface_timer.hh"
diff --git a/src/SConscript b/src/SConscript
index 50842c9..14c1e8a 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,80 +1,81 @@
 Import('main_env')
 
 import os
 
 env = main_env
 
 # Lib roughcontact
 generator_list = """
 surface_generator.cpp
 surface_generator_voss.cpp
 surface_generator_crenel.cpp
 surface_generator_ellipsoid.cpp
 surface_generator_filter.cpp
 surface_generator_filter_bessel.cpp
 surface_generator_filter_fft.cpp
 """.split()
 #env.SharedLibrary('Generator',generator_list)
 
 # Lib SURFACE
 surface_list = """
 map_2d.cpp
 map_2d_square.cpp
 surface.cpp
 surface_timer.cpp
 tamaas.cpp
 """.split()
 #env.SharedLibrary('Surface',surface_list)
 
 # Lib PERCOLATION
 percolation_list = """
 cluster_grow.cpp
 contact_area.cpp
 contact_cluster.cpp
 contact_cluster_collection.cpp
+flood_fill.cpp
 """.split()
 #env.SharedLibrary('Percolation',percolation_list)
 
 # BEM PERCOLATION
 bem_list = """
 tamaas_info.cpp
 bem_kato.cpp
 bem_polonski.cpp
 bem_gigi.cpp
 bem_gigipol.cpp
 bem_penalty.cpp
 bem_uzawa.cpp
 bem_fft_base.cpp
 functional.cpp
 meta_functional.cpp
 elastic_energy_functional.cpp
 exponential_adhesion_functional.cpp
 squared_exponential_adhesion_functional.cpp
 maugis_adhesion_functional.cpp
 complimentary_term_functional.cpp
 fftransform.cpp
 fftransform_fftw.cpp
 fft_plan_manager.cpp
 grid.cpp
 grid_hermitian.cpp
 types.cpp
 bem_grid.cpp
 bem_grid_polonski.cpp
 bem_grid_kato.cpp
 bem_grid_teboulle.cpp
 bem_grid_condat.cpp
 """.split()
 #env.SharedLibrary('BEM',bem_list)
 
 rough_contact_list = generator_list + surface_list + percolation_list + bem_list
 current_files = set(os.listdir(Dir('.').srcnode().abspath))
 unexpected_files = list(current_files - set(rough_contact_list + ['SConscript']))
 
 unexpected_headers = [ f for f in unexpected_files if os.path.splitext(f)[1] == '.hh' ]
 
 unexpected_files = [ f for f in unexpected_files if not os.path.splitext(f)[1] == '.hh' ]
 
 #print "Unexpected files : {}".format(unexpected_files)
 #print "Include path : {}".format([str(x) for x in env['CPPPATH']])
 #print unexpected_headers
 env.SharedLibrary('Tamaas',rough_contact_list)
diff --git a/src/flood_fill.cpp b/src/flood_fill.cpp
new file mode 100644
index 0000000..d259c4f
--- /dev/null
+++ b/src/flood_fill.cpp
@@ -0,0 +1,105 @@
+/**
+ *
+ * @author Lucas Frérot <lucas.frerot@epfl.ch>
+ *
+ * @section LICENSE
+ *
+ * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
+ * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
+ * Solides)
+ *
+ * Tamaas is free  software: you can redistribute it and/or  modify it under the
+ * terms  of the  GNU Lesser  General Public  License as  published by  the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
+ * details.
+ *
+ * You should  have received  a copy  of the GNU  Lesser General  Public License
+ * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#include "flood_fill.hh"
+#include <stack>
+/* -------------------------------------------------------------------------- */
+__BEGIN_TAMAAS__
+
+Cluster::Cluster(const Cluster & other):
+  points(other.points),
+  area(other.area)
+{}
+
+Cluster::Cluster(Point start,
+		 const Grid<bool, 2> & contact,
+		 Grid<bool, 2> & visited,
+		 bool diagonal):area(0) {
+  // Visiting stack
+  std::stack<Point> visiting;
+  visiting.push(start);
+  // Contact sizes
+  const UInt * n = contact.sizes();
+
+  while (!visiting.empty()) {
+    Int i = visiting.top().first;
+    Int j = visiting.top().second;
+
+    if (visited(i, j)) {
+      visiting.pop();
+      continue;
+    }
+
+    visited(i, j) = true;
+    points.push_back(visiting.top());
+    area++;
+    visiting.pop();
+
+    std::list<Point> neighbors;
+    neighbors.push_back(std::make_pair(i+1, j));
+    neighbors.push_back(std::make_pair(i-1, j));
+    neighbors.push_back(std::make_pair(i, j-1));
+    neighbors.push_back(std::make_pair(i, j+1));
+      
+    if (diagonal) {
+      neighbors.push_back(std::make_pair(i+1, j+1));
+      neighbors.push_back(std::make_pair(i-1, j-1));
+      neighbors.push_back(std::make_pair(i-1, j+1));
+      neighbors.push_back(std::make_pair(i+1, j-1));
+    }
+      
+    for (auto & p : neighbors) {
+      Int im = p.first % n[0];
+      Int jm = p.second % n[1];
+      if (!visited(im, jm) && contact(im, jm)) {
+	visiting.push(p);
+      }
+    }
+  }
+}
+
+Cluster::~Cluster()
+{}
+
+std::list<Cluster> FloodFill::getClusters(const Grid<bool, 2> & contact,
+					  bool diagonal) {
+  const UInt * n = contact.sizes();
+  Grid<bool, 2> visited(n, contact.lengths(), 1);
+  visited = false;
+  std::list<Cluster> clusters;
+
+  for (UInt i = 0 ; i < n[0] ; ++i) {
+    for (UInt j = 0 ; j < n[1] ; ++j) {
+      if (contact(i, j) && !visited(i, j)) {
+	clusters.push_back(Cluster(std::make_pair(i, j), contact, visited, diagonal));
+      }
+    }
+  }
+
+  return clusters;
+}
+
+
+__END_TAMAAS__
diff --git a/src/flood_fill.hh b/src/flood_fill.hh
new file mode 100644
index 0000000..d8dcad7
--- /dev/null
+++ b/src/flood_fill.hh
@@ -0,0 +1,73 @@
+/**
+ *
+ * @author Lucas Frérot <lucas.frerot@epfl.ch>
+ *
+ * @section LICENSE
+ *
+ * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
+
+ * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
+ * Solides)
+ *
+ * Tamaas is free  software: you can redistribute it and/or  modify it under the
+ * terms  of the  GNU Lesser  General Public  License as  published by  the Free
+ * Software Foundation, either version 3 of the License, or (at your option) any
+ * later version.
+ *
+ * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
+ * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+ * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
+ * details.
+ *
+ * You should  have received  a copy  of the GNU  Lesser General  Public License
+ * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#ifndef FLOOD_FILL_H
+#define FLOOD_FILL_H
+/* -------------------------------------------------------------------------- */
+#include "grid.hh"
+#include <list>
+/* -------------------------------------------------------------------------- */
+__BEGIN_TAMAAS__
+/* -------------------------------------------------------------------------- */
+typedef std::pair<Int, Int> Point;
+typedef Grid<bool, 2> Grid2dBool;
+
+class Cluster {
+
+public:
+  /// Constructor
+  Cluster(Point start,
+	  const Grid<bool, 2> & contact,
+	  Grid<bool, 2> & visited,
+	  bool diagonal);
+  /// Copy constructor
+  Cluster(const Cluster & other);
+  /// Default constructor
+  Cluster() {};
+  /// Destructor
+  virtual ~Cluster();
+  /// Get area of cluster
+  UInt getArea() {return area;}
+
+private:
+  std::list<Point> points;
+  UInt area;
+};
+
+/* -------------------------------------------------------------------------- */
+
+class FloodFill {
+public:
+  /// Returns a list of clusters from a boolean contact map
+  static std::list<Cluster> getClusters(const Grid<bool, 2> & contact,
+					bool diagonal);
+};
+
+
+/* -------------------------------------------------------------------------- */
+__END_TAMAAS__
+/* -------------------------------------------------------------------------- */
+#endif // FLOOD_FILL_H
diff --git a/src/grid.hh b/src/grid.hh
index f68415b..d47e507 100644
--- a/src/grid.hh
+++ b/src/grid.hh
@@ -1,311 +1,342 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_HH__
 #define __GRID_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "array.hh"
 #include <cstddef>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /**
  * @brief Multi-dimensional & multi-component array class
  *
  * This class is a container for multi-component data stored on a multi-
  * dimensional grid.
  *
  * The access function is the parenthesis operator. For a grid of dimension d,
  * the operator takes d+1 arguments: the first d arguments are the position on
  * the grid and the last one is the component of the stored data.
  *
  * It is also possible to use the access operator with only one argument, it is
  * then considering the grid as a flat array, accessing the given cell of the
  * array.
  */
 template <typename T, UInt dim>
 class Grid {
 /* -------------------------------------------------------------------------- */
 /* Constructors                                                               */
 /* -------------------------------------------------------------------------- */
 public:
   /// Constructor by default (empty array)
   Grid();
   /// Constructor
   Grid(const UInt n[dim], const Real L[dim], UInt nb_components);
   /// Destructor
   virtual ~Grid();
 
 /* -------------------------------------------------------------------------- */
 /* Iterator class                                                             */
 /* -------------------------------------------------------------------------- */
   struct iterator {
     /// constructor
     iterator(T * start, ptrdiff_t offset):ptr(start), offset(offset) {}
     /// destructor
     ~iterator() {}
     /// pre-increment
     iterator & operator++() { ptr += offset; return *this; }
     /// comparison
     bool operator<(const iterator & a) { return ptr < a.ptr; }
     /// increment with given offset
     iterator & operator+=(ptrdiff_t a) { ptr += a*offset; return *this;}
     /// dereference iterator
     T * operator*() { return ptr; }
     /// needed for OpenMP range calculations
     Int operator-(const iterator & a) const { return (ptr - a.ptr)/offset; }
   private:
     T * ptr;
     ptrdiff_t offset;
   };
 
 public:
   /// Begin iterator with n components
   iterator begin(UInt n = 1) { return iterator(data.getPtr(), n); }
   /// End iterator with n components
   iterator end(UInt n = 1) { return iterator(data.getPtr()+computeSize(), n); }
 
 public:
 /* -------------------------------------------------------------------------- */
 /* Common operations                                                          */
 /* -------------------------------------------------------------------------- */
   /// Resize array
   void resize(const UInt n[dim]);
   /// Compute size
   inline UInt computeSize() const;
   /// Print
   void printself(std::ostream & str) const;
   /// Get sizes
   const UInt * sizes() const {return n;}
   /// Get total size
   UInt dataSize() const {return data.size(); }
   /// Get number of points
   UInt getNbPoints() const {return dataSize() / getNbComponents(); }
   /// Get lengths
   const Real * lengths() const {return L;}
   /// Get number of components
   UInt getNbComponents() const {return nb_components;}
   /// Set number of components
   void setNbComponents(UInt n) {nb_components = n;}
   /// Get internal data pointer (const)
   const T * getInternalData() const {return data.getPtr();}
   /// Get internal data pointer (non-const)
   T * getInternalData() {return data.getPtr();}
   /// Set components
   void uniformSetComponents(const Grid<T, 1> & vec);
 
 /* -------------------------------------------------------------------------- */
 /* Access operators                                                           */
 /* -------------------------------------------------------------------------- */
+  /// Variadic access operator (non-const)
   template <typename... T1>
   inline T & operator()(T1... args);
+  /// Variadic access operator
   template <typename... T1>
   inline const T & operator()(T1... args) const;
 
+  /// Tuple index access operator
   inline T & operator()(UInt tuple[dim+1]);
   inline const T & operator()(UInt tuple[dim+1]) const;
 
 private:
+  /// Unpacking the arguments of variadic () operator when nargs == dim+1
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index, T1... rest) const;
+  /// End case for recursion
   template <typename... T1>
   inline UInt unpackOffset(UInt offset, UInt index) const;
+  template <typename... T1>
+  /// Unpacking the arguments of variadic () operator when nargs == dim (for nb_components == 1)
+  inline UInt unpackOffsetReduced(UInt offset, UInt index, T1... rest) const;
+  template <typename... T1>
+  /// End case for recursion
+  inline UInt unpackOffsetReduced(UInt offset, UInt index) const;
+  /// Computing offset for a tuple index
   inline UInt computeOffset(UInt tuple[dim+1]) const;
 
 /* -------------------------------------------------------------------------- */
 /* Arithmetic operators                                                       */
 /* -------------------------------------------------------------------------- */
 public:
 #define GRID_VEC_OPERATOR(op)					\
   template <typename T1>					\
   void operator op(const Grid<T1, dim> & other)			\
 
   GRID_VEC_OPERATOR(+=);
   GRID_VEC_OPERATOR(*=);
   GRID_VEC_OPERATOR(-=);
   GRID_VEC_OPERATOR(/=);
 #undef GRID_VEC_OPERATOR
 
 #define GRID_SCALAR_OPERATOR(op)		\
   void operator op(const T & e)			\
 
   GRID_SCALAR_OPERATOR(+=);
   GRID_SCALAR_OPERATOR(*=);
   GRID_SCALAR_OPERATOR(-=);
   GRID_SCALAR_OPERATOR(/=);
   GRID_SCALAR_OPERATOR(=);
 #undef GRID_SCALAR_OPERATOR
 
   // = operator
   void operator=(const Grid<T, dim> & other);
   // = operator (not const input, otherwise gcc is confused)
   void operator=(Grid<T, dim> & other);
 
   // Copy data from another grid
   template <typename T1>
   void copy(const Grid<T1, dim> & other);
 
 /* -------------------------------------------------------------------------- */
 /* Member variables                                                           */
 /* -------------------------------------------------------------------------- */
 protected:
   UInt n[dim];
   Real L[dim];
   Array<T> data;
   UInt nb_components;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline function definitions                                                */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline UInt Grid<T, dim>::computeSize() const {
   UInt size = 1;
   for (UInt i = 0 ; i < dim ; i++) size *= n[i];
   size *= nb_components;
   return size;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index, T1... rest) const {
   UInt size = sizeof...(T1);
   offset += index;
   offset *= (size == 1) ? this->nb_components : this->n[dim-size+1];
   return unpackOffset(offset, rest...); // tail-rec bb
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline UInt Grid<T, dim>::unpackOffset(UInt offset, UInt index) const {
   return offset + index;
 }
 
+template <typename T, UInt dim>
+template <typename... T1>
+inline UInt Grid<T, dim>::unpackOffsetReduced(UInt offset, UInt index, T1... rest) const {
+  UInt size = sizeof...(T1);
+  offset += index;
+  offset *= this->n[dim-size-1];
+  return unpackOffsetReduced(offset, rest...);
+}
+
+template <typename T, UInt dim>
+template <typename... T1>
+inline UInt Grid<T, dim>::unpackOffsetReduced(UInt offset, UInt index) const {
+  return unpackOffset(offset, index);
+}
+
 template <typename T, UInt dim>
 template <typename... T1>
 inline T & Grid<T, dim>::operator()(T1... args) {
-  TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1,
+  TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1 ||
+		(sizeof...(T1) == dim && this->nb_components == 1),
 		"operator() does not match dimension");
-  UInt offset = unpackOffset(0, args...);
+  UInt offset = (sizeof...(T1) == dim) ?
+    unpackOffsetReduced(0, args...) : unpackOffset(0, args...);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline const T & Grid<T, dim>::operator()(T1... args) const {
-  TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1,
+  TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1 ||
+		(sizeof...(T1) == dim && this->nb_components == 1),
 		"operator() does not match dimension");
-  UInt offset = unpackOffset(0, args...);
+  UInt offset = (sizeof...(T1) == dim) ?
+    unpackOffsetReduced(0, args...) : unpackOffset(0, args...);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 inline UInt Grid<T, dim>::computeOffset(UInt tuple[dim+1]) const {
   UInt offset = 0;
   for (UInt i = 0 ; i < dim-1 ; i++) {
     offset += tuple[i];
     offset *= n[i+1];
   }
 
   offset += tuple[dim-1];
   offset *= nb_components;
   return offset + tuple[dim];
 }
 
 template <typename T, UInt dim>
 inline T & Grid<T, dim>::operator()(UInt tuple[dim+1]) {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 template <typename T, UInt dim>
 inline const T & Grid<T, dim>::operator()(UInt tuple[dim+1]) const {
   UInt offset = computeOffset(tuple);
   return this->data[offset];
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define GRID_VEC_OPERATOR_IMPL(op)					\
   template <typename T, UInt dim>					\
   template <typename T1>						\
   inline void Grid<T, dim>::operator op(const Grid<T1, dim> & other) {	\
     TAMAAS_ASSERT(other.computeSize() == data.size(), "surface size does not match");\
     _Pragma("omp parallel for")						\
     for (UInt i = 0 ; i < this->data.size() ; i++) {			\
       data[i] op other(i);						\
     }									\
   }									\
 
 GRID_VEC_OPERATOR_IMPL(+=);
 GRID_VEC_OPERATOR_IMPL(-=);
 GRID_VEC_OPERATOR_IMPL(*=);
 GRID_VEC_OPERATOR_IMPL(/=);
 #undef GRID_VEC_OPERATOR
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::operator=(const Grid<T, dim> & other) {
   this->copy(other);
 }
 
 template <typename T, UInt dim>
 void Grid<T, dim>::operator=(Grid<T, dim> & other) {
   this->copy(other);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename T1>
 void Grid<T, dim>::copy(const Grid<T1, dim> & other) {
   resize(other.sizes());
 #pragma omp parallel for
   for (UInt i = 0 ; i < data.size() ; i++)
     data[i] = other(i);
 }
 
 /* -------------------------------------------------------------------------- */
 /* Stream operator                                                            */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 inline std::ostream & operator<<(std::ostream & stream,
 				 const Grid<T, dim> & _this) {
   _this.printself(stream);
   return stream;
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #endif // __GRID_HH__
diff --git a/tests/SConscript b/tests/SConscript
index e1fb27b..c6fcf60 100644
--- a/tests/SConscript
+++ b/tests/SConscript
@@ -1,24 +1,25 @@
 from os.path import join
 
 Import('main_env')
 
 test_files = Split("""
 run_tests.sh
 test_hertz_pressure.py
 test_westergaard.py
 test_surface.py
 test_autocorrelation.py
 test_hertz_disp.py
 test_hertz_kato.py
 test_hertz_adhesion.py
 test_saturated_pressure.py
 test_fftransform.py
 test_bem_grid.py
+test_flood_fill.py
 """)
 
 src_dir = "#/tests"
 build_dir = 'build-' + main_env['build_type'] + '/tests'
 
 for file in test_files:
     source = join(src_dir, file)
     Command(file, source, Copy("$TARGET", "$SOURCE"))
diff --git a/tests/test_flood_fill.py b/tests/test_flood_fill.py
new file mode 100644
index 0000000..e1929c3
--- /dev/null
+++ b/tests/test_flood_fill.py
@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+# coding: utf-8
+# -----------------------------------------------------------------------------
+# @author Lucas Frérot <lucas.frerot@epfl.ch>
+#
+# @section LICENSE
+#
+# Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
+# Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
+# Solides)
+#
+# Tamaas is free  software: you can redistribute it and/or  modify it under the
+# terms  of the  GNU Lesser  General Public  License as  published by  the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
+# details.
+#
+# You should  have received  a copy  of the GNU  Lesser General  Public License
+# along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
+# -----------------------------------------------------------------------------
+from __future__ import print_function, division
+
+import sys
+import tamaas as tm
+import numpy as np
+
+field = np.zeros((18, 18))
+
+# Cluster of permeter 18 area 13
+field[2:4, 3:6] = 1.
+field[4  , 4:6] = 1.
+field[5  , 5  ] = 1.
+field[3:5, 6:8] = 1.
+
+# Cluster of perimeter 8 area 4
+field[14:16, 3:5] = 1.
+
+# Cluster of perimeter 20 area 11
+field[9:11, 8:11 ] = 1.
+field[7:9 , 9    ] = 1.
+field[10  , 11:14] = 1.
+
+# Cluster of perimeter 18 area 9
+field[3:5, 14:16] = 1.
+field[5:10, 15] = 1.
+
+bool_field = np.zeros_like(field, np.bool)
+bool_field[field > 0] = True
+
+clusters = tm.FloodFill.getClusters(bool_field, False)
+
+
+# Dummy class for clusters
+class dummy:
+    def __init__(self, area, perimeter):
+        self.area = area
+        self.perimeter = perimeter
+
+
+# Solution
+ref = [dummy(13, 18), dummy(9, 18), dummy(11, 20), dummy(4, 8)]
+
+for x, y in zip(clusters, ref):
+    if x.getArea() != y.area:
+        sys.exit(1)
diff --git a/tests/test_surface.py b/tests/test_surface.py
index aff8852..f7efd80 100644
--- a/tests/test_surface.py
+++ b/tests/test_surface.py
@@ -1,123 +1,147 @@
+#!/usr/bin/env python
+# coding: utf-8
+# -----------------------------------------------------------------------------
+# @author Lucas Frérot <lucas.frerot@epfl.ch>
+#
+# @section LICENSE
+#
+# Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
+# Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
+# Solides)
+#
+# Tamaas is free  software: you can redistribute it and/or  modify it under the
+# terms  of the  GNU Lesser  General Public  License as  published by  the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
+# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
+# A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
+# details.
+#
+# You should  have received  a copy  of the GNU  Lesser General  Public License
+# along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
+# -----------------------------------------------------------------------------
 from __future__ import print_function, division
 
 import sys
 import numpy as np
 import tamaas as tm
 
 
 def compute_mse(values, mean):
     return np.mean((values - mean)**2) / values.size
 
 
 def main():
     tm.initialize(1)
 
     # Spectrum parameters in wavelengths
     disctretization_size = 1.
     lambda_s = 8.
     lambda_r = 64.
     lambda_l = 64.
     box_size = 512
 
     # Independent parameters
     Hurst = 0.8
     hrms = 1.
 
     # Spectrum parameters in wavenumbers
     k_l = int(box_size // lambda_l)
     k_r = int(box_size // lambda_r)
     k_s = int(box_size // lambda_s)
     N = int(box_size // disctretization_size)
 
     # Analytical moments, cf. Yastrebov et al. (2015)
     # "From infinitesimal to full contact between rough surfaces: Evolution
     # of the contact area", appendix A
     xi = k_l / k_r
     zeta = k_s / k_r
     H = Hurst
     T = {0: 2*np.pi,
          2: np.pi,
          4: 3*np.pi/4}
 
     m = {}
     for q in [0, 2, 4]:
         m[q] = T[q] * k_r**(q-2*H)*((1-xi**(q+2))/(q+2) + (zeta**(q-2*H)-1)/(q-2*H))
 
     # Some spectral parameters
     alpha = m[0]*m[4]/m[2]**2
     rms_slopes = np.sqrt(np.pi * m[2] / 2)
 
     # Output setup
     realizations = 100
     hrms_out = np.zeros(realizations)
     alpha_out = np.zeros(realizations)
     rms_slopes_out = {"spectral": np.zeros(realizations),
                       "geometric": np.zeros(realizations)}
 
     for i in range(realizations):
         # Surface generator setup
         sg = tm.SurfaceGeneratorFilterFFT()
         sg.getGridSize().assign(N)
         sg.getHurst().assign(Hurst)
         sg.getRMS().assign(hrms)
         sg.getQ0().assign(k_l)
         sg.getQ1().assign(k_r)
         sg.getQ2().assign(k_s)
         sg.getRandomSeed().assign(i)
         sg.Init()
         # Generating surface
         surface = sg.buildSurface() / sg.analyticalRMS()
 
         # Computing RMS heights
         # We know surface.mean() = 0, so we get an unbiased estimator here
         hrms_out[i] = np.sqrt(1 / (N**2-1) * np.sum(surface**2))
 
         # Computing spectrum bandwidth
         surface /= box_size**2  # Normalize so that the computed PSD does not contain the box size
         moments = tm.SurfaceStatistics.computeMoments(surface)
         m0 = hrms_out[i]**2 / box_size**4
         m2 = moments[0]
         m4 = moments[1]
         alpha_out[i] = m0*m4/m2**2
 
         # Computing RMS slopes
         rms_slopes_out["spectral"][i] = tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
         rms_slopes_out["geometric"][i] = tm.SurfaceStatistics.computeRMSSlope(surface)
 
     # print(rms_slopes_out)
     # print(rms_slopes)
     # Computing Mean Square Error
     print("""========= Analytical values =========
 hrms = {}
 alpha = {}
 m0, m2, m4 = {}
 rms_slopes = {}""".format(hrms, alpha, m, rms_slopes))
     print("========= MSE Calculations  =========")
 
     data_association = [
         (hrms_out, hrms, "hrms"),
         (alpha_out, alpha, "alpha"),
         (rms_slopes_out["spectral"], rms_slopes, "RMS Slopes Spectral"),
         (rms_slopes_out["geometric"], rms_slopes, "RMS Slopes Geometric")
     ]
 
     # Computing and printing MSE values
     mse_vals = []
     for out, analytical, string in data_association:
         mse = np.sqrt(compute_mse(out, analytical)) / analytical
         print("Normalized error (MSE) for {}: {}".format(string, mse))
         mse_vals.append(mse)
 
     # Checking tolerances
     tols = [5e-3, 5e-3, 1e-1, 1e-1]
     for mse, tol in zip(mse_vals, tols):
         if mse >= tol:
             print("{} >= {}".format(mse, tol))
             return 1
 
     tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())