Skip to content
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

Apply fixes to skimage.transform scheduled for scikit-image 0.19.2 #208

Merged
merged 2 commits into from
Mar 25, 2022

Conversation

grlee77
Copy link
Contributor

@grlee77 grlee77 commented Jan 26, 2022

Closes #199

Corresponds to:
scikit-image/scikit-image#6207
scikit-image/scikit-image#6211
scikit-image/scikit-image#6214

Additionally, the affines contained within the PiecewiseAffineTransform will have parameters in CuPy arrays when the inputs to estimate are CuPy arrays. This is one to make it consistent with the other transform classes.

@chrisroat, can you confirm if this fixes the issue for you?

@grlee77 grlee77 requested a review from a team as a code owner January 26, 2022 21:10
@grlee77
Copy link
Contributor Author

grlee77 commented Jan 26, 2022

@chrisroat, Can you provide some context of the sizes of the src and dst arrays that are being used with the estimate method in your application? I did not spend much time thinking about this _geometry.py module on the cuCIM side and it is probably one area where we are unlikely to get a benefit over staying on the CPU given the small sizes of the affine matrices. PiecewiseAffineTransform has a particular difficulty in that it relies on methods from scipy.spatial.Delaunay that do not have a corresponding GPU implementation in CuPy, requiring us to transfer to/from the host in multiple places!

For that reason this _geometry.py is possibly the only place in cuCIM where the classes accept both NumPy or CuPy arrays as inputs. Most of the rest of the library only allows CuPy arrays. This is probably something we should discuss further to decide on the desired behavior long term.

@chrisroat
Copy link

Thanks for consolidating everything, and adding the conversion to cupy arrays. My incremental approach was too slow. :)

The application in question is warping a chunk of data to a new coordinate system. The code basically looks like the following. The shape is as large as 128x512x512, which gives coords roughly 100M elements.

    ndim = chunk.ndim
    shape = chunk.shape

    coords = np.indices(shape, np.int32)
    coords = coords.reshape(ndim, -1).T
    coords = tform(coords)
    coords = coords.T.reshape((-1,) + shape)
    ...
    return ndimage.map_coordinates(chunk, coords, ...)

@grlee77 grlee77 added breaking Introduces a breaking change bug Something isn't working labels Jan 26, 2022
grlee77 and others added 2 commits January 26, 2022 16:41
corresponds to scikit-image gh-6207 and gh-6214

Co-authored-by: Chris Roat <1053153+chrisroat@users.noreply.github.com>
incorporates the change from scikit-image gh-6211

Also updates PiecewiseAffineTransform to be consistent with the other transform types
in keeping self.params as a CuPy array when the inputs are CuPy arrays. Note that
performance of this transform will be very slow due to host/device transfers to
access the CPU-only scipy.spatial.Delaunay.

Co-authored-by: Chris Roat <1053153+chrisroat@users.noreply.github.com>
@chrisroat
Copy link

chrisroat commented Jan 26, 2022

Oh, sorry -- I realize now that you were asking about the estimate function, not the __call__. Yes, the arrays are small that calculate the transform -- typically less than 1k elements. I could use CPU skimage to calculate the transform, and then switch over to cucium based affines for the use of the transform and the warp.

Copy link
Member

@jakirkham jakirkham left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks Greg! 😄

Had a couple questions below

Comment on lines +1062 to +1063
# TODO: update if spatial.Delaunay become available in CuPy
self._tesselation = spatial.Delaunay(cp.asnumpy(src))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we open an issue about this?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We could, it is probably not very easy to address though. We would need a Delaunay triangulation class with a vertices attribute and find_simplex method for the GPU (There are several other attributes/methods on this class in SciPy, but they are not used here). On the SciPy side the implementation relies on the QHull C library.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm sure there are GPU implementations of Delaunay triangulation, but hadn't looked into whether there is a compatibly licensed one we could reuse.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we did find one, ideally it would be available via upstream cupyx.scipy.spatial

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Idk maybe it could happen. There has been work in the past like CudaHull


# determine triangle index for each coordinate
simplex = self._tesselation.find_simplex(coords)
# coords must be on host for calls to _tesselation methods
simplex = self._tesselation.find_simplex(cp.asnumpy(coords))
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why do we coerce to a NumPy array here?

Note: This comes up again below, but only commenting here for simplicity.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

self._tesselation is of class Delaunay from SciPy so this find_simplex method is wrapping some CPU implementation from QHull under the hood.

@grlee77
Copy link
Contributor Author

grlee77 commented Jan 27, 2022

Oh, sorry -- I realize now that you were asking about the estimate function, not the call. Yes, the arrays are small that calculate the transform -- typically less than 1k elements. I could use CPU skimage to calculate the transform, and then switch over to cucium based affines for the use of the transform and the warp.

Yeah, that is probably the best approach here.

@jakirkham
Copy link
Member

So RAPIDS is releasing next week. Should we retarget this for the next release?

@jakirkham jakirkham changed the base branch from branch-22.02 to branch-22.04 February 22, 2022 20:13
@jakirkham
Copy link
Member

Retargeted to branch-22.04 per our discussion in the meeting today. Though ok to include PyPI 22.02.01

@jakirkham
Copy link
Member

@gpucibot merge

@rapids-bot rapids-bot bot merged commit 7ead774 into rapidsai:branch-22.04 Mar 25, 2022
@jakirkham
Copy link
Member

Thanks Greg! 😄

Let's follow up on anything else in a new PR 🙂

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
breaking Introduces a breaking change bug Something isn't working
Projects
None yet
Development

Successfully merging this pull request may close these issues.

[BUG] PiecewiseAffineTransform can trigger TypeError
3 participants