Implementação de convolução em GPU

altor:Fernando Paolieri Neto
Data:22 de outubro de 2010

introdução

Todos os filtros lineares podem ser implementados na forma de uma convolução isso a faz uma das mais importantes operações de processamento de imagem, trivialmente são implementados em c usando processamento sequencial com o surgimento das arquiteturas many-cores (existente nas placas de vídeo) houve a possibilidade de torna-los muito mais eficientes, por sua implementação agora ser paralela, e assim prover um ambiente de visualização de imagens em tempo real e aplicações que levariam horas sejam possíveis em segundos.

Arquitetura CUDA

Esta sessão tenta explicar como a arquitetura cuda funciona.

Arquitetura Multicore e Manycore

/media/Attachments/courseIA366F2S2010/fp100356P1/CudaVsGpu.png

Figura 1: Comparação entre GPU e CPU

As arquiteturas existentes hoje na CPU são chamadas arquiteturas multicore pois existem vários núcleos (em torno de 3 em computadores pessoais) e todos eles trabalham de forma independente, isto significa que cada unidade tem sua unidade de controle, suas unidades logicas e aritméticas e um cache muito grande. Estas características fazem ela ser eficiente em operações sequenciais.

No caso da GPU existem vários processadores de fluxo (em torno de 15 em placas de médio custo) cada um independente, com sua unidade de controle, um pequeno cache e em torno de 20 unidades logicas e aritméticas, assim as GPU's são capazes de processão uma quantidade grande de dados em paralelo esta característica faz com que ela seja muito eficiente em operações paralelas, pois são em torno de 300 unidades logicas e aritméticas em paralelo, por causa desse paralelismo ao extremo, esta arquitetura é chamada de arquitetura manycore.

Organização lógica da GPU

/media/Attachments/courseIA366F2S2010/fp100356P1/GridBlockThread.png

Figura 2: Organização lógica da GPU

Cada stream processor executa um conjunto de threads, as threads são organizadas em uma estrutura chamada bloco (que são estruturas bi dimensionais de threads), os blocos são organizados em um grid (estruturas tri dimensionaisde blocos)

Por conta da organização em blocos a GPU é altamente escalável entre diversos tipos de placas de video, pois o no caso de uma GPU com muitos stream processors executará mais blocos simultaneamente e no caso de menos, menos blocos simultaneamente.

Organização da memória

/media/Attachments/courseIA366F2S2010/fp100356P1/MemoryHierarqui.png

Figura 3: Hierarquia de memórias

Todas as threads tem acesso a mesma memoria global que é a memória mais devagar de todas. As threads de mesmo bloco tem acesso a uma memoria chamada de memória compartilhada que é mais rápida que a memoria global porem só é visível no escopo do bloco. E cada thread tem acesso aos seus registradores que tem o menor tempo de acesso de todos.

Estratégias para uso da memória compartilhada

Usualmente é feito apenas um acesso a memória global para carregar o valor para memória compartilhada, a partir daí é feito todo o processamento usando apenas a memória compartilhada e após o processamento é feito o retorno dos dados para a memória global.

A convolução 2D

A convolução 2D é uma técnica que necessita de uma matriz de entrada e uma mascara de convolução, esta mascara tem seu sistema de referência definido (não necessariamente na posição 0,0) uma de suas aplicações é na area de processamento de imagem onde a matriz de entrada representa uma imagem de entrada, e a mascara de convolução é uma matriz retangular (pode ser quadrada) que representa o comportamento do filtro.

Supondo um pedaço de uma imagem:

\begin{vmatrix}
2 & 1 & 4 \\
1 & 4 & 8 \\
8 & 2 & 1
\end{vmatrix}

e a mascara com origem no centro

\begin{vmatrix}
1 & 1 & 1 \\
1 & 1 & 1 \\
1 & 1 & 1
\end{vmatrix}/9

O resultado esperado para o pixel central do pedaço da imagem de entrada é calculado sobrepondo a mascara sobre o pixel que deseja-se obter o resultado. Assim o resultado esperado é: 2 * 1 + 1 * 1 + 4 * 1 + 1 * 1 + 4 * 1 + 8 * 1 + 8 * 2 + 1 * 1 = 3.44

Analisando o resultado obtido, a mascara escolhida faz uma média entre os pixels vizinhos e o atribui para o pixel central.

Tratamento de bordas

No caso da mascara de convolução estiver em uma borda, ela não terá um valor válido de imagem para os coordenadas da borda, neste caso, varias estratégias podem ser adotadas, neste artigo a estratégia adotada é imaginar que a imagem é infinita e periódica, esta estratégia pode ser aplicada apenas usando a fórmula:

novoIndice = (indice + maximoValorIndice) \% maximoValorIndice

Assim aplicando esta fórmula todos os valores serão válidos e ainda a convolução terá a sua relação com a transformada de fourier.

Implementações de convolução

Esta sessão exemplifica como pode ser implementado diversas implementações de convolução começando com uma implementação simples na CPU passando por implementações na GPU.

Utilizando a CPU

A CPU implementa o algoritmo de forma sequencial, ou seja, um processador é responsável por processar todos os pixels, ele processa o primeiro, depois o segundo, e assim por diante.

A formula da convolução normalmente utilizada para implementação em CPU é:

g(i) = (f*h)(i) = \sum^m \sum^n f(i-n,j-m) h(n,m)

A implementação trivial de convolução é dada pelo algoritmo a baixo:

\begin{algorithmic}[1]

    \For{$y_f = 0$ até $H_f$}
        \For{$x_f = 0$ até $W_f$}

            \State $V \leftarrow 0$

            \For{$y_m = 0$ até $H_m$}
                \For{$x_m = 0$ até $W_m$}

                    \State Calcular $id_m$ a partir de $y_m$ e $x_m$
                    \State Calcular $id_f$ baseado em $x_f + x_m$ e $y_f + y_m$

                    \State $V \leftarrow V + f[id_f] m[id_m]$

                \EndFor
            \EndFor

            \State Calcular $id_f$ baseado em $x_f$ e $y_f$
            \State $g[id_f] \leftarrow V$

        \EndFor
    \EndFor



\end{algorithmic}

Algoritmo 1: Convolução implementada na CPU

Utilizando a GPU

No caso da implementação da GPU, funciona da mesma forma que a implementação na CPU, só que neste caso é implementado um kernel que é executado por cada nucleo da GPU e cada nucleo é responsavel por um pixel. Assim o algoritmo é identico ao da CPU, com a diferença de que não é necessário mais um laço para varrer cada pixel da imagem, pois este processamento será feito em paralelo.

Um exemplo de implementação trivial em GPU

\begin{algorithmic}[1]

    \State $x_m \leftarrow$ indice X do pixel que esta thread é responsavel
    \State $y_m \leftarrow$ indice Y do pixel que esta thread é responsavel
    \State $V \leftarrow 0$

    \For{$y_m = 0$ até $H_m$}
        \For{$x_m = 0$ até $W_m$}

            \State Calcular $id_m$ a partir de $y_m$ e $x_m$
            \State Calcular $id_f$ baseado em $x_f + x_m$ e $y_f + y_m$

            \State $V \leftarrow V + f[id_f] m[id_m]$

        \EndFor
    \EndFor

    \State Calcular $id_f$ baseado em $x_f$ e $y_f$
    \State $g[id_f] \leftarrow V$
\end{algorithmic}

Algoritmo 2: Convolução trivial em GPU

Este mesmo algoritmo foi implementado em C, segue o codigo fonte:

 1 #include <cuda.h>
 2 #define BLOCK_SIZE 16
 3 
 4 // in = Imagem de entrada
 5 // out= Imagem de saída
 6 // sizex = Numero de colunas
 7 __global__ void convolution_kernel(float * in, float* out, int sizex,int sizey, float mask[], int msizex,int msizey) {
 8     int x = blockIdx.x * blockDim.x + threadIdx.x;
 9     int y = blockIdx.y * blockDim.y + threadIdx.y;
10 
11     float v = 0.0f;
12 
13     if (x >= sizex || y >= sizey) return;
14 
15     for (int i=0; i<msizey; i++) {
16         int posy = y+i -(msizey/2);
17             posy = (posy+sizey)%sizey;
18         for (int j=0; j<msizex; j++) {
19             int posx = x+j -(msizex/2);
20                 posx = (posx + sizex)%sizex;
21 
22             float maskVal = mask[ (msizey-i-1) * msizex + (msizex-j-1)];
23             if (maskVal){
24                 v +=  maskVal * in[posx + posy * sizex];
25             }
26         }
27     }
28     out[y * sizex + x] = v;
29 }

É possível notar que o código da GPU implementado é bem mais complexo, pois exige muita manipulação de indice, principalmente porque é assumido que a origem da mascara de convolução é sempre no centro.

Utilizando a GPU com memória compartilhada

Como a implementação trivial se baseia apenas no uso da memória global, esta segunda implementação tenta ganhar eficiencia com o uso de memória compartilhada e com a simplificação da convolução, onde na hora que a imagem for carregada para a memória compartilhada ela já é carregada com uma translação, assim o resultado não será diferente do resultado obtido quando a origem da mascara se encontra no centro da mascara. A partir do momento que a mascara é carregada ela é processada como se a mascara tivesse sua origem no canto superior esquerdo.

Assim implementação deste algoritmo consiste de 3 etapas:

  • Carregar os valores da memória global para a memória compartilhada;
  • Processar os valores na memória compartilhada;
  • Gravar o resultado na memória global;

Como no caso da convolução o valor de um pixel é dependente de seus vizinhos (dependendo das dimensões da mascara) então é necessário que sejam carregados mais valores para a memória compartilhada. A estratégia adotada por este algoritmo foi que cada thread carregue mais de um pixel, assim todos os valores necessários para realizar a conta estarão na memória compartilhada

Supondo o uma convolução de uma imagem 10x10 e que cada bloco tenham dimensões de 5x5 e a mascara de convolução seja 3x3, para calcular o valor de cada thread precisam para calcular mais 8 pixels, entre eles pixels a direita e pixels a baixo (por conta da mascara ter a origem no canto acima a esquerda), assim faltaram valores de pixels para calcular a convolução, neste contexto foi criado a umbra. A umbra são pixels que são carregados a mais para que o calculo seja possível.

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraP0.png

Figura 4: Funcionamento do bloco 0,0. Parte vermelha são os pixels cuja saída são de responsabilidade do bloco. Parte laranja representa a umbra do bloco.

O mesmo ocorre com os outros blocos de forma periodica, como mostra as imagens a baixo

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraP1.png

Figura 5: Funcionamento do bloco 1,0. Parte vermelha são os pixels cuja saída são de responsabilidade do bloco. Parte laranja representa a umbra do bloco.

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraP2.png

Figura 6: Funcionamento do bloco 0,1. Parte vermelha são os pixels cuja saída são de responsabilidade do bloco. Parte laranja representa a umbra do bloco.

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraP3.png

Figura 7: Funcionamento do bloco 1,1. Parte vermelha são os pixels cuja saída são de responsabilidade do bloco. Parte laranja representa a umbra do bloco.

Desta forma a imagem na memória compartilhada ficará sempre desta forma:

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraShared.png

Figura 8: Organização dos pixels da imagem na memória compartilhada para qualquer bloco.

Com isso é possivel calcular o valor do pixel sobrepondo uma mascara de convolução e multiplicando a mascara pelos elementos da imagem que ela sobrepõem e depois somando todos os elementos resultantes. Isso deve ser feito para cada thread, segue alguns exemplos de sobreposição:

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraShared00.png

Figura 9: Sobreposição da mascara para o calculo do valor resultante da thread 0,0.

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraShared10.png

Figura 10: Sobreposição da mascara para o calculo do valor resultante da thread 5,0.

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraShared01.png

Figura 11: Sobreposição da mascara para o calculo do valor resultante da thread 0,5.

/media/Attachments/courseIA366F2S2010/fp100356P1/gridUmbraShared11.png

Figura 12: Sobreposição da mascara para o calculo do valor resultante da thread 5,5.

Assim foi implementado um código que implementa a técnica:

 1 #include <string.h>
 2 #include <math.h>
 3 #include <cuda.h>
 4 
 5 #define KerShElem 1024*4
 6 #define KerCoElem 64*64
 7 
 8 #define BlockDimAxis 16
 9 
10 __constant__ float Md[KerCoElem];
11 
12 #define Bx        (blockDim.x*blockIdx.x + threadIdx.x)
13 #define By        (blockDim.y*blockIdx.y + threadIdx.y)
14 #define Bid       (By * Awidth + Bx)
15 #define Swidth    (blockDim.x+Mwidth -1)
16 #define Sheight   (blockDim.y+Mheight-1)
17 
18 // Ad = Imagem de entrada
19 // Aheight = numero de linhas da Imagem de entrada
20 // Awidth  = numero de colunas Imagem de entrada
21 // Md = mascara
22 // Mheight = Numero de linhas da mascara
23 // Mwidth  = Numero de colunas da mascara
24 // Bd = Imagem de saida
25 __global__ void conv2D_kernel(unsigned char *Ad, int Aheight, int Awidth, int Mheight,int Mwidth, float * Bd)
26 {
27     __shared__ unsigned char S[KerShElem];
28 
29     int Ax,Ay;
30     int Sx,Sy,Sid;
31     int i,j;
32 
33     // Etapa 1 - carrega os valores na memória compartilhada transladada
34     for (j=0;j<=Sheight/blockDim.y;j++){
35         for (i=0;i<=Swidth/blockDim.x;i++){
36 
37             Sx = threadIdx.x + i*blockDim.x;
38             Sy = threadIdx.y + j*blockDim.y;
39 
40             // se nao for interessante, continua;
41             if (Sx < Swidth && Sy < Sheight) {
42                 Ax= (Bx-Mwidth /2  +i*blockDim.x+Awidth )% Awidth;
43                 Ay= (By-Mheight/2 +j*blockDim.y+Aheight)% Aheight;
44 
45                 S[Sy*Swidth + Sx] = Ad[Ay * Awidth + Ax];
46             }
47         }
48     }
49     __syncthreads();
50 
51     if ( Bx >= Awidth || By >= Aheight) return;
52 
53 
54     // Etapa 2 - Processa os valores na memória compartilhada
55     float v = 0.0f;
56 
57     Sx = threadIdx.x;
58     Sy = threadIdx.y;
59     Sid = Sx + Sy*Swidth;
60 
61     for (j=0;j<Mheight;j++){
62         for (i=0;i<Mwidth;i++){
63             v += Md[ (Mheight-j-1)*Mwidth + (Mwidth-i-1) ]
64                  * S[Sid + i + j*Swidth];
65         }
66     }
67 
68     // Etapa 3 - Coloca o resultado na memória global
69     Bd[Bid] = v;
70 }

Comparação dos resultados obtidos

O teste foi realizado com o uso de uma imagem de 256 e 256 feita de números aleatórios e como as três implementações que serão testadas fazem o mesmo tratamento de borda, assim é esperado o resultado igual.

CPU: 4.7861328125 ms
GPU: 0.776123046875 ms
GPU (com memória compartilhada): 0.591064453125 ms

É possível notar que a implementação usando memória compartilhada é mais rapida que a implementação mais trivial utilizando GPU e a implementação mais trivial usando GPU teve um ganho significativo com relação a implementação feita com a CPU, este resultado é esperado pelo alto poder de processamento da GPU.

Outros tipos de implementação

Existem outras estratégias de implementação de convolução utilizando a GPU, estas estratégias serão abordadas nesta sessão.

Implementação que uma thread carrega apenas um pixel.

Nesta implementação cada thread inicia seu processamento, carrega seu pixel para a memória compartilhada, e caso ele possa calcular o valor do pixel que ele está associado, ele o processa, caso contrário, ele não realiza tarefa alguma.

Esta implementação é bem problemática, pois vão haver diversas threads que após carregar o pixels na memória compartilhada, não irão mais fazer nada até o termino da execução deste nucleo.

Por este motivo podemos dizer que este algoritimo irá subutilizar as threads de processamento.

Incremental.

Nesta implementação cada thread inicia carregando um pedaço da imagem original para a memória compartilhada, após esta etapa é feito o processamento de um dos elementos do nucleo para todos os pixels da imagem e este processo se repete para cada um dos elementos da mascara de convolução.

Esta implementação possuí um problema simples, que ela necessita de muitas barreiras de sincronismo que pode acabar causando uma serialização na execução do código, tirando isto, é um algoritmo eficiente pois não subutiliza as threads de processamento e nem sobrecarrega a memória compartilhada

Estratégias secundárias

Além de alguns tipos diferentes de implementação existem mais duas estratégias que podem otimizar o processo do calculo da convolução, entre elas são o uso da memória de texturas e o processo de mais de um pixel por thread.

Uso de textura

A biblioteca fornecida pela NVIDIA permite que seja utilizado uma memória especial usada para textura, esta memória possui um cache e fornece a possibilidade de realizar tarefas via hardware como a interpolação bi-linear.

Ela é altamente recomendada pela NVIDIA caso o acesso a memória não seja coalecente (ou seja, cada thread leia uma posição de memória de forma seqüencial) ou seja feito o uso dos recursos fornecidos via hardware como a interpolação bi-linear.

Processamento de mais de um pixel por thread

O processamento de mais de um pixel por thread pode trazer beneficio para a convolução por diminuir o numero de leituras da memória global para a memória compartilhada, se existir um numero menor de leituras, o algoritmo irá ser mais rápido assim esta é uma pratica desejável.

Análise dos resultados

Para realizar a comparação de tempo dos resultados obtidos foi utilizado uma imagem de 256 por 256 e uma mascara de convolução 3 por 3, como cada um dos algoritmos tratou a borda de uma maneira diferente, será realizado um pad de duas linhas acima e duas a baixo da imagem e duas colunas a esquerda e duas a direita e a comparação será feita ignorando as linhas adicionadas para que o tratamento da borda não influencie no resultado.

A tabela a seguir ilustra na primeira coluna se a comparação foi ou não bem sucedida, a segunda mostra o tempo de execução e a tarceira mostra a técnica utilizada.

OK 0.5891113281 ms => GPU (Mais de um pixel por thread e memória compartilhada)
FL 0.6340332031 ms => GPU (Textura, memória compartilhada e 4 pixels por thread)
FL 0.7866210938 ms => GPU (Textura e memória compartilhada)
OK 0.7919921875 ms => GPU (Apenas memória global)
FL 0.8046875000 ms => GPU (Incremental)
FL 0.9101562500 ms => GPU (Textura)
OK 0.9799804688 ms => GPU (Memória compartilhada e 2 pixels por thread)
OK 4.9682617188 ms => CPU

É possível observar que algumas técnicas não deram o resultado esperado, como por exemplo a técnica incremental e as técnicas que fazem uso da textura isto pode se dar ao fato de que a técnica usa um elemento estruturante que não está com a origem da mascara em seu centro.

Mas observando apenas os tempos de execução é possível notar que processar mais de um pixel por vez nem sempre traz um bom resultado, e que o uso de textura é benéfico por conta de seu cache.

Conclusão

Existem dezenas de técnicas capazes de realizar a tarefa da convolução, mas para implementação de um código eficiente em cuda, exige uma implementação com alto grau de paralelismo e um conhecimento de como o hardware funciona, para que possa tomar melhor proveito da hierarquia de memória.

Sempre existe uma maneira de otimizar o código da GPU, mas nem sempre vale a pena, pois cada otimização exige um conhecimento elevado de como o hardware funciona mas a cada dia que passa o hardware fica melhor e menos conhecimento é exegido na hora de otimizar o código.

Referências

[KirkHwu]David B. Kirk and Wen-mei W. Hwu, ``Programming Massively Parallel Processors'', 1st ed. Burlington, MA: Morgan Kaufmann, 2010.
[CUDA]NVIDIA CUDA C Programming Guide Version 3.2. NVIDIA.