diff --git a/src/percolation/flood_fill.cpp b/src/percolation/flood_fill.cpp
index 69781ee..c4507b4 100644
--- a/src/percolation/flood_fill.cpp
+++ b/src/percolation/flood_fill.cpp
@@ -1,271 +1,270 @@
 /*
  *  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>
+#include <queue>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
-Int wrap_pbc(Int i, Int n) {
-  if (i < 0)
-    return i + n;
-  if (i >= n)
-    return i - n;
-  return i;
-}
+// Both i and n should be signed, otherwise weird things happen
+Int wrap_pbc(Int i, Int n) { return ((i % n) + n) % n; }
 
 template <std::size_t dim>
 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] = 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] = 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<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<Int, 1>& /*p*/) {
   return std::vector<Point>();
 }
 
 template <>
 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<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<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<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>
 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 queue: a queue here prioritizes visiting local neighbors
+  // instead of growing a stack with each new visited point it gets rid of older
+  // neighbors first this prevents the points positions to be shifted
+  // indefinitely across periodic boundaries
+  std::queue<Point> visiting;
   visiting.push(start);
+
   // Contact sizes
   const auto& n = cast_int(map.sizes());
 
   while (not visiting.empty()) {
-    auto p = visiting.top();
+    auto p = visiting.front();
     auto p_wrapped = wrap_pbc(p, n);
 
     if (visited(p_wrapped)) {
       visiting.pop();
       continue;
     }
 
     visited(p_wrapped) = true;
-    points.push_back(visiting.top());
+    points.push_back(visiting.front());
     visiting.pop();
 
     auto process = [&](const std::vector<Point>& neighbors,
                        const bool do_perimeter) {
       for (const auto& p : neighbors) {
         const auto index = wrap_pbc(p, n);
         if (not visited(index) and map(index)) {
           visiting.push(p);
         }
 
         perimeter += do_perimeter and (not map(index));
       }
     };
 
     process(getNextNeighbors(p), true);
 
     if (diagonal) {
       process(getDiagonalNeighbors(p), false);
     }
   }
 }
 
 typename FloodFill::List<1> FloodFill::getSegments(const Grid<bool, 1>& map) {
   auto n = map.sizes();
   Grid<bool, 1> visited(n, 1);
   visited = false;
   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);
     }
   }
 
   return clusters;
 }
 
 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;
   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;
 }
 
 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;
   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);
         }
       }
     }
   }
 
   return clusters;
 }
 
 template class Cluster<1>;
 template class Cluster<2>;
 template class Cluster<3>;
 
 }  // namespace tamaas
diff --git a/tests/test_flood_fill.py b/tests/test_flood_fill.py
index 2fe8547..7173033 100644
--- a/tests/test_flood_fill.py
+++ b/tests/test_flood_fill.py
@@ -1,101 +1,104 @@
 # -*- 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[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
 
     # 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, extent=None):
             self.area = area
             self.perimeter = perimeter
             self.extent = extent
 
     # Solution
     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)
+    field[:] = True
+    field[2:4, 2:4] = False
+    clusters = tm.FloodFill.getClusters(field, False)
+    assert clusters[0].extent == [18, 18]
 
 
 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