Implementação da Convolução Bi-Dimensional em ambiente CUDA

Autor: Adriano Ferreira
Data: 21/10/2010

Resumo

Este artigo trata de implementações da Convolução Bi-Dimensional utilizando CUDA. É abordado o funcionamento dos algoritmos e suas implementações na GPU, com foco nas estratégias utilizadas para fazer uso dos diferentes componentes contidos na placa GPU, especificamente o uso dos diferentes tipos de memória, e na busca de ganhos em desempenho através de um paralelismo eficiente programado na linguagem CUDA-C. Resultados são apresentados assim como uma conclusão sobre os desempenhos gerados pelas implementações.

Introdução

Convolução

A Convolução é uma técnica amplamente utilizada em aplicações de processamento de imagens. Neste sentido, a Convolução caracteriza-se pela filtragem do domínio espacial, sendo este o próprio plano da imagem, ou seja, o conjunto de pontos que compõe a imagem, com a finalidade de alterar as características desta. A filtragem é realizada por uma matriz denominada máscara ou filtro, que possui em seus elementos valores conhecidos como peso ou coeficiente. No processo de convolução, cada ponto f(x,y) da imagem é convoluído e a computação neste ponto dependerá do seu valor e dos valores dos pontos vizinhos. A máscara é então centrada sobre o ponto f(x,y) da imagem, isto é, seu ponto central na horizontal e na vertical sobrepõe o ponto f(x,y) da imagem, e os outros pontos do filtro sobrepõem os pontos vizinhos à este ponto da imagem, como mostra a figura 1.

/media/Attachments/courseIA366F2S2010/as41664_p1/dominioEspacial.PNG

Figura 1: Exemplo de filtragem no domínio espacial. A máscara, de dimensão 3x3, sobrepõe uma região da imagem. O ponto central do filtro alinha-se sobre o pixel a ser convoluído.

Em geral, seleciona-se um filtro com número ímpar de elementos para que seu centro esteja localizado sobre o ponto a ser computado. A máscara deve sofre uma reflexão (ou rotação de 180º) antes de ser aplicada na imagem, como mostra a figura abaixo.

/media/Attachments/courseIA366F2S2010/as41664_p1/filtroReflexao.PNG

Figura 2: Máscara de dimensão 3x3 refletida.

Após a reflexão e centralização do filtro sobre o pixel a ser processado, cada elemento da máscara multiplica o elemento sobreposto correspondente na imagem. Depois de multiplicados, todos os valores são então somados. O resultado gerado será passado ao elemento da imagem de saída de mesma coordenada do ponto convoluído da imagem de entrada, como é ilustrado abaixo:

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicao.PNG

Figura 3: À esquerda, o filtro sobrepõe uma região da imagem. O centro, em destaque, corresponde ao ponto a ser convoluído, os elementos em vermelho, à esquerda, correspondem aos pesos da máscara e os elementos em preto, à direita, correspondem aos valores dos elementos naquela região da imagem. A figura à direita mostra a saída na mesma coordenada do ponto convoluído.

A vizinhança necessária para a convolução em um determinado ponto da imagem depende do tamanho do filtro. No exemplo mostrado, a região sobreposta pela máscara tem dimensão 3x3 e os elementos vizinhos correspondem aos elementos nas coordenadas (x + 1, y), (x - 1, y), (x, y + 1) e (x, y – 1) (vizinhos horizontais e verticais) e às coordenadas (x – 1), (y – 1), (x – 1, y + 1), (x + 1, y – 1) e (x + 1, y + 1) (vizinhos diagonais) da imagem. Sendo assim, a computação executada para a geração do resultado no exemplo acima e considerando I e Y as imagens de entrada e saída, respectivamente, e W o filtro, é: Y[1, 1] = I[0, 0] * W[1, 1] + I[1, 0] * W[0, 1] + I[2, 0] * W[-1, 1] + I[0, 1] * W[1, 0] + I[1, 1] * W[0, 0] + I[2, 1] * W[-1, 0] + I[0, 2] * W[1, -1] + I[1, 2] * W[0, -1] + I[2, 2] * W[-1, -1] = 1 * 1 + 2 * 2 + 3 * 1 + 4 * 0 + 5 * 0 + 6 * 0 + 7 * (-1) + 8 * (-2) + 9 * (-1) = -24. Assim, a convolução é definida matematicamente como a ilustração seguinte.

Onde e . As incógnitas "m" e "n" indicam a largura e altura do filtro.

Um fator importante no mecanismo de filtragem é o tratamento de borda. Quando a máscara centraliza-se sobre pontos das extremidades da imagem, parte do filtro pode se localizar fora mas dimensões da imagem. Nessas condições, cuidados especiais são requeridos e diversas estratégias podem ser usadas para contornar este problema.

/media/Attachments/courseIA366F2S2010/as41664_p1/borda.PNG

Figura 4: Parte da máscara (região pontilhada) localizada fora das dimensões da imagem.

Nos algoritmos presentes neste trabalho, usou-se estratégias como zerar valores que correspondem aos elementos fora do limite da imagem (gerando uma imagem maior), indexação circular (períodica) e sobreposição. Os detalhes serão mostrados nos algoritmos que implementaram as respectivas estratégias.

Considerações sobre o algoritmo em CUDA

Em um ambiente de hardware semelhante à CPU (Central Processing Unit), o algoritmo da convolução é executado seqüencialmente, percorrendo-se a imagem de entrada e aplicando o filtro em um elemento por vez. Na GPU, este processamento é feito paralelamente, possibilitado por sua arquitetura. De uma forma geral, a GPU, composta de muitos núcleos (cores), está focada em paralelizar a execução de múltiplos dados através de uma única instrução. Com a disponibilidade de centenas desses núcleos, a GPU caracteriza-se por um rico sistema multithread, onde é possível alocar uma thread para cada elemento em uma determinada estrutura de dados. Na convolução, pode-se atribuir uma thread por pixel na imagem de entrada, executando a computação de cada pixel paralelamente. Como a quantidade de elementos da imagem pode ser maior que a quantidade de núcleos disponíveis na GPU, o CUDA encarrega-se de alocar as threads para uso dos processadores conforme estes forem liberando suas execuções. Para isto, o CUDA organiza internamente as threads dentro de blocos, que podem se cooperar na sincronização de suas execuções e dividir dados eficientemente a baixo custo. No entanto, threads de diferentes blocos não podem se cooperar. A figura abaixo ilustra esta organização em CUDA.

/media/Attachments/courseIA366F2S2010/as41664_p1/organizacaoThreads.png

Figura 5: Organização das threads em blocos. Blocos e threads são indexados em 2 dimensões, com o objetivo de manter índices similares às matrizes eventualmente processadas.

Como mostrado na figura acima, o CUDA cria uma grid de threads e as organiza dentro de blocos sempre quando é invocado uma execução na GPU. Nota-se no lado esquerdo da figura uma região identificada como host, que indica que uma parte da execução de um programa em CUDA inicia-se não diretamente na GPU (por vezes chamado de device), mas sim no host (CPU). Quando a parte do device precisa ser executada, isto é feito através da invocação de uma função, chamada de kernel, onde efetivamente ocorre o paralelismo. Portanto, quando o kernel é invocado, é criada a grid com blocos de threads, e os dados são então passados do host para o device. Dentro deste, existe uma hierarquia de memória, configurando um escopo global e local para os dados. Quando os dados são copiados do host para o device, eles assumem um escopo global, pois seu armazenamento pode ser feito nas memórias global, textura ou constante. Neste sentido, todos as threads de todos blocos podem acessar estes dados. Por outro lado, pode-se fazer uso da memória compartilhada (shared memory) dada a velocidade de acesso a esta memória, que tem seu escopo por bloco, como mostra a ilustração seguinte.

/media/Attachments/courseIA366F2S2010/as41664_p1/memoriasCuda.png

Figura 6: Tipos de memória no device.

O algoritmo pode sofrer algumas mudanças na implementação dependendo da memória utilizada. As particularidades de cada tipo de memória podem exigir códigos adicionais e algoritmos estruturados para executar nessas memórias. Neste artigo, foram implementados algoritmos usando as memórias global, textura e compartilhada. Ressalta-se que um determinado algoritmo não restringe-se em apenas um tipo de memória em sua implementação, como se pode notar nos códigos que serão apresentados.

Implementação

Memória Global

A imagem de entrada e a máscara são carregados do host para a memória global da GPU. Então todo o processamento é feito nesta memória. Portanto cada thread precisa acessar a memória global sempre que precisar ler ou escrever nesta. Como esta memória é fisicamente longe dos núcleos, os acessos à esta memória tornam-se custosos. Na implementação da Convolução Bi-Dimensional na memória global, o tratamento das bordas da imagem de entrada é feita de maneira periódica, deslocando o índice que não possui coordenada dentro das dimensões da imagem para a outra extremidade da imagem, como ilustra a figura abaixo.

/media/Attachments/courseIA366F2S2010/as41664_p1/periodica.PNG

Figura 7: Indexação circular (periódica) para tratamento das bordas da imagem. As regiões pontilhadas azul e amarela correspondem à elementos fora das dimensões da imagem.

A seguir o código que garante a indexação circular. Soma-se o índice da linha (ou coluna) com a largura (ou altura) da imagem e tira-se o módulo pela própria largura ou altura.

rowIdx = (rowIdx + width) % width;
colIdx = (colIdx + height) % height;

Sendo que width e height são as dimensões da imagem. O centro do filtro é calculado para posicioná-lo sobre o pixel a ser computado, como ilustrado a seguir.

int maskCenterWidth = maskWidth / 2;
int maskCenterHeight = maskHeight / 2;

A rotação da máscara é feita no trecho de código a seguir:

int i = maskWidth - 1 - maskCol;
int j = maskHeight - 1 - maskRow;

Sendo que maskWidth e maskHeight são as dimensões da máscara e maskCol e maskRow são os índices do filtro.

Finalmente, os elementos do filtro e da imagem de entrada sobrepostos são multiplicados e então passados à variável que armazena a somatória das multiplicações.

sum += input[colIdx + width * rowIdx] * mask[i + (maskWidth * j)];

Depois de armazenado a somatória das multiplicações, o valor que representa o pixel convoluído é passado para a imagem de saída.

output[offset] = sum;

Sendo assim, cada thread representada pelo índice "offset" escreve o resultado de sua execução no elemento de mesma coordenada na imagem de saída.

Memória de Textura

A memória de textura caracteriza-se por ser apenas de leitura e pode melhorar o desempenho da aplicação, reduzindo tráfico de memória em acessos a dados que não possuem endereços de memória aritmeticamente consecutivos. A implementação é aparentemente simples, com alguns detalhes de diferença em relação à implementação feita em memória global. No algoritmo da Convolução Bi-Dimensional, o uso da memória de textura foi feito da seguinte maneira: Primeiro, declara-se o objeto que será alocado na memória de textura. Neste algoritmo, apenas a imagem de entrada é passada à memória de textura.

texture<float, 2> texImg;

Dentro do Kernel, é necessário indicar para a memória de textura os índices que garantem a manipulação do objeto alocado. Os índices são passados através da função tex2D(), que tem como argumento o objeto alocado e seus índices.

sum += tex2D(texImg, colIdx, rowIdx) * mask[i + maskWidth * j];

O restante do código dentro do kernel é semelhante ao código elaborado para a memória global. A única diferença é que o kernel não receberá mais como parâmetro a imagem de entrada. A imagem de entrada será passada para a memória de textura através da função cudaBindTexture2D().

cudaBindTexture2D(NULL, texImg, dInput, desc, width, height, width * sizeof(float));

O argumento "desc" é do tipo cudaChannelFormatDesc e deve ser declarado sempre quando trabalha-se com matrizes bi-dimensionais. Após a utilização da memória de textura, faz-se necessário desalocar o objeto da memória, através da função cudaUnbindTexture(), passando como argumento o objeto a ser desalocado.

Memória Compartilhada

A memória compartilhada reside fisicamente na GPU e devido a isso a latência para acessá-la tende a ser bem menor do que as outras memórias. O contexto da shared memory é o bloco. Portanto qualquer thread dentro do bloco divide a memória e pode acessar ou modificar uma variável dentro desta. Neste sentido, as threads dentro do bloco podem se comunicar e colaborar-se entre si. A estratégia para solucionar a Convolução Bi-Dimensional com a memória compartilhada é o carregamento de pedaços da imagem dentro da memória compartilhada. Sendo assim, as threads carregam os valores da memória global para a memória compartilhada, acessando, assim, apenas vez a memória global. Uma vez que os dados são carregados para a shared memory, toda a computação é realizada dentro dela. Depois, o resultado é escrito na memória global. A estratégia ganha complexidade ao tratar as bordas da imagem. A implementação proposta sugere zerar os valores para os índices fora das dimensões da imagem total, como segue na figura e trecho de código abaixo.

/media/Attachments/courseIA366F2S2010/as41664_p1/limiteImagem.PNG

Figura 8: Os elementos em verde representam as bordas da imagem inteira. A região pontilhada representa os elementos fora das dimensões da imagem, que recebem o valor zero.

if (idx < 0 || idy < 0 || idx >= width || idy >= height)
{
      sPartialImage[tx][ty] = 0;
}
else
{
      sPartialImage[tx][ty] = input[offset];
}

 __syncthreads();

O objeto sPartialImage do código acima está contido na memória compartilhada e foi definido anteriormente, no ínicio do kernel (como ilustra o trecho de código abaixo). Nesta sintaxe, a memória compartilhada é também carregada. Depois do carregamento dos dados, é necessário sincronizar as threads de modo que nenhuma outra thread tente ler um dado ainda não carregado. A linha "_syncthreads()" garante a sincronização das threads.

__shared__ float sPartialImage[TILE_WIDTH][TILE_WIDTH];

É importante notar que a dimensão do objeto criado na shared memory é de mesma dimensão do bloco. Além das bordas da imagem total, as bordas de cada bloco deverão ser tratadas, lembrando-se que cada bloco tem seu sub-conjunto de pixels a serem processados e que esses pixels são privados para o bloco. Portanto, as bordas de cada bloco não seriam tratadas se os trechos da imagem de entrada fossem exatamente divididos entre os blocos e depois carregados dentro da memória compartilhada. Nestas condições, o conjunto de pixels que compõe todas as bordas não seriam convoluídos, como podemos mostrar na ilustração seguinte.

/media/Attachments/courseIA366F2S2010/as41664_p1/bordaTotal.PNG

Figura 9: Apenas os pixels em azul são convoluídos se as bordas não forem tratadas. As bordas nas extremidades serão tratadas atribuindo-se zero aos elementos fora das dimensões da imagem de entrada.

Para o tratamento da borda, o algoritmo então usa de uma sobreposição dos blocos, afim de garantir que cada borda de cada bloco tenha os valores necessários para a computação adequada de seus elementos. Assim, as bordas superior, inferior e laterais de um determinado bloco serão sobrepostos por outros blocos com a intenção de copiar os valores das bordas destes outros blocos e convoluí-los. Segue exemplo ilustrativo da abordagem de sobreposição de blocos.

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicaoBlocos1.PNG

Figura 10: Quando a máscara sobrepõe os pixels no limite à direita da sub-imagem no Bloco(0,0), são necessários os pixels contidos no limite à esquerda do Bloco (1,0) para efetuar a convolução e vice-versa.

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicaoBlocos2.PNG

Figura 11: Os pixels do limite inferior do Bloco(0,0) precisarão dos pixels do limite superior do Bloco(0,1).

Para o bloco(0,0), os cantos superior e esquerdo são tratados zerando-se os valores dos elementos para índices fora das dimensões da imagem. Porém o canto inferior e direito são bordas que devem ser tratadas através do uso de dados que estão em outros blocos.

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicaoBlocos3.PNG

Figura 12: Sobreposição Horizontal dos Blocos para tratamento de borda.

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicaoBlocos4.PNG

Figura 13: Sobreposição Vertical dos Blocos para tratamento de borda

Neste sentido, aplica-se a sobreposição entre os blocos para que seja possível adquirir estes valores. A sobreposição cobrirá 2 linhas inferiores e 2 colunas mais a direita do bloco(0,0) através dos blocos(0,1) e (1,0), respectivamente, para satisfazer os elementos de borda de ambos blocos.

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicaoBlocos6.PNG

Figura 14: Sobreposição Horizontal dos Blocos para tratamento de borda

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicaoBlocos8.PNG

Figura 15: Sobreposição Vertical dos Blocos para tratamento de borda

Note que as 2 últimas linhas do bloco(0,0) são idênticas às 2 linhas superiores do bloco(0,1), assim como as 2 primeiras colunas do bloco(1,0) são idênticas às colunas mais a direita do bloco(0,0). Portanto os blocos(0,1) e (1,0) estão sobrepostos nos elementos correspondentes do bloco(0,0). Para a garantia desta sobreposição, subtraí-se duas vezes o raio da máscara das dimensões do bloco, que é proporcional ao número de bordas que deseja-se convoluir. O raio corresponde a metade da dimensão da máscara. Esta, ao centralizar-se sobre um pixel de borda, sobrepõe elementos fora das dimensões da borda em uma dimensão proporcional ao raio do filtro. Multipica-se o raio por dois pois considera-se as bordas em pares, seja na vertical ou na horizontal.

int radius = maskWidth / 2;

int idx = blockIdx.x * (blockDim.x - (2 * radius)) + tx;
int idy = blockIdx.y * (blockDim.y - (2 * radius)) + ty;

Um detalhe importante desta abordagem é que apesar da borda ser pertencente à um determinado bloco, o bloco que faz a sobreposição é que calculará a convolução para aquela borda.

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicaoBlocos5.PNG

Figura 16: Blocos sobrepostos verticalmente.

/media/Attachments/courseIA366F2S2010/as41664_p1/sobreposicaoBlocos7.PNG

Figura 17: Blocos sobrepostos horizontalmente.

Além disso, a sobreposição faz com que um certo conjunto de threads tornam-se subutilizadas. Estas threads correspondem àquelas que foram sobrepostas, dado que na configuração do Grid, dizemos que cada bloco tem TILE_WIDTH * TILE_WIDTH threads. Dada que a indexação proposta pela sobreposição é em função da dimensão do bloco subtraído por duas vezes o raio da máscara, então é necessário calcular o múmero de blocos baseado nessa indexação.

width / (TILE_WIDTH - (2 * radius))

Sabendo-se que width é a largura total da imagem e que TILE_WIDTH é a dimensão do bloco, podemos criar o número adequado de blocos baseado na expressão acima. Assim teremos blocos suficientes para realizar a sobreposição. Outro detalhe é que a indexação dentro da shared memory usa apenas os índices das threads, pois a manipulação é por bloco.

int rowIdx = ty + (maskRow - maskCenterHeight);
int colIdx = tx + (maskCol - maskCenterWidth);

sum +=  sPartialImage[colIdx][rowIdx] * mask[i + maskWidth * j];

__syncthreads();

As variáveis ty e ty representam os índices locais da thread dentro do bloco. Após o cálculo dos índices usados pelo objeto "sPartialImage", as multiplicações são efetuadas, seguida da somatória. Novamente, "_synthreads()" aparece para garantir que nenhuma outra thread pegue um valor desatualizado da variável sum.

Testes de Desempenho

Os testes a seguir são resultados das implementações para as memórias global, textura e compartilhada. Diferentes tamanhos de imagens foram testados bem como variados tamanhos da máscara.

-Memória Global:

Mascara: 3x3

Lado GPU CPU Speedup Compare
177 0.547 3.505 6.408 False
366 1.365 13.803 10.112 False
655 3.312 45.996 13.888 False
1076 8.303 132.556 15.965 False
1766 20.822 366.342 17.594 False
2027 26.864 490.424 18.256 False

Mascara: 5x5

Lado GPU CPU Speedup Compare
177 0.662 9.503 14.358 False
366 1.744 37.775 21.660 False
655 4.611 124.227 26.941 False
1076 11.504 363.796 31.624 False
1766 29.734 981.974 33.025 False
2027 38.524 1303.957 33.848 False

Mascara: 9x9

Lado GPU CPU Speedup Compare
177 1.048 30.539 29.138 False
366 3.085 121.713 39.451 False
655 8.473 398.296 47.008 False
1076 22.348 1160.716 51.938 False
1766 56.275 3147.960 55.939 False
2027 74.761 4173.974 55.831 False

-Memória de Textura:

Mascara: 3x3

Lado GPU CPU Speedup Compare
177 0.544 3.484 6.404 False
366 1.329 13.786 10.372 False
655 3.325 45.865 13.794 False
1076 8.016 134.221 16.744 False
1766 20.400 367.041 17.992 False
2027 26.764 489.019 18.272 False

Mascara: 5x5

Lado GPU CPU Speedup Compare
177 0.631 9.513 15.080 False
366 1.667 37.782 22.664 False
655 4.318 124.403 28.810 False
1076 10.681 364.054 34.085 False
1766 27.538 984.288 35.743 False
2027 36.090 1300.624 36.038 False

Mascara: 9x9

Lado GPU CPU Speedup Compare
177 0.922 30.583 33.172 False
366 2.670 121.720 45.587 False
655 7.577 399.486 52.724 False
1076 19.473 1161.720 59.658 False
1766 50.945 3146.789 61.768 False
2027 66.899 4145.378 61.965 False

-Memória Compartilhada:

Mascara: 3x3

Lado GPU CPU Speedup Compare
177 0.529 3.484 6.588 False
366 1.320 13.800 10.455 False
655 3.259 45.683 14.017 False
1076 7.853 132.655 16.892 False
1766 20.096 366.653 18.245 False
2027 26.147 489.292 18.713 False

Mascara: 5x5

Lado GPU CPU Speedup Compare
177 0.616 9.519 15.451 False
366 1.708 37.803 22.133 False
655 4.462 123.688 27.720 False
1076 11.088 363.344 32.769 False
1766 28.805 985.432 34.210 False
2027 37.439 1302.082 34.779 False

Mascara: 9x9

Lado GPU CPU Speedup Compare
177 1.105 30.545 27.641 False
366 3.446 121.688 35.312 False
655 9.833 397.190 40.394 False
1076 25.473 1159.400 45.515 False
1766 65.876 3150.707 47.828 False
2027 85.945 4152.594 48.317 False

Conclusões

Como mostra os resultados acima apresentados, as 3 abordagens de resolução da Convolução Bi-dimensional apresentaram ganhos significativos de desempenho em relação às implementações executadas na CPU, como mostra as tabelas acima, com variações na velocidade de execução para as diferentes abordagens. A implementação da Convolução Bi-Dimensional pode ser melhorada como um todo através da Convolução separável, que transforma o filtro de 2 dimensões em 2 filtros de 1 dimensão cada. A implementação da sobreposição de blocos, utilizada na memória compartilhada, pode ser melhorada também, através da utilização de outras estratégias como uso de outros padrões de carregamento e acesso aos elementos. Nota-se que a sobreposição de blocos, ao aumentar a máscara, não apresenta muitos ganhos em relação a máscara menor, provavelmente devido ao fato de desperdiçar-se threads e acessos de memória não coalescidos. De uma forma geral, ganhos em desempenho são proporcionais ao tamanho do filtro. A memória global apresentou speedups com ganhos significativos em relação à CPU, para uma máscara 9 x 9 e imagem de lado 2027. A memória de textura alcançou resultado significativo se comparado às memórias global e compartilhada, atingindo um speedup em torno de para uma imagem de lado 1766 e máscara 9 x 9. Esta implementação mostrou-se mais simples do que a de memória compartilhada, sendo uma alternativa de implementação simples e eficiente. De um modo geral, conclui-se que a utilização do CUDA para a implementação da Convolução Bi-Dimensional pode acelerar significativamente sua execução, trazendo ganhos de desempenho muito relevantes através de um paralelização de exucação das aplicações.

Referências

David B. Kirk; Wen-mei W. Hwu, Programming Massively Parallel Processors, 1st ed. Burlington, MA: Morgan Kaufmann, 2010.

Pedrini, Helio; Schwartz, William Robson; Análise de Imagens Digitais, Principios, Algoritmos e Aplicações. Editora Thomson Pioneira, 2008.

Jason Sanders, Edward Kandrot, CUDA by Example: An Introduction to General-Purpose GPU Programming, Addison-Wesley, 2010.