\section{Introduction}
Après avoir conçu des  filtres médians aux peformances élevées, nous avons cherché à en appliquer les principes à d'autres types d'algorithmes de filtrage.
Les filtres de convolution, par la diversité des traitements qu'ils permettent de réaliser et leur universalité, nous ont semblé être un objectif particulièrement intéressant.

Les principes et formulation de la convolution sont présentés au chapitre \ref{sec-op-base} et nous nous attacherons uniquement dans les paragraphes qui suivent à détailler les solutions et expérimentations qui permettent de concevoir des filtres de convolution performants sur GPU. Nous faisons l'hypothèse que les fonctions de convolution sont à support carré de taille impaire, permettant ainsi de considérer un \textit{pixel central}. Cette hypothèse ne constitue pas une restriction en termes de traitement car tout support non carré peut être étendu à un support carré même si dans ce cas de figure, l'exécution impliquera plus d'opérations que nécessaire et ne sera ainsi plus optimale. 

L'étude la plus complète et qui montre les performances les plus élevées émane du constructeur Nvidia lui-même dans \cite{convolutionsoup}. Nous l'avons présentée dans le chapitre \ref{} et nous rappellons simplement ici qu'elle a utilisé des modèles à architecture GT200 (GTX280) dont nous disposons également et qu'elle a choisi comme traitement de référence une convolution non-séparable de masque 5$\times$5 sur une image en profondeur 8 bits de 2048$\times$2048 pixels.
Leur implémentation la plus rapide effectue cette opération en 1.4~ms et permet un débit global (incluant les temps de transfert des données) de 945 millions de pixels à la seconde (MP/s). Elle est prise comme référence pour nos implémentations. 

\section{Implémentation générique de la convolution non séparable}

L'implémentation GPU de la  convolution non-séparable d'une fonction image $I$ par une fonction masque $h$ définie sur un support $\Omega$ peut-être décrite comme dans l'algorithme \ref{algo-convo-gene}. Pour le cas où la somme $S_h$ des valeurs du masque est différente de 1, l'image résultante $I'$ est obtenue après une normalisation nécessaire pour ne pas modifier l'intensité moyenne de l'image. Par exemple, pour une profondeur de 8 bits :
\begin{enumerate}
\item Si $S_h > 0$ alors $I' = I_{\Omega}/S_h$
\item Si $S_h = 0$ alors $I' = I_{\Omega} + 128$
\item Si $S_h < 0$ alors $I' = I_{\Omega} + 255$
\end{enumerate}


\begin{algorithm}
\caption{Convolution générique sur GPU}   
\label{algo-convo-gene}
  \ForEach(\tcc*[f]{\textbf{en parallèle}}){pixel $(x, y)$}{
    Lire les niveaux de gris $I(x, y)$ des voisins sur $\Omega$ \;
    Calculer la somme \( I_\Omega(x, y) = \sum_{(j,i) \in \Omega}I(x-j, y-j).h(j,i) \) \;
    Normaliser $I_{\Omega}(x, y)$ pour obtenir $I'(x, y)$ \;
    Mémoriser $I'(x, y)$ \;
  }
\end{algorithm}

Il est tout à fait possible d'envisager ici l'application brute des principes mis en \oe uvre pour les filtres médians. Cela conduit au code du listing \ref{lst-convo-gene3reg8} où les coefficients du masque (moyenneur 3$\times$3) sont fixés et mémorisés chacun dans un registre, le calcul de la somme s'effectuant également dans un registre. Pour éviter des opérations coûteuses comme la division, on remarque que la normalisation est évitée et pré-effectuée au niveau des coefficients du masque dont la somme est ainsi égale à 1.
Par ailleurs, pour des raisons de lisibilité de ce premier code, chaque thread ne traite ici qu'un seul pixel.

\lstinputlisting[label={lst-convo-gene3reg8},caption={Kernel réalisant la convolution par un masque moyenneur 3$\times$3 dont les coefficients normalisés sont codés \textit{en dur}, dans les registres du GPU.}]{Chapters/chapter6/code/convoGene3Reg8.cu}

Les performances de cette implémentation directe ont été regroupées dans les tableaux \ref{tab-convo-gene3reg8-480} et \ref{tab-convo-gene3reg8-2070} où l'on peut immédiatement constater que la solution optimale Nvidia demeure plus rapide. 
L'analyse plus détaillée nous apprend aussi que le modèle GTX280 exécute le kernel plus vite que le plus récent C2070, en raison d'un plus grand nombre de registres disponibles. Malgré tout, lorsqu'on prend en compte les temps de transfert des données, l'avantage va au C2070 qui réalise ce traitement à 875~MP/s.

\begin{table}
\centering
{\normalsize
\begin{tabular}{clrrrr}
\toprule
&&\multicolumn{4}{c}{Taille d'image}\\
Masque&&$\mathbf{512\times 512}$&$\mathbf{1024\times 1024}$&$\mathbf{2048\times 2048}$&$\mathbf{4096\times 4096}$\\
\midrule
\multirow{2}*{3$\times$3}& temps exéc. (ms)   & 0.077 & 0.297 & 1.178 & 4.700 \\
                         & débit global (MP/s)& 1165  & 1432  & 1549  & 1585  \\
\midrule
\multirow{2}*{5$\times$5}& temps exéc. (ms)   & 0.209 & 0.820 & {\bf 3.265} & 13.050\\
                         & débit global (MP/s)& 559   & 836   & {\bf 875}   & 533   \\
\midrule
\multirow{2}*{7$\times$7}& temps exéc. (ms)   & 0.407 & 1.603 & 6.398 & 25.560\\ 
                         & débit global (MP/s)& 472   & 515   & 529   & 533 \\
\bottomrule
\end{tabular}
}  
\caption{Performances des kernels effectuant la convolution non-séparable sur le modèle du listing \ref{lst-convo-gene3reg8}, sur GPU C2070. Le temps d'exécution correspond à la seule exécution du kernel. Le débit global intègre les temps de transfert. Les valeurs en gras correspondent au traitement de référence.}
\label{tab-convo-gene3reg8-2070}
\end{table} 

\begin{table}
\centering
{\normalsize
\begin{tabular}{clrrrr}
\toprule
&&\multicolumn{4}{c}{Taille d'image}\\
Masque&&$\mathbf{512\times 512}$&$\mathbf{1024\times 1024}$&$\mathbf{2048\times 2048}$&$\mathbf{4096\times 4096}$\\
\midrule
\multirow{2}*{3$\times$3}& temps exéc. (ms)   & 0.060 & 0.209 & 0.801 & 3.171 \\
                         & débit global (MP/s)& 1186  & 1407  & 1092  & 1075  \\
\midrule
\multirow{2}*{5$\times$5}& temps exéc. (ms)   & 0.148 & 0.556 & {\bf 2.189} & 8.7200\\
                         & débit global (MP/s)& 848   & 960   & {\bf 802}   & 793   \\
\midrule
\multirow{2}*{7$\times$7}& temps exéc. (ms)   & 0.280 & 1.080 & 4.278 & 17.076\\ 
                         & débit global (MP/s)& 594   & 649   & 573   & 569 \\
\bottomrule
\end{tabular}
}  
\caption{Performances des kernels effectuant la convolution non-séparable sur le modèle du listing \ref{lst-convo-gene3reg8}, sur GPU GTX280. Le temps d'exécution correspond à la seule exécution du kernel. Le débit global intègre les temps de transfert. Les valeurs en gras correspondent au traitement de référence.}
\label{tab-convo-gene3reg8-480}
\end{table}

\section{Implémentation optimisée de la convolution non séparable}

Les coefficients du masque de convolution sont indépendants et il est donc impossible, sauf cas particulier, de réduire le nombre de registres nécessaire. Pour cette même raison, multiplier le nombre de pixels traités par chaque thread ne permet pas d'économiser des registres au niveau bloc comme il a été possible de la faire pour les médians. 

De surcroît, autant il était envisageable de concevoir un kernel par taille de masque lorsqu'il s'agissait de filtres médians car ils ne comportent qu'un seul paramètre, autant cela devient inconcevable pour les filtres de convolution et l'immense variété de paramétrage qu'ils recouvrent. La contrainte de définir les valeurs des coefficients du masque de manière littérale à l'intérieur du kernel doit donc être levée pour permettre de rendre toute leur souplesse aux opérations de convolution.

Parmi les types de mémoire disponibles, nous avons opté pour le stockage des coefficients du masque en mémoire constante (\textit{symbol memory}) en raison de ses performances et du petit volume requis. L'abandon des registres permet aussi d'adopter un style de codage beaucoup plus conventionnel utilisant des structures de contrôle classiques (itérations et tableaux). 

L'augmentation du nombre de pixels traités par chaque thread est alors de nouveau envisageable, puisque l'utilisation de la mémoire constante pour les coefficients libère autant de registres. Il faut cependant organiser les calculs de manière à réduire autant que possible les accès en lecture aux valeurs de l'image en texture ; en d'autres termes, on cherche à exploiter le recouvrement entre positions voisines du masque de convolution et n'effectuer qu'une seule lecture par pixel de l'image pour en distribuer la valeur sur l'ensemble des calculs de convolutions en cours dans le thread. Cela complique quelque peu les expressions des sommes partielles mais réalise l'objectif opérationnel de la mémoire partagée sans en subir le coût ni les contraintes d'accès et permet ainsi d'envisager de meilleures performances.

Multiplier les pixels traités par un même thread impose également de faire un choix sur la forme de ce que l'on appellera un \textit{paquet} de pixels (centraux, par opposition aux pixels des voisinages, même si un pixel a successivement l'un et l'autre des statuts). La contrainte de contiguité des accès en mémoire globale pour la mémorisation des valeurs de sortie font que seule l'organisation \textit{en ligne} des paquets de pixels est bénéfique, bien que n'étant pas celle qui présente systématiquement les recouvrements les plus importants. 

Multiplier les pixels traités par un même thread impose également de faire un choix sur la forme de ce que l'on appellera un \textit{paquet} de pixels (centraux, par opposition aux pixels des voisinages, même si un pixel adopte successivement l'un et l'autre des statuts.

\begin{figure}[ht]
  \centering
  \subfigure[Cas d'un masque de taille 3$\times$3 ($k=1$) où l'on dénombre 6 colonnes centrales, soit 18 pixels de multiplicité maximale 3.]{\includegraphics[width=0.4\linewidth]{Chapters/chapter6/img/convoOverlap1.png}}\quad
\subfigure[Cas d'un masque de taille 5$\times$5 ($k=2$) où l'on dénombre 4 colonnes centrales, soit 20 pixels de multiplicité maximale 5.]{\includegraphics[width=0.5\linewidth]{Chapters/chapter6/img/convoOverlap2.png}}
\caption{Multiplicité des implications des pixels de la zone d'intérêt d'un thread dans les calculs de convolution. Le nombre de calculs dans lequel est impliqué un pixel est inscrit en son centre. Le premier pixel du paquet, ou pixel de base, est repéré par ses coordonnées $(x, y)$ ; le dernier a pour coordonnées $(x+7,y)$}
\label{fig-convo-overlap}
\end{figure}

Une valeur de 8 pixels comme taille des paquets, déterminée expérimentalement, s'est avérée optimale sur les deux types d'architecture GPU et pour toutes les tailles de masques soumis au mesures. Cela signifie que chaque thread conduit simultanément les calculs de convolution attachés à chacun des  8 pixels du paquet qu'il traite. La somme partielle de chaque convolution est mémorisée dans un registre. Sur cette base, on a schématisé à la figure \ref{fig-convo-overlap}, l'implication de chaque pixel de la zone d'intéret d'un thread découlant du recouvrement des 8 positions du masque. Pour chaque pixel, cette implication est figurée par une valeur de \textit{multiplicité} représentant le nombre de convolutions différentes dans lesquelles il est impliqué au sein d'un même thread. Tous les pixels d'une colonne partagent la même multiplicité et chaque pixel étant au moins impliqué dans un des 8 calculs, les valeurs de cette multiplicité varient de 1 à k, si k est le \textit{rayon} du masque tel que $n=2k+1$.

On peut dénombrer globalement les multiplicités comme suit :
\begin{itemize}
\item les $(8-2k)$ colonnes centrales de la zone d'intéret, soient $(8-2k)(2k+1)$ pixels sont impliqués dans $k$ calculs.
\item les paires de colonnes symétriques par rapport au bloc des colonnes centrales précédentes ont leurs $2(2k+1)$ pixels impliqués dans $(k-1-e)$ calculs, si $e$ représente l'éloignement avec le bloc de colonnes centrales.
\item Les deux colonnes extérieures ont ainsi leurs pixels impliqués chacun dans un seul calcul de convolution.   
\end{itemize}

Le listing \ref{lst-convo-8x8pL3} présente pour exemple, le code implémentant ces solutions pour les masques de taille 3$\times$3 et l'ensemble des mesures de performance associées, sur C2070,  est regroupé dans le tableau \ref{tab-convo-8x8p}. Cette implémentation atteint des débits supérieurs aux précédentes, mais aussi et surtout surpasse la solution Nvidia avec une exécution du traitement de référence en 1.21~ms sur GTX280, soit une accélération de plus de 14\%. Le gain au niveau du débit reste modeste car les transferts représentent à eux seuls plus de 72\% du temps total. Le modèle GTX280 traite ainsi 962~MP à la seconde, soit un gain de seulement 1.7\% par rapport à la solution de référence. 
Sur C2070, grâce à une bande passante mémoire supérieure, les débits mesurés peuvent dépasser les 2100~MP/s, pour une convolution 3$\times$3 sur une image de 4096$^2$ pixels. Le traitement de référence quant à lui est effectué en 0.987~ms pour un débit de 1666~MP/s. 

\begin{table}
\centering
{\normalsize
\begin{tabular}{clrrrr}
\toprule
&&\multicolumn{4}{c}{Taille d'image}\\
Masque&&$\mathbf{512\times 512}$&$\mathbf{1024\times 1024}$&$\mathbf{2048\times 2048}$&$\mathbf{4096\times 4096}$\\
\midrule
\multirow{2}*{3$\times$3}& temps exéc. (ms)   & 0.036 & 0.128 & 0.495 & 1.964 \\
                         & débit global (MP/s)& 1425  & 1862  & 2071  & 2138  \\
\midrule
\multirow{2}*{5$\times$5}& temps exéc. (ms)   & 0.069 & 0.253 & {\bf 0.987} & 3.926\\
                         & débit global (MP/s)& 1208  & 1524  & {\bf 1666}  & 1711  \\
\midrule
\multirow{2}*{7$\times$7}& temps exéc. (ms)   & 0.110 & 0.413 & 1.615 & 6.416\\ 
                         & débit global (MP/s)& 1016  & 1237  & 1334  & 1364\\
\bottomrule
\end{tabular}
}  
\caption{Performances des kernels effectuant la convolution non-séparable sur le modèle du listing \ref{lst-convo-8x8pL3}, sur GPU C2070. Le temps d'exécution correspond à la seule exécution du kernel. Le débit global intègre les temps de transfert. Les valeurs en gras correspondent au traitement de référence }
\label{tab-convo-8x8p}
\end{table}

\lstinputlisting[label={lst-convo-8x8pL3},caption={Kernel réalisant la convolution par un masque 3$\times$3 dont les coefficients normalisés sont en mémoire constante.}]{Chapters/chapter6/code/convoGene8x8pL3.cu}

\section{Cas de la convolution séparable}

Dans la pratique, les traitements appliqués aux images par des opérations de convolution à deux dimensions reposent souvent sur des masques présentant une ou plusieurs symétries. Lorsqu'un tel masque $h$ peut s'écrire comme le produit de 2 vecteurs $h_v$ et $h_h$, comme dans l'exemple ci-dessous, alors on dit que la convolution 2D est séparable et peut donc être effectuée en deux opérations de convolution 1D de masques respectifs $h_v$ et $h_h$.

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

Une convolution séparable $n\times n$ est donc moins coûteuse en nombre d'opérations arithmétiques, avec seulement $2n$ paires addition/multiplication par pixel contre $n^2$ pour une convolution non séparable. Cela représente un gain de 60\% du nombre d'opérations pour un masque 5$\times$5 et nous laisse entrevoir des performances supérieures à celles de la convolution non séparable.

Il faut cependant considérer qu'effectuer un traitement en 2 exécutions de kernel(s) consécutives implique de multiplier aussi les écritures en mémoire globale, ce qui a un coût. La plupart des implémentations séquentielles de la convolution séparable utilisent la même fonction pour réaliser les 2 passes, horizontale et verticale, la première mémorisant la transposée de l'image de sortie pour qu'elle soit traitée directement par la seconde passe. Sur GPU, cette solution se heurte aux contraintes de contiguïté dans les accès à la mémoire globale, il faut donc préférer deux kernels distincts : un pour la convolution verticale, l'autre par l'horizontale. La mémorisation de l'image intermédiaire est effectuée en mémoire globale, qui est ensuite recopiée en texture. Nos mesures (tableau \ref{tab-convo-memcpy}) montrent que le coût de la copie en texture est largement compensé par le gain apporté par le cache 2D de la texture pour les lectures des valeurs des pixels.

\begin{table}
\centering
{\normalsize
\begin{tabular}{cr}
\toprule
\textbf{Image}& Temps (ms)\\
\midrule
$\mathbf{512\times 512}$  & 0.029\\
$\mathbf{1024\times 1024}$& 0.101\\
$\mathbf{2048\times 2048}$& 0.387\\
$\mathbf{4096\times 4096}$& 1.533\\
\bottomrule
\end{tabular}
}  
\caption{Coût, en ms, de la copie effectuée entre les deux phases de convolution 1D, sur C2070.}
\label{tab-convo-memcpy}
\end{table}

En revanche, les latences d'accès aux textures ne sont plus compensées par les distributions des valeurs sur plusieurs des calculs menés par un thread pour un paquet de pixels. En effet, dans une convolution séparable, il n'y a pas de recouvrement entre les différentes positions du masque associées aux pixels d'un paquet. Aucun gain n'est donc possible de ce côté et il s'avère même que l'utilisation de la mémoire partagée est ici la solution la plus performante.

Pour chacune des convolutions 1D, la zone d'intérêt d'un bloc de threads ne s'étend que dans une direction et l'on peut donc appliquer une version simplifiée du cadre général d'emploi de la mémoire partagée présenté au paragraphe \ref{sec-bilateral}. Les listings \ref{lst-convo-1Dv} et \ref{lst-convo-1Dh} détaillent la mise en \oe uvre complète des kernels de convolutions verticale et horizontale, pour des paquets de 8 pixels, qui demeure la taille optimale dans le cas séparable. 
 
\lstinputlisting[label={lst-convo-1Dv},caption={Kernel réalisant la convolution verticale 3$\times$1.}]{Chapters/chapter6/code/convoSepShV.cu}
\lstinputlisting[label={lst-convo-1Dh},caption={Kernel réalisant la convolution horizontale 1$\times$3.}]{Chapters/chapter6/code/convoSepShH.cu}

Les temps d'exécution et débits effectifs globaux de cette implémentation sont détaillés dans le tableau \ref{tab-convo-sep} jusqu'à la taille 13$\times$13. Les temps d'excécution des deux kernels étant très voisins, le tableau présente la somme des temps des deux et de la copie mémoire afin de disposer d'une base de comparaison claire avec la convolution non séparable. 
L'analyse des valeurs nous confirme que la complexité réduite de la convolution séparable permet une moindre dépendance à la taille du masque. Elle confirme également que le coût de la copie intérmédiaire n'est pas amorti pour les petites tailles de masque et l'implémentation optimisée de la convolution non séparable demeure plus rapide que la séparable pour le tailles 3$\times$3 et 5$\times$5.

\begin{table}
\centering
{\normalsize
\begin{tabular}{clrrrr}
\toprule
&&\multicolumn{4}{c}{Taille d'image}\\
Masque&&$\mathbf{512\times 512}$&$\mathbf{1024\times 1024}$&$\mathbf{2048\times 2048}$&$\mathbf{4096\times 4096}$\\
\midrule
\multirow{2}*{3$\times$3}& temps exéc. (ms)   & 0.080 & 0.306 & 1.094 & 4.262 \\
                         & débit global (MP/s)& 1150  & 1415  & 1598  & 1654  \\
\midrule
\multirow{2}*{5$\times$5}& temps exéc. (ms)   & 0.087 & 0.333 & 1.191 & 4.631\\ 
                         & débit global (MP/s)& 1116  & 1365  & 1541  & 1596  \\
\midrule
\multirow{2}*{7$\times$7}& temps exéc. (ms)   & 0.095 & 0.333 & 1.260 & 5.000\\ 
                         & débit global (MP/s)& 1079  & 1365  & 1503  & 1542\\
\midrule
\multirow{2}*{9$\times$9}& temps exéc. (ms)   & 0.108 & 0.378 & 1.444 & 5.676\\ 
                         & débit global (MP/s)& 1024  & 1290  & 1410  & 1452\\
\midrule
\multirow{2}*{11$\times$11}& temps exéc. (ms)   & 0.115 & 0.404 & 1.545 & 6.105\\ 
                         & débit global (MP/s)  & 997  & 1250  & 1364  & 1400\\
\midrule
\multirow{2}*{13$\times$13}& temps exéc. (ms)   & 0.126 & 0.468 & 1.722 & 6.736\\ 
                         & débit global (MP/s)  & 957  & 1169  & 1290  & 1330\\
\bottomrule
\end{tabular}
}  
\caption{Performances des kernels effectuant la convolution séparable sur le modèle des listings \ref{lst-convo-1Dv} et \ref{lst-convo-1Dh}, sur GPU C2070. Le temps d'exécution correspond à l'exécution des 2 kernels et de la copie intérmédiaire. Le débit global intègre les temps de transfert.}
\label{tab-convo-sep}
\end{table}

\section{Conclusion}

L'architecture des GPU et le modèle de programmation CUDA permettent d'implémenter efficacement les opérations de convolution, séparable ou non séparable. 
Nous avons transposé les principes appliqués aux filtres médians et montré qu'ils n'étaient pas tous pertinents dans le cas de la convolution. Nous avons cependant proposé des solutions adaptées qui ont permis d'atteindre des performances encore inégalées sur GPU Nvidia avec jusqu'à 2138 millions de pixels traités à la seconde, transferts inclus.
Les expérimentations conduites sur les kernels de convolution tendent également à confirmer dans un cadre plus large ce que les travaux sur les filtres médians avaient fait apparaître : l'usage de la mémoire partagée ne représente souvent pas la solution apportant les meilleure performances. Cela peut cependant être les cas, en particulier lorsque les voisinages des pixels d'un même paquet ne se recouvrent pas, rendant sans objet toute optimisation liée à ces recouvrements, comme la distribution des valeurs sur les calculs multiples.

Conscients du manque de souplesse découlant de l'optimisation de ces kernels et pour que cela ne soit pas un frein à l'utilisation de ces solutions, nous avons enfin proposé une application en ligne qui génère à la demande les codes des kernels médians et de convolution d'après les critères indiqués par l'utilisateur, qui peut alors télécharger un ensemble suffisant et immédiatement fonctionnel comprenant un fichier kernel GPU, un fichier main.c, un Makefile et une image de test. Il est accessible à l'adresse http://info.iut-bm.univ-fcomte.fr/staff/perrot/convomed.
 