Rolling ball algorithm scales badly with radius

[ctrueden]

Description:

I was playing with rolling ball background subtraction, and naively tried to invoke skimage.restoration.rolling_ball(img, radius=50) on a single-plane 582x200 RGB image. I noticed it eating CPU like crazy, so I did some further testing. Here is what I found:

radius	time (s)	s/r
5	0.627	0.1254
10	3.83	0.383
15	10	0.6666666667
20	20.7	1.035
25	39.4	1.576
30	67	2.233333333

Conversely, the following Groovy script in Fiji takes 17 ms to complete on my system with radius=50 on the same image:

#@ ImagePlus imp
#@ Double(value=50) radius

import ij.plugin.filter.BackgroundSubtracter

start = System.currentTimeMillis()
BackgroundSubtracter bs = new BackgroundSubtracter()
bs.rollingBallBackground(imp.getProcessor(), radius, false, false, false, false, false)
end = System.currentTimeMillis()
println(end - start)

imp.updateAndDraw()

Is this a known limitation of the algorithm as currently implemented in skimage?

Way to reproduce:

img = np.random.random([200, 582, 3])
from skimage.restoration import rolling_ball
rolling_ball(img, radius=50)

Version information:

3.11.6 | packaged by conda-forge | (main, Oct  3 2023, 10:40:35) [GCC 12.3.0]
Linux-6.5.0-28-generic-x86_64-with-glibc2.38
scikit-image version: 0.22.0
numpy version: 1.26.2

I also experienced the same issue using skimage within Pyodide (not sure which versions, sorry).

Given that the default radius for the rolling_ball routine as currently written is 100, I feel like I must be doing something wrong?

[stefanv]

Thanks for the report, @ctrueden! The skimage version clearly needs some work. Despite starting from a faster baseline, do you also see quadratic scaling with radius for the Fiji implementation? If not, it hints that we may need a different implementation.

[ctrueden]

@stefanv Good question! Reading the source, yes, the ImageJ implementation is still quadratic complexity. It has four nested loops: y in [-radius, height+radius], x in [-radius, width+radius], yp in [-radius, radius], and xp in [-radius, radius]. So it's O(r**2 * w * h).

Here are empirical results on a 1024x1024 float32 image:

(I'm not sure what the artifacts are at smaller radii—at first I thought it could be the JDK's JIT, but rerunning the same script yields the same artifacts, so maybe not...)

I looked briefly at the skimage implementation, and it looks reasonable at a glance, using Cython for speed. I see in #4851 that there is even a benchmark. @FirefoxMetzger Does the data above match the expected performance in your experience, or is possible that a regression might have been introduced?

[jni]

@ctrueden these shrinkFactors seem pertinent:

https://github.com/imagej/ImageJ/blob/7c00dc5295ffda22fa98e504a1d6cb228441bbe3/ij/plugin/filter/BackgroundSubtracter.java#L779-L795

😂

[jni]

(It is used before background substraction here and then the image is upscaled again.)

[FirefoxMetzger]

We had an issue like this before. Last time it was different defaults between our implementation and the Fiji (ImageJ ?) implementation for RGB images.

Our version can roll an N-dimensional ball in the image, i.e., it can process RGB images directly. Theirs could - at least at the time - only process grayscale. It's been a while, so I don't remember the exact algorithm but I do recall that they converted the image (HSV perhaps?) applied rolling ball to a single channel, and then converted back to RGB.

We should have the previous issue in our history somewhere ... let's see if we can find it.

[FirefoxMetzger]

Here is the old issue: #5193

Could this be the source of the performance degradation? (aka, it's the curse of dimensionality)

[jni]

Looks like that was a different case of Fiji being much more magical than us. I'm not 100% sure whether it is in scope for skimage to automatically downscale. But we should probably change the default radius at least, 😅 and maybe warn for larger kernels/images that downscaling the image then upscaling the resulting background would be a good approach.

Or we can add relevant kwargs to the function, dunno.

[jni]

(And there is yet a third way Fiji is clever, btw, which is that by default Fiji doesn't use rolling ball but rather what I call "rolling umbrella spokes" 😂, see this imagesc thread.)

[ctrueden]

@FirefoxMetzger @jni Thanks for the links to the prior discussions! Very helpful.

Perhaps it would make sense to put these links into the docstring on the rolling_ball, with an explanation of the limitations of the algorithm itself, as well as a list of alternatives? That way we avoid being "magical" but at the same time educate people on how to "design their own preferred magic" that accomplishes their goals.

Perhaps something like this?

diff --git a/skimage/restoration/_rolling_ball.py b/skimage/restoration/_rolling_ball.py
index 2f8c69e11..c349bf6c5 100644
--- a/skimage/restoration/_rolling_ball.py
+++ b/skimage/restoration/_rolling_ball.py
@@ -53,10 +53,18 @@ def rolling_ball(image, *, radius=100, kernel=None, nansafe=False, num_threads=N
     noise). If this is a problem in your image, you can apply mild
     gaussian smoothing before passing the image to this function.
 
+    This algorithm runs in quadratic time to the size of the ball radius,
+    meaning it can take a long time as the radius grows beyond 30 or so [2, 3].
+    It is an exact N-dimensional calculation; if all you need is an
+    approximation, another option to consider is a top-hat filter [4].
+
     References
     ----------
     .. [1] Sternberg, Stanley R. "Biomedical image processing." Computer 1
            (1983): 22-34. :DOI:`10.1109/MC.1983.1654163`
+    .. [2] https://github.com/scikit-image/scikit-image/issues/5193
+    .. [3] https://github.com/scikit-image/scikit-image/issues/7423
+    .. [4] https://forum.image.sc/t/59267/7
 
     Examples
     --------

[jni]

It's polynomial in the radius, with exponent = image.ndim.

I would like the suggestion to suggest either approximating with top-hat or downscaling-then-upscaling. I think for something generally smooth like background intensity, that's actually a good approach.

Otherwise, improving the docs is always a good idea! 🙏