Um Estudo sobre Algoritmos de Convolução Digital em GPUs

Victor Oliveira

Introdução

Neste trabalho avaliaremos diversas formas de implementação de algoritmos de convolução digital em placas gráficas.

A convolução é uma operação entre duas funções e que para o caso de e serem sinais discretos bidimensionais pode ser descrita da forma:

Se considerarmos e como imagens, teremos a convolução de uma imagem por uma máscara .

Como computadores são máquinas finitas, precisamos remover a infinitude da equação acima, isso pode ser feito de duas maneiras:

  • Consideramos que e são sinais periódicos finitos de tamanho , desse modo em cada dimensão.
  • O mesmo que 1., mas e são considerados como zero fora dos limites de sua dimensão.

Vamos a um exemplo de aplicação:

Filtragem Gaussiana

A convolução de uma imagem por uma máscara , em que os pesos de são a discretização dos valores de uma distribuição gaussiana tem como resultado o "espalhamento" dos pontos de , dando uma impressão de suavização na imagem.

 1 import ia636
 2 import numpy as np
 3 
 4 image = adread('cameraman.pgm').astype('float32')
 5 
 6 mask = 1.0/16 * np.array([
 7                 [1, 2, 1],
 8                 [2, 4, 2],
 9                 [1, 2, 1]
10                 ], dtype=np.float32)
11 
12 conv = ia636.iapconv(image, mask)
13 
14 adshow(image)
15 adshow(conv)

Convolução em GPUs

GPUs (placas gráficas) são hardware extremamente otimizado para tarefas paralelas e em que o acesso a memória é localizado.

Assim, vemos que a convolução é apropriada para implementação em GPU, pois o resultado de cada pixel em é independente, tornando a operação bastante paralelizável. Além disso, temos que há a avaliação dos pixels em torno de uma vizinhança para cada . O que caracteriza um acesso localizado à memória.

Neste trabalho usaremos a definição de convolução em que o sinal é zero fora dos limites e a máscara é centrada em . Essa definição pode ser facilmente alterada para tratar os outros casos.

Implementação

1 - Memória global

Começaremos com uma implementação básica, que usa acesso direto à memória global. Ela será a base para comparação com as outras implementações.

No modelo CUDA de programação, cada thread na GPU recebe um índice. Esse índice pode ser linear, bidimensional ou tridimensional, isto é, as threads podem ser agrupadas como sequências de linhas, retângulos ou paralelepípedos.

/media/Attachments/courseIA366F2S2010/vm072589_p1/relation.png

Figura 1: Vemos aqui um agrupamento 2D das threads, isso cria uma relação natural com uma imagem, pois o índice x e y da thread pode ser considerado como o pixel (x,y) na imagem.

 1 #include <string.h>
 2 #include <math.h>
 3 #include <cuda.h>
 4 
 5 #include "simple_arrays.h"
 6 #define BLOCK_X 16
 7 #define BLOCK_Y 16
 8 
 9 __global__ void convolution_simple_kernel(float* dImg, float* dOut, int imgW, int imgH, float dMask[], int maskW, int maskH) {
10     int x = blockIdx.x * blockDim.x + threadIdx.x;
11     int y = blockIdx.y * blockDim.y + threadIdx.y;
12     float v = 0.0f;
13 
14     if (x < imgW && y < imgH) {
15         for (int dy=0; dy<maskH; dy++) {
16             for (int dx=0; dx<maskW; dx++) {
17                 int px = x+dx;
18                 int py = y+dy;
19                 if (px>=0 && px<imgW && py>=0 && py<imgH)
20                     v += dMask[(maskH-dy-1)*maskW+(maskW-dx-1)] * dImg[py*imgW+px];
21             }
22         }
23         dOut[y*imgW+x] = v;
24     }
25 
26 }
27 
28 void stub_convolution_simple(float* gImg, float* gOut, int imgW, int imgH, float *gMask, int maskW, int maskH) {
29     float* dImg;
30     float* dOut;
31     float* dMask;
32 
33     cudaMalloc((void **)&dImg,  imgW*imgH*sizeof(float));
34     cudaMalloc((void **)&dOut,  imgW*imgH*sizeof(float));
35     cudaMalloc((void **)&dMask, maskW*maskH*sizeof(float));
36 
37     cudaMemcpy(dImg, gImg, imgW*imgH*sizeof(float), cudaMemcpyHostToDevice);
38     cudaMemcpy(dMask, gMask, maskW*maskH*sizeof(float), cudaMemcpyHostToDevice);
39 
40     dim3 gridDim((imgW+BLOCK_X-1)/BLOCK_X, (imgH+BLOCK_Y-1)/BLOCK_Y);
41     dim3 blockDim(BLOCK_X, BLOCK_Y);
42 
43     convolution_simple_kernel <<< gridDim, blockDim >>>
44         (dImg, dOut, imgW, imgH, dMask, maskW, maskH);
45 
46     cudaMemcpy(gOut, dOut, imgW*imgH*sizeof(float), cudaMemcpyDeviceToHost);
47 
48     cudaFree(dImg);
49     cudaFree(dOut);
50     cudaFree(dMask);
51 }
1 void convolution_simple(int h, int w, float in[], int my, int mx, float *mask,
2                         float **out, int *out_y, int *out_x)
3 {
4     *out = new float[w*h];
5     *out_y = h;
6     *out_x = w;
7 
8     stub_convolution_simple(in, *out, w, h, mask, mx, my);
9 }

Como podemos ver no código, cada thread verifica seu pixel correspondente na imagem (em oposição à máscara) e seus vizinhos para fazer a operação de soma dos valores ponderados pelo valor na máscara. Podemos ver que basicamente fizemos uma paralelização dos dois somatórios presentes na equação da convolução.

2 - Memória global com máscara na memória constante

Nessa implementação, apenas colocamos a máscara na memória constante, que é uma memória pequena, mas visível a todas as threads e de acesso muito rápido.

Esse uso é feito através da declaração __constant__, temos uma limitação aqui que colocamos o tamanho máximo do kernel 9x9, isso pode ser alterado no código com um #DEFINE evidentemente. Mas a posição que seguiremos nesse artigo é que para máscaras grandes não vale a pena usar a convolução, mas sim métodos , como a Fast Fourier Transform.

 1 #define MASK_MAX 20
 2 
 3 __constant__ float cMask[MASK_MAX*MASK_MAX];
 4 
 5 __global__ void convolution_const_kernel(float* dImg, float* dOut, int imgW, int imgH, int maskW, int maskH) {
 6     int x = blockIdx.x * blockDim.x + threadIdx.x;
 7     int y = blockIdx.y * blockDim.y + threadIdx.y;
 8     float v = 0.0f;
 9 
10     if (x < imgW && y < imgH) {
11         for (int dy=0; dy<maskH; dy++) {
12             for (int dx=0; dx<maskW; dx++) {
13                 int px = x+dx;
14                 int py = y+dy;
15                 if (px>=0 && px<imgW && py>=0 && py<imgH)
16                     v += cMask[(maskH-1-dy)*maskW+(maskW-1-dx)] * dImg[py*imgW+px];
17             }
18         }
19         dOut[y*imgW+x] = v;
20     }
21 
22 }
23 
24 void stub_convolution_const(float* gImg, float* gOut, int imgW, int imgH, float *gMask, int maskW, int maskH) {
25     float* dImg;
26     float* dOut;
27 
28     cudaMalloc((void **)&dImg,  imgW*imgH*sizeof(float));
29     cudaMalloc((void **)&dOut,  imgW*imgH*sizeof(float));
30 
31     cudaMemcpy(dImg, gImg, imgW*imgH*sizeof(float), cudaMemcpyHostToDevice);
32     cudaMemcpyToSymbol(cMask, gMask, maskW*maskH*sizeof(float), 0, cudaMemcpyHostToDevice);
33 
34     dim3 gridDim((imgW+BLOCK_X-1)/BLOCK_X, (imgH+BLOCK_Y-1)/BLOCK_Y);
35     dim3 blockDim(BLOCK_X, BLOCK_Y);
36 
37     convolution_const_kernel <<< gridDim, blockDim >>>
38         (dImg, dOut, imgW, imgH, maskW, maskH);
39 
40     cudaMemcpy(gOut, dOut, imgW*imgH*sizeof(float), cudaMemcpyDeviceToHost);
41 
42     cudaFree(dImg);
43     cudaFree(dOut);
44 }

3 - Uso de Texturas

Nesta implementação, usaremos a memória de textura, que é uma forma de cache sobre a memória global otimizada para operações locais bidimensionais, como filtragem.

Basicamente, carregamos os dados da imagem em um cuArray e habilitamos o suporte a textura em tex. No kernel, só há a alteração do acesso a imagem, que agora é feito através da função tex2D, que pega o pixel (x,y) na textura tex.

/media/Attachments/courseIA366F2S2010/vm072589_p1/texture.png

Figura 2: Cada acesso à memória global agora passa por um cache que é otimizado para acesso 2D localizado. Portanto, quanto mais o uso da memória divergir desse modelo, menos o uso de texturas vai ajudar no desempenho.

 1 texture<float, 2, cudaReadModeElementType> tex;
 2 
 3 __global__ void convolution_texture_kernel(float* dOut, int imgW, int imgH, int maskW, int maskH) {
 4     int x = blockIdx.x * blockDim.x + threadIdx.x;
 5     int y = blockIdx.y * blockDim.y + threadIdx.y;
 6     float v = 0.0f;
 7 
 8     if (x < imgW && y < imgH) {
 9         for (int dy=0; dy<maskH; dy++) {
10             for (int dx=0; dx<maskW; dx++) {
11                 int px = x+dx;
12                 int py = y+dy;
13                 if (px>=0 && px<imgW && py>=0 && py<imgH)
14                     v += cMask[(maskH-1-dy)*maskW+(maskW-1-dx)] * tex2D(tex, px, py);
15             }
16         }
17         dOut[y*imgW+x] = v;
18     }
19 
20 }
21 
22 void stub_convolution_texture(float* gImg, float* gOut, int imgW, int imgH, float *gMask, int maskW, int maskH) {
23     float* dOut;
24 
25     cudaMalloc((void **)&dOut,  imgW*imgH*sizeof(float));
26     cudaMemcpyToSymbol(cMask, gMask, 9*sizeof(float), 0, cudaMemcpyHostToDevice);
27 
28     cudaChannelFormatDesc floatdesc =
29         cudaCreateChannelDesc<float>();
30 
31     cudaArray* cuArray;
32     cudaMallocArray(&cuArray, &floatdesc, imgW, imgH);
33     cudaMemcpyToArray(cuArray, 0, 0, gImg, imgW*imgH*sizeof(float), cudaMemcpyHostToDevice);
34 
35     cudaBindTextureToArray(tex, cuArray, floatdesc);
36 
37     dim3 gridDim((imgW+BLOCK_X-1)/BLOCK_X, (imgH+BLOCK_Y-1)/BLOCK_Y);
38     dim3 blockDim(BLOCK_X, BLOCK_Y);
39 
40     convolution_texture_kernel <<< gridDim, blockDim >>>
41         (dOut, imgW, imgH, maskW, maskH);
42 
43     cudaMemcpy(gOut, dOut, imgW*imgH*sizeof(float), cudaMemcpyDeviceToHost);
44 
45     cudaFreeArray(cuArray);
46     cudaFree(dOut);
47 }

4 - Memória compartilhada

Nesta impĺementação, faremos com que cada pixel carregue seus vizinhos na memória compartilhada, de modo que cada bloco carrega na imagem os pixels correspondentes aos seus mais uma borda. Esse carregamente é feito de forma ordenada (coalesced), assim temos o melhor padrão de acesso possível à memória global.

Convolução usando a memória compartilhada (baseado no código do André Körbes):

  1 #include <string.h>
  2 #include <math.h>
  3 #include <cuda.h>
  4 
  5 #define CHECK_ERROR(msg) if(checkError(msg)) return
  6 #define TILE_WIDTH   16
  7 #define MAX_MASK     2500
  8 
  9 //
 10 // Error Checking
 11 //
 12 int checkError(char *msg)
 13 {
 14     cudaError_t status = cudaGetLastError();
 15     if(status != cudaSuccess) {
 16         pyprintf("*** %s [%s] ***\n", msg, cudaGetErrorString(status));
 17         return 1;
 18     }
 19     return 0;
 20 }
 21 
 22 __constant__ float c_m[MAX_MASK];
 23 
 24 //
 25 // CUDA Kernel
 26 //
 27 __global__ void convolution_shared_kernel(float *g_f, int dimyf, int dimxf, int dimym, int dimxm, float *g_g)
 28 {
 29     __shared__ float s_f[TILE_WIDTH][TILE_WIDTH];
 30 
 31     int posxf = TILE_WIDTH*blockIdx.x + threadIdx.x;
 32     int posyf = TILE_WIDTH*blockIdx.y + threadIdx.y;
 33     int ixf;
 34     int posxsf = threadIdx.x;
 35     int posysf = threadIdx.y;
 36 
 37     int dimxTotal = dimxm + TILE_WIDTH;
 38     int dimyTotal = dimym + TILE_WIDTH;
 39 
 40     int dimxTiles = (int)(dimxTotal/TILE_WIDTH) + (dimxTotal % TILE_WIDTH == 0 ? 0 : 1);
 41     int dimyTiles = (int)(dimyTotal/TILE_WIDTH) + (dimyTotal % TILE_WIDTH == 0 ? 0 : 1);
 42 
 43     int os, ot;
 44 
 45     float v = 0;
 46     int oxf = 0, oxm = dimxm-1, oyf = 0, oym = dimym-1;
 47     int oxmBegin = dimxm-1;
 48     int oxfBegin = 0;
 49 
 50     int oymBegin = dimym-1;
 51     int oyfBegin = 0;
 52 
 53     for (int t = 0; t < dimyTiles; t++)
 54     {
 55         // calcula offset em T
 56         ot = t * TILE_WIDTH;
 57         for (int s = 0; s < dimxTiles; s++)
 58         {
 59             // calcula offset em S
 60             os = s * TILE_WIDTH;
 61 
 62             // calcula indice na imagem (periodico)
 63             ixf = ((posyf + ot) % dimyf)*dimxf + ((posxf + os) % dimxf);
 64 
 65             // carrega para smem (provavel problema de coalescencia aqui, a borda periodica nao deveria ser carregada assim)
 66             // poderia colocar uma condicao de soh carregar se posxf + os < TILE_WIDTH + dimxm, mas pode ser que seja pior
 67             s_f[posysf][posxsf] = g_f[ixf];
 68 
 69             // aguarda carregamento
 70             __syncthreads();
 71 
 72             oyf = oyfBegin;
 73             oym = oymBegin;
 74 
 75             // calcula o valor parcial do pixel
 76             for (int cy = 0; cy < TILE_WIDTH; cy++)
 77             {
 78                 if (oyf < dimym && (posysf + oyf) - ot < TILE_WIDTH)
 79                 {
 80                     oxf = oxfBegin;
 81                     oxm = oxmBegin;
 82                 }
 83 
 84                 for (int cx = 0; cx < TILE_WIDTH; cx++)
 85                 {
 86                     if (oxf < dimxm && (posxsf + oxf) - os < TILE_WIDTH)
 87                     {
 88                         v += s_f[(posysf + oyf) - ot][(posxsf + oxf) - os] * c_m[oym*dimym + oxm];
 89                         oxf++;
 90                         oxm--;
 91                     }
 92                 }
 93 
 94                 if (oyf < dimym && (posysf + oyf) - ot < TILE_WIDTH)
 95                 {
 96                     oyf++;
 97                     oym--;
 98                 }
 99             }
100 
101             oxfBegin = oxf;
102             oxmBegin = oxm;
103 
104         }
105 
106         oyfBegin = oyf;
107         oymBegin = oym;
108 
109         oxfBegin = 0;
110         oxmBegin = dimxm-1;
111 
112     }
113 
114     if (posxf >= dimxf || posyf >= dimyf)
115         return;
116 
117     ixf = posyf*dimxf + posxf;
118 
119     g_g[ixf] = v;
120 
121 }
122 
123 //
124 // Kernel Stub
125 //
126 void stub_convolution_shared(float *h_f, int dimyf, int dimxf, float *h_m, int dimym, int dimxm, float *h_g)
127 {
128     float *d_f;
129     float *d_g;
130 
131     int size = dimxf*dimyf;
132 
133     if (dimym * dimxm >= MAX_MASK)
134     {
135         pyprintf("Erro: mask is too big");
136         return;
137     }
138 
139 
140     // Allocate device memory
141     cudaMalloc((void **)&d_f, size);
142     CHECK_ERROR("cudaMalloc d_f");
143 
144     cudaMalloc((void **)&d_g, size*sizeof(float));
145     CHECK_ERROR("cudaMalloc d_g");
146 
147     cudaMemset(d_g, 0, size*sizeof(float));
148 
149     // Copy data to device
150     cudaMemcpy(d_f, h_f, size, cudaMemcpyHostToDevice);
151     CHECK_ERROR("cudaMemcpy h_f -> d_f");
152 
153     cudaMemcpyToSymbol(c_m, h_m, dimym*dimxm*sizeof(float), 0, cudaMemcpyHostToDevice);
154     CHECK_ERROR("cudaMemcpy h_m -> c_m");
155 
156     // Invoke kernel
157     int gridW = (int)(dimxf/TILE_WIDTH) + (dimxf % TILE_WIDTH == 0 ? 0 : 1);
158     int gridH = (int)(dimyf/TILE_WIDTH) + (dimyf % TILE_WIDTH == 0 ? 0 : 1);
159 
160     dim3 gridDim(gridW, gridH);
161     dim3 blockDim(TILE_WIDTH, TILE_WIDTH);
162 
163     convolution_shared_kernel<<<gridDim, blockDim>>>(d_f, dimyf, dimxf, dimym, dimxm, d_g);
164     CHECK_ERROR("kernel invoke");
165 
166     // Copy data from device
167     cudaMemcpy(h_g, d_g, size*sizeof(float), cudaMemcpyDeviceToHost);
168     CHECK_ERROR("cudaMemcpy d_g -> h_g");
169 
170     // Release allocated device memory
171     cudaFree(d_f);
172     CHECK_ERROR("cudaFree d_f");
173 
174     cudaFree(d_g);
175     CHECK_ERROR("cudaFree d_g");
176 }
177 
178 //
179 // Exported Function
180 //
181 Image32F *convolution_shared(Image32F *A, Image32F *mask)
182 {
183     Image32F *output = new Image32F(A->nd, A->dims);
184     stub_convolution_shared((float*)A->raster, A->dims[0], A->dims[1], (float*)mask->raster, mask->dims[0], mask->dims[1], (float*)output->raster);
185     return output;
186 }
OK

Resultados

 1 void convolution_simple(int DIM1, int DIM2, float *IN_ARRAY2,
 2                         int DIM1, int DIM2, float *IN_ARRAY2,
 3                         float **ARGOUT_ARRAY2, int *DIM1, int *DIM2);
 4 
 5 
 6 void convolution_const(int DIM1, int DIM2, float *IN_ARRAY2,
 7                        int DIM1, int DIM2, float *IN_ARRAY2,
 8                        float **ARGOUT_ARRAY2, int *DIM1, int *DIM2);
 9 
10 void convolution_texture(int DIM1, int DIM2, float *IN_ARRAY2,
11                          int DIM1, int DIM2, float *IN_ARRAY2,
12                          float **ARGOUT_ARRAY2, int *DIM1, int *DIM2);
13 
14 Image32F *convolution_shared(Image32F *A, Image32F *mask);
OK [C/C++ extension is up-to-date]
 1 def run():
 2     import vm072589_p1 as vm
 3 
 4     import numpy as np
 5     from ia636 import iaconv
 6     import time
 7 
 8     image = np.ones((4000, 4000)).astype('float32')
 9 
10     masks = [3, 5, 7, 9, 20]
11 
12     print 'Speedup considerando o tamanho da máscara e formas de implementação:'
13     print ''
14     print "%20s %20s %20s %20s %20s" % (('='*20,)*5)
15     print "%-20s %-20s %-20s %-20s %-20s" % ("..", " Simple", " Constant mem", " Texture mem", " Shared mem")
16     print "%20s %20s %20s %20s %20s" % (('='*20,)*5)
17 
18     for kk, mask_s in enumerate(masks):
19 
20         mask = np.ones((mask_s, mask_s), 'float32')/(mask_s*mask_s)
21 
22         v = []
23         for i in range(1):
24             t1 = 1000.0*time.time()
25             iaconv(image, mask)
26             t2 = 1000.0*time.time()
27             v.append(t2-t1)
28 
29         t_serial = min(v)
30 
31         v = []
32         for i in range(1):
33             t1 = 1000.0*time.time()
34             vm.convolution_simple(image, mask)
35             t2 = 1000.0*time.time()
36             v.append(t2-t1)
37 
38         t_simple = min(v)
39 
40         v = []
41         for i in range(1):
42             t1 = 1000.0*time.time()
43             vm.convolution_const(image, mask)
44             t2 = 1000.0*time.time()
45             v.append(t2-t1)
46 
47         t_const = min(v)
48 
49         v = []
50         for i in range(1):
51             t1 = 1000.0*time.time()
52             vm.convolution_texture(image, mask)
53             t2 = 1000.0*time.time()
54             v.append(t2-t1)
55 
56         t_texture = min(v)
57 
58         v = []
59         for i in range(1):
60             t1 = 1000.0*time.time()
61             vm.convolution_shared(image, mask)
62             t2 = 1000.0*time.time()
63             v.append(t2-t1)
64 
65         t_shared = min(v)
66 
67         print "%20s %20.2f %20.2f %20.2f %20.2f" % ("%dx%d" % (mask_s, mask_s), t_serial/t_simple, t_serial/t_const, t_serial/t_texture, t_serial/t_shared)
68         if mask_s != masks[-1]: print "%20s %20s %20s %20s %20s" % (('-'*20,)*5)
69 
70     print "%20s %20s %20s %20s %20s" % (('='*20,)*5)
71     print
72 
73 run()

Speedup considerando o tamanho da máscara e formas de implementação:

Simple Constant mem Texture mem Shared mem
3x3 1.10 11.47 11.86 1.75
5x5 18.17 24.65 26.36 4.57
7x7 24.05 36.81 40.75 8.44
9x9 30.78 46.48 52.87 13.06
20x20 49.17 68.77 83.54 25.40

Conclusão

Nossa análise foi feita para máscaras de tamanho 3x3, 5x5, 7x7, 9x9 e 20x20, ou seja, máscaras pequenas. Isso foi motivado pelo fato que para máscaras maiores, existem algoritmos mais eficientes para fazer a convolução, como usando a Transformada de Fourier.

Podemos ver que o desempenho do código usando texturas (código 3) foi melhor que o utilizando memória compartilhada (código 4), sendo que houve basicamente a mudança de um comando (tex2D) em relação à versão simples.

Além disso, o código 3 é muito mais limpo e claro que o 4 e é independente de especificações (como tamanho das memórias) que podem mudar com o tempo e quebrar o código.

Portanto, nem sempre o uso da memória compartilhada é a melhor forma de resolver um problema em GPUs, conforme se der o avanço do hardware e dos compiladores, a tendência é que os códigos em GPU se aproximem mais da forma 3 (e talvez até mesmo da 1 e 2) do que da 4. Apenas códigos extremamente otimizados vão se preocupar com detalhes de implementação como conflitos em memória e acesso coordenado (coalesced), que é exatamente a situação em que se encontram as CPUs (com tempo muito maior de existência no mercado).