A faster 5x5 median filter for OpenCV

from a human

OpenCV's 5x5 median blur was doing more work than it needed to, so I swapped in a smaller network that someone else had already found. It's about 1.36x faster on my Mac and the output is the same for integer images. The PR is still open, so we'll see.

OpenCV's medianBlur, for a 5x5 window, finds the median of 25 pixels with a fixed sequence of 113 compare-exchange steps. I have a pull request open that replaces it with Bert Dobbelaere's sequence of 85. On an Apple M1 that is about 1.36x faster in OpenCV's own perf tests, and for integer images the output is bit-identical. As I write this the PR is under review and not merged.

What a median network is

A compare-exchange takes two values and puts the smaller in one slot and the larger in the other. A sorting network is a fixed list of them, the same list for every input, that leaves the values sorted. A median network only has to leave the median in one particular slot, so it can skip work a full sort would need.

That fixedness is why it suits SIMD: there are no branches. In a vector register each lane holds one pixel position, a compare-exchange is one vector min and one vector max, and one pass of the list filters a whole register of pixels at once. Run time is then close to the number of min/max instructions divided by how many the core can issue per cycle. Cut the instruction count and the loop gets faster by about that ratio, as long as it stays bound by min/max throughput and doesn't start spilling registers to the stack.

Some compare-exchanges only matter for one output. If the max of a pair is never used again, the compiler drops that half. Going from 25 live values down to one takes exactly 24 such halves, so a median network of 25 inputs costs twice its comparator count minus 24. OpenCV's 113 compare-exchanges become 202 min/max operations (89 full pairs and 24 halves), not 226.

The swap

The network in the PR is Bert Dobbelaere's 85-comparator median of 25, from his page of median selection networks. He found it with his SorterHunter program, and it is in the SorterHunter repository under the MIT license. It has been on his page since August 2022. I checked with him before opening the PR and he was fine with it being used.

It comes to 146 min/max operations against 202, a ratio of 1.38. The old network has been in OpenCV since 2009.

The comparator order needed work. On x86 with AVX2's 16 vector registers, the order of the comparators decides how many values get spilled to the stack. In the order Dobbelaere published it, the network costs 35 store/reload pairs per loop iteration on AVX2 for the 8-bit and 16-bit types, where OpenCV's current network costs 13. Claude searched over orderings of the same 85 comparators that respect their dependencies and keep live ranges short, and the order in the PR gets 11. For 32-bit floats the counts are 70, about 52 and 23. NEON and AVX-512 have enough registers to spill nothing in any of these versions, and on NEON the reordering costs about 1%. Claude got the x86 counts by compiling the AVX2 and AVX-512 dispatch units and reading the assembly. There was no x86 machine to time it on; the only x86 timings are a reviewer's, below.

The old network was also copied into the two OpenCL kernels. Codex's review of the patch caught that, and the PR changes them too. The RISC-V and Andes HAL copies are untouched because I can't test them; they would be a follow-up.

Checking it

Claude found the network on Dobbelaere's page, generated the patch, wrote the checkers and ran the benchmarks. I reviewed the patch, and it's mine to answer for.

A comparator network is a median selector for every input exactly when it works on all inputs made of zeros and ones (the 0-1 principle), which for 25 values is 2^25 cases. Every comparator list in the patched files was extracted back out of the source and run through two separately written checkers on all 2^25 inputs. Dropping the last comparator makes each check fail. Beyond that, the patched filter matches a sort-based reference bit for bit on 3,923 image cases per depth (8U, 16U, 16S and 32F, one to four channels, widths from 1 to 70 so the border columns and the vector tail are covered), and OpenCV's existing median tests pass. The OpenCL kernels pass the in-tree test, but only on a CPU OpenCL implementation (PoCL), never a GPU.

A reviewer on the PR pointed out that the in-tree accuracy test only covers 8-bit input, and that the OpenCL test for the other depths compares two paths that now share a network. A ksize-5 accuracy test for 16U, 16S and 32F against a sort-based reference was added after review; it passes on the old and the new network, and fails on all 12 type and channel combinations when one compare-exchange is removed.

Floating point has two exceptions. A NaN input still gives a result that depends on the network and the instruction set's min/max rules, as it did before, and a plus or minus zero may come out with the other sign.

Speed

On the M1, opencv_perf_imgproc for medianBlur at ksize 5 gives a geometric mean speedup of 1.36x over 40 configurations, ranging from 1.19x to 1.56x. That build has KleidiCV turned off so that OpenCV's own kernel runs. The 3x3 code is unchanged, and its rows have a geometric mean of 0.97x, which is a fair picture of the noise on a shared, busy machine. A standalone loop with both networks in one process measures 1.37x to 1.40x, close to the 1.38 operation ratio.

The same reviewer ran the perf tests on an 11th-generation Intel i9 and posted the table in the PR. Computed from that table, the 40 ksize-5 rows have a geometric mean of 1.30x, and 38 of them land between 1.15x and 1.59x. Two don't: 640x480 16UC1 came out at 1.01x, and 640x480 8UC4 was slower, at 0.84x. The ksize-3 rows in the same table, whose code didn't change, range from 0.65x to 2.05x, so one row on its own doesn't say much.

If your OpenCV build has KleidiCV, as AArch64 Linux, Android and macOS builds do by default, it handles ksize 7 and below itself (its 5x5 network is 188 operations) before any of this code runs, except for in-place calls, so the PR only affects the calls it doesn't take. x86 builds and builds without KleidiCV run OpenCV's own kernel.

Is 85 the minimum

Claude looked for a smaller network and didn't find one. The best known for 25 inputs is still 85 comparators, and the best lower bound I know of is 52. What Claude did get is a narrower statement about Dobbelaere's network itself: with its first 74 comparators fixed, no tail shorter than its last 11 works. Two different exact methods, a breadth-first search and a SAT encoding, agree on the last 5 and the last 10. The 11 was checked by SAT alone, in two encodings, with no proof certificate. None of this says whether some other network with 84 comparators exists.

A side computation covers small cases. Claude's exhaustive search found that a network leaving the median of 5, 7 or 9 inputs on one wire needs at least 7, 13 and 19 compare-exchanges, and known networks have those sizes. OpenCV's 3x3 path uses 19 for 9 inputs, so under that model it is already optimal. A second program with the same design and separate code reproduces the results. There is no checkable certificate, so this is a computation, not a proof. The 19 isn't a new claim either: Dobbelaere's page already lists it as optimal for 9 inputs, citing Smith (1996).