Yesterday, when I came home after work, I stumbled upon a Github repository with the source code of a graphical user interface for TomoPy. TomoPy is a well-known open-source Python package for tomographic data processing and image reconstruction. I thought the code there looked interesting and in no way related to what I do at work.
Namely, the following code jumped out at me (see here):
if median_filter:
_prj_gpu = ndi_cp.median_filter(_prj_gpu, size=(1, 5, 5))
_sim_gpu = ndi_cp.median_filter(_sim_gpu, size=(1, 5, 5))
if mask_sim:
_sim_gpu = cp.where(_prj_gpu < 1e-7, 0, _sim_gpu)
What is happening here? Let’s assume that both median_filter and mask_sim are true. Basically, we apply a median filter to input data (_prj_gpu is a 2-dimensional array of floats), and then compare the result to some small ‘limit value’ (1e-7), and put either 0 or some other value in _sim_gpu, depending on the result. As Wikipedia says, median filter is a popular way to get rid of noise in an image. Suppose you have an input image. We’ll assume it’s grayscale, so it has only one channel (typically, you can process each channel separately). Choose some odd number, say, 5 (like in the code above). Now, to get the output value at some pixel, imagine a 5×5 pixel square centered at this pixel in the input image. The 25 pixels in this square give you 25 values. Calculate the median of these values—this is the output value.
In other words, every pixel is replaced by the median of the pixels that are not too far from it. This is similar to more traditional blur filters, where every pixel is replaced by the average of the pixels that are not too far from it, but the median version is kind of ‘less blurry’: for example, every output pixel value is actually present somewhere in the original array. In the extreme case, when the input array is black-and-white (all values are either 0.0 or 1.0), the output image after a median filter is also black-and-white, while averaging typically produces lots of intermediate values between 0.0 and 1.0.
Another variation is to change the shape of the region used to calculate median or average (take a circle instead of a square), or use some weighted average (like Gaussian blur), but here we will consider the simplest case of a square median filter. Also, let’s not bother with defining what happens at the pixels close to the edges, where you can’t fit a 5×5 square around a given pixel. The simplest way to deal with that is to assume the value 0 for the pixels outside the input image, and still calculate the median of 25 values.
However, finding the median can be expensive. One simple observation is that finding the median is easier than sorting. You can see it on the example of probably the most popular sorting algorithm: quicksort (you can find its description on Wikipedia).
Let’s assume for simplicity that the numbers in our array are distinct. In general terms, quicksort sorts the array by choosing an element in it (called a pivot), rearranging the values so that pivot comes after all smaller elements and before all bigger elements, and recursing into those two parts: before the pivot and after the pivot. Note that after the rearranging (called partitioning) the pivot appears in its final position: it won’t ever be moved again. Also, the elements before it and after it will be independently sorted, but they won’t mix with each other.
Now imagine that instead of sorting, you just need to find the median. You start applying the quicksort algorithm: choose a pivot, and do the partitioning. Now the array looks like this: first, there are some elements smaller than pivot, then the pivot itself, and then some elements that are bigger than pivot. We are only interested in the median, which is the element that will be in the middle after all the sorting is done. But we can already tell which of the two parts contains the median: it’s the one that has more elements. Of course, it could happen that the pivot divides the array into two equal parts, and then the pivot is the median, and we are done; otherwise, if the first part turns out to have more elements than the second (so the pivot is located to the right of the center), then we can ignore the second part completely: it can’t contain the median. Similarly, if the first part has less elements than the second, the we can ignore it. In any case, it’s sufficient to recurse into only one of the two parts.
There’s an additional complication here: when you recurse, you are no longer looking for the median in the chosen sub-array, but for some other order statistics, meaning, for the kth largest element in that sub-array for some k. So we have to make our algorithm more generic and say that it should be able to find kth largest element in an input array for any k given as additional input. But the procedure described in the previous paragraph still works and allows one to do partitioning and recursion into one of the two sub-arrays. The resulting algorithm is called quickselect (and it’s not esoteric at all, see Wikipedia).
I guess I should talk about asymptotic complexity of quickselect here, but, frankly, it’s not that important. The worst case for quickselect is O(n2), as for quicksort, but it’s obvious from the description that quickselect is faster, because it does less work. Then you can apply some techniques to choose the pivot in a better way, and say some incantations about the ‘average’ case. These techniques for quickselect are mostly the same as for quicksort, though.
So, here comes the twist. In the code above, the very next step in the algorithm takes the results of the median filter processing and compares every pixel value with some fixed constant (1e-7). And then, around 16:43 local time, sitting on my sofa at home, I suddenly realized: we don’t need to calculate the median value, we just need to know if it’s less than 1e-7.
But if you don’t need to know the actual value, here’s what you can do: in a given 5×5 square, go through all the values, and count how many of them are less than 1e-7, and how many are not. If the majority of the values are less than 1e-7, then the median is also less than 1e-7; otherwise, it’s greater or equal to 1e-7. So you can just do one sweep, keep two counters (how many values are greater than 1e-7 and how many are not), and stop when one of the counters reaches 13 (more than half of 25). It might happen that you don’t even need to look at all the values in the square!
So, here are some benchmarks on a 800×600 image and 5×5 median filter. Three methods are tested: using array sorting to calculate the median and then compare with the limit value; same with quickselect instead of sorting; and, finally, the method described above that avoids calculation of the median altogether. As you can see, the quickselect version is almost 1.3 times faster than quicksort, and the last version is more than 3 times faster than quickselect.
| Method | Mean | Error | StdDev |
|----------------- |----------:|---------:|---------:|
| UsingSort | 181.76 ms | 1.677 ms | 1.569 ms |
| UsingQuickSelect | 144.36 ms | 2.089 ms | 1.954 ms |
| WithoutMedian | 42.80 ms | 0.046 ms | 0.043 ms |
The lesson here, believe it or not, sounds trivial: don’t calculate what you don’t need. Or, if you wish, don’t solve the general problem, solve your particular case. I’m sure one can also extract some other lessons from this story. For example, are we as hobbyist developers, who are just playing with random Python code after we are done with our daily jobs, capable of noticing the possibility of fusing median calculation with comparison? I believe that some blame for not noticing things like that should be put on modern software development practices, like the so-called ‘single responsibility principle’: developers are taught to split problems into smaller tasks and solve those separately. In fact, it’s not even hard to imagine a situation where the two tasks (calculation of the median and comparison with a limit) are given to different developers or even different teams (or get deployed to different microservices). It sounds absurd, but some friends told me that it’s a grim reality for many software developers today. And you can see how it’s the total opposite of what is needed to achieve tremendous improvements.
Back to top