The problem, of course, is that along the edges you don't want to check for neighbours outside of the image boundary. I found no good way to do this for an n-dimensional image, after a few attempts of writing very complicated code. In the end, I ended up padding the whole image with a boundary of zeros, iterating over the padded image with an "if(current_voxel)" conditional that skipped checking for lit neighbours around the boundary, and when checking for lit neighbours at the original image's boundaries would give no neighbours at the padded zero boundaries.
The code was cleaner, but I incurred a big realloc, because N-dimensional images in Octave are stored as contiguous Fortran-order (column-major) arrays. I'm still looking for a better solution to this problem.
So, how do you do this cleanly?