feat: Continuous value CCL for direct application to microscopy images (#80)

* wip: experimenting with continuous value CCL on 4 & 8 connected

* feat: add in wip implementation of 6,18,26 connected

* fix: missing allocating provisional labels

* fix: remove debug statement

* chore: add new file to MANIFEST.in

* test: demonstrate that continuous values minimally works

* feat: add support for float and double input types

* fix: compilation error on MSVC

* test: add tests for floats and doubles

* docs: show how to use delta mode

* docs: show how to use delta

* docs: make delta visible in function sig

* test: add a randomized test to see if more voxels are joined

Author: william-silversmith

MANIFEST.in

@@ -1,3 +1,4 @@
 include cc3d.hpp
 include cc3d_graphs.hpp
+include cc3d_continuous.hpp
 include LICENSE

README.md

@@ -3,7 +3,7 @@
 cc3d: Connected Components on Multilabel 3D Images
 =======================
 
-Implementation of connected components in three dimensions using a 26, 18, or 6 connected neighborhood in 3D or 4 and 8-connected in 2D. This package uses a 3D variant of the two pass method by Rosenfeld and Pflatz augmented with Union-Find and a decision tree based on the 2D 8-connected work of Wu, Otoo, and Suzuki. This implementation is compatible with images containing many different labels, not just binary images. It can be used with 2D or 3D images. 
+Implementation of connected components in three dimensions using a 26, 18, or 6 connected neighborhood in 3D or 4 and 8-connected in 2D. This package uses a 3D variant of the two pass method by Rosenfeld and Pflatz augmented with Union-Find and a decision tree based on the 2D 8-connected work of Wu, Otoo, and Suzuki. This implementation is compatible with images containing many different labels, not just binary images. It also supports continuously valued images such as grayscale microscope images with an algorithm that joins together nearby values. It can be used with 2D or 3D images. 
 
 I wrote this package because I was working on densely labeled 3D biomedical images of brain tissue (e.g. 512x512x512 voxels). Other off the shelf implementations I reviewed were limited to binary images. This rendered these other packages too slow for my use case as it required masking each label and running the connected components algorithm once each time. For reference, there are often between hundreds to thousands of labels in a given volume. The benefit of this package is that it labels all connected components in one shot, improving performance by one or more orders of magnitude. 
 
@@ -51,6 +51,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're working with continuously valued images like microscopy
+# images you can use cc3d to perform a very rough segmentation. 
+# If delta = 0, standard high speed processing. If delta > 0, then
+# neighbor voxel values <= delta are considered the same component.
+# The algorithm can be 2-10x slower though.
+labels_out = cc3d.connected_components(labels_in, delta=10)
+
 # You can extract the number of labels (which is also the maximum 
 # label value) like so:
 labels_out, N = cc3d.connected_components(labels_in, return_N=True) # free

automated_test.py

@@ -7,6 +7,7 @@
 TEST_TYPES = [
   np.int8, np.int16, np.int32, np.int64,
   np.uint8, np.uint16, np.uint32, np.uint64,
+  np.float32, np.float64
 ]
 
 OUT_TYPES = [ np.uint16, np.uint32, np.uint64 ]
@@ -848,3 +849,69 @@ def test_statistics(order):
     "centroids": None 
   }
 
+@pytest.mark.parametrize("connectivity", (8, 18, 26))
+@pytest.mark.parametrize("dtype", TEST_TYPES)
+@pytest.mark.parametrize("order", ("C", "F"))
+def test_continuous_ccl_diagonal(order, dtype, connectivity):
+  labels = np.zeros((2,2), dtype=dtype, order=order)
+  labels[0,0] = 1
+  labels[1,0] = 2
+  labels[0,1] = 3
+  labels[1,1] = 4
+
+  out = cc3d.connected_components(labels, delta=0, connectivity=connectivity)
+  assert np.all(np.unique(labels) == [1,2,3,4])
+
+  out = cc3d.connected_components(labels, delta=1, connectivity=connectivity)
+  assert np.all(out == 1)
+
+@pytest.mark.parametrize("connectivity", (4, 6))
+@pytest.mark.parametrize("dtype", TEST_TYPES)
+@pytest.mark.parametrize("order", ("C", "F"))
+def test_continuous_ccl_4_6(order, dtype, connectivity):
+  labels = np.zeros((2,2), dtype=dtype, order=order)
+  labels[0,0] = 1
+  labels[1,0] = 2
+  labels[0,1] = 3
+  labels[1,1] = 4
+
+  out = cc3d.connected_components(labels, delta=0, connectivity=connectivity)
+  assert np.all(np.unique(labels) == [1,2,3,4])
+
+  out = cc3d.connected_components(labels, delta=1, connectivity=connectivity)
+  assert np.all(out == np.array([
+    [1, 2],
+    [1, 2],
+  ]))
+
+@pytest.mark.parametrize("dtype", TEST_TYPES)
+@pytest.mark.parametrize("connectivity", (4, 8))
+@pytest.mark.parametrize("order", ("C", "F"))
+def test_continuous_blocks(dtype, connectivity, order):
+  mask = np.random.randint(0,5, size=(64,64)).astype(dtype)
+  img = np.zeros((512,512), dtype=dtype)
+  img[64:128, 64:128] = 50 + mask
+  img[200:264, 64:128] = 70 + mask
+
+  img = np.ascontiguousarray(img)
+  if order == "F":
+    img = np.asfortranarray(img)
+
+  out = cc3d.connected_components(
+    img, connectivity=connectivity, delta=0
+  )
+  assert np.unique(out).size > 1000
+
+  out = cc3d.connected_components(
+    img, connectivity=connectivity, delta=1
+  )
+  assert np.unique(out).size > 3
+
+  out = cc3d.connected_components(
+    img, connectivity=connectivity, delta=5
+  )
+  assert np.all(np.unique(out)== [0,1,2])
+
+
+
+