feat(voxel_connectivity_graph): extracts directed voxel connectivity … (#48)

* feat(voxel_connectivity_graph): extracts directed voxel connectivity graph

Each voxel of the output is a bitfield that represents the allowed
transit directions.

* docs: added description of bitfields to python help

* refactor: move region_graph tests to bottom

* test: add simple tests for 4 and 8 graph connectivity

* test: skip big memory test for 32-bit windows

* fix: voxel size + fix compiler warnings

* test: adds test for 3d voxel connectivity graph

* fix: the masks dont make sense in row-major format so just drop it for now

* docs: update README with voxel_connectivity_graph

* docs: describe how to use the voxel graph in C++

* docs: correction to data width of black cube

Author: william-silversmith

.gitignore

@@ -140,4 +140,4 @@ saved_images
 save_images
 
 .DS_Store
-callgrind*
+callgrind*

README.md

@@ -50,7 +50,7 @@ import numpy as np
 labels_in = np.ones((512, 512, 512), dtype=np.int32)
 labels_out = cc3d.connected_components(labels_in) # 26-connected
 
-connectivity = 6 # only 26, 18, and 6 are allowed
+connectivity = 6 # only 4,8 (2D) and 26, 18, and 6 (3D) are allowed
 labels_out = cc3d.connected_components(labels_in, connectivity=connectivity)
 
 # You can adjust the bit width of the output to accomodate
@@ -73,7 +73,15 @@ for segid in range(1, N+1):
 
 # We also include a region adjacency graph function 
 # that returns a set of undirected edges.
-graph = cc3d.region_graph(labels_out, connectivity=connectivity) 
+edges = cc3d.region_graph(labels_out, connectivity=connectivity) 
+
+# You can also generate a voxel connectivty graph that encodes
+# which directions are passable from a given voxel as a bitfield.
+# This could also be seen as a method of eroding voxels fractionally
+# based on their label adjacencies.
+# See help(cc3d.voxel_connectivity_graph) for details.
+graph = cc3d.voxel_connectivity_graph(labels, connectivity=connectivity)
+
 ```
 
 If you know approximately how many labels you are going to generate, you can save some memory by specifying a number a safety factor above that range. The max label ID in your input labels must be less than `max_labels`.
@@ -106,11 +114,22 @@ uint16_t* cc_labels = cc3d::connected_components3d<int, uint16_t>(
   /*connectivity=*/18 // default is 26 connected
 );
 
+#include "cc3d_graphs.hpp"
+
 // edges is [ e11, e12, e21, e22, ... ]
 std::vector<uint64_t> edges = cc3d::extract_region_graph<uint64_t>(
   labels, /*sx=*/512, /*sy=*/512, /*sz=*/512, 
   /*connectivity=*/18 // default is 26 connected
 );
+
+// graph is a series of bitfields that describe inter-voxel
+// connectivity based on adjacent labels. See "cc3d_graphs.hpp"
+// for details on the bitfield. 
+uint32_t* graph = extract_voxel_connectivity_graph<T>(
+  labels, /*sx=*/512, /*sy=*/512, /*sz=*/512, 
+  /*connectivity=*/6 // default is 26 connected
+);
+
 ```
 
 ## 26-Connected CCL Algorithm

automated_test.py

@@ -1,4 +1,5 @@
 import pytest
+import sys
 
 import cc3d
 import numpy as np
@@ -469,6 +470,55 @@ def test_compare_scipy_6(sparse):
 
   assert np.all(cc3d_labels == scipy_labels)
 
+@pytest.mark.skipif("sys.maxsize <= 2**33")
+@pytest.mark.xfail(raises=MemoryError, reason="Some build tools don't have enough memory for this.")
+def test_sixty_four_bit():
+  input_labels = np.ones((1626,1626,1626), dtype=np.uint8)
+  cc3d.connected_components(input_labels, max_labels=3)  
+
+@pytest.mark.parametrize("size", (255,256))
+def test_stress_upper_bound_for_binary_6(size):
+  labels = np.zeros((size,size,size), dtype=np.bool)
+  for z in range(labels.shape[2]):
+    for y in range(labels.shape[1]):
+      off = (y + (z % 2)) % 2
+      labels[off::2,y,z] = True
+
+  out = cc3d.connected_components(labels, connectivity=6)
+  assert np.max(out) + 1 <= (256**3) // 2 + 1
+
+@pytest.mark.parametrize("size", (255,256))
+def test_stress_upper_bound_for_binary_8(size):
+  labels = np.zeros((size,size), dtype=np.bool)
+  labels[0::2,0::2] = True
+
+  out = cc3d.connected_components(labels, connectivity=8)
+  assert np.max(out) + 1 <= (256**2) // 4 + 1
+
+  for _ in range(10):
+    labels = np.random.randint(0,2, (256,256), dtype=np.bool)
+    out = cc3d.connected_components(labels, connectivity=8)
+    assert np.max(out) + 1 <= (256**2) // 4 + 1    
+
+@pytest.mark.parametrize("size", (255,256))
+def test_stress_upper_bound_for_binary_18(size):
+  labels = np.zeros((size,size,size), dtype=np.bool)
+  labels[::2,::2,::2] = True
+  labels[1::2,1::2,::2] = True
+
+  out = cc3d.connected_components(labels, connectivity=26)
+  assert np.max(out) + 1 <= (256**3) // 4 + 1
+
+  for _ in range(10):
+    labels = np.random.randint(0,2, (256,256,256), dtype=np.bool)
+    out = cc3d.connected_components(labels, connectivity=26)
+    assert np.max(out) + 1 <= (256**2) // 4 + 1
+
+@pytest.mark.parametrize("size", (255,256))
+def test_stress_upper_bound_for_binary_26(size):
+  labels = np.zeros((size,size,size), dtype=np.bool)
+  labels[::2,::2,::2] = True
+
 @pytest.mark.parametrize("connectivity", (8, 18, 26))
 @pytest.mark.parametrize("dtype", TEST_TYPES)
 @pytest.mark.parametrize("out_dtype", OUT_TYPES)
@@ -595,53 +645,89 @@ def test_region_graph_6():
     (1,2), (1,3), (1,4), (1,5), (1,6), (1,7)
   ])
 
-@pytest.mark.parametrize("size", (255,256))
-def test_stress_upper_bound_for_binary_6(size):
-  labels = np.zeros((size,size,size), dtype=np.bool)
-  for z in range(labels.shape[2]):
-    for y in range(labels.shape[1]):
-      off = (y + (z % 2)) % 2
-      labels[off::2,y,z] = True
-
-  out = cc3d.connected_components(labels, connectivity=6)
-  assert np.max(out) + 1 <= (256**3) // 2 + 1
-
-@pytest.mark.parametrize("size", (255,256))
-def test_stress_upper_bound_for_binary_8(size):
-  labels = np.zeros((size,size), dtype=np.bool)
-  labels[0::2,0::2] = True
-
-  out = cc3d.connected_components(labels, connectivity=8)
-  assert np.max(out) + 1 <= (256**2) // 4 + 1
-
-  for _ in range(10):
-    labels = np.random.randint(0,2, (256,256), dtype=np.bool)
-    out = cc3d.connected_components(labels, connectivity=8)
-    assert np.max(out) + 1 <= (256**2) // 4 + 1    
-
-@pytest.mark.parametrize("size", (255,256))
-def test_stress_upper_bound_for_binary_18(size):
-  labels = np.zeros((size,size,size), dtype=np.bool)
-  labels[::2,::2,::2] = True
-  labels[1::2,1::2,::2] = True
-
-  out = cc3d.connected_components(labels, connectivity=26)
-  assert np.max(out) + 1 <= (256**3) // 4 + 1
-
-  for _ in range(10):
-    labels = np.random.randint(0,2, (256,256,256), dtype=np.bool)
-    out = cc3d.connected_components(labels, connectivity=26)
-    assert np.max(out) + 1 <= (256**2) // 4 + 1
+def test_voxel_graph_2d():
+  labels = np.ones((3,3), dtype=np.uint8)
+  graph = cc3d.voxel_connectivity_graph(labels, connectivity=4)
+  assert graph.dtype == np.uint8
+  assert np.all(graph)
+
+  graph = cc3d.voxel_connectivity_graph(labels, connectivity=8)
+  assert graph.dtype == np.uint8
+  assert np.all(graph)
+
+  labels[1,1] = 0
+  graph = cc3d.voxel_connectivity_graph(labels, connectivity=4)
+  gt = np.array([
+    [0x0f,       0b00001011, 0x0f      ],
+    [0b00001110, 0x00,       0b00001101],
+    [0x0f,       0b00000111, 0x0f      ]
+  ], dtype=np.uint8)
+  assert np.all(gt.T == graph)
 
-@pytest.mark.parametrize("size", (255,256))
-def test_stress_upper_bound_for_binary_26(size):
-  labels = np.zeros((size,size,size), dtype=np.bool)
-  labels[::2,::2,::2] = True
+  graph = cc3d.voxel_connectivity_graph(labels, connectivity=8)
+  gt = np.array([
+    [0b11101111, 0b11111011, 0b11011111],
+    [0b11111110, 0x00,       0b11111101],
+    [0b10111111, 0b11110111, 0b01111111]
+  ], dtype=np.uint8)
+  assert np.all(gt.T == graph)
+
+def test_voxel_graph_3d():
+  labels = np.ones((3,3,3), dtype=np.uint8)
+  graph = cc3d.voxel_connectivity_graph(labels, connectivity=26)
+  assert graph.dtype == np.uint32
+  assert np.all(graph)
+
+  E = {
+    (-1,-1,-1): 0b01111111111111111111111111,
+    (+1,-1,-1): 0b10111111111111111111111111,
+    (-1,+1,-1): 0b11011111111111111111111111,
+    (+1,+1,-1): 0b11101111111111111111111111,
+    (-1,-1,+1): 0b11110111111111111111111111,
+    (+1,-1,+1): 0b11111011111111111111111111,
+    (-1,+1,+1): 0b11111101111111111111111111,
+    (+1,+1,+1): 0b11111110111111111111111111,
+    ( 0,-1,-1): 0b11111111011111111111111111,
+    ( 0,+1,-1): 0b11111111101111111111111111,
+    (-1, 0,-1): 0b11111111110111111111111111,
+    (+1, 0,-1): 0b11111111111011111111111111,
+    ( 0,-1,+1): 0b11111111111101111111111111,
+    ( 0,+1,+1): 0b11111111111110111111111111,
+    (-1, 0,+1): 0b11111111111111011111111111,
+    (+1, 0,+1): 0b11111111111111101111111111,
+    (-1,-1, 0): 0b11111111111111110111111111,
+    (+1,-1, 0): 0b11111111111111111011111111,
+    (-1,+1, 0): 0b11111111111111111101111111,
+    (+1,+1, 0): 0b11111111111111111110111111,
+    ( 0, 0,-1): 0b11111111111111111111011111,
+    ( 0, 0,+1): 0b11111111111111111111101111,
+    ( 0,-1, 0): 0b11111111111111111111110111,
+    ( 0,+1, 0): 0b11111111111111111111111011,
+    (-1, 0, 0): 0b11111111111111111111111101,
+    (+1, 0, 0): 0b11111111111111111111111110,
+  }
+
+
+  labels[1,1,1] = 0
+  graph = cc3d.voxel_connectivity_graph(labels, connectivity=26)
+  gt = np.array([
+    [
+      [ E[(1, 1,1)], E[(0, 1,1)], E[(-1, 1,1)] ],
+      [ E[(1, 0,1)], E[(0, 0,1)], E[(-1, 0,1)] ],
+      [ E[(1,-1,1)], E[(0,-1,1)], E[(-1,-1,1)] ],
+    ],
+    [
+      [ E[(1, 1,0)], E[(0,1,0)], E[(-1, 1,0)] ],
+      [ E[(1, 0,0)], 0x00,       E[(-1, 0,0)] ],
+      [ E[(1,-1,0)], E[(0,-1,0)], E[(-1,-1,0)] ]    
+    ],
+    [
+      [ E[(1,1,-1)], E[(0,1,-1)], E[(-1,1,-1)] ],
+      [ E[(1,0,-1)], E[(0,0,-1)], E[(-1,0,-1)] ],
+      [ E[(1,-1,-1)], E[(0,-1,-1)], E[(-1,-1,-1)] ]
+    ]
+  ], dtype=np.uint32)
+  
+  assert np.all(gt.T == graph)
 
-  out = cc3d.connected_components(labels, connectivity=26)
-  assert np.max(out) + 1 <= (256**3) // 8 + 1
 
-  for _ in range(10):
-    labels = np.random.randint(0,2, (256,256,256), dtype=np.bool)
-    out = cc3d.connected_components(labels, connectivity=26)
-    assert np.max(out) + 1 <= (256**2) // 8 + 1

benchmarks/README.md

@@ -136,7 +136,7 @@ Here there's a slight difference in the memory usage. SciPy uses about 850 MB wh
 
 <p style="font-style: italics;" align="center">
 <img height=384 src="https://raw.githubusercontent.com/seung-lab/connected-components-3d/master/benchmarks/cc3d_sparse_black.png" alt="Fig. 4: Different configurations run against a uint64 512x512x512 black cube using 26-connectivity. (black) SciPy 1.5.2 (blue) cc3d 1.14.0 (red) cc3d 1.14.0 in sparse mode." /><br>
-Fig. 4: Different configurations run against a uint32 512x512x512 black cube using 26-connectivity. (black) SciPy 1.5.2 (blue) cc3d 1.14.0 (red) cc3d 1.14.0 in sparse mode.
+Fig. 4: Different configurations run against a uint64 512x512x512 black cube using 26-connectivity. (black) SciPy 1.5.2 (blue) cc3d 1.14.0 (red) cc3d 1.14.0 in sparse mode.
 </p>   
 
 Sometimes empty data shows up in your pipeline. Sometimes a lot of it. How do your libraries handle it? At full speed? Slower? Faster than normal?