diff --git a/cellfinder/core/detect/filters/volume/structure_splitting.py b/cellfinder/core/detect/filters/volume/structure_splitting.py index f0018df8..7480e8ec 100644 --- a/cellfinder/core/detect/filters/volume/structure_splitting.py +++ b/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) diff --git a/tests/core/test_integration/test_detection_structure_splitting.py b/tests/core/test_integration/test_detection_structure_splitting.py index 6fa83d07..31f04672 100644 --- a/tests/core/test_integration/test_detection_structure_splitting.py +++ b/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