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

\chapterauthor{Gilles Perrot}{Femto-ST Institute, University of Franche-Comte, France}

\chapter{Implementing an efficient convolution operation on GPU}





\section{Overview}
In this chapter, after dealing with GPU median filter implementations,
we propose to explore how convolutions\index{convolution}  can be implemented on modern
GPUs. Widely used in digital image processing filters, the \emph{convolution
operation} basically consists of taking the sum of products of elements
from two 2D 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 addressing the matter of separable convolutions. We shall refer
to $I$ as an $H\times L$ pixel gray-level image and to $I(x,y)$ as the gray-level
value of each pixel of coordinates $(x,y)$.


\clearpage
\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 often 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 not a 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.]{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\;
    Output 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 > 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 techniques 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 significant 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 parameters; hard-coding
the elements of the mask 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=\frac{1}{9}\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 outputs 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 gray-level value of 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 multiplying  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$, and $4096\times 4096$ run on a C2070 card with $32\times 8$ thread blocks.


\begin{table}[htbp]
\centering
{\normalsize
\begin{tabular}{|c||r|r||r|r||r|r|}
\hline
\textbf{Mask size}$\rightarrow$&\multicolumn{2}{c||}{$\mathbf{3\times 3}$}&\multicolumn{2}{c||}{$\mathbf{5\times 5}$}&\multicolumn{2}{c|}{$\mathbf{7\times 7}$}\\
\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 MP/s) of one register-only nonseparable convolution kernel, for small mask sizes of $3\times 3$, $5\times 5$, and $7\times 7$ pixels, on a C2070 card.]{Timings (time) and throughput values (TP in MPx/s) of one register-only nonseparable 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}. The bold value points out the result obtained in the reference situation.}
\label{tab:convoNonSepReg1}
\end{table} 







Table \ref{tab:convoNonSepReg3} shows timings and global throughput values achieved by those convolution masks on an NVIDIA GT200 Tesla architecture (GTX280 card) with $16\times 8$ thread blocks. This measurement has been done in order to make a relevant comparison with a reference given by NVIDIA in \cite{convolutionsoup} in which they state that their fastest kernel achieves a $5\times 5$ convolution of an 8-bit  $2048\times 2048$ pixel image in $1.4~ms$, leading to a throughput value of 945~MP/s. In all the result tables, the values associated to this reference will be presented in boldface.
Our current value of 802~MP/s, though not unsatisfactory, remains lower to the one reached by the manufacturer's 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 as opposed to 128 for Tesla GT200).

\begin{table}[htbp]
\centering
{\normalsize
\begin{tabular}{|c||r|r||r|r||r|r|}
\hline
\textbf{Mask size}$\rightarrow$&\multicolumn{2}{c||}{$\mathbf{3\times 3}$}&\multicolumn{2}{c||}{$\mathbf{5\times 5}$}&\multicolumn{2}{c|}{$\mathbf{7\times 7}$}\\
\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 MP/s) of one register-only nonseparable convolution kernel, for small mask sizes of $3\times 3$, $5\times 5$, and $7\times 7$ pixels, on a GTX280.]{Timings (time) and throughput values (TP in MP/s) of one register-only nonseparable 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}. The bold value points out the result obtained in the reference situation.}
\label{tab:convoNonSepReg3}
\end{table}

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 proportionality factor, that we call \textit{slope},  is $3.14.10^{-8}$~ms/pix on C2070 in this first implementation. 
As a reminder, Table \ref{tab:memcpy1} details the data transfer costs that helped in computing throughput values.
\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}

\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 this process, involving
the CUDA function \emph{cudaMemcpyToSymbol()}.

\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.10^{-8}$~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, e.g., 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 2D 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: (a) a mask of radius 1 ($3\times 3$) applied to a packet of 8 pixels in a row; (b) 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}[htbp]
\centering
   \subfigure[$3\times 3$ mask: there are 18 pixels (out of 30) involved in 3 computations.]{ \includegraphics[width=5.8cm]{Chapters/chapter4/img/convoOverlap1.png}}\\
   \subfigure[$5\times 5$ mask: only 20 pixels (out of 60) are involved in 5 computations.]{ \includegraphics[width=7cm]{Chapters/chapter4/img/convoOverlap2.png}}
   \caption[Mask window overlapping when processing a packet of 8 pixels per thread.]{Mask window overlapping when processing a packet of 8 pixels per thread. The dark gray pixels are the center pixels, while light gray pixels belong to the halo. The number in each pixel box is the convolution count in which it is involved. (a) $3\times 3$ mask; (b) $5\times 5$ mask.}
   \label{fig:convoOverlap1}
\end{figure}

Although we actually have written 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} reproduces 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, in the packet, to be processed, the last one being $(x+7,y)$. 

In this particular case of a $3\times 3$ mask, each pixel value is used in 3 different convolution sums, except for pixels located near both ends of the packet, whose values are used in fewer sums.
The general rule, when performing an $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}[htbp]
\centering
{\normalsize
\begin{tabular}{|c||r|r||r|r||r|r|}
\hline
\textbf{Mask size}$\rightarrow$&\multicolumn{2}{c||}{$\mathbf{3\times 3}$}&\multicolumn{2}{c||}{$\mathbf{5\times 5}$}&\multicolumn{2}{c|}{$\mathbf{7\times 7}$}\\
\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 MP/s) of our generic fixed mask size convolution kernel run on a C2070 card.]{Timings (time) and throughput values (TP in MP/s) of our generic fixed mask size convolution kernel run on a C2070 card. Data transfer durations are those of Table \ref{tab:memcpy1}. The bold value points out the result obtained in the reference situation.}
\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 writing one kernel per mask size, which can be seen as a major constraint. To make it easier to use this method, we are working on a kernel code generator that is currently under development and will be made available in the near future. 

\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}

\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 because 
\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}[htbp]
\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$.]{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).    
Tables \ref{tab:convoGeneSh1} and \ref{tab:convoGeneSh2} detail timing results and throughput values 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}[htbp]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\shortstack{\textbf{Mask size}$\rightarrow$\\\textbf{Image size$\downarrow$}}&$\mathbf{3\times 3}$&$\mathbf{5\times 5}$&$\mathbf{7\times 7}$&$\mathbf{9\times 9}$&$\mathbf{11\times 11}$&$\mathbf{13\times 13}$\\\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. Data transfers duration are not included.}
\label{tab:convoGeneSh1}
\end{table}
\begin{table}[htbp]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\shortstack{\textbf{Mask size}$\rightarrow$\\\textbf{Image size$\downarrow$}}&$\mathbf{3\times 3}$&$\mathbf{5\times 5}$&$\mathbf{7\times 7}$&$\mathbf{9\times 9}$&$\mathbf{11\times 11}$&$\mathbf{13\times 13}$\\\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.]{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 us to replace a generic 2D convolution operation by two consecutive stages of a 1D 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 an $H\times L$ image basically represents $HLn^2$ multiplications and as many additions, while two consecutive $n\times 1$ convolutions represents only $2HLn$ of each, e.g.,  60\% operations are saved per pixel of the image for a $5\times 5$ mask.

However, besides 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 1D convolution stages. To do so, this function reads the input image and actually ouputs the transposed filtered image. 
Applying this 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 the 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 1D convolution stages.
This, 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} contains only 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 1D 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 both Listings (\ref{lst:convoSepShV} and \ref{lst:convoSepShH}), in the comment lines before its declaration.

\begin{table}[h]
\centering
{\normalsize
\begin{tabular}{|c||r|r|r|r|r|r|}
\hline
\shortstack{\textbf{Mask size}$\rightarrow$\\\textbf{Image size$\downarrow$}}&$\mathbf{3\times 3}$&$\mathbf{5\times 5}$&$\mathbf{7\times 7}$&$\mathbf{9\times 9}$&$\mathbf{11\times 11}$&$\mathbf{13\times 13}$\\\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 1D convolution kernels using shared memory, run  on a C2070 card.]{Performances, in milliseconds, of our generic 8 pixels per thread 1D 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$}}&$\mathbf{3\times 3}$&$\mathbf{5\times 5}$&$\mathbf{7\times 7}$&$\mathbf{9\times 9}$&$\mathbf{11\times 11}$&$\mathbf{13\times 13}$\\\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 1D convolution kernel using shared memory, run on a C2070 card.]{Throughput values, in MegaPixel per second, of our generic 8 pixels per thread 1D convolution kernel using shared memory, run on a C2070 card. Bold values correspond to situations where separable-convolution kernels run faster than non separable ones (data transfer durations are those of Table \ref{tab:memcpy1}).}
\label{tab:convoSepSh2}
\end{table} 
\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 1D convolution stages, on a C2070 cards (in milliseconds).}
\label{tab:cpyToArray}
\end{table}
\lstinputlisting[label={lst:convoSepSh},caption=data copy between the calls to 1D convolution kernels achieving a 2D separable convolution operation]{Chapters/chapter4/code/convoSepSh.cu}
\lstinputlisting[label={lst:convoSepShV},caption=CUDA kernel achieving a horizontal 1D convolution operation after a preloading \index{prefetching} of data into shared memory]{Chapters/chapter4/code/convoSepShV.cu}
\lstinputlisting[label={lst:convoSepShH},caption=CUDA kernel achieving a vertical 1D convolution operation after a preloading of data into 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 in 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. This is the case for generic convolution kernels.
\item coding such optimized data fetching is not straightforward. Consequently, we are currently developing 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~MP/s including data transfers, which comes close to the absolute maximum throughput value allowed by the Fermi architecture. The next GPU generation (called 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 in 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]