feat: aggressive fill holes (#9)

* fix: closing and opening are idempotent, no need to provide iterations parameters

* feat(fill_holes): add fix_borders and morphological_close params

* fix: rename to morphological_closing and add pixel change counts

* fix: label may not be present in tiny cutouts after erosion

* test: check fix borders and morphological closing work

* ci: add crackle-codec to test run

* docs: show how to use new options

Author: william-silversmith

.github/workflows/run_tests.yml

@@ -28,7 +28,7 @@ jobs:
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
-          python -m pip install pytest pybind11 setuptools wheel
+          python -m pip install pytest pybind11 setuptools wheel crackle-codec
 
       - name: Compile
         run: python setup.py develop

README.md

@@ -10,7 +10,7 @@ We provide the following multithreaded (except where noted) operations:
 - Grayscale Stenciled Dilation, Erosion, Opening, Closing
 - Multi-Label Spherical Erosion
 - Binary Spherical Dilation, Opening, and Closing
-- Multi-Label Fill Voids (single threaded)
+- Multi-Label Fill Voids (mostly single threaded)
 
 Highlights compared to other libraries:
 
@@ -71,6 +71,17 @@ morphed = fastmorph.spherical_erode(labels, radius=1, parallel=2, anisotropy=(1,
 # 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)
+
+# If the holes in your segmentation are imperfectly sealed, consider
+# using the following options.
+filled_labels = fastmorph.fill_holes(
+	labels, 
+	# runs 2d fill on the sides of the cube for each binary image
+	fix_borders=True, 
+	# does a dilate and then an erode after filling holes
+	morphological_closing=True,
+)
+
 ```
 
 ## Performance

automated_test.py

@@ -69,6 +69,46 @@ def test_fill_holes():
 	assert res[1] == 1
 	assert res[2] == set()
 
+def test_fill_fix_borders():
+	labels = np.ones([100,100,100], dtype=np.uint8)
+	labels[40:60,40:60,:] = 0
+	labels[40:60,:,40:60] = 0
+	labels[:,40:60,40:60] = 0
+
+	res, cts = fastmorph.fill_holes(
+		labels, 
+		fix_borders=False,
+		return_fill_count=True,
+		morphological_closing=False,
+	)
+
+	assert cts[1] == 0
+	assert np.all(res == labels)
+
+	res, cts = fastmorph.fill_holes(
+		labels, 
+		fix_borders=True,
+		return_fill_count=True,
+		morphological_closing=False,
+	)
+
+	assert cts[1] == np.count_nonzero(labels == 0)
+	assert np.all(res == 1)
+
+def test_complex_fill():
+	import crackle
+	labels = crackle.load("test_data/complex_soma.ckl.gz")
+	ans = crackle.load("test_data/complex_soma_result.ckl.gz")
+
+	res = fastmorph.fill_holes(
+		labels, 
+		remove_enclosed=True,
+		fix_borders=True,
+		morphological_closing=True,
+	)
+
+	assert np.all(res == ans)
+
 def test_spherical_open_close_run():
 	labels = np.zeros((10,10,10), dtype=bool)
 	res = fastmorph.spherical_open(labels, radius=1)

fastmorph/__init__.py

@@ -121,74 +121,56 @@ def opening(
   background_only:bool = True,
   parallel:int = 0,
   mode:Mode = Mode.multilabel,
-  iterations:int = 1,
   erode_border:bool = True,
 ) -> np.ndarray:
   """Performs morphological opening of labels.
 
+  This operator is idempotent with the exception
+  of boundary effects.
+
   background_only is passed through to dilate.
     True: Only evaluate background voxels for dilation.
     False: Allow labels to erode each other as they grow.
   parallel: how many pthreads to use in a threadpool
   mode: 
     Mode.multilabel: are all surrounding pixels the same?
     Mode.grey: use grayscale image dilation (min value)
-
-  iterations: number of times to iterate the result
-
   erode_border: if True, the border is treated as background,
     else it is regarded as a value that would preserve the
     current value. Only has an effect for multilabel erosion.
   """
-  if iterations < 0:
-    raise ValueError(f"iterations ({iterations}) must be a positive integer.")
-  elif iterations == 0:
-    return np.copy(labels, order="F")
-
-  output = labels
-  for i in range(iterations):
-    output = dilate(
-      erode(output, parallel, mode, iterations=1, erode_border=erode_border),
-      background_only, parallel, mode, iterations=1
-    )
-  return output
+  return dilate(
+    erode(labels, parallel, mode, iterations=1, erode_border=erode_border),
+    background_only, parallel, mode, iterations=1
+  )
 
 def closing(
   labels:np.ndarray, 
   background_only:bool = True,
   parallel:int = 0,
   mode:Mode = Mode.multilabel,
-  iterations:int = 1,
   erode_border:bool = True,
 ) -> np.ndarray:
   """Performs morphological closing of labels.
 
+  This operator is idempotent with the exception
+  of boundary effects.
+
   background_only is passed through to dilate.
     True: Only evaluate background voxels for dilation.
     False: Allow labels to erode each other as they grow.
   parallel: how many pthreads to use in a threadpool
   mode: 
     Mode.multilabel: are all surrounding pixels the same?
     Mode.grey: use grayscale image dilation (min value)
-
-  iterations: number of times to iterate the result
-
   erode_border: if True, the border is treated as background,
     else it is regarded as a value that would preserve the
     current value. Only has an effect for multilabel erosion.
   """
-  if iterations < 0:
-    raise ValueError(f"iterations ({iterations}) must be a positive integer.")
-  elif iterations == 0:
-    return np.copy(labels, order="F")
-
-  output = labels
-  for i in range(iterations):
-    output = erode(
-      dilate(output, background_only, parallel, mode, iterations=1), 
-      parallel, mode, iterations=1, erode_border=erode_border,
-    )
-  return output
+  return erode(
+    dilate(labels, background_only, parallel, mode, iterations=1), 
+    parallel, mode, iterations=1, erode_border=erode_border,
+  )
 
 def spherical_dilate(
   labels:np.ndarray, 
@@ -283,6 +265,8 @@ def fill_holes(
   return_fill_count:bool = False,
   remove_enclosed:bool = False,
   return_removed:bool = False,
+  fix_borders:bool = False,
+  morphological_closing:bool = False,
 ) -> np.ndarray:
   """
   For fill holes in toplogically closed objects.
@@ -298,7 +282,7 @@ def fill_holes(
   Return value: (filled_labels, fill_count (if specified), removed_set (if specified))
   """
   assert np.issubdtype(labels.dtype, np.integer) or np.issubdtype(labels.dtype, bool), "fill_holes is currently only supported for integer or binary images."
-  if np.issubdtype(labels.dtype, bool):
+  if np.issubdtype(labels.dtype, bool) and not fix_borders and not morphological_closing:
     filled_labels, filled_ct = fill_voids.fill(labels, return_fill_count=True)
     ret = [ filled_labels ]
     if return_fill_count:
@@ -327,10 +311,36 @@ def fill_holes(
       continue
 
     binary_image = (cc_labels[slices] == label)
-    binary_image, pixels_filled = fill_voids.fill(
+
+    pixels_filled = 0
+
+    if morphological_closing:
+      dilated_binary_image = dilate(binary_image)
+      pixels_filled += np.sum(dilated_binary_image != binary_image)
+      binary_image = dilated_binary_image
+      del dilated_binary_image
+
+    if fix_borders:
+      binary_image[:,:,0], pf1 = fill_voids.fill(binary_image[:,:,0], return_fill_count=True)
+      binary_image[:,:,-1], pf2 = fill_voids.fill(binary_image[:,:,-1], return_fill_count=True)
+      binary_image[:,0,:], pf3 = fill_voids.fill(binary_image[:,0,:], return_fill_count=True)
+      binary_image[:,-1,:], pf4 = fill_voids.fill(binary_image[:,-1,:], return_fill_count=True)
+      binary_image[0,:,:], pf5 = fill_voids.fill(binary_image[0,:,:], return_fill_count=True)
+      binary_image[-1,:,:], pf6 = fill_voids.fill(binary_image[-1,:,:], return_fill_count=True)
+      pixels_filled += pf1 + pf2 + pf3 + pf4 + pf5 + pf6
+
+    binary_image, pf7 = fill_voids.fill(
       binary_image, in_place=True, 
       return_fill_count=True
     )
+    pixels_filled += pf7
+
+    if morphological_closing:
+      eroded_binary_image = erode(binary_image, erode_border=False)
+      pixels_filled -= np.sum(eroded_binary_image != binary_image)
+      binary_image = eroded_binary_image
+      del eroded_binary_image
+
     fill_counts[label] = pixels_filled
     output[slices][binary_image] = mapping[label]
 
@@ -339,7 +349,7 @@ def fill_holes(
 
     sub_labels = fastremap.unique(cc_labels[slices][binary_image])
     sub_labels = set(sub_labels)
-    sub_labels.remove(label)
+    sub_labels.discard(label)
     sub_labels.discard(0)
     if not remove_enclosed and sub_labels:
       sub_labels = [ int(l) for l in sub_labels ]