Hello, and welcome to
Rush to deadline: AP 186 edition
and we’ll be discussing about the
Fourier Transform Model of Image Formation
Part 1. The Fourier transform
Fourier transforms are one of the most useful transforms known to man. Snapchat has this mathematical beauty to appreciate, else they’d have to find another way to modify filters. For 2-D, it follows the general form
It is an integral transform that extracts a function in x to be a function in inverse space (or 1/x). Along with photo apps, it is extensively used in signal processing, image processing, and physics (as if we don’t already know that).
For pretty much the rest of the images to be used here, we’ll be going with 128×128 resolution to be kind to our fan-favorite fft2 function. Contrary to usual situations, smaller resolutions are actually better here to better visualize the pixel-large delta spikes for the sinusoids to be produced later in this activity.
Part 2. Familiarization with the 2-D discrete fast Fourier transform
circle = imread('AP186_circle.png');
fft_circle = fft2(double(circle)); //fast fourier transform. converts matrix values from integer to double
We start with a circle of radius 0.1. Load it up as imread, then we apply the function fft2 on the double. As with the manual, fft2 works best for pixel resolutions with base 2. An application of the double function is used to ensure that the resulting FFT is a 2-D matrix, and not just a single list (which happens too often as I am forgetful).
A3 = mat2gray(abs(fft_circle)); //mat2gray solves this: //imshow(uint8(255*abs(fftshift(fft_circle))/max(abs(fftshift(fft_circle))))) //absolute value of the fft of the circle imwrite(A3,'5A3.png');
I feel that I need some explaining to do for mat2gray, as I will happily use it for future activities dealing with FFT and images. Mat2gray is basically your one-stop-fix-your-problems-shorten-your-code function that normalizes the matrix and then multiplies it by 255, the maximum brightness value for an unsigned 8 bit integer (which mat2gray also converts).
Long story short, it converts a matrix to a grayscale image.
Additionally, we need to get the absolute value since the FFT operation of a real function or matrix yields a complex matrix. As such, we only need the modulus of the value.
From the manual, we know that the quadrants are diagonally reflected, the resulting image will be presented the wrong way.
A4 = mat2gray(abs(fftshift(fft_circle))); imwrite(A4,'5A4.png');
Much better
With the fftshift (2-D) mode, here we see the familiar Airy pattern when the Fourier transform. It looks like a ripple of water radiating from the center point. This will be important for the convolution parts later. Remember:
A5 = mat2gray(abs(fft2(fft_circle))); //FFT it back
imwrite(A5,'5A5.png');
Finally, we apply the Fourier Transform once more to generate the “same” image once again. Due to the symmetry of the image along the horizontal center line, it is hard to determine whether a double application of fft2 will really provide the same result.
Now, let’s move on to an assymetric figure.
Part 2b. Letter A
I used paint to generate this baby.
(I’ll refrain from putting the code here as it is just similar to 2a).
The next three pictures show the same method done in 2a, but for the letter A (would’ve done the character more justice if it were part 2a).
Looking at the FFT-shifted photograph, we see that a cross has been formed for the letter A, and also signifying that the FFT was centered on the middle thanks to shifting. An additional point to note is that the letter A is asymmetrical to both the horizontal and vertical center lines. As such, the shifted FFT is also expected to be asymmetric along both axes. My guess is that a symmetrical character along the vertical would show symmetry.
After applying another FFT to the FFT of the character, we see that the letter A is inverting.
Anarchy.
I like it as it shows that I am not just copying the original photo. For people obsessing over resolving matters, change fft2 to the second application to ifft.
fft2(fft_matrix); //inverts A ifft(fft_matrix); //does not invert A
Part 2c. Sinusoid along x (corrugated roof)
Similar to 2a and 2b, I will refrain from posting similar codes for 2c.
However, there is a nuance to be addressed when plotting functions with negative values as greyscale images. Take the standard sine function for example
B = sin(X*25) //typical sine function
It would still look like a sine function, with the exception that the bright lines are thinner than the dark lines. For a trigonometric function with no vertical biasing, this should not be the case.
B_vertical_bias = sin(X*25) + bayas
The reason for this is that in converting a grayscale image from a 2-D matrix, the negative values are treated as zeros. But if we simply just add 1 to make the minimum value of the sine function to 0, the maximum value would be 2. This is another red flag.
To avoid these flags, we perform vertical shrinking and biasing via
B = sin(X*25)*0.5 + 0.5 //with vertical shrinking and biasing
With these changes, I would be able to make a sine matrix with values ranging from 0 to 1.
The resulting unshifted FFT looks like a completely black picture with the exception of a few incomplete dots on each corner of the image. With the fftshift, we see three dots in the center vertical line, with the center dot being the brightest. Analytically, this means that the center dot represents the vertical biasing performed a while ago, and the two dots equidistant from the center represent the frequencies of the image. As expected, the FFT of a sine function is a 
Part 2d. Simulated double slit
Next part deals with the double slit aperture.
The double slit was generated via this code. Thankfully, no biasing needed:
C = zeros(nx,ny)
C(find(X<0.05 & X>0.025)) = 1; //vary to change the width and separation
C(find(X<-0.025 & X>-0.05)) = 1;
If we want thicker lines in the original space, we just adjust the parameters wherein the elements in the zero matrix is 1. The FFT results are quite interesting.
No one loves you, unfftshifted photograph.
For the FT shifted photograph, we see a familiar fringe pattern found in the famous double slit experiment. Here I found the realization that
The Fourier transform of the double slit aperture is the interference pattern.
For this shape, I decided to use inverse FT just because.
Part 2e. Square function
Next section, we deal with the square aperture. After obtaining the fftshifted image, we obtain an interesting crosshair-like pattern. Projecting the image to either the center horizontal or vertical line, we can see the sinc function. This sinc function is actually the FT of the rectangle function, which is basically the 1-D version of the quadrilateral. The FT looks symmetric due to the fact that the square is symmetric both vertically and horizontally.
Part 2f. 2-D Gaussian bell curve
For the next part, we deal with the peculiar Gaussian function. Peculiar in the sense that
…the FT of a Gaussian function, is… well…, another Gaussian function.
Again, nothing too interesting for the unfftshifted photograph, but look at when it’s shifted. Everything becomes clear.
…it’s another Gaussian. How surprising! This verifies that what I generated originally is actually a Gaussian function.
And here we see that returning it back to its original space leads to another Gaussian function.Very reminiscient of old memes
Part 3. Simulation of an imaging device
We now approach the convolution section of this activity. I am tasked to convolve the photo of Ma’am Jing’s glorious lab’s initials to a circular aperture.
These are the two photos to be convolved. VIP image as the object, and the circle as the aperture. We already take the circle to be in frequency space, so we just fftshift it as seen in the code below.
vip = imread('5B_vip.bmp'); //already grayscaled for my convenience
ape = imread('5B_circle.bmp');
F_ape = fftshift(double(ape)); //aperture already in fourier plane
F_vip = fft2(double(vip));
Short note: remember that the FT of the circle is the airy pattern. So in the original space, we are convolving the image to this Airy pattern. In frequency space, we are multiplying the Fourier transform of the VIP image to the circle.
We use the following code for the convolution process itself:
FAP = F_ape.*(F_vip);
IAP = ifft(FAP); //Inverse FT
FImage = abs(IAP);
FImage1 = mat2gray(FImage);
imwrite(FImage1,'5B_final.png');
Here we see that the Airy pattern has been casted t the VIP lettering, thus completing the convolution. Now it looks like it’s the letters that are forming the ripples.
Part 4. Template matching using correlation
Correlations are usually needed to determine repetitive patterns, such as repeating letters in certain statements, such as the one below
One of the dominant letters in this sentence is the letter ‘A’ (seems to be a recurring theme this blog). We then turn our attention to the code:
spain = imread('5C_text.png');
arial_a = imread('5C_a.png');
g_spain = rgb2gray(spain);
g_arial_a = rgb2gray(arial_a);
fft_spain = fftshift(fft2(double(g_spain)));
fft_a = fftshift(fft2(double(g_arial_a)));
product = fft_a.*conj(fft_spain); //Correlation deals with conjugate multiplication
final_5C = abs(fftshift(fft2(product)));
final_5C1 = mat2gray(final_5C)
imwrite(final_5C1,'5C_final.png');
If the reader would kindly notice the sharp spikes of intensity in the right photograph, then he/she/it/they would notice that these spikes are coincident to the locations where the capital letter A was placed. As such, the ctrl+F function of the correlation function was put to test with these photographs.
.
Part 5. Edge detection using the convolution integral
We’re not going to do away with this image.
For the last part, we will trace the edges of the VIP letters via these 3 x 3 matrices. It is important to note that the elements of these matrices are zero-sum, and produce edge elements similar to their shapes.
Code:
//Part 5D
p_hori = [-1 -1 -1; 2 2 2; -1 -1 -1];
p_vert = [-1 2 -1; -1 2 -1; -1 2 -1];
p_spot = [-1 -1 -1; -1 8 -1; -1 -1 -1];
p_mist = [-1 1 1; -1 1 1; -1 1 1]
I deliberately
imwrite(conv2(double(vip),p_hori),'5D_horizontal.png');
The convolution here is horizontal, and the edge elements could be seen as horizontal as well. The vertical edges cannot be seen.
imwrite(conv2(double(vip),p_vert),'5D_vertical.png');
With here, the edge elements are now vertical. The horizontal edges cannot be seen, sadly.
imwrite(conv2(double(vip),p_spot),'5D_spot.png');
We use a central spot for this convolution. It seems to be the best edging convolution so far as it encapsulates the whole lettering.
imwrite(conv2(double(vip),p_mist),'5D_mistake.png')
I added a non-zero sum matrix just to see what it does. If the sum isn’t zero, it serves to not really change the image much. I would want to investigate why is this so given that I spend more time on this eye-opening activity.
Part 6. Evaluation
Peformance-wise: 8/10. I could have done a lot more with the activity but given my self-inflicted time constraint, my performance was clipped. I could have provided a better explanation for the correlation, or the edge-tracing capabilities of convolutions to zero sum matrices. Despite that, I believe that I was able to perform the activity well.
References
- Soriano, M. Fourier Transform Model of Image Formation. Applied Physics 186 Activity Manual.
- Sinc and rect function photograph taken from https://www.chem.purdue.edu/courses/chm621/text/ft/basiset/rect/rectangle.gif.
- Quicklatex.com for the matrix images.
- Carlo Solibet, Robbie Esperanza, Roland Romero, and Dr. Maricor Soriano for their invaluable assistance in helping me finish the activity.
Apol's blog 