feat: voxel connectivity graph to labeled image (#123)

* feat: add color connectivty graph for 4 and 6 connectivities

* chore: whitespace

* feat: add 6 and 8 color connectivity

* feat: added colorization for 8-bit 8-connected

* feat: add support for 8-connected

* feat: add support for 26 connected vcg coloring

* fix: start numbering from 1, iterate over z

* test: check that 6-way color connectivity works

* fix: handle different array sizes

* fix: two recurring issues in 6 and 8 connected

* test: check 4 and 8 connectivity for 4 connected logic

* test: remove debugging statements

* fix: mixed up uint8 and uint32 8-connected

* fix: label incrementing for 6 and 26 connected

* test: a basic test for 8-connected uint8

* test: basic test for 8-way uint32

* test: add check for 6-connectivity

* test: 6-connected won't accept uint8

* fix: limit print out of DisjointSet to actual size

* fix: incorrect limitation

* fix: segmentation fault

* fix: more segmentation faults

* fix: segfault

* test: check 26-connected supports 6-connected

* test: check for simple 26-connectedness

* test: check connecting center

* perf: reduce memory usage for 26-connected

* fix: segmentation fault

* fix: segmentation fault

* ensure self-contact is handled properly

automated_test.py

@@ -1260,3 +1260,167 @@ def test_periodic_boundary_6():
   assert N == 2
 
 
+@pytest.mark.parametrize("dtype", [np.uint8, np.uint32])
+@pytest.mark.parametrize("connectivity", [4,8])
+def test_color_connectivity_graph_4(dtype, connectivity):
+  vcg = np.array([[]], dtype=dtype)
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+  assert cc_labels.size == 0
+
+  vcg = np.zeros([2,1], dtype=dtype)
+  vcg[:] = 0x0f
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+  ans = np.array([1,1], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+  vcg = np.zeros([2,1], dtype=dtype)
+  vcg[:] = 0b1100
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+  ans = np.array([[1],[2]], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+  vcg = np.zeros([1,2], dtype=dtype)
+  vcg[:] = 0b1100
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+  ans = np.array([1,1], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+  vcg[:] = 0b0011
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+
+  ans = np.array([1,2], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+  vcg = np.zeros([2,2], dtype=dtype)
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+  ans = np.array([[1, 3],[2, 4]], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+  vcg[:] = 0b1111
+  vcg[0,0] = 0b1110
+  vcg[1,0] = 0b1101
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+  ans = np.array([[1, 1],[1, 1]], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+  vcg[:] = 0b1111
+  vcg[0,0] = 0b1110
+  vcg[1,0] = 0b1101
+  vcg[1,1] = 0b0111
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+  ans = np.array([[1, 1],[2, 1]], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+  vcg[:] = 0b1111
+  vcg[0,0] = 0b1010
+  vcg[1,0] = 0b1101
+  vcg[0,1] = 0b0111
+  vcg[1,1] = 0b0111
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=connectivity)
+  ans = np.array([[1, 3],[2, 3]], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+def test_color_connectivity_graph_8_uint8():
+  vcg = np.zeros([2,2], dtype=np.uint8)
+  vcg[0,0] = 0b10000
+  vcg[1,1] = 0b10000000
+  vcg[1,0] = 0b100000
+  vcg[0,1] = 0b1000000
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=8)
+  ans = np.array([[1, 2],[2, 1]], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+def test_color_connectivity_graph_8_uint32():
+  vcg = np.zeros([2,2], dtype=np.uint32)
+  vcg[0,0] = 0b1000000
+  vcg[1,1] = 0b1000000000
+  vcg[1,0] = 0b10000000
+  vcg[0,1] = 0b100000000
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=8)
+  ans = np.array([[1, 2],[2, 1]], dtype=np.uint32)
+  assert np.all(cc_labels == ans)
+
+def test_color_connectivity_graph_6():
+
+  with pytest.raises(ValueError):
+    vcg = np.zeros([2,2,2], dtype=np.uint8, order="F")
+    cc3d.color_connectivity_graph(vcg, connectivity=6)
+
+
+  vcg = np.zeros([10,10,10], dtype=np.uint32, order="F")
+  vcg[:,:,:] = 0b11111111111111111111111111
+
+  vcg[:,:,4] = vcg[:,:,4] & 0b101111
+  vcg[:,:,5] = vcg[:,:,5] & 0b011111
+
+  vcg[:,4,:] = vcg[:,4,:] & 0b111011
+  vcg[:,5,:] = vcg[:,5,:] & 0b110111
+
+  vcg[4,:,:] = vcg[4,:,:] & 0b111110
+  vcg[5,:,:] = vcg[5,:,:] & 0b111101
+
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=6)
+
+  out = np.zeros(vcg.shape, dtype=np.uint32, order="F")
+  out[:5,:5,:5] = 1
+  out[5:,:5,:5] = 2
+  out[:5,5:,:5] = 3
+  out[5:,5:,:5] = 4
+  out[:5,:5,5:] = 5
+  out[5:,:5,5:] = 6
+  out[:5,5:,5:] = 7
+  out[5:,5:,5:] = 8
+
+  assert np.all(out == cc_labels)
+
+def test_color_connectivity_graph_26():
+
+  with pytest.raises(ValueError):
+    vcg = np.zeros([2,1,2], dtype=np.uint8, order="F")
+    cc3d.color_connectivity_graph(vcg, connectivity=26)
+
+  vcg = np.zeros([10,10,10], dtype=np.uint32, order="F")
+  vcg[:,:,:] = 0b11111111111111111111111111
+
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=26)
+  assert np.all(cc_labels == 1)
+
+  vcg[:,:,4] = vcg[:,:,4] & 0b101111
+  vcg[:,:,5] = vcg[:,:,5] & 0b011111
+
+  vcg[:,4,:] = vcg[:,4,:] & 0b111011
+  vcg[:,5,:] = vcg[:,5,:] & 0b110111
+
+  vcg[4,:,:] = vcg[4,:,:] & 0b111110
+  vcg[5,:,:] = vcg[5,:,:] & 0b111101
+
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=26)
+
+  out = np.zeros(vcg.shape, dtype=np.uint32, order="F")
+  out[:5,:5,:5] = 1
+  out[5:,:5,:5] = 2
+  out[:5,5:,:5] = 3
+  out[5:,5:,:5] = 4
+  out[:5,:5,5:] = 5
+  out[5:,:5,5:] = 6
+  out[:5,5:,5:] = 7
+  out[5:,5:,5:] = 8
+
+  assert np.all(out == cc_labels)
+
+  vcg[4,4,4] = 0b11111111111111111111000000
+  vcg[5,4,4] = 0b11111111111111111111000000
+  vcg[4,5,4] = 0b11111111111111111111000000
+  vcg[5,5,4] = 0b11111111111111111111000000
+
+  vcg[4,4,5] = 0b11111111111111111111000000
+  vcg[5,4,5] = 0b11111111111111111111000000
+  vcg[4,5,5] = 0b11111111111111111111000000
+  vcg[5,5,5] = 0b11111111111111111111000000
+
+  cc_labels = cc3d.color_connectivity_graph(vcg, connectivity=26)
+  assert np.all(cc_labels == 1)
+
+
+
+

cc3d.hpp

@@ -137,7 +137,8 @@ class DisjointSet {
   }
 
   void print() {
-    for (int i = 0; i < 15; i++) {
+    int size = std::min(static_cast<int>(length), 15);
+    for (int i = 0; i < size; i++) {
       printf("%d, ", ids[i]);
     }
     printf("\n");
@@ -1082,7 +1083,7 @@ OUT* connected_components2d_8(
       else {
         next_label++;
         out_labels[loc] = next_label;
-        equivalences.add(out_labels[loc]);        
+        equivalences.add(out_labels[loc]);      
       }
     }
   }

cc3d.pyx

@@ -93,6 +93,13 @@ cdef extern from "cc3d_graphs.hpp" namespace "cc3d":
     vector[cpp_pair[size_t, size_t]] all_runs,
     T* labels, size_t voxels
   ) except +
+  cdef OUT* color_connectivity_graph_N[T,OUT](
+    T* vcg,
+    int64_t sx, int64_t sy, int64_t sz,
+    int connectivity,
+    OUT* out_labels,
+    uint64_t& N
+  ) except +
 
 ctypedef fused UINT:
   uint8_t
@@ -781,6 +788,86 @@ def _statistics_helper(
 
   return output
 
+@cython.binding(True)
+def color_connectivity_graph(
+  vcg, 
+  connectivity = 26, 
+  return_N = False, 
+) -> np.ndarray:
+  """
+  Given a voxel connectivity graph following the same bit convention as 
+  cc3d.voxel_connectivity_graph (see docstring), assuming an undirected
+  graph (the format supports directed graphs, but this is not implemented
+  for the sake of efficiency), this function will return a uint32 image
+  that contains connected components labeled according to the boundaries 
+  described in the voxel connectivity graph (vcg).
+  """
+  cdef int dims = len(vcg.shape)
+  if dims not in (2,3):
+    raise DimensionError("Only 2D, and 3D arrays supported. Got: " + str(dims))
+
+  if dims == 2 and connectivity not in [4,8,6,26]:
+    raise ValueError(f"Only 4 and 8 connectivity is supported for 2D images. Got: {connectivity}")
+  elif dims != 2 and connectivity not in [6,26]:
+    raise ValueError(f"Only 6 and 26 connectivity are supported for 3D images. Got: {connectivity}")
+
+  dtype = vcg.dtype
+  if dtype not in [np.uint8, np.uint32]:
+    raise ValueError(f"Only uint8 and uint32 are supported. Got: {vcg.dtype}")
+
+  if vcg.size == 0:
+    return np.zeros([0] * dims, dtype=np.uint32, order="F")
+
+  while vcg.ndim < 3:
+    vcg = vcg[..., np.newaxis]
+
+  vcg = np.asfortranarray(vcg)
+
+  shape = vcg.shape
+
+  cdef int sx = shape[0]
+  cdef int sy = shape[1]
+  cdef int sz = shape[2]
+
+  if connectivity in [6, 26] and sz > 1 and dtype != np.uint32:
+    raise ValueError(f"Only uint32 is supported for 3d connectivites. Got: {vcg.dtype}")
+
+  cdef uint8_t[:,:,:] arr_memview8u
+  cdef uint32_t[:,:,:] arr_memview32u
+  cdef uint32_t[:,:,:] out_labels32
+
+  cdef int64_t voxels = <int64_t>sx * <int64_t>sy * <int64_t>sz
+  out_labels = np.zeros( (sx,sy,sz,), dtype=np.uint32, order='F' )
+  out_labels32 = out_labels
+
+  cdef size_t N = 0
+
+  if dtype == np.uint32:
+    arr_memview32u = vcg
+    color_connectivity_graph_N[uint32_t, uint32_t](
+      &arr_memview32u[0,0,0],
+      sx, sy, sz,
+      connectivity,
+      &out_labels32[0,0,0],
+      N
+    )
+  else:
+    arr_memview8u = vcg
+    color_connectivity_graph_N[uint8_t, uint32_t](
+      &arr_memview8u[0,0,0],
+      sx, sy, sz,
+      connectivity,
+      &out_labels32[0,0,0],
+      N
+    )
+
+  while out_labels.ndim > dims:
+    out_labels = out_labels[...,0]
+
+  if return_N:
+    return (out_labels, N)
+  return out_labels
+
 @cython.binding(True)
 def voxel_connectivity_graph(
   data:np.ndarray,