In the first five notebooks of this workshop session, we've reviewed the fundamentals tools of image processing: numerical arrays, convolutions, point clouds and gradient descent.

Now, for those of you who went *really* fast through these first steps, here's a little bonus: an introduction to the JPEG (1992) and JPEG2000 compression standards, which are ubiquitous in cameras and cinemas. These two formats rely on the **Fourier** and **Wavelet transforms**, two mathematical tools that can be thought of as **precursors to Convolutional Neural Networks**.

**References, going further.**
These two bonus notebooks are based on the Numerical Tour on approximation with orthogonal bases, written by Gabriel PeyrÃ©. For additional references, you may have a look at:

- The chapter 3 (Harmonic analysis) of my course of
*Mathematics for non-mathematicians*. - Gabriel's handbook + my workshop notes, written for maths students at the ENS.

First, we re-import our libraries:

In [1]:

```
%matplotlib inline
import center_images # Center our images
import matplotlib.pyplot as plt # Display library
import numpy as np # Numerical computations
from imageio import imread # Load .png and .jpg images
```

Re-define custom display routines:

In [2]:

```
def display(im): # Define a new Python routine
"""
Displays an image using the methods of the 'matplotlib' library.
"""
plt.figure(figsize=(8,8)) # Square blackboard
plt.imshow( im, cmap="gray", vmin=0, vmax=1) # Display 'im' using a gray colormap,
# from 0 (black) to 1 (white)
def display_2(im_1, title_1, im_2, title_2):
"""
Displays two images side by side; typically, an image and its Fourier transform.
"""
plt.figure(figsize=(12,6)) # Rectangular blackboard
plt.subplot(1,2,1) ; plt.title(title_1) # 1x2 waffle plot, 1st cell
plt.imshow(im_1, cmap="gray") # Auto-equalization
plt.subplot(1,2,2) ; plt.title(title_2) # 1x2 waffle plot, 2nd cell
plt.imshow(im_2, cmap="gray", vmin=-7, vmax=15)
```

And load, once again, a simple axial slice:

In [3]:

```
I = imread("data/aortes/1.jpg", as_gray=True) # Import as a grayscale array
I = I / 255 # Normalize intensities in the [0,1] range
I = I[::2,::2] # Subsample the image, for convenience
display(I) # Let's display our slice:
```

Compression algorithms rely on **transforms** `f`

, which turn an image `I`

into
a new array `f(I)`

that is supposed to be easier to handle.
The most fundamental of these "helpers" is the Fourier Transform (click, it's great!), which decomposes a signal or an image as a **superposition of harmonics** (just like a piano note, really), with weights encoded in the array

In [4]:

```
# The numpy packages provides a Fast Fourier Transform in 2D,
# and its inverse (the iFFT). FFTshift and iFFTshift
# are just there to get nicer, centered plots:
from numpy.fft import fft2, ifft2, fftshift, ifftshift
fI = fft2(I) # Compute the Fourier transform of our slice
# Display the logarithm of the amplitutde of Fourier coefficients.
# The "fftshift" routine allows us to put the zero frequency in
# the middle of the spectrum, thus centering the right plot as expected.
display_2( I, "Image",
fftshift( np.log(1e-7 + abs(fI)) ), "Fourier Transform" )
```

To get an intuition of this new object, the simplest thing to do is to take our Fourier Transform `fI`

, edit it, and see what the image `ifft2( edit( fI ))`

looks like:

In [5]:

```
def Fourier_bandpass(fI, fmin, fmax) :
"""
Truncates a Fourier Transform fI, before reconstructing a bandpassed image.
"""
Y, X = np.mgrid[:fI.shape[0], :fI.shape[1]] # Horizontal and vertical gradients
radius = (X - fI.shape[0]/2) ** 2 \
+ (Y - fI.shape[1]/2) ** 2 # Squared distance to the middle point
radius = ifftshift( np.sqrt(radius) ) # Reshape to be fft-compatible
fI_band = fI.copy() # Create a copy of the Fourier transform
fI_band[ radius <=fmin ] = 0 # Remove all the low frequencies
fI_band[ radius > fmax ] = 0 # Remove all the high frequencies
I_band = np.real(ifft2(fI_band)) # Invert the new transform...
display_2( I_band, "Image", # And display
fftshift( np.log(1e-7 + abs(fI_band)) ), "Fourier Transform" )
```

As evidenced below, **Fourier coefficients that are close to the center encode the low frequencies of the signal**:

In [6]:

```
Fourier_bandpass(fI, 0, 10)
```

As we add more coefficients, we see that **details** start to appear:

In [7]:

```
Fourier_bandpass(fI, 0, 50)
```

We can also keep specific frequencies and compute "detail-only" images at a cheap numerical cost. **Convolutions**, that we presented in the 2nd notebook, can all be implemented this way:

In [8]:

```
Fourier_bandpass(fI, 50, 100)
```

Our image can then be simply expressed as a sum of low, medium and high frequencies:

In [9]:

```
Fourier_bandpass(fI, 0, 100) # = Sum of the last two images
```