Wednesday, April 8th, 2026

Last class, we continued work on applying the mean/median filter to remove noise from a grayscale image.

import numpy as np
import matplotlib.pyplot as plt

Project 4: Image denoising

If you have a working sp_noise function, then you can use any grayscale image you'd like and add salt and pepper noise to generate a corresponding noisy_img array. Since I've not written such a function in class, I will use the pre-noised image (and the corresponding clean version) from the project page.

img = np.mean(plt.imread('everton.png')[:,:,:3], axis=2)
noisy_img = np.mean(plt.imread('noisy_img.png')[:,:,:3], axis=2)

We previously developed some code that would apply the mean filter. If you've already written a mean_filter function, you should use that instead.

nrows, ncols = noisy_img.shape

mean_filtered_img = noisy_img.copy()
for i in range(1,nrows-1):
    for j in range(1,ncols-1):
        grid = noisy_img[i-1:i+2, j-1:j+2]
        mean_filtered_img[i,j] = np.mean(grid)

The median filter works in much the same way. If you've already written a median_filter function, you should use that instead.

nrows, ncols = noisy_img.shape

median_filtered_img = noisy_img.copy()
for i in range(1,nrows-1):
    for j in range(1,ncols-1):
        grid = noisy_img[i-1:i+2, j-1:j+2]
        median_filtered_img[i,j] = np.median(grid)

Last class, we also discussed using a padded version of our image array to deal with pixels along the edge.

nrows, ncols = noisy_img.shape

padded_img = np.zeros((nrows + 2, ncols + 2))   # Initialize padded array with two extra rows/columns
padded_img[1:-1, 1:-1] = noisy_img              # Set the noisy_img array into the center of padded_img

median_filtered_img = noisy_img.copy()
for i in range(nrows):
    for j in range(ncols):
        grid = padded_img[i:i+3, j:j+3]
        median_filtered_img[i,j] = np.median(grid)

Let's see how this filter peforms.

plt.figure(figsize=(8,4))

plt.subplot(1,2,1)
plt.imshow(noisy_img, cmap='gray', vmin=0, vmax=1)
plt.axis('off')
plt.title('Noisy image')

plt.subplot(1,2,2)
plt.imshow(median_filtered_img, cmap='gray', vmin=0, vmax=1)
plt.axis('off')
plt.title('Median filtered image')

plt.tight_layout()

Reminder: The exercises from previous classes ask you to create functions to accomplish the above tasks. For example:

• The sp.noise function should take in a grayscale image array and a noise parameter and create a copy of the image array with salt and pepper noise added (with approximately $100\cdot$noise% of the pixels being hit with salt/pepper noise).
• The get_padded_img function should take in a grayscale array and a desired pad size and create a padded array with pad additional rows at the top/bottom and pad additional columns at the left/right of the image array.
• The mean_filter function should take in a grayscale image array and apply the mean filter, accounting for edge pixels (using the get_padded_img function).
• The median_filter function should take in a grayscale image array and apply the median filter, accounting for edge pixels (using the get_padded_img function).
• You may choose to combine the mean_filter and median_filter functions into a single function, which takes in an additional argument filter_func that controls the type of filter (e.g. filter_func = np.mean or filter_func = np.median).

Handling different grid sizes

So far, we've been applying the mean/median filters using specifically a 3 by 3 grid centered on each pixel. What if we wan to use grids of different sizes, say s=3, s=5, s=7,...?

Let's think through how we would change our code to treat 5 by 5 grids. For now, let's again ignore the filtering of edge pixels.

nrows, ncols = noisy_img.shape

median_filtered_img = noisy_img.copy()
for i in range(2,nrows-2):
    for j in range(2,ncols-2):
        grid = noisy_img[i-2:i+2+1, j-2:j+2+1]
        median_filtered_img[i,j] = np.median(grid)

plt.imshow(median_filtered_img, cmap='gray', vmin=0, vmax=1)
plt.axis('off')

(np.float64(-0.5), np.float64(699.5), np.float64(494.5), np.float64(-0.5))

What about 7 by 7 grids?

nrows, ncols = noisy_img.shape

median_filtered_img = noisy_img.copy()
for i in range(3,nrows-3):
    for j in range(3,ncols-3):
        grid = noisy_img[i-3:i+3+1, j-3:j+3+1]
        median_filtered_img[i,j] = np.median(grid)

plt.imshow(median_filtered_img, cmap='gray', vmin=0, vmax=1)
plt.axis('off')

(np.float64(-0.5), np.float64(699.5), np.float64(494.5), np.float64(-0.5))

Note: The grid size s and the pad size pad are related by pad = (s-1)//2 or s = 2*pad + 1.

Exercise: Modify the examples above to take in a parameter s, which is an odd integer, and apply the mean/median filter using an s by s grid. You can skip filtering any rows/columns along the edge as necessary.

Putting things together

We've discussed how to deal with edge pixels (for the special case of 3 by 3 grids) and how to deal with varying grid sizes. We want to be able to handle both at the same time.

Exercise: Combine the two discussions (dealing with different grid sizes and dealing with edge pixels) to write functions mean_filter(img, s) and median_filter(img,s) that apply the mean/median filters to an input 2D array img using s by s grids, including all edge pixels.

Remaining tasks to be addressed for Project 4

Once you have working mean/median filters that deal with edge pixels using any odd-sized grid, you can explore how well they work on images with varying levels of noise. Does the noise level affect which filter works better? What about the grid size? Does the "best" grid size change when the amount of noise changes?

Additionally, try to think of other ways the images could be filtered. You can start from the mean/median filters and think of ways to modify or improve them, or come up with something new altogether.

Speed - "NumPy-ifying" the process

For large images, the median filter code above is somewhat slow. It would be nice if we could use NumPy to avoid the nested for loop iteration through every row and column.

For reference, let's time how long it takes to apply the median filter to an image. I will use the sample code at the top of the notebook, but you should use your median_filter function instead.

from time import time

t0 = time()

nrows, ncols = noisy_img.shape

padded_img = np.zeros((nrows + 2, ncols + 2))   # Initialize padded array with two extra rows/columns
padded_img[1:-1, 1:-1] = noisy_img              # Set the noisy_img array into the center of padded_img

median_filtered_img = noisy_img.copy()
for i in range(nrows):
    for j in range(ncols):
        grid = padded_img[i:i+3, j:j+3]
        median_filtered_img[i,j] = np.median(grid)

t1 = time()

print('Execution time:', t1-t0)

Execution time: 13.136903047561646

Idea: Let's construct nine versions of our array. In particular, we want the following variations of the noisy_img array:

1. The noisy_img array,
2. The noisy_img array shifted down by one pixel,
3. The noisy_img array shifted up by one pixel,
4. The noisy_img array shifted left by one pixel,
5. The noisy_img array shifted right by one pixel,
6. The noisy_img array shifted up by one pixel and left by one pixel,
7. The noisy_img array shifted up by one pixel and right by one pixel,
8. The noisy_img array shifted down by one pixel and left by one pixel, and
9. The noisy_img array shifted down by one pixel and right by one pixel.

Then if we collect the [i,j]th entry of each of these shifted versions, we have all the values that makeup our 3 by 3 grid to compute a median.

How can we construct each of these shifted versions? We will again want to use a padded version of our array to make room for each of these shifts. We can then just take various slices of the padded_img array to construct these shifted versions.

nrows, ncols = noisy_img.shape
padded_img = np.zeros((nrows + 2, ncols + 2))
padded_img[1:-1, 1:-1] = noisy_img

pnrows, pncols = padded_img.shape

shift1 = padded_img[1:pnrows-1, 1:pncols-1]  # Unshifted noisy_img
shift2 = padded_img[0:pnrows-2, 1:pncols-1]  # Shifted down by 1
shift3 = padded_img[2:pnrows-0, 1:pncols-1]  # Shifted up by 1
shift4 = padded_img[1:pnrows-1, 2:pncols-0]  # Shifted left by 1
shift5 = padded_img[1:pnrows-1, 0:pncols-2]  # Shifted right by 1
shift6 = padded_img[2:pnrows-0, 2:pncols-0]  # Shifted up by 1, left by 1
shift7 = padded_img[2:pnrows-0, 0:pncols-2]  # Shifted up by 1, right by 1
shift8 = padded_img[0:pnrows-2, 2:pncols-0]  # Shifted down by 1, left by 1
shift9 = padded_img[0:pnrows-2, 0:pncols-2]  # Shifted down by 1, right by 1

Now that we have all of these shifted versions, we need some way to compute medians element-wise through these nine shifts. The np.stack function can stack several arrays together into a new dimension. The desired new dimension can be specified with the axis input argument. For example, calling np.stack with axis=-1 will place the new dimension after any existing dimensions.

In particular, we can use np.stack with axis=-1 to combine these nine shifted arrays into a 3-dimensional array with shape (num_rows,num_cols, 9), and then use np.median to compute the median through axis=-1.

shifts = [shift1, shift2, shift3, shift4, shift5, shift6, shift7, shift8, shift9]

stacked_shifts = np.stack(shifts, axis=-1)
print(stacked_shifts.shape)

(495, 700, 9)

noisy_img.shape

(495, 700)

median_filtered_img = np.median(stacked_shifts, axis=-1)

Let's time it:

t0 = time()

nrows, ncols = noisy_img.shape
padded_img = np.zeros((nrows + 2, ncols + 2))
padded_img[1:-1, 1:-1] = noisy_img

pnrows, pncols = padded_img.shape

shift1 = padded_img[1:pnrows-1, 1:pncols-1]  # Unshifted noisy_img
shift2 = padded_img[0:pnrows-2, 1:pncols-1]  # Shifted down by 1
shift3 = padded_img[2:pnrows-0, 1:pncols-1]  # Shifted up by 1
shift4 = padded_img[1:pnrows-1, 2:pncols-0]  # Shifted left by 1
shift5 = padded_img[1:pnrows-1, 0:pncols-2]  # Shifted right by 1
shift6 = padded_img[2:pnrows-0, 2:pncols-0]  # Shifted up by 1, left by 1
shift7 = padded_img[2:pnrows-0, 0:pncols-2]  # Shifted up by 1, right by 1
shift8 = padded_img[0:pnrows-2, 2:pncols-0]  # Shifted down by 1, left by 1
shift9 = padded_img[0:pnrows-2, 0:pncols-2]  # Shifted down by 1, right by 1

shifts = [shift1, shift2, shift3, shift4, shift5, shift6, shift7, shift8, shift9]
stacked_shifts = np.stack(shifts, axis=-1)

median_filtered_img = np.median(stacked_shifts, axis=-1)

t1 = time()

print('Execution time:', t1-t0)

Execution time: 0.13090229034423828

We see that this is way faster than our previous approach!

Note: In this implementation, we manually wrote out each of the nine shifted arrays that were used. Can we do this in a more algorithmic way? That is, can we have Python automatically generate these shifts? One idea is to again use a pair of nested for loops, but this time we iterate through lists of horizontal and vertical shifts (instead of rows and columns).

# shift9 = padded_img[0:pnrows-2, 0:pncols-2]  # Shifted down by 1, right by 1
# shift2 = padded_img[0:pnrows-2, 1:pncols-1]  # Shifted down by 1
# shift8 = padded_img[0:pnrows-2, 2:pncols-0]  # Shifted down by 1, left by 1

# shift5 = padded_img[1:pnrows-1, 0:pncols-2]  # Shifted right by 1
# shift1 = padded_img[1:pnrows-1, 1:pncols-1]  # Unshifted noisy_img
# shift4 = padded_img[1:pnrows-1, 2:pncols-0]  # Shifted left by 1

# shift7 = padded_img[2:pnrows-0, 0:pncols-2]  # Shifted up by 1, right by 1
# shift3 = padded_img[2:pnrows-0, 1:pncols-1]  # Shifted up by 1
# shift6 = padded_img[2:pnrows-0, 2:pncols-0]  # Shifted up by 1, left by 1

shifts = []
for row_shift in [-1, 0, 1]:
    for col_shift in [-1, 0, 1]:
        shifts.append(padded_img[row_shift + 1:row_shift + pnrows-1, col_shift + 1:col_shift + pnrows-1])

stacked_shifts = np.stack(shifts, axis=-1)

t0 = time()

nrows, ncols = noisy_img.shape
padded_img = np.zeros((nrows + 2, ncols + 2))
padded_img[1:-1, 1:-1] = noisy_img

pnrows, pncols = padded_img.shape

shifts = []
for row_shift in [-1, 0, 1]:
    for col_shift in [-1, 0, 1]:
        shifts.append(padded_img[row_shift + 1:row_shift + pnrows-1, col_shift + 1:col_shift + pnrows-1])

stacked_shifts = np.stack(shifts, axis=-1)

median_filtered_img = np.mean(stacked_shifts, axis=-1)

t1 = time()

print('Execution time:', t1-t0)

Execution time: 0.05407381057739258

Note: This implementation is hard-coded to use 3 by 3 grid sizes. Can we modify it to work with any odd-sized grid?