Convolução bidimensional com GPU’s e CUDA: uma abordagem utilizando memória compartilhada

Introdução

A convolução é uma operação de vizinhança utilizada em processamento de imagens para filtragem no domínio espacial, como detecção de contornos, aplicação de filtros passa-alta, passa-baixa etc. Quando falamos em operação de vizinhança, significa que o processamento de cada pixel depende do valor dos pixels que estão ao seu redor dentro de uma vizinhança fixa. Uma imagem f e uma máscara w são dadas como entrada, sendo que o tamanho de w determina a vizinhança utilizada no processamento. O processamento de um pixel da imagem independe do processamento de outros pixels, desta forma, este processamento pode ser realizado de forma simultânea (paralela).

CUDA é uma arquitetura de programação paralela criada pela NVIDIA que possibilita aproveitar o potencial de processamengo da GPU (unidade de processamento gráfico). Nesta arquitetura, diversas threads são executadas ao mesmo tempo, o que permite reduzir consideravelmente o tempo de processamento de programas que podem ser paralelizáveis.

O objetivo do presente trabalho é apresentar a implementação de um algoritmo de convolução utilizando a memória compartilhada da arquitetura CUDA. Para tal, ates será apresentado um algoritmo de convolução utilizando apenas a memória global e destacando as principais etapas do desenvolvimento de um algoritmo paralelo.

Convolução

A convolução é uma operação de vizinhança utilizada em processamento de imagens para filtragem no domínio espacial, operando diretamente nos pixels da imagem. Nesta operação, cada pixel da imagem é processado tomando através de seu valor, vizinhança e uma subimagem (também chamada de filtro, janela ou kernel) com o mesmo tamanho da vizinhança utilizada.

Em cada ponto da imagem, o resultado da filtragem é dado pela soma do produto dos coeficientes da máscara com os valores dos pixels correpondentes na área coberta pela máscara. Ou seja, considerando uma máscara w de tamanho mxn, dado um ponto (x,y) da imagem f, o resultado da filtragem neste ponto é dado por:

Onde e .

A figura figConvolucao ilustra a aplicação do filtro utilizando uma máscara de tamanho 3x3.

/media/Attachments/courseIA366F2S2010/ra033072p1/convolucao.png

Uma consideração importante sobre a implementação da convolução ou qualquer outra operação que considere a vizinhança de um pixel, é definir o tratamento quando o filtro se aproxima da borda da imagem e. Suponha uma máscara de tamanho nxn, a borda da mascara coincide com a borda da imagem quando esta está a uma distância de (n-1)/2 pixels da borda da imagem, desta forma, se o centro da máscara se move em direção à borda da imagem, alguns pixels da máscara estarão fora do plano da imagem.

Existem diversas maneiras de contornar este problema, uma delas, e a mais simples, é limitar o uso da máscara a pixels com distância superior a (n-1)/2 pixels em relação à borda, resultando numa imagem menor que a imagem original. Outra abordagem, para manter o tamanho original da imagem, é trabalhar apenas com os pixels da máscara que estão contidos no domínio da imagem, assim alguns pixels da máscara são descartados durante o processamento da borda da imagem. Pode-se também ampliar o tamanho da imagem, preenchendo o excedente com o valor 0, ou replicando linhas e colunas da imagem de entrada.

Algoritmo de convolução utilizando GPU

Um programa é paralelizável se operações aritméticas podem ser executadas simultaneamente numa estrutura de dados. No caso da convolução, descrita na seção anterior, o processamento de cada pixel da imagem não depende do processamento dos demais pixels, caracterizando um problema que pode ser paralelizado. Cada elemento da imagem resultante g(x,y) é formado pelo elemento f(x,y) da imagem de entrada, sua vizinhança e os pixels da mascara w, como mostra a figura figCalculoPixel, desta forma, os elementos da matriz resultante g podem ser computados simultaneamente, isto é, o valor obtido em cada pixel da imagem resultante independe do valor obtido em outros pixels.

/media/Attachments/courseIA366F2S2010/ra033072p1/calculopixel.png

O modelo de memória utilizado pelo CUDA é mostrado na figura figCUDAmodel. Como pode-se observar neste modelo, existem três tipos de memória: constante (constant), global e compartilhada (shared). Os dados são transmitidos da CPU (host) para a memória global ou constante do dispositivo (device), ou GPU, sendo que as duas memórias são compartilhadas por todas as threads, porém as threads não podem “escrever” na memória constante, ao contrário do que ocorre com a memória global. Além disso, cada bloco possui uma memória compartilhada que é visível apenas às threads do bloco.

/media/Attachments/courseIA366F2S2010/ra033072p1/CUDAmodel.png

Tendo em vista os diferentes tipos de memória da arquitetura CUDA, a seguir serão apresentadas duas formas de implementação para convolução de uma mascara de tamanho nxn numa imagem de tamanho HxW: a primeira utilizando apenas a memória global e limitando o uso da máscara a uma distância de (n-1)/2 pixels da borda da imagem; e a segunda utilizando também a memória compartilhada, ampliando a imagem de entrada e preenchendo com o valor 0 os pixels excedentes à borda da imagem.

Utilizando memória global

Em CUDA, para que o processamento ocorra em paralelo, devemos chamar, dentro do programa main(), uma função kernel que especifica o código que será executado por todas as threads durante a fase paralela. Porém, antes de chamar esta função, devemos definir no host a imagem que será processada (f) e a máscara (w), além da imagem de saída (g).

Nesta implementação, consideraremos que a imagem tem tamanho HxW, a máscara tem tamanho nxn e a imagem resultante terá o mesmo tamanho da imagem original, porém apenas a área de (H-(n-1)/2)x(W-(n-1)/2) pixels da imagem original serão tratados (isto é, apenas o interior da imagem, sem a borda de (n-1)/2 pixels) e a borda da imagem resultante terá valor 0.

As matrizes f,w e g serão as matrizes utilizadas para a convolução utilizando a GPU, desta forma, estas matrizes precisam ser copiadas para o dispositivo. A alocação da memória global do dispositivo é realizada através do comando:

cudaMalloc((void**)&fd, fd.size);
...
cudaFree(fd);

sendo fd um ponteiro para a memória global, onde será alocada a matriz f e fd.size é o tamanho da memória a ser alocada. Neste caso, como f é uma imagem do tipo float, então fd.size = H*W*sizeof(float). Após o processamento paralelo, o comando cudaFree() é executado para liberar o espaço reservado para a matriz f no dispositivo.

Com o espaço para as matrizes f,w g alocados no dispositivo, a função cudaMemcpy() é utilizada para transferir os dados da CPU para a GPU:

cudaMemcpy(fd, f, fd.size, cudaMemcpyHostToDevice);

com esta função, a matriz f é copiada para a matriz fd do dispositivo. O comando cudaMemcpyHostToDevice define o fluxo da operação, isto é, estamos copiando dados da CPU (host) para o dispositivo (device). Após o processamento na GPU, o caminho inverso é adotado através da função

cudaMemcpy(g, gd, gd.size, cudaMemcpyDeviceToHost);

isto é, a matriz resultante da operação na GPU, matriz gd, é copiada para a matriz g do host (cudaMemcpyDeviceToHost).

Quando um kernel é chamado pela função main, dois parâmetros são setados: dimensão do grid e do bloco (block). Estes parâmetros são declarados como variáveis do tipo dim3:

dim3 gridDim(W/TILE_WIDTH+1,H/TILE_WIDTH+1);
dim3 blockDim(TILE_WIDTH, TILE_WIDTH);

Neste caso, TILE_WIDTH é uma constante de tamanho 16, esta constante deve-se ao fato que o processamento no bloco é mais eficiente quando o número de threads é múltiplo de 32 (devido ao processamento de wraps), porém o tamanho do bloco é limitado a 512 threads, sendo assim, 16*16=256 é o maior bloco quadro que não ultrapassa o limite de 512 threads pro bloco.

As threads são organizadas em dois níveis: bloco e grid. Cada grid consiste num grupo de um ou mais blocos de threads, sendo que o grid tem no máximo 3 dimensões e o bloco 3. Assim, variável dimBlock descreve o número de threads em cada bloco e dimGrid descreve o número de blocos utilizados. Deve-se atentar ao fato que o número de blocos deve ser suficiente para processar todos os pixels da imagem, desta forma somamos uma unidade em gridDim para o caso em que a dimensão da imagem não é múltipla de TILE_WIDTH.

Finalmente, a chamada do kernel é realizada através da função

Convolucaokernel2D<<< dimGrid , dimBlock>>>(fd,wd,gd,H,W,n);

A função kernel define o código que será executado em paralelo e é identificada pela palavra chave “__global__” em frente à declaração de Convolucaokernel2D(). Esta palavra chave identifica que a função declarada é uma função kernel do CUDA:

__global__ void Convolucaokernel2D(float *fd, float *gd, float *w, int H, int W, int n)

Dentro do kernel, todas as threads executam a função ao mesmo tempo e, para serem identificadas, as variáveis blockIdx.x (blockIdx.y) e threadIdx.x (threadIdx.y) fornecidas pelo CUDA especificam as coordenadas da thread. A figura figIdx mostra esta organização para blockDim(3, 3) e gridDim(4,3), onde cada quadrado pequeno representa uma thread. A thread destacada em vermelho tem coordenadas threadIdx.x = 2, threadIdx.y = 0 , blockIdx.x =1 e blockIdx.y =1.

/media/Attachments/courseIA366F2S2010/ra033072p1/idx.png

Sabendo a posição de cada thread, deve-se definir a função que todas as threads irão executar. No caso do algoritmo de convolução proposto, cada thread será responsável por calcular o valor final de um pixel na matriz resultante gd. Suponha que a imagem fd seja de tamanho H=8 e W=10, isto é, tem 8 linhas e 10 colunas, e blockDim(3, 3) e gridDim(4,3). Neste caso, a thread localizada nas coordenadas threadIdx.x = 0, threadIdx.y = 0 , blockIdx.x =0 e blockIdx.y =0 irá processar a posição 0 do vetor fd. De forma análoga, a thread localizada nas coordenadas threadIdx.x = 1, threadIdx.y = 0 , blockIdx.x =0 e blockIdx.y =0 irá processar a posição 1 do vetor fd e assim sucessivamente, como mostra a figura figProcfd.

/media/Attachments/courseIA366F2S2010/ra033072p1/procfd.png

Deve-se então encontrar uma relação entre as posições das threads e as posições dos valores a serem processados no vetor fd e armazenados no vetor gd.

De forma genérica,

x = blockDim.x*blockIdx.x + threadIdx.x;
y = blockDim.y*blockIdx.y + threadIdx.y;

nos fornece a coluna(x) e a linha(y) da matriz fd que a thread blockIdx.x, blockIdx.y, blockIdx.x, blockIdx.y processa. Porém, como a matriz fd está vetorizada, o índice que devemos acessar o vetor fd e o vetor gd na memória global é dado por:

idx = y*W + x;

O trecho de código abaixo mostra a função kernel

__global__ void Convolucaokernel2D(float *fd, float *gd, float *wd, int H, int W, int n)
{
    int i;
    int j;
    int t;
    int k;
    int x = blockDim.x*blockIdx.x + threadIdx.x;
    int y = blockDim.y*blockIdx.y + threadIdx.y;
    int idx = y*W + x;
    float Pvalue = 0;


    t= n*n-1;
    k = n/2;

    if (x >= W-k|| y >= H-k|| x <= k-1 || y <= k-1) return;
    for(i=0; i <= n-1; i++){
        for(j=0; j<=n-1; j++){
           Pvalue += wd[t-((j)*n +i)]*fd[(idx-k+i)+ W*(j-k)];
        }
    }
    gd[idx]=Pvalue;

}

A condição “if” do kernel testa se o índice calculado pertence à borda da imagem; neste caso, a thread não calcula Pvalue. Os dois “loops” seguintes servem para percorrer todos os valores da máscara para que estes sejam multiplicados pelos valores correspondentes no vetor fd e o resultado estocado na variável Pvalue. Os índices de wd dentro do loop são calculados de forma análoga ao idx, tomando o cuidado de inverter a máscara (por isso calculamos t menos o indice de wd, onde t é o índice máximo da máscara vetorizada). O índice de fd também precisa variar de acordo com a posição na máscara por isso um novo cálculo de índice é introduzido.

Finalmente, atribuímos à posição processada de gd o valor de Pvalue. Esta última linha poderia ser substituída fazendo a alteração:

gd[idx] += wd[t-((j)*n +i)]*fd[(idx-k+i)+ W*(j-k)];

entretanto, esta alteração supõe a atualização de valores na memória global, assim a utilização de Pvalue reduz os acessos a esta memória.

Utilizando memória compartilhada

Se observarmos o problema da convolução, podemos notar que para o processamento de pixels vizinhos na imagem, são utilizados alguns pixels em comum. A figura figPixelsComum mostra dois pixels vizinhos e, em destaque, os pixels que serão utilizados em comum para o processamento de cada um deles no caso de uma máscara de tamanho 3x3.

/media/Attachments/courseIA366F2S2010/ra033072p1/pixelscomum.png

Cada entrada do vetor fd será acessado por 3 ou mais threads (no caso da máscara 3x3, o número de acessos chega a ate 8 vezes, número que cresce de acordo com o tamanho da máscara). Mara melhorar o tempo de acesso à memória, esta segunda implementação introduz o uso da memória compartilhada, permitindo que cada bloco tenha acesso a uma parte do vetor fd.

Como cada memória compartilhada é visível apenas para as threads no respectivo bloco e a convolução é uma operação que considera os vizinhos dos pixels processados, cada memória compartilhada deve conter os pixels que serão processados pelas threads do bloco mais os pixels vizinhos, além disso, um mesmo pixel pode ser armazenado em memórias compartilhadas diferentes. Suponha que a imagem a ser processada, por uma máscara de tamanho 3x3, seja de tamanho 8x10 e queremos processar 9 pixels em cada bloco, como mostra a figura figDivisaoBloco. Para que estes pixels sejam processados utilizando a memória compartilhada, os pixels ao redor do bloco processado também precisam ser armazenados na memória, além disso, existem pixels que serão armazenados em mais de uma memória compartilhada. Se o tamanho da memória compartilhada é de 5x5 pixels, como mostra a figura, então a quantidade de pixels processada será de (5-2)x(5-2), ou 3x3. De maneira geral, se a memória compartilhada tem tamanho TILE_WIDTH x TILE_WIDTH, então a quantidade de pixels a ser processada por aquele bloco será (TILE_WIDTH-2n) x (TILE_WIDTH-2n), onde n é o tamanho da máscara utilizada, isto é, deve-se armazenar na memória compartilhada não apenas os pixels a serem processados, como também sua vizinhança.

/media/Attachments/courseIA366F2S2010/ra033072p1/divisaobloco.png

Desta maneira, considerando que cada bloco processa apenas (TILE_WIDTH-2n) x (TILE_WIDTH-2n) pixels, as dimensões de gridDim e blockDim são dadas por:

dim3 gridDim(gridW, gridH)
dim3 blockDim(TILE_WIDTH, TILE_WIDTH)

onde

gridW = W/(TILE_WIDTH-2*n) +1

e

gridH = H/(TILE_WIDTH-2*n) +1

Assim, o kernel é chamado através da função

Convolucaokernel2D<<<gridDim, blockDim>>>(fd, gd, wd, H, W, n);

Dentro da função kernel, temos blocos de tamanho TILE_WIDTH x TILE_WIDTH, onde cada thread do bloco é identificada por threadIdx e blockIdx. Calculando o índice idx, como no algoritmo anterior, teríamos a situação apresentada pela figura figBloco. As threads do bloco (1,0), por exemplo, não acessam todos os pixels de cor azul, pois estes serão alocados na memória compartilhada do bloco (0,0). Para que os pixels sejam corretamente alocados nas memórias compartilhadas, a posição dos blocos em relação à imagem precisa ser transladada e os blocos precisam se sobrepor, formando a situação ilustrada pela figura figDivisaoBloco.

/media/Attachments/courseIA366F2S2010/ra033072p1/bloco.png

A relação entre os índices da imagem vetorizada fd e o índice dos blocos é dada por:

int x = blockIdx.x*(TILE_WIDTH - 2*k)+threadIdx.x - k;
int y = blockIdx.y*(TILE_WIDTH - 2*k)+threadIdx.y - k;
int idx = y*W + x;

onde k = (n-1)/2.

Agora, a thread com coordenadas threadIdx.x = 1, threadIdx.y = 1 , blockIdx.x =0 e blockIdx.y =0, com uma máscara de tamanho 3x3, irá acessar o pixel da posição 0 do vetor fd.

A alocação da memória compartilhada é feita através da função

__shared__ float Sh[TILE_WIDTH][TILE_WIDTH];

Isto é, alocamos um espaço na memória compartilhada de tamanho TILE_WIDTH x TILE_WIDTH. Os valores dos pixels são passados para a memória compartilhada através de

if (x < 0 || y < 0 || x >= W || y >= H){
    Sh[threadIdx.y][threadIdx.x] = 0;}
else{
    Sh[threadIdx.y][threadIdx.x] = fd[idx];
}
__syncthreads();

A condição “if” trata quando o indice do loco está associado a uma posição de fd que não existe, como no caso da thread com coordenadas threadIdx.x = 0, threadIdx.y = 0 , blockIdx.x =0 e blockIdx.y =0, com uma máscara de tamanho 3x3, tentará acessar a posição –W-1 do vetor fd. Neste caso, o valor será atribuído a 0. Cada elemento do vetor fd é associado a uma ou mais memórias compartilhadas e tais elementos serão acessados através da variável Sh. Ao final, o parâmetro “__syncthreads()” é utilizado para sincronizar todas as threads.

O trecho de código abaixo mostra o kernel da função Convolucaokernel2D() utilizando memória compartilhada:

__global__ void Convolucaokernel2D(float *fd, float *gd, float *wd, int H, int W, int n)
{

    int t = n*n;
    int k = n/2;
    float Pvalue= 0;

    __shared__ float Sh[TILE_WIDTH][TILE_WIDTH];


    //Indices no bloco deslocado
    int x = blockIdx.x*(TILE_WIDTH - 2*k)+threadIdx.x - k;
    int y = blockIdx.y*(TILE_WIDTH - 2*k)+threadIdx.y - k;
    int idx = y*W + x;

    //Bordas são consideradas iguais a 0
    if (x < 0 || y < 0 || x >= W || y >= H){
        Sh[threadIdx.y][threadIdx.x] = 0;}
    else{
        Sh[threadIdx.y][threadIdx.x] = fd[idx];
    }


    __syncthreads();

    if(x < 0 || y < 0 || x >= W || y >= H ||
       threadIdx.x < k || threadIdx.x >= TILE_WIDTH - k ||
       threadIdx.y < k || threadIdx.y >= TILE_WIDTH - k)
       return;


    for(int i=0; i <= n-1; i++){
        for(int j=0; j<=n-1; j++){

            Pvalue+= wd[t-(i*n +j)-1]*Sh[threadIdx.y+i-k][threadIdx.x+j-k];
        }
    }

    gd[idx] = Pvalue;

}

Nesta implementação, Pvalue é obtido através dos valores de wd e Sh, sendo o último um acesso à memória compartilhada. Finalmente o valor de Pvalue é escrito no vetor resultante gd.

Avaliação de Desempenho

As tabelas abaixo mostram o desempenho dos algoritmos utilizando uma máscara 3x3 e uma máscara 5x5 em imagens de tamanho 177 x 177, 1076 x 1076 e 1577 x 1577. Os testes forem executados na GPU e comparados aos resultados obtidos através função iaconv(), pertencente à biblioteca ia636 do python. Para cada tamanho de imagem, a tabela mostra o tempo de execução na CPU, na GPU, a taxa de speedup, sendo esta dada pela razão (tempo de processamento na CPU) (tempo de processamento na GPU). Além disso, a ultima coluna da tabela mostra se o resultado obtido pelo processamento na GPU é próximo ao resultado gerado pela função iaconv().

Tabela 1: Utilizando memória gobal e máscara de tamanho 3x3

Lado GPU CPU Speedup Compare
177 0.579 1.728 2.985 True
1076 9.166 67.298 7.342 True
1577 18.915 145.031 7.668 True

Tabela 2: Utilizando memória compartilhada e máscara de tamanho 3x3

Lado GPU CPU Speedup Compare
177 0.577 1.728 2.995 True
1076 9.076 65.840 7.254 True
1577 18.675 142.917 7.653 True

Conclusão

O presente trabalho apresentou dois algoritmos para convolução utilizando processaento paralelo e arquitetura CUDA. A primeira abordagem tratou da utlização da memória global, enquanto a segunda tratou da utilização da memória compartilhada. Os resultados mostraram que ao tempo de execução dos dois algoritmos foi menor que o tempo de execução na CPU, entretanto, a abordagem utilizando memória compartilhada mostrou-se mais vantajosa, com taxas de speedup speriores às mesmas taxas obtidas através do primeiro algoritmo. Além disso, os resultados obtidos motivam a utilização de algoritmos paralelos em problemas de processamento de imagens, como exemplo a convolção.