Monday, March 31st

Image denoising

Last week, we developed some code to apply the mean and median filter using 3 by 3 grids and ignoring the left/right/top/bottom edges of the image:

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

We read in a sample image and converted to grayscale by averaging the RGB values:

In [2]:
hamburg = plt.imread('hamburg.png')

hamburg_gray = hamburg.mean(axis=2)
plt.imshow(hamburg_gray,cmap='gray',vmin=0,vmax=1)
Out[2]:
<matplotlib.image.AxesImage at 0x1951102e610>

We copied the hamburg_gray array before adding noise. We used Boolean arrays to identify which pixels will be hit with salt noise and which will be hit with pepper noise.

In [3]:
noisy_hamburg = hamburg_gray.copy()
num_rows, num_cols = hamburg_gray.shape
random_array = np.random.random((num_rows, num_cols))

# Let's generate an Boolean array with True in the locations that we want
# to add salt noise and False everywhere else
salt_mask = (random_array > .9)

# Then, we can use `salt_mask` as a slice on the `noisy_hamburg` array:
noisy_hamburg[salt_mask] = 1

# Do the same thing for pepper noise:
pepper_mask = (random_array < .1)
noisy_hamburg[pepper_mask] = 0

Note: The cell above is a useful skeleton for writing the sp_noise function described in the project page. It requires some modification to take inputs img and noise, where img is the image array that will have noise added to it, and noise measures the total proportion of pixels that will be hit by salt or pepper noise.

In [4]:
plt.imshow(noisy_hamburg, cmap='gray',vmin=0,vmax=1)
Out[4]:
<matplotlib.image.AxesImage at 0x195111be580>

We then applied the mean filter using 3 by 3 grids, ignoring the left/right/top/bottom edges of the image:

In [5]:
filtered_hamburg = noisy_hamburg.copy()

# For now, we need to skip the first row and last row
for i in range(1,num_rows-1):
    # We also need to skip the first column and last column
    for j in range(1, num_cols-1):
        grid = noisy_hamburg[i-1:i+2, j-1:j+2]
        mean = np.mean(grid)

        filtered_hamburg[i,j] = mean
In [6]:
plt.imshow(filtered_hamburg,cmap='gray',vmin=0,vmax=1)
plt.title('Mean filter')
Out[6]:
Text(0.5, 1.0, 'Mean filter')

We also applied the median filter:

In [7]:
filtered_hamburg = noisy_hamburg.copy()

# For now, we need to skip the first row and last row
for i in range(1,num_rows-1):
    # We also need to skip the first column and last column
    for j in range(1, num_cols-1):
        grid = noisy_hamburg[i-1:i+2, j-1:j+2]
        median = np.median(grid)
        
        filtered_hamburg[i,j] = median
In [8]:
plt.imshow(filtered_hamburg,cmap='gray',vmin=0,vmax=1)
plt.title('Median filter')
Out[8]:
Text(0.5, 1.0, 'Median filter')

Some updates:

Let's write a mean_filter function based on the mean filter code above:

In [9]:
def mean_filter(noisy_img):
    filtered_img = noisy_img.copy()
    num_rows, num_cols = noisy_img.shape

    # For now, we need to skip the first row and last row
    for i in range(1,num_rows-1):
        # We also need to skip the first column and last column
        for j in range(1, num_cols-1):
            grid = noisy_img[i-1:i+2, j-1:j+2]
            mean = np.mean(grid)

            filtered_img[i,j] = mean
            
    return filtered_img

Can we write a median_filter function as well?

In [12]:
def median_filter(noisy_img):
    filtered_img = noisy_img.copy()
    num_rows, num_cols = noisy_img.shape

    # For now, we need to skip the first row and last row
    for i in range(1,num_rows-1):
        # We also need to skip the first column and last column
        for j in range(1, num_cols-1):
            grid = noisy_img[i-1:i+2, j-1:j+2]
            median = np.median(grid)

            filtered_img[i,j] = median
            
    return filtered_img
In [13]:
filtered_hamburg = median_filter(noisy_hamburg)
In [14]:
plt.imshow(filtered_hamburg, cmap='gray', vmin=0, vmax=1)
Out[14]:
<matplotlib.image.AxesImage at 0x19512ea57f0>

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

The filter_func input should be a function that takes in an 2D array and returns a float.

In [16]:
def img_filter(noisy_img, filter_func=np.mean):
    filtered_img = noisy_img.copy()
    num_rows, num_cols = noisy_img.shape

    # For now, we need to skip the first row and last row
    for i in range(1,num_rows-1):
        # We also need to skip the first column and last column
        for j in range(1, num_cols-1):
            grid = noisy_img[i-1:i+2, j-1:j+2]
            filtered_img[i,j] = filter_func(grid)
            
    return filtered_img
In [17]:
mean_filtered_hamburg = img_filter(noisy_hamburg, np.mean)
median_filtered_hamburg = img_filter(noisy_hamburg, np.median)
In [20]:
plt.figure(figsize=(10,10))

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

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

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

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

plt.tight_layout()

Edge pixels

Until now, we've simply ignored filtering of edge pixels. Let's try to deal with them now.

The project page discusses adding extra rows/columns to our array so that we are able to construct a 3 by 3 grid centered at all pixels of our original image.

If we are using 3 by 3 grids, we need to add one extra row/column on all sides.

In [22]:
num_rows, num_cols = noisy_hamburg.shape

# Define a larger array with an extra on top/bottom
# and an extra column on left/right
padded_img = np.zeros((num_rows+2, num_cols+2))

# Place the noisy array in the center of the padded array
padded_img[1:-1, 1:-1] = noisy_hamburg

Now, we want to apply our mean filter to this padded_img array. We need to be careful about our indexing:

In [25]:
filtered_img = np.zeros((num_rows, num_cols))

for i in range(num_rows):
    for j in range(num_cols):
        grid = padded_img[i:i+3, j:j+3]
        filtered_img[i,j] = np.mean(grid)
In [28]:
plt.figure(figsize=(10,10))
plt.imshow(filtered_img, cmap='gray',vmin=0,vmax=1)
plt.axis('off')
Out[28]:
(-0.5, 899.5, 599.5, -0.5)

Note: We have added black border to the padded_img (since we used np.zeros), which propogates into the filtered_img when we take any mean that includes padded pixels.

We could instead add a white border (by using np.ones), but that's no better off. We could split the difference by using np.ones(...)/2 to use a gray border.

Note: To implement this strategy of dealing with edge pixels, we will need update our mean_filter/median_filter/img_filter functions. It will helpful to first a get_padded_img function that takes in an array and returns a padded version.

Grid size

What needs to change if we want to use grids of different sizes, say s=3, s=5, s=7, ...?

Let's try 5 by 5 grids:

In [34]:
filtered_hamburg = noisy_hamburg.copy()

# For now, we need to skip the first row and last row
for i in range(2,num_rows-2):
    # We also need to skip the first column and last column
    for j in range(2, num_cols-2):
        grid = noisy_hamburg[i-2:i+3, j-2:j+3]
        mean = np.mean(grid)

        filtered_hamburg[i,j] = mean
In [35]:
plt.figure(figsize=(10,10))
plt.imshow(filtered_hamburg, cmap='gray',vmin=0,vmax=1)
plt.axis('off')
Out[35]:
(-0.5, 899.5, 599.5, -0.5)

What about 7 by 7 grids?

In [36]:
filtered_hamburg = noisy_hamburg.copy()

# For now, we need to skip the first row and last row
for i in range(3,num_rows-3):
    # We also need to skip the first column and last column
    for j in range(3, num_cols-3):
        grid = noisy_hamburg[i-3:i+4, j-3:j+4]
        mean = np.mean(grid)

        filtered_hamburg[i,j] = mean
In [37]:
plt.figure(figsize=(10,10))
plt.imshow(filtered_hamburg, cmap='gray',vmin=0,vmax=1)
plt.axis('off')
Out[37]:
(-0.5, 899.5, 599.5, -0.5)

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

Exercise: Once you've done this, combine it with the earlier discussion on filtered edge pixels.

Wednesday, April 2nd

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

Let's look at some other considerations when dealing with image denoising (not things that are necessary for the project).

What about "salt and pepper" noise with color images?

In [7]:
color_img = plt.imread('Pilotwings.jpg') / 255 # Convert the integer RGB data to floating point

plt.imshow(color_img)
Out[7]:
<matplotlib.image.AxesImage at 0x13dc105f5e0>

Let's add some salt and pepper noise:

In [11]:
noisy_img = color_img.copy()

num_rows, num_cols = noisy_img.shape[:2]

random_array = np.random.random((num_rows, num_cols))

salt_mask = random_array > .9
pepper_mask = random_array < .1

noisy_img[salt_mask] = 1
noisy_img[pepper_mask] = 0
In [10]:
plt.imshow(noisy_img)
Out[10]:
<matplotlib.image.AxesImage at 0x13dc06e8190>

Let's try to apply our median filter:

In [12]:
filtered_img = noisy_img.copy()

# For now, we need to skip the first row and last row
for i in range(1,num_rows-1):
    # We also need to skip the first column and last column
    for j in range(1, num_cols-1):
        grid = noisy_img[i-1:i+2, j-1:j+2]
        median = np.median(grid)
        
        filtered_img[i,j] = median
In [13]:
plt.imshow(filtered_img)
Out[13]:
<matplotlib.image.AxesImage at 0x13dc06e0550>

With our old code: for each pixel to be filtered, we've constructed a grid consisting of 3 rows, 3 columns, and 3 color channels (RGB). We then the np.median of this 3 by 3 by 3 slice, and set the RGB value for the filtered pixel to this median.

It might be better to separately compute the median red value, median green value, and median blue value.

In [15]:
#help(np.median)

We can use the keyword argument axis=[0,1] inside np.median to only compute medians throught the rows and columns, but the colors:

In [16]:
filtered_img = noisy_img.copy()

# For now, we need to skip the first row and last row
for i in range(1,num_rows-1):
    # We also need to skip the first column and last column
    for j in range(1, num_cols-1):
        grid = noisy_img[i-1:i+2, j-1:j+2]
        median = np.median(grid, axis=[0,1])
        
        filtered_img[i,j] = median
In [17]:
plt.imshow(filtered_img)
Out[17]:
<matplotlib.image.AxesImage at 0x13dc1099eb0>

In the code above, we've hit all color channels together with salt/pepper noise (if they get hit at all). It could also be the case that the color channels are independently affected by "salt"/"pepper" noise.

In [18]:
noisy_img = color_img.copy()

num_rows, num_cols = noisy_img.shape[:2]

random_array = np.random.random((num_rows, num_cols,3))

salt_mask = random_array > .9
pepper_mask = random_array < .1

noisy_img[salt_mask] = 1
noisy_img[pepper_mask] = 0
In [19]:
plt.imshow(noisy_img)
Out[19]:
<matplotlib.image.AxesImage at 0x13dc0741280>
In [23]:
filtered_img = noisy_img.copy()

# For now, we need to skip the first row and last row
for i in range(1,num_rows-1):
    # We also need to skip the first column and last column
    for j in range(1, num_cols-1):
        grid = noisy_img[i-1:i+2, j-1:j+2]
        median = np.median(grid, axis=[0,1])
        
        filtered_img[i,j] = median
In [21]:
plt.imshow(filtered_img)
Out[21]:
<matplotlib.image.AxesImage at 0x13dc11af7f0>

Speed

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

For simplicity, let's return to working with grayscale images:

In [30]:
img = (plt.imread('Pilotwings.jpg') / 255).mean(axis=2)
In [31]:
noisy_img = img.copy()

num_rows, num_cols = noisy_img.shape

random_array = np.random.random((num_rows, num_cols))

salt_mask = random_array > .9
pepper_mask = random_array < .1

noisy_img[salt_mask] = 1
noisy_img[pepper_mask] = 0
In [32]:
plt.imshow(noisy_img, cmap='gray', vmin=0, vmax=1)
Out[32]:
<matplotlib.image.AxesImage at 0x13dc31113d0>
In [34]:
from time import time
In [44]:
t0 = time()

filtered_img = noisy_img.copy()

# For now, we need to skip the first row and last row
for i in range(1,num_rows-1):
    # We also need to skip the first column and last column
    for j in range(1, num_cols-1):
        grid = noisy_img[i-1:i+2, j-1:j+2]
        median = np.median(grid)
        
        filtered_img[i,j] = median
        
t1 = time()

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

Execution time: 9.303221464157104

Idea: Let's construct nine versions of our array, in the following manner:

  1. The original array
  2. Shift the original array down by one pixel
  3. Shift the original array up by one pixel
  4. Shift left by one
  5. Shift right by one
  6. Shift up by one, left by one
  7. Shift up by one, right by one
  8. Shift down by one, left by one
  9. Shift down by one, right by one

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

In [36]:
num_rows, num_cols = noisy_img.shape

shift1 = noisy_img.copy()

# Shift down by one
shift2 = np.zeros((num_rows, num_cols))
shift2[1:, :] = noisy_img[:-1, :]

# Shift up by one
shift3 = np.zeros((num_rows, num_cols))
shift3[:-1, :] = noisy_img[1:, :]

# Shift left by one
shift4 = np.zeros((num_rows, num_cols))
shift4[:, :-1] = noisy_img[:, 1:]

# Shift right by one
shift5 = np.zeros((num_rows, num_cols))
shift4[:, 1:] = noisy_img[:, :-1]

# Shift up by one, left by one
shift6 = np.zeros((num_rows, num_cols))
shift6[:-1,:-1] = noisy_img[1:, 1:]

# Shift up by one, right by one
shift7 = np.zeros((num_rows, num_cols))
shift7[:-1,1:] = noisy_img[1:, :-1]

# Shift down by one, left by one
shift8 = np.zeros((num_rows, num_cols))
shift8[1:,:-1] = noisy_img[:-1, 1:]

# Shift down by one, right by one
shift9 = np.zeros((num_rows, num_cols))
shift9[1:,1:] = noisy_img[:-1, :-1]

The np.stack function can stack several arrays together into a new dimension:

In [39]:
stacked_img = np.stack([shift1, shift2, shift3, 
                        shift4, shift5, shift6, 
                        shift7, shift8, shift9])

print(stacked_img.shape)

(9, 662, 868)

Now compute median across axis 0:

In [41]:
filtered_img = np.median(stacked_img, axis=0)
In [42]:
plt.imshow(filtered_img, cmap='gray',vmin=0, vmax=1)
Out[42]:
<matplotlib.image.AxesImage at 0x13dc1395940>

Let's go back and time this new code.

In [55]:
t0 = time()

num_rows, num_cols = noisy_img.shape

shift1 = np.zeros((num_rows, num_cols))
shift1[:,:] = noisy_img[:,:]

# Shift down by one
shift2 = np.zeros((num_rows, num_cols))
shift2[1:, :] = noisy_img[:-1, :]

# Shift up by one
shift3 = np.zeros((num_rows, num_cols))
shift3[:-1, :] = noisy_img[1:, :]

# Shift left by one
shift4 = np.zeros((num_rows, num_cols))
shift4[:, :-1] = noisy_img[:, 1:]

# Shift right by one
shift5 = np.zeros((num_rows, num_cols))
shift5[:, 1:] = noisy_img[:, :-1]

# Shift up by one, left by one
shift6 = np.zeros((num_rows, num_cols))
shift6[:-1,:-1] = noisy_img[1:, 1:]

# Shift up by one, right by one
shift7 = np.zeros((num_rows, num_cols))
shift7[:-1,1:] = noisy_img[1:, :-1]

# Shift down by one, left by one
shift8 = np.zeros((num_rows, num_cols))
shift8[1:,:-1] = noisy_img[:-1, 1:]

# Shift down by one, right by one
shift9 = np.zeros((num_rows, num_cols))
shift9[1:,1:] = noisy_img[:-1, :-1]

stacked_img = np.stack([shift1, shift2, shift3, 
                        shift4, shift5, shift6, 
                        shift7, shift8, shift9])

filtered_img2 = np.median(stacked_img, axis=0)

t1 = time()

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

Execution time: 0.0906379222869873

This new version is around ~100 times faster.

Can we clean up our code a little bit?

In [ ]:
num_rows, num_cols = noisy_img.shape


for row_shift in [-1,0,1]:
    for column_shift in [-1,0,1]:
        shift = np.zeros((num_rows,num_cols))
        shift[row_shift:, :] = noisy_img[:-row_shift]
        
        
        # Needs some thought after class...