3 ms·

The example there about edges on an array is something I've had to directly deal with myself when I implemented a multidimensional image-processing function for

by jordigh 10y ago

The example there about edges on an array is something I've had to directly deal with myself when I implemented a multidimensional image-processing function for GNU Octave. The problem is to find connected components in a binary image, where voxels are either 0 or 1. You want to find all the islands of ones, and you want to be flexible if diagonals count as being connected or not. 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?

arghbleargh 10y ago
Isn't it simple enough to make three bounds checks to essentially do the same thing you did with the zero padding? You can even make it a function is_in_bounds(x, y, z)
- jordigh 10y ago
  It's n-dimensional, not 3-dimensional. Also, it's really slow and wasteful to be doing bounds-checking for each voxel of an image that is almost completely not a boundary, whereas you can't avoid checking if each voxel is lit. :-)
ThatGeoGuy 10y ago
I think there are three generally accepted ways to do this: 1) Pad with zeros until you fill your extent 2) Pad with an arbitrary constant of some kind 3) Pretend that spilling over one edge loops back to the opposite edge, as if the image was just a view over a repeating signal in that dimension. The first two are naive and work well enough for simple tasks. The third actually has some basis to it, although it does seem odd. The third option is actually how the image is treated when you transform the image using a 2D Fourier transform. Instead of cutting off the "signal" and padding with zeros, it assumes that the signal is continuous and repeating (i.e. it is a 2D or N-dimensional Fourier series). I know MATLAB uses this as the default for some (irritatingly not _all_) of it's operations, especially those concerned with Fourier transforms. I can't speak to Octave, but I believe they try to aim for bug compatibility with MATLAB where possible. One nice property of this kind of behaviour is that it avoids introducing some artefacts in the image when performing some sort of kernel filter or convolution. Note that the interesting bit here is that Python + Numpy actually makes this whole debacle a bit easier to deal with. You can implement solution 3) easier than 1) or 2) with Python, since you can use negative indices in any iterable. So for an image `I`, you can call `I[xdim, ydim, -1, wdim]` for any N-dimensional image and the wrapping works well. Of course you then have to ensure you don't go past the "extent" of your N-dimensional space, but that's simple enough to enforce with 0 and `I.shape`. A lot of people complain that negative indices make no sense and are an annoying language flaw, but this is one scenario where I'm really glad they're available. Numpy views would also probably alleviate your realloc problem, since you can use views to simulate extents within larger matrices, without requiring copying / separate allocation.

by jordigh 10y ago

arghbleargh 10y ago

Isn't it simple enough to make three bounds checks to essentially do the same thing you did with the zero padding? You can even make it a function is_in_bounds(x, y, z)

jordigh 10y ago
It's n-dimensional, not 3-dimensional. Also, it's really slow and wasteful to be doing bounds-checking for each voxel of an image that is almost completely not a boundary, whereas you can't avoid checking if each voxel is lit. :-)

ThatGeoGuy 10y ago

I think there are three generally accepted ways to do this: 1) Pad with zeros until you fill your extent 2) Pad with an arbitrary constant of some kind 3) Pretend that spilling over one edge loops back to the opposite edge, as if the image was just a view over a repeating signal in that dimension. The first two are naive and work well enough for simple tasks. The third actually has some basis to it, although it does seem odd. The third option is actually how the image is treated when you transform the image using a 2D Fourier transform. Instead of cutting off the "signal" and padding with zeros, it assumes that the signal is continuous and repeating (i.e. it is a 2D or N-dimensional Fourier series). I know MATLAB uses this as the default for some (irritatingly not _all_) of it's operations, especially those concerned with Fourier transforms. I can't speak to Octave, but I believe they try to aim for bug compatibility with MATLAB where possible. One nice property of this kind of behaviour is that it avoids introducing some artefacts in the image when performing some sort of kernel filter or convolution. Note that the interesting bit here is that Python + Numpy actually makes this whole debacle a bit easier to deal with. You can implement solution 3) easier than 1) or 2) with Python, since you can use negative indices in any iterable. So for an image `I`, you can call `I[xdim, ydim, -1, wdim]` for any N-dimensional image and the wrapping works well. Of course you then have to ensure you don't go past the "extent" of your N-dimensional space, but that's simple enough to enforce with 0 and `I.shape`. A lot of people complain that negative indices make no sense and are an annoying language flaw, but this is one scenario where I'm really glad they're available. Numpy views would also probably alleviate your realloc problem, since you can use views to simulate extents within larger matrices, without requiring copying / separate allocation.