CUDA - 从 3D 数组中提取图层

Question

我有一个 3D 矩阵，其中 x-y 平面代表图像，z 平面代表图像层。
问题是当我尝试使用 idz 提取第一层（或其他层）时，我没有得到预期的结果。看起来数组，一旦在 CUDA 中，x、y 或 z 的索引与我预期的不同（如在 pycuda 中）。我通过下面的结果数组看到了这一点。

以下是这个迷你示例的分步过程（我使用通用 int 数字来表示我的图像以保存上传图像和整个代码）！
我在这里导入库并定义图像大小和图层...

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy
from pycuda.gpuarray import to_gpu

row = 10
column = 10
depth = 5

然后我定义我的输入 3D 数组和我的输出 2D 数组...

#--==== Input 3D Array ====---
arrayA = numpy.full((row, column, depth), 0)

#populate each layer with fixed values
for i in range(depth):
    arrayA[:,:,i] = i + 1

arrayA = arrayA.astype(numpy.uint16)
arrayA_gpu = cuda.mem_alloc(arrayA.nbytes)
cuda.memcpy_htod(arrayA_gpu, arrayA)
arrayA_Answer = numpy.empty_like(arrayA)

#--==== Output 2D array container ====---
arrayB = numpy.zeros([row, column], dtype = numpy.uint16)
arrayB_gpu = cuda.mem_alloc(arrayB.nbytes)
cuda.memcpy_htod(arrayB_gpu, arrayB)
arrayB_Answer = numpy.empty_like(arrayB)

接下来我在pycuda中定义CUDA内核和函数

mod = SourceModule("""
    __global__ void getLayer(int *arrayA, int *arrayB)
    {
        int idx = threadIdx.x + (blockIdx.x * blockDim.x); // x coordinate (numpy axis 2) 
        int idy = threadIdx.y + (blockIdx.y * blockDim.y); // y coordinate (numpy axis 1)
        int idz = 0; //The first layer, this can set in range from 0-4 
        int x_width = (blockDim.x * gridDim.x); 
        int y_width = (blockDim.y * gridDim.y); 

        arrayB[idx + (x_width * idy)] = arrayA[idx + (x_width * idy) + (x_width * y_width) * idz];
    }
    """)

func = mod.get_function("getLayer")
func(arrayA_gpu, arrayB_gpu, block=(row, column, 1), grid=(1,1))

使用标准 pycuda 命令，我提取结果（不是我所期望的）
arrayA[:,:,0] = 填充 1 的 10x10 矩阵（好）

print(arrayA_Answer[:,:,0])
[[1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]]

arrayB[:,:] = 10x10 矩阵填充以下（错误），预计等于 arrayA[:,:,0]...

print(arrayB_Answer)
[[1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]
 [1 2 3 4 5 1 2 3 4 5]]

Answer 1

正如here 所讨论的，numpy 3D 存储顺序模式是“z”（即“3rd”）索引是快速变化的索引，因为您在内存中线性前进。您的代码假定第一个索引 ("x") 是快速变化的索引。

由于您的内核已经针对高效（“合并”）加载/存储行为进行了组织，因此您可以通过在 numpy.这是一个有效的例子：

$ cat t10.py
from __future__ import print_function
import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy
from pycuda.gpuarray import to_gpu

row = 5
column = 10
depth = 10

#--==== Input 3D Array ====---
arrayA = numpy.full((row, column, depth), 0)
my_slice=numpy.int32(3)  # choose the layer
#populate each layer with fixed values
for i in range(row):
    arrayA[i,:,:] = i + 1

arrayA = arrayA.astype(numpy.int32)
arrayA_gpu = cuda.mem_alloc(arrayA.nbytes)
cuda.memcpy_htod(arrayA_gpu, arrayA)
arrayA_Answer = numpy.empty_like(arrayA)

#--==== Output 2D array container ====---
arrayB = numpy.zeros([column, depth], dtype = numpy.int32)
arrayB_gpu = cuda.mem_alloc(arrayB.nbytes)
cuda.memcpy_htod(arrayB_gpu, arrayB)
arrayB_Answer = numpy.empty_like(arrayB)

mod = SourceModule("""
    __global__ void getLayer(int *arrayA, int *arrayB, int slice)
    {
        int idx = threadIdx.x + (blockIdx.x * blockDim.x); // x coordinate (numpy axis 2)
        int idy = threadIdx.y + (blockIdx.y * blockDim.y); // y coordinate (numpy axis 1)
        int idz = slice; //The "layer"
        int x_width = (blockDim.x * gridDim.x);
        int y_width = (blockDim.y * gridDim.y);

        arrayB[idx + (x_width * idy)] = arrayA[idx + (x_width * idy) + (x_width * y_width) * idz];
    }
    """)

func = mod.get_function("getLayer")
func(arrayA_gpu, arrayB_gpu, my_slice, block=(depth, column, 1), grid=(1,1))
cuda.memcpy_dtoh(arrayB_Answer,arrayB_gpu)

print(arrayA[my_slice,:,:])

print(arrayB_Answer[:,:])
$ python t10.py
[[4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]]
[[4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]
 [4 4 4 4 4 4 4 4 4 4]]
$

请注意，我还将您对uint16 的使用更改为int32，以匹配内核类型int。