Region Properties Performance Overhaul - Part 4: Moment-Based Properties #846
Conversation
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
The functions introduced here are not being added to the public API. They will be used behind the scenes from `regionprops_table` to enable orders of magnitude faster computation of region properties for all labels in an image. The basic approach here is to compute a property for all labels in an image from a single CUDA kernel call. This is in contrast to the approach from the `RegionProperties` class which first splits the full image into small sub-images corresponding to each region and then loops over these small sub-images, computing the requested property for each small region in turn. That approach is not amenable to good acceleration on the GPU as individual regions are typically small. Provides batch implementation that computes the following properties for all properties in a single kernel call: - bbox - label_filled (creates version of label_image with all holes filled) - num_pixels - num_pixels_filled - num_perimeter_pixels (number of pixels at perimeter of each labeled region) - num_boundary_pixels (number of pixels touching the image boundary for each region) The following properties are simple transformations of the properties above and have negligable additional cost to compute: - area - area_bbox - area_filled - equivalent_diameter_area - equivalent_spherical_perimeter (as in ITK) - extent - perimeter_on_border_ratio (as in ITK) - slice The following split the label image into a list of sub-images or subsets of coordinates where each element in the list corresponds to a label. The background of the label image has value 0 and is not represented in the sequences. Sequence entry `i` corresponds to label `i + 1`. In most cases, these will not be needed as properties are now computed for all regions at once from the labels image, but they are provided for completeness and to match the scikit-image API. - coords - coords_scaled - image (label mask subimages) - image_convex (convex label mask subimages) - image_intensity (intensity_image subimages) - image_filled (subimages of label mask but with holes filled) - label (sequence of integer label ids) Test cases are added that compare the results of these batch computations to results from scikit-image `regionprops_table`.
This function operates similarly to `regionprops_table`. In a future commit, once all properties have been supported, it will be used within the existing regionprops_table function so that it will provide much higher performance.
- intensity_mean - intensity_std - intensity_min - intensity_max Both single and multi-channel intensity images are supported
grlee77
force-pushed
the
regionprops_part4_moments
branch
from
15280b1 to
a7f5274
Compare
grlee77
added
improvement
gforsyth
approved these changes
Mar 3, 2025
grlee77
force-pushed
the
regionprops_part4_moments
branch
from
a7f5274 to
df640e3
Compare
gigony
approved these changes
Apr 2, 2025
python/cucim/src/cucim/skimage/measure/_regionprops_gpu_moments_kernels.py
Outdated
Show resolved
Hide resolved
python/cucim/src/cucim/skimage/measure/_regionprops_gpu_moments_kernels.py
Outdated
Show resolved
Hide resolved
python/cucim/src/cucim/skimage/measure/_regionprops_gpu_moments_kernels.py
Outdated
Show resolved
Hide resolved
| tmp2 = m00 - m11; | ||
| tmp2 *= tmp2; | ||
| tmp2 += tmp1; | ||
| tmp2 = sqrt(tmp2); |
There was a problem hiding this comment.
Not sure if this is needed but suggesting here.
Suggested change
| tmp2 = sqrt(tmp2); | |
| // add numerical stability check | |
| tmp2 = max(tmp2, 0.0); // ensure non-negative before sqrt | |
| tmp2 = sqrt(tmp2); | |
| // more robust handling of small values | |
| const double eps = 1e-10; | |
| if (fabs(tmp1) < eps) { | |
| tmp1 = 0.0; | |
| } |
There was a problem hiding this comment.
I think tmp2 and tmp1 individually cannot be negative due to squares that happen above (e.g. tmp1 = m01 * m01 and tmp2 *= tmp2).
I think the problematic case is when the magnitude of tmp1 is within numerical precision of tmp2 (e.g. for a perfectly circular region we would expect equal eigenvalues). The max(tmp1 - tmp2, 0.0) is to handle that case. The other max in max(tmp1 + tmp2, 0.0) seems unnecessary and could be removed.
…ity_kernels.py Co-authored-by: Gigon Bae <gigony@gmail.com>
|
Auto-sync is disabled for draft pull requests in this repository. Workflows must be run manually. Contributors can view more details about this message here. |
These properties are computed based on the image_convex subimages: - area_convex - feret_diameter_max - solidity
…moments_analytical.py
Co-authored-by: Gigon Bae <gigony@gmail.com>
grlee77
force-pushed
the
regionprops_part4_moments
branch
from
433d166 to
aa3a89a
Compare
|
/merge |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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.
Please review #843 first as that explains the general approach in more detail.
Overview
This MR implements many properties based on image moments. These include centroids, normalized/central moments and measurements based on the inertia-tensor. These come in both unweighted and "weighted" variants where the weighted versions rely on values in a corresponding
intensity_image.The primary kernels that run on all image pixels are the computation of
momentsormoments_weighted. Other properties are typically derived from these raw moments and are very fast to compute once the moments of the required order have been computed.Unweighted properties implemented are:
Weighted properties implemented are:
And a couple of properties that are useful, but no currently in scikit-image
axis_major_lengthis axis_lengths[:, 0] while -axis_minor_lengthis axis_lengths[:, -1])Benchmarks
Performance vs. Image Size (with # regions fixed)
The following show performance for a small fixed number of label regions at different spatial scale in both 2D and 3D
In 2D, there are 16 labeled regions for shapes ranging from (64, 64) up to (8192, 8192)
In 3D, there are 8 labeled regions for shapes ranging from (32, 32) up to (512, 512, 512)
Note that "multi-moment" is the time to compute the following list of region properties rather than just a single property
and `"multi-moment-weighted" corresponds to
Performance vs. Region Size (with image shape fixed)
Here a single large 2D image (7680, 4320) is used, but with varying numbers of labeled regions within it. The total % of foreground vs. background voxels remains similar (i.e. regions are larger when there are fewer of them). The number of regions range from 4 up through 16,384.
Here a single large 3D image (384, 384, 384) is used, but with varying numbers of labeled regions within it. The number of regions range from 8 up through 1,728.
Benchmark conclusions
Note: The results for the older GPU-based regionprops from cuCIM are not shown here. However, for many properties that implementation became much slower as the number of regions increased. We can see that for the new implementation proposed here, performance does not continuously decline once a certain number of objects are reached. There is still better acceleration for a small number of objects, but acceleration holds even for very large numbers of objects.