Image Filtering¶
The main objectives of this module are:
- Implement point filtering with look-up tables.
- Implement region filtering with kernel convolution, morphological operations and edge detection.
- Understand spatial frequency information with the Fourier transform.
1. Look-up tables¶
Look-up tables are a very useful tool to perform simple pixel-level filtering of an image.
Numpy arrays indexing (see reference here) is a very powerful tool, but can be a bit counter-intuitive at first. We will use it to quickly create and use look-up tables to transform an image.
Look at the short example below.
We construct an image with 4 possible values for each pixel (-> $I(x,y) \in \{0,1,2,3\}$). We then build the following LUT:
| Input value | Output value |
|---|---|
| 0 | 2 |
| 1 | 3 |
| 2 | 1 |
| 3 | 0 |
And we apply it on the image.
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
%matplotlib inline
# Let's create a 3x3 "image" with 4 possible pixel values (in [0,3])
im = np.array([[1,0,2],[2,1,3],[3,3,1]])
# look-up table is a vector of size 4 mapping the new value for the pixels of m
lut = np.array([2,3,1,0]) # maps 0 -> 2, 1 -> 3, 2 -> 1, 3 -> 0
plt.figure(figsize=(8,2))
plt.plot(lut)
plt.title('Look-up table')
plt.show()
# we can use numpy vector indexing to apply the look-up table:
new_im = lut[im] # This will create a new array with the same size as m
plt.figure(figsize=(8,4))
plt.subplot(1,2,1) #supblot is useful to display images side-by-side
plt.imshow(im, cmap=cm.gray)
plt.title('Original')
plt.subplot(1,2,2)
plt.imshow(new_im, cmap=cm.gray)
plt.title('LUT(Original)')
plt.show()
Starting from this example, create the following LUTs and apply them to the "walking.jpg" image:
- For inverting an image (so that $0 \rightarrow 255$, $255 \rightarrow 0$)
- To keep only graylevel such that $g \ge t$, setting all other values to 0.
- To reduce the number of gray levels from 256 to 8 (so that $255 \rightarrow 7$, $0 \rightarrow 0$).
- To "stretch" the histogram so that, given a minimum value $T_{min}$ and a maximum value $T_{max}$, we have :
- If $I(x,y) < T_{min}$, the new value is set to 0
- If $I(x,y) > T_{max}$, the new value is set to 255
- Values between $T_{min}$ and $T_{max}$ are stretched to cover the entire histogram.
- to perform an equalization of the cameraman image (so that the histogram becomes "as flat as possible")
from skimage.io import imread,imshow,imsave
im = imread('walking.jpg')
print(im.shape, im.dtype)
plt.figure(figsize=(10,10))
imshow(im)
plt.show()
## -- Your code here -- ##
(799, 640) uint8
Need more help? You can check the following videos:
2. Kernel convolution¶
Kernal convolution allows us to filter an image based on the values of a neighborhood.
Write a program that applies a 3x3 kernel convolution on an image. Use it on the "walking" image to perform a mean filter.
def conv(im, kernel):
pass # TO DO
im = imread('walking.jpg')
## -- Your code here -- ##
Using the median filter and the mean filter from scikit-image, compare the behaviour of the two filters with neighborhood of increasing sizes on the "noisy" astronaut image:
from skimage.filters.rank import mean, median
im = imread('astronaut_noisy.jpg')
imshow(im)
## -- Your code here -- ##
<matplotlib.image.AxesImage at 0x7ff12cc685d0>
Need more help? You can check the following videos:
3. 2D Fourier transform¶
The 2D Fourier transform allows us to get a representation of the spatial frequencies present in an image. A very powerful way of filtering images is to modify those frequencies directly by modifying the "Fourier image", and to use the inverse transform to get the pixel values of the filtered image.
Let's illustrate that with a fake example. We start with a completely random signal and compute the Fourier transform. Note that the Fourier image is complex, so we use its amplitude. As the range of amplitudes is so high that we don't see much in the resulting image, we display the log of the amplitude.
from numpy.fft import fft2,ifft2,fftshift,ifftshift
# fft2 -> Fourier transform
# ifft2 -> Inverse transform
# fftshift & ifftshift -> reorganize the "Fourier image" to make it more easily interpretable.
im = np.random.random((512,512))
f = fftshift(fft2(im))
amplitude = np.sqrt(np.real(f)**2+np.imag(f)**2)
plt.figure(figsize=[10,5])
plt.subplot(1,2,1)
plt.imshow(im,cmap = plt.cm.gray)
plt.title('$f(x,y)$')
plt.subplot(1,2,2)
plt.imshow(np.log(amplitude))# show log so as to see more than just the global maximum
plt.title('$log(sqrt(|F(u,v)|^2))$')
plt.show()
We could now, for instance, decide to add an horizontal frequency by creating peaks in the Fourier image on the horizontal axis. The further away from the center we put the peaks, the higher the frequency.
f2 = f.copy()
f2[250:262,262:272] *= 1000
f2[250:262,240:250] *= 1000
amplitude = np.sqrt(np.real(f2)**2+np.imag(f2)**2)
im2 = ifft2(ifftshift(f2)).real
plt.figure(figsize=[10,5])
plt.subplot(1,2,1)
plt.imshow(im2,cmap = plt.cm.gray)
plt.title('$f(x,y)$')
plt.subplot(1,2,2)
plt.imshow(np.log(amplitude))# show log so as to see more than just the global maximum
plt.title('$log(sqrt(|F(u,v)|^2))$')
plt.show()
Now we could also decide to mask the central region, which would remove the peaks that we just created and put all those values at zero. The image will once again look random:
f3 = f2.copy()
f3[250:262,240:272] = 0
amplitude = np.sqrt(np.real(f3)**2+np.imag(f3)**2)
im3 = ifft2(ifftshift(f3)).real
plt.figure(figsize=[10,5])
plt.subplot(1,2,1)
plt.imshow(im3,cmap = plt.cm.gray)
plt.title('$f(x,y)$')
plt.subplot(1,2,2)
plt.imshow(np.log(amplitude, where=amplitude>0))
plt.title('$log(sqrt(|F(u,v)|^2))$')
plt.show()
Starting from those examples and the code below:
- Use the Fourier transform to reduce the dithering of the moire.png image by building a low-pass filter.
- Build a high-pass filter using Fourier transform and apply the filter to the "walking" image.
im = imread('moire1.png').astype(np.float)
f = fftshift(fft2(im)) # shift Fourier image so that the center corresponds to low frequencies
amplitude = np.sqrt(np.real(f)**2+np.imag(f)**2)
plt.figure(figsize=[10,5])
plt.subplot(1,2,1)
plt.imshow(im,cmap = plt.cm.gray)
plt.title('$f(x,y)$')
plt.subplot(1,2,2)
plt.imshow(np.log(amplitude))# show log so as to see more than just the global maximum
plt.title('$log(|F(u,v)|^2)$')
plt.show()
## -- Your code here -- ##
Need more help? You can check the following videos:
4. Morphological operations¶
- For morphology functions, see the skimage documentation
Starting from the example below:
- Using the "opening" operation with a disk structuring element of increasing size, determine how many circles of any given radius there are in the circles.png image.
- Build a morphological filter that eliminates one size of circles.
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from skimage.io import imread
from skimage.morphology import disk, erosion, dilation
im = imread('circles1.png')==0 #to be sure objects are = 1
plt.imshow(im,cmap=plt.cm.gray);
# use the local maximum and local minimum for dilation and erosion
eroded_image = erosion(im,disk(5))
dilated_image = dilation(im,disk(5))
plt.figure(figsize=[10,8])
plt.subplot(1,2,1)
plt.imshow(eroded_image)
plt.title('erosion')
plt.subplot(1,2,2)
plt.imshow(dilated_image)
plt.title('dilation')
# note that "erosion" and "dilation" are the same as the "local minimum" and "local maximum" :
from skimage.filters.rank import minimum,maximum
eroded_image = minimum(im.astype('uint8'), disk(5))
dilated_image = maximum(im.astype('uint8'), disk(5))
plt.figure(figsize=[10,8])
plt.subplot(1,2,1)
plt.imshow(eroded_image)
plt.title('erosion')
plt.subplot(1,2,2)
plt.imshow(dilated_image)
plt.title('dilation')
Text(0.5, 1.0, 'dilation')
Need more help? You can check the following videos:
5. Edge detection¶
Using the convolve2d function, develop the Sobel filter. The Sobel operator is given by $$\mathbf{G} = \sqrt{ {\mathbf{G}_x}^2 + {\mathbf{G}_y}^2 }$$ where $G_x$ and $G_y$ are images respectively obtained by a convolution with the given kernels $$\mathbf{g}_y = \begin{bmatrix} -1 & -2 & -1 \\ 0 & 0 & 0 \\ +1 & +2 & +1 \end{bmatrix} \quad \mbox{and} \quad \mathbf{g}_x = \begin{bmatrix} -1 & 0 & +1 \\ -2 & 0 & +2 \\ -1 & 0 & +1 \end{bmatrix} $$
Apply the Sobel filter to find the edges in the "road" image.
Compare with the results of the Canny edge filter with different "sigma" values (see the skimage documentation)
from skimage.data import camera
from skimage.io import imread
im = imread('road.jpg', as_gray=True)
plt.figure(figsize=(15,15))
plt.imshow(im, cmap=plt.cm.gray)
plt.show()
## -- Your code here -- ##
Use the Hough transform to find the main straight lines in the road image, using the results from the Canny edge detector.
## -- Your code here -- ##
Need more help? You can check the following videos:
Coding project - Picture enhancement¶
Write code that automatically enhances a photograph.
"Enhance" may mean a lot of different things, and we encourage you to be creative in which enhancements you want to implement.
Some possibilities include (not an exhaustive list):
- Noise reduction
- Auto-level
- Gamma correction (with gamma provided by the user or automatically determined from the image histogram)
- Increase colour saturation
- ...
# -- Your code here -- #