1

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

现在，我正在尝试制作一个平均滤波器，它取 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 帖子，我看不出我在做什么与提出原始问题的人不同......

将 Numpy 数组传递给 C 函数以进行输入和输出

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

4

1 回答 1

1

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

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

于 2015-06-13T23:45:11.857 回答