Repository: seung-lab/connected-components-3d (language: python)

Issue:
Periodic Boundary Conditions
> Hi,

> It means that D(x)=D(x+L), where L is the length of the box along x, so that in the following case:

> data = np.array([[3,3,0,3,3]])
> labels_out = cc3d.connected_components(data, connectivity=6)

> Thereby instead of array([[1,1,0,2,2]], dtype=unit16) it will output array([[1,1,0,1,1]], dtype=unit16) if periodic boundary condition is applied.
Maybe this [wiki](https://en.wikipedia.org/wiki/Periodic_boundary_conditions) would be helpful in explaining what I meant.


Discussion / hints:
A simple solution could be to post-process the labels at the boundaries after labelling.

Fix patch:
diff --git a/README.md b/README.md
index f5f6a8f..a68fb7f 100644
--- a/README.md
+++ b/README.md
@@ -69,6 +69,13 @@ labels_out = cc3d.connected_components(labels_in) # 26-connected
 connectivity = 6 # only 4,8 (2D) and 26, 18, and 6 (3D) are allowed
 labels_out = cc3d.connected_components(labels_in, connectivity=connectivity)
 
+# If you need the borders to wrap around (useful for simulations)
+# specify periodic_boundary=True, currently only supported for
+# 4 (2d) and 6 (3d) connectivities.
+labels_out = cc3d.connected_components(
+  labels_in, connectivity=connectivity, periodic_boundary=True
+)
+
 # If you need a particular dtype you can specify np.uint16, np.uint32, or np.uint64
 # You can go bigger, not smaller, than the default which is selected
 # to be the smallest that can be safely used. This can save you the copy
diff --git a/cc3d.hpp b/cc3d.hpp
index 06da001..c03b70f 100644
--- a/cc3d.hpp
+++ b/cc3d.hpp
@@ -49,10 +49,11 @@
 #include <cmath>
 #include <cstdio>
 #include <cstdint>
+#include <limits>
+#include <memory>
 #include <stdexcept>
 #include <unordered_set>
 #include <vector>
-#include <limits>
 
 namespace cc3d {
 
@@ -252,7 +253,7 @@ OUT* relabel(
   }
 
   OUT label;
-  OUT* renumber = new OUT[num_labels + 1]();
+  std::unique_ptr<OUT[]> renumber(new OUT[num_labels + 1]());
   OUT next_label = 1;
 
   for (int64_t i = 1; i <= num_labels; i++) {
@@ -279,7 +280,6 @@ OUT* relabel(
     }
   }
 
-  delete[] renumber;
   return out_labels;
 }
 
@@ -366,7 +366,9 @@ OUT* connected_components3d_26(
   
   DisjointSet<OUT> equivalences(max_labels);
 
-  const uint32_t *runs = compute_foreground_index(in_labels, sx, sy, sz);
+  const std::unique_ptr<uint32_t[]> runs(
+    compute_foreground_index(in_labels, sx, sy, sz)
+  );
      
   /*
     Layout of forward pass mask (which faces backwards). 
@@ -590,10 +592,7 @@ OUT* connected_components3d_26(
     }
   }
   
-
-  out_labels = relabel<OUT>(out_labels, sx, sy, sz, next_label, equivalences, N, runs);
-  delete[] runs;
-  return out_labels;
+  return relabel<OUT>(out_labels, sx, sy, sz, next_label, equivalences, N, runs.get());
 }
 
 template <typename T, typename OUT = uint32_t>
@@ -621,7 +620,9 @@ OUT* connected_components3d_18(
 
   DisjointSet<OUT> equivalences(max_labels);
 
-  const uint32_t *runs = compute_foreground_index(in_labels, sx, sy, sz);
+  const std::unique_ptr<uint32_t[]> runs(
+    compute_foreground_index(in_labels, sx, sy, sz)
+  );
      
   /*
     Layout of forward pass mask (which faces backwards). 
@@ -751,9 +752,7 @@ OUT* connected_components3d_18(
     }
   }
 
-  out_labels = relabel<OUT>(out_labels, sx, sy, sz, next_label, equivalences, N, runs);
-  delete[] runs;
-  return out_labels;
+  return relabel<OUT>(out_labels, sx, sy, sz, next_label, equivalences, N, runs.get());
 }
 
 template <typename T, typename OUT = uint32_t>
@@ -761,7 +760,8 @@ OUT* connected_components3d_6(
     T* in_labels, 
     const int64_t sx, const int64_t sy, const int64_t sz,
     size_t max_labels, 
-    OUT *out_labels = NULL, size_t &N = _dummy_N
+    OUT *out_labels = NULL, size_t &N = _dummy_N,
+    bool periodic_boundary = false
   ) {
 
   const int64_t sxy = sx * sy;
@@ -781,7 +781,9 @@ OUT* connected_components3d_6(
 
   DisjointSet<OUT> equivalences(max_labels);
 
-  const uint32_t *runs = compute_foreground_index(in_labels, sx, sy, sz);
+  const std::unique_ptr<uint32_t[]> runs(
+    compute_foreground_index(in_labels, sx, sy, sz)
+  );
 
   /*
     Layout of forward pass mask (which faces backwards). 
@@ -860,9 +862,34 @@ OUT* connected_components3d_6(
     }
   }
 
-  out_labels = relabel<OUT>(out_labels, sx, sy, sz, next_label, equivalences, N, runs);
-  delete[] runs;
-  return out_labels;
+  if (periodic_boundary) {
+    for (int64_t z = 0; z < sz; z++) {
+      for (int64_t y = 0; y < sy; y++) {
+        loc = sx * (y + sy * z);
+        if (in_labels[loc] != 0 && in_labels[loc] == in_labels[loc + sx - 1]) {
+          equivalences.unify(out_labels[loc], out_labels[loc + sx - 1]);
+        }
+      }
+    }
+    for (int64_t z = 0; z < sz; z++) {
+      for (int64_t x = 0; x < sx; x++) {
+        loc = x + sxy * z;
+        if (in_labels[loc] != 0 && in_labels[loc] == in_labels[loc + sx * (sy - 1)]) {
+          equivalences.unify(out_labels[loc], out_labels[loc + sx * (sy - 1)]);
+        }
+      }
+    }
+    for (int64_t y = 0; y < sy; y++) {
+      for (int64_t x = 0; x < sx; x++) {
+        loc = x + sx * y;
+        if (in_labels[loc] != 0 && in_labels[loc] == in_labels[loc + sxy * (sz - 1)]) {
+          equivalences.unify(out_labels[loc], out_labels[loc + sxy * (sz - 1)]);
+        }
+      }
+    }
+  }
+
+  return relabel<OUT>(out_labels, sx, sy, sz, next_label, equivalences, N, runs.get());
 }
 
 
@@ -875,7 +902,8 @@ OUT* connected_components2d_4(
     T* in_labels, 
     const int64_t sx, const int64_t sy, 
     size_t max_labels, 
-    OUT *out_labels = NULL, size_t &N = _dummy_N
+    OUT *out_labels = NULL, size_t &N = _dummy_N,
+    const bool periodic_boundary = false
   ) {
 
   const int64_t voxels = sx * sy;
@@ -894,7 +922,9 @@ OUT* connected_components2d_4(
 
   DisjointSet<OUT> equivalences(max_labels);
 
-  const uint32_t *runs = compute_foreground_index(in_labels, sx, sy, /*sz=*/1);
+  const std::unique_ptr<uint32_t[]> runs(
+    compute_foreground_index(in_labels, sx, sy, /*sz=*/1)
+  );
     
   /*
     Layout of forward pass mask. 
@@ -945,9 +975,21 @@ OUT* connected_components2d_4(
     }
   }
 
-  out_labels = relabel<OUT>(out_labels, sx, sy, /*sz=*/1, next_label, equivalences, N, runs);
-  delete[] runs;
-  return out_labels;
+  if (periodic_boundary) {
+    for (int64_t x = 0; x < sx; x++) {
+      if (in_labels[x] != 0 && in_labels[x] == in_labels[x + sx * (sy - 1)]) {
+        equivalences.unify(out_labels[x], out_labels[x + sx * (sy - 1)]);
+      }
+    }
+    for (int64_t y = 0; y < sy; y++) {
+      loc = sx * y;
+      if (in_labels[loc] != 0 && in_labels[loc] == in_labels[loc + (sx - 1)]) {
+        equivalences.unify(out_labels[loc], out_labels[loc + (sx - 1)]);
+      }
+    }
+  }
+
+  return relabel<OUT>(out_labels, sx, sy, /*sz=*/1, next_label, equivalences, N, runs.get());
 }
 
 // K. Wu, E. Otoo, K. Suzuki. "Two Strategies to Speed up Connected Component Labeling Algorithms". 
@@ -979,7 +1021,9 @@ OUT* connected_components2d_8(
   
   DisjointSet<OUT> equivalences(max_labels);
 
-  const uint32_t *runs = compute_foreground_index(in_labels, sx, sy, /*sz=*/1);
+  const std::unique_ptr<uint32_t[]> runs(
+    compute_foreground_index(in_labels, sx, sy, /*sz=*/1)
+  );
 
   /*
     Layout of mask. We start from e.
@@ -1042,9 +1086,7 @@ OUT* connected_components2d_8(
     }
   }
 
-  out_labels = relabel<OUT>(out_labels, sx, sy, /*sz=*/1, next_label, equivalences, N, runs);
-  delete[] runs;
-  return out_labels;
+  return relabel<OUT>(out_labels, sx, sy, /*sz=*/1, next_label, equivalences, N, runs.get());
 }
 
 template <typename T, typename OUT = uint32_t>
@@ -1052,7 +1094,8 @@ OUT* connected_components3d(
     T* in_labels, 
     const int64_t sx, const int64_t sy, const int64_t sz,
     size_t max_labels, const int64_t connectivity,
-    OUT *out_labels = NULL, size_t &N = _dummy_N
+    OUT *out_labels = NULL, size_t &N = _dummy_N, 
+    bool periodic_boundary = false
   ) {
 
   if (connectivity == 26) {
@@ -1070,7 +1113,7 @@ OUT* connected_components3d(
   else if (connectivity == 6) {
     return connected_components3d_6<T, OUT>(
       in_labels, sx, sy, sz, 
-      max_labels, out_labels, N
+      max_labels, out_labels, N, periodic_boundary
     );
   }
   else if (connectivity == 8) {
@@ -1088,7 +1131,7 @@ OUT* connected_components3d(
     }
     return connected_components2d_4<T, OUT>(
       in_labels, sx, sy, 
-      max_labels, out_labels, N
+      max_labels, out_labels, N, periodic_boundary
     );
   }
   else {
@@ -1100,11 +1143,16 @@ template <typename T, typename OUT = uint32_t>
 OUT* connected_components3d(
     T* in_labels, 
     const int64_t sx, const int64_t sy, const int64_t sz,
-    const int64_t connectivity=26, size_t &N = _dummy_N
+    const int64_t connectivity=26, size_t &N = _dummy_N,
+    const bool periodic_boundary = false
   ) {
   const size_t voxels = sx * sy * sz;
   size_t max_labels = std::min(estimate_provisional_label_count(in_labels, sx, voxels), voxels);
-  return connected_components3d<T, OUT>(in_labels, sx, sy, sz, max_labels, connectivity, NULL, N);
+  return connected_components3d<T, OUT>(
+    in_labels, sx, sy, sz, 
+    max_labels, connectivity, 
+    NULL, N, periodic_boundary
+  );
 }
 
 };
diff --git a/cc3d.pyx b/cc3d.pyx
index c5ef2a9..82cf9aa 100644
--- a/cc3d.pyx
+++ b/cc3d.pyx
@@ -68,7 +68,7 @@ cdef extern from "cc3d_continuous.hpp" namespace "cc3d":
     T* in_labels, 
     int64_t sx, int64_t sy, int64_t sz,
     int64_t max_labels, int64_t connectivity, T delta,
-    U* out_labels, size_t &N
+    U* out_labels, size_t &N, native_bool periodic_boundary
   ) except +
 
 cdef extern from "cc3d_graphs.hpp" namespace "cc3d":
@@ -231,10 +231,14 @@ def estimate_provisional_labels(data:np.ndarray) -> Tuple[int,int,int]:
 
 @cython.binding(True)
 def connected_components(
-  data:np.ndarray, int64_t max_labels=-1, 
-  int64_t connectivity=26, native_bool return_N=False,
-  delta:Union[int,float] = 0, out_dtype:Optional[Any] = None,
-  out_file:Optional[Union[str, BinaryIO]] = None
+  data:np.ndarray, 
+  int64_t max_labels=-1, 
+  int64_t connectivity=26, 
+  native_bool return_N=False,
+  delta:Union[int,float] = 0, 
+  out_dtype:Optional[Any] = None,
+  out_file:Optional[Union[str, BinaryIO]] = None,
+  periodic_boundary:bool = False,
 ) -> np.ndarray:
   """
   Connected components applied to 3D images with 
@@ -262,6 +266,7 @@ def connected_components(
       operation by outputting the appropriately sized type instead.
     out_file: If specified, the output array will be an mmapped
       file. Can be a file-name or a file-like object.
+    periodic_boundary: the boundary edges wrap around
 
   let OUT = 1D, 2D or 3D numpy array remapped to reflect
     the connected components sequentially numbered from 1 to N. 
@@ -285,6 +290,11 @@ def connected_components(
   elif dims != 2 and connectivity not in (6, 18, 26):
     raise ValueError("Only 6, 18, and 26 connectivities are supported for 3D images. Got: " + str(connectivity))
 
+  if periodic_boundary and connectivity not in (4, 6):
+    raise ValueError(f"periodic_boundary is not yet implemented for {connectivity}-connectivity.")
+  if periodic_boundary and delta != 0:
+    raise ValueError(f"periodic_boundary is not yet implemented continuous data.")
+
   if data.size == 0:
     out_labels = np.zeros(shape=(0,), dtype=data.dtype)
     if return_N:
@@ -416,19 +426,19 @@ def connected_components(
         connected_components3d[uint64_t, uint16_t](
           &arr_memview64u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint16_t*>&out_labels16[0], N
+          <uint16_t*>&out_labels16[0], N, periodic_boundary
         )
       elif out_dtype == np.uint32:
         connected_components3d[uint64_t, uint32_t](
           &arr_memview64u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint32_t*>&out_labels32[0], N
+          <uint32_t*>&out_labels32[0], N, periodic_boundary
         )
       elif out_dtype == np.uint64:
         connected_components3d[uint64_t, uint64_t](
           &arr_memview64u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint64_t*>&out_labels64[0], N
+          <uint64_t*>&out_labels64[0], N, periodic_boundary
         )
     elif dtype in (np.uint32, np.int32):
       arr_memview32u = data.view(np.uint32)
@@ -436,19 +446,19 @@ def connected_components(
         connected_components3d[uint32_t, uint16_t](
           &arr_memview32u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint16_t*>&out_labels16[0], N
+          <uint16_t*>&out_labels16[0], N, periodic_boundary
         )
       elif out_dtype == np.uint32:
         connected_components3d[uint32_t, uint32_t](
           &arr_memview32u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint32_t*>&out_labels32[0], N
+          <uint32_t*>&out_labels32[0], N, periodic_boundary
         )
       elif out_dtype == np.uint64:
         connected_components3d[uint32_t, uint64_t](
           &arr_memview32u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint64_t*>&out_labels64[0], N
+          <uint64_t*>&out_labels64[0], N, periodic_boundary
         )
     elif dtype in (np.uint16, np.int16):
       arr_memview16u = data.view(np.uint16)
@@ -456,19 +466,19 @@ def connected_components(
         connected_components3d[uint16_t, uint16_t](
           &arr_memview16u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint16_t*>&out_labels16[0], N
+          <uint16_t*>&out_labels16[0], N, periodic_boundary
         )
       elif out_dtype == np.uint32:
         connected_components3d[uint16_t, uint32_t](
           &arr_memview16u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint32_t*>&out_labels32[0], N
+          <uint32_t*>&out_labels32[0], N, periodic_boundary
         )
       elif out_dtype == np.uint64:
         connected_components3d[uint16_t, uint64_t](
           &arr_memview16u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint64_t*>&out_labels64[0], N
+          <uint64_t*>&out_labels64[0], N, periodic_boundary
         )
     elif dtype in (np.uint8, np.int8, bool):
       arr_memview8u = data.view(np.uint8)
@@ -476,19 +486,19 @@ def connected_components(
         connected_components3d[uint8_t, uint16_t](
           &arr_memview8u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint16_t*>&out_labels16[0], N
+          <uint16_t*>&out_labels16[0], N, periodic_boundary
         )
       elif out_dtype == np.uint32:
         connected_components3d[uint8_t, uint32_t](
           &arr_memview8u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint32_t*>&out_labels32[0], N
+          <uint32_t*>&out_labels32[0], N, periodic_boundary
         )
       elif out_dtype == np.uint64:
         connected_components3d[uint8_t, uint64_t](
           &arr_memview8u[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint64_t*>&out_labels64[0], N
+          <uint64_t*>&out_labels64[0], N, periodic_boundary
         )
     elif dtype == np.float32:
       arr_memviewf = data
@@ -496,19 +506,19 @@ def connected_components(
         connected_components3d[float, uint16_t](
           &arr_memviewf[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint16_t*>&out_labels16[0], N
+          <uint16_t*>&out_labels16[0], N, periodic_boundary
         )
       elif out_dtype == np.uint32:
         connected_components3d[float, uint32_t](
           &arr_memviewf[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint32_t*>&out_labels32[0], N
+          <uint32_t*>&out_labels32[0], N, periodic_boundary
         )
       elif out_dtype == np.uint64:
         connected_components3d[float, uint64_t](
           &arr_memviewf[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint64_t*>&out_labels64[0], N
+          <uint64_t*>&out_labels64[0], N, periodic_boundary
         )
     elif dtype == np.float64:
       arr_memviewd = data
@@ -516,19 +526,19 @@ def connected_components(
         connected_components3d[double, uint16_t](
           &arr_memviewd[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint16_t*>&out_labels16[0], N
+          <uint16_t*>&out_labels16[0], N, periodic_boundary
         )
       elif out_dtype == np.uint32:
         connected_components3d[double, uint32_t](
           &arr_memviewd[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint32_t*>&out_labels32[0], N
+          <uint32_t*>&out_labels32[0], N, periodic_boundary
         )
       elif out_dtype == np.uint64:
         connected_components3d[double, uint64_t](
           &arr_memviewd[0,0,0],
           sx, sy, sz, max_labels, connectivity, delta,
-          <uint64_t*>&out_labels64[0], N
+          <uint64_t*>&out_labels64[0], N, periodic_boundary
         )
     else:
       raise TypeError("Type {} not currently supported.".format(dtype))
diff --git a/cc3d_continuous.hpp b/cc3d_continuous.hpp
index 7c272ae..3dd047f 100644
--- a/cc3d_continuous.hpp
+++ b/cc3d_continuous.hpp
@@ -118,7 +118,9 @@ OUT* connected_components3d_continuous(
 
   DisjointSet<OUT> equivalences(max_labels);
 
-  const uint32_t *runs = compute_foreground_index(in_labels, sx, sy, sz);
+  const std::unique_ptr<uint32_t[]> runs(
+    compute_foreground_index(in_labels, sx, sy, sz)
+  );
 
 
   int64_t loc = 0;
@@ -183,9 +185,7 @@ OUT* connected_components3d_continuous(
     }
   }
 
-  out_labels = relabel<OUT>(out_labels, sx, sy, sz, next_label, equivalences, N, runs);
-  delete[] runs;
-  return out_labels;
+  return relabel<OUT>(out_labels, sx, sy, sz, next_label, equivalences, N, runs.get());
 }
 
 template <typename T, typename OUT = uint32_t>
@@ -212,7 +212,9 @@ OUT* connected_components2d_4(
 
   DisjointSet<OUT> equivalences(max_labels);
 
-  const uint32_t *runs = compute_foreground_index(in_labels, sx, sy, /*sz=*/1);
+  const std::unique_ptr<uint32_t[]> runs(
+    compute_foreground_index(in_labels, sx, sy, /*sz=*/1)
+  );
 
   /*
     Layout of forward pass mask. 
@@ -265,12 +267,10 @@ OUT* connected_components2d_4(
     }
   }
 
-  out_labels = relabel<OUT>(
+  return relabel<OUT>(
     out_labels, sx, sy, /*sz=*/1, next_label, 
-    equivalences, N, runs
+    equivalences, N, runs.get()
   );
-  delete[] runs;
-  return out_labels;
 }
 
 template <typename T, typename OUT = uint32_t>
@@ -297,7 +297,9 @@ OUT* connected_components2d_8(
   
   DisjointSet<OUT> equivalences(max_labels);
 
-  const uint32_t *runs = compute_foreground_index(in_labels, sx, sy, /*sz=*/1);
+  const std::unique_ptr<uint32_t[]> runs(
+    compute_foreground_index(in_labels, sx, sy, /*sz=*/1)
+  );
 
   T gmax = in_labels[0];
   T gmin = in_labels[0];
@@ -392,9 +394,7 @@ OUT* connected_components2d_8(
     }
   }
 
-  out_labels = relabel<OUT>(out_labels, sx, sy, /*sz=*/1, next_label, equivalences, N, runs);
-  delete[] runs;
-  return out_labels;
+  return relabel<OUT>(out_labels, sx, sy, /*sz=*/1, next_label, equivalences, N, runs.get());
 }
 
 template <typename T, typename OUT = uint32_t>
@@ -402,7 +402,8 @@ OUT* connected_components3d(
     T* in_labels, 
     const int64_t sx, const int64_t sy, const int64_t sz,
     size_t max_labels, const int64_t connectivity, const T delta,
-    OUT *out_labels = NULL, size_t &N = _dummy_N
+    OUT *out_labels = NULL, size_t &N = _dummy_N, 
+    const bool periodic_boundary = false
   ) {
 
   // for performance, shouldn't be "more correct"
@@ -410,10 +411,14 @@ OUT* connected_components3d(
     return connected_components3d<T,OUT>(
       in_labels, sx, sy, sz, 
       max_labels, connectivity, 
-      out_labels, N
+      out_labels, N, periodic_boundary
     );
   }
 
+  if (periodic_boundary) {
+    throw std::runtime_error("periodic_boundary is not currently supported for continuous data.");
+  }
+
   if (connectivity == 26 || connectivity == 18 || connectivity == 6) {
     return connected_components3d_continuous<T, OUT>(
       in_labels, sx, sy, sz, 


Test patch:
diff --git a/automated_test.py b/automated_test.py
index 642ead5..5c11613 100644
--- a/automated_test.py
+++ b/automated_test.py
@@ -1132,8 +1132,70 @@ def test_largest_k(k):
 
   assert np.all(lbls == retained_labels)
 
+def test_periodic_boundary_4():
+  labels = np.zeros((10,10), dtype=np.uint32, order="F")
 
+  labels[:,0] = 1
+  labels[:,-1] = 1
 
+  out, N = cc3d.connected_components(labels, connectivity=4, return_N=True, periodic_boundary=False)
+  assert N == 2
 
+  out, N = cc3d.connected_components(labels, connectivity=4, return_N=True, periodic_boundary=True)
+  assert N == 1
+
+  labels = np.zeros((10,10), dtype=np.uint32, order="F")
+
+  labels[0,:] = 1
+  labels[-1,:] = 1
+
+  out, N = cc3d.connected_components(labels, connectivity=4, return_N=True, periodic_boundary=False)
+  assert N == 2
+
+  out, N = cc3d.connected_components(labels, connectivity=4, return_N=True, periodic_boundary=True)
+  assert N == 1
+
+def test_periodic_boundary_6():
+  labels = np.zeros((10,10,10), dtype=np.uint32, order="F")
+
+  labels[:,:,0] = 1
+  labels[:,:,-1] = 1
+
+  out, N = cc3d.connected_components(labels, connectivity=6, return_N=True, periodic_boundary=False)
+  assert N == 2
+
+  out, N = cc3d.connected_components(labels, connectivity=6, return_N=True, periodic_boundary=True)
+  assert N == 1
+
+  labels[:,:,:] = 0
+  labels[:,0,:] = 1
+  labels[:,-1,:] = 1
+
+  out, N = cc3d.connected_components(labels, connectivity=6, return_N=True, periodic_boundary=False)
+  assert N == 2
+
+  out, N = cc3d.connected_components(labels, connectivity=6, return_N=True, periodic_boundary=True)
+  assert N == 1
+
+  labels[:,:,:] = 0
+  labels[0,:,:] = 1
+  labels[-1,:,:] = 1
+
+  out, N = cc3d.connected_components(labels, connectivity=6, return_N=True, periodic_boundary=False)
+  assert N == 2
+
+  out, N = cc3d.connected_components(labels, connectivity=6, return_N=True, periodic_boundary=True)
+  assert N == 1
+
+
+  labels[:,:,:] = 0
+  labels[1,:,:] = 1
+  labels[-2,:,:] = 1
+
+  out, N = cc3d.connected_components(labels, connectivity=6, return_N=True, periodic_boundary=False)
+  assert N == 2
+
+  out, N = cc3d.connected_components(labels, connectivity=6, return_N=True, periodic_boundary=True)
+  assert N == 2