This post is the second of the Fourier-transform for time series, check the first here:
Fourier transform for time-series: fast convolution explained with numpy
Quick review of the previous post
In the first post, I explained how the Fourier-transform can be used to convolve signals very efficiently. I showed that Convolution using the Fourier-transform in numpy is many orders of magnitude faster that the standard algebraic approach, and that it corresponds to a certain type of convolution called circular convolution.
In this post, I want to emphasize what the circular convolution means and how it all applies to images. Images are also a good way to extend the 1-dimension intuition into 2 dimensions.
All images were made by the author.
Image convolution with scipy
If you’ve ever worked with images for image processing, you most likely have encountered functions to apply convolution. Convoluting images is used everywhere – image enhancement, denoising, segmentation, feature extraction, compression – and is at the base of Convolutionnal Neural Networks, the gold standard of deep learning model to process visual data.
In Python, image convolution can be done quite simply using Scipy and its ndimage subpackage. At this point, I recommend taking a quick look at the documentation of the convolve
function, and then come back here.
The use is very simple: you can pass two images to convolve them together. Let’s see an example:
Note that scipy proposes several ways to handle the boundaries using the parameter ‘mode‘: as we will see below, the mode ‘wrap‘ corresponds to circular convolution and hence to convolution using a Fourier-transform approach. Other approaches exist, like ‘reflect‘ that reflects the images inside-out, or ‘constant‘ that repeats the outermost value. Notice also how ‘wrap‘ works: it repeats the whole signal, as if it was periodic.
Convolution of 2D images
Let’s start coding to see the differences between different convolution modes.
First, we create a class to represent 2D periodic images: remember from the previous post that when using Fourier-transform tool, the signal are considered to be periodic. This class is just syntactic sugar to plot such 2d periodic arrays.
We show the "base" image in the [0, V, 0, H] rectangle, as well as its 8 first replicas around. As stated in the previous post, the signal is considered periodic hence with infinite support, but we only need and use a single period.
Let’s now create a sample image to play with: it shall contain random noise, a sinusoidal pattern, a slope pattern, and a few square spots. We also create the periodic version of this sample image: it represents the periodic image that the Fourier-transform considers when applying its operators:
Let’s now create a kernel to use for the convolution: we’ll use a simple constant kernel, also called averaging kernel since the convolution with this kernel just gives the local average of the input image.
We then start playing with scipy convolution function and its differents modes to handle the boundaries, and wrap the result as a periodic array for easy plotting: notice how the middle of the convoluted image is always the same whatever the mode used, but the boundaries vary.
Now we can use a Fourier-transform approach to compute the convolution: as shown in the previous post, we just need to take the inverse Fourier-transform of the product of the Fourier-transform of both signals, the image and the kernel:
Comparing the result with the "wrap" mode of scipy, we can see that the results look a lot alike, just with a slight shift:
This is just a matter of indexing, and we can get the exact same results using a shifted-centered kernel:
Using proper centering, we then got identical results between scipy‘s convolution with mode=’wrap‘, and theFourier-transform approach.
Out of curiosity, let’s see which approach is faster:
Again, the Fourier-transform approach was faster, and in this case faster than a scipy function, which is nice.
Wrap up
We have seen in this post how the circular convolution translates to images, and how it is equivalent to scipy convolution function using mode=’wrap’.
In the next post, we’ll dive in the use of window functions in the context of Fourier-transform to reduce spectral leakage and improve spectral analysis.
Subscribe to get future posts about Fourier transform directly onto your feed!
Also, check out my other posts and if you like any of them, please subscribe it helps me a lot to reach my goal of 100 subscribers:
Fourier transform for time-series: fast convolution explained with numpy
PCA/LDA/ICA: a components analysis algorithms comparison
PCA-whitening vs ZCA-whitening: a numpy 2d visual
300-times faster resolution of Finite-Difference Method using numpy
If those articles seem interesting to you, remember to follow me, the new articles will appear on your feed.