add ch4

[book_gpu.git] / BookGPU / Chapters / chapter4 / ch4.tex

\chapterauthor{Gilles Perrot}{FEMTO-ST Institute}

%\newcommand{\kl}{\includegraphics[scale=0.6]{Chapters/chapter4/img/kernLeft.png}~}
%\newcommand{\kr}{\includegraphics[scale=0.6]{Chapters/chapter4/img/kernRight.png}}

%% \lstset{
%%   language=C,
%%   columns=fixed,
%%   basicstyle=\footnotesize\ttfamily,
%%   numbers=left,
%%   firstnumber=1,
%%   numberstyle=\tiny,
%%   stepnumber=5,             
%%   numbersep=5pt,              
%%   tabsize=3,                  
%%   extendedchars=true,         
%%   breaklines=true,       
%%   keywordstyle=\textbf,
%%   frame=single,         
%%   % keywordstyle=[1]\textbf,   
%%   %identifierstyle=\textbf,
%%   commentstyle=\color{white}\textbf,
%%   stringstyle=\color{white}\ttfamily,
%%   % xleftmargin=17pt,
%%   % framexleftmargin=17pt,
%%   % framexrightmargin=5pt,
%%   % framexbottommargin=4pt,
%%   backgroundcolor=\color{lightgray},
%%   }


\chapter{Implementing an efficient convolution \index{Convolution} operation on GPU}
\section{Overview}
In this chapter, after dealing with GPU median filter implementations,
we propose to explore how convolutions can be implemented on modern
GPUs. Widely used in digital image processing filters, the \emph{convolution
operation} basically consists in taking the sum of products of elements
from two 2-D functions, letting one of the two functions move over
every element of the other, producing a third function that is typically
viewed as a modified version of one of the original functions. To
begin with, we shall examine non-separable or generic convolutions,
before adressing the matter of separable convolutions. We shall refer
to $I$ as an H x L pixel gray-level image, and to $I(x,y)$ as the gray-level
value of each pixel of coordinates $(x,y)$.
%dire qqs mots sur le filtrage IIR/FIR ?


\section{Definition}
Within a digital image $I$, the convolution operation is performed between
image $I$ and convolution mask \emph{h} (To avoid confusion with other
GPU functions referred to as kernels, we shall use\emph{ convolution
mask} instead of \emph{convolution kernel}) is defined by:
\begin{equation}
I'(x, y) = \left(I * h\right) = \sum_{i < H} \sum_{j < L}I(x-j, y-j).h(j,i)
\label{convoDef}
\end{equation}
While processing an image, function \emph{h} is bounded by a square
window of size \emph{k = 2r + 1}, i.e an uneven number, to ensure
there is a center. We shall also point out that, as stated earlier,
the square shape is no limiting factor to the process, as any shape
can be inscribed into a square. In the case of a more complex shape,
the remaining space is filled by null values (padding).


\section{Implementation}
The basic principle of computing a convolution between one $I$ picture
and one \emph{h} convolution mask defined on domain $\Omega$ is given
by algorithm \ref{algo_genconv} and illustrated by Figure \ref{fig:convoPrinciple}, which mainly shows how gray-level values of the center pixel's neighborhood are combined with the convolution mask values to compute the output value.  
For more readability, only part of the connecting lines are shown.
 \begin{figure}
\centering
   \includegraphics[width=11cm]{Chapters/chapter4/img/convo1.png}
   \caption{Principle of a generic convolution implementation. The center pixel is represented with a black background and the pixels of its neighborhood are denoted $I_{p,q}$ where $(p,q)$ is the relative position of the neighbor pixel. Elements $h_{t,u}$ are the values of the convolution mask.}
   \label{fig:convoPrinciple}
\end{figure}
\begin{algorithm}
\caption{generic convolution}   
\label{algo_genconv}
  \ForEach{pixel at position $(x, y)$}{
    Read all gray-level values $I(x, y)$ in the neighborhood\\\;
    Compute the weighted sum \( I_\Omega = \sum_{(j,i) \in \Omega}I(x-j, y-j).h(j,i) \)\\\;
    Normalize $I'(x, y)$ value\\\;
    Outputs the new gray-level value 
  }
\end{algorithm}

The gray-level value of each pixel of output image $I'$ is the weighted
sum of pixels included in the neighborhood defined by $\Omega$ around
the corresponding pixel in the input image. It has to be noted that,
in case the sum $S$ of all coefficients in the mask is not 1, the original
brightness of the image will be altered and a normalization stage
has to take place, as, for example, in the case of an 8-bit coded
image:
\begin{enumerate}
\item if $S \ge 0$ then $I' = I_\Omega / S$
\item if $S = 0$ then $I' = I_\Omega + 128$
\item if $S < 0$ then $I' = I_\Omega + 255$
\end{enumerate}
In case one, normalizing means performing a division operation for
each pixel, which will be quite time-costly when performed on a GPU. A simple work-around is to normalize mask values before using them in GPU kernels.


\subsection{First test implementation}
This first implementation consists of a rather naive application to
convolutions of the tuning recipes applied to median filters in the
previous chapter, as a reminder : texture memory used with incoming
data, pinned memory with output data, optimized use of registers
while processing data and multiple output per thread\index{Multiple output per thread}. 
One signifcant difference lies in the fact
that the median filter uses only one parameter, the size of the window mask,
which can be hard-coded, while a convolution mask requires referring to several; hard-coding
its elements would lead to severe lack of flexibility (one function
per filter, no external settings) so we will just use it as a starting
point in our approach. 

Let us assume that we are planning to implement the convolution defined by the following $3\times 3$ mask (low-pass filter or averaging filter):
$$h=\begin{bmatrix}1&1&1\\1&1&1\\1&1&1\end{bmatrix}$$ 
The kernel code presented in Listing \ref{lst:convoGene3Reg8} implements the convolution operation and applies all above optimizations except, for clarity reasons, multiple output per thread.
In the particular case of a generic convolution, it is important to note how mask coefficients are applied to image pixels in order to fit the definition of equation \ref{convoDef}: if the coordinates of the center pixel had been set to (0,0), then the pixel of coordinates $(i,j)$ would have been multiplied by the element $(-i,-j)$ of the mask, which, transposed in our kernel code, leads to multiply the $p^{th}$ pixel of the window by the $(n-p)^{th}$ element of the convolution mask.

\lstinputlisting[label={lst:convoGene3Reg8},caption=Generic CUDA kernel achieving a convolution operation with hard-coded mask values]{Chapters/chapter4/code/convoGene3Reg8.cu}

Table \ref{tab:convoNonSepReg1} shows kernel timings and throughput values for such a low-pass filter extended to $5\times 5$ and $7\times 7$ masks applied on 8-bit coded gray-level
images of sizes $512\times 512$, $1024\times 1024$, $2048\times 2048$, $4096\times 4096$ and run on a C2070 card with $32\times 8$ thread blocks.
As a reminder, Table \ref{tab:memcpy1} details the data transfer costs that helped computing throughput values.


\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\textbf{Mask size$\rightarrow$}&\multicolumn{2}{|c|}{\textbf{3x3}}&\multicolumn{2}{|c|}{\textbf{5x5}}&\multicolumn{2}{|c|}{\textbf{7x7}}\\
\textbf{Image size$\downarrow$}&time (ms)&TP&time (ms)&TP&time (ms)&TP\\\hline\hline
$\mathbf{512\times 512}$  &0.077&1165 &0.209&559  &0.407   &472 \\\hline
$\mathbf{1024\times 1024}$&0.297&1432 &0.820&836  &1.603   &515 \\\hline
$\mathbf{2048\times 2048}$&1.178&1549 &\bf 3.265&\bf 875 &6.398&529 \\\hline
$\mathbf{4096\times 4096}$&4.700&1585 &13.05&533     &25.56&533 \\\hline
\end{tabular}
}  
\caption{Timings ($time$) and throughput values ($TP$ in Mpix/s) of one register-only non separable convolution kernel, for small mask sizes of $3\times 3$, $5\times 5$ and $7\times 7$ pixels, on a C2070 card (fermi architecture). Data transfer duration are those of Table \ref{tab:memcpy1}.}
\label{tab:convoNonSepReg1}
\end{table} 

\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\textbf{Mask size$\rightarrow$}&\multicolumn{2}{|c|}{\textbf{3x3}}&\multicolumn{2}{|c|}{\textbf{5x5}}&\multicolumn{2}{|c|}{\textbf{7x7}}\\
\textbf{Image size$\downarrow$}&time (ms)&TP&time (ms)&TP&time(ms)&TP\\\hline\hline
$\mathbf{512\times 512}$  &0.060&1186 &0.148&848 &0.280&594 \\\hline
$\mathbf{1024\times 1024}$&0.209&1407 &0.556&960 &1.080&649 \\\hline
$\mathbf{2048\times 2048}$&0.801&1092 &\bf 2.189&\bf 802 &4.278&573 \\\hline
$\mathbf{4096\times 4096}$&3.171&1075 &8.720&793 &17.076&569 \\\hline
\end{tabular}
}  
\caption{Timings ($time$) and throughput values ($TP$ in Mpix/s) of one register-only non separable convolution kernel, for small mask sizes of $3\times 3$, $5\times 5$ and $7\times 7$ pixels, on a GTX280 (GT200 architecture). Data transfer duration are those of Table \ref{tab:memcpy1}.}
\label{tab:convoNonSepReg3}
\end{table}

\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|}
\hline
\shortstack{\textbf{GPU card$\rightarrow$}\\\textbf{Image size$\downarrow$}}&\textbf{C2070}&\textbf{GTX280}\\\hline\hline
$\mathbf{512\times 512}$  &0.148 &0.161 \\\hline
$\mathbf{1024\times 1024}$&0.435 &0.536 \\\hline
$\mathbf{2048\times 2048}$&1.530 &3.039 \\\hline
$\mathbf{4096\times 4096}$&5.882 &12.431 \\\hline
\end{tabular}
}  
\caption{Time cost of data transfers between CPU and GPU memories, on C2070 and GTX280 cards (in milliseconds).}
\label{tab:memcpy1}
\end{table}

Table \ref{tab:convoNonSepReg3} shows timings and global throughput values achieved by those convolution masks on an Nvidia GT200 Tesla architecture (GTX480 card) with $16x8$ thread blocks. This measurement has been done in order to make a relevant comparison with a reference given by Nvidia in \cite{convolutionsoup} where they state that their fastest kernel achieves a $5\times5$ convolution of an 8-bit  $2048\times 2048$ pixelimage in $1.4~ms$, which lead to a throughput value of 945~Mpix/s.
Our current value of 802~Mpix/s, though not unsatisfactory, remains lower to the one reached by the manufacturer'own coding. 
Tested in the same conditions, the newer Fermi architecture of
Nvidia's GPUs proved slower (3.3 ms, see Table \ref{tab:convoNonSepReg1}) due to the lower maximum
register count allowed (63, against 128 for Tesla GT200).

It is interesting to note that, as long as each thread processes one single pixel, kernel execution time is ruled in proportion
with the number of pixels in the image multiplied by that of the mask. 
The slope in this first implementaion is $3.14.10^{-8}~ms/pix$ on C2070. 

\subsection{Using parameterizable masks}

To further improve the above implementation, it becomes necessary
to free ourselves from the hard-coding constraint. To achieve this,
as was the case with input image storing, several memory options are
available, but, since the amount of data involved in processing a
mask is quite small and constant, we considered it relevant to copy data
into \emph{symbol memory}. Listing \ref{lst:symbolmem} details the process, involving
the Cuda function \emph{CudaMemCopyToSymbol()}.

\lstinputlisting[label={lst:symbolmem},caption=code snippet showing how to setup a mask in GPU symbol memory]{Chapters/chapter4/code/maskInSymbol.cu}

In parallel, giving up the register-only constraint allows a more
conventional coding practice (loops). Listing \ref{lst:convoGene8r} presents
a generic convolution kernel, whose code immediately
appears both simple and concise. Its global time
performance, however, is comparatively lower than the register-only
process, due to the use of constant memory and of the \emph{r} parameter
(radius of the mask). The average slope amounts to $3.81~ms/pix$ on C2070,
which means a time-cost increase of around $20~\%$.

\lstinputlisting[label={lst:convoGene8r},caption=Generic CUDA kernel achieving a convolution operation with the mask in symbol memory and its radius passed as a parameter]{Chapters/chapter4/code/convoGene8r.cu}

\subsection{Increasing the number of pixels processed by each thread}
Much in the same way as we did with the Median Filter, we shall now
attempt to reduce the average latency due to writes into global memory
by having each thread process more than one output value. As the basic
structure of the above GPU kernel uses only 14 registers per thread, regardless
of the size of the convolution mask, one can envisage processing 2
or more pixels per thread while keeping safely within the 63-per-thread
rule.

However, when doing so, \textit{eg} processing what we shall call a \textit{packet} of pixels, window mask overlapping has to be taken into account
to avoid multiple texture fetches of each pixel's gray-level value, while benefiting from the 2-D cache.
In that case, both mask size and pixel packet shape determine the number of texture fetches to be performed for each pixel value.
Figure \ref{fig:convoOverlap1} illustrates two different situations: on top, a mask of radius 1 ($3\times 3$) applied to a packet of 8 pixels in row; at bottom, a mask of radius 2 ($5\times 5$).
The dark gray pixels are the center pixels (pixels of the packet), while light gray pixels belong to the halo around the packet. The number in each pixel box corresponds to the convolution count in which it is involved. 
There would be little interest in using different \textit{packet} shapes, as the final global memory writes would not be coalescent; generating multiple latencies.  
 \begin{figure}
\centering
   \subfigure[$3\times 3$ mask: there are 18 center pixels (out of 30) involved in 3 computations.]{ \includegraphics[width=5.8cm]{Chapters/chapter4/img/convoOverlap1.png}}\\
   \subfigure[$5\times 5$ mask: only 20 center pixels (out of 60), involved in 5 computations.]{ \includegraphics[width=7cm]{Chapters/chapter4/img/convoOverlap2.png}}
   \caption{Mask window overlapping when processing 8 pixels per thread. Top: $3\times 3$ mask. Bottom: $5\times 5$ mask.}
   \label{fig:convoOverlap1}
\end{figure}

Altough we actually wrote GPU kernels able to process 2, 4, 8 and 16 pixels per thread, only the one that processes 8 pixels per thread is presented below, as it proved to be the fastest one. Listing \ref{lst:convoGene8x8pL3} reproduce the source code of the kernel for $3\times 3$ masks.
The bottom line is that each thread is associated with one base pixel of coordinates $(x,y)$ which is the first of the packet to be processed, the last one being $(x+7,y)$. 
\lstinputlisting[label={lst:convoGene8x8pL3},caption=CUDA kernel achieving a $3\times 3$ convolution operation with the mask in symbol memory and direct data fetches in texture memory]{Chapters/chapter4/code/convoGene8x8pL3.cu}

In this particular case of a $3\times 3$ mask, each pixel value is used in 3 different convolution sums, except pixels located near both ends of the packet, whose values are used in fewer sums.
The general rule, when performing a $n\times n$ convolution (radius $k$) by 8-pixel packets is that each of the $(8-2k).(2k+1)$ \textit{center} pixels of the halo is used in $k$ sums, while the $4k.(2k+1)$ remaining pixels, located around the ends of the packet are used in fewer sums, from $k-1$ to $1$ ($2.(2k+1)$ pixels each).         
\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\textbf{Mask size$\rightarrow$}&\multicolumn{2}{|c|}{\textbf{3x3}}&\multicolumn{2}{|c|}{\textbf{5x5}}&\multicolumn{2}{|c|}{\textbf{7x7}}\\
\textbf{Image size$\downarrow$}&time (ms)&TP&time (ms)&TP&time (ms)&TP\\\hline\hline
$\mathbf{512\times 512}$  &0.036&1425 &0.069&1208 &0.110&1016 \\\hline
$\mathbf{1024\times 1024}$&0.128&1862 &0.253&1524 &0.413&1237 \\\hline
$\mathbf{2048\times 2048}$&0.495&2071 &\bf 0.987&1666 &1.615&1334 \\\hline
$\mathbf{4096\times 4096}$&1.964&2138 &3.926&1711 &6.416&1364 \\\hline
\end{tabular}
}  
\caption{Timings ($time$) and throughput values ($TP$ in Mpix/s) of our generic fixed mask size convolution kernel run on a C2070 card. Data transfer durations are those of Table \ref{tab:memcpy1}.}
\label{tab:convoGene8x8p}
\end{table}
 
Timing results and throughput values are shown in Table \ref{tab:convoGene8x8p}, and show that this solution now outperforms Nvidia references. 
It is important to remember that the above kernels have been optimized for the Fermi architecture, unlike those mentioned earlier, which were more efficient on the GT200 architecture.  
However, our technique requires to write one kernel per mask size, which can be seen as a major constraint. To make it easier to use this method, we shall propose a kernel code generator that will be available in the near future. 

\subsection{Using shared memory to store prefetched data\index{Prefetching}.}
 \index{memory~hierarchy!shared~memory}
A more convenient way of coding a convolution kernel is to use shared memory to perform a prefetching stage of the whole halo before computing the convolution sums.
This proves to be quite efficient and more versatile, but it obviously generates some overhead as:
\begin{itemize}
\item Each pixel value has to be read at least twice, first from texture memory into shared memory and then one or several more times from shared memory to be used in convolution computations.
\item Reducing the number of times a single pixel value is read from shared memory is bound to generate bank conflicts, hence once again performance loss.    
\end{itemize}
 \begin{figure}
\centering
   \includegraphics[width=12cm]{Chapters/chapter4/img/convoShMem.png}
   \caption{Organization of the prefetching stage of data, for a $5\times 5$ mask and a thread block size of $8\times 4$. Threads in both top corners of the top figure are identified either by a circle or by a star symbol. The image tile, loaded into shared memory includes the pixels to be updated by the threads of the block, as well as its 2-pixel wide halo. Here, circle and star symbols in the image tile show which pixels are actually loaded into one shared memory vector by its corresponding thread. }
   \label{fig:ShMem1}
\end{figure}
Still, we also implemented this method, in a similar manner as Nvidia did in its SDK sample code.
Some improvement has been obtained by increasing the number of pixels processed by each thread, to an optimum 8 pixels per thread.
The principle is to prefetch all pixel values involved in the computations performed by all threads of a block, including 8 pixels per thread plus the halo of radius $r$ (the radius of the convolution mask). As this obviously represents more values than the thread count in one block, some threads have to load more than one value.
The general organization is reproduced in Figure \ref{fig:ShMem1} for $5\times 5$ mask and a $8\times 4$ thread block, while Listing \ref{lst:convoGeneSh1} gives the details of the implementation with its two distinct code blocks: preload in shared memory (Lines 20 to 42) and convolution computations (Lines 45 to 57).    
Table \ref{tab:convoGeneSh1} details timing results of this implementation ($16\times 8$ threads/block), up to $13\times 13$ masks, that will serve as a reference in the next section, devoted to separable convolution. 
\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\shortstack{\textbf{Mask size$\rightarrow$}\\\textbf{Image size$\downarrow$}}&\textbf{3x3}&\textbf{5x5}&\textbf{7x7}&\textbf{9x9}&\textbf{11x11}&\textbf{13x13}\\\hline\hline
$\mathbf{512\times 512}$  &0.040 &0.075 &0.141    &0.243&0.314&0.402\\\hline
$\mathbf{1024\times 1024}$&0.141 &0.307 &0.524    &0.917&1.192&1.535\\\hline
$\mathbf{2048\times 2048}$&0.543 &\bf 1.115&2.048 &3.598&4.678&6.037\\\hline
$\mathbf{4096\times 4096}$&2.146 &4.364 &8.156    &14.341&18.652&24.020\\\hline
\end{tabular}
}  
\caption{Performances, in milliseconds, of our generic 8 pixels per thread kernel using shared memory, run  on a C2070 card.}
\label{tab:convoGeneSh1}
\end{table}
\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\shortstack{\textbf{Mask size$\rightarrow$}\\\textbf{Image size$\downarrow$}}&\textbf{3x3}&\textbf{5x5}&\textbf{7x7}&\textbf{9x9}&\textbf{11x11}&\textbf{13x13}\\\hline\hline
$\mathbf{512\times 512}$  &1394 &1176 &907      &670&567&477\\\hline
$\mathbf{1024\times 1024}$&1820 &1413 &1093     &776&644&532\\\hline
$\mathbf{2048\times 2048}$&2023 &\bf 1586 &1172 &818&676&554\\\hline
$\mathbf{4096\times 4096}$&2090 &1637 &1195     &830&684&561\\\hline
\end{tabular}
}  
\caption{Throughput values, in MegaPixel per second, of our generic 8 pixels per thread kernel using shared memory, run on a C2070 card. Data transfer durations are those of Table \ref{tab:memcpy1}.}
\label{tab:convoGeneSh2}
\end{table} 
\lstinputlisting[label={lst:convoGeneSh1},caption=CUDA kernel achieving a generic convolution operation after a preloading of data in shared memory.]{Chapters/chapter4/code/convoGeneSh1.cu}

\section{Separable convolution}
A convolution operation is said separable when its masks $h$ is the product of 2 vectors $h_v$ and $h_h$, as is the case in the following example:
$$h = h_v \times h_h = \begin{bmatrix}1\\2\\1\end{bmatrix} \times \begin{bmatrix}-1&2&-1\end{bmatrix} = \begin{bmatrix}
-1&2&-1\\
-2&4&-2\\
-1&2&-1
\end{bmatrix}$$
Such a mask allows to replace a generic 2-D convolution operation by two consecutive stages of a 1-D convolution operation: a vertical of mask $h_v$ and a horizontal of mask $h_h$.
This saves a lot of arithmetic operations, as a generic $n\times n$ convolution applied on a $H\times L$ image basically represents $H.L.n^2$ multiplications and as many additions, while two consecutive $n\times 1$ convolutions only represents $2.H.L.n$ of each, \textit{eg} 60\% operations are saved per pixel of the image for a $5\times 5$ mask.\\
However, beside reducing the operation count, performing a separable convolution also means writing an intermediate image into global memory.
CPU implementations of separable convolutions often use a single function to perform both 1-D convolution stages. To do so, this function reads the input image and actually ouputs the transposed filtered image. 
Applying that principle to GPUs is not efficient, as outputting the transposed image means non-coalescent writes into global memory, generating severe performance loss. Hence the idea of developing two different kernels, one for each of both vertical and horizontal convolutions.

Here, the use of Shared memory is the best choice, as there is no overlapping between neighbor windows and thus no possible optimization.
Moreover, to ensure efficiency, it is important to read the input image from texture memory, which implies an internal GPU data copy between both 1-D convolution stages.
Which, even if it is faster than CPU/GPU data transfer, makes separable convolutions slower than generic convolutions for small mask sizes. On C2070, the lower limit is $7\times 7$ pixels ($9\times 9$ for $512\times 512$ images).

Both vertical and horizontal kernels feature similar runtimes: Table \ref{tab:convoSepSh1} only contains their average execution time, including the internal data copy stage, while Table \ref{tab:convoSepSh2} shows the  achieved global throughput values. Timings of the data copy stage are given in Table \ref{tab:cpyToArray}. 
Listings \ref{lst:convoSepShV} and \ref{lst:convoSepShH} detail the implementation of both 1-D kernels, while Listing \ref{lst:convoSepSh} shows how to use them in addition with the data copy function in order to achieve a whole separable convolution. The shared memory size is dynamically passed as a parameter at kernel call time. Its expression is given in the comment line before its declaration.
\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|}
\hline
\textbf{Image size}&\textbf{C2070}\\\hline\hline
$\mathbf{512\times 512}$  &0.029 \\\hline
$\mathbf{1024\times 1024}$&0.101 \\\hline
$\mathbf{2048\times 2048}$&0.387 \\\hline
$\mathbf{4096\times 4096}$&1.533 \\\hline
\end{tabular}
}  
\caption{Time cost of data copy between the vertical and the horizontal 1-D convolution stages, on a C2070 cards (in milliseconds).}
\label{tab:cpyToArray}
\end{table}
\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\shortstack{\textbf{Mask size$\rightarrow$}\\\textbf{Image size$\downarrow$}}&\textbf{3x3}&\textbf{5x5}&\textbf{7x7}&\textbf{9x9}&\textbf{11x11}&\textbf{13x13}\\\hline\hline
$\mathbf{512\times 512}$  &0.080 &0.087 &0.095 &\bf 0.108&\bf 0.115&\bf 0.126\\\hline
$\mathbf{1024\times 1024}$&0.306 &0.333 &\bf 0.333 &\bf 0.378&\bf 0.404&\bf 0.468\\\hline
$\mathbf{2048\times 2048}$&1.094 &1.191 &\bf 1.260 &\bf 1.444&\bf 1.545&\bf 1.722\\\hline
$\mathbf{4096\times 4096}$&4.262 &4.631 &\bf 5.000 &\bf 5.676&\bf 6.105&\bf 6.736\\\hline
\end{tabular}}  
\caption{Performances, in milliseconds, of our generic 8 pixels per thread 1-D convolution kernels using shared memory, run  on a C2070 card. Timings include data copy. Bold values correspond to situations where separable-convolution kernels run faster than non separable ones.}
\label{tab:convoSepSh1}
\end{table}
\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\shortstack{\textbf{Mask size$\rightarrow$}\\\textbf{Image size$\downarrow$}}&\textbf{3x3}&\textbf{5x5}&\textbf{7x7}&\textbf{9x9}&\textbf{11x11}&\textbf{13x13}\\\hline\hline
$\mathbf{512\times 512}$  &1150 &1116 &1079 &\bf 1024&\bf 997 &\bf 957\\\hline
$\mathbf{1024\times 1024}$&1415 &1365 &\bf 1365 &\bf 1290&\bf 1250&\bf 1169\\\hline
$\mathbf{2048\times 2048}$&1598 &1541 &\bf 1503 &\bf 1410&\bf 1364&\bf 1290\\\hline
$\mathbf{4096\times 4096}$&1654 &1596 &\bf 1542 &\bf 1452&\bf 1400&\bf 1330\\\hline
\end{tabular}
}  
\caption{Throughput values, in MegaPixel per second, of our generic 8 pixels per thread 1-D convolution kernel using shared memory, run on a C2070 card. Data transfer durations are those of Table \ref{tab:memcpy1}.}
\label{tab:convoSepSh2}
\end{table} 

\lstinputlisting[label={lst:convoSepSh},caption=data copy between the calls to 1-D convolution kernels achieving a 2-D separable convolution operation.]{Chapters/chapter4/code/convoSepSh.cu}
\lstinputlisting[label={lst:convoSepShV},caption=CUDA kernel achieving a horizontal 1-D convolution operation after a preloading \index{Prefetching} of data in shared memory.]{Chapters/chapter4/code/convoSepShV.cu}
\lstinputlisting[label={lst:convoSepShH},caption=CUDA kernel achieving a vertical 1-D convolution operation after a preloading of data in shared memory.]{Chapters/chapter4/code/convoSepShH.cu}
 
\section{Conclusion}
Extensively detailing the various techniques that may be applied when designing a median or a convolution operation on GPU has enabled us determine that:
\begin{itemize}
\item the use of registers with direct data fetching from texture often allows kernels to run faster than those which use the more conventionnal way of prefetching data from texture memory and storing them into shared memory.
\item increasing the pixel count processed by each thread brings important speedups. In this case, if neighboring windows overlap, optimized direct data fetching from texture will likely outperform the shared memory prefetching technique. That is the case for generic convolution kernels.
\item coding such optimized data fetching is not straightforward. Consequently, we are planning to provide a kernel code generator that will make our kernels more accessible by GPU users.
\end{itemize}
The presented kernels, optimized for a C2070 card, achieve up to 2138~Mpix/s including data transfers, which comes close to the absolute maximum throughput value allowed by the Fermi architecture. The next GPU generation (Kepler) may allow us not only to benefit from  new dynamic parallelism capability to increase kernel paralelism level, but also to take advantage of an increase of the register count allowed per thread block which would allow us, for example, to extend our register-only median filter technique to larger mask sizes.

\putbib[Chapters/chapter4/biblio4]