1 \chapterauthor{Gilles Perrot}{Femto-ST Institute, University of Franche-Comte, France}
3 \chapter{Implementing an efficient convolution operation on GPU}
10 In this chapter, after dealing with GPU median filter implementations,
11 we propose to explore how convolutions\index{convolution} can be implemented on modern
12 GPUs. Widely used in digital image processing filters, the \emph{convolution
13 operation} basically consists of taking the sum of products of elements
14 from two 2D functions, letting one of the two functions move over
15 every element of the other, producing a third function that is typically
16 viewed as a modified version of one of the original functions. To
17 begin with, we shall examine nonseparable or generic convolutions,
18 before addressing the matter of separable convolutions. We shall refer
19 to $I$ as an $H\times L$ pixel gray-level image and to $I(x,y)$ as the gray-level
20 value of each pixel of coordinates $(x,y)$.
25 Within a digital image $I$, the convolution operation is performed between
26 image $I$ and convolution mask \emph{h} (To avoid confusion with other
27 GPU functions referred to as kernels, we shall use\emph{ convolution
28 mask} instead of \emph{convolution kernel}) is defined by
30 I'(x, y) = \left(I * h\right) = \sum_{(i < H)} \sum_{(j < L)}I(x-j, y-j)h(j,i)
33 While processing an image, function \emph{h} is often bounded by a square
34 window of size \emph{k = 2r + 1}, i.e., an uneven number, to ensure
35 there is a center. We shall also point out that, as stated earlier,
36 the square shape is not a limiting factor to the process, as any shape
37 can be inscribed into a square. In the case of a more complex shape,
38 the remaining space is filled by null values (padding).
41 \section{Implementation}
42 The basic principle of computing a convolution between one $I$ picture
43 and one \emph{h} convolution mask defined on domain $\Omega$ is given
44 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.
45 For more readability, only part of the connecting lines are shown.
48 \includegraphics[width=11cm]{Chapters/chapter4/img/convo1.png}
49 \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.}
50 \label{fig:convoPrinciple}
53 \caption{generic convolution}
55 \ForEach{pixel at position $(x, y)$}{
56 Read all gray-level values $I(x, y)$ in the neighborhood\;
57 Compute the weighted sum \( I_\Omega = \sum_{(j,i) \in \Omega}I(x-j, y-j)h(j,i) \)\;
58 Normalize $I'(x, y)$ value\;
59 Output the new gray-level value
63 The gray-level value of each pixel of output image $I'$ is the weighted
64 sum of pixels included in the neighborhood defined by $\Omega$ around
65 the corresponding pixel in the input image. It has to be noted that,
66 in case the sum $S$ of all coefficients in the mask is not 1, the original
67 brightness of the image will be altered and a normalization stage
68 has to take place, as, for example, in the case of an 8-bit coded
71 \item if $S > 0$ then $I' = I_\Omega / S$
72 \item if $S = 0$ then $I' = I_\Omega + 128$
73 \item if $S < 0$ then $I' = I_\Omega + 255$
75 In case one, normalizing means performing a division operation for
76 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.
79 \subsection{First test implementation}
80 This first implementation consists of a rather naive application to
81 convolutions of the techniques applied to median filters in the
82 previous chapter, as a reminder: texture memory used with incoming
83 data, pinned memory with output data, optimized use of registers
84 while processing data and multiple output per thread\index{multiple output per thread}.
85 One significant difference lies in the fact
86 that the median filter uses only one parameter, the size of the window mask,
87 which can be hard-coded, while a convolution mask requires referring to several parameters; hard-coding
88 the elements of the mask would lead to severe lack of flexibility (one function
89 per filter, no external settings) so we will just use it as a starting
90 point in our approach.
92 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):
93 $$h=\frac{1}{9}\begin{bmatrix}1&1&1\\1&1&1\\1&1&1\end{bmatrix}$$
94 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.
95 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.
97 \lstinputlisting[label={lst:convoGene3Reg8},caption=generic CUDA kernel achieving a convolution operation with hard-coded mask values]{Chapters/chapter4/code/convoGene3Reg8.cu}
99 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
100 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.
106 \begin{tabular}{|c||r|r||r|r||r|r|}
108 \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}$}\\
109 \textbf{Image size}$\downarrow$&time (ms)&TP&time (ms)&TP&time (ms)&TP\\\hline\hline
110 $\mathbf{512\times 512}$ &0.077&1165 &0.209&559 &0.407 &472 \\\hline
111 $\mathbf{1024\times 1024}$&0.297&1432 &0.820&836 &1.603 &515 \\\hline
112 $\mathbf{2048\times 2048}$&1.178&1549 &\bf 3.265&\bf 875 &6.398&529 \\\hline
113 $\mathbf{4096\times 4096}$&4.700&1585 &13.05&533 &25.56&533 \\\hline
116 \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.}
117 \label{tab:convoNonSepReg1}
126 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.
127 Our current value of 802~MP/s, though not unsatisfactory, remains lower to the one reached by the manufacturer's own coding.
128 Tested in the same conditions, the newer Fermi architecture of
129 NVIDIA's GPUs proved slower (3.3 ms, see Table \ref{tab:convoNonSepReg1}) due to the lower maximum
130 register count allowed (63 as opposed to 128 for Tesla GT200).
135 \begin{tabular}{|c||r|r||r|r||r|r|}
137 \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}$}\\
138 \textbf{Image size}$\downarrow$&time (ms)&TP&time (ms)&TP&time(ms)&TP\\\hline\hline
139 $\mathbf{512\times 512}$ &0.060&1186 &0.148&848 &0.280&594 \\\hline
140 $\mathbf{1024\times 1024}$&0.209&1407 &0.556&960 &1.080&649 \\\hline
141 $\mathbf{2048\times 2048}$&0.801&1092 &\bf 2.189&\bf 802 &4.278&573 \\\hline
142 $\mathbf{4096\times 4096}$&3.171&1075 &8.720&793 &17.076&569 \\\hline
145 \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.}
146 \label{tab:convoNonSepReg3}
149 It is interesting to note that, as long as each thread processes one single pixel, kernel execution time is ruled in proportion
150 with the number of pixels in the image multiplied by that of the mask.
151 The proportionality factor, that we call \textit{slope}, is $3.14.10^{-8}$~ms/pix on C2070 in this first implementation.
152 As a reminder, Table \ref{tab:memcpy1} details the data transfer costs that helped in computing throughput values.
156 \begin{tabular}{|c||r|r|}
158 \shortstack{\textbf{GPU card}$\rightarrow$\\\textbf{Image size$\downarrow$}}&\textbf{C2070}&\textbf{GTX280}\\\hline\hline
159 $\mathbf{512\times 512}$ &0.148 &0.161 \\\hline
160 $\mathbf{1024\times 1024}$&0.435 &0.536 \\\hline
161 $\mathbf{2048\times 2048}$&1.530 &3.039 \\\hline
162 $\mathbf{4096\times 4096}$&5.882 &12.431 \\\hline
165 \caption{Time cost of data transfers between CPU and GPU memories, on C2070 and GTX280 cards (in milliseconds).}
169 \subsection{Using parameterizable masks}
170 To further improve the above implementation, it becomes necessary
171 to free ourselves from the hard-coding constraint. To achieve this,
172 as was the case with input image storing, several memory options are
173 available, but, since the amount of data involved in processing a
174 mask is quite small and constant, we considered it relevant to copy data
175 into \emph{symbol memory}. Listing \ref{lst:symbolmem} details this process, involving
176 the CUDA function \emph{cudaMemcpyToSymbol()}.
178 \lstinputlisting[label={lst:symbolmem},caption=code snippet showing how to setup a mask in GPU symbol memory]{Chapters/chapter4/code/maskInSymbol.cu}
180 In parallel, giving up the register-only constraint allows a more
181 conventional coding practice (loops). Listing \ref{lst:convoGene8r} presents
182 a generic convolution kernel, whose code immediately
183 appears both simple and concise. Its global time
184 performance, however, is comparatively lower than the register-only
185 process, due to the use of constant memory and of the \emph{r} parameter
186 (radius of the mask). The average slope amounts to $3.81.10^{-8}$~ms/pix on C2070,
187 which means a time-cost increase of around 20~\%.
189 \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}
191 \subsection{Increasing the number of pixels processed by each thread}
192 Much in the same way as we did with the Median Filter, we shall now
193 attempt to reduce the average latency due to writes into global memory
194 by having each thread process more than one output value. As the basic
195 structure of the above GPU kernel uses only 14 registers per thread, regardless
196 of the size of the convolution mask, one can envisage processing 2
197 or more pixels per thread while keeping safely within the 63-per-thread
200 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
201 to avoid multiple texture fetches of each pixel's gray-level value, while benefiting from the 2D cache.
202 In that case, both mask size and pixel packet shape determine the number of texture fetches to be performed for each pixel value.
203 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$).
204 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.
205 There would be little interest in using different \textit{packet} shapes, as the final global memory writes would not be coalescent, generating multiple latencies.
208 \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}}\\
209 \subfigure[$5\times 5$ mask: only 20 pixels (out of 60) are involved in 5 computations.]{ \includegraphics[width=7cm]{Chapters/chapter4/img/convoOverlap2.png}}
210 \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.}
211 \label{fig:convoOverlap1}
214 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.
215 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)$.
217 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.
218 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).
222 \begin{tabular}{|c||r|r||r|r||r|r|}
224 \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}$}\\
225 \textbf{Image size}$\downarrow$&time (ms)&TP&time (ms)&TP&time (ms)&TP\\\hline\hline
226 $\mathbf{512\times 512}$ &0.036&1425 &0.069&1208 &0.110&1016 \\\hline
227 $\mathbf{1024\times 1024}$&0.128&1862 &0.253&1524 &0.413&1237 \\\hline
228 $\mathbf{2048\times 2048}$&0.495&2071 &\bf 0.987&1666 &1.615&1334 \\\hline
229 $\mathbf{4096\times 4096}$&1.964&2138 &3.926&1711 &6.416&1364 \\\hline
232 \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.}
233 \label{tab:convoGene8x8p}
236 Timing results and throughput values are shown in Table \ref{tab:convoGene8x8p}, and show that this solution now outperforms NVIDIA references.
237 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.
238 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.
240 \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}
242 \subsection{Using shared memory to store prefetched data\index{prefetching}}
243 \index{memory hierarchy!shared memory}
244 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.
245 This proves to be quite efficient and more versatile, but it obviously generates some overhead because
247 \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.
248 \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.
252 \includegraphics[width=12cm]{Chapters/chapter4/img/convoShMem.png}
253 \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. }
256 Still, we also implemented this method, in a similar manner as NVIDIA did in its SDK sample code.
257 Some improvement has been obtained by increasing the number of pixels processed by each thread, to an optimum 8 pixels per thread.
258 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.
259 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).
260 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.
264 \begin{tabular}{|c||r|r|r|r|r|r|}
266 \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
267 $\mathbf{512\times 512}$ &0.040 &0.075 &0.141 &0.243&0.314&0.402\\\hline
268 $\mathbf{1024\times 1024}$&0.141 &0.307 &0.524 &0.917&1.192&1.535\\\hline
269 $\mathbf{2048\times 2048}$&0.543 &\bf 1.115&2.048 &3.598&4.678&6.037\\\hline
270 $\mathbf{4096\times 4096}$&2.146 &4.364 &8.156 &14.341&18.652&24.020\\\hline
273 \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.}
274 \label{tab:convoGeneSh1}
279 \begin{tabular}{|c||r|r|r|r|r|r|}
281 \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
282 $\mathbf{512\times 512}$ &1394 &1176 &907 &670&567&477\\\hline
283 $\mathbf{1024\times 1024}$&1820 &1413 &1093 &776&644&532\\\hline
284 $\mathbf{2048\times 2048}$&2023 &\bf 1586 &1172 &818&676&554\\\hline
285 $\mathbf{4096\times 4096}$&2090 &1637 &1195 &830&684&561\\\hline
288 \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}.}
289 \label{tab:convoGeneSh2}
291 \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}
293 \section{Separable convolution}
294 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:
295 $$h = h_v \times h_h = \begin{bmatrix}1\\2\\1\end{bmatrix} \times \begin{bmatrix}-1&2&-1\end{bmatrix} = \begin{bmatrix}
300 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$.
301 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.
303 However, besides reducing the operation count, performing a separable convolution also means writing an intermediate image into global memory.
304 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.
305 Applying this principle to GPUs is not efficient, as outputting the transposed image means noncoalescent 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.
307 Here, the use of shared memory is the best choice, as there is no overlapping between neighbor windows and thus no possible optimization.
308 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.
309 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).
311 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}.
312 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.
317 \begin{tabular}{|c||r|r|r|r|r|r|}
319 \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
320 $\mathbf{512\times 512}$ &0.080 &0.087 &0.095 &\bf 0.108&\bf 0.115&\bf 0.126\\\hline
321 $\mathbf{1024\times 1024}$&0.306 &0.333 &\bf 0.333 &\bf 0.378&\bf 0.404&\bf 0.468\\\hline
322 $\mathbf{2048\times 2048}$&1.094 &1.191 &\bf 1.260 &\bf 1.444&\bf 1.545&\bf 1.722\\\hline
323 $\mathbf{4096\times 4096}$&4.262 &4.631 &\bf 5.000 &\bf 5.676&\bf 6.105&\bf 6.736\\\hline
325 \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 nonseparable ones.}
326 \label{tab:convoSepSh1}
331 \begin{tabular}{|c||r|r|r|r|r|r|}
333 \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
334 $\mathbf{512\times 512}$ &1150 &1116 &1079 &\bf 1024&\bf 997 &\bf 957\\\hline
335 $\mathbf{1024\times 1024}$&1415 &1365 &\bf 1365 &\bf 1290&\bf 1250&\bf 1169\\\hline
336 $\mathbf{2048\times 2048}$&1598 &1541 &\bf 1503 &\bf 1410&\bf 1364&\bf 1290\\\hline
337 $\mathbf{4096\times 4096}$&1654 &1596 &\bf 1542 &\bf 1452&\bf 1400&\bf 1330\\\hline
340 \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 nonseparable ones (data transfer durations are those of Table \ref{tab:memcpy1}).}
341 \label{tab:convoSepSh2}
346 \begin{tabular}{|c||r|}
348 \textbf{Image size}&\textbf{C2070}\\\hline\hline
349 $\mathbf{512\times 512}$ &0.029 \\\hline
350 $\mathbf{1024\times 1024}$&0.101 \\\hline
351 $\mathbf{2048\times 2048}$&0.387 \\\hline
352 $\mathbf{4096\times 4096}$&1.533 \\\hline
355 \caption{Time cost of data copy between the vertical and the horizontal 1D convolution stages, on a C2070 cards (in milliseconds).}
356 \label{tab:cpyToArray}
358 \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}
359 \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}
360 \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}
363 Extensively detailing the various techniques that may be applied when designing a median or a convolution operation on GPU has enabled us determine that
365 \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.
366 \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.
367 \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.
369 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.
371 \putbib[Chapters/chapter4/biblio4]
