3D numpy array iteration using shared C libraries

Question

I'm currently working on an image processing project. I'm using Python, the SimpleITK, numpy, and a couple of other libraries to take a stack of DICOM images, turn them into a 3D numpy array, and then do some image processing operations using the SITK or other mathematical techniques (masking, etc.)

Right now, I'm trying to make an averaging filter that takes the average of a 3x3 neighborhood and replaced the center pixel of that neighborhood with that average value. The result is just a blurred image. Since Python's not really good at looping through 300x300x400 pixels really fast, I'm trying to use a C library to do it for me. The problem is, I'm not good at C. (or python for that matter...)

Below is my C code:

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);
            }
        }
    }
}

And here's some of the python code I'm using. poststack is a large 3D numpy array.

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

I should mention that I looked at this StackOverflow post and I don't see what I'm doing different from the guy who asked the original question...

Passing Numpy arrays to a C function for input and output

I know I'm making a stupid mistake, but what is it?

Answer 1

The problem is that the C code produce a segmentation fault, because it tries to access kern[m*N + n*N] where indexes are outside of the the allocated array boundaries.

The indexing of your multidimensional arrays is wrong. For an array X of shape (n, m) the equivalent of X[i][j] for a flattened array in C is X[i*m + j], not the way you were using it in the code above.