-
Notifications
You must be signed in to change notification settings - Fork 59
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Add Euclidean distance transform for images/volumes #318
Merged
rapids-bot
merged 12 commits into
rapidsai:branch-22.08
from
grlee77:distance-transform-edt
Jul 28, 2022
Merged
Changes from all commits
Commits
Show all changes
12 commits
Select commit
Hold shift + click to select a range
8e9b5a5
add license info for the PBA+ algorithm
grlee77 0462a2c
Add distance_transform_edt implementation based on PBA+
grlee77 6a9997f
add tests for 2d distance transform
grlee77 3d44899
Update 3D kernels to allow use of size > 1024 per axis
grlee77 8509bff
expand tests and extend to 3D cases
grlee77 0933d73
test error on non-uniform sampling
grlee77 e98b611
flake8/isort fixes
grlee77 e5a37a9
minor update to kernel headers
grlee77 e83d9a7
update outdated comment
grlee77 c67b7c4
add fallback implementation of math.lcm for Python < 3.9
grlee77 cb06d78
make sure fallback lcm function is actually used if needed
grlee77 2a190b3
Apply suggestions from code review
grlee77 File filter
Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,21 @@ | ||
MIT License | ||
|
||
Copyright (c) 2019 School of Computing, National University of Singapore | ||
|
||
Permission is hereby granted, free of charge, to any person obtaining a copy | ||
of this software and associated documentation files (the "Software"), to deal | ||
in the Software without restriction, including without limitation the rights | ||
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell | ||
copies of the Software, and to permit persons to whom the Software is | ||
furnished to do so, subject to the following conditions: | ||
|
||
The above copyright notice and this permission notice shall be included in all | ||
copies or substantial portions of the Software. | ||
|
||
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR | ||
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, | ||
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE | ||
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER | ||
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, | ||
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE | ||
SOFTWARE. |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
3 changes: 3 additions & 0 deletions
3
python/cucim/src/cucim/core/operations/morphology/__init__.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,3 @@ | ||
from ._distance_transform import distance_transform_edt | ||
|
||
__all__ = ["distance_transform_edt"] |
182 changes: 182 additions & 0 deletions
182
python/cucim/src/cucim/core/operations/morphology/_distance_transform.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -0,0 +1,182 @@ | ||
import numpy as np | ||
|
||
from ._pba_2d import _pba_2d | ||
from ._pba_3d import _pba_3d | ||
|
||
# TODO: support sampling distances | ||
# support the distances and indices output arguments | ||
# support chamfer, chessboard and l1/manhattan distances too? | ||
|
||
|
||
def distance_transform_edt(image, sampling=None, return_distances=True, | ||
return_indices=False, distances=None, indices=None, | ||
*, block_params=None, float64_distances=False): | ||
"""Exact Euclidean distance transform. | ||
|
||
This function calculates the distance transform of the `input`, by | ||
replacing each foreground (non-zero) element, with its shortest distance to | ||
the background (any zero-valued element). | ||
|
||
In addition to the distance transform, the feature transform can be | ||
calculated. In this case the index of the closest background element to | ||
each foreground element is returned in a separate array. | ||
|
||
Parameters | ||
---------- | ||
image : array_like | ||
Input data to transform. Can be any type but will be converted into | ||
binary: 1 wherever image equates to True, 0 elsewhere. | ||
sampling : float, or sequence of float, optional | ||
Spacing of elements along each dimension. If a sequence, must be of | ||
length equal to the image rank; if a single number, this is used for | ||
all axes. If not specified, a grid spacing of unity is implied. | ||
return_distances : bool, optional | ||
Whether to calculate the distance transform. | ||
return_indices : bool, optional | ||
Whether to calculate the feature transform. | ||
distances : float32 cupy.ndarray, optional | ||
An output array to store the calculated distance transform, instead of | ||
returning it. `return_distances` must be True. It must be the same | ||
shape as `image`. | ||
indices : int32 cupy.ndarray, optional | ||
An output array to store the calculated feature transform, instead of | ||
returning it. `return_indicies` must be True. Its shape must be | ||
`(image.ndim,) + image.shape`. | ||
|
||
Other Parameters | ||
---------------- | ||
block_params : 3-tuple of int | ||
The m1, m2, m3 algorithm parameters as described in [2]_. If None, | ||
suitable defaults will be chosen. Note: This parameter is specific to | ||
cuCIM and does not exist in SciPy. | ||
float64_distances : bool, optional | ||
If True, use double precision in the distance computation (to match | ||
SciPy behavior). Otherwise, single precision will be used for | ||
efficiency. Note: This parameter is specific to cuCIM and does not | ||
exist in SciPy. | ||
|
||
Returns | ||
------- | ||
distances : float64 ndarray, optional | ||
The calculated distance transform. Returned only when | ||
`return_distances` is True and `distances` is not supplied. It will | ||
have the same shape as `image`. | ||
indices : int32 ndarray, optional | ||
The calculated feature transform. It has an image-shaped array for each | ||
dimension of the image. See example below. Returned only when | ||
`return_indices` is True and `indices` is not supplied. | ||
|
||
Notes | ||
----- | ||
The Euclidean distance transform gives values of the Euclidean distance:: | ||
|
||
n | ||
y_i = sqrt(sum (x[i]-b[i])**2) | ||
i | ||
|
||
where b[i] is the background point (value 0) with the smallest Euclidean | ||
distance to input points x[i], and n is the number of dimensions. | ||
|
||
Note that the `indices` output may differ from the one given by | ||
`scipy.ndimage.distance_transform_edt` in the case of input pixels that are | ||
equidistant from multiple background points. | ||
|
||
The parallel banding algorithm implemented here was originally described in | ||
[1]_. The kernels used here correspond to the revised PBA+ implementation | ||
that is described on the author's website [2]_. The source code of the | ||
author's PBA+ implementation is available at [3]_. | ||
|
||
References | ||
---------- | ||
..[1] Thanh-Tung Cao, Ke Tang, Anis Mohamed, and Tiow-Seng Tan. 2010. | ||
Parallel Banding Algorithm to compute exact distance transform with the | ||
GPU. In Proceedings of the 2010 ACM SIGGRAPH symposium on Interactive | ||
3D Graphics and Games (I3D ’10). Association for Computing Machinery, | ||
New York, NY, USA, 83–90. | ||
DOI:https://doi.org/10.1145/1730804.1730818 | ||
.. [2] https://www.comp.nus.edu.sg/~tants/pba.html | ||
.. [3] https://github.com/orzzzjq/Parallel-Banding-Algorithm-plus | ||
|
||
Examples | ||
-------- | ||
>>> import cupy as cp | ||
>>> from cucim.core.operations import morphology | ||
>>> a = cp.array(([0,1,1,1,1], | ||
... [0,0,1,1,1], | ||
... [0,1,1,1,1], | ||
... [0,1,1,1,0], | ||
... [0,1,1,0,0])) | ||
>>> morphology.distance_transform_edt(a) | ||
array([[ 0. , 1. , 1.4142, 2.2361, 3. ], | ||
[ 0. , 0. , 1. , 2. , 2. ], | ||
[ 0. , 1. , 1.4142, 1.4142, 1. ], | ||
[ 0. , 1. , 1.4142, 1. , 0. ], | ||
[ 0. , 1. , 1. , 0. , 0. ]]) | ||
|
||
With a sampling of 2 units along x, 1 along y: | ||
|
||
>>> morphology.distance_transform_edt(a, sampling=[2,1]) | ||
array([[ 0. , 1. , 2. , 2.8284, 3.6056], | ||
[ 0. , 0. , 1. , 2. , 3. ], | ||
[ 0. , 1. , 2. , 2.2361, 2. ], | ||
[ 0. , 1. , 2. , 1. , 0. ], | ||
[ 0. , 1. , 1. , 0. , 0. ]]) | ||
|
||
Asking for indices as well: | ||
|
||
>>> edt, inds = morphology.distance_transform_edt(a, return_indices=True) | ||
>>> inds | ||
array([[[0, 0, 1, 1, 3], | ||
[1, 1, 1, 1, 3], | ||
[2, 2, 1, 3, 3], | ||
[3, 3, 4, 4, 3], | ||
[4, 4, 4, 4, 4]], | ||
[[0, 0, 1, 1, 4], | ||
[0, 1, 1, 1, 4], | ||
[0, 0, 1, 4, 4], | ||
[0, 0, 3, 3, 4], | ||
[0, 0, 3, 3, 4]]]) | ||
|
||
""" | ||
if distances is not None: | ||
raise NotImplementedError( | ||
"preallocated distances image is not supported" | ||
) | ||
if indices is not None: | ||
raise NotImplementedError( | ||
"preallocated indices image is not supported" | ||
) | ||
scalar_sampling = None | ||
if sampling is not None: | ||
sampling = np.unique(np.atleast_1d(sampling)) | ||
if len(sampling) == 1: | ||
scalar_sampling = float(sampling) | ||
sampling = None | ||
else: | ||
raise NotImplementedError( | ||
"non-uniform values in sampling is not currently supported" | ||
) | ||
|
||
if image.ndim == 3: | ||
pba_func = _pba_3d | ||
elif image.ndim == 2: | ||
pba_func = _pba_2d | ||
else: | ||
raise NotImplementedError( | ||
"Only 2D and 3D distance transforms are supported.") | ||
|
||
vals = pba_func( | ||
image, | ||
sampling=sampling, | ||
return_distances=return_distances, | ||
return_indices=return_indices, | ||
block_params=block_params | ||
) | ||
|
||
if return_distances and scalar_sampling is not None: | ||
vals = (vals[0] * scalar_sampling,) + vals[1:] | ||
|
||
if len(vals) == 1: | ||
vals = vals[0] | ||
|
||
return vals |
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Thank you for updating license texts! :)