1

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?

Amit Joshi
  • 15,448
  • 21
  • 77
  • 141
SciurusDoomus
  • 113
  • 2
  • 14
  • can you finish your question first before publishing! – Samer Jun 12 '15 at 20:24
  • Sorry! I accidentally hit Enter. All done now. – SciurusDoomus Jun 12 '15 at 20:26
  • 3
    `scipy.ndimage.filters.generic_filter(data, function=numpy.mean, size=3)` should do that rather efficiently (you can also specify the `footprint` parameter). – rth Jun 12 '15 at 20:35
  • @rth I've looked at the scipy filters but I'd like the iterator to work since the filtering isn't the only thing I want it to do. I'd like to find the standard deviation of each neighborhood and look at that, perhaps take the Laplacian and try and find edges, etc. The reason I used the mean is because it was the simplest to code in C. – SciurusDoomus Jun 12 '15 at 20:39
  • 2
    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 Jun 12 '15 at 20:45
  • @rth Thanks. I'll probably actually use these for those options and I'll dig through ndimage a little. I do still need the iterator (my boss has a laundry list of uses for it), but that's helpful. – SciurusDoomus Jun 12 '15 at 20:49
  • Actually what's the problem with the above code? Does it raise an error or is the result wrong? – rth Jun 12 '15 at 20:59
  • 1
    @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 Jun 12 '15 at 21:01

1 Answers1

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.

rth
  • 10,680
  • 7
  • 53
  • 77