Very fast implementation of a median filter using graphics hardware (GPU)

February 8, 2013

Abstract

Median filtering is a well-known method used in a wide range of application frameworks as well as a standalone filter especially for salt and pepper denoising. It is able to highly reduce power of noise without blurring edges too much. It has been widely studied and efficient algorithms and implementations have been proposed. In this paper, we propose a very fast GPU implementation of fixed-size kernel median filters, able to output more than one billion pixels per second in the case of a 3 × 3 pixel kernel. It is based on a Branchless Vectorized Median class algorithm and implemented through a fine tuning of memory and register use allowing us to obtain, as far as we know, the fastest median filter to date.

1 Introduction

First introduced by Tukey in [?], median filtering has been widely studied since then, and many researchers have proposed efficient implementations of it, adapted to various hypothesis, architectures and processors. Originally, its main drawbacks were its compute complexity, its non linearity and its datadependent runtime. Several researchers have adress these issues and designed, for example, efficient histogram-based median filter with predictible runtime [?, ?]. More recently, the advent of GPUs opened new perspectives in terms of image processing performance, and some researchers managed to take advantage of the new graphic capabilities: in that respect, we can cite the Branchless Vectorized Median filter (BVM) [?, ?] which allows very interesting runtimes on CUDA-enabled devices but, as far as we know, the fastest implementation was the sequential, histogram-based, CCDS median filter [?]. When using GPUs as general purpose computing processor, the issue of data transfers brings up, especially when kernel runtime is fast and/or when large data sets are processed. The truth is that, in certain cases, data transfers between GPU and CPU are slower than the actual computation on GPU. It remains that global runtime can still be faster than similar processes run on CPU. Therefore, to fully optimize global runtimes, it is important to pay attention to how memory transfers are done. This leads us to propose, in the following section, an overall code structure to be used with our median kernels. So as to propose concise and more readable code, we will assume the following limitations: 8 or 16 bit-coded gray-level input images whose height (H) and width (W) are both multiples of 512 pixels. Moreover, the following implementation is targeted on Nvidia Tesla GPU (Fermi architecture, compute capability 2.x) but may easily be adapted to other models e.g. those of compute capability 1.3.

2 General structure

Algorithm 1 describes how data are handled in our code. The use of texture memory to store input image data in GPU memory is usual and allows to benefit from the 2-D caching mechanism. After kernel execution, copying output image back to CPU memory is done by use of pinned-memory, which drastically accelerates data transfer.

Algorithm 1: Global memory management on CPU and GPU sides.

1: allocate and populate CPU memory h in; 2: allocate CPU pinned-memory h out; 3: allocate GPU global memory d out; 4: declare GPU texture reference tex img in; 5: allocate GPU array array img in; 6: bind array img in to texture tex img in; 7: copy data from h in to array img in;

8: kernel  gridDim,blockDim  /* to

  d out */;

9: copy data from d out to h out ;

3 Implementing a fast median filter

3.1 Basic principles

Designing a 2-D median filter basically consists in defining a square window H(i, j) for each pixel I(i, j) of the input image, containing n = k×k pixels and centered on I(i, j). The output value I′(i, j) is the median value of the gray level values of the k×k pixels of H(i, j). Figure 1 illustrates this principle with an example of a 5x5 median filter applied on pixel I(5, 6).

Obviously, one key issue is the selection method that identifies the median value. But, as shown in figure 1, since two neighboring pixels share part of the values to be sorted, a second key issue is how to rule redundancy

(a) 5 × 5 median filtering applied on pixel of coordinates (5,6)

(b) window overlapping in 5 × 5 median filtering

Figure 1: Illustration of 5 × 5 median filtering

between consecutive positions of the running window H(i, j). Selection of the median value can be performed by more than one technique, using either histogram-based or sorting methods, each of them having its own benefits and drawbacks.

3.2 Using registers

The aim of using registers to store temporary data results from the perfomance of register access, which is at least 20 times faster than all other types of memory available on the GPU. On fermi architecture, each individual thread can use a maximum of 63 registers, within the limit of 32K per thread block which is the actual peformance limiting factor. That lead us to use the iterative forgetfull selection al-Figure 2: Median value determination by Forgetfull selection process, in a 3 × 3 neighborhood window.

gorithm which allows to save registers while keeping a high pixel throughput value. The principle of the forgetfull selection is, at each iteration, to identify then to eliminate both elements showing the maximum and the minimum values among the current list, then to introduce one of the remaining candidate elements. When there are only three elements in the lists, a simple sort makes it easy to select the median value, which is actually the global median value. Figure 2 illustrates the forgetfull selection process applied to a 3 × 3 pixel median filter. For clarity reasons, the nine values have been represented in a row. The minimum register count R_{n} needed to start this iterative selection process among k × k values is given by:

R_{n} = ⌈\frac{n}{2}⌉ + 1

The selection of both maxima is achieved by a simple incomplete bubble sort implemented through an elementary 2-element swapping function. This ensures that the GPU kernel code is free of divergent branches which would highly impact performance.

3.3 Hiding Latencies

Optimizing a GPU kernel also means hiding latencies. The massive thread parallelism helps in doing so transparently but depending on the actual computation performed by each thread and the grid dimensions, it may be useful to investigate further. We first maximized the Instruction Level Parallelism, especially inside the forgetfull selection method by arranging the instruction sequence in order to reduce the data dependency of consecutive instructions. This helps hiding the latency of arithmetic instruction calls. We also divided the thread block size by 2 while having each thread perform the computation of two neighbor input pixels instead of only one. This keeps the grid size unchanged while reducing the effect of global memory access latency. However, it imposes to manage the window overlapping illustrated by Figure 1 in order to minimize the increase of register count per thread brought by such a new computing organization. Two k × k consecutive neighboring windows share S_{n} = n−\sqrt{n} pixels, this is more than the number of registers needed to perform the median selection, i.e. R_{n} (or equal for 3 × 3 median filter). The (S_{n} − R_{n} + 1) first selection stages can be considered common to both windows, leaving only the k non-shared pixels of each window being processed separately. The above technique saves k + 1 registers for each pair of input pixels, which means that each thread block now uses less registers while processing the same pixel count, thus allowing a higher parallelism level at warp level. Figure 3 illustrates this by representing the different classes of pixels in the 5 × 5 median example.

3.4 Compute complexity

The incomplete sorting (forgetful selection) and the redundancy management performed when outputing two pixels per thread lead to the below instruction count. Arithmetic instructions and texture fetches are easy to count but the number of element swap needed to select the median value is data dependent and only its maximum can be evaluated.

• 5 integer multiplications and 5 integer additions to compute thread indexes and ouput pixel coordinates.

• 2n−1 additions to compute neighbor pixel coordinates.

• n − \sqrt{n} texture fetches to load gray-level values.

• within a m-element vector and according to our selection method, the maximum number of element swapping needed to move the minimum value at first position and the maximum at last is given by:

sc(m) = −2 + ⌈\frac{3.m}{2}⌉

Consequently, the number of element swapping needed by the entire selection of two median values performed by one thread is less or equal to:

( ∑_{i=⌈ \frac{n}{2} ⌉+1}^{n−\sqrt{n}} sc(i)) + 2 ( ∑_{i=n−\sqrt{n}+1}^{n} sc(i))

The technique used allows us to keep runtime in O(n)log_{10}(n) for small window sizes as confirmed by Figure ?? which shows the evolution of GPU kernel runtime when the window size varies.

4 Experiments

Runtimes have been obtained, by averaging 1000 executions, on C2070 GPU card hosted by a system with one Xeon E56202.40GHz processor running a linux kernel 2.6.18. Each kernel has been run on images of size 512×512, 1024×1024, 2048×2048 and 4096×4096. As many authors, we computed the pixel throughput value as our main performance indicator. It includes kernel runtime as well as transfer times, to and from GPU. We also measured Figure 3: Reducing register count in a 5×5 register-only median kernel processing 2 input pixels. The first 7 forgetful selection stages are common to both processed center pixels.

the maximum effective pixel throughput that our GPU/host couple could achieve. This has been done by running a identity kernel that only fetches the gray level of each pixel in texture memory and outputs it in global memory without any other process. Knowing this peak value allow us to evaluate the potential for improvement of each kernel and helps in deciding wether its worth improving further or not. Those peak values are gathered in Table 2, which shows that the larger the image, the higher the expected throughput.

5 Results

The main contribution of this work is to show how to tune a GPU implementation in order to achieve high pixel throughput values. Table 3 shows throughput values of several implementations compared to ours and applied to small window median filters on almost identical GPU cards. Due to the lack of available source code, comparison is based on the most recent published results [?, ?]. In this case, it can be misleading to talk about the relative performance of one algorithm compared to another without knowing much about the implementations. We actually implement one similar algorithm as the one in ArrayFire but the fine tuning made around GPU’s specificities brings important speedup and leads to the fastest implementation of a GPU median filter known to date with 1805 MPix/s. It’s important to notice that this maximum throughput value, obtain by a 3 × 3 median filter comes quite close to the peak effective pixel throughput value of 2444 Mpix/s allowed by our developpement platform and that whatever algorithm run on this GPU would only bring a meager speedup.

Gray-level format→ image size↓ 8 bits 16 bits

512×512  1598  975
1024×1024  2101  1200
2048×2048  2359  1308
4096×4096  2444  1335

Table 2: Maximum effective pixel throughput values (in MPixel per second) on C2070 GPU for each image size.

Gray level format  8 bits  16 bits

time costs→ image size (pixels)↓

{to GPU \atop (ms)}

from GPU (ms)

{Total \atop (ms)}

\begin{matrix}→GPU \\ (ms) \end{matrix}

\begin{matrix}→GPU \\ (ms) \end{matrix}

{Total \atop (ms)}

512×512  0.08  0.06  0.14  0.14  0.10  0.24
1024×1024  0.24  0.19  0.43  0.45  0.35  0.80
2048×2048  0.85  0.68  1.53  1.59  1.32  2.91
4096×4096  3.27  2.61  5.88  6.21  5.21 11.42

Table 1: Time cost of data transfer for each image size and gray-level format on C2070 GPU.

Window size→

image size - perf.↓  3×3  5×5  7×7
512^{2}t(ms)  19.25  14  11
TP (Mpix/s)  18.49  14  11
E (%TP)  20.27  13  10
1024^{2}  75.49  14  11
73.88  14  11
77.40  13  10
2048^{2}  296.18  14  11
294.55  14  11
303.48  14  11
4096^{2}  1184.16  14  11
1158.26  14  11
1213.55  14  11

Table 3: Performance results of the fast median kernels applied on 8 and 16-bit-coded gray-level images.