From 4373a172edeaadd63991b76291152b8975918d16 Mon Sep 17 00:00:00 2001 From: Gregory Lee Date: Sat, 8 Feb 2025 09:00:21 -0500 Subject: [PATCH] add cpu_fallback_threshold to just run very small images via scikit-image --- .../cucim/skimage/morphology/convex_hull.py | 62 +++++++++++++- .../morphology/tests/test_convex_hull.py | 85 +++++++++++++++---- 2 files changed, 126 insertions(+), 21 deletions(-) diff --git a/python/cucim/src/cucim/skimage/morphology/convex_hull.py b/python/cucim/src/cucim/skimage/morphology/convex_hull.py index c0cd35c9..293faf90 100644 --- a/python/cucim/src/cucim/skimage/morphology/convex_hull.py +++ b/python/cucim/src/cucim/skimage/morphology/convex_hull.py @@ -64,7 +64,7 @@ def _check_coords_in_hull( mask_temp = cp.ones((n_coords,), dtype=bool) dot_out = cp.empty((n_coords,), dtype=float_dtype) for idx in range(n_hull_equations): - cp.dot(hull_equations[idx, :ndim], gridcoords, out=dot_out) + cp.dot(hull_equations[idx, :ndim], gridcoords, ou4t=dot_out) dot_out += hull_equations[idx, ndim:] cp.less(dot_out, tolerance, out=mask_temp) mask *= mask_temp @@ -80,6 +80,7 @@ def convex_hull_image( *, omit_empty_coords_check=False, float64_computation=True, + cpu_fallback_threshold=None, ): """Compute the convex hull image of a binary image. @@ -100,19 +101,52 @@ def convex_hull_image( some points erroneously being classified as being outside the hull. include_borders: bool, optional If ``False``, vertices/edges are excluded from the final hull mask. + + + Extra Parameters + ---------------- omit_empty_coords_check : bool, optional If ``True``, skip check that there are not any True values in `image`. + float64_computation : bool, optional + If False, allow use of 32-bit float during the postprocessing stage + that determines whether each pixel falls within the convex hull. + cpu_fallback_threshold : non-negative int or None + Number of pixels in an image before convex_hull_image will fallback + to pure CPU implementation. Returns ------- hull : (M, N) array of bool Binary image with pixels in convex hull set to True. + Notes + ----- + The parameters listed under "Extra Parameters" above are present only + in cuCIM and not in scikit-image. + References ---------- .. [1] https://blogs.mathworks.com/steve/2011/10/04/binary-image-convex-hull-algorithm-notes/ """ + if cpu_fallback_threshold is None: + # Fallback to scikit-image implementation of total number of pixels + # is less than this + cpu_fallback_threshold = 30000 if image.ndim == 2 else 13000 + + if image.size < cpu_fallback_threshold: + # Fallback to pure CPU implementation + from skimage import morphology as morphology_cpu + + return cp.asarray( + morphology_cpu.convex_hull_image( + cp.asnumpy(image), + offset_coordinates=offset_coordinates, + tolerance=tolerance, + include_borders=include_borders, + ) + ) + if not scipy_available: raise ImportError( "This function requires SciPy, but it could not import: " @@ -194,7 +228,13 @@ def convex_hull_image( return mask -def convex_hull_object(image, *, connectivity=2): +def convex_hull_object( + image, + *, + connectivity=2, + float64_computation=False, + cpu_fallback_threshold=None, +): r"""Compute the convex hull image of individual objects in a binary image. The convex hull is the set of pixels included in the smallest convex @@ -216,6 +256,14 @@ def convex_hull_object(image, *, connectivity=2): | / | \ [ ] [ ] [ ] [ ] + Extra Parameters + ---------------- + float64_computation : bool, optional + If False, allow use of 32-bit float during the postprocessing stage + cpu_fallback_threshold : non-negative int or None + Number of pixels in an image before convex_hull_image will fallback + to pure CPU implementation. + Returns ------- hull : ndarray of bool @@ -228,6 +276,9 @@ def convex_hull_object(image, *, connectivity=2): these regions with logical OR. Be aware the convex hulls of unconnected objects may overlap in the result. If this is suspected, consider using convex_hull_image separately on each object or adjust ``connectivity``. + + The parameters listed under "Extra Parameters" above are present only + in cuCIM and not in scikit-image. """ if connectivity not in tuple(range(1, image.ndim + 1)): raise ValueError("`connectivity` must be between 1 and image.ndim.") @@ -238,7 +289,12 @@ def convex_hull_object(image, *, connectivity=2): max_label = int(labeled_im.max()) for i in range(1, max_label + 1): - convex_obj = convex_hull_image(labeled_im == i) + convex_obj = convex_hull_image( + labeled_im == i, + omit_empty_coords_check=True, + float64_computation=float64_computation, + cpu_fallback_threshold=cpu_fallback_threshold, + ) convex_img = cp.logical_or(convex_img, convex_obj) return convex_img diff --git a/python/cucim/src/cucim/skimage/morphology/tests/test_convex_hull.py b/python/cucim/src/cucim/skimage/morphology/tests/test_convex_hull.py index aa3adacb..848d7d87 100644 --- a/python/cucim/src/cucim/skimage/morphology/tests/test_convex_hull.py +++ b/python/cucim/src/cucim/skimage/morphology/tests/test_convex_hull.py @@ -37,13 +37,17 @@ def test_basic(): dtype=bool, ) - assert_array_equal(convex_hull_image(image), expected) + assert_array_equal( + convex_hull_image(image, cpu_fallback_threshold=0), expected + ) def test_empty_image(): image = cp.zeros((6, 6), dtype=bool) with expected_warnings(["entirely zero"]): - assert_array_equal(convex_hull_image(image), image) + assert_array_equal( + convex_hull_image(image, cpu_fallback_threshold=0), image + ) def test_qhull_offset_example(): @@ -181,7 +185,9 @@ def test_qhull_offset_example(): image[nonzeros] = True image = cp.asarray(image) expected = image.copy() - assert_array_equal(convex_hull_image(image), expected) + assert_array_equal( + convex_hull_image(image, cpu_fallback_threshold=0), expected + ) def test_pathological_qhull_example(): @@ -193,7 +199,9 @@ def test_pathological_qhull_example(): [[0, 0, 0, 1, 1, 1, 0], [0, 1, 1, 1, 1, 1, 1], [1, 1, 1, 1, 1, 0, 0]], dtype=bool, ) - assert_array_equal(convex_hull_image(image), expected) + assert_array_equal( + convex_hull_image(image, cpu_fallback_threshold=0), expected + ) @pytest.mark.skipif(True, reason="include_borders option not implemented") @@ -208,7 +216,9 @@ def test_pathological_qhull_labels(): dtype=bool, ) - actual = convex_hull_image(image, include_borders=False) + actual = convex_hull_image( + image, include_borders=False, cpu_fallback_threshold=0 + ) assert_array_equal(actual, expected) @@ -244,7 +254,8 @@ def test_object(): ) assert_array_equal( - convex_hull_object(image, connectivity=1), expected_conn_1 + convex_hull_object(image, connectivity=1, cpu_fallback_threshold=0), + expected_conn_1, ) expected_conn_2 = cp.array( @@ -263,35 +274,42 @@ def test_object(): ) assert_array_equal( - convex_hull_object(image, connectivity=2), expected_conn_2 + convex_hull_object(image, connectivity=2, cpu_fallback_threshold=0), + expected_conn_2, ) with pytest.raises(ValueError): - convex_hull_object(image, connectivity=3) + convex_hull_object(image, connectivity=3, cpu_fallback_threshold=0) - out = convex_hull_object(image, connectivity=1) + out = convex_hull_object(image, connectivity=1, cpu_fallback_threshold=0) assert_array_equal(out, expected_conn_1) def test_non_c_contiguous(): # 2D Fortran-contiguous image = cp.ones((2, 2), order="F", dtype=bool) - assert_array_equal(convex_hull_image(image), image) + assert_array_equal( + convex_hull_image(image, cpu_fallback_threshold=0), image + ) # 3D Fortran-contiguous image = cp.ones((2, 2, 2), order="F", dtype=bool) - assert_array_equal(convex_hull_image(image), image) + assert_array_equal( + convex_hull_image(image, cpu_fallback_threshold=0), image + ) # 3D non-contiguous image = cp.transpose(cp.ones((2, 2, 2), dtype=bool), [0, 2, 1]) - assert_array_equal(convex_hull_image(image), image) + assert_array_equal( + convex_hull_image(image, cpu_fallback_threshold=0), image + ) def test_consistent_2d_3d_hulls(): from cucim.skimage.measure.tests.test_regionprops import SAMPLE as image image3d = cp.stack((image, image, image)) - chimage = convex_hull_image(image) + chimage = convex_hull_image(image, cpu_fallback_threshold=0) chimage[8, 0] = True # correct for single point exactly on hull edge - chimage3d = convex_hull_image(image3d) + chimage3d = convex_hull_image(image3d, cpu_fallback_threshold=0) assert_array_equal(chimage3d[1], chimage) @@ -309,14 +327,18 @@ def test_few_points(): ) image3d = cp.stack([image, image, image]) with assert_warns(UserWarning): - chimage3d = convex_hull_image(image3d, offset_coordinates=False) + chimage3d = convex_hull_image( + image3d, offset_coordinates=False, cpu_fallback_threshold=0 + ) assert_array_equal(chimage3d, cp.zeros(image3d.shape, dtype=bool)) # non-zero when using offset_coordinates # (This is an improvement over skimage v0.25 implementation due to how # initial points are determined without a separate ConvexHull call before # the addistion of the offset coordinates) - chimage3d = convex_hull_image(image3d, offset_coordinates=True) + chimage3d = convex_hull_image( + image3d, offset_coordinates=True, cpu_fallback_threshold=0 + ) chimage3d.sum() > 0 @@ -345,11 +367,12 @@ def test_diamond( offset_coordinates=offset_coordinates, omit_empty_coords_check=omit_empty_coords_check, float64_computation=float64_computation, + cpu_fallback_threshold=0, ) if offset_coordinates: assert_array_equal(chimage, expected) else: - # may not be an exact match if offset coordinates are used + # may not be an exact match if offset coordinates are not used num_mismatch = cp.sum(chimage != expected) percent_mismatch = 100 * num_mismatch / expected.sum() if float64_computation: @@ -378,14 +401,40 @@ def test_octahedron( offset_coordinates=offset_coordinates, omit_empty_coords_check=omit_empty_coords_check, float64_computation=float64_computation, + cpu_fallback_threshold=0, ) if offset_coordinates: assert_array_equal(chimage, expected) else: - # may not be an exact match if offset coordinates are used + # may not be an exact match if offset coordinates are not used num_mismatch = cp.sum(chimage != expected) percent_mismatch = 100 * num_mismatch / expected.sum() if float64_computation: assert percent_mismatch < 5 else: assert percent_mismatch < 20 + + +@pytest.mark.parametrize("radius", [1, 10, 100]) +@pytest.mark.parametrize("offset_coordinates", [False, True]) +@pytest.mark.parametrize("cpu_fallback_threshold", [0, None, 1000000000]) +def test_cpu_fallback(radius, offset_coordinates, cpu_fallback_threshold): + expected = diamond(radius) + + # plus sign should become a diamond once convex + image = cp.zeros_like(expected) + image[:, radius] = True + image[radius, :] = True + + chimage = convex_hull_image( + image, + offset_coordinates=offset_coordinates, + cpu_fallback_threshold=cpu_fallback_threshold, + ) + if offset_coordinates: + assert_array_equal(chimage, expected) + else: + # may not be an exact match if offset coordinates are not used + num_mismatch = cp.sum(chimage != expected) + percent_mismatch = 100 * num_mismatch / expected.sum() + assert percent_mismatch < 5