Indexación 3D del núcleo CUDA para el filtrado de imágenes?

Tengo una matriz de características de imagen A es la matriz n * m * 31 acolchada para filtrar y tengo B como un filtro de objetos k * l * 31 . Quiero obtener una matriz de salida C es p * r * 31 con el tamaño de la imagen A sin relleno . Intento escribir un código CUDA para ejecutar el filtro B sobre A y obtener C.

Supongo que cada operación de filtrado sobre A con filtro B ocupada por un bloque de hilos, por lo que habrá una operación k * l dentro de cada bloque de hilos. Y cada operación de filtrado desplazado se realizará en diferentes bloques de hilos. Para A (0,0) el filtrado estará en thread_block (0,0) y para A (0,1) estará en thread_block (1,0), etc. También tengo una tercera dimensión como 31. Cada espacio en la tercera dimensión se calculará en sí mismo. Por lo tanto, con la correcta indexación en 3D sobre las matrices, podría tener toda la operación en forma muy paralela.

Entonces la operación es

A n*m*31 XB k*l*31 = C p*r*31 

¿Cómo podría hacer indexación kernel para las operaciones de manera eficiente?

Aquí hay un código que he escrito para hacer más o menos lo que describes.

Crea un conjunto de datos 3D (llamado cell , pero sería comparable a su matriz A ), lo rellena con datos aleatorios y luego calcula una matriz resultante de salida 3D (llamada node , pero sería comparable a la matriz C ) basada en los datos en A El tamaño de datos de A es mayor que el tamaño de C (“acolchado” como lo llama) para permitir el paso de la función B sobre los elementos de contorno de A En mi caso, la función B es simplemente encontrar el mínimo dentro del volumen cúbico 3D de A que está asociado con el tamaño de B (lo que llamo WSIZE , crear una región 3D WSIZE x WSIZE x WSIZE ) y almacenar el resultado en C

Este código particular intenta explotar la reutilización de datos copiando una cierta región de la entrada A en la memoria compartida para cada bloque. Cada bloque calcula múltiples puntos de salida (es decir, calcula B una cantidad de veces para completar una región de C) a fin de explotar la oportunidad de reutilización de datos de los cómputos vecinos de B.

Esto podría ayudarlo a comenzar. Obviamente tendrías que reemplazar B (mi código de búsqueda mínima) con lo que sea tu función B deseada. Además, necesitarías modificar el dominio B de uno cúbico en mi caso a cualquier tipo de prisma rectangular que corresponda a tus dimensiones B. Esto también afectará la operación de memoria compartida, por lo que es posible que desee prescindir de la memoria compartida para su primera iteración solo para que las cosas funcionen correctamente, luego agregue la optimización de memoria compartida para ver qué beneficio puede obtener.

 #include  #include  // these are just for timing measurments #include  // Computes minimum in a 3D volume, at each output point // To compile it with nvcc execute: nvcc -O2 -o grid3d grid3d.cu //define the window size (cubic volume) and the data set size #define WSIZE 6 #define DATAXSIZE 100 #define DATAYSIZE 100 #define DATAZSIZE 20 //define the chunk sizes that each threadblock will work on #define BLKXSIZE 8 #define BLKYSIZE 8 #define BLKZSIZE 8 // for cuda error checking #define cudaCheckErrors(msg) \ do { \ cudaError_t __err = cudaGetLastError(); \ if (__err != cudaSuccess) { \ fprintf(stderr, "Fatal error: %s (%s at %s:%d)\n", \ msg, cudaGetErrorString(__err), \ __FILE__, __LINE__); \ fprintf(stderr, "*** FAILED - ABORTING\n"); \ return 1; \ } \ } while (0) // device function to compute 3D volume minimum at each output point __global__ void cmp_win(int knode[][DATAYSIZE][DATAXSIZE], const int kcell[][DATAYSIZE+(WSIZE-1)][DATAXSIZE+(WSIZE-1)]) { __shared__ int smem[(BLKZSIZE + (WSIZE-1))][(BLKYSIZE + (WSIZE-1))][(BLKXSIZE + (WSIZE-1))]; int tempnode, i, j, k; int idx = blockIdx.x*blockDim.x + threadIdx.x; int idy = blockIdx.y*blockDim.y + threadIdx.y; int idz = blockIdx.z*blockDim.z + threadIdx.z; if ((idx < (DATAXSIZE+WSIZE-1)) && (idy < (DATAYSIZE+WSIZE-1)) && (idz < (DATAZSIZE+WSIZE-1))){ smem[threadIdx.z][threadIdx.y][threadIdx.x]=kcell[idz][idy][idx]; if ((threadIdx.z > (BLKZSIZE - WSIZE)) && (idz < DATAZSIZE)) smem[threadIdx.z + (WSIZE-1)][threadIdx.y][threadIdx.x] = kcell[idz + (WSIZE-1)][idy][idx]; if ((threadIdx.y > (BLKYSIZE - WSIZE)) && (idy < DATAYSIZE)) smem[threadIdx.z][threadIdx.y + (WSIZE-1)][threadIdx.x] = kcell[idz][idy+(WSIZE-1)][idx]; if ((threadIdx.x > (BLKXSIZE - WSIZE)) && (idx < DATAXSIZE)) smem[threadIdx.z][threadIdx.y][threadIdx.x + (WSIZE-1)] = kcell[idz][idy][idx+(WSIZE-1)]; if ((threadIdx.z > (BLKZSIZE - WSIZE)) && (threadIdx.y > (BLKYSIZE - WSIZE)) && (idz < DATAZSIZE) && (idy < DATAYSIZE)) smem[threadIdx.z + (WSIZE-1)][threadIdx.y + (WSIZE-1)][threadIdx.x] = kcell[idz+(WSIZE-1)][idy+(WSIZE-1)][idx]; if ((threadIdx.z > (BLKZSIZE - WSIZE)) && (threadIdx.x > (BLKXSIZE - WSIZE)) && (idz < DATAZSIZE) && (idx < DATAXSIZE)) smem[threadIdx.z + (WSIZE-1)][threadIdx.y][threadIdx.x + (WSIZE-1)] = kcell[idz+(WSIZE-1)][idy][idx+(WSIZE-1)]; if ((threadIdx.y > (BLKYSIZE - WSIZE)) && (threadIdx.x > (BLKXSIZE - WSIZE)) && (idy < DATAYSIZE) && (idx < DATAXSIZE)) smem[threadIdx.z][threadIdx.y + (WSIZE-1)][threadIdx.x + (WSIZE-1)] = kcell[idz][idy+(WSIZE-1)][idx+(WSIZE-1)]; if ((threadIdx.z > (BLKZSIZE - WSIZE)) && (threadIdx.y > (BLKYSIZE - WSIZE)) && (threadIdx.x > (BLKXSIZE - WSIZE)) && (idz < DATAZSIZE) && (idy < DATAYSIZE) && (idx < DATAXSIZE)) smem[threadIdx.z+(WSIZE-1)][threadIdx.y+(WSIZE-1)][threadIdx.x+(WSIZE-1)] = kcell[idz+(WSIZE-1)][idy+(WSIZE-1)][idx+(WSIZE-1)]; } __syncthreads(); if ((idx < DATAXSIZE) && (idy < DATAYSIZE) && (idz < DATAZSIZE)){ tempnode = knode[idz][idy][idx]; for (i=0; i>>(d_node, d_cell); cudaCheckErrors("Kernel launch failure"); // copy output data back to host cudaMemcpy(node, d_node, ((nx*ny*nz)*sizeof(int)), cudaMemcpyDeviceToHost); cudaCheckErrors("CUDA memcpy3 failure"); t2 = clock(); t2sum = ((double)(t2-t1))/CLOCKS_PER_SEC; printf(" Device compute took %3.2f seconds. Beginning host compute.\n", t2sum); // now compute the same result on the host for (u=0; u cell[i+u][j+v][k+w]) temphnode = cell[i+u][j+v][k+w]; hnode[u][v][w] = temphnode; } t3 = clock(); t3sum = ((double)(t3-t2))/CLOCKS_PER_SEC; printf(" Host compute took %3.2f seconds. Comparing results.\n", t3sum); // and compare for accuracy for (i=0; i