Thursday, July 29, 2010

AP 186 Act8 - Enhancement in the Frequency Domain

This activity is mostly about applying what we've learned in the last two activities about Fourier Transforms (FT). We'll be applying the techniques to images, making them of better quality. :D

Part A. Convolution Theorem
For this part we're first asked to make certain patterns and take their FT's.

Two dots along x-axis symmetric about the center:

pattern (left) and its FT (right).

Two circles along x-axis symmetric about the center:
Increasing the radius...

Note that by just adding a small finite radius to the dots, the FT immediately changes such that it is now a combination of an Airy disk with a sinusoid running along x-axis. This is because the FT of two peaks symmetric along the x-axis is a sinusoid while the FT of a circle is the Airy disk. The pattern is a convolution of the two peaks with the circle. Therefore, the FT modulus seen is a product of the sinusoid with the Airy disk. Furthermore, as we increase the radius of the dots, the Airy disk becomes smaller and smaller.

Two squares along x-axis symmetric about the center:

pattern (left) and its FT (right).

Increasing the width...

This pattern is again a convolution of two dirac delta's or peaks with a square aperture. Therefore, the FT modulus should be the product of a sinusoid with the cross-like rectangular patterns (FT of square aperture) and this is what we got. Increasing the width of the squares would make the rectangular pattern smaller (a property of FT's).

Two gaussians along x-axis symmetric about the center:
pattern (left) and its FT (right)

Increasing the variance...

This pattern is a convolution of two dirac delta's and a gaussian. The peaks are spaced futher apart therefore the lines are much closer together. The FT of a gaussian is still a gaussian so the FT modulus is the product of the sinusoid and the gaussian as seen in the figure. Increasing the variance of the gaussian leads to a smaller gaussian in the FT modulus.

Next, we're asked to create a 200 x 200 array with ten 1's randomly located. Also, we made a 3x3 pattern of

[0 0 1]
0 1 0
1 0 0]

and convolved the two. The result is:

(left) array with random 1's. (right) convolved image

The dots are now replaced with the 3x3 pattern.

Then, we create a 200 x 200 array with equally spaced 1's along the x and y axis. We vary the spacing of the 1's and then see what happens in the FT modulus:

grid pattern (left) and its FT (right). width = 5

grid pattern (left) and its FT (right). width = 20

grid pattern (left) and its FT (right). width = 40

It is observed that the closer the 1's are spaced, the farther apart they are in frequency space, which is consisted with the inverse property of FT.


Part B. Fingerprints: Ridge Enhancement
For this part we'll be enhancing the ridges of the fingerprints by filtering in frequency space.

I'm going to use this fingerprint that I got from the manual:

fingerprint

The fingerprint reminds me of sine waves that are propagating in all direction (2D) and it's like a ripple. So I thought: If you rotate the two peaks (FT of sinusoid) in all directions then you would get an annulus.

FT of fingerprint (left). Filtering mask (middle). Product of mask and FT (right).

It's FT is the left figure and information is contained mostly in the bright annulus. So I designed the filter to let this part pass through as well as the center spot. Removing the center spot will invert the image (black --> white and white --> black) since 0's and 1's are interpreted as black and white respectively in scilab. Multiplying this filter or mask to the FT of the fingerprint and then inverse FT-ing results in:

original picture (left). Filtered fingerprint picture (right).


The black blotches that were present in the original (left) are now greatly reduced in the frequency filtered image (right). Putting a threshold or converting this to binary, we can see the effect of filtering more clearly:

original fingerprint (left). filtered fingerprint (right).

We can see that the blotches are removed and the ridges are enhanced. :)


Part C. Lunar Landing Scanned Pictures: Line Removal
Basically we're asked to remove the vertical lines in the picture caused by piecing together individual pictures into a composite one.

lunar picture (left) and its FT (right).


It's FT is shown on the right and since the vertical lines are equally spaced, we can think of this as a sinusoid. Therefore, the frequency components that make up these vertical lines belong to the horizontal bright line in the FT. Therefore I designed my filter to be:


and applying this, I get:


and the vertical lines are gone! The left one is the original while the right figure is the filtered one. :)


Part D. Canvas weave modeling and Removal
For this part, we're asked to remove the canvas weave patterns by filtering in the frequency domain.

Oil painting from UP Vargas Museum Collection.


FT of the painting (left) and with a filtering mask applied to it (right).

The FT of the picture is shown in the left figure. Note that there are bright horizontal and vertical lines as well as some bright peaks. The canvas weaves are like grids so they can be thought of sine waves in the horizontal and vertical so I designed my mask to "cover" them as shown in the right. And after applying the filter, what I got is:

Original painting (left) and filtered image of painting (right).

Voila! The canvas weave patterns are gone! :D :D Also, the brushstrokes in the painting seems to have been enhanced as compared to the original.

Lastly, we invert the filtering mask and then take it's inverse FT. What I get is:

Inverted mask (left) and its inverse FT (right).


Cool pattern! Not only that, it looks like the canvas weave pattern in some way because it has grids too. :D


Grading myself, I give a 10/10 for understanding the lesson and producing the required outputs (whew!)
Score: 10/10

Finally, I would like to acknowledge Dr. Soriano, Tisza Trono, Arvin Mabilangan, BA Racoma and Joseph Bunao for the helpful discussions. :D

References:
1. M. Soriano, "A8 - Enhancement in the Frequency Domain"

Appendix: (Code) (There are missing semi-colons and slash because i removed them. This blog sometimes just won't display the codes correctly.)

nx = 256;
ny = 256;
x = linspace(-128,128,nx);
y = linspace(-128,128,ny);
[X,Y] = ndgrid(x,y);

// two dots
A = ones(nx, ny);
A(find(Y>20)) = 0;
A(find(Y<19))>1)) = 0;
A(find(X<-1)) = 0; two circles r= sqrt(X.^2 + Y.^2); A = zeros (nx,ny); A(find((X.^2 + (Y+20).^2) < a =" ones(nx,">25)) = 0;
A(find(Y<15))>5)) = 0;
A(find(X<-5)) = 0; two gaussians nx = 256; ny = 256; x = linspace(-1,1,nx); y = linspace(-1,1,ny); [X,Y] = ndgrid(x,y); r= sqrt((X.^2 + (Y+0.4).^2)/0.0001); A1 = exp(-r); r= sqrt((X.^2 + (Y-0.4).^2)/0.0001); A2 = exp(-r); A = A1 + A2; grid A = zeros(200, 200); width = 5; for i = 1:width:200 for j = 1: width:200 A(i, j) = 1; end end For Part B: nx = 419; ny = 581; x = linspace(-1,1,nx); y = linspace(-1,1,ny); [X,Y] = ndgrid(x,y); r= sqrt(X.^2 + Y.^2); finger = gray_imread("C:\fingerprint.jpg"); FTfinger = fftshift(fft2(finger)); mask = ones(nx, ny); mask = zeros(nx, ny); mask(find(r<0.20)) a =" FTfinger" b =" im2bw(finger," c =" im2bw(abs(ifft(A))," a =" gray_imread(" fta =" fftshift(fft2(A));" nx =" 480;" ny =" 640;" x =" linspace(-1,1,nx);" y =" linspace(-1,1,ny);" mask =" ones(nx,">0.10)) = 0;
mask(find(abs(X)>0.045)) = 1;

C = FTA .* mask;

scf();
imshow(log(abs(C)+0.1), []);

scf();
imshow(A, []);

scf();
imshow(abs(fft2(C)),[]);


For Part D:

nx = 419; ny = 581; //defines the number of elements along x and y

x = linspace(-1,1,nx); //defines the range
y = linspace(-1,1,ny);
[X,Y] = ndgrid(x,y);
r= sqrt(X.^2 + Y.^2);


mask1 = ones(nx, ny);
mask1(find(abs(Y)<0.02)) = 0;
mask1(find(abs(X)<0.04)) = 1;

mask2 = ones(nx, ny);
mask2(find(abs(X)<0.02)) = 0;
mask2(find(abs(Y)<0.04)) = 1;


mask = ones(nx, ny);
mask(find((Y.^2 + (X+0.43).^2) < 0.0004)) = 0;
mask(find((Y.^2 + (X-0.43).^2) < 0.0004)) = 0;
mask(find((Y.^2 + (X+0.205).^2) < 0.0004)) = 0;
mask(find((Y.^2 + (X-0.205).^2) < 0.0004)) = 0;

mask(find(((Y+0.31).^2 + X.^2) < 0.0004)) = 0;
mask(find(((Y-0.31).^2 + X.^2) < 0.0004)) = 0;
mask(find(((X+0.1).^2 + (Y-0.16).^2) < 0.0004)) = 0;
mask(find(((X-0.1).^2 + (Y-0.16).^2) < 0.0004)) = 0;
mask(find(((X+0.1).^2 + (Y+0.16).^2) < 0.0004)) = 0;
mask(find(((X-0.1).^2 + (Y+0.16).^2) < 0.0004)) = 0;


//scf();
//imshow(log(abs(FTA)+0.1),[]);

C = FTA .* mask2 .* mask1 .* mask;

scf();
imshow(log(abs(C)+0.1), []);

scf();
imshow(A, []);

scf();
imshow(abs(ifft(C)),[]);

Wednesday, July 21, 2010

AP 186 Act7 - Properties of the 2D Fourier Transform

For this activity, we'll explore some properties of the 2D Fourier Transform (FT):

Part A: Familiarization with FT of different 2D patterns
In this part, we created 2D patterns and then took their FT's and see what it looks like.

Square:

square pattern and its FT

Annulus (Donut):

Annulus pattern and its FT

Square annulus:
square annulus and its FT

Two slits along the x-axis symmetric about the center:

double slit pattern along x-axis and its FT

Two dots along x-axis symmetric about the center:

Two dots (1 pixel) along x-axis pattern and its FT

Note: If the dots are not of 1 pixel, a whole different FT pattern will appear:
different FT pattern if dots are not of 1 pixel


Part B: Anamorphic property of the Fourier Transform

First, we're asked to make a 2D sinusoid in the x direction (like a corrugated roof) and take it's FT:
sine wave (f=2Hz) and its FT

Next we change the frequency of the sinusoid:

f = 4 Hz

f = 8 Hz

f = 16 Hz

As the frequency of the sine wave is increased or decreased, the two peaks in the frequency space moves farther apart or closer together, respectively. This is because low frequencies are close to the center (zero) in frequency space and as you move away from the center, the frequency increases.

Adding a bias to the sinusoid and taking the FT:

sinusoid with bias and its FT

By adding the bias, the sinusoid now oscillates with non-negative values. However, we can see that there's an additional peak right in the middle of its FT. Digital images have no negative values so for example, what I can do to find the actual frequencies of the interferogram from Young's double slit experiment is to DC-filter it first. DC filtering will remove the bias and hence the middle peak in the FT. If the bias is non-constant, e.g. low frequency sinusoids, I can still get the frequencies of the interferogram by applying a high pass filter. This way, the low frequencies (the noise in this case), will be removed.

Next, we rotate the sinusoid and take its FT:

rotated 30 degrees

rotated 45 degrees

Rotating the sinusoid also amounts to rotating their FT by the same amount. This is a unique property of the 2D FT compared to the 1D FT.

Next, we created a pattern which is a combination of sinusoids in the X and Y direction, take its FT and see what it looks like:

combination of sinusoids in the X and Y and its FT

Then we add several rotated sinusoids in the pattern. I added two sinusoids: 1 rotated by 30 degrees and another one rotated by 45 degrees. Since I know that rotating the sinusoids will also rotate their FT, the FT's of the added sinusoids individually are shown in the last part. At first I predict that adding the sinusoids would result in a rotation of the FT of the "eggcart". But remembering that the FT is a linear transformation, the FT of the result should be the addition of the FT's of the components. Therefore, the resulting FT should be something like 4 peaks on the corners of a square, then two peaks tilted at 30 degrees and another 2 peaks tilted 45 degrees. Checking this prediction:


Indeed, the resulting FT is the same as the prediction! :D


Rating myself, I would give myself a 10/10 for understanding the lesson and producing the required outputs.

Score: 10/10

Lastly, I would like to acknowledge Dr. Soriano, BA Racoma and Androphil Polinar
for the helpful discussions.

- Dennis

References:
1. M. Soriano, "A7 - Properties of the 2D Fourier Transform"

Monday, July 19, 2010

AP 186 Act6 – Fourier Transform Model of Image Formation

In this activity we will be dealing with Fourier Transforms (Fast Fourier Transforms) and their applications in image processing. Cooley and Tukey made the FFT algorithm and it’s widely used in image as well as signal processing.

Part A. Familiarization with FFT

First we’re asked to make a white circle centered on a black background:

(Note that this circle can also represent a circular aperture. I’ve made my circle smaller since this will result in a bigger Airy disk as will be seen later.)

Then we converted this to grayscale and performed the Fast Fourier Transform (FFT). Abs() was used to get the intensity values since the FFT will return an array of complex numbers.

I = imread('C:\circle.PNG');
Igray = im2gray(I);
FIgray = fft2(Igray);
imshow(abs(FIgray), []);
And this is the result:

FFT of the white circle

It looks like a pure black picture but take a look at the four corners. :p

This is a result of the FFT. It’s fast, but it also results in the reversal of the quadrants diagonally. Therefore, applying the fftshift will give:

imshow(fftshift(abs(FIgray)), []);
shifted FFT of the circle

Now that’s better. :) That’s the Fourier transform of the white circle.

Applying the FFT again on the transformed image (or consequently, applying the FFT twice on the original image) will theoretically result in an inverted original image:

imshow(abs(fft2(FIgray)));

circle FFT'ed twice

Well, we got the original image alright. But we can’t tell if this is inverted since an inverted circle is still a circle. So we try this on the letter “A”:

The (unshifted) FFT of A:

I = imread('A.bmp');
Igray = im2gray(I);
FIgray = fft2(Igray);
imshow(abs(FIgray), []);

unshifted FFT of "A"

And applying fftshift:

imshow(fftshift(abs(FIgray)), []);

shifted FFT of "A"

And that’s the Fourier transform of “A”.

Now applying the FFT twice… we get:

imshow(abs(fft2(FIgray)));

"A" FFT'ed twice

Which is indeed, the inverted image of the letter “A”. :D An inverted image results in the application of FFT twice because doing so will result in the introduction of a negative sign in the x’s and y’s. If you want to get your original image back from a Fourier transformed image, do the inverse Fourier transform. :)

Part B. Simulation of an imaging device

We’re asked to create an image of the letters “VIP” and this will be our object to be imaged. I used Arial font because according to the activity manual, sans serif fonts are sharp.

object

We also created a white circle on a black background and this represents the aperture of a lens.

aperture

The image of the “VIP” can be obtained by convolving the object with the transfer function of the imaging system. Convolution is done by taking the product of the FFT of the grayscaled object with the fftshifted aperture and then inverse fourier transforming. The aperture image doesn’t need to be FFT-ed because it is already in the Fourier Plane.

T = gray_imread('C:\VIP.bmp');
A = gray_imread('C:\circle.PNG');
Ar = fftshift(A);
Ta = fft2(T);
FRA = Ta.*Ar;
IRA = fft2(FRA);
imshow(abs(IRA), []);
The output is:

r = 0.7

Note that it is inverted and this is true for real images.

Making the size of the aperture smaller, we get…

r = 0.4

r = 0.25

r = 0.1

The image quality reduces when the lens aperture becomes smaller and smaller because less and less rays are able to pass through it. The more rays the lens can collect, the better the image quality.

Part C. Template Matching using correlation

Template matching is basically a type of pattern recognition in which it will determine matches between the test object and the template.

For this part, we created the white texts in black background of: “THE RAIN IN SPAIN STAYS MAINLY IN THE PLAIN.”

text

and our template will be the letter “A” with the same font style and size as in our text:

template

We get the correlation of these two by multiplying the FT of A and the conjugate of the FT of the text and then inverse Fourier transforming by using ifft().

T = gray_imread('C:\text.PNG');
A = gray_imread('C:\A.bmp');
FT = fft2(T);
FA = fft2(A);
FTconj = conj(FT);
FRA = FA.*FTconj;
IRA = ifft(FRA);
imshow(abs(IRA), []);
This is what I got:

correlation of text and template using ifft

It looked like the quadrants were reversed diagonally so I applied the fftshift after inverse FT-ing:

imshow(fftshift(abs(IRA)), []);

shifted correlation

Now it looks inverted… so I tried using the fft2 instead of ifft:

IRA = fft2(FRA);
imshow(abs(IRA), []);
correlation of text with template using ifft2

and fftshifting…


imshow(fftshift(abs(IRA)), []);

shifted correlation

Now this looks more correct. The original text could be made out from this image but there seems to be bright spots were the letter “A” is present. Therefore the bright spots represent a match since these are the positions were the correlation is highest.

I thought… what if the template is a “P” and there are “B”’s in the text… what would happen? So I made my own text:

text 2

And my own template:

template 2

And the correlation is…

limitation of method

It can be seen that despite there is only one capital “P”, there are multiple bright spots. This is because the other letters like “B” looks very similar to P except that it has an extra part at the bottom. The other letters also have parts that are similar to “P” so they have a high correlation. Therefore this shows the limitation of the technique.


Part D. Edge detection using the convolution integral

We were asked to make a 3x3 matrix pattern of edges such that the sum of all the elements is zero. Then we convolved this with the "VIP" image using imcorrcoef():

For the matrix:

[-1 -1 -1

2 2 2

-1 -1 -1]

We get:

convolved image

We can see that the horizontal edges of the image are brighter. Also note that the matrix has the same values horizontally. What this does it that it scans the image and tries to match its edges with the matrix.

For the vertical matrix:
[-1 2 -1
-1 2 -1
-1 2 -1]

We get:

convolved image

This time, the vertical edges of the image are brighter :D

For the spot pattern:
[-1 -1 -1
-1 8 -1
-1 -1 -1] :

convolved image

This time, all of the edges of the image are bright. This is because in the pattern, all of the edges have the same value.

The code for this is:
T = gray_imread("C:\VIP.bmp");

pattern = [-1 -1 -1; 2 2 2; -1 -1 -1];

F = imcorrcoef(T, pattern);

imshow(F, []);

And that wraps up the activity. :)

For the score, I give myself a 10/10 for producing the required output and understanding the lesson. Plus 2 points for testing the limitations of the technique of correlation, which gives a total of...

Score: 12/10

Lastly, I would like to acknowledge Dr. Soriano for the helpful discussions. :D


- Dennis

References:
1. M. Soriano, "A6 - Fourier Transform Model of Image Formation"