使用共享 C 库的 3D numpy 数组迭代

Question

我目前正在从事一个图像处理项目。我正在使用 Python、SimpleITK、numpy 和其他几个库来获取一堆 DICOM 图像，将它们转换为 3D numpy 数组，然后使用 SITK 或其他数学技术（掩蔽、等等）

现在，我正在尝试制作一个平均过滤器，它采用 3x3 邻域的平均值，并用该平均值替换该邻域的中心像素。结果只是一个模糊的图像。由于 Python 并不擅长快速循环遍历 300x300x400 像素，因此我正在尝试使用 C 库来为我做这件事。问题是，我不擅长 C。（或者 Python 就这点...）

下面是我的 C 代码：

int i, j, k, m, n, p;
double kernsum;

void iter(double *data, int N, int height, int width, int depth, double *kernavg){
    double kern[N*N];
    for (k = 0; k < depth; k++){
        for (i = (N - 1)/2; i < height - 1; i++){
            for (j = (N - 1)/2; j < width - 1; j++){                
                for (m = i - (N - 1)/2; m < i + (N - 1)/2; m++){
                    for (n = j - (N - 1)/2; n < j + (N - 1)/2; n++){
                        kern[m + n*N] = data[i + j*width + k*width*depth];
                    }
                }       
                kernsum = 0;
                for (p = 0; p < N*N; p++){
                    kernsum += kern[p];
                }
                kernavg[i + j*width + k*width*depth] = kernsum/(N*N);
            }
        }
    }
}

这是我正在使用的一些 python 代码。 poststack 是一个大型 3D numpy 数组。

height = poststack.shape[1]
width = poststack.shape[2]
depth = poststack.shape[0]
N = 3

kernavgimg = np.zeros(poststack.shape, dtype = np.double)
lib = ctypes.cdll.LoadLibrary('./iter.so')
iter = lib.iter

iter(ctypes.c_void_p(poststack.ctypes.data), ctypes.c_int(N), 
    ctypes.c_int(height), ctypes.c_int(width), ctypes.c_int(depth),
    ctypes.c_void_p(kernavgimg.ctypes.data))

print kernavgimg

pyplot.imshow(kernavgimg[0, :, :], cmap = 'gray')
pyplot.show()
image.imsave('/media/sd/K/LabCode/python_code/dump/test.png', kernavgimg.data[0, :, :], cmap = 'gray')

pyplot.imshow(poststack[0, :, :], cmap = 'gray')
pyplot.show()
image.imsave('/media/sd/K/LabCode/python_code/dump/orig.png', poststack[0, :, :], cmap = 'gray')

print kernavgimg[0, :, :] == poststack[0, :, :]
print kernavgimg.shape
print poststack.shape

我应该提到，我查看了这个 StackOverflow 帖子，我看不出我在做什么与提出原始问题的人不同......

Passing Numpy arrays to a C function for input and output

我知道我犯了一个愚蠢的错误，但这是什么？

Answer 1

问题在于 C 代码会产生分段错误，因为它会尝试访问 kern[m*N + n*N]，其中索引位于分配的数组边界之外。

您的多维数组的索引是错误的。对于形状为(n, m) 的数组X，C 中扁平数组的X[i][j] 的等价物是X[i*m + j]，而不是您在上面的代码中使用它的方式。