diff --git a/python/wrap/percolation.cpp b/python/wrap/percolation.cpp
index fd1aa79..cfff21c 100644
--- a/python/wrap/percolation.cpp
+++ b/python/wrap/percolation.cpp
@@ -1,93 +1,95 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "flood_fill.hh"
 #include "wrap.hh"
 /* -------------------------------------------------------------------------- */
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 namespace wrap {
 
 /* -------------------------------------------------------------------------- */
 using namespace py::literals;
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 void wrapCluster(py::module& mod) {
   auto name = makeDimensionName("Cluster", dim);
   py::class_<Cluster<dim>>(mod, name.c_str())
       .def(py::init<>())
       .def_property_readonly("area", &Cluster<dim>::getArea, "Area of cluster")
       .def_property_readonly("points", &Cluster<dim>::getPoints,
                              "Get list of points of cluster")
       .def_property_readonly("perimeter", &Cluster<dim>::getPerimeter,
                              "Get perimeter of cluster")
       .def_property_readonly("bounding_box", &Cluster<dim>::boundingBox,
                              "Compute the bounding box of a cluster")
+      .def_property_readonly("extent", &Cluster<dim>::extent,
+                             "Compute the extents of a cluster")
 
       .def(TAMAAS_DEPRECATE_ACCESSOR(getArea, Cluster<dim>, "area"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getPoints, Cluster<dim>, "points"))
       .def(TAMAAS_DEPRECATE_ACCESSOR(getPerimeter, Cluster<dim>, "perimeter"))
 
       .def("__str__", [](const Cluster<dim>& cluster) {
         std::stringstream sstr;
         sstr << '{';
         auto&& points = cluster.getPoints();
         auto print_point = [&sstr](const auto& point) {
           sstr << '(';
           for (UInt i = 0; i < point.size() - 1; ++i)
             sstr << point[i] << ", ";
           sstr << point.back() << ')';
         };
 
         auto point_it = points.begin();
         for (UInt p = 0; p < points.size() - 1; ++p) {
           print_point(*(point_it++));
           sstr << ", ";
         }
         print_point(points.back());
         sstr << "}";
         return sstr.str();
       });
 }
 
 void wrapPercolation(py::module& mod) {
   wrapCluster<1>(mod);
   wrapCluster<2>(mod);
   wrapCluster<3>(mod);
 
   py::class_<FloodFill>(mod, "FloodFill")
       .def_static("getSegments", &FloodFill::getSegments,
                   "Return a list of segments from boolean map", "contact"_a)
       .def_static("getClusters", &FloodFill::getClusters,
                   "Return a list of clusters from boolean map", "contact"_a,
                   "diagonal"_a)
       .def_static("getVolumes", &FloodFill::getVolumes,
                   "Return a list of volume clusters", "map"_a, "diagonal"_a);
 }
 
 }  // namespace wrap
 
 }  // namespace tamaas
diff --git a/src/percolation/flood_fill.cpp b/src/percolation/flood_fill.cpp
index 0228f05..69781ee 100644
--- a/src/percolation/flood_fill.cpp
+++ b/src/percolation/flood_fill.cpp
@@ -1,258 +1,271 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "flood_fill.hh"
 #include "partitioner.hh"
 #include <algorithm>
 #include <limits>
 #include <stack>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
-Int unsigned_modulo(Int i, UInt n) { return (i % n + n) % n; }
+Int wrap_pbc(Int i, Int n) {
+  if (i < 0)
+    return i + n;
+  if (i >= n)
+    return i - n;
+  return i;
+}
 
 template <std::size_t dim>
-std::array<UInt, dim> unsigned_modulo(const std::array<Int, dim>& t,
-                                      const std::array<UInt, dim>& n) {
+std::array<UInt, dim> wrap_pbc(const std::array<Int, dim>& t,
+                               const std::array<Int, dim>& n) {
   std::array<UInt, dim> index;
   for (UInt i = 0; i < dim; ++i)
-    index[i] = unsigned_modulo(t[i], n[i]);
+    index[i] = static_cast<UInt>(wrap_pbc(t[i], n[i]));
   return index;
 }
 
 template <std::size_t dim>
 std::array<Int, dim> cast_int(const std::array<UInt, dim>& a) {
   std::array<Int, dim> r;
   for (UInt i = 0; i < dim; ++i)
-    r[i] = (Int)a[i];
+    r[i] = static_cast<Int>(a[i]);
   return r;
 }
 
 template <UInt dim>
 Cluster<dim>::Cluster(const Cluster& other)
     : points(other.points), perimeter(other.perimeter) {}
 
 template <>
-auto Cluster<1>::getNextNeighbors(const std::array<UInt, 1>& p) {
+auto Cluster<1>::getNextNeighbors(const std::array<Int, 1>& p) {
   std::vector<Point> neighbors(2);
   Int i = p[0];
   neighbors[0] = {i - 1};
   neighbors[1] = {i + 1};
   return neighbors;
 }
 
 template <>
-auto Cluster<1>::getDiagonalNeighbors(const std::array<UInt, 1>& /*p*/) {
+auto Cluster<1>::getDiagonalNeighbors(const std::array<Int, 1>& /*p*/) {
   return std::vector<Point>();
 }
 
 template <>
-auto Cluster<2>::getNextNeighbors(const std::array<UInt, 2>& p) {
+auto Cluster<2>::getNextNeighbors(const std::array<Int, 2>& p) {
   std::vector<Point> neighbors(4);
   Int i = p[0], j = p[1];
   neighbors[0] = {i + 1, j};
   neighbors[1] = {i - 1, j};
   neighbors[2] = {i, j - 1};
   neighbors[3] = {i, j + 1};
   return neighbors;
 }
 
 template <>
-auto Cluster<2>::getDiagonalNeighbors(const std::array<UInt, 2>& p) {
+auto Cluster<2>::getDiagonalNeighbors(const std::array<Int, 2>& p) {
   std::vector<Point> neighbors(4);
   Int i = p[0], j = p[1];
   neighbors[0] = {i + 1, j + 1};
   neighbors[1] = {i - 1, j - 1};
   neighbors[2] = {i - 1, j + 1};
   neighbors[3] = {i + 1, j - 1};
   return neighbors;
 }
 
 template <>
-auto Cluster<3>::getNextNeighbors(const std::array<UInt, 3>& p) {
+auto Cluster<3>::getNextNeighbors(const std::array<Int, 3>& p) {
   std::vector<Point> neighbors(6);
   Int i = p[0], j = p[1], k = p[2];
   neighbors[0] = {i + 1, j, k};
   neighbors[1] = {i - 1, j, k};
   neighbors[2] = {i, j - 1, k};
   neighbors[3] = {i, j + 1, k};
   neighbors[4] = {i, j, k - 1};
   neighbors[5] = {i, j, k + 1};
   return neighbors;
 }
 
 template <>
-auto Cluster<3>::getDiagonalNeighbors(const std::array<UInt, 3>& p) {
+auto Cluster<3>::getDiagonalNeighbors(const std::array<Int, 3>& p) {
   std::vector<Point> neighbors(20);
   Int i = p[0], j = p[1], k = p[2];
   // 8 corners
   neighbors[0] = {i + 1, j + 1, k + 1};
   neighbors[1] = {i + 1, j + 1, k - 1};
   neighbors[2] = {i + 1, j - 1, k + 1};
   neighbors[3] = {i + 1, j - 1, k - 1};
   neighbors[4] = {i - 1, j + 1, k + 1};
   neighbors[5] = {i - 1, j + 1, k - 1};
   neighbors[6] = {i - 1, j - 1, k + 1};
   neighbors[7] = {i - 1, j - 1, k - 1};
 
   // 4 diagonals in i = 0
   neighbors[8] = {i, j + 1, k + 1};
   neighbors[9] = {i, j + 1, k - 1};
   neighbors[10] = {i, j - 1, k + 1};
   neighbors[11] = {i, j - 1, k - 1};
 
   // 4 diagonals in j = 0
   neighbors[12] = {i + 1, j, k + 1};
   neighbors[13] = {i + 1, j, k - 1};
   neighbors[14] = {i - 1, j, k + 1};
   neighbors[15] = {i - 1, j, k - 1};
 
   // 4 diagonals in k = 0
   neighbors[16] = {i + 1, j + 1, k};
   neighbors[17] = {i + 1, j - 1, k};
   neighbors[18] = {i - 1, j + 1, k};
   neighbors[19] = {i - 1, j - 1, k};
   return neighbors;
 }
 
 template <UInt dim>
-std::pair<std::array<Int, dim>, std::array<Int, dim>>
-Cluster<dim>::boundingBox() const {
-  // TODO make sure periodic boundaries are correctly unwrapped
+typename Cluster<dim>::BBox Cluster<dim>::boundingBox() const {
   std::array<Int, dim> mins, maxs;
   mins.fill(std::numeric_limits<Int>::max());
   maxs.fill(std::numeric_limits<Int>::min());
 
   for (const auto& p : points)
     for (UInt i = 0; i < dim; ++i) {
       mins[i] = std::min(mins[i], p[i]);
       maxs[i] = std::max(maxs[i], p[i]);
     }
 
   return std::make_pair(mins, maxs);
 }
 
+template <UInt dim>
+std::array<Int, dim> Cluster<dim>::extent() const {
+  const auto bb = boundingBox();
+  std::array<Int, dim> extent;
+  std::transform(bb.first.begin(), bb.first.end(), bb.second.begin(),
+                 extent.begin(),
+                 [](auto min, auto max) { return max - min + 1; });
+  return extent;
+}
+
 template <UInt dim>
 Cluster<dim>::Cluster(Point start, const Grid<bool, dim>& map,
                       Grid<bool, dim>& visited, bool diagonal) {
   // Visiting stack
   std::stack<Point, std::vector<Point>> visiting;
   visiting.push(start);
   // Contact sizes
-  const auto& n = map.sizes();
+  const auto& n = cast_int(map.sizes());
 
   while (not visiting.empty()) {
-    auto p = unsigned_modulo(visiting.top(), n);
+    auto p = visiting.top();
+    auto p_wrapped = wrap_pbc(p, n);
 
-    if (visited(p)) {
+    if (visited(p_wrapped)) {
       visiting.pop();
       continue;
     }
 
-    visited(p) = true;
+    visited(p_wrapped) = true;
     points.push_back(visiting.top());
     visiting.pop();
 
     auto process = [&](const std::vector<Point>& neighbors,
                        const bool do_perimeter) {
-      for (auto& p : neighbors) {
-        auto index = unsigned_modulo(p, n);
+      for (const auto& p : neighbors) {
+        const auto index = wrap_pbc(p, n);
         if (not visited(index) and map(index)) {
-          visiting.push(cast_int(index));
+          visiting.push(p);
         }
 
-        if (do_perimeter and (not map(index)))
-          ++perimeter;
+        perimeter += do_perimeter and (not map(index));
       }
     };
 
     process(getNextNeighbors(p), true);
 
     if (diagonal) {
       process(getDiagonalNeighbors(p), false);
     }
   }
 }
 
-std::list<Cluster<1>> FloodFill::getSegments(const Grid<bool, 1>& map) {
+typename FloodFill::List<1> FloodFill::getSegments(const Grid<bool, 1>& map) {
   auto n = map.sizes();
   Grid<bool, 1> visited(n, 1);
   visited = false;
-  std::list<Cluster<1>> clusters;
+  List<1> clusters;
 
   for (UInt i = 0; i < n[0]; ++i) {
     if (map(i) && !visited(i)) {
-      clusters.emplace_back(std::array<Int, 1>{(Int)i}, map, visited,
-                            false);
+      clusters.emplace_back(std::array<Int, 1>{(Int)i}, map, visited, false);
     }
   }
 
   return clusters;
 }
 
-std::list<Cluster<2>> FloodFill::getClusters(const Grid<bool, 2>& map,
-                                             bool diagonal) {
+typename FloodFill::List<2> FloodFill::getClusters(const Grid<bool, 2>& map,
+                                                   bool diagonal) {
   auto global_contact = Partitioner<2>::gather(map);
   auto n = global_contact.sizes();
   Grid<bool, 2> visited(n, 1);
   visited = false;
-  std::list<Cluster<2>> clusters;
+  List<2> clusters;
 
   if (mpi::rank() == 0)
     for (UInt i = 0; i < n[0]; ++i) {
       for (UInt j = 0; j < n[1]; ++j) {
         if (global_contact(i, j) && !visited(i, j)) {
           clusters.emplace_back(std::array<Int, 2>{(Int)i, (Int)j},
                                 global_contact, visited, diagonal);
         }
       }
     }
 
   return clusters;
 }
 
-std::list<Cluster<3>> FloodFill::getVolumes(const Grid<bool, 3>& map,
-                                            bool diagonal) {
+typename FloodFill::List<3> FloodFill::getVolumes(const Grid<bool, 3>& map,
+                                                  bool diagonal) {
   auto n = map.sizes();
   Grid<bool, 3> visited(n, 1);
   visited = false;
-  std::list<Cluster<3>> clusters;
+  List<3> clusters;
 
   for (UInt i = 0; i < n[0]; ++i) {
     for (UInt j = 0; j < n[1]; ++j) {
       for (UInt k = 0; k < n[2]; ++k) {
         if (map(i, j, k) && !visited(i, j, k)) {
-          clusters.emplace_back(std::array<Int, 3>{(Int)i, (Int)j, (Int)k},
-                                map, visited, diagonal);
+          clusters.emplace_back(std::array<Int, 3>{(Int)i, (Int)j, (Int)k}, map,
+                                visited, diagonal);
         }
       }
     }
   }
 
   return clusters;
 }
 
 template class Cluster<1>;
 template class Cluster<2>;
 template class Cluster<3>;
 
 }  // namespace tamaas
diff --git a/src/percolation/flood_fill.hh b/src/percolation/flood_fill.hh
index c760b9d..a7a58c1 100644
--- a/src/percolation/flood_fill.hh
+++ b/src/percolation/flood_fill.hh
@@ -1,77 +1,85 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef FLOOD_FILL_H
 #define FLOOD_FILL_H
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 #include <list>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 class Cluster {
   using Point = std::array<Int, dim>;
+  using BBox = std::pair<std::array<Int, dim>, std::array<Int, dim>>;
 
 public:
   /// Constructor
   Cluster(Point start, const Grid<bool, dim>& map, Grid<bool, dim>& visited,
           bool diagonal);
   /// Copy constructor
   Cluster(const Cluster& other);
   /// Default constructor
   Cluster() = default;
   /// Get area of cluster
   UInt getArea() const { return getPoints().size(); }
   /// Get perimeter of cluster
   UInt getPerimeter() const { return perimeter; }
   /// Get contact points
-  const std::list<Point>& getPoints() const { return points; }
+  const auto& getPoints() const { return points; }
   /// Get bounding box
-  std::pair<std::array<Int, dim>, std::array<Int, dim>> boundingBox() const;
+  BBox boundingBox() const;
+  /// Get bounding box extent
+  std::array<Int, dim> extent() const;
   /// Assign next neighbors
-  auto getNextNeighbors(const std::array<UInt, dim>& p);
+  auto getNextNeighbors(const std::array<Int, dim>& p);
   /// Assign diagonal neighbors
-  auto getDiagonalNeighbors(const std::array<UInt, dim>& p);
+  auto getDiagonalNeighbors(const std::array<Int, dim>& p);
 
 private:
-  std::list<Point> points;
+  /// List of points in the cluster
+  std::vector<Point> points;
+  /// Perimeter size (number of segments)
   UInt perimeter = 0;
 };
 
 /* -------------------------------------------------------------------------- */
 
 class FloodFill {
+  template <UInt dim>
+  using List = std::vector<Cluster<dim>>;
+
 public:
-  static std::list<Cluster<1>> getSegments(const Grid<bool, 1>& map);
-  /// Returns a list of clusters from a boolean contact map
-  static std::list<Cluster<2>> getClusters(const Grid<bool, 2>& map,
-                                           bool diagonal);
-  static std::list<Cluster<3>> getVolumes(const Grid<bool, 3>& map,
-                                          bool diagonal);
+  /// Return a list of connected segments
+  static List<1> getSegments(const Grid<bool, 1>& map);
+  /// Return a list of connected areas
+  static List<2> getClusters(const Grid<bool, 2>& map, bool diagonal);
+  /// Return a list of connected volumes
+  static List<3> getVolumes(const Grid<bool, 3>& map, bool diagonal);
 };
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 #endif  // FLOOD_FILL_H
diff --git a/tests/test_flood_fill.py b/tests/test_flood_fill.py
index 37b2b60..2fe8547 100644
--- a/tests/test_flood_fill.py
+++ b/tests/test_flood_fill.py
@@ -1,70 +1,101 @@
 # -*- coding: utf-8 -*-
 #
 # Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
 # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 # Copyright (©) 2020-2022 Lucas Frérot
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Affero General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # This program 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 Affero General Public License for more details.
 #
 # You should have received a copy of the GNU Affero General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 from __future__ import print_function, division
 
 import tamaas as tm
 import numpy as np
 
 
 def test_flood_fill2d():
     field = np.zeros((18, 18), np.bool)
 
     # Cluster of permeter 18 area 13
     field[2:4, 3:6] = True
-    field[4  , 4:6] = True  # noqa
-    field[5  , 5  ] = True  # noqa
+    field[4, 4:6] = True  # noqa
+    field[5, 5] = True  # noqa
     field[3:5, 6:8] = True  # noqa
 
     # Cluster of perimeter 8 area 4
     field[14:16, 3:5] = True
 
     # Cluster of perimeter 20 area 11
-    field[9:11, 8:11 ] = True  # noqa
-    field[7:9 , 9    ] = True  # noqa
-    field[10  , 11:14] = True  # noqa
+    field[9:11, 8:11] = True  # noqa
+    field[7:9, 9] = True  # noqa
+    field[10, 11:14] = True  # noqa
 
     # Cluster of perimeter 18 area 9
     field[3:5, 14:16] = True
     field[5:10, 15] = True
 
+    # Cluster spanning periodic boundary
+    field[0:3, 12:14] = True
+    field[-2:, 12:14] = True
+
+    # Percolating cluster
+    field[12, :] = True
+
+    if False:
+        # Show contact map
+        import matplotlib.pyplot as plt
+        plt.imshow(field, interpolation='none', cmap='gray_r', origin='lower')
+        plt.xlim(0, 17)
+        plt.ylim(0, 17)
+        plt.show()
+
     clusters = tm.FloodFill.getClusters(field, False)
 
     # Dummy class for clusters
     class dummy:
-        def __init__(self, area, perimeter):
+
+        def __init__(self, area, perimeter, extent=None):
             self.area = area
             self.perimeter = perimeter
+            self.extent = extent
 
     # Solution
-    ref = [dummy(13, 18), dummy(9, 18), dummy(11, 20), dummy(4, 8)]
+    ref = [
+        dummy(10, 14, [5, 2]),   # cluster spanning pbc
+        dummy(13, 18),
+        dummy(9, 18),
+        dummy(11, 20),
+        dummy(18, 36, [1, 18]),  # percolating cluster
+        dummy(4, 8)
+    ]
+
+    assert len(ref) == len(clusters)
 
     for x, y in zip(clusters, ref):
         assert x.area == y.area
         assert x.perimeter == y.perimeter
 
+        if y.extent is not None:
+            assert x.extent == y.extent
+
+    print(clusters[0].extent, clusters[4].extent)
+
 
 def test_flood_fill3d():
     field = np.zeros((5, 5, 5), np.bool)
     field[2:4, 3, 1:3] = True
 
     clusters = tm.FloodFill.getVolumes(field, False)
 
     assert len(clusters) == 1
     assert clusters[0].area == 4