CUDA中的三線性插值

Question

我需要在CUDA中執行三線性插值。這是問題的定義：

點考慮三個矢量：x[nx]，y[ny]，z[nz]和功能的矩陣值func[nx][ny][nz]，我想在x，y和z範圍之間的一些隨機點，找到函數值。

我可以在CUDA中編寫我自己的插值內核，但我想知道是否有一個已經存在的工作。

謝謝！

Answer 1

正如@Farzad所述，您可以使用紋理過濾在CUDA中執行三線性插值。 simpleTexture3D示例提供了有關如何使用它的完整示例。然而，現在它可能不是直接使用，因爲它涉及使用像OpenGL和glut這樣的庫以及其他外部依賴項，如cutil.h。

因此，我發現將上述代碼簡化爲顯示該概念的「最小尺寸」示例非常有用。正如你將會看到的那樣，代碼會載入位於名爲的文件中的外部數據，這些文件是我從上面鏈接的github頁面「借用」的。

下面的代碼將一個長方體數據插入位於其中央切片中的常規笛卡爾網格。如果一切順利，您將重建的圖像將在下面報告。

的代碼如下：

#include <stdio.h> 
#include <stdlib.h> 
#include <fstream> 

#include <cuda_runtime.h> 
#include <cuda.h> 

typedef unsigned char uchar; 

#define BLOCKSIZE 16 

float w = 0.5; // texture coordinate in z 

/********************/ 
/* CUDA ERROR CHECK */ 
/********************/ 
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); } 
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true) 
{ 
    if (code != cudaSuccess) 
    { 
     fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line); 
     if (abort) { getchar(); exit(code); } 
    } 
} 

typedef unsigned int uint; 
typedef unsigned char uchar; 

texture<uchar, 3, cudaReadModeNormalizedFloat> tex; // 3D texture 

cudaArray *d_volumeArray = 0; 

uint *d_output = NULL; 
uint *h_output = NULL; 

/************************************************/ 
/* TEXTURE-BASED TRILINEAR INTERPOLATION KERNEL */ 
/************************************************/ 
__global__ void 
d_render(uint *d_output, uint imageW, uint imageH, float w) 
{ 
    uint x = __umul24(blockIdx.x, blockDim.x) + threadIdx.x; 
    uint y = __umul24(blockIdx.y, blockDim.y) + threadIdx.y; 

    float u = x/(float) imageW; 
    float v = y/(float) imageH; 

    // read from 3D texture 
    float voxel = tex3D(tex, u, v, w); 

    if ((x < imageW) && (y < imageH)) { 
     // write output color 
     uint i = __umul24(y, imageW) + x; 
     d_output[i] = voxel*255; 
    } 
} 

void main() { 

    int N = 32; 
    int imageH = 512; 
    int imageW = 512; 

    const char* filename = "Bucky.raw"; 

    // --- Loading data from file 
    FILE *fp = fopen(filename, "rb"); 
    if (!fp) { fprintf(stderr, "Error opening file '%s'\n", filename); getchar(); return; } 

    uchar *data = (uchar*)malloc(N*N*N*sizeof(uchar)); 
    size_t read = fread(data, 1, N*N*N, fp); 
    fclose(fp); 

    printf("Read '%s', %lu bytes\n", filename, read); 

    gpuErrchk(cudaMalloc((void**)&d_output, imageH*imageW*sizeof(uint))); 

    // --- Create 3D array 
    const cudaExtent volumeSize = make_cudaExtent(N, N, N); 

    cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<uchar>(); 
    gpuErrchk(cudaMalloc3DArray(&d_volumeArray, &channelDesc, volumeSize)); 

    // --- Copy data to 3D array (host to device) 
    cudaMemcpy3DParms copyParams = {0}; 
    copyParams.srcPtr = make_cudaPitchedPtr((void*)data, volumeSize.width*sizeof(uchar), volumeSize.width, volumeSize.height); 
    copyParams.dstArray = d_volumeArray; 
    copyParams.extent = volumeSize; 
    copyParams.kind  = cudaMemcpyHostToDevice; 
    gpuErrchk(cudaMemcpy3D(&copyParams)); 

    // --- Set texture parameters 
    tex.normalized = true;      // access with normalized texture coordinates 
    tex.filterMode = cudaFilterModeLinear;  // linear interpolation 
    tex.addressMode[0] = cudaAddressModeWrap; // wrap texture coordinates 
    tex.addressMode[1] = cudaAddressModeWrap; 
    tex.addressMode[2] = cudaAddressModeWrap; 

    // --- Bind array to 3D texture 
    gpuErrchk(cudaBindTextureToArray(tex, d_volumeArray, channelDesc)); 

    // --- Launch the interpolation kernel 
    const dim3 blockSize(BLOCKSIZE, BLOCKSIZE, 1); 
    const dim3 gridSize(imageW/blockSize.x, imageH/blockSize.y); 
    d_render<<<gridSize, blockSize>>>(d_output, imageW, imageH, w); 
    gpuErrchk(cudaPeekAtLastError()); 
    gpuErrchk(cudaDeviceSynchronize()); 

    // --- Copy the interpolated data to host 
    h_output = (uint*)malloc(imageW*imageH*sizeof(uint)); 
    gpuErrchk(cudaMemcpy(h_output,d_output,imageW*imageH*sizeof(uint),cudaMemcpyDeviceToHost)); 

    std::ofstream outfile; 
    outfile.open("out_texture.dat", std::ios::out | std::ios::binary); 
    outfile.write((char*)h_output, imageW*imageH*sizeof(uint)); 
    outfile.close(); 

    getchar(); 

} 

代碼以二進制格式保存在一個out_texture.dat結果。可以從MATLAB根據

fd = fopen('out_texture.dat','r'); 
U = fread(fd,imageH*imageW,'unsigned int'); 
fclose(fd);