IMPLEMENTAÇÃO DE CONVOLUÇÃO BIDIMENSIONAL DE IMAGENS UTILIZANDO PROCESSAMENTO PARALELO ATRAVÉS DA LINGUAGEM CUDA

Wendell Fioravante da Silva Diniz

Departamento de Mecânica Computacional

Faculdade de Engenharia Mecânica

Universidade Estadual de Campinas

e-mail: wdiniz[at]fem[dot]unicamp[dot]br

Resumo: Dentre o domínio do Processamento de Imagens, existe uma ampla variedade de técnicas que se utilizam do processo da convolução, destacadamente os processos de filtragem. Os algoritmos clássicos que aplicam a convolução são relativamente fáceis de implementar computacionalmente, porém seu custo computacional cresce bastante à medida que a máscara utilizada aumenta. O uso de técnicas tais como processamento paralelo visa minimizar este custo. Neste texto, discute-se a implementação de uma versão da convolução utilizando processamento paralelo, especificamente utilizando a linguagem CUDA-C. Será apresentada a metodologia utilizada bem como os resultados alcançados.

1.INTRODUÇÃO

A operação de convolução pode ser considerada uma das mais importantes técnicas para o processamento de imagens, haja vista a variada gama de soluções em que pode ser utilizada. Pode-se definir a convolução como a filtragem de uma imagem de forma que o valor de cada pixel seja determinado por uma média ponderada dos valores dos pixels vizinhos, sendo o valor da ponderação determinado por um operador matricial. Os valores do operador definem o tipo de filtragem a que a imagem será submetida.

Matematicamente, a convolução pode ser definida como a ponderação de infinitas cópias invertidas de uma função sobre outra, sobre um deslocamento qualquer, ou seja:

De forma geral, as imagens digitais são estruturas de dados discretos e finitos. Em termos práticos, a convolução será limitada pelas dimensões da imagem e da máscara. Faz-se necessário então, recorrer à forma discreta da operação de convolução, que se define como:

Outra característica das imagens digitais é que elas são estruturas de duas dimensões. Faz-se necessário definir a convolução bidimensional. Assim, substituindo por , representando a imagem de saída, por , para a imagem de entrada e para o operador e fazendo os índices dos pixels das imagens nas linhas e colunas respectivamente e , e e os índices das linhas e colunas do operador, a convolução bidimensional de uma imagem fica definida desta forma:

O processo tradicional consiste em varrer cada pixel da imagem, sobrepondo o centro da máscara ao pixel corrente e realizando o somatório dos pixels vizinhos multiplicados pelos valores correspondentes sobrepostos à máscara e escrevendo o resultado no pixel da resposta.

/media/Attachments/courseIA366F2S2010/wd101245_artigo/fig1.png

Figura 1: Exemplo de cálculo da convolução de uma máscara 3x3 em uma imagem 10x10. Considerando o primeiro pixel da imagem como corrente e origem da máscara no primeiro elemento.

Uma consideração importante a ser levado em conta é o tratamento da borda da imagem. Uma vez que a convolução originalmente é definida como uma operação periódica num domínio infinito, ao operar em um domínio discreto, é necessário tratar esse problema. Várias abordagens podem ser utilizadas, sendo a mais comum, considerar a imagem infinita sendo que as regiões fora da borda tem valor zero. Considerando o exemplo anterior utilizando esta abordagem, o valor resultante do pixel será 15.

Observando a operação, nota-se que o número de operações aritméticas (multiplicação-acumulação) será igual ao número de elementos da máscara multiplicado pelo número de pixels da imagem. As estratégias de processamento paralelo visam diminuir o tempo de processamento calculando vários pixels simultaneamente, reduzindo o número de passos para a varredura da imagem de entrada.

2.CARACTERÍSTICAS DE IMPLEMENTAÇÃO EM CUDA

A linguagem CUDA (Compute Unified Device Architecture) surgiu em meados de 2006 como uma alternativa para o desenvolvimento de algoritmos paralelos. A linguagem consiste predominantemente em um conjunto de extensões da linguagem C que permitem utilizar os recursos de uma unidade de processamento gráfico (GPU) para aumentar sensivelmente o desempenho de programas através da paralelização do processamento. A GPU então, passa a se comportar como um poderoso co-processador, auxiliando a unidade de processamento central (CPU) em tarefas que exigem operações massivas sobre grandes porções de dados. A GPU tem seu próprio espaço de memória, coletivamente chamado de device memory. Mais adiante, veremos como este espaço de memória está subdividido.

O modelo adotado em CUDA permite uma abstração do hardware da GPU através do mapeamento dos dados a serem processados em unidades de processamento mínimas chamadas threads. As threads são agrupadas em unidades maiores denominadas blocos e os blocos, por sua vez, são organizados em uma estrutura chamada grid. Cada thread executa um bloco de programa sobre um determinado conjunto de dados. Este bloco de programa recebe o nome de kernel. O trabalho de programação em CUDA consiste em mapear o espaço de dados que cada bloco deve processar. Cada thread realiza a operação definida no kernel sobre esse conjunto de dados de forma independente. Em um dado instante, várias threads diferentes de blocos distintos estarão em execução de forma assíncrona. O hardware se encarrega do agendamento de execução dos vários blocos nas unidades de processamento. Como esse processo é transparente ao desenvolvedor, a forma como é feito o mapeamento entre a abstração em software e hardware não será abordada neste trabalho. Na camada de software, o programador deve definir a configuração na qual o kernel irá rodar através de uma sintaxe específica. Nessa chamada, é informado o número de blocos que o grid irá conter, assim como o número de threads em cada bloco. Detalhes sobre esta divisão e sua chamada serão vistos mais adiante.

Na abordagem utilizada pela linguagem CUDA, o código é divido em duas partes, o host code e o device code. O host code consiste na parte “serial” do programa, que irá ser executado pela CPU. As porções de código que serão executadas na GPU são o que chamados de device code. Quando se faz uma chamada a uma rotina em device code, antes de mais nada, é necessário preparar os dados a serem executados. Esta preparação consiste na transferência destes dados para o espaço de memória da GPU, para que esta possa processá-los. A forma como esses dados serão utilizados pelo device é um aspecto importante que deve ser observado, em virtude das características da GPU, para que bons resultados possam ser alcançados. A GPU possui uma hierarquia de memória bastante particular. Num código CUDA, a memória RAM do computador recebe o nome de host memory e é tratado como um espaço de acesso exclusivo do host code. A GPU não tem acesso a esses dados. Então, toda operação na GPU depende de uma transferência dos dados na host memory para a device memory.

A device memory recebe uma nomenclatura de acordo com certas finalidades determinadas pela arquitetura de hardware desenvolvida pelo fabricante. O principal espaço de memória é a chamada memória global (global memory). A memória global é o maior espaço de memória disponível no dispositivo, normalmente conhecida como VRAM (Video RAM). Num primeiro momento, é para onde devem ser transferidos os dados que serão processados. A memória global possui dois espaços menores, chamados de memória de textura (texture memory) e memória de constantes (constant memory). Dados residentes nestes espaços de memória estão fisicamente na memória global, porém o acesso a elas se dá através de cache, então a velocidade de acesso a esses espaços é significativamente maior que um acesso direto à memória global. Dados nesses espaços de memória podem ser acessados, para leitura e escrita, por todos os blocos. Existem ainda, dois importantes espaços de memória que podem ser utilizados. Esses espaços são bem menores que a memória global, mas têm velocidades de acesso muito maiores, garantindo tempos reduzidos de transferência. O primeiro espaço são os registradores. Este espaço geralmente é utilizado para armazenar valores temporários, como por exemplo variáveis auxiliares, acumuladores e outros. Este espaço é definido para cada thread, ou seja, os dados armazenados ali não podem ser acessados por outras threads do bloco nem por threads de outros blocos. O outro espaço recebe o nome de memória compartilhada (shared memory). Como seu nome indica, trata-se de um espaço de memória compartilhado entre as threads de um mesmo bloco. Ou seja, este espaço é utilizado para cooperação e/ou comunicação entre as threads de um mesmo bloco. Como trata-se de memória on chip, seu uso é extremamente importante para alcançar bom desempenho.

As threads e blocos recebem índices únicos que os identificam. Esses índices podem ser consultados através de variáveis automáticas. As threads em um bloco podem assumir configurações em uma, duas ou três dimensões e os blocos em um grid podem assumir uma ou duas dimensões. A forma como serão divididas as threads e blocos é determinante para a disciplina de acesso aos dados, uma vez que seus índices são usados para compor o índice que acessará os ponteiros para os dados.

A seguir, serão apresentadas algumas propostas de implementação do algoritmo de convolução bidimensional em sua forma paralela, utilizando os recursos da linguagem CUDA.

3.PROPOSTAS DE IMPLEMENTAÇÃO

3.1.PRIMEIRA ABORDAGEM: OPERAÇÃO EM MEMÓRIA GLOBAL

Esta primeira abordagem utiliza apenas a memória global e pode ser considerada a forma mais trivial de implementação. A imagem de entrada é dividida em blocos contíguos de 16 por 16 pixels, tantos quantos forem necessários para cobrir toda a extensão da imagem nas duas dimensões. Cada thread processa um pixel da imagem de saída, portanto, o bloco será configurado em duas dimensões, assim como o grid. Desta forma, cada bloco será responsável por 256 pixels.

A seguir, a ideia da implementação em pseudo-código. As variáveis e representam os índices para acesso calculados a partir dos índices da thread e bloco. A imagem de entrada é representada por e a de saída por . A máscara de filtragem é representada por .

Algoritmo 1: Implementação em Memória Global

Esta implementação é praticamente uma transcrição direta da forma serial, realizando a inversão da máscara e o resultado de cada pixel calculado como a ponderação da imagem de entrada sobre a máscara utilizada. Porém, em vez de se calcular cada pixel a cada varredura do laço, estes são calculados simultaneamente em blocos. Neste ponto, já é possível observar uma melhora no desempenho, visto que a arquitetura de memória da GPU é bastante otimizada em relação às CPUs sequenciais, devido a sua disciplina de acesso por caches e largura de dados maiores. Porém, devido a própria arquitetura da GPU, existem memórias com velocidades ainda maiores que a própria memória global. O uso destas estruturas garante um desempenho ainda melhor, porém o código correspondente torna-se mais complexo. Implementações que usam estes espaços de forma correta tem o desempenho consideravelmente melhorado. Porém, a própria tendência da evolução da arquitetura do hardware, indica que no futuro, o uso destes espaços de memória se dará de forma transparente.

Uma variante deste código, utiliza a memória de textura para armazenamento da imagem de entrada. Os dados continuam residentes na memória global, porem como o acesso à memória de textura se dá por um cache bastante otimizado, inclusive com endereçamento em 2D por hardware, a velocidade de acesso melhora consideravelmente. A memória de textura também apresenta algumas vantagens, como por exemplo o tratamento da borda. Por padrão, a memória de textura atua num modo chamado clamping, em que os valores acessados além do limite da imagem retornam o valor do pixel mais próximo, ou seja, é realizada a replicação da borda. A textura também pode ser configurada para trabalhar no modo wrap, que considera a imagem como uma parede de ladrilhos em que a imagem de entrada se repete indefinidamente. Neste modo, os índices de memória são normalizados para valores entre 0 e 1.

3.2.SEGUNDA ABORDAGEM: OPERAÇÃO EM MEMÓRIA COMPARTILHADA COM LEITURA ATRAVÉS DE TEXTURA E CONSTANTE

Esta segunda abordagem, faz o uso da memória compartilhada para diminuir o número de acessos à memória global, aumentando a velocidade de execução do código. Para a implementação desta estratégia, foram feitas algumas considerações. A dimensão máxima da máscara foi fixada em um raio de 16, ou seja, um quadro de 21 por 21 elementos. A memória compartilhada foi declarada para conter os dados necessários para processar este caso máximo. Considerando a superposição da máscara, será necessário carregar um bloco de dados de tamanho igual a DIMENSÃO DO BLOCO + 2 x RAIO DA MÁSCARA em cada sentido, conforme exemplificado na figura abaixo. A região em verde mostra os pixels do primeiro bloco. A região azul, mostra a vizinhança necessária para o processamento dos pixels da borda e a região em vermelho representa a máscara.

/media/Attachments/courseIA366F2S2010/wd101245_artigo/fig2.png

Figura 2: Sobreposição da máscara ao conjunto de dados mostrando os blocos úteis e a borda

A execução do kernel, ficou então dividida em duas fases. Na primeira fase, as threads cooperam para carregar os dados para execução. No pior caso (kernel de raio 16), cada thread carrega 9 dados da memória global para a compartilhada. Este carregamento é feito através de dois laços, um em cada dimensão, que fazem com que uma janela com as dimensões do bloco percorra a região de interesse da imagem, carregando os dados.

A figura abaixo exemplifica a técnica, supondo um bloco de processamento de lado 3 e uma máscara de raio 1. Os blocos em azul escuro mostram os pixels que serão processados. Os dados em azul claro representam a vizinhança que deve ser carregada para o processamento dos pixels da borda. O quadro vermelho representa a janela de carregamento. Nesta situação, serão necessários quatro passos para o carregamento dos dados. É possível notar que, devido às dimensões do bloco e da janela, uma parte das threads não carregará nenhum dado, representado pela parte vermelha mais clara.

/media/Attachments/courseIA366F2S2010/wd101245_artigo/fig3.png

Figura 3: Carregamento dos dados. Uma janela do tamanho do bloco movimenta-se sobre o conjunto de dados

Uma vez carregados os dados, tem início a segunda fase da execução, durante a qual o valor do pixel correspondente na imagem de saída será calculado através do método convencional. A vantagem em se utilizar a memória compartilhada está no fato de é feito uma operação de leitura e uma operação de escrita na memória global apenas. Uma vez carregados estes dados, o cálculo irá acessar a memória compartilhada, tirando proveito de sua proximidade e velocidade. A imagem de entrada foi armazenada na memória de textura e a máscara de convolução foi armazenada na memória de constantes.

Segue o pseudocódigo referente a esta implementação, utilizando a seguinte nomenclatura:

  • F : imagem de entrada;
  • G : imagem de saída;
  • H : máscara de filtragem;
  • lH : lado da máscara (2 * raio + 1);
  • sData : estrutura de dados na memória compartilhada;
  • raio : raio da máscara;
  • bloco : tamanho do bloco de execução;
  • w e h : largura e altura das imagens de entrada e de saída;
  • idx e idy : índices associados às imagens de entrada e saída;
  • tx e ty : índices da thread no bloco;
  • step : número de passos necessários para carregar os dados;

Algoritmo 2: Implementação em Memória Global

4.COMPARAÇÃO DE RESULTADOS

As tabelas abaixo mostram os resultados obtidos com a implementação da segunda abordagem (memória compartilhada, textura e constante), para tamanhos diferentes de máscara e imagens. A função de referência utilizada para comparação foi a iapconv, implementada em python e disponível neste ambiente como parte da biblioteca ia636.

Teste com máscara 3x3

Lado GPU CPU Speedup Compare
256 0.560 7.027 12.552 False
366 0.865 13.798 15.946 False
512 1.479 27.024 18.276 False
1055 5.430 128.944 23.745 False
2073 19.653 581.778 29.602 False

Teste com máscara 5x5

Lado GPU CPU Speedup Compare
256 0.593 19.127 32.240 False
366 0.979 37.775 38.584 False
512 1.684 73.699 43.755 False
1055 6.296 350.447 55.665 False
2073 21.347 1552.823 72.742 False

Teste com máscara 7x7

Lado GPU CPU Speedup Compare
256 0.676 37.301 55.186 False
366 1.141 73.786 64.686 False
512 1.992 143.774 72.176 False
1055 7.604 686.296 90.251 False
2073 26.320 3007.575 114.271 False

Teste com máscara 9x9

Lado GPU CPU Speedup Compare
256 0.783 61.570 78.589 False
366 1.350 121.712 90.172 False
512 2.402 237.514 98.882 False
1055 9.293 1127.953 121.378 False
2073 32.919 4957.392 150.595 False

A última coluna indica a similaridade das respostas obtidas com o algoritmo desenvolvido e com a função de referência. O resultado não foi idêntico, devido à diferença de tratamento das bordas. O algoritmo implementado neste trabalho usa a memória de textura no modo clamping, enquanto a função de referência faz o tratamento da borda considerando a imagem um sinal periódico. Desta forma, os valores dos pixels nas bordas são ligeiramente diferentes, mas numa escala muito pequena, conforme pode ser conferido nas imagens de saída.

Pode-se observar que a variação do tempo de execução ao aumentar o tamanho da máscara é muito maior na versão serial. Com isso, quanto maior a máscara, maior a vantagem da GPU em relação à CPU.

5.CONCLUSÕES

Neste trabalho, pôde-se observar claramente uma grande diferença nos tempos de execução da operação de convolução bidimensional com vantagem para a versão que utiliza CUDA. A linguagem permite que a transição entre os algoritmo serial e paralelo seja feito sem grande dificuldade, com desempenho bastante encorajador.

O código poderia ser otimizado ainda mais, observando características específicas do hardware, porém esta otimização implicaria numa degradação da legibilidade do código. A manutenção do código também seria muito difícil, já que a implementação ficaria muito dependente da arquitetura da placa. A substituição por uma placa de arquitetura mais nova implicaria na reescrita do código, para se adaptar às novas características. Deste modo, é desejável gerar código mais independente da arquitetura. A implementação sugerida conseguiu um desempenho considerável, mantendo o código bastante genérico, podendo ser executado com pouca ou nenhuma alteração nas placas de geração mais atual.

Uma otimização a ser implementada é a mudança do modo de leitura da memória de textura para o modo normalizado, o que permitiria considerar a imagem como periódica, obtendo um resultado idêntico à função de referência.

6.ANEXO

Código C/CUDA da solução para convolução bidimensional com máscara genérica

__global__ void cudaConvolution2DKernel ( float *dG,
                                int wi,
                                int hi,
                                int radius   )
{

    const int idx = IMUL(blockIdx.x, blockDim.x) + threadIdx.x;
    const int idy = IMUL(blockIdx.y, blockDim.y) + threadIdx.y;

    if(idx > wi || idy > hi) return;

    __shared__ unsigned char sData[2 * MAX_RADIUS + BLOCK_SIZE_X][2 * MAX_RADIUS + BLOCK_SIZE_Y ];

    const int step = iDivUp((2 * radius) + BLOCK_SIZE_X, BLOCK_SIZE_X),
              kLength = IMUL(2, radius) + 1;
    float pixel = 0.0f;

    #pragma unroll
    for(int i= 0; i < step; i++)
        #pragma unroll
       for (int j = 0; j < step; ++j)
       {
           sData[threadIdx.x + IMUL(i, BLOCK_SIZE_X)][threadIdx.y + IMUL(j, BLOCK_SIZE_Y)] = tex2D(tF, idx + IMUL(i, BLOCK_SIZE_X), idy + IMUL(j, BLOCK_SIZE_X));
       }

    __syncthreads();

    for(int i = 0 ; i < kLength; i++) // loop rows
    {
       #pragma unroll
       for(int j = 0 ; j < kLength; j++) // loop cols
       {
           pixel +=  cH[IMUL(kLength, (kLength-j))-(1+i)] * sData[threadIdx.x + i][threadIdx.y + j];
       }
    }

    dG[IMUL(idy, wi) + idx] =  pixel;
}