feat(fill_holes_v2): true multi-label fill holes

.github/workflows/build_wheels.yml

@@ -28,7 +28,7 @@ jobs:
         uses: docker/setup-qemu-action@v1
 
       - name: Build wheels
-        uses: pypa/cibuildwheel@v3.1.4
+        uses: pypa/cibuildwheel@v3.2.0
         # to supply options, put them in 'env', like:
         env:
           CIBW_ARCHS_LINUX: ${{matrix.arch}}

.github/workflows/run_tests.yml

@@ -32,7 +32,7 @@ jobs:
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
-          python -m pip install pytest pybind11 setuptools wheel crackle-codec 
+          python -m pip install pytest pybind11 setuptools wheel crackle-codec scipy
           python -m pip install connected-components-3d edt fastremap numpy tqdm fill-voids
 
       - name: Compile

README.md

@@ -64,17 +64,53 @@ morphed = fastmorph.spherical_close(labels, radius=1, parallel=2, anisotropy=(1,
 # The rest support multilabel images.
 morphed = fastmorph.spherical_erode(labels, radius=1, parallel=2, anisotropy=(1,1,1))
 
+
+# Rapid multilabel hole filling. There are two versions that use different techniques
+# and have different interfaces for their "aggressive" modes. Both modes fill
+# holes appropriately by default.
+#
+# Generally speaking, fill_holes_v2 will be much faster. v2 uses a 
+# mostly linear time contact graph analysis. v1 analyzes a sequence 
+# of binary images. v2 exhibits much better scaling behavior and supports 
+# returning the filled and hole labels as  CrackleArray compressed objects 
+# to save memory.
+# 
+# The main advantage of v1 is that it includes a morphological closure mode
+# that operates on voxels for closing small holes. The downside is that this
+# can modify the surface of the object.
+#
+# v2 allows merging holes that are less than 100% closed, but if this
+# threshold is set too high, holes won't be closed. If it is too low,
+# improper merging can occur. 
+# 
+# In both methods, objects that contact the sides or more than one side
+# (in the case of fix_borders) cannot be merged.
+
+filled_labels, hole_labels = fastmorph.fill_holes_v2(labels)
+# requires: pip install crackle-codec
+# returns as compressed CrackleArrays that have speedy access to labels 
+# in the compressed state (often hundreds of times smaller than the full array)
+filled_labels, hole_labels = fastmorph.fill_holes_v2(labels, return_crackle=True)
+
+# fix_borders runs hole filling for each object on the edge to reduce edge contacts
+filled_labels, hole_labels = fastmorph.fill_holes_v2(labels, fix_borders=True)
+
+# merge_threshold (range 0.0 - 1.0) controls how much surface area can be 
+# "exposed" for a hole to still be filled. The default (1.0) means a hole
+# must be perfectly sealed (typical for hole filling algorithms).
+filled_labels, hole_labels = fastmorph.fill_holes_v2(labels, merge_threshold=0.97)
+
 # Note: for boolean images, this function will directly call fill_voids
 # and return a scalar for ct 
 # For integer images, more processing will be done to deal with multiple labels.
 # A dict of { label: num_voxels_filled } for integer images will be returned.
 # Note that for multilabel images, by default, if a label is totally enclosed by another,
 # a FillError will be raised. If remove_enclosed is True, the label will be overwritten.
-filled_labels, ct = fastmorph.fill_holes(labels, return_fill_count=True, remove_enclosed=False)
+filled_labels, ct = fastmorph.fill_holes_v1(labels, return_fill_count=True, remove_enclosed=False)
 
 # If the holes in your segmentation are imperfectly sealed, consider
 # using the following options.
-filled_labels = fastmorph.fill_holes(
+filled_labels = fastmorph.fill_holes_v1(
 	labels, 
 	# runs 2d fill on the sides of the cube for each binary image
 	fix_borders=True, 

automated_test.py

@@ -1,6 +1,9 @@
+import time
+
 import pytest
 import numpy as np
 import fastmorph
+import scipy.ndimage
 
 def test_spherical_dilate():
 	labels = np.zeros((10,10,10), dtype=bool)
@@ -109,6 +112,88 @@ def test_complex_fill():
 
 	assert np.all(res == ans)
 
+def test_fill_holes_v2():
+	labels = np.zeros((10,10,10), dtype=np.uint64)
+	filled, holes = fastmorph.fill_holes_v2(labels)
+	assert not np.any(filled)
+	assert not np.any(holes)
+
+	labels = np.ones((10,10,10), dtype=np.uint8)
+	labels[:,:,:5] = 2
+
+	labels[5,5,2] = 0
+	labels[5,5,7] = 0
+
+	assert np.count_nonzero(labels) == 998
+	filled, holes = fastmorph.fill_holes_v2(labels)
+	assert np.count_nonzero(filled) == 1000
+	assert list(np.unique(filled)) == [1,2]
+
+	assert filled[5,5,2] == 2
+	assert filled[5,5,7] == 1
+
+	labels = np.ones((10,10,10), dtype=np.uint32)
+	labels[5,5,2] = 777
+
+	filled, holes = fastmorph.fill_holes_v2(labels, return_crackle=True)
+	assert set(holes.labels()) == set([0,777])
+
+	labels = np.ones((10,10,10), dtype=bool)
+	labels[5,5,2] = 0
+
+	res, holes = fastmorph.fill_holes_v2(labels)
+	assert np.all(res)
+	assert not np.any(holes)
+
+def test_fill_v2_fix_borders():
+	labels = np.ones([100,100,100], dtype=np.uint8)
+	labels[40:60,40:60,:] = 2
+	labels[40:60,:,40:60] = 2
+	labels[:,40:60,40:60] = 2
+
+	filled, holes = fastmorph.fill_holes_v2(
+		labels, 
+		fix_borders=False,
+	)
+	assert np.all(filled == labels)
+
+	filled, holes = fastmorph.fill_holes_v2(
+		labels, 
+		fix_borders=True,
+	)
+	assert np.all(filled == 1)
+	assert np.count_nonzero(holes) == 20*20*100*3 - 20*20*20*2
+
+	labels = np.ones([100,100,100], dtype=np.uint8)
+	labels[40:60,40:60,0] = 2
+	labels[40:60,0,40:60] = 2
+	labels[0,40:60,40:60] = 2
+
+	filled, holes = fastmorph.fill_holes_v2(
+		labels, 
+		fix_borders=True,
+	)
+
+	assert np.all(filled == 1)
+	assert np.count_nonzero(holes) == 20*20*3
+
+def test_fill_holes_v2_data():
+	import crackle
+	labels = crackle.load("connectomics.npy.ckl.gz")
+
+	s = time.time()
+	filled, holes = fastmorph.fill_holes_v2(labels)
+	print(f"{time.time() - s:.3f}")
+
+	uniq = np.unique(filled)
+	for lbl in uniq[:10]:
+		print(lbl)
+		if lbl == 0:
+			continue
+		binary_image_fm = filled == lbl
+		binary_image_scipy = scipy.ndimage.binary_fill_holes(labels == lbl)
+		assert np.all(binary_image_fm == binary_image_scipy)
+
 def test_spherical_open_close_run():
 	labels = np.zeros((10,10,10), dtype=bool)
 	res = fastmorph.spherical_open(labels, radius=1)