-Median filtering technique is often used to remove additive white, salt and pepper noise from a signal or a source image. This filtering method is essential for the processing of digital data representing analog signals in real time. The median filter considers each pixel in the image in turn and looks at its nearby neighbors to determine whether or not it is representative of its surroundings. It replaces the pixel value with the median of neighboring pixel values. The median is calculated by first sorting all the pixel values from the surrounding neighborhood into numerical order and then replacing the pixel being considered with the middle pixel value. We have used graphics processing units (GPUs) to implement the post-processing, performed by NVIDIA Compute Unified Device Architecture (CUDA). Such a system is faster than the CPU version, or other traditional computing, for processing medical applications such as echography or Doppler. This paper shows the effect of the Median Filtering and a comparison of the performance of the CPU and GPU in terms of response time. Index Terms-GPU programming, performances, filters, images
The noise and distortion are the main factors that limit the
capacity of data transmission in telecommunications affecting
the accuracy of the results in the signal measurement systems,
in communications and signal processing [
The salt and pepper noise or impulsive noise is visible in images and is similar to many black and white spots and is generated during image acquisition. It is characterized by positive or negative peaks that saturate the pixels of the image. An example of salt and pepper noise in images is shown in figure 6b. An example of white noise audio signal is shown in figure 7a. A white noise is also a sound that contains every frequency within the range of human hearing and it can be generated on a sound synthesizer. Sound designers can use this sound, with some processing and filtering, to create a multitude of effects such as wind, surf, space whooshes, and rumbles. The analog voltage signal produced by the phototransistors is converted into a digital or numeric value, the ADC adds broad spectrum noise or white noise. The Median filter is a nonlinear digital filtering technique often used to remove salt and pepper noise.
Computing the neighborhood requires more computation time in a post processing step. More recently, the opened perspectives have been offered by modern GPUs that allow us to develop CUDA-based filters, which will be shown as the fastest median filter implementation known. The surgical applications such as echography or Doppler to biomedical research and clinical medicine for real-time image request image brightness and exceptional image quality. In this paper, the first section explains the filtering problem and mathematical model. In CUDA section, we have presented the algorithm implemented using CUDA of median filter comparing the performance of CPU implementations of filtering effects, analysing one dimensional (1D) audio signal.
II. ANALOG TO DIGITAL CONVERSION
In Signal Theory, Sampling is a technique that allows the
conversion of a continuous signal (for example audio) to a
discrete-time signal evaluating the amplitude of that signal at
regular time interval (sample time, ts) [
Shannon Theorem establishes the minimum sampling
frequency fs so that the continuous time signal can be
reconstructed from the discrete time samples without distortion and
error [
Digital filters [
The main types of digital filters are: Infinite Impulse Response (IIR) filters, and Finite Impulse Response (FIR) filters. The first class of digital filters is distinguished by having an impulse response which continues indefinitely. A FIR filter is a causal LTI system whose impulse response is of finite duration; the transfer function of a FIR filter turns out to be a polynomial in z 1. Essentially, the filtering process is a function in which the value that assumes the output sample is determined by h(n) that operates on neighbors too. The filter operates on a window w of N elements with a width smaller than the signal duration. For this reason it is called local operator. Two classes of filter are distinguished: Linear and Non Linear. The linear filtering operation returns an output that is a linear combination of the input. A non linear filter is a filter whose output is not a linear function of its input and it is based on rank function for operation of erosion and dilation. The Fourier transform is a reversible, linear transform used only for linear filters.
In an invariant-time system, in order to filter a signal we need to implement the convolution of the input signal x(t) with the filter transfer function h[t] y(t) = x(t) h(t) (2)
A discrete-time linear system transforms the input succession x[n] 2 Z to the output succession y[n] 2 Z across the discrete convolution with the discrete transfer function h[n] y[n] = 1 X x[k] h[n
k] = so that for a discrete-time filter the output is a weighted sum of the values taken from the input at the current time and at earlier times. This is described by the following operation: y[n] = h0x[n]+h1x[n 1]+ +hN x[n
N ] = (4) where hi are the filter coefficients that determine the impulsive response and N is the filter order.
N X hix[n i] i=0
A simple example of FIR filter is the mobile mean filter, the coefficients b0; ; bN with N=2 , are determined by equation (5) bi =
1 (N + 1) (5)
Median filtering [
Fig. 4: Example of Median filtering operation on a window of 5 samples. After sorted, the central element is the Median, for this reason the window elements must be odd.
Two examples of median filtering are shown in figure 6 for an image (2D domain) and in figure 7 for an audio signal (1D domain). It should be recalled that if m is the median among an ordered set of numbers, the same amount C of numbers will precede and follow that median. If the cardinality of that set is N then it follows that
C = (N 1)=2 (6) then C represent the width of the domain boundary. The condition for the 1D signal is that the first term of the samples C and the last term C can not be filtered because they have no neighbours. These samples will be copied without any processing. The novelty consists in determining the median element without ordering the window. For each mi 2 w it is useful to count the number of elements equal to, greater and lower than each current element in the window, mi 2 w is the median if 9 exactly C elements > mi 9 at least C + 1 elements = mi 9 exactly C elements < mi (7a) (7b) (7c) where w denotes the window and mi represents the i-th element of w, i = f1 : : : N g. Of course, all the said conditions must be satisfied.
The computational complexity (CC) defines the number of
operations required for the extraction of the median from
the window of N elements. The CC in an efficient sorting
algorithm is estimated as N ln(N ) [
Counting algorithm equals sorting algorithm in terms of computational complexity when N=10.
The developer still programs in the familiar C, C++ and
incorporates extensions of these languages in the form of
a few basic keywords [
In this section the main differences between the CPU and the GPU are analysed. Originally GPUs—Graphic Processing Units—were used to accelerate the drawing of lines, circles and arcs, rectangles, and character bitmaps, then were also used to perform computations. Instead the CPU architecture allows performing computations directly and it is more evolved and complex with respect to the GPU. GPUs have far more processor cores than CPUs, but because each GPU core runs significantly slower than a CPU core and does not have the features needed for modern operating systems, it is not appropriate for performing most of the processing in everyday computing. GPUs are more suited to compute-intensive operations such as video processing and physics simulations. A GPU is simpler than a CPU constituted by full adders, multipliers and comparators. A program is executed in an element called thread. The GPU has multiple hardware threads per core comparing to CPU that it is constituted by multicore (2,4,8 or 16 ) and single thread architecture. The CPU architecture has multi-core and multi-thread at the same time.
Computing threads on a GPU are organized into thread blocks and grid. Each thread belonging to the same block communicates with others via a shared memory and can be synchronized manually by the user. CUDA contains special C function called kernels, which are simply C code that is executed on the GPU on fixed number of threads. Prior to launching a kernel, all the data required for the computation, must be transferred from the host (CPU) memory to the device (GPU) memory across the PCI Express bus. This operation takes a considerable time, because the PCI Express bus is managed by the operation system according to the system requirements and not the single job. If the amount of data is not great, the total time (the computation time plus the data transfer time) may be larger than the execution time of the serial program which does not require data transfer because it operate on the RAM locations directly.
Threads organization in grids and blocks is a logical structure that the user defines at the time of the call to the global function. In the next subsection IV-C we will see how to call a global function and how define the logical structure during call.
To write a program in a parallel paradigm we need to follow
some rules reported in listing 1 [
The keyword global indicates a function that runs on the device and is called from host code. Such a function is called kernel and runs in multiple instances on several blocks and threads. The number of blocks and threads is determined by using the <<<; >>> symbols, with bkdims and thdims indicate the number of blocks and threads to execute.
cudaMalloc(. . . ) is used to allocate variables in the device memory, it takes the parameter (void**) devPtr, a variable address to allocate; size is an integer requested to allocate size in bytes, generally it is used N*sizeof(. . . ) to calculate it.
cudaMemcp(. . . ) is used to copy data from host to device end vice versa; the needed parameters are: dst destination memory address; src source memory address; size the requested allocation size in byte and kind the type of transfer and is a specific token cudaMemcpyHostToDevice cudaMemcpyDeviceToHost the first copies the data through the PCI Express bus from HOST to DEVICE, instead the second from DEVICE to HOST.
cudaFree(. . . ) is used to free device memory, being limited, when the computation is ended. In this manner, we indicate the device which are the useful variables to keep in memory for other calculation and which are not.
cudaError: each cuda function returns an enumerator type,
named cudaError or simply “cuda status”. It indicates if the
operation in the GPU was successful or not. If all is well,
it returns the cuda status: cudaSuccess, in case of error the
cuda status will be a keyword describing the type of error. A
complete list of cuda errors can be found on NVIDIA site in
“related page/cudaError” [
Indexes linearization is an important issue in cuda programming, so great attention is required. Indexes linearization technique allows us to express a multidimensional index as one dimensional. Let us express the concept with some example.
1) Example of indexes linearization, from 2 indexes to 1: Assume that we have a squared sheet and y is the row index, row is the row size, x is the column index and col is the column size; the element (y; x) = (2; 5) in a linear index is determined by equation (8)
Id = y row + x (8) as shown in figure 9; note that the index starts from 0.
2) Example of indexes linearization, from n indexes to 1: This example extends the previous, i.e. assume to have a new dimension where k is the sheet index and sheet is the block size; now an element in the position (k; y; x) = (2; 2; 5) has the linear index determined by equation (9)
Id = k row col + y row + x (9)
If we have one more dimension, for example the ream index r, the linear index Id is determined by equation (10) Id = r sheet row col + k row col + y row + x (10) Fig. 9: example of indexes linearization from two to one index. Fig. 10: Example of indexes linearization from n to one index
In general, assume that I Zn with n 2 N is a discrete and limited set, so that for each dimension it exists a maximum Mi and a minimum mi:
I =< m1; M1 >< m2; M2 >
n = Y < mi; Mi >
i=1 the i-th volume is defined by < mn; Mn > so index Id is determined by equation (13)
Vi = (Mi
mi) n Id = x1 + X xj j=2 j 1 Y Vi i=1
This section explains the implementation of an efficient sequential and parallel algorithm for a median filtering. The 1D audio signal is called track in. In a sequential approach, to determinate the median from a window of N elements, first of all, the window must be sorted in ascendant or descendant order, then the median element can be taken. This element is situated in the middle of the sorted window, for this reason N must be an odd number, as shown in figure 4 with N = 5. The Serial Median Filter (SMF) is shown in Listing 2, relating to input audio track in. Then, C is calculated using equation (6) and the window of N elements is created. Moreover, the first C samples of the input track (track in) are copied from it to the output track (track out) because they are in the boundary of (11) (12) (13) the track and also the last C samples. Finally, the real Median filter is applied, and in order to filter all samples of track in the window must select from the beginning to the end of it, so that a scrollable index k is implemented for the loop (k starts from C because the first track in element starts from 0 and it ends with the samples C that represent the total number of samples in filter); therefore the window is filled with the samples as shown on the top of figure 4 and the median function is on the window. Median function sorts the window and returns the C-th element that is the median in the window. These procedures are repeated for each sample of interval [C; sample C] // serial median filter SMF void SMF(int* track_in, int* track_out) { int c = (N-1)/2; int window[N] = {0}; //copy the first c boundary elements from TRACK_IN to
TRACK_OUT for (int i = 0; i < c; i++) track_out[i] = track_in[i]; //copy the last c elements from TRACK_IN to TRACK_OUT for (int i=samples-c; i<samples; i++) track_out[i] = track_in[i]; //scroll the input track for (int k=c; k<samples-c; k++) { //fill the window for (int i=0; i<N; i++) window[i] = track_in[k-c+i]; track_out[k] = Median(window, N, c); } } // extract the median element from input array w int Median(int *w, int N, int c) { bool exchange = true; int last = N-1, i=0; while (exchange) { exchange = false; for (i = 0; i<last; i++) { if (w[i]>w[i+1]) { int tmp = w[i]; w[i] = w[i+1]; w[i+1] = tmp; exchange = true; } } } last--; } return w[c]; Listing 2: Serial algorithm used to implement the Median filter.
We have used the property described in equation (7) to determine the median. An example of application is shown in figure 11 and we can find on the top an input array, on the bottom the output array, instead, in the center three array: major (M), equal (e) and minor (m). These arrays have the same dimensions of the window and take into account the number of elements greater, equals or minor than the current element. In the example of figure 11, the number 6 has two elements much greater in the window, one element equals to it and two minor elements. In this case the number 6 has all requirements to be the median element as specified in equation (7).
The parallel algorithm (CMF—Median Filter with CUDA extensions) is shown in listing 4, that is the global function, furthermore the cuda model is shown in the listing 1 Fig. 11: Example of median filtering operation with a window of 5 elements executed in the GPU with the CMF method with memory allocations and data transfer. It is important to determine the number of threads for the simulation of the filter. Each sample must be compared to all samples of the N elements of the window to obtain the total number of threads: Ttot = samples
N (14) the maximum number of threads per core in the GPU (1024) will be indicated as nthread. So, as shown in figure 12, each block will represent the dimension bkarea = bkd N where bkd is gived by
nthread bkd = (15)
N finally,we have assumed a mono grid of dimension grd 1 1, the total block with dimension bkd N will be: gnd = Ttot + 1 = samples + 1 (16)
bkarea bkd
The global CMF function is in listing 4. While in the sequential algorithm k was a scrollable index, now in the parallel version it is a global index as in the indexes linearization, example IV-C1, therefore the k-th sample will be assigned to the k-th thread with a correspondence 1 : 1. The term C is calculated using equation (6) and takes into account the boundary condition, defining another index i. The column index is used to create the window. For each sample, arrays with major, equal and minor elements are created in order to determine the median mi in the window as shown in figure 11. Each thread executes two comparison to determine the major, minor or equal sample. The median can be found by observing the three arrays vertically. The median is present if a column complies with the specifications (7). Note that the sum of number in this column is always the window width of N. //CMF = Median Filter with CUDA extension __global__ void CMF(int* in, int* out, int c) { // k = input array index int k = c+blockDim.x*blockIdx.x+threadIdx.x; int M[N] = {0}; //major array int e[N] = {0}; //equal array int m[N] = {0}; //minor array //window index int i = k-c+threadIdx.y; for (int j=k-c; j<k-c+N; j++) { if (in[j]>in[i]) M[threadIdx.y]+=1; else if (in[j]<in[i]) m[threadIdx.y]+=1; else e[threadIdx.y]+=1; } for (int j = 0; j<N; j++) { if (M[j]==c || m[j]==c || e[j]>=c) { out[k] = in[k-c+j]; return; } } } Listing 4: Parallel algorithm used to implement the Median filter in CUDA.
For computing results, two different GPUs and an Intel CPU, in two different OSs, have been used: the NVIDIA’s GeForce 820M in a node with Windows 10, the NVIDIA’s GeForce GTX 480 in a node with Linux MIND 14.4 and the Intel Core i7-4510U in a node with Windows 10. The 820M is considered an entry level GPU, instead the GTX 480 is a midrange GPU. The main features of the three devices are given in Table I. the realization of real-time systems in image processing. These results have been obtained by filtering the audio signals as shown in figure 7 with different duration and several sample numbers. During the experiments different window sizes have been used. The work on 1D data can be easily extended to 2D data.
The characteristic of the different audio tracks are reported in Table II and Table III
Figure 13 shows the simulation results: the numbers of samples to process are represented in the x axis; the y axis represents the response time of the SMF and CMF expressed in milliseconds. Figure 13a, 13b and 13c show the obtained results varying the window width from 5 to 9 and then to 15, respectively. Data have been expressed in logarithmic scales and the results highlight the behavior of few samples during filtering. The CPU graphic has a computational complexity, when the filtering of the samples increases, the time increase exponentially. For a few samples the GPU uses more time than the CPU. The reason is the data transfer from the host to device and vice versa through the PCI Express bus, this time isn’t negligible. The PCI Express bus is managed by the OS, the data flow through the bus when it isn’t busy. A dedicated bus will make the computation more easy. In figure 7 we can see the median effect of the filter in noise reduction. The noise, clearly, can’t be totally removed but is reduced. In fact, the spectral analysis still reveals the presence of noise. The same data of figure 13 are reported in the table IV with a more clear result.
To highlight the behavior of a lot of samples during filtering, figure 13c reports a linear scale in figure 13d. The speed gain is about 2 for the 820M and 4 for the GTX 480. In figure 13e and 13f are reported the ratio between the CPU performance and GPU. The trend of the three graphs is ascendant, this demonstrates that GPU computing is more powerful than CPU computing.
To measure the execution time was used the following code. #include <time.h> float duration = 0; start = clock(); /*your codes here*/ end = clock(); //time in milliseconds duration = (float)(end-start)*1000... .../CLOCKS_PER_SEC;
Listing 5: Execution time measure method
In this section we evaluate the extension of this work to a 2D domain in the real-time scenario. The total number of samples is 100 106 s = d fs = d
Ts where s represents samples and d represents the duration. So that the samples per sec (sample rate = sr) are given by
s sr = = fs (18)
d Supposing negligible the data transfer time for many samples, the elaboration capability of the GeForce GTX 480 is evaluated by: srGT X480 = s tr = n i a G 1:5 0:5 2 1 0 4 3 2 1 0 103 104
105 106 Number of samples 107
108 (c) Case N=15 in logarithmic scale 103 104
105 106 Number of samples 107
108 (d) Case N=15 in linear scale N=5 N=9 N=15
N=5 N=9
N=15 i7-4510U
GeForce 820M GeForce GTX 480 103 104
105 106 Number of samples 107
108 (a) Case N=5 in logarithmic scale
i7-4510U
GeForce 820M GeForce GTX 480 103 104
105 106 Number of samples 107
108 (b) Case N=9 in logarithmic scale 105
i7-4510U
GeForce 820M GeForce GTX 480 s 105 d n o se 104 c i l l i m in 103 e m iT 102 e s n sp 101 o e R
100 sd 105 n o c is 104 e l l i m in 103 e m i T 102 e s n sp 101 o e R
100 n i a G 1:5 0:5 1 0 103 104
105 106 Number of samples 107 108 103 104
105 106 Number of samples 107 108 (e) Speed Gain of GeForce 820M (f) Speed Gain of GeForce GTX 480.
Fig. 13: Simulation results: Comparison between sequential algorithm and parallel algorithm in terms of response time. a) case N=5; b) case N=9; c) case N=15; e and f) Representation of Speed Gain: CP UResponseT ime=GP UResponseT ime.
1 106[samples] (20)
Applying the CMF takes about 766 milliseconds. If we have a video at 480p@24fps, it means we have the following sample rate 22:1 106h samples i sec
(21) In this case the GeForce GTX 480 can not be used for a realtime filtering process on a video at 480p@24fps so that a more powerful GPU is required.
Other similar works are “High Performance Median Filter”
[
In conclusion, the method reported by equation (7) is effective for the CUDA implementation. The simulation results demonstrate that is possible to obtain gain in response time with an entry level GPU, allowing real-time image and audio filtering. However, the bottleneck of these systems is the PCI Express bus, for dedicated and direct bus through GPU/RAM and CPU/ GPU the response time is reduced.