\newcommand{\kl}{\includegraphics[scale=0.6]{Chapters/chapter3/img/kernLeft.png}~}
\newcommand{\kr}{\includegraphics[scale=0.6]{Chapters/chapter3/img/kernRight.png}}
+
\chapter{Setting up the environnement.}
Image processing using a GPU often means using it as a general purpose computing processor, which soon brings up the issue of data transfers, 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.
8 or 16~bit-coded gray-level input images whose dimensions $H\times W$ are multiples of 512 pixels.
\begin{algorithm}
- \SetNlSty{textbf}{}{:}
+\SetNlSty{}{}{:}
allocate and populate CPU memory \textbf{h\_in}\;\\
allocate CPU pinned-memory \textbf{h\_out}\;\\
allocate GPU global memory \textbf{d\_out}\;\\
declare GPU texture reference \textbf{tex\_img\_in}\;\\
allocate GPU array in global memory \textbf{array\_img\_in}\;\\
bind GPU array \textbf{array\_img\_in} to texture \textbf{tex\_img\_in}\;\\
- copy data from \textbf{h\_in} to \textbf{array\_img\_in}\label{algo:memcopy:H2D}\; \\
+ copy data from \textbf{h\_in} to \textbf{array\_img\_in}\label{algo:memcopy:H2D}\;\\
kernel\kl gridDim,blockDim\kr()\tcc*[f]{outputs to d\_out}\label{algo:memcopy:kernel}\;\\
copy data from \textbf{d\_out} to \textbf{h\_out} \label{algo:memcopy:D2H}\;\\
\caption{Global memory management on CPU and GPU sides.}
\label{algo:memcopy}
- \end{algorithm}
+\end{algorithm}
\section{Data transfers, memory management.}
This section deals with the following issues:
\begin{enumerate}
-\item data transfer from CPU memory to GPU global memory: several GPU memory areas are available as destination memory but the 2-D caching mechanism of texture memory, specifically designed for fetching neighboring pixels, is currently the fastest way to fetch gray-level pixel values inside a kernel computation. This has lead us to choose \textbf{texture memory} as primary GPU memory area for images.
-\item data fetching from GPU global memory to kernel local memory: as said above, we use texture memory. Depending on which process is run, texture data is used either by direct fetching in kernel local memory or through a prefetching in thread block shared memory.
-\item data outputting from kernels to GPU memory: there is actually no alternative to global memory, as kernels can not directly write into texture memory and as copying from texture to CPU memory would not be faster than from simple global memory.
-\item data transfer from GPU global memory to CPU memory: it can be drastically accelerated by use of \textbf{pinned memory}, keeping in mind it has to be used sparingly.
+\item Data transfer from CPU memory to GPU global memory: several GPU memory areas are available as destination memory but the 2-D caching mechanism of texture memory, specifically designed for fetching neighboring pixels, is currently the fastest way to fetch gray-level pixel values inside a kernel computation. This has lead us to choose \textbf{texture memory} as primary GPU memory area for images.
+\item Data fetching from GPU global memory to kernel local memory: as said above, we use texture memory. Depending on which process is run, texture data is used either by direct fetching in kernel local memory or through a prefetching in shared memory.
+\item Data outputting from kernels to GPU memory: there is actually no alternative to global memory, as kernels can not directly write into texture memory and as copying from texture to CPU memory would not be faster than from simple global memory.
+\item Data transfer from GPU global memory to CPU memory: it can be drastically accelerated by use of \textbf{pinned memory}, keeping in mind it has to be used sparingly.
\end{enumerate}
Algorithm \ref{algo:memcopy} summarizes all the above considerations and describe how data are handled in our examples. For more information on how to handle the different types of GPU memory, we suggest to refer to CUDA programmer's guide.
-At debug stage, for simplicity's sake, we use the \textbf{cutil} library supplied by the NVidia developpement kit (SDK). Thus, in order to easily implement our examples, we suggest readers download download and install the latest NVidia-SDK (ours is SDK4.0), create a new directory \textit{SDK-root-dir/C/src/fast\_kernels} and adapt the generic \textit{Makefile} present in each sub-directory of \textit{SDK-root-dir/C/src/}. Then, only two more files will be enough to have a fully operational environnement: \textit{main.cu} and \textit{fast\_kernels.cu}.
+At debug stage, for simplicity's sake, we use the \textbf{cutil} library supplied by the NVidia developpement kit (SDK). Thus, in order to easily implement our examples, we suggest readers download and install the latest NVidia-SDK (ours is SDK4.0), create a new directory \textit{SDK-root-dir/C/src/fast\_kernels} and adapt the generic \textit{Makefile} that can be found in each sub-directory of \textit{SDK-root-dir/C/src/}. Then, only two more files will be enough to have a fully operational environnement: \textit{main.cu} and \textit{fast\_kernels.cu}.
Listings \ref{lst:main1}, \ref{lst:fkern1} and \ref{lst:mkfile} implement all the above considerations minimally, while remaining functional.
The main file of Listing \ref{lst:main1} is a simplified version of our actual main file.
It has to be noticed that cutil functions \texttt{cutLoadPGMi} and \texttt{cutSavePGMi} only operate on unsigned integer data. As data is coded in short integer format for performance reasons, the use of these functions involves casting data after loading and before saving. This may be overcome by use of a different library. Actually, our choice was to modify the above mentioned cutil functions.
-Listing \ref{lst:fkern1} gives a minimal kernel skeleton that will serve as the basis for all other kernels. Lines 5 and 6 determine the coordinates $(i, j)$ of the pixel to be processed. Each pixel is associated with one thread.
+Listing \ref{lst:fkern1} gives a minimal kernel skeleton that will serve as the basis for all other kernels. Lines 5 and 6 determine the coordinates $(i, j)$ of the pixel to be processed, each pixel being associated to one thread.
The instruction in line 8 combines writing the output gray-level value into global memory and fetching the input gray-level value from 2-D texture memory.
The Makefile given in Listing \ref{lst:mkfile} shows how to adapt examples given in SDK.
\section{Performance measurements}
As our goal is to design very fast implementations of basic image processing algorithms, we need to make quite accurate time-measurements, within the order of magnitude of $0.01~ms$. Again, the easiest way of doing so is to use the helper functions of the cutil library. As usual, as the durations we are measuring are short and possibly suject to non neglectable variations, a good practice is to measure multiple executions and issue the mean runtime. All time results given in this chapter have been obtained through 1000 calls to each kernel.
-Listing \ref{lst:chronos} shows how to use the dedicated cutil functions. Timer declaration and creation only need to be performed once while reset, start and stop can be used as often as necessary. Synchronization is mandatory before stopping the timer (Line 7), to avoid runtime measure being biased.
+Listing \ref{lst:chronos} shows how to use the dedicated cutil functions. Timer declaration and creation only need to be performed once while reset, start and stop can be used as often as necessary. Synchronization is mandatory before stopping the timer (Line 7), to avoid runtime measurement being biased.
\lstinputlisting[label={lst:chronos},caption=Time measurement technique using cutil functions]{Chapters/chapter3/code/exChronos.cu}
-In an attempt to provide relevant speedup values, we either implemented CPU versions of the algorithms studied, or used the values found in existing literature. Still, the large number and diversity of hardware platforms and GPU cards make it impossible to benchmark every possible combination and significant differences may occur between the speedups we announce and those obtained with different devices. As a reference, our developing platform details as follows:
+In an attempt to provide relevant speedup values, we either implemented CPU versions of the algorithms studied, or used the values found in existing literature. Still, the large number and diversity of hardware platforms and GPU cards makes it impossible to benchmark every possible combination and significant differences may occur between the speedups we announce and those obtained with different devices. As a reference, our developing platform details as follows:
\begin{itemize}
\item CPU codes run on:
\begin{itemize}
- \item Quad Core Xeon E31245 at 3.3GHz-8GByte RAM running Linux kernel 3.2
- \item Quad Core Xeon E5620 at 2.40GHz-12GByte RAM running Linux kernel 2.6.18
+ \item \textbf{Xeon}: a recent and very efficient Quad Core Xeon E31245 at 3.3GHz-8GByte RAM running Linux kernel 3.2.
\end{itemize}
\item GPU codes run on:
\begin{itemize}
- \item Nvidia Tesla C2070 hosted by a PC QuadCore Xeon E5620 at 2.4GHz-12GByte RAM, running Linux kernel 2.6.18
- \item NVidia GeForce GTX 280 hosted by a PC QuadCore Xeon X5482 at 3.20GHz-4GByte RAM, running Linux kernel 2.6.32
+ \item \textbf{C2070}: Nvidia Tesla C2070 hosted by a PC QuadCore Xeon E5620 at 2.4GHz-12GByte RAM, running Linux kernel 2.6.18
+ \item \textbf{GTX280}: NVidia GeForce GTX 280 hosted by a PC QuadCore Xeon X5482 at 3.20GHz-4GByte RAM, running Linux kernel 2.6.32
\end{itemize}
\end{itemize}
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) \cite{5402362, chen09} which allows very interesting runtimes on CUDA-enabled devices but, as far as we know, the fastest implementation to date is the histogram-based CCDS median filter \cite{6288187}.
-Some of the following implementations, feature very fast runtimes. They are 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.
+Some of the following implementations feature very fast runtimes; They are 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.
The fastest ones are based on one efficient parallel implementation of the BVM algorithm described in \cite{mcguire2008median}, improving its performance through fine tuning of its implementation.
\section{Median filtering}
\subsection{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\times n$ pixels and centered on $I(i,j)$. The output value $I'(i,j)$ is the median value of the gray level values of the $n\times n$ pixels of $H(i,j)$. Figure \ref{fig:median_1} illustrates this principle with an example of a 5x5 median filter applied on pixel $I(5,6)$. The output value is the median value of the 25 values of the dark gray window centered on pixel $I(5,6)$.
- The generic filtering method is given by Algorithm \ref{algo_median_generic}. After the data transfer stage of line \ref{algo_median_generic:memcpyH2D} which copies data from CPU memory to GPU texture memory, the actual median computing occurs between lines \ref{algo_median_generic:cptstart} and lines \ref{algo_median_generic:cptend}, before the final transfer which copies data back to CPU memory at line \ref{algo_median_generic:memcpyD2H}. Obviously, on key issue is the selection method that identifies the median value. But, as shown in figure \ref{fig:median_overlap}, since two neighboring pixels share part of the values to be sorted, a second key issue is how to rule redundancy between consecutive positions of the running window $H(i,j)$.
-As mentioned earlier, the selection of the median value can be performed by por than one technique, using either histogram-based or sorting methods, each of them having its own benefits and drawbacks as will be discussed further down.
+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\times n$ pixels and centered on $I(i,j)$. The output value $I'(i,j)$ is the median value of the gray level values of the $n\times n$ pixels of $H(i,j)$. Figure \ref{fig:median_1} illustrates this principle with an example of a 5x5 median filter applied on pixel $I(5,6)$. The output value is the median value of the 25 values of the dark gray window centered on pixel $I(5,6)$.
+Figure \ref{fig:sap_examples} shows an example of a $512\times 512$ pixel image, corrupted by a \textit{salt and pepper} noise and the denoised versions, output respectively by a $3\times 3$, a $5\times 5$ and a 2 iterations $3\times 3 $ median filter.
+\begin{figure}
+\centering
+ \subfigure[Airplane image, corrupted by salt and pepper noise of density 0.25]{\label{img:sap_example_ref} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25.png}}\qquad
+ \subfigure[Image denoised by a $3\times 3$ median filter]{\label{img:sap_example_med3} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med3.png}}\\
+ \subfigure[Image denoised by a $5\times 5$ median filter]{\label{img:sap_example_med5} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med5.png}}\qquad
+ \subfigure[Image denoised by 2 iterations of a $3\times 3$ median filter]{\label{img:sap_example_med3_it2} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_med3_it2.png}}\\
+ \caption{Exemple of median filtering, applied to salt \& pepper noise reduction.}
+ \label{fig:sap_examples}
+\end{figure}
+
+ The generic filtering method is given by Algorithm \ref{algoMedianGeneric}. After the data transfer stage of the first line, which copies data from CPU memory to GPU texture memory, the actual median computing occurs, before the final transfer which copies data back to CPU memory at the last line. Obviously, one key issue is the selection method that identifies the median value. But, as shown in figure \ref{fig:median_overlap}, since two neighboring pixels share part of the values to be sorted, a second key issue is how to rule redundancy between consecutive positions of the running window $H(i,j)$.
+As mentioned earlier, the 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 as will be discussed further down.
\subsection{A naive implementation}
As a reference, Listing \ref{lst:medianGeneric} gives a simple, not to say simplistic implementation of a CUDA kernel (\texttt{kernel\_medianR}) achieving generic $n\times n$ histogram-based median filtering. Its runtime has a very low data dependency, but this implementation does not suit very well GPU architecture. Each pixel loads the whole of its $n\times n$ neighborhood meaning that one pixel is loaded multiple times inside one single thread block, and above all, the use of a local vector (histogram[]) considerably downgrades performance, as the compiler automatically stores such vectors in local memory (slow).
Table \ref{tab:medianHisto1} displays measured runtimes of \texttt{kernel\_medianR} and pixel throughputs for each GPU version and for both CPU and GPU implementations. Usual window sizes of $3\times 3$, $5\times 5$ and $7\times 7$ are shown. Though some specific applications require larger window sizes and dedicated algorithms , such small square window sizes are most widely used in general purpose image processing. GPU runtimes have been obtained with a grid of 64-thread blocks. This block size, is a good compromise in this case.
The first observation to make when analysing results of Table \ref{tab:medianHisto1} is that, on CPU, window size has almost no influence on the effective pixel throughput.
-Since inner loops that fill the histogram vector contain very few fetching instructions (from 9 to 49, depending on the window size), it is not surprising to note neglectable runtime compared to the runtime of outer loops that fetch image pixels (from 256k to 16M instructions).
-One could be tempted to claim that CPU has no chance to win, which is not so obvious as it highly depends on what kind of algorithm is run and above all, how it is implemented. Despite a maximum effective throughput potential that is almost five times higher, measured GTX280 throughput values sometimes prove slower than CPU values, as shown in Table \ref{tab:medianHisto1}.
+Since inner loops that fill the histogram vector contain very few fetching instructions (from 9 to 49, depending on the window size), it is not surprising to note their neglectable impact compared to outer loops that fetch image pixels (from 256k to 16M instructions).
+One could be tempted to claim that CPU has no chance to win, which is not so obvious as it highly depends on what kind of algorithm is run and above all, how it is implemented. To illustrate this, we can notice that, despite a maximum effective throughput potential that is almost five times higher, measured GTX280 throughput values sometimes prove slower than CPU values, as shown in Table \ref{tab:medianHisto1}.
-On the GPU's side, we note high dependence on window size due to the redundancy induced by the multiple fetches of each pixel inside each block, becoming higher with the window size as illustrated by Figure \ref{fig:median_overlap}. On C2070 card, thanks to a more efficient caching mechanism, this effect is lesser. On GPUs, dependency over image size is low, due to slightly more efficient data transfers when copying larger data amounts. Thus transferring a 4096$\times$4096 pixel image (32~MBytes) is a bit faster than transferring 64 times a 512$\times$512 pixel image (0.5~MBytes).
+On the GPU's side, we note high dependence on window size due to the redundancy induced by the multiple fetches of each pixel inside each block, becoming higher with the window size as illustrated by Figure \ref{fig:median_overlap}. On C2070 card, thanks to a more efficient caching mechanism, this effect is lesser. On GPUs, dependency over image size is low, and due to slightly more efficient data transfers when copying larger data amounts, pixel throughputs increases with image size. As an example, transferring a 4096$\times$4096 pixel image (32~MBytes) is a bit faster than transferring 64 times a 512$\times$512 pixel image (0.5~MBytes).
%% mettre l'eau à la bouche
\end{figure}
\begin{algorithm}
- \SetNlSty{textbf}{}{:}
- copy data from CPU to GPU texture memory\label{algo_median_generic:memcpyH2D}\;
+ %\SetNlSty{}{}{:}
+ % \SetLine
+ %\linesnumbered
+ copy data from CPU to GPU texture memory\label{algoMedianGeneric:memcpyH2D}\;\\
\ForEach(\tcc*[f]{in parallel}){pixel at position $(x, y)$}{
- Read gray-level values of the n$\times$n neighborhood\label{algo_median_generic:cptstart}\;
- Selects the median ($(n^2/2)^{th}$) value among those n$\times$n values\;
- Outputs the new gray-level value \label{algo_median_generic:cptend}\;
+ Read gray-level values of the n$\times$n neighborhood\label{algoMedianGeneric:cptstart}\;\\
+ Selects the median value among those n$\times$n values\;\\
+ Outputs the new gray-level value \label{algoMedianGeneric:cptend}\;\\
}
-copy data from GPU global memory to CPU memory\label{algo_median_generic:memcpyD2H}\;
-\caption{generic n$\times$n median filter}
-\label{algo_median_generic}
+copy data from GPU global memory to CPU memory\label{algoMedianGeneric:memcpyD2H}\;\\
+\caption{\label{algoMedianGeneric}generic n$\times$n median filter}
\end{algorithm}
\begin{figure}
\hline
\multicolumn{2}{|l||}{Processor} & \multicolumn{3}{c|}{\textbf{GTX280}} & \multicolumn{3}{c|}{\textbf{C2070}} & \multicolumn{3}{c|}{\textbf{Xeon}} \\ \hline
\multicolumn{2}{|l||}{\shortstack{Performances$\rightarrow$\\sizes (pixels)$\downarrow$}} & \shortstack{t\\(ms)}& \shortstack{output\\(MP/s)}& \shortstack{rate\\\% }&\shortstack{t\\(ms)}& \shortstack{output\\(MP/s)}& \shortstack{rate\\\% }&\shortstack{t\\(ms)}& \shortstack{output\\(MP/s)}& \shortstack{rate\\\% } \\ \whline
-\multirow{3}{*}{\rotatebox{90}{512$^2$}} &3$\times$3&11.50 &22 &2.2 &7.58 &33 &3.4 & 19.25& 14&-\\
- &5$\times$5&19.10 &14 &1.3 &8.60 &30 &3.0 &18.49 &14 &-\\
- &7$\times$7&31.30 &8 &0.8 &10.60 &24 &2.5 &20.27 &13 &-\\\whline
-\multirow{3}{*}{\rotatebox{90}{1024$^2$}}&3$\times$3&44.50 &23 &2.3 &29.60 &34 &3.5 &75.49 &14 &-\\
- &5$\times$5&71.10 &14 &1.4 &33.00 &31 &3.2 &73.88 &14 &-\\
- &7$\times$7&114.50 &9 &0.9 &39.10 &26 &2.7 &77.40 &13 &-\\\whline
-\multirow{3}{*}{\rotatebox{90}{2048$^2$}}&3$\times$3&166.00 &24 &2.4 &115.20 &36 &3.6 &296.18&14 &-\\
- &5$\times$5&261.00&16 &1.5 &128.20&32 &3.3 &294.55&14 &-\\
- &7$\times$7&411.90 &10&1.0 &143.30&28 &2.8 &303.48&14 &-\\\whline
-\multirow{3}{*}{\rotatebox{90}{4096$^2$}}&3$\times$3&523.80 &31 &3.0 &435.00 &38 &3.9 &1184.16&14 &-\\
- &5$\times$5&654.10&25 &2.4 &460.20&36 &3.7 &1158.26&14 &-\\
- &7$\times$7&951.30 &17&1.7 &509.60&32 &3.3 &1213.55&14 &-\\\whline
+\multirow{3}{*}{\rotatebox{90}{512$^2$}} &3$\times$3&11.50 &22 &2.2 &7.58 &33 &3.4 & 19.25& 14&11\\
+ &5$\times$5&19.10 &14 &1.3 &8.60 &30 &3.0 &18.49 &14 &11\\
+ &7$\times$7&31.30 &8 &0.8 &10.60 &24 &2.5 &20.27 &13 &10\\\whline
+\multirow{3}{*}{\rotatebox{90}{1024$^2$}}&3$\times$3&44.50 &23 &2.3 &29.60 &34 &3.5 &75.49 &14 &11\\
+ &5$\times$5&71.10 &14 &1.4 &33.00 &31 &3.2 &73.88 &14 &11\\
+ &7$\times$7&114.50 &9 &0.9 &39.10 &26 &2.7 &77.40 &13 &10\\\whline
+\multirow{3}{*}{\rotatebox{90}{2048$^2$}}&3$\times$3&166.00 &24 &2.4 &115.20 &36 &3.6 &296.18&14 &11\\
+ &5$\times$5&261.00&16 &1.5 &128.20&32 &3.3 &294.55&14 &11\\
+ &7$\times$7&411.90 &10&1.0 &143.30&28 &2.8 &303.48&14 &11\\\whline
+\multirow{3}{*}{\rotatebox{90}{4096$^2$}}&3$\times$3&523.80 &31 &3.0 &435.00 &38 &3.9 &1184.16&14 &11\\
+ &5$\times$5&654.10&25 &2.4 &460.20&36 &3.7 &1158.26&14 &11\\
+ &7$\times$7&951.30 &17&1.7 &509.60&32 &3.3 &1213.55&14 &11\\\whline
\end{tabular}}
\caption{Performance results of \texttt{kernel medianR}. }
\end{table}
\section{NVidia GPU tuning recipes}
-When designing GPU code, besides thinking of the actual data computing process, one must choose the memory type into which to store temporary data. Three type of GPU memory are available:
+When designing GPU code, besides thinking of the actual data computing process, one must choose the memory type into which to store temporary data. Three types of GPU memory are available:
\begin{enumerate}
\item \textbf{Global memory, the most versatile:}\\Offers the largest storing space and global scope but is slowest (400 cycles latency). \textbf{Texture memory} is physically included into it, but allows access through an efficient 2-D caching mechanism.
-\item \textbf{Registers, the fastest:}\\Allows access wihtout latency, but only 63 registers are available per thread (thread scope), with a maximum of 32K per Symetric Multiprocessor (SM).
+\item \textbf{Registers, the fastest:}\\Allow access wihtout latency, but only 63 registers are available per thread (thread scope), with a maximum of 32K per Symetric Multiprocessor (SM).
\item \textbf{Shared memory, a complex compromise:}\\All threads in one block can access 48~KBytes of shared memory, which is faster than global memory (20 cycles latency) but slower than registers.
-However, bank conflicts can occur if two threads of a warp try to access data stored in one single memory bank. In such cases, the parallel process is re-serialized which may cause significant performance decrease. One easy way to avoid it is to ensure that two consecutive threads in one block always access 32 bit data at two consecutive adresses.
+However, bank conflicts can occur if two threads of a warp try to access data stored in one single memory bank. In such cases, the parallel process is re-serialized which may cause significant performance decrease. One easy way to avoid it is to ensure that two consecutive threads in one block always access 32-bit data at two consecutive adresses.
\end{enumerate}
\noindent As observed earlier, designing a median filter GPU implementation using only global memory is fairly straightforward, but its performance remains quite low even if it is faster than CPU.
-To overcome this, the most frequent choice made in efficient implementations found in literature is to use shared memory. Such option implies prefetching data prior to doing the actual computations, a relevant choice, as each pixel of an image belongs to n$\times$n different neighborhoods. Thus, it can be expected that fetching each gray-level value from global memory only once should be more efficient than do it each time it is required. One of the most efficient implementations using shared memory is presented in \cite{5402362}. In the case of the generic kernel of Listing \ref{lst:medianGeneric}, using shared memory without further optimization would not bring valuable speedup because that would just move redundancy from texture to shared memory fetching and would generate bank conflicts. For information, we wrote such a version of the generic median kernel and our measurements showed a speedup of around 3\% (for example: 32ms for 5$\times$5 median on a 1024$^2$ pixel image).
+To overcome this, the most frequent choice made in efficient implementations found in literature is to use shared memory. Such option implies prefetching data prior to doing the actual computations, a relevant choice, as each pixel of an image belongs to n$\times$n different neighborhoods. Thus, it can be expected that fetching each gray-level value from global memory only once should be more efficient than doing it each time it is required. One of the most efficient implementations using shared memory is presented in \cite{5402362}. In the case of the generic kernel of Listing \ref{lst:medianGeneric}, using shared memory without further optimization would not bring valuable speedup because that would just move redundancy from texture to shared memory fetching and would generate bank conflicts. For information, we wrote such a version of the generic median kernel and our measurements showed a speedup of around 3\% (as an example: 32ms for 5$\times$5 median on a 1024$^2$ pixel image, i.e. 33~Mpixel/s ).
As for registers, designing a generic median filter that would only use that type of memory seems difficult, due to the above mentioned 63 register-per-thread limitation.
Yet, nothing forbids us to design fixed-size filters, each of them specific to one of the most popular window sizes. It might be worth the effort as dramatic increase in performance could be expected.
\section{A 3$\times$3 median filter: using registers }
Designing a median filter dedicated to the smallest possible square window size is a good challenge to start using registers.
-One first issue is that the exclusive use of registers forbids us to implement a naive histogram-based method. In a \textit{8-bit gray level pixel per thread} rule, each histogram requires one 256-element vector to store its values, e.g. four times the maximum register count allowed per thread (63). Considering a 3$\times$3 median filter involves only 9 pixel values per thread, it seem obvious they can be sorted within the 63-register limit.
+One first issue is that the exclusive use of registers forbids us to implement a naive histogram-based method. In a \textit{8-bit gray level pixel per thread} rule, each histogram requires one 256-element vector to store its values, i.e. four times the maximum register count allowed per thread (63). Considering that a 3$\times$3 median filter involves only 9 pixel values per thread, it seem obvious they can be sorted within the 63-register limit.
\subsection{The simplest way}
-In the case of a 3$\times$3 median filter, the simplest solution consists in associating one register to each gray-level value, then sorting those 9 values and selecting the fifth one, e.g. the median value. For such a small amount of data to sort, a simple selection method is well indicated. As shown in Listing \ref{lst:kernelMedian3RegTri9} (\texttt{kernelMedian3RegTri9()}), the constraint of only using registers leads to adopt an unusual manner of coding. However, results are persuasive: runtimes are divided by around 120 on GTX280 and 80 on C2070, while only reduced by a 3.5 factor on CPU.
-The diagram of Figure \ref{fig:compMedians1} summarizes these first results. Only C2070 throughputs are shown and compared to CPU results. We included the maximum effective pixel throughput in order to see the improvement potential of the different implementations. We also introduced throughputd achieved by \textit{libJacket}, a commercial implementation, as it was the fastest known implementation of 3$\times$3 median filter to date, as illustrated in \cite{chen09}. One of the authors of libJacket kindly posted the CUDA code of its 3$\times$3 median filter, that we inserted into our own coding structure. The algorithm itself is quite similar to ours, but running it in our own environement produced higher throughput values than those published in \cite{chen09}, not due to different hardware capabilities between our GTX280 and the GTX260 used in the paper, but to the way we perform memory transfers and to our register-only method of storing temporary data.
+In the case of a 3$\times$3 median filter, the simplest solution consists in associating one register to each gray-level value, then sorting those 9 values and selecting the fifth one, i.e. the median value. For such a small amount of data to sort, a simple selection method is well indicated. As shown in Listing \ref{lst:kernelMedian3RegTri9} (\texttt{kernelMedian3RegTri9()}), the constraint of only using registers leads to adopt an unusual manner of coding. However, results are persuasive: runtimes are divided by around 120 on GTX280 and 80 on C2070, while only reduced by a 3.5 factor on CPU.
+The diagram of Figure \ref{fig:compMedians1} summarizes these first results. Only C2070 throughputs are shown and compared to CPU results. We included the maximum effective pixel throughput in order to see the improvement potential of the different implementations. We also introduced throughputd achieved by \textit{libJacket}, a commercial implementation, as it was the fastest known implementation of a 3$\times$3 median filter to date, as illustrated in \cite{chen09}. One of the authors of libJacket kindly posted the CUDA code of its 3$\times$3 median filter, that we inserted into our own coding structure. The algorithm itself is quite similar to ours, but running it in our own environement produced higher throughput values than those published in \cite{chen09}, not due to different hardware capabilities between our GTX280 and the GTX260 used in the paper, but to the way we perform memory transfers and to our register-only method of storing temporary data.
\lstinputlisting[label={lst:kernelMedian3RegTri9},caption= 3$\times$3 median filter kernel using one register per neighborhood pixel and bubble sort]{Chapters/chapter3/code/kernMedianRegTri9.cu}
\begin{figure}
\centering
- \includegraphics[width=11cm]{Chapters/chapter3/img/debitPlot1.png}
+ \includegraphics[width=11cm]{Chapters/chapter3/img/debitPlot1.pdf}
\caption{Comparison of pixel throughputs on GPU C2070 and CPU for generic median, in 3$\times$3 median register-only and \textit{libJacket}.}
\label{fig:compMedians1}
\end{figure}
\subsubsection{Reducing register count}
Our current kernel (\texttt{kernelMedian3RegTri9}) uses one register per gray-level value, which amounts to 9 registers for the entire 3$\times$3 window.
-This count can be reduced by use of an iterative sorting process called \textit{forgetful selection}, where both \textit{extrema} are eliminated at each sorting stage, until only 3 elements remain. The question is to find out the minimal register count $k_n$ that allows the selection of the median amoung $n^2$ values. The answer can be evaluated considering that, when eliminating the maximum and the minimum values, one has to make sure not to eliminate the global median value. Such a situation is illustrated in Figure \ref{fig:forgetful_selection} for a 3$\times$3 median filter. For better comprehension, the 9 elements of the 3$\times$3 pixel window have been represented in a row.
+This count can be reduced by use of an iterative sorting process called \textit{forgetful selection}, where both \textit{extrema} are eliminated at each sorting stage, until only 3 elements remain. The question is to find out the minimal register count $k_{n^2}$ that allows the selection of the median amoung $n^2$ values. The answer can be evaluated considering that, when eliminating the maximum and the minimum values, one has to make sure not to eliminate the global median value. Such a situation is illustrated in Figure \ref{fig:forgetful_selection} for a 3$\times$3 median filter. For better comprehension, the 9 elements of the 3$\times$3 pixel window have been represented in a row.
\begin{figure}
\centering
\includegraphics[width=6cm]{Chapters/chapter3/img/forgetful_selection.png}
\caption{Forgetful selection with the minimal element register count. Illustration for 3$\times$3 pixel window represented in a row and supposed sorted.}
\label{fig:forgetful_selection}
\end{figure}
-We must remember that, in the fully sorted vector, the median value will have the middle index e.g. $\lfloor n^2/2\rfloor$.
+We must remember that, in the fully sorted vector, the median value will have the middle index i.e. $\lfloor n^2/2\rfloor$.
Moreover, assuming that both \textit{extrema} are eliminated from the first $k$ elements and that the global median is one of them would mean that:
\begin{itemize}
\item if the global median was the minimum among the $k$ elements, then at least $k-1$ elements would have a higher index. Considering the above median definition, at least $k-1$ elements should also have a lower index in the entire vector.
Therefore, the number $k$ of elements that are part of the first selection stage can be defined by the condition
$$n^2-k \leq \lfloor \frac{n^2}{2} \rfloor -1$$
which leads to:
-$$k_n=\lceil \frac{n^2}{2}\rceil+1 $$
-This rule can be applied to the first eliminating stage and remains true with the next ones as each stage suppresses exactly two values.
+$$k_{n^2}=\lceil \frac{n^2}{2}\rceil+1 $$
+This rule can be applied to the first eliminating stage and remains true with the next ones as each stage suppresses exactly two values, one above and one below the median value.
In our 3$\times$3 pixel window example, the minimum register count becomes $k_9=\lceil 9/2\rceil+1 = 6$.
The \textit{forgetful selection} method, used in \cite{mcguire2008median} does not imply full sorting of values, but only selecting minimum and maximum values, which, at the price of a few iteration steps ($n^2-k$), reduces arithmetic complexity.
\subsubsection{More data output per thread}
In the case of a kernel achieving an effective memory throughput value far from the GPU peak value, and if enough threads are run, another technique may help hiding memory latency and thus leverage performance: make sure that each thread generates multiple pixel outputs.\\
Attentive readers could remark that it would increase the register count per thread, which can be compensated by dividing thread block size accordingly, thus allowing to keep the same register count per block.
-Moreover, it is now possible to take advantage of window overlapping, first illustrated In Figure \ref{fig:median_overlap}, and further detailed in Figure \ref{fig:median3_overlap}. As the selection is first processed on the first 6 gray-level values, e.g. exactly the number of pixels that overlap between the neighborhoods of two adjacent center pixels, 6 texture fetches and one \texttt{minmax6} selection per thread can be saved. There again, some speedup can be expected through our modified kernel source code presented in Listing \ref{lst:medianForget2pix3}. One important difference lies in the way pixel coordinates are computed from thread indexes. As each thread has to process two pixels, the number of threads in each block is divided by 2, while the grid size remains unchanged. Consequently, in our kernel code, each thread whose block-related coordinates are $(tx, ty)$ will be in charge of processing pixels of block-related coordinates $(2tx, ty)$ and $(2tx+1, ty)$; lines 5 and 6 implement this.
+Moreover, it is now possible to take advantage of window overlapping, first illustrated In Figure \ref{fig:median_overlap}, and further detailed in Figure \ref{fig:median3_overlap}. As the selection is first processed on the first 6 gray-level values, i.e. exactly the number of pixels that overlap between the neighborhoods of two adjacent center pixels, 6 texture fetches and one \texttt{minmax6} selection per thread can be saved. There again, some speedup can be expected through our modified kernel source code presented in Listing \ref{lst:medianForget2pix3}. One important difference with previous versions lies in the way pixel coordinates are computed from thread indexes. As each thread has to process two pixels, the number of threads in each block is divided by 2, while the grid size remains unchanged. Consequently, in our kernel code, each thread whose block-related coordinates are $(tx, ty)$ will be in charge of processing pixels of block-related coordinates $(2tx, ty)$ and $(2tx+1, ty)$; lines 5 and 6 implement this.
\begin{figure}
\centering
\begin{figure}
\centering
- \includegraphics[width=11cm]{Chapters/chapter3/img/debitPlot2.png}
+ \includegraphics[width=11cm]{Chapters/chapter3/img/debitPlot2.pdf}
\caption{Comparison of pixel throughput on GPU C2070 for the different 3$\times$3 median kernels.}
\label{fig:compMedians2}
\end{figure}
The next two sections will first detail the particular case of the 5$\times$5 median through register-only method and eventually a generic kernel for larger window sizes.
\subsection{A register-only 5$\times$5 median filter \label{sec:median5}}
-The minimum register count required to apply the forgetful selection method to a 5$\times$5 median filter is $k_{25}=\lceil 25/2\rceil+1 = 14$. Moreover, two adjacent overlapping windows share 20 pixels ($n^2-one\_column$) so that, when processing 2 pixels simultaneously, a count of 6 common selection stages can be carried out from the first selection stage with 14 common values to the processing of the last common value. That allows to limit register count to 14+8=22 per thread. Figure \ref{fig:median5overlap} describes the distribution of overlapping pixels, implemented in Listing \ref{lst:medianForget2pix5}: common selection stages take place from line 25 to line 37, while the remaining separate selection stages occur between lines 45 and 62 after the separation of line 40.
+The minimum register count required to apply the forgetful selection method to a 5$\times$5 median filter is $k_{25}=\lceil 25/2\rceil+1 = 14$. Moreover, two adjacent overlapping windows share 20 pixels ($n^2-one\_column$) so that, when processing 2 pixels simultaneously, a count of 7 common selection stages can be carried out from the first selection stage with 14 common values to the processing of the last common value. That allows to limit register count to 22 per thread. Figure \ref{fig:median5overlap} describes the distribution of overlapping pixels, implemented in Listing \ref{lst:medianForget2pix5}: common selection stages take place from line 25 to line 37, while the remaining separate selection stages occur between lines 45 and 62 after the separation of line 40.
\begin{figure}
\centering
\includegraphics[width=6cm]{Chapters/chapter3/img/median5_overlap.png}
- \caption{Reducing register count in a 5$\times$5 register-only median kernel outputting 2 pixels simultaneously. The first 6 forgetful selection stages are common to both processed center pixels. Only the last 5 selections have to be done separately.}
+ \caption{Reducing register count in a 5$\times$5 register-only median kernel outputting 2 pixels simultaneously. The first 7 forgetful selection stages are common to both processed center pixels. Only the last 5 selections have to be done separately.}
\label{fig:median5overlap}
\end{figure}
\end{table}
\subsection{Fast approximated n$\times$n median filter }
-Large window median filters are less widespread and used in more specific fields, such as digital microscopy: for example in \cite{paper_bio_background}, an estimation of the background gray-level intensity is achieved through a $111\times 111$ median filtering (processed in around 2s for a 4MPixel image). In such cases, a possible technique is to split median selection into two separate 1-D stages: one in the vertical direction and the other in the horizontal direction. Image processing specialists may object that this method does not select the actual median value. This is true but, in the case of large window sizes and \textit{real-life} images, the so selected value is statistically near the actual median value and often represents an acceptable approximation. Such a filter is sometimes called \textit{smoother}.
+Large window median filters are less widespread and used in more specific fields, such as digital microscopy where, for example, background estimation of images is achieved through $64\times 64$ or $128\times 128$ median filters \cite{Wu2010}. In such cases, a possible technique is to split median selection into two separate 1-D stages: one in the vertical direction and the other in the horizontal direction. Image processing specialists may object that this method does not select the actual median value. This is true but, in the case of large window sizes and \textit{real-life} images, the so selected value is statistically near the actual median value and often represents an acceptable approximation. Such a filter is sometimes called \textit{smoother}.
As explained earlier in this section, the use of large window median filters rules out register-only implementation,
which suggests to privilege the use of shared memory. The 1-D operation almost completely avoids bank conflicts in shared memory accesses.
Furthermore, the above-described forgetful selection method cannot be used anymore, as too many registers would be required.\\Instead, the Torben Morgensen sorting algorithm is used, as its required register count is both low and constant, and avoids the use of a local vector, unlike histogram-based methods.
Listing \ref{lst:medianSeparable} presents a kernel code that implements the above considerations and achieves a 1-D vertical $n \times 1$ median filter. The shared memory vector is declared as \texttt{extern} (Line 16) as its size is determined at runtime and passed to the kernel call as an argument. Lines 20 to 29 perform data prefetching, including the $2n$-row halo ($n$ at the bottom and $n$ at the top of each block). Then one synchronization barrier is mandatory (line 31) to ensure that all needed data is ready prior to its use by the different threads.
-Torben Morgensen sorting takes place between lines 37 and 71 and eventually, the transposed output value is stored in global memory at line 73. Outputting the transposed image in global memory saves time and allows to re-use the same kernel to achieve the second step, e.g 1-D horizontal $n \times 1$ median filtering. The final transpose is done at transfer time, when copying data from GPU to CPU memory, which once more saves time while actually generating the expected image.
-It has to be noticed that this smoother, unlike the technique we proposed for fixed-size median filters, can not be considered as a state-of-the-art technique, as for example the one presented in \cite{4287006}. However, it may be considered as a good, easy to use and efficient alternative as confirmed by the results presented in Table \ref{tab:medianSeparable}. Pixel throughput values achieved by our kernel, though not constant with window size, remain very competitive if window size is kept under $120\times 120$ pixels (in \cite{4287006}, pixel throughput is around 7MP/s).
+Torben Morgensen sorting takes place between lines 37 and 71 and eventually, the transposed output value is stored in global memory at line 73. Outputting the transposed image in global memory saves time and allows to re-use the same kernel to achieve the second step, e.g 1-D horizontal $n \times 1$ median filtering. The final transpose is done at transfer time, when copying data from GPU to CPU memory, which once more saves time while actually generates the expected image.
+It has to be noticed that this smoother, unlike the technique we proposed for fixed-size median filters, can not be considered as a state-of-the-art technique, as for example the one presented in \cite{4287006}. However, it may be considered as a good, easy to use and efficient alternative as confirmed by the results presented in Table \ref{tab:medianSeparable}. Pixel throughput values achieved by our kernel, though not constant with window size, remain very competitive if window size is kept under $120\times 120$ pixels, especially when outputting 2 pixels per thread (in \cite{4287006}, pixel throughput is around 7MP/s).
+Figure \ref{fig:sap_examples2} shows an example of a $512\times 512$ pixel image, corrupted by a \textit{salt and pepper} noise and the denoised versions, output respectively by a $3\times 3$, a $5\times 5$ and a $55\times 55 $ separable smoother.
+\begin{figure}
+ \subfigure[Airplane image, corrupted with by salt and pepper noise of density 0.25]{\label{img:sap_example_ref} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25.png}}\qquad
+ \subfigure[Image denoised by a $3\times 3$ separable smoother]{\label{img:sap_example_sep_med3} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med3.png}}\\
+ \subfigure[Image denoised by a $5\times 5$ separable smoother]{\label{img:sap_example_sep_med5} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med5.png}}\qquad
+ \subfigure[Image background estimation by a $55\times 55$ separable smoother]{\label{img:sap_example_sep_med3_it2} \includegraphics[width=5cm]{Chapters/chapter3/img/airplane_sap25_sep_med111.png}}\\
+ \caption{Exemple of separable median filtering (smoother), applied to salt \& pepper noise reduction.}
+ \label{fig:sap_examples2}
+\end{figure}
+
\begin{table}[h]
%\newlength\savedwidth
\newcommand\whline{\noalign{\global\savedwidth
\lstinputlisting[label={lst:medianSeparable},caption= generic pseudo median kernel.]{Chapters/chapter3/code/kernMedianSeparable.cu}
-\section{Glossary}
-\begin{Glossary}
-\item[CUDA] Compute Unified Device Architecture.
-\end{Glossary}
+% \section{Glossary}
+% \begin{Glossary}
+% \item[CUDA] Compute Unified Device Architecture.
+% \end{Glossary}
\putbib[Chapters/chapter3/biblio3]
