Monday, October 28th, 2026¶

Last class, we continued to work on creating tartan pattern images.

In [1]:
import numpy as np
import matplotlib.pyplot as plt

Project 4: Image denoising¶

In recent weeks, we've been working with images using (m,n,3)-shaped NumPy RGB arrays and visualizing them with plt.imshow. The cell below will download an example image file from the course webpage into your Jupyter notebook directory with the file name noisy_img.png.

In [2]:
import requests

with open('noisy_img.png','wb') as f:
    data = requests.get(r'https://jllottes.github.io/_images/noisy_img.png').content
    f.write(data)

Let's read this image in as a NumPy array using plt.imread and plot it using plt.imshow:

In [3]:
noisy_img = plt.imread('noisy_img.png')
plt.imshow(noisy_img)
Out[3]:
<matplotlib.image.AxesImage at 0x168e525b770>
No description has been provided for this image

Notice that the image looks to have been affected by some kind of noise. In particular, some pixels seem to have mistakenly been colored white, while others have been colored black. Let's take a closer look to see this more clearly.

In [4]:
zoom_noisy_img = noisy_img[:200, 200:400]
In [5]:
plt.imshow(zoom_noisy_img)
Out[5]:
<matplotlib.image.AxesImage at 0x168e6354a50>
No description has been provided for this image

This type of noise is called salt and pepper noise, where some pixels have been replaced with white (salt) pixels and other pixels have been replaced with black (pepper) pixels. For our next project, we are interested in exploring strategies to try to filter and remove this type of noise.

Loading in sample images¶

In order to explore filtering salt and pepper noise, it's useful to first be able to generate this kind of noise ourselves. That is, it will be helpful to have some Python code that can take an image and randomly assign some pixels to be salt noise and others to be pepper noise.

There are several example images available on the project page that can be used for this exploration. Pick whichever you like, download it, and place it into your Jupyter notebook directory. Then, read it into a NumPy array using plt.imread.

For example, I've downloaded circuit.png into my Jupyter notebook directory. Now, I will read it into a NumPy array:

In [7]:
circuit = plt.imread('circuit.png')

Let's use plt.imshow to see what this image looks like.

In [8]:
plt.imshow(circuit)
Out[8]:
<matplotlib.image.AxesImage at 0x168e7286850>
No description has been provided for this image

Before proceeding, let's take a look at the structure of the img array. For example, does it contain RGB values or RGBA values (with transparency)? Does it contain floats between 0 and 1 or does it contain integers between 0 and 255?

In [9]:
print(circuit.shape)   # Check the resolution and color channels

(426, 640, 3)

We see that the resolution is 640 by 426 and there are three color channels (i.e. RGB, not RGBA values).

In [11]:
print(circuit.dtype)
print(circuit[0,0])

float32
[0.12156863 0.12156863 0.12156863]

We see that the circuit.png image contains floating point RGB values (without transparency). If your chosen image has a different format, transform it so that it also contains floating point RGB values. As a reminder:

  • img[:,:,:3] will produce a slice consisting of only the RGB parts of the image and will drop any transparency information, and
  • img / 255 produce a floating point RGB array from an integer RGB array.

Grayscale images¶

For this project, we will be working with grayscale images. That is, each image contains only shades of gray.

Note: In an RGB image, a given pixel will be some shade of gray when the red, green, and blue values match. Let's confirm that our image file contains only gray pixels.

In [12]:
red_channel = circuit[:, :, 0]
green_channel = circuit[:, :, 1]
blue_channel = circuit[:, :, 2]
In [17]:
print(np.min(red_channel == green_channel))

True

Since the minimum of red_channel == green_channel is True, it must have True in every position. That is, the red and green channels match.

In [19]:
print(np.min(red_channel == blue_channel))

True

The red channel also matches the blue channel, so all three RGB values match, and the image is a grayscale image.

We can represent grayscale images as (m,n) arrays (rather than (m,n,3)) without any loss of information, since all three color channels give the same (m,n) array. In the case of the img array, we could take img[:,:,0] (the red channel) to represent the grayscale values.

In [20]:
gray_circuit = circuit[:,:,0]   # Take the red channel for grayscale values

Let's confirm that there's no loss of information:

In [21]:
plt.imshow(gray_circuit)
Out[21]:
<matplotlib.image.AxesImage at 0x168e73dc910>
No description has been provided for this image

Recall: when plt.imshow is passed a grayscale array, it defaults to the viridis colormap. Let's override this choice to use the gray colormap. It would also be helpful to specify vmin=0 and vmax=1 (or vmax=255 if working with integer data) to ensure that plt.imshow is correctly assigning colors to the values in the array.

In [23]:
plt.imshow(gray_circuit, cmap='gray', vmin=0, vmax=1)
Out[23]:
<matplotlib.image.AxesImage at 0x168e7db02d0>
No description has been provided for this image

What could we do if we started from a color image that we want to convert to a grayscale image, like mario.png? Note: If you did not previous download this image from the course page, you can run the cell below to do so.

In [24]:
with open('noisy_img.png','wb') as f:
    data = requests.get(r'https://jllottes.github.io/_static/downloads/mario.png').content
    f.write(data)
In [25]:
mario = plt.imread('mario.png')
plt.imshow(mario)
Out[25]:
<matplotlib.image.AxesImage at 0x168e803fed0>
No description has been provided for this image
In [37]:
print(mario.shape)

(224, 256, 4)

This array includes an alpha channel, so let's remove it.

In [38]:
mario = mario[:,:,:3]

The above image does not contain purely grayscale data, since the red/green/blue values do not match for at least some of the pixels. Can we construct a reasonable grayscale version of this image?

In [39]:
red_channel = mario[:, :, 0]
plt.imshow(red_channel, cmap='gray', vmin=0, vmax=1)
Out[39]:
<matplotlib.image.AxesImage at 0x168e8373ed0>
No description has been provided for this image
In [40]:
green_channel = mario[:, :, 1]
plt.imshow(green_channel, cmap='gray', vmin=0, vmax=1)
Out[40]:
<matplotlib.image.AxesImage at 0x168e840f110>
No description has been provided for this image
In [41]:
blue_channel = mario[:, :, 2]
plt.imshow(blue_channel, cmap='gray', vmin=0, vmax=1)
Out[41]:
<matplotlib.image.AxesImage at 0x168e848a490>
No description has been provided for this image

The np.mean function can be used to compute averages of an array.

In [42]:
np.mean(mario)
Out[42]:
np.float32(0.57220167)
In [43]:
print(mario[10, 20])  # RGB triple in row 10, column 20
print(np.mean(mario[10,20]))   # Mean of the RGB triple in row 10, column 20

[0.36078432 0.5803922  0.9882353 ]
0.6431373
In [ ]:
In [44]:
nrows, ncols = mario.shape[:2]

gray_mario = np.zeros((nrows, ncols))

for row in range(nrows):
    for col in range(ncols):
        gray_mario[row, col] = np.mean(mario[row, col])    # Store the mean of the the RGB triple
In [45]:
plt.imshow(gray_mario, cmap='gray', vmin=0, vmax=1)
Out[45]:
<matplotlib.image.AxesImage at 0x168e8511950>
No description has been provided for this image

For our purposes, we don't want to take the average value of every value in the mario array. Instead, we want to average the RGB triples separately for each pixel in the image.

We can include an optional argument axis to specify the dimension through which means should be calculated. In particular, including the argument axis=2 tells NumPy to compute the average through the color channel only (i.e. axis 2 of the (m,n,3) array). In other words, this will average the RGB triple separately for each pixel.

In [47]:
#help(np.mean)

In our previous previous calculation with nested for loops, we went pixel-by-pixel and separately computed the mean of the RGB triple for each pixel. In other words, we are computing the mean only through the color channel, which is axis 2 of our 3D array.

In [49]:
gray_mario = np.mean(mario, axis = 2)
In [50]:
plt.imshow(gray_mario, cmap='gray', vmin=0, vmax=1)
Out[50]:
<matplotlib.image.AxesImage at 0x168e858ce10>
No description has been provided for this image

Adding salt and pepper noise¶

Before adding noise to our image, let's make a copy of the array to preserve the original image array. Let's call this copy noisy_img (we will come back to add noise to this image later).

In [60]:
noisy_img = circuit.copy()
In [61]:
print(noisy_img.shape)

(426, 640, 3)
In [62]:
noisy_img = np.mean(noisy_img, axis = 2)   # Convert to a 2D grayscale array

We would like a method for randomly selecting pixels of the noisy_img array to which we can add salt/pepper noise. The np.random.random function will be useful for this task, which works as follows:

  • np.random.random() will return a randomly selected float in the interval $[0,1)$.
  • np.random.random(shape) will return an array with shape shape, where each entry contains a randomly selected float in the interval $[0,1)$.

Strategy: let's generate an array of the same shape as our image filled with randomly drawn values between 0 and 1. If the randomly drawn value is greater than, say .9, we'll add salt noise in that position.

In [53]:
np.random.random()
Out[53]:
0.033091847824884724
In [54]:
np.random.random((10,5))
Out[54]:
array([[0.87203521, 0.20884608, 0.2828869 , 0.6131427 , 0.08446156],
       [0.33186108, 0.60462635, 0.07329707, 0.75846063, 0.98428425],
       [0.84741076, 0.00695036, 0.81554276, 0.39958771, 0.26426598],
       [0.12184005, 0.35203429, 0.46816805, 0.80483077, 0.34337356],
       [0.24568637, 0.61783906, 0.47797448, 0.30085885, 0.72488346],
       [0.85804667, 0.75388118, 0.49241381, 0.69991157, 0.79295514],
       [0.42242437, 0.44042547, 0.87137992, 0.9174853 , 0.52198081],
       [0.11223625, 0.05886283, 0.7613832 , 0.47763146, 0.44576682],
       [0.01611859, 0.69527523, 0.03983186, 0.22832442, 0.74305388],
       [0.25466495, 0.65366248, 0.63906875, 0.36768453, 0.18897103]])
In [64]:
nrows, ncols = noisy_img.shape

for row in range(nrows):
    for col in range(ncols):
        if np.random.random() > .9:
            noisy_img[row, col] = 1
In [65]:
plt.imshow(noisy_img, cmap='gray', vmin=0, vmax=1)
Out[65]:
<matplotlib.image.AxesImage at 0x168e8d10f50>
No description has been provided for this image

What if we want to add pepper noise as well?

In [66]:
nrows, ncols = noisy_img.shape

for row in range(nrows):
    for col in range(ncols):
        r = np.random.random()
        if np.random.random() > .9:
            noisy_img[row, col] = 1   # Add salt noise
        if np.random.random() < .1:
            noisy_img[row, col] = 0   # Add pepper noise
In [67]:
plt.imshow(noisy_img, cmap='gray', vmin=0, vmax=1)
Out[67]:
<matplotlib.image.AxesImage at 0x168e91d8e10>
No description has been provided for this image
In [ ]:

Note: the implementation above gives a preference to pepper noise, since it's possible for salt pixels to be overwritten as pepper (but not the reverse). Instead, let's use a single random array, but select select different subintervals between 0 and 1 to identify salt/pepper noise:

In [74]:
noisy_img = np.mean(circuit.copy(), axis=2)   # Initialize noisy_img as a 2D grayscale array
nrows, ncols = noisy_img.shape

for row in range(nrows):
    for col in range(ncols):
        r = np.random.random()
        if r > .9:
            noisy_img[row, col] = 1   # Add salt noise
        if r < .1:
            noisy_img[row, col] = 0   # Add pepper noise
In [75]:
plt.imshow(noisy_img, cmap='gray', vmin=0, vmax=1)
Out[75]:
<matplotlib.image.AxesImage at 0x168e90abed0>
No description has been provided for this image

Note: In the example code above, we iterated through every row and column and decided for each entry whether or not to add salt and pepper noise. This is perfect case to use Boolean arrays.

Idea: create a salt_mask Boolean array containing True in each entry that we wish to hit with salt noise. That is, salt_mask should contain True everywhere that the random_array is less than, say, 0.1.

In [79]:
np.random.random((4,6)) < .5
Out[79]:
array([[False, False, False,  True,  True, False],
       [ True, False,  True,  True, False, False],
       [ True, False, False, False, False, False],
       [ True,  True, False, False, False,  True]])
In [81]:
noisy_img = np.mean(circuit.copy(), axis=2)   # Initialize noisy_img as a 2D grayscale array
nrows, ncols = noisy_img.shape

r = np.random.random((nrows, ncols))  # Array of the same shape as noisy_img full of random floats
                                      # between 0 and 1.

salt_mask =   r > .9 # Identify 10% of entries that will be hit with salt noise
pepper_mask = r < .1 # Identify 10% of entries that will be hit with pepper noise

noisy_img[salt_mask]   = 1 # Add salt noise
noisy_img[pepper_mask] = 0 # Add pepper noise
In [82]:
plt.imshow(noisy_img, cmap='gray', vmin=0, vmax=1)
Out[82]:
<matplotlib.image.AxesImage at 0x168e912bb10>
No description has been provided for this image

Exercise. Use the above sample code to write a function sp_noise(img, noise) that adds salt and pepper noise to an image. Its first argument img should be a 2-dimensional numpy array representing the image and the second argument noise should be the fraction of pixels that are to be replaced by noise (for example, with noise = 0.05 about 5% of pixels should be noise, consisting in roughly equal parts of white and black pixels). The function should return a 2-dimensional numpy array representing the original image with noise added.

In [ ]:
In [ ]:
In [ ]:

A description of the mean filter is provided on the project page. Let's try to apply the mean filter.

  • For right now, we will use a 3 by 3 grid centered centered at each pixel to compute the mean.
  • For right now, we will ignore all pixels at the edges of the image.

For now, I will use the pre-noised image file that was downloaded from the course webpage. If you have a function sp_noise function, you should generate your own noisy image to use.

In [ ]:
In [ ]:
In [ ]:

Again, let's make a copy filtered_img of the noisy array that we will filter so that we do not modify the noisy array. This way, we will be able to compare how our filtered image looks with how the noisy image looked.

In [ ]:
In [ ]:

If we are considering the pixel in the ith row and jth column (i.e. noisy_img[i,j], we want a 3 by 3 slice centered on the [i,j]th pixel. That is, we want a slice that includes the i-1, i, and i+1 rows and the j-1, j, and j+1 columns.

In [ ]:
In [ ]:
In [ ]:

Now let's go pixel by pixel and take the mean of each 3 by 3 grid to create our filtered image array. As a reminder, for now we will ignore any problematic pixels along the edge of the image.

In [ ]:
In [ ]:
In [ ]:

What about the median filter? The np.median function will compute the median of an array.

In [ ]:
In [ ]:
In [ ]:

Exercise: Use the code above to write a function simple_mean_filter that applies the mean filter to non-edge pixels using a 3 by 3 grid.

In [ ]:
In [ ]:
In [ ]:

Exercise: Modify the simple_mean_filter function to write a simple_median_filter that instead applies the median filter.

In [ ]:
In [ ]:
In [ ]:

Note: We could combine these two functions into a single simple_filter function that takes in a noisy image array (noisy_img) along with some function (filter_func) that will be applied to each grid.

The input filter_func should be a Python function that takes in an 2D array and returns a float (e.g. np.mean or np.median).

In [ ]:
In [ ]:
In [ ]: