Coding a circular filter in Python

Question

I found a code snippet for making a circular filter using scipy and I'd like to understand how it works. I know there's a better one in skimage, but I'm interested in what's going on in this one.

 from scipy.ndimage.filters import generic_filter as gf

 #   Define physical shape of filter mask
 def circular_filter(image_data, radius):
     kernel = np.zeros((2*radius+1, 2*radius+1))
     y, x = np.ogrid[-radius:radius+1, -radius:radius+1]
     mask = x**2 + y**2 <= radius**2
     kernel[mask] = 1                
     filtered_image = gf(image_data, np.median, footprint = kernel)
     return filtered_image

But I'm not sure I understand perfectly what's going on. In particular, what exactly do the lines

     y, x = np.ogrid[-radius:radius+1, -radius:radius+1]
     mask = x**2 + y**2 <= radius**2
     kernel[mask] = 1

do?

I posted this as an answer to one of my previous questions, but it wasn't replied to, so I'm posting it as a new question.

Answer 1

I'm not familiar with SciPy but I'll give it a shot trying to explain the basic concepts.

This entire function's purpose is to alter the original image by applying a filter. This filter could do a lot of things, from changing the contrast of the image, or adding special effects, etc.

Let's go through the different lines:

     kernel = np.zeros((2*radius+1, 2*radius+1))

In this line, a copy of the image data is being created, but with all the data being zeros (hence the zeros function is being used). This is so the mask can be applied later onto it.

     y, x = np.ogrid[-radius:radius+1, -radius:radius+1]

This is creating what is known as a "meshgrid" or a multi-dimensional grid. This is to create the circular "mask". Just like how on a graph, x and y axes have evenly spaced scaling, the same is necessary here in the meshgrid. The x and y variables in this case store evenly spaced values that serve as the axes' scaling.

     mask = x**2 + y**2 <= radius**2

Here, a "mask" is being created. A mask will serve as the region in the image to be protected from the filter, so as to not alter any original data. Notice how x and y variables are used here in a Pythagorean inequality (important to see that it's not just a circle but a disk), just like how they would be in a mathematical sense. This will create a disk with the given radius that is now considered the mask. The mask variable now contains all coordinates (x,y) where the original data values should not be altered.

     kernel[mask] = 1 

This is where the mask is now applied to the copy of the image that was created earlier. Now, there is a perfect copy of the image (i.e. same dimensions) but with a disk-like "mask" that "protects" the original data from being altered. This is why all the points covered by the disk is set to 1. Also, notice how the dimensions of kernel and mask match. Both are multi-dimensional. The rest of values in the image copy are still set to zero, as was done in the first line.

     filtered_image = gf(image_data, np.median, footprint = kernel)

This is final part where everything is pieced together. There is the original data stored in image_data and there is the kernel, which is the image copy with the mask applied on it indicating where the data should not be altered. Both of them are passed as parameters into the actual filter function gf (stands for generic filter) and the output is a new filtered image.

This is a core concept in image filtering and if you want to learn more about it, I suggest starting out by learning basic signal processing concepts. Signal processing courses cover the mathematics of how these concepts work, but are usually explained in really abstract mathematics because this concept can be applied to numerous different examples.

Answer 2

Looking at your code in detail:

 kernel = np.zeros((2*radius+1, 2*radius+1))
 y, x = np.ogrid[-radius:radius+1, -radius:radius+1]
 mask = x**2 + y**2 <= radius**2
 kernel[mask] = 1                

The first line:

 kernel = np.zeros((2*radius+1, 2*radius+1))

creates a 2-d array of zeros, with a center point and "radius" points on either side. For radius = 2, you would get:

# __r__ +1 __r__
[ 0, 0, 0, 0, 0, ] #\
[ 0, 0, 0, 0, 0, ] #_} r
[ 0, 0, 0, 0, 0, ] #   +1
[ 0, 0, 0, 0, 0, ] #\
[ 0, 0, 0, 0, 0, ] #_} r

Next, you get two arrays from the open mesh grid created by numpy.ogrid. Mesh grids are a "trick" in numpy that involves storing a "parallel" array or matrix that holds the x or y coordinate of a particular cell at the location of that cell.

For example, a y-mesh grid might look like this:

 [ 0, 0, 0, 0, 0, ]
 [ 1, 1, 1, 1, 1, ]
 [ 2, 2, 2, 2, 2, ]
 [ 3, 3, 3, 3, 3, ]
 [ 4, 4, 4, 4, 4, ]

And an x-mesh grid might look like this:

 [ 0, 1, 2, 3, 4, ]
 [ 0, 1, 2, 3, 4, ]
 [ 0, 1, 2, 3, 4, ]
 [ 0, 1, 2, 3, 4, ]
 [ 0, 1, 2, 3, 4, ]

If you look at them, you'll realize that Y_grid[x][y] == y and X_grid[x][y] == x which is so often useful that it has more than one numpy function to support it. ;-)

An open mesh grid is similar to a closed one, except that it only has "one dimension." That is, instead of a pair of (for example) 5x5 arrays, you get a 1x5 array and a 5x1 array. That's what ogrid does - it returns two open grids. The values are from -radius to radius+1, according to python rules (meaning the radius+1 is left out):

 y, x = np.ogrid[-radius:radius+1, -radius:radius+1]

So y is a numpy array storing from e.g., -2..2 (inclusive), and x is an array from -2..2 inclusive. The next step is to build a boolean mask - that is, an array full of boolean values. As you know, when you operate on a numpy array, you get another numpy array. So involving two arrays in an expression with a constant produces another array:

 mask = x**2 + y**2 <= radius**2

The value of mask is going to be a 2-color bitmap, where one color is "True" and the other color is "False." The bitmap will describe a solid circle or disk. (Because of the <= relation. Remember that x and y contain -2..2, not 0..4.)

Finally, you convert from type Boolean to int by using the masking array as an overlay on the kernel array (of zeroes), setting the zeroes to ones whenever the mask is "True":

 kernel[mask] = 1

At this point, kernel looks like:

# __r__ +1 __r__
[ 0, 0, 1, 0, 0, ] #\
[ 0, 1, 1, 1, 0, ] #_} r
[ 1, 1, 1, 1, 1, ] #   +1
[ 0, 1, 1, 1, 0, ] #\
[ 0, 0, 1, 0, 0, ] #_} r