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?

scipy.ndimage.filters.generic_filter(data, function=numpy.mean, size=3) should do that rather efficiently (you can also specify the footprint parameter). — rth
– rth, Commented Jun 12, 2015 at 20:35
You can calculate the local standard deviation (using function=numpy.std) with the same approach, and there is a laplacian function in scipy.ndimage. Anyway, hope you find your issue with ctypes. Most of the classical image processing operators are already efficiently implemented in scipy, scikit-image, etc. — rth
– rth, Commented Jun 12, 2015 at 20:45
@rth The code fails silently (I think) after the iter function call and does not print kernavgimg in Python nor does it display the images, save them to disk, etc. — SciurusDoomus
– SciurusDoomus, Commented Jun 12, 2015 at 21:01

rth · Accepted Answer · 2015-07-06 08:42:51Z

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.

edited Jul 6, 2015 at 8:42

answered Jun 13, 2015 at 23:45

rth

11.3k7 gold badges58 silver badges88 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

SciurusDoomus Over a year ago

I've changed my code in the OP, but I'm still getting the same error... Any ideas?

Collectives™ on Stack Overflow

3D numpy array iteration using shared C libraries

1 Answer 1

1 Comment

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Linked

Related