Add Euclidean distance transform for images/volumes #318

grlee77 · 2022-06-28T16:46:07Z

Overview

closes #319

This PR adds an implementation of scipy.ndimage.distance_transform_edt. For each foreground pixel in a binary image, this function computes the minimal Euclidean distance to reach a background pixel. This is a SciPy function rather than a scikit-image one so I have put it under cucim.core instead of cucim.skimage. Longer term we could move this upstream to CuPy, but I think we want to have a copy here so that we can use it in the near term to implement some of the missing scikit-image functionality. This function is used by the following functions in the scikit-image API, so this PR will help us implement these in future PRs:

skimage.morphology.medial_axis (update: see #342)
skimage.segmentation.expand_labels (update: see #341)
skimage.segmentation.chan_vese (update: see #343)
It is also used in examples related to skimage.segmentation.watershed and is needed for implementation of ITK's SignedMaurerDistanceMapImageFilter for itk-cucim.

The algorithm used here is adapted from the C++ code in the PBA+ repository (MIT-licensed)

extensions made to PBA+ kernel source

larger sizes can be supported due to runtime code generation of the integer datatypes and encode/decode functions.

current known limitations

2D and 3D only (this should cover the most common use cases)
The sampling argument is not fully supported. This can likely be done with a bit of effort, but will require modifications to the CUDA kernels
User-specified output arguments are not currently supported. We can potentially add this in the future
The indices returned by return_indices=True are equally valid, but not always identical to those returned by SciPy. This is because there can be multiple indices with an identical distance, so which one gets chosen in that case is implementation-dependent.

initial benchmarks relative to SciPy

Here % true is the percentage of the image that corresponds to foreground pixels. Even for fairly small images there is some benefit, with the benefit becoming two orders of magnitude at larger sizes.

shape	% true	cuCIM	SciPy	acceleration
(256, 256)	5	0.00108	0.00353	3.25
(512, 512)	5	0.00108	0.01552	14.42
(2048, 2048)	5	0.00252	0.34434	136.86
(4096, 4096)	5	0.00765	1.58948	207.88
(32, 32, 32)	5	0.00103	0.00305	2.98
(128, 128, 128)	5	0.00153	0.30103	196.26
(256, 256, 256)	5	0.00763	3.17872	416.37
(384, 384, 384)	5	0.02460	14.28779	580.89
(256, 256)	50	0.00107	0.00430	4.01
(512, 512)	50	0.00109	0.01878	17.30
(2048, 2048)	50	0.00299	0.39304	131.60
(4096, 4096)	50	0.00896	1.84686	206.19
(32, 32, 32)	50	0.00102	0.00361	3.53
(128, 128, 128)	50	0.00163	0.31657	194.66
(256, 256, 256)	50	0.00914	3.35914	367.49
(384, 384, 384)	50	0.03005	13.83219	460.30
(256, 256)	95	0.00107	0.00344	3.22
(512, 512)	95	0.00108	0.01638	15.17
(2048, 2048)	95	0.00314	0.36996	117.81
(4096, 4096)	95	0.00943	1.90475	202.05
(32, 32, 32)	95	0.00102	0.00314	3.07
(128, 128, 128)	95	0.00180	0.28843	159.80
(256, 256, 256)	95	0.01073	3.23450	301.41
(384, 384, 384)	95	0.03577	12.40526	346.82

other comments

can likely reduce memory overhead and improve performance a bit by refactoring some of the pre/post-processing code into elementwise kernels. (e.g. encode3d/decode3d, etc.)
(JK) This may be able to leverage CuPy in the future ( scipy.ndimage.distance_transform_edt (Euclidean distance transform) cupy/cupy#6919 )

fix automated setting of m1, m2, m3 parameters for 2D PBA+

When switching to 64-bit int can go up to size 2**20 on each axis.

grlee77 · 2022-06-28T18:54:20Z

The failure here is due to use of math.lcm in a helper function. Apparently this function was only recently introduced in Python 3.9, so will need to introduce a workaround for that.

grlee77 · 2022-07-26T14:07:47Z

I am pretty confident that this is working as expected based on the tests here as well as use in other pending feature PRs.

It would be good to review and merge this soon so I can submit a series of smaller new feature PRs that rely on this functionality (issues #334, #335, #336).

gigony

I don't understand detailed algorithms so I feel I am not capable of reviewing this patch in detail.
Overall, it looks good to me.

gigony · 2022-07-27T01:30:39Z

LICENSE-3rdparty.md

@@ -281,3 +281,9 @@ StainTools
 - https://github.com/Peter554/StainTools/blob/master/LICENSE.txt
 - Copyright: Peter Byfield
 - Usage: reference for stain color normalization algorithm
+
+PBA+


Thank you for updating license texts! :)

python/cucim/src/cucim/core/operations/morphology/tests/test_distance_transform.py

remove irrelevant text (copy/paste error)

grlee77 · 2022-07-27T19:42:03Z

I don't understand detailed algorithms so I feel I am not capable of reviewing this patch in detail.
Overall, it looks good to me.

Thanks @gigony. I think the kernels are pretty well tested both here and in the original PBA+ repo, so probably doesn't need to much review in the details of the kernels in the .h files. I can hopefully revisit adding support for non-uniform values for sampling in a follow-up PR (for a future release).

This PR is a bit unique for cuCIM in that I kept the bulk of the kernel code in the form of CUDA C++ header files so that it is easier to view than if it was stored within a Python string. The file gets loaded into a string for purposes of runtime compilation, but keeping it in a separate file makes it easier to view with proper code highlighting either on GitHub or in an IDE.

jakirkham · 2022-07-27T22:19:08Z

Made an edit to the OP to xref the related upstream CuPy issue. Hope that is ok 🙂

jakirkham

Thanks Greg! 🙏

Would suggest converting any TODOs into issues so we can track them better

jakirkham · 2022-07-28T06:46:39Z

@gpucibot merge

closes #334 This adds a simple function that can be used to expand a set of labels in a label image. It should be reviewed after #318 is merged. The new commits only start from 6241bd2. Authors: - Gregory Lee (https://github.com/grlee77) Approvers: - https://github.com/jakirkham URL: #341

closes #335 This PR adds a function for segmentation based on active-contours. Test cases are directly from scikit-image, but the code itself was heavily refactored to make extensive use of kernel fusion via `cupy.fuse()`. I will add benchmark results here soon. It should be reviewed after #318 is merged. The new commits only start from 502f41a Authors: - Gregory Lee (https://github.com/grlee77) Approvers: - https://github.com/jakirkham URL: #343

closes #336 This PR adds a function for skeletonization of 2D images via the medial axis transform. It should be reviewed after #318 is merged. The new commits only start from 19a6fed. There is one sequential component to this algorithm that still must be run on the CPU, but the majority of the computations are on the GPU and acceleration is good. I will add benchmark results here soon. Authors: - Gregory Lee (https://github.com/grlee77) - https://github.com/jakirkham Approvers: - https://github.com/jakirkham - Gigon Bae (https://github.com/gigony) URL: #342

jakirkham · 2023-02-13T20:35:53Z

Filed upstream issue ( cupy/cupy#7400 ) about moving this to CuPy

grlee77 added 9 commits June 27, 2022 11:17

add license info for the PBA+ algorithm

8e9b5a5

Add distance_transform_edt implementation based on PBA+

0462a2c

add tests for 2d distance transform

6a9997f

fix automated setting of m1, m2, m3 parameters for 2D PBA+

Update 3D kernels to allow use of size > 1024 per axis

3d44899

When switching to 64-bit int can go up to size 2**20 on each axis.

expand tests and extend to 3D cases

8509bff

test error on non-uniform sampling

0933d73

flake8/isort fixes

e98b611

minor update to kernel headers

e5a37a9

update outdated comment

e83d9a7

grlee77 added feature request New feature or request non-breaking Introduces a non-breaking change labels Jun 28, 2022

grlee77 requested a review from a team as a code owner June 28, 2022 16:46

grlee77 mentioned this pull request Jun 28, 2022

[FEA] add euclidean distance transform #319

Closed

add fallback implementation of math.lcm for Python < 3.9

c67b7c4

grlee77 mentioned this pull request Jun 29, 2022

Implement distance function filters InsightSoftwareConsortium/itk_cucim#13

Open

make sure fallback lcm function is actually used if needed

cb06d78

This was referenced Jul 26, 2022

Add cucim.skimage.segmentation.expand_labels #341

Merged

Add cucim.skimage.morphology.medial_axis #342

Merged

Add segmentation with the Chan-Vese active contours method #343

Merged

gigony approved these changes Jul 27, 2022

View reviewed changes

grlee77 commented Jul 27, 2022

View reviewed changes

python/cucim/src/cucim/core/operations/morphology/tests/test_distance_transform.py Outdated Show resolved Hide resolved

grlee77 commented Jul 27, 2022

View reviewed changes

python/cucim/src/cucim/core/operations/morphology/tests/test_distance_transform.py Outdated Show resolved Hide resolved

Apply suggestions from code review

2a190b3

remove irrelevant text (copy/paste error)

grlee77 mentioned this pull request Jul 27, 2022

scipy.ndimage.distance_transform_edt (Euclidean distance transform) cupy/cupy#6919

Closed

jakirkham approved these changes Jul 28, 2022

View reviewed changes

rapids-bot bot merged commit 1d784ab into rapidsai:branch-22.08 Jul 28, 2022

grlee77 mentioned this pull request Sep 1, 2022

Implement isotropic versions of binary morphologic operators scikit-image/scikit-image#6492

Merged

6 tasks

jakirkham mentioned this pull request Feb 13, 2023

Adding distance_transform_edt (from scipy.ndimage) cupy/cupy#7400

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add Euclidean distance transform for images/volumes #318

Add Euclidean distance transform for images/volumes #318

grlee77 commented Jun 28, 2022 •

edited

Loading

grlee77 commented Jun 28, 2022

grlee77 commented Jul 26, 2022

gigony left a comment

gigony Jul 27, 2022

grlee77 commented Jul 27, 2022

jakirkham commented Jul 27, 2022

jakirkham left a comment

jakirkham commented Jul 28, 2022

jakirkham commented Feb 13, 2023

Add Euclidean distance transform for images/volumes #318

Add Euclidean distance transform for images/volumes #318

Conversation

grlee77 commented Jun 28, 2022 • edited Loading

Overview

extensions made to PBA+ kernel source

current known limitations

initial benchmarks relative to SciPy

other comments

grlee77 commented Jun 28, 2022

grlee77 commented Jul 26, 2022

gigony left a comment

Choose a reason for hiding this comment

gigony Jul 27, 2022

Choose a reason for hiding this comment

grlee77 commented Jul 27, 2022

jakirkham commented Jul 27, 2022

jakirkham left a comment

Choose a reason for hiding this comment

jakirkham commented Jul 28, 2022

jakirkham commented Feb 13, 2023

grlee77 commented Jun 28, 2022 •

edited

Loading