Prevent zero underflow in volume

cellfinder/core/detect/filters/volume/structure_splitting.py

@@ -74,6 +74,7 @@ def ball_filter_imgs(
     good_tiles_mask = np.ones((1, 1, volume.shape[2]), dtype=np.bool_)
 
     plane_width, plane_height = volume.shape[:2]
+    current_z = ball_z_size // 2
 
     bf = BallFilter(
         plane_width,
@@ -86,9 +87,7 @@ def ball_filter_imgs(
         threshold_value=threshold_value,
         soma_centre_value=soma_centre_value,
     )
-    cell_detector = CellDetector(
-        plane_width, plane_height, start_z=ball_z_size // 2
-    )
+    cell_detector = CellDetector(plane_width, plane_height, start_z=current_z)
 
     # FIXME: hard coded type
     ball_filtered_volume = np.zeros(volume.shape, dtype=np.uint32)
@@ -98,7 +97,11 @@ def ball_filter_imgs(
         if bf.ready:
             bf.walk()
             middle_plane = bf.get_middle_plane()
-            ball_filtered_volume[:, :, z] = middle_plane[:]
+
+            # first valid middle plane is the current_z, not z
+            ball_filtered_volume[:, :, current_z] = middle_plane[:]
+            current_z += 1
+
             # DEBUG: TEST: transpose
             previous_plane = cell_detector.process(
                 middle_plane.copy(), previous_plane
@@ -134,7 +137,10 @@ def iterative_ball_filter(
         vol, cell_centres = ball_filter_imgs(
             vol, threshold_value, soma_centre_value
         )
-        vol -= 1
+
+        # vol is unsigned, so can't let zeros underflow to max value
+        vol[:, :, :] = np.where(vol != 0, vol - 1, 0)
+
         n_structures = len(cell_centres)
         ns.append(n_structures)
         centres.append(cell_centres)

tests/core/test_integration/test_detection_structure_splitting.py

@@ -8,9 +8,13 @@
 
 import os
 
+import numpy as np
 import pytest
 from brainglobe_utils.IO.image.load import read_with_dask
 
+from cellfinder.core.detect.filters.volume.structure_splitting import (
+    split_cells,
+)
 from cellfinder.core.main import main
 
 data_dir = os.path.join(
@@ -46,3 +50,34 @@ def test_structure_splitting(signal_array, background_array):
         voxel_sizes,
         n_free_cpus=0,
     )
+
+
+def test_underflow_issue_435():
+    # two cells centered at (9, 10, 10), (19, 10, 10) with radius 5
+    p1 = np.array([9, 10, 10])
+    p2 = np.array([19, 10, 10])
+    radius = 5
+
+    bright_voxels = np.zeros((30, 20, 20), dtype=np.bool_)
+
+    pos = np.empty((30, 20, 20, 3))
+    pos[:, :, :, 0] = np.arange(30).reshape((-1, 1, 1))
+    pos[:, :, :, 1] = np.arange(20).reshape((1, -1, 1))
+    pos[:, :, :, 2] = np.arange(20).reshape((1, 1, -1))
+
+    dist1 = pos - p1.reshape((1, 1, 1, 3))
+    dist1 = np.sqrt(np.sum(np.square(dist1), axis=3))
+    inside1 = dist1 <= radius
+    dist2 = pos - p2.reshape((1, 1, 1, 3))
+    dist2 = np.sqrt(np.sum(np.square(dist2), axis=3))
+    inside2 = dist2 <= radius
+
+    bright_voxels[np.logical_or(inside1, inside2)] = True
+    bright_indices = np.argwhere(bright_voxels)
+
+    centers = split_cells(bright_indices)
+
+    # for some reason, same with pytorch, it's shifted by 1. Probably rounding
+    expected = {(10, 11, 11), (20, 11, 11)}
+    got = set(map(tuple, centers.tolist()))
+    assert expected == got