Projeto e Implementação de Convolução usando CUDA

Data:22 de outubro de 2010
Autor:André Körbes
Apresentação: Apresentação

Introdução

A operação de convolução tem grande importância no processamento de imagens, provendo a possibilidade de se aplicar técnicas de filtragem, oriundas do processamento de sinais unidimensionais e da filtragem no tempo, para se obter elementos filtrados no espaço, como bordas a partir das altas frequências, ou reduzir ruídos, como filtros gaussianos. Além destas características, o processamento em janelas - ou vizinhança - é utilizado em diversos outros algoritmos, especialmente na morfologia matemática. Desta forma, o domínio de técnicas de implementação de convolução, buscando a otimização, tem efeito também na implementação de outros filtros não lineares, como medianas, e em processamentos mais complexos, como a reconstrução morfológica.

Em termos matemáticos, a convolução 2-D é definida de forma contínua como a integração em e . Ao aplicar esta teoria sobre imagens, vale sua forma discreta, sob a forma de somatórios de multiplicações dos pontos da imagem. Entretanto, sua definição discreta, com parâmetros e variando entre e , assume a periodicidade, tanto da máscara quanto da imagem . Em termos práticos, e serão limitados pelas dimensões da imagem e da máscara.

Observando-se a equação da convolução, percebe-se que as operações são efetuadas para cada pixel de forma independente, ou seja, o resultado final de um não depende do resultado de outros, e dentro de uma janela. Estas duas características sugerem forte paralelismo no cálculo, com possibilidade de se aproveitar as arquiteturas de memória diferenciadas de modo a otimizar o acesso dentro da janela.

Esta mesma equação poderia ser reescrita de forma diferente, para refletir uma outra possibilidade de interpretação da convolução. Visto que a convolução é comutativa, poderia-se trocar a imagem pela máscara, e obter o resultado a partir de translações da imagem, ponderadas pelos valores da máscara, e somadas. Esta possibilidade não foi explorada e não será tratada neste trabalho.

Considerando-se a transformada de Fourier, onde a convolução é dada por uma multiplicação de sinais, de forma periódica, é indiferente onde considera-se o centro da máscara, pois o resultado final será uma translação periódica de qualquer posição possível do centro. Desta forma pode-se optar por implementações mais simplificadas, que evitem índices negativos, facilitando o cálculo da periodicidade.

Propostas de Implementação

Neste trabalho são avaliadas três propostas de implementação da convolução: (1) máscara genérica em memória global, (2) máscara fixa em memória compartilhada, (3) máscara genérica em blocos na memória compartilhada. Por máscara genérica entende-se máscaras de filtragem de tamanho arbitrário, ao contrário de máscara fixa, onde é pré-determinado um tamanho para estas, por exemplo 3x3. O tipo de memória utilizada na implementação é de extrema importância, pois afeta o desempenho do programa e sua forma de programação. A proposta por máscara genérica em memória global não foi implementada, sendo apenas ilustrativa para compreensão melhor do problema.

Máscara Genérica em Memória Global

Esta abordagem oferece a maior facilidade de implementação. A transição de um código sequencial para esta forma seria praticamente direta. Em relação à organização das threads CUDA, cada pixel seria calculado por uma thread, varrendo a máscara em sentido contrário. O seguinte pseudocódigo apresenta o algoritmo executado em cada kernel, onde representa um índice rasterizado, e os subscritos e indicam a imagem de entrada e a máscara de filtragem, respectivamente.

Algoritmo 1: Máscara Genérica em Memória Global

Nota-se que este pseudocódigo é voltado inteiramente para o algoritmo de cálculo da convolução, executando a rotação de sobre a máscara, e calculando em o somatório ponderado e atribuindo este ao resultado final. Entretanto, a simplicidade de código e no uso do hardware tem como característica a dependência do próprio hardware no que diz respeito ao desempenho. Isto ocorre pois a melhora das arquiteturas reduziu consideravelmente o tempo de acesso à memória global através do uso de caches. Entretanto, é de se esperar que o uso de máscaras de tamanho grande vejam piora no desempenho.

Reforça-se então que a grande, e talvez única, desvantagem desta implementação, seja a limitação do desempenho à capacidade do hardware de gerenciar memória e cache. Todavia, esta é a mesma limitação imposta à programação em CPUs sequenciais modernas, pois não se tem acesso à arquitetura de memória, porém, dada a capacidade de gerenciamento e ao tamanho do cache, isto não se mostra como um problema. Assim, de certa forma, apesar de ser a implementação mais simples, é plausível imaginar que em novos dispositivos seja ainda mais competitiva com implementações que utilizam diferentes formas de acesso à memória e otimizações de cálculo.

Máscara Fixa em Memória Compartilhada

De forma a melhorar o desempenho da implementação da convolução em memória global, um grande fator de modificação é utilizar a memória compartilhada. Esta memória tem maior banda de acesso, e é compartilhada entre as threads executando em um mesmo bloco. Entretanto, apesar dos diversos benefícios, é extremamente limitada em tamanho, 16Mb por processador de fluxo (com CUDA compute capability 1.3), além de ser fator limitante para a alocação de blocos ao processador [CUDA]. Assim, por ser um recurso limitado, deve ser utilizada de modo a não prejudicar o desempenho geral.

Uma primeira abordagem, útil ao se computar convoluções onde a máscara pode ser decomposta em máscaras menores, é fixar um tamanho padrão, 3x3 neste caso, e utilizar a memória compartilhada associada às threads do bloco, para realizar o carregamento de um bloco de tamanho conhecido. Assim, temos blocos de 16x16 threads, que efetuarão o carregamento da memória global para a memória de textura. Considerando-se os dados disponíveis, pode-se obter o número de pixels que serão efetivamente processados pelo bloco de threads. Neste caso, como é necessária uma borda com 2 pixels, o resultado serão blocos reduzidos, de cada 14x14 pixels. A Fig. figMascaraFixa apresenta esta ideia, onde os pixels com cores mais fortes são o resultado do bloco, e o contorno - da mesma cor, porém mais transparente - indica os pixels que são carregados na memória compartilhada.

/media/Attachments/courseIA366F2S2010/ak079907p1/mascara-fixa.png

Figure 1: Sobreposição de blocos para carregamento de bordas com bloco 6x6 e máscara 3x3

Assim, pode-se especificar um algoritmo, que, dada uma máscara de tamanho 3x3, utilize as threads para carregar para a memória compartilhada, e, após isto, apenas as threads associadas aos pixels válidos continuam a processar, e calculam efetivamente a convolução. Portanto, a quantidade de blocos necessários depende não mais do número de pixels apenas, mas do número de pixels em relação à quantidade de pixels que cada bloco será capaz de processar. Claramente, são necessárias mais threads do que o número de pixels, havendo assim uma sobrecarga devido à criação e gerenciamento destas, e ao maior número de blocos necessários. Entretanto, o uso da memória compartilhada compensa esta sobrecarga por ser de acesso muito mais rápido, e de fato compartilhando os valores necessários dentro do bloco. Além deste fato, a arquitetura CUDA e intrínseca do hardware gráfico utiliza threads leves, cujo gerenciamento traz sobrecarga quase desprezível.

Algoritmo 2: Máscara Fixa em Memória Compartilhada

Há diversos aspectos neste algoritmo a serem considerados. A introdução da memória compartilhada, representada por , se dá através de uma matriz, indexada pelos índices da thread, e . Pode-se interpretar esta proposta como em dois passos distintos, sendo o primeiro o carregamento, sincronização e término antecipado das threads que não estão associadas a pixels dentro do bloco reduzido; e um segundo passo onde ocorre o cálculo da convolução utilizando os dados da memória compartilhada. A implementação desta proposta é dada pelo seguinte código C/CUDA.

__global__ void conv2D_kernel(unsigned char *a, int height, int width, float *mask, float * out)
{
    __shared__ unsigned char data[TILE_WIDTH][TILE_WIDTH];

    int x = REAL_TILE*blockIdx.x + threadIdx.x - 1;
    int y = REAL_TILE*blockIdx.y + threadIdx.y - 1;
    int idx = y*width + x;

    if (x < 0 || y < 0 || x >= width || y >= height)
    {
        data[threadIdx.y][threadIdx.x] = 0; // zeroed border
    }
    else
    {
        data[threadIdx.y][threadIdx.x] = a[idx];
    }

    __syncthreads();

    if (x < 0 || y < 0 || x >= width || y >= height ||
        threadIdx.x == 0 || threadIdx.x == TILE_WIDTH - 1 ||
        threadIdx.y == 0 || threadIdx.y == TILE_WIDTH - 1)
        return;

    float value = 0;

    for (int dx = -1, mx = 1; dx < 2; dx++, mx--)
    {
        for (int dy = -1, my = 1; dy < 2; dy++, my--)
        {
            value += data[threadIdx.y + dy][threadIdx.x + dx] * mask[(my+1)*3 + (mx+1)];
        }
    }

    out[idx] = value;


}

Uma diferença entre o código apresentado e a proposta apresentada inicialmente diz respeito aos valores de borda da imagem. Neste código, a imagem não é considerada como periódica, e utiliza-se a borda com valor fixo em zero. Além disto, continua-se utilizando o centro da máscara em sua coordenada (0,0), o que, para efeito de comparação com outras implementações, também requer que estas utilizem borda zerada, e uma translação periódica da imagem resultante.

Porém esta implementação não reflete o propósito final deste artigo, que é a implementação da convolução utilizando máscaras genéricas. Para este propósito, são possíveis diversas alternativas, onde um ou outro recurso limitante será desperdiçado, sejam threads utilizadas apenas para o carregamento, espaço em memória compartilhada por se utilizar menos espaço do que o alocado, ou a alocação máxima de blocos por processador, por estar limitado devido ao tamanho da memória compartilhada necessária. A próxima proposta apresentada busca minimizar estes dois efeitos, utilizando todas as threads do bloco para processamento de pixels na imagem final, e um espaço em memória compartilhada delimitado pelo bloco de threads.

Máscara Genérica em Blocos na Memória Compartilhada

Ao tentar encontrar uma alternativa viável para o processamento de máscaras extremamente grandes, uma opção busca inspiração no método de multiplicação de matrizes, apresentado no livro de Kirk e Hwu [KirkHwu]. Nesta proposta, todas as threads do bloco são utilizadas para carregar um bloco de dados para a memória compartilhada, realizar o processamento parcial sobre estes dados, e movimentar o bloco sobre a imagem de entrada, carregando e processando em passos, de forma a cobrir toda a área necessária para o bloco de pixels de saída, do mesmo tamanho do bloco de threads utilizado. A Fig. figMascaraGenerica apresenta o esquema de carregamento dos blocos, onde cada cor representa os pixels carregados na memória em um determinado passo do processamento. O bloco resultante deste processamento é o bloco de 6x6 pixels vermelho. Os pixels em vermelho que são cobertos pelas áreas verde, azul e amarela são representações da borda necessária para processar uma máscara 3x3. Entretanto, esta máscara poderia ser estendida até 7x7, sem ser necessário carregar mais dados.

/media/Attachments/courseIA366F2S2010/ak079907p1/mascara-generica.png

Figure 2: Carregamento de 4 blocos 6x6 para processamento de um bloco de pixels com máscara até 7x7

Esta proposta, apesar de buscar otimizar o uso de memória e das threads, tem como desvantagem a complexidade adicional na implementação. Esta complexidade é refletida no cálculo dos índices da máscara e da imagem, em cada passo. Entretanto, esta complexidade permite grande flexibilidade em relação ao tamanho das máscaras, fato desejável quando se expande o horizonte de uso desta técnica para algoritmos morfológicos, onde o elemento estruturante, que determina a vizinhança, nem sempre é quadrado. De forma abstrata, o pseudocódigo é apresentado a seguir.

Algoritmo 3: Máscara Genérica em Blocos na Memória Compartilhada

Claramente, são introduzidos dois novos laços, para controle de qual bloco está sendo visitado. Entretanto, os cálculos dos índices tanto da imagem quanto da máscara dependem da posição inicial na imagem e do bloco atual. Isto implica em uma implementação mais intrincada, e sujeita a erros. Além disto, tais erros tornam-se difíceis de serem rastreados, pelo fato do algoritmo executar em vários passos iterativos. Desta forma, o código-fonte abaixo apresenta o programa C/CUDA que implementa esta proposta.

//
// CUDA Kernel
//
__global__ void conv2D_kernel(unsigned char *g_f, int dimyf, int dimxf, int dimym, int dimxm, float *g_g)
{
    __shared__ unsigned char s_f[TILE_WIDTH][TILE_WIDTH];

    int posxf = TILE_WIDTH*blockIdx.x + threadIdx.x;
    int posyf = TILE_WIDTH*blockIdx.y + threadIdx.y;
    int ixf;
    int posxsf = threadIdx.x;
    int posysf = threadIdx.y;

    int dimxTotal = dimxm + TILE_WIDTH;
    int dimyTotal = dimym + TILE_WIDTH;

    int dimxTiles = (int)(dimxTotal/TILE_WIDTH) + (dimxTotal % TILE_WIDTH == 0 ? 0 : 1);
    int dimyTiles = (int)(dimyTotal/TILE_WIDTH) + (dimyTotal % TILE_WIDTH == 0 ? 0 : 1);

    int os, ot;

    float v = 0;
    int oxf = 0, oxm = dimxm-1, oyf = 0, oym = dimym-1;
    int oxmBegin = dimxm-1;
    int oxfBegin = 0;

    int oymBegin = dimym-1;
    int oyfBegin = 0;

    for (int t = 0; t < dimyTiles; t++)
    {
        // calcula offset em T
        ot = t * TILE_WIDTH;
        for (int s = 0; s < dimxTiles; s++)
        {
            // calcula offset em S
            os = s * TILE_WIDTH;

            // calcula indice na imagem (periodico)
            ixf = ((posyf + ot) % dimyf)*dimxf + ((posxf + os) % dimxf);

            // carrega para smem (provavel problema de coalescencia aqui, a borda periodica nao deveria ser carregada assim)
            // poderia colocar uma condicao de soh carregar se posxf + os < TILE_WIDTH + dimxm, mas pode ser que seja pior
            s_f[posysf][posxsf] = g_f[ixf];

            // aguarda carregamento
            __syncthreads();

            oyf = oyfBegin;
            oym = oymBegin;

            // calcula o valor parcial do pixel
            for (int cy = 0; cy < TILE_WIDTH; cy++)
            {
                if (oyf < dimym && (posysf + oyf) - ot < TILE_WIDTH)
                {
                    oxf = oxfBegin;
                    oxm = oxmBegin;
                }

                for (int cx = 0; cx < TILE_WIDTH; cx++)
                {
                    if (oxf < dimxm && (posxsf + oxf) - os < TILE_WIDTH)
                    {
                        v += s_f[(posysf + oyf) - ot][(posxsf + oxf) - os] * c_m[oym*dimym + oxm];
                        oxf++;
                        oxm--;
                    }
                }

                if (oyf < dimym && (posysf + oyf) - ot < TILE_WIDTH)
                {
                    oyf++;
                    oym--;
                }
            }

            oxfBegin = oxf;
            oxmBegin = oxm;

        }

        oyfBegin = oyf;
        oymBegin = oym;

        oxfBegin = 0;
        oxmBegin = dimxm-1;

    }

    if (posxf >= dimxf || posyf >= dimyf)
        return;

    ixf = posyf*dimxf + posxf;

    g_g[ixf] = v;

}

Pode-se ver claramente a diferença entre as propostas até mesmo pelo número de linhas necessário para cada implementação. Nesta implementação, o cálculo dos índices é realizado de forma incremental, armazenando-se os últimos índices válidos utilizados no bloco anterior. Desta forma, cada bloco é varrido completamente e utilizam-se apenas os pixels válidos. Esta implementação, apesar de buscar otimizar o acesso à memória e uso do processador, executa mais instruções, e dependendo do caso possui warps divergentes.

Comparação de Desempenho

A comparação de desempenho executa a segunda proposta de algoritmo, para matrizes aleatórias de diversos tamanhos, com máscaras de lado 3, 5, 7 e 9. As máscaras são calculadas como a média dos pixels, mas seu conteúdo é irrelevante para o teste, bem como o das imagens. O algoritmo tem seu resultado conferido com a implementação iapconv da toolbox ia636.

 1 from funcTimer import *
 2 import scipy.signal as sp
 3 import ak079907_8_2 as m2
 4 import ak079907_4_cod as m1
 5 
 6 import numpy.random as rp
 7 import time
 8 from ia636 import iapconv, iaptrans
 9 
10 mmfreedom(2)
11 
12 f = adreadgray('cameraman.tif')
13 
14 dims = [177, 2027, 366, 655,1766,1076]#,2021,1305, 400, 493,1964, 939, 946, 139,1499, 659, 431,1046, 298,1110,1577, 117,1500,1747,2006,1629]
15 dims.sort()
16 
17 def Compare(a,b):
18     c = abs(a - b)
19     c = c > 0.001
20     return not c.any()
21 
22 
23 def AvgRatio(d, h):
24     f = rp.randint(1000, size=(d,d)).astype('uint8')
25 
26     tGPU2,rGPU2 = funcTimer( lambda: m2.conv2D( f, h ) )
27 
28     tCPU2,rCPU2 = funcTimer( lambda: iapconv( f, h) )
29 
30     c2 = Compare(iaptrans(rGPU2, list(array(h.shape)/2)), rCPU2)
31 
32     return d, tGPU2, tCPU2, tCPU2 / tGPU2, c2
33 
34 def MakeTable(hsize):
35     h = (ones((hsize,hsize))/(hsize**2)).astype('float32')
36     print 'Mascara: %dx%d' % (hsize,hsize)
37     print
38     print '==== ========== ========== ========== ========'
39     print 'Lado GPU        CPU        Speedup    Compare '
40     print '==== ========== ========== ========== ========'
41 
42     for d in dims:
43         print '%4d %10.3f %10.3f %10.3f %s' % AvgRatio(d,h)
44 
45     print '==== ========== ========== ========== ========'
46     print
47 
48 MakeTable(3)
49 MakeTable(5)
50 MakeTable(7)
51 MakeTable(9)

Mascara: 3x3

Lado GPU CPU Speedup Compare
177 1.604 3.489 2.175 True
366 4.942 13.791 2.790 True
655 14.837 45.784 3.086 True
1076 39.973 132.786 3.322 True
1766 104.935 367.498 3.502 True
2027 137.100 491.266 3.583 True

Mascara: 5x5

Lado GPU CPU Speedup Compare
177 1.665 9.508 5.711 True
366 5.193 37.810 7.281 True
655 15.591 124.021 7.955 True
1076 41.642 368.184 8.842 True
1766 109.803 992.262 9.037 True
2027 143.980 1309.769 9.097 True

Mascara: 7x7

Lado GPU CPU Speedup Compare
177 1.758 18.558 10.556 True
366 5.604 73.759 13.162 True
655 16.551 241.138 14.569 True
1076 44.352 713.195 16.080 True
1766 116.889 1926.918 16.485 True
2027 152.550 2536.331 16.626 True

Mascara: 9x9

Lado GPU CPU Speedup Compare
177 1.891 30.582 16.173 True
366 6.064 121.708 20.070 True
655 17.977 398.479 22.166 True
1076 48.207 1168.826 24.246 True
1766 126.575 3166.020 25.013 True
2027 166.074 4179.244 25.165 True

Como se pode ver nas tabelas, o algoritmo tem desempenho melhor quanto maior a máscara utilizada. Isto vem do fato que o paralelismo começa a ter efeitos maiores, e que os dados adicionais carregados pelos blocos também passam a ser utilizados no cômputo do bloco de pixels. O tempo da operação em CPU aumenta pois é necessário um maior número de operações sobre a matriz, conforme o tamanho da máscara. Não é possível comparar de forma sistemática a primeira proposta, pois aquela é fixada em máscaras de 3x3. Entretanto, em uma comparação apenas com este tamanho, aquela proposta é significativamente melhor. Todavia, não é possível afirmar se esse comportamento é escalável para máscaras maiores.

Conclusões

Neste trabalho foram apresentadas três propostas de implementação da convolução utilizando CUDA e explorando sua arquitetura de memória. Estas propostas foram baseadas no entendimento mais comum da convolução, como uma soma ponderada sobre um conjunto de dados. As propostas partem da abordagem trivial, onde não se utilizam os recursos de memória compartilhada, para uma proposta onde são utilizadas mais threads do que o número de pixels na imagem, e por último uma proposta baseada no carregamento de blocos, e cálculos parciais do resultado final.

Cada proposta traz vantagens e desvantagens. Este fato, combinado com a constante evolução do hardware, aponta para o uso desta arquitetura buscando algoritmos mais simples, mas que aproveitem a natureza paralela do problema. Pelos resultados obtidos, percebe-se uma tendência a utilizar a memória compartilhada de forma simplificada, utilizando threads adicionais quando necessário, devido à baixa sobrecarga de gerenciamento destas. Ainda neste sentido, pode-se também vislumbrar o uso das 2 abordagens, que utilizam esta memória, em conjunto, de forma complementar, extraindo o maior desempenho dependendo do caso, tendo como parâmetros o tamanho da máscara e da imagem.

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.