Convolución unidimensional

La convolución es un patrón de cálculo muy común. Su método se basa en que cada elemento de los datos de salida es la suma ponderada del elemento de entrada en la misma posición y sus elementos circundantes. La convolución utiliza un segmento de datos llamado kernel o máscara, que especifica los pesos de los elementos vecinos. Para evitar confusiones, nos referiremos a él como máscara. El cálculo se ilustra a continuación:

El vector M es la máscara, el vector N es la entrada y el vector P es la salida. Por ejemplo, P[2] = N[0] * M[0] + ... + N[4] * M[4]. El cálculo de la convolución debe considerar los límites. Como se muestra, con una máscara de longitud 5, al calcular los primeros y últimos dos elementos, las posiciones de los datos de entrada necesarios pueden superar el límite.

Para manejar esto, se deben agregar manualmente elmeentos conocidos como Celdas Fantasma, cuyo valor depende del caso. Para simplificar, aquí se utiliza el valor 0. En algunas aplicaciones, este valor se establece como P[0].

Código del kernel para convolución 1D

__global__ void convolucion1D_kernel(float* inputSignal, float* maskData, float* outputSignal, int maskLen, int signalLen) {
    int globalIdx = blockIdx.x * blockDim.x + threadIdx.x;
    float resultado = 0.0f;
    int inicio = globalIdx - (maskLen / 2);
    for(int offset = 0; offset < maskLen; offset++) {
        int posActual = inicio + offset;
        if(posActual >= 0 && posActual < signalLen) {
            resultado += inputSignal[posActual] * maskData[offset];
        }
    }
    outputSignal[globalIdx] = resultado;
}

Convolución bidimensional

La convolución bidimensional es una extensión de la unidimensional. Las Celdas Fantasma resultantes también son bidimensionales.

Uso de Tile para calcular la convolución

Optimización 1D con tile

En la convolución 1D, al obtener valores de N, algunos se utilizen múltiples veces. Por ejemplo:

Para calcular P4 se necesitan N2, N3, N4, N5, N6. Para P5 se necesitan N3, N4, N5, N6, N7. Para P7 se necesitan N5, N6, N7, N8, N9.

Esto revela mucha duplicación en los cálculos. Si se divide en tiles de longitud 4, como se muestra, para calcular P4, P5, P6, P7 se requieren N2 a N9.

Almacenar estos valores en memoria compartida reduce significativamente el número de accesos a memoria global.

Implementación del kernel 1D con tile

#define TAM_TILE 4
#define TAM_MASCARA 5
#define EXTENSION TAM_MASCARA / 2
#define TAM_BLOQUE (TAM_TILE + TAM_MASCARA - 1)

// Configuración de la malla de ejecución
dim3 dimensionBloque(TAM_BLOQUE, 1, 1);
dim3 dimensionMalla((senalLong - 1) / TAM_TILE + 1, 1, 1);

// Fragmento del kernel
__global__ void convolucion1D_tile(float* senal, const float* __restrict__ mascara, float* resultado, int mascaraLen, int senalLen) {
    __shared__ float cacheDatos[TAM_BLOQUE];
    int idxLocal = threadIdx.x;
    int idxGlobal = blockIdx.x * TAM_TILE + idxLocal - EXTENSION;
    
    // Carga de datos a memoria compartida
    if(idxGlobal >= 0 && idxGlobal < senalLen) {
        cacheDatos[idxLocal] = senal[idxGlobal];
    } else {
        cacheDatos[idxLocal] = 0.0f;
    }
    __syncthreads();
    
    // Cálculo de la convolución para el tile actual
    if(idxLocal < TAM_TILE) {
        float acumulador = 0.0f;
        for(int k = 0; k < mascaraLen; k++) {
            acumulador += mascara[k] * cacheDatos[k + idxLocal];
        }
        resultado[blockIdx.x * TAM_TILE + idxLocal] = acumulador;
    }
}

Optimización 2D con tile

La convolución bidimensional se basa en la unidimensional. La máscara M, al ser constante, puede almacenarse en memoria constante para mejorar la velocidad de acceso. La memoria constante se distribuye en un warp mediante broadcast a todos los hilos.

La convolución 2D se aplica comúnmente en procesameinto de imágenes.

Primera implementación 2D con tile

#include <wb.h>

#define ANCHO_MASCARA 5
#define RADIO_MASCARA (ANCHO_MASCARA / 2)
#define ANCHO_TILE 16
#define ANCHO_BLOQUE (ANCHO_TILE + ANCHO_MASCARA - 1)
#define LIMITAR(x) (min(max((x), 0.0), 1.0))

__global__ void convolucion2D_kernel(float* imagen, const float* __restrict__ mascara, float* salida,
                                     int canales, int ancho, int alto) {
    __shared__ float bloqueDatos[ANCHO_BLOQUE][ANCHO_BLOQUE];
    int tx = threadIdx.x;
    int ty = threadIdx.y;
    int filaSalida = blockIdx.y * ANCHO_TILE + ty;
    int colSalida = blockIdx.x * ANCHO_TILE + tx;
    int filaEntrada = filaSalida - RADIO_MASCARA;
    int colEntrada = colSalida - RADIO_MASCARA;
    
    for(int canal = 0; canal < canales; canal++) {
        // Carga de datos considerando límites
        if(filaEntrada >= 0 && filaEntrada < alto && colEntrada >= 0 && colEntrada < ancho) {
            bloqueDatos[ty][tx] = imagen[(filaEntrada * ancho + colEntrada) * canales + canal];
        } else {
            bloqueDatos[ty][tx] = 0.0f;
        }
        __syncthreads();
        
        // Cálculo de la convolución para el tile actual
        float valorSalida = 0.0f;
        if(ty < ANCHO_TILE && tx < ANCHO_TILE) {
            for(int i = 0; i < ANCHO_MASCARA; i++) {
                for(int j = 0; j < ANCHO_MASCARA; j++) {
                    valorSalida += mascara[j * ANCHO_MASCARA + i] * bloqueDatos[i + ty][j + tx];
                }
            }
            if(filaSalida < alto && colSalida < ancho) {
                salida[(filaSalida * ancho + colSalida) * canales + canal] = LIMITAR(valorSalida);
            }
        }
        __syncthreads();
    }
}

int main(int argc, char* argv[]) {
    wbArg_t argumentos;
    int mascFilas, mascColumnas;
    int imgCanales, imgAncho, imgAlto;
    char* archivoImagen;
    char* archivoMascara;
    wbImage_t imagenEntrada, imagenSalida;
    float *datosImgCPU, *datosImgSalidaCPU, *datosMascaraCPU;
    float *datosImgGPU, *datosImgSalidaGPU, *datosMascaraGPU;
    
    argumentos = wbArg_read(argc, argv);
    archivoImagen = wbArg_getInputFile(argumentos, 0);
    archivoMascara = wbArg_getInputFile(argumentos, 1);
    
    imagenEntrada = wbImport(archivoImagen);
    datosMascaraCPU = (float*)wbImport(archivoMascara, &mascFilas, &mascColumnas);
    
    assert(mascFilas == 5 && mascColumnas == 5);
    
    imgAncho = wbImage_getWidth(imagenEntrada);
    imgAlto = wbImage_getHeight(imagenEntrada);
    imgCanales = wbImage_getChannels(imagenEntrada);
    
    imagenSalida = wbImage_new(imgAncho, imgAlto, imgCanales);
    datosImgCPU = wbImage_getData(imagenEntrada);
    datosImgSalidaCPU = wbImage_getData(imagenSalida);
    
    // Asignación de memoria en GPU
    cudaMalloc(&datosImgGPU, imgAncho * imgAlto * imgCanales * sizeof(float));
    cudaMalloc(&datosImgSalidaGPU, imgAncho * imgAlto * imgCanales * sizeof(float));
    cudaMalloc(&datosMascaraGPU, mascFilas * mascColumnas * sizeof(float));
    
    // Copia de datos a GPU
    cudaMemcpy(datosImgGPU, datosImgCPU, imgAncho * imgAlto * imgCanales * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(datosMascaraGPU, datosMascaraCPU, mascFilas * mascColumnas * sizeof(float), cudaMemcpyHostToDevice);
    
    // Ejecución del kernel
    dim3 dimMalla(ceil((float)imgAncho / ANCHO_TILE), ceil((float)imgAlto / ANCHO_TILE));
    dim3 dimBloque(ANCHO_BLOQUE, ANCHO_BLOQUE, 1);
    convolucion2D_kernel<<<dimmalla dimbloque="">>>(datosImgGPU, datosMascaraGPU, datosImgSalidaGPU,
                                                  imgCanales, imgAncho, imgAlto);
    cudaDeviceSynchronize();
    
    // Copia de resultados a CPU
    cudaMemcpy(datosImgSalidaCPU, datosImgSalidaGPU, imgAncho * imgAlto * imgCanales * sizeof(float), cudaMemcpyDeviceToHost);
    
    wbSolution(argumentos, imagenSalida);
    
    // Liberación de memoria
    cudaFree(datosImgGPU);
    cudaFree(datosImgSalidaGPU);
    cudaFree(datosMascaraGPU);
    free(datosMascaraCPU);
    wbImage_delete(imagenSalida);
    wbImage_delete(imagenEntrada);
    
    return 0;
}</dimmalla></wb.h>

Segunda implementación 2D con tile

La diferencia con la implementación anterior radica en la configuración del bloque. Aquí, el tamaño del bloque coincide con el del tile, por lo que la carga de datos de entrada requiere dos etapas, ya que la cantidad de datos necesarios excede la capacidad de un solo bloque.

#include <wb.h>

#define ANCHO_MASCARA 5
#define RADIO_MASCARA (ANCHO_MASCARA / 2)
#define TAM_TILE 16
#define TAM_BLOQUE (TAM_TILE + ANCHO_MASCARA - 1)
#define RESTRINGIR(x) (min(max((x), 0.0), 1.0))

__global__ void convolucion2D_avanzada(float* imagen, const float* __restrict__ mascara, float* salida,
                                       int canales, int ancho, int alto) {
    __shared__ float datosCompartidos[TAM_BLOQUE][TAM_BLOQUE];
    
    for(int canal = 0; canal < canales; canal++) {
        // Primera etapa de carga
        int idxCarga1 = threadIdx.y * TAM_TILE + threadIdx.x;
        int filaDest1 = idxCarga1 / TAM_BLOQUE;
        int colDest1 = idxCarga1 % TAM_BLOQUE;
        int filaOrig1 = blockIdx.y * TAM_TILE + filaDest1 - RADIO_MASCARA;
        int colOrig1 = blockIdx.x * TAM_TILE + colDest1 - RADIO_MASCARA;
        
        if(filaOrig1 >= 0 && filaOrig1 < alto && colOrig1 >= 0 && colOrig1 < ancho) {
            datosCompartidos[filaDest1][colDest1] = imagen[(filaOrig1 * ancho + colOrig1) * canales + canal];
        } else {
            datosCompartidos[filaDest1][colDest1] = 0.0f;
        }
        
        // Segunda etapa de carga
        int idxCarga2 = idxCarga1 + TAM_TILE * TAM_TILE;
        int filaDest2 = idxCarga2 / TAM_BLOQUE;
        int colDest2 = idxCarga2 % TAM_BLOQUE;
        int filaOrig2 = blockIdx.y * TAM_TILE + filaDest2 - RADIO_MASCARA;
        int colOrig2 = blockIdx.x * TAM_TILE + colDest2 - RADIO_MASCARA;
        
        if(filaDest2 < TAM_BLOQUE) {
            if(filaOrig2 >= 0 && filaOrig2 < alto && colOrig2 >= 0 && colOrig2 < ancho) {
                datosCompartidos[filaDest2][colDest2] = imagen[(filaOrig2 * ancho + colOrig2) * canales + canal];
            } else {
                datosCompartidos[filaDest2][colDest2] = 0.0f;
            }
        }
        __syncthreads();
        
        // Cálculo de la convolución
        float resultado = 0.0f;
        for(int fy = 0; fy < ANCHO_MASCARA; fy++) {
            for(int fx = 0; fx < ANCHO_MASCARA; fx++) {
                resultado += datosCompartidos[threadIdx.y + fy][threadIdx.x + fx] * mascara[fy * ANCHO_MASCARA + fx];
            }
        }
        
        int filaFinal = blockIdx.y * TAM_TILE + threadIdx.y;
        int colFinal = blockIdx.x * TAM_TILE + threadIdx.x;
        if(filaFinal < alto && colFinal < ancho) {
            salida[(filaFinal * ancho + colFinal) * canales + canal] = RESTRINGIR(resultado);
        }
        __syncthreads();
    }
}

int main(int argc, char* argv[]) {
    wbArg_t arg;
    int mascFilas, mascColumnas;
    int imgCanales, imgAncho, imgAlto;
    char* archivoImgEntrada;
    char* archivoMascara;
    wbImage_t imgEntrada, imgSalida;
    float *datosImgCPU, *datosSalidaCPU, *datosMascaraCPU;
    float *datosImgGPU, *datosSalidaGPU, *datosMascaraGPU;
    
    arg = wbArg_read(argc, argv);
    archivoImgEntrada = wbArg_getInputFile(arg, 0);
    archivoMascara = wbArg_getInputFile(arg, 1);
    
    imgEntrada = wbImport(archivoImgEntrada);
    datosMascaraCPU = (float*)wbImport(archivoMascara, &mascFilas, &mascColumnas);
    
    assert(mascFilas == 5 && mascColumnas == 5);
    
    imgAncho = wbImage_getWidth(imgEntrada);
    imgAlto = wbImage_getHeight(imgEntrada);
    imgCanales = wbImage_getChannels(imgEntrada);
    
    imgSalida = wbImage_new(imgAncho, imgAlto, imgCanales);
    datosImgCPU = wbImage_getData(imgEntrada);
    datosSalidaCPU = wbImage_getData(imgSalida);
    
    // Asignación de memoria en GPU
    cudaMalloc(&datosImgGPU, imgAncho * imgAlto * imgCanales * sizeof(float));
    cudaMalloc(&datosSalidaGPU, imgAncho * imgAlto * imgCanales * sizeof(float));
    cudaMalloc(&datosMascaraGPU, mascFilas * mascColumnas * sizeof(float));
    
    // Transferencia de datos a GPU
    cudaMemcpy(datosImgGPU, datosImgCPU, imgAncho * imgAlto * imgCanales * sizeof(float), cudaMemcpyHostToDevice);
    cudaMemcpy(datosMascaraGPU, datosMascaraCPU, mascFilas * mascColumnas * sizeof(float), cudaMemcpyHostToDevice);
    
    // Ejecución del kernel
    dim3 dimMalla(ceil((float)imgAncho / TAM_TILE), ceil((float)imgAlto / TAM_TILE));
    dim3 dimBloque(TAM_TILE, TAM_TILE, 1);
    convolucion2D_avanzada<<<dimmalla dimbloque="">>>(datosImgGPU, datosMascaraGPU, datosSalidaGPU,
                                                   imgCanales, imgAncho, imgAlto);
    
    // Transferencia de resultados a CPU
    cudaMemcpy(datosSalidaCPU, datosSalidaGPU, imgAncho * imgAlto * imgCanales * sizeof(float), cudaMemcpyDeviceToHost);
    
    wbSolution(arg, imgSalida);
    
    // Liberación de memoria
    cudaFree(datosImgGPU);
    cudaFree(datosSalidaGPU);
    cudaFree(datosMascaraGPU);
    free(datosMascaraCPU);
    wbImage_delete(imgSalida);
    wbImage_delete(imgEntrada);
    
    return 0;
}</dimmalla></wb.h>