ENH: ndimage: 1D rank filter speed up

[gideonKogan]

Is your feature request related to a problem? Please describe.

During working on the Hampel filter, I implemented a median filter and tested it vs the current implementation. Obviously, there are some tests that I did not implement. Nevertheless, seems like for any window bigger than 5, the new implementation is faster. I checked vs the faster implementation scipy.ndimage.median_filter rather than the slower signal.medfilt

I wonder if I am missing something. Otherwise, is in the community's interest that I will continue the development and replace the current implementation?
Why do I think that the current implementation is slower? I think that they did not use the median heap. Obviously, in case I am right, it might lead to new implementations of rank filter too.

The testing code and the implementation are here:

Describe the solution you'd like.

No response

Describe alternatives you've considered.

No response

Additional context (e.g. screenshots, GIFs)

No response

[ev-br]

There certainly is interest, I'd say!

Ideally the perf improvement happens in ndimage --- the duplication between scipy.ndimage and scipy.signal filtering routines is mostly historical, and would be nice to reduce where reasonable (see e.g. [1]). This is not set in stone, however, so if there's a reason to keep duplicating things, let's hear them.

[1] https://mail.python.org/archives/list/scipy-dev@python.org/thread/Z3Y6CZHGJ2G5GVA3EKLGATDCJX6PDW5P/#T53S74F46GCNEJMSUSAAC4BCR4VQFO55

[gideonKogan]

Well, for both rank filter and median filter, a reasonable 1D use case in signal processing can have a 21 window length. Meaning, that for each step we will replace 1 element from 21. I am not familiar with the image processing but it seems to me that for a window of 6x6=36, every step will have new 6 elements. Therefore a ratio between the new and the current number of elements is smaller and the current code is reasonable (in my tests the advantage os from a widow size of 7 which is larger than 6).
Practiaclly, I was working in the direction of Hampel filter and there is no big effort for me to cover the median and the rank filter for the 1d case. 2D case is kind of changing the direction. I prefer to do 1 step a time :)

[gideonKogan]

@ev-br , in your opinion, assuming I am covering all the boundary conditions (symetric, anti-symetric, zero-pad, constant, etc), what are the additional tests required for the function?

[ev-br]

1 step at a time is good. So the story along the lines of "scipy.signal.medfilt does better in 1D and delegates to ndimage otherwise" sounds reasonable. As a first step at least :-).

Re tests. Good question. Both ndimage and scipy.signal versions are documented to operate on general N-D arrays, so ideally there are tests to clarify and fix the behavior.
Having kernels incommensurate with the array lengths would be great, I guess this is going to be covered in what you call boundary conditions.
Also from a cursory glance at the scipy.signal code, it looks like it matters whether the array and/or kernel size is even or odd. Would be great to test all combinations, too.

Also the question of backwards compatibility is going to come up, so you'll need to be able to show it's preserved (or clearly show where it is not). I'm not entirely sure the current coverage is stellar, so maybe think of how to help in this.
Off the cuff, maybe structure the commits into first tests which pass on main, then tests which may not (are any? hopefully it's all transparent for a user?), then your new implementation.

[gideonKogan]

Seems like the median filters in both modules are running on top of the rank filter. I suggest renaming the issue to rank filter speed up and handle this function for the 1D case by integrating to the source code and show improvement in performance. Is it ok with you @ev-br ?

Also from a cursory glance at the scipy.signal code, it looks like it matters whether the array and/or kernel size is even or odd. Would be great to test all combinations, too.

I figured out that though it though the documentation says that the kernal should be odd, it accepts even values too. Assuming we move to handling the rank filter it will not bother us and I will modify the code to consider it.

Also the question of backwards compatibility is going to come up, so you'll need to be able to show it's preserved (or clearly show where it is not). I'm not entirely sure the current coverage is stellar, so maybe think of how to help in this.
Off the cuff, maybe structure the commits into first tests which pass on main, then tests which may not (are any? hopefully it's all transparent for a user?), then your new implementation.

Assuming we are moving to handle the background functions only, it will be transparent to the user, ok?

[gideonKogan]

PS: a second thought about when it should be more efficient in 2D: I assume that it will outperform the current implementation in the cases where ratio of the computation time is higher than the width of the kernel. For example: all the kernels with a width of 2 and a total kernel size larger than 11 should be beneficial. Or all the kernels with a width of 3 and a total kernel size larger than 19.

[gideonKogan]

Finished handling all the modes/boundary coditions and implemented all the required modifications to handle the rank_filter. Current results are attached. Now I will start moving from ctypes API to Python/C API.
I don't have an experience with this API, meaning that you will probably have to go over my code and make some corrections at the end of the process :(

implementation and testing code

[ev-br]

A couple of quick superficial comments:

1. For the build, you'll need to use meson, see ndimage/meson.build file.
2. It may (or may not be) easier to use cython for the python <-> C glue. Or maybe if it's a single function which has no python user API, a manual C extension is easier. My point is I guess, not everything which ndimage does needs replicating.
3. At some point we'll ask for an ASV benchmark to show the speedup. See https://github.com/scipy/scipy/tree/main/benchmarks/benchmarks
4. Would be nice to document the source of https://github.com/ggkogan/hampel/blob/main/rank_filter/_rank_filter_1d.c

[gideonKogan]

Thanks for your support.

1. I did not manage yet to handle the build
2. To support all the data types I figured out that I have to move to c++ (templates) + cython. C+python/numpy API seems to not support polymorphism. Here I figured out that the current version does not support longdouble and produce wrong results for np.uint64:
   ndimage.rank_filter(np.arange(2**62, 2**62+10).astype('uint64'), 2, size=5) -> np.ones(10, dtype='uint64') * 2**62
3. moving to cython + optimization flag produced further improvement. It seems that now we improve all 1d cases, with respect to the current implementation.
4. added the link to the source.

[ggkogan]

@ev-br I have created a pull request. Will you please review it and tell me what you think of it?
PS: I figured out that the bug associated with longdouble is probably related to the compilation rather than to the code as now it does not work with my implementation too.