In the previous notebook, we introduced the Fourier Transform that represents images as harmonics or waves of varying frequencies. Using a local (block-wise) Cosine Transform, we were able to mimick the JPEG standard and compress images... at a significant cost: the introduction of blocking artifacts.
To alleviate this problem, engineers and mathematicians from the 80's came up with a new idea: the Wavelet Transform, that we are now going to present. But first, we need to re-import our standard libraries:
%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
And re-define our
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=-1,vmax=1) # Auto-equalization
Once again, we work on a simple axial slice:
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:
What if "Fourier" coefficients referred to well-located wavelets instead of pure, oscillating waves than span across the full domain? After all, when we write down music, we tend to use notes that have a well-defined frequency and duration... The Fast Wavelet Transform, introduced by Stéphane Mallat in 1989, is all about computing such a "music score" for images through the iteration of convolutions and subsamplings.
In the picture below, which is adapted from one of Mallat's papers, the input image is displayed on the right. At each layer, we compute highpass features through convolutions with oriented edge detectors, and a lowpass thumbnail that is subsampled before being fed to the same procedure once again. We stop when the "sub-sub-...-sub-sampled" residual is reduced to a single pixel, and the "Wavelet Transform" of the input image is then defined as the concatenation of all the feature maps, computed at different resolutions:
For image compression, we typically use three "highpass" filters at each scale, thus ending up with a nice set of thumbnails that fit into each other like Russian dolls:
# We use Gabriel Peyre's implementation of 2D wavelet transforms: from nt_toolbox.general import * from nt_toolbox.signal import * from nt_toolbox.compute_wavelet_filter import * def wavelet_transform(I, Jmin=1, h = compute_wavelet_filter("Daubechies",6)) : """ 2D-Wavelet decomposition, using Mallat's algorithm. By default, the convolution filters are those proposed by Ingrid Daubechies in her landmark 1988-1992 papers. """ wI = perform_wavortho_transf(I, Jmin, + 1, h) return wI def iwavelet_transform(wI, Jmin=1, h = compute_wavelet_filter("Daubechies",6)) : """ Invert the Wavelet decomposition by rolling up the operations above. """ I = perform_wavortho_transf(wI, Jmin, - 1, h) return I wI = wavelet_transform(I) # Compute the wavelet transform plt.figure(figsize = (8,8)) # Square blackboard plot_wavelet(wI, 1) # Custom display routine, with a red tesselation plt.title('Wavelet coefficients - with added recursive grid') plt.show()
Multiscale pyramid. Here, the low frequencies sit in the upper-left corner and represent the image at coarser and coarser scales. To extract "bandpass" images as in the previous notebook, we simply have to zero specific bands in our Wavelet Transform
def Wavelet_bandpass(wI, depthMin, depthMax) : wI_band = wI.copy() # Create a copy of the Wavelet transform # Remove the low frequencies: wI_band[ :int(2**(depthMin-1)), :int(2**(depthMin-1))] = 0 wI_band[ 2**depthMax:, :] = 0 # Remove the high frequencies wI_band[ :, 2**depthMax:] = 0 # Remove the high frequencies I_band = iwavelet_transform(wI_band) display_2( I_band, "Image", # And display wI_band, "Wavelet Transform" )
If we just keep the coarser scales, we retrieve a "blurred" image:
Wavelet_bandpass(wI, 0, 4) # Scales 0, 1, 2, 3, 4
Wavelet decompositions produce features that are well-located:
Wavelet_bandpass(wI, 5, 6) # Scales 5 and 6
In sharp contrast to the Fourier bandpass images of the sixth notebook, these wavelet decompositions do not really suffer from ringing artifacts: oscillations remain close to the shapes' supports:
Wavelet_bandpass(wI, 7, 8) # Scales 7 and 8
Crucially, Wavelet Transforms perform very well in compression tasks at a cheap numerical cost. To illustrate this, we follow the same steps as in Part 6:
def Wavelet_threshold(wI, threshold) : # Re-implement a thresholding routine wI_thresh = wI.copy() # Create a copy of the Wavelet transform wI_thresh[ abs(wI) < threshold ] = 0 # Remove all the small coefficients I_thresh = iwavelet_transform(wI_thresh) # Invert the new transform... display_2( I_thresh, "Image", # And display wI_thresh, "Wavelet Transform" ) return I_thresh
Thresholding operations tend to preserve the sharpness of geometric features:
Wavelet_threshold(wI, .2) ;
If we compare our results to the 1:33 compressions of the previous notebook, the Wavelet Transform looks good!
Abs_values = np.sort( abs(wI.ravel()) ) # Sort the coefficients' magnitudes print(Abs_values) # in ascending order... cutoff = Abs_values[-2000] # And select the 2000th largest value JPEG2000 = Wavelet_threshold(wI, cutoff) # as a cutoff display(JPEG2000)
[ 0. 0. 0. ... 12.37511 12.566159 42.94246 ]
To remember: Multiscale image decompositions (pyramids), computed through a succession of convolutions and sub/supsamplings, have a long history.
Easy to compute (convolutions are cheap), these representations encode a multiscale prior: features at scale
k should be computed locally as aggregations of finer descriptors at scale
From the 70's to the turn of the XXIst century, a great amount of effort was put into the engineering of convolution filters with optimal properties. The Daubechies filters or the CDF wavelets (used in the actual JPEG2000 standard) are prime examples of this research, which powers digital cinemas or MRI reconstruction softwares.
Convolutional "Neural Networks". But since 2012, the combination of GPUs and high-level automatic differentiation libraries (Theano, TensorFlow, PyTorch...) has changed the rules of the game: as we'll see on Saturday, filter weights can now be tuned automatically using gradient descent, and don't have to be hand-picked anymore.
This major breakthrough has allowed engineers to best fit their models to the available data, and to explore new families of convolutional architectures that were previously "too hard to tune" to be useful. This is impressive research.
But don't be fooled by the hand-waving, pseudo-biological intuitions that are often bundled with these methods. As should be clear by now, the "Neural Networks" in Convolutional Neural Networks is little more than a cool, but very misleading piece of mathematical jargon:
Final words. When marketing people put a biological neuron in their slides, they sell you a dream, not a product: "AI" simply doesn't exist. You should always be cautious (but optimistic) when reading about progress in the field: