Activity 7 - Properties of the 2D Fourier Transform

In the first part of the activity, (a) different 2D patterns were created and their (b) fourier transforms were taken (see Figure 1-5).

(a)

(b)

Figure 1.A 128x128 px white square

(a)

(b)

Figure 2. A circular annulus with inner radius 32 px and outer radius 64 px

(a)

(b)

Figure 3. A square annulus with inner square side of 64 px and outer square side of 128 px.

(a)

(b)

Figure 4. Two horizontal slits 128 pixels apart.

(a)

(b)

Figure 5. Two dots horizontally aligned and 128 pixels apart.

For the analysis of the anamorphic property of the Fourier Transform, 2D sinusoids were created in Scilab. Figure 6 shows sinusoids with varying frequencies and their corresponding Fourier transforms. Also, images of the FT modulus are shown with a constant bias added to the original sinusoids. This addition of constant is for simulating real images, that is, digital images with no negative values. It has been observed that as the frequency increases, the distance between the peaks increases also. Also, when the bias was added, a peak at f=0 appeared.

Figure 6. Horizontal sinusoids [z = sin(2*%pi*f*X)] with frequency (a) f = 2, (b) f = 4,

(c) f = 8, (d) f = 16 and (e) f = 32.

Second and third rows display the FT modulus of the sinusoids and biased sinusoids.

In Figure 7 and Figure 8, rotating the sinusoid by 30 degrees resulted to the rotation of the line containing the two peak frequencies.

Figure 7. 30-degree rotated sinusoids [z = sin(2*%pi*f*(Y*sin(theta) + X*cos(theta)))] with frequency

(a) f = 2, (b) f = 4, (c) f = 8, (d) f = 16 and (e) f = 32.

Second and third rows display the FT modulus of the sinusoids and biased sinusoids.

Figure 8. 45-degree rotated sinusoids [z = sin(2*%pi*f*(Y*sin(theta) + X*cos(theta)))] with frequency

(a) f = 2, (b) f = 4, (c) f = 8, (d) f = 16 and (e) f = 32.

Second and third rows display the FT modulus of the sinusoids and biased sinusoids.

Figure 9 shows sinusoids that resulted from the product of two corrugated roofs, one running in the x-direction and the other in y-direction.

Figure 9. corrugated roofs or sinusoids [z = sin(2*%pi*f*X).*sin(2*%pi*f*Y)] with frequency

(a) f = 2, (b) f = 4, (c) f = 8, (d) f = 16 and (e) f = 32.

Second and third rows display the FT modulus of the sinusoids and biased sinusoids.

Figure 10 and 11 show what happened when rotated sinusoids were added to the sinusoids in Figure 9. Ans as expected, as the frequency increases, distance between peaks increases and also, a rotation of the line connecting the peaks of the rotated sinusoids resulted.

Figure 10. 30-degree rotated sinusoids added to sinusoids in Figure 9

[z = sin(2*%pi*f*X).*sin(2*%pi*f*Y) + sin(2*%pi*f*(Y*sin(theta) + X*cos(theta)))]

with frequency (a) f = 2, (b) f = 4, (c) f = 8, (d) f = 16 and (e) f = 32.

Second and third rows display the FT modulus of the sinusoids and biased sinusoids.

Figure 11. 45-degree rotated sinusoids added to sinusoids in Figure 9

[z = sin(2*%pi*f*X).*sin(2*%pi*f*Y) + sin(2*%pi*f*(Y*sin(theta) + X*cos(theta)))]

with frequency (a) f = 2, (b) f = 4, (c) f = 8, (d) f = 16 and (e) f = 32.

Second and third rows display the FT modulus of the sinusoids and biased sinusoids.

The following code was used to create the images above.

  
//Activity 7.B Anamorphic property of the Fourier Transform
nx = 100; ny = 100;
x = linspace(-1,1,nx);
y = linspace(-1,1,ny);
[X,Y] = ndgrid(x,y);
 
name = 'sinusoid';
for f = [2 4 8 16 32] //frequency
  k = 'a_';
  z = sin(2*%pi*f*X);
  //theta = 45;
  //z = sin(2*%pi*f*(Y*sin(theta) + X*cos(theta)));
  //z = sin(2*%pi*f*X).*sin(2*%pi*f*Y);
  //z = sin(2*%pi*f*X).*sin(2*%pi*f*Y) + sin(2*%pi*f*(Y*sin(theta) + X*cos(theta)));
  //z = sin(2*%pi*f*X).*sin(2*%pi*f*Y) + sin(2*%pi*2*3*f*X).*sin(2*%pi*2*3*f*Y);
  //scf(); imshow(z,[]);

  Im = z;
  Im = (Im - min(Im))/(max(Im)-min(Im));
  imwrite(Im, strcat([PATH,k,name,string(f),'.bmp']));
  output = abs(fftshift(fft2(z)));
  Im = output;
  Im = (Im - min(Im))/(max(Im)-min(Im));
  imwrite(Im, strcat([PATH,k,'fft_',name,string(f),'.bmp']));

  //bias
  z2 = z + abs(min(z));
  Im = z2;
  Im = (Im - min(Im))/(max(Im)-min(Im));
  imwrite(Im, strcat([PATH,k,name,'_bias_',string(f),'.bmp']));
  output = abs(fftshift(fft2(z2)));
  Im = output;
  Im = (Im - min(Im))/(max(Im)-min(Im));
  imwrite(Im, strcat([PATH,k,'fft_',name,'_bias_',string(f),'.bmp']));
end
  

Reference:

[1] M.N. Soriano, "Applied Physics 186 - Properties of the 2D Fourier Transform", 2010.

Credits:

Self Evaluation:

Activity 6 - Fourier Transform Model of Image Formation

To become familiarized with discrete FFT, a 128 x 128 image of a white centered-circle with black background was created (see Figure 1.a.). This image was converted to grayscale and then the fourier transform of this 2D image was taken. The intensity values were calculated and are plotted as shown in Figure 1.b. The fftshift() function was used to shift the peak onto the center of the image as shown in Figure 1.c Lastly, the fft2() was applied to the image in Figure 1.c and obtained an image as shown in Figure 1.d, which is similar to the original image. This method was done the same with a 128 x 128 image of a letter "A" as shown Figure 2. Observe that the image shown in Figure 2.d is an inverted image of Figure 2.a.

  
stacksize(100000000);
PATH = 'C:\Users\Micent\Documents\Applied Physics\acad 2010-2011 1st sem\App Phy 186 - Instrumentation Physics II\Activity 6\';
 
//centered circle
x = linspace(-64,64,128);
[X,Y] = meshgrid(x);
r = sqrt(X.^2 + Y.^2);
circle = zeros(128,128);
circle = zeros(size(X,1),size(X,1));
circle(find(r<=16)) = 1.0;
//imshow(circle,[]);
imwrite(circle,strcat([PATH,'circle_16.bmp']));
  

(a)

(b)

(c)

(d)

Figure 1. Results from the Simulation of an imaging devive for a circle.

//image

filename = 'circle'

//filename = 'letterA'

Igray = gray_imread(strcat([PATH,filename,'.bmp']));

  
 //fast fourier transform for 2D
FIgray = fft2(Igray);
//scf(1); imshow(abs(FIgray),[]);
//xs2bmp(1,strcat([PATH,filename,'_FIgray','.bmp']));
Im = abs(FIgray);
Im = (Im - min(Im))/(max(Im)-min(Im));
imwrite(Im,strcat([PATH,filename,'_FIgray','.bmp']));
  

//shift the FT
shiftFIgray = fftshift(abs(FIgray));
//scf(2); imshow(shiftFIgray), []);
//xs2bmp(2,strcat([PATH,filename,'_FIgray_shift','.bmp']));
Im = shiftFIgray;
Im = (Im - min(Im))/(max(Im)-min(Im));
imwrite(Im,strcat([PATH,filename,'_FIgray_shift','.bmp']));
  
  
//fast fourier transform
FIgray2 = fft2(FIgray);
//scf(3); imshow(abs(FIgray2),[]);
//xs2bmp(3,strcat([PATH,filename,'_FIgray2','.bmp']));
Im = abs(FIgray2);
Im = (Im - min(Im))/(max(Im)-min(Im));
imwrite(Im,strcat([PATH,filename,'_FIgray2','.bmp']));
  

(a)

(b)

(c)

(d)

Figure 2. Results from the Simulation of an imaging devive for a letter A.

In the second part of the activity, convolution was used for the simulation of an imaging device. Given two 2D functions f and g, the convolution between this functions is

Equation 1

Equation 1 can be similarly written in Equation 2 with "*" as the shorthand notation for convolution.

Equation 2

Following the convolution theorem, f and g can be recasted by linear transformations, such as Laplace or Fourier transform, and with H, F, and G as the transforms of h, f and g, respectively, Equation 3 can be obtained

Equation 3

Figure 3. Image of the letters "VIP" in a 128 x 128 black background.

Figure 3 shows a 128 x 128 image of the letters "VIP", which served as an object for the simulation. A 128 x 128 image of a white centered-circle was created over a black background. This circle represented the aperture of a circular lens. Both image were converted to grayscale. The 2D fourier transform of the object (VIP image) was taken and fftshift() was only applied to the aperture since it is already in the Fourier Plane. The product of their FFT is calculated and then the inverse of the product was determined to obtain the convolved image.

The process was repeated for different radii of the white circle such as 16 px, 32 px, 48 px and 64 px. Figure 4 shows the different outputs as the radius of the circular lens was varied. The radius of the aperture of lens shows that its size affects the quality of the image formed. The smallest radius produced a low quality image because the lens can only gather a limited number of rays.

  
filename1 = 'circle_16'
filename2 = 'VIP';
I = gray_imread(strcat([PATH,filename1,'.bmp'])); //circle
VIP = gray_imread(strcat([PATH,filename2,'.bmp'])); //vip
  
FI = fftshift(I);
FVIP = fft2(VIP);

FRVIP = FI.*FVIP;
IRVIP = fft2(FRVIP); //inverse FFT
VIPimage = abs(IRVIP);
Im = VIPimage;
Im = (Im - min(Im))/(max(Im)-min(Im));
imwrite(Im,strcat([PATH,filename2,'_JImage_16','.bmp']));

(a)

(b)

(c)

(d)

Figure 4. The convolution of the circular apertures with the VIP image

resulted to the images on the right column.

For the third part this actvity, correlation function was used in pattern recognition technique called template matching in which identical patterns in a scene are being looked for. The correlation between two 2D functions f and g si given by

Equation 4

or in shorthand notation

Equation 5

Equation 6 shows how Equation 5 is related to the convolution integral, where the superscript asterisk at function f means complex conjugation.

Equation 6

Equation 7 was derived from the correlation theorem which holds for the linear transforms of functions f and g. P, F, and G ar the Fourier transforms of p, f, and g, respectively, and the the superscript asterisk stands for complex conjugation.

Equation 7

The correlation measures the degree of similarity betwee two functions. That is, the more identical they are at certain positions (x,y), the higher correlation value.

Figure 5. Image of a sentence over a white background.

Figure 6. Pattern need to be matched in the image in Figure 5.

  
filename1 = 'sentence';
filename2 = 'A'
  

S = gray_imread(strcat([PATH,filename1,'.bmp'])); //sentence
A = gray_imread(strcat([PATH,filename2,'.bmp'])); //letter A
  

FA = fft2(A); //G
FS = fft2(S); //F
R = FA.*conj(FS); //F*G
output = fftshift(fft2(R));
Im = abs(output);
Im = (Im - min(Im))/(max(Im)-min(Im));
imwrite(Im,strcat([PATH,filename1,'_output','.bmp']));
  

FA = fftshift(fft2(A));
FS = fftshift(fft2(S));
R = FA.*conj(FS);
output = (fft2(R));
Im = abs(output);
Im = (Im - min(Im))/(max(Im)-min(Im));
imwrite(Im,strcat([PATH,filename1,'_output2','.bmp']));
  

Figure 5 shows a 128 x 128 image with white background was created with the text "THE RAIN IN SPAIN STAYS MAINLY IN THE PLAIN." Another image of size 128 x 128 was created with letter "A" of the same font and size in the text in Figure 6. Both these images were imported in Scilab and converted to grayscale. Their Fourier transforms were determined. To calculate the correlation of the two transforms, element-per-element multiplication of the FT of "A" and the conjugate of the FT of the text image was performed. Lastly, the inverse FFT of the correlation was determined and the output is diplayed in Figure 7.

Figure 7. Result of the correlation of image in Figure 5 and 6.

Comparing Figure 5 and Figure 7, it can be observed that the peaks locations show where the letter "A" can be found.

In the last part of the activity, convolution integral was applied in the template matching of an edge pattern with an image. A 3x3 matrix pattern of an edge was created such that the total sum of the elements is zero. The image formed was convolved with the VIP image using the imcorrcoef() function of Scilab. Figure 8 shows the result of using different patterns and their corresponding convolved image.

  
filename2 = 'VIP';
VIP = gray_imread(strcat([PATH,filename2,'.bmp'])); //vip
  

pattern1 = [-1 -1 -1; 2 2 2; -1 -1 -1];
pattern2 = [-1 2 -1; -1 2 -1; -1 2 -1];
pattern3 = [-1 -1 -1; -1 8 -1; -1 -1 -1];
output = imcorrcoef(VIP,pattern3);
Im = abs(output);
Im = (Im - min(Im))/(max(Im)-min(Im));
imwrite(Im,strcat([PATH,'convolve_output3','.bmp']));
  

(a)

(b)

(c)

Figure 8. Edge detected using three differe nt matrix patterns.

For Figure 8.a., matrix pattern [-1 -1 -1; 2 2 2; -1 -1 -1] was used that is why horizontal edges brighten up. In Figure 8.b, the convolution of matrix pattern [-1 2 -1; -1 2 -1; -1 2 -1] with the VIP image resulted to the higlighting of the vertical edges. The use of matrix patter [-1 -1 -1; -1 8 -1; -1 -1 -1] resulted to the detection of the entire edges of the VIP image.

Reference:

[1] M.N. Soriano, "Applied Physics 186 - Fourier Transform Model of Image Formation",2010.

Credits:

Self Evaluation:

Activity 5 - Enhancement by Histogram Manipulation

In this activity, the histogram of a grayscale image was studied. The graylevel probability distribution function (PDF) of an image is equivalent to the normalized histogram of the grayscale image. The grayscale PDF of an image can be modified by back projecting the grayscale values using the cumulative distribution function (CDF) of a desired PDF. In this histogram manipulation, the quality of an image is further improved, certain image features are more enhanced, or the response of different imaging systems can be mimicked.

stacksize(100000000);

PATH = 'C:\Users\Micent\Documents\Applied Physics\acad 2010-2011 1st sem\App Phy 186 - Instrumentation Physics II\Activity 5\';

//PATH = 'C:\Users\Micent\Documents\Applied Physics\acad 2010-2011 1st sem\App Phy 186 - Instrumentation Physics II\Activity 5\images\';

filename = 'cousins';

imagefile = strcat([PATH,filename,'.jpg']);

I = gray_imread(imagefile);

//imwrite(I,strcat([PATH,'cousins_gray.bmp']));

Figure 1. A dark image taken by a Digital SLR

Figure 1 shows a dark looking image. The image was imported in Scilab using the gray_imread() function. The grayscale histogram of the image was tabulated and plotted as shown in Figure 2. The PDF of the image was determined by normalizing the histogram, that is, dividing by the total number of pixels. Afterwhich, the CDF was determined from the PDF. This was done by performing a cumulative sum on the frequencies of the graylevel. Figures 3 and 4 shows the PDF and CDF of the image, respectively.

//plot the histogram of the grayscale image
scf(0); histplot(256,I);
title('Histogram of a grayscale image');
xlabel('grayscale value'); ylabel('frequency');
xs2bmp(0,strcat([PATH,'histplot.bmp'])); //sends the graphics into a file
  
 
//tabulate the grayscale values
[H] = tabul(I,"i");
gs = H(:,1); //saves the grayscale value
f = H(:,2); //saves the frequency for each gray level
scf(1); plot(gs,f);
title('Histogram of a grayscale image');
xlabel('grayscale value'); ylabel('frequency');
xs2bmp(1,strcat([PATH,'histogram.bmp'])); //sends the graphics into a file
  
  
//calculating the Probability Distribution Function (PDF)
pdf = f/sum(f); //normalize
scf(2); plot(gs,pdf);
title('PDF (normalized histogram) of a grayscale image');//histogram plot of the grayscale image
xlabel('grayscale value'); ylabel('frequency');
xs2bmp(2,strcat([PATH,'pdf.bmp'])); //sends the graphics into a file
 
  
//calculating the Cumulative Distribution Function (CDF)
F = cumsum(pdf); //cumulative sum of the frequency
cdf = F/max(F);
scf(3); plot(gs,cdf); //cumulative sum plot of the pdf plot
title('CDF of the PDF of the grayscale image');
xlabel('grayscale value'); ylabel('frequency');
xs2bmp(3,strcat([PATH,'cdf.bmp'])); //sends the graphics into a file

Figure 2. Histogram of the grayscale image.

Figure 3. The probability distribution function (PDF) of the grayscale image.

Figure 4. The Cumulative distribution function (CDF) determined from the PDF.

In this activity, I used three different CDF, namely uniform distribution (linear), parabolic and logarithmic. For each CDF, the x-axis is the range of grayscale values of the original untreated image [0-1] and the y-axis is equal to the value of the defined function (linear, parabolic or logarithmic).

//Pixel-per-Backprojection

[r,c] = size(I);

for i = 1:r

for j = 1:c

n = find(gs==I(i,j)); //index

N = cdf(n,1); //value on the index

M = 255*N; //uniform

// M = exp(N); //logarithmic

// M = -sqrt(N); //parabolic

I(i,j) = M;

end

end

I = (I - min(I))/(max(I)-min(I));

imwrite(I,strcat([PATH,filename,'_parabolic.bmp']));

Figure 5. Modified image using a linear CDF.

Figure 6. Modified image using a parabolic CDF.

Figure 7. Modified image using a logarithmic CDF.

The histogram manipulation was done pixel-per-pixel, backprojecting the image pixel values by looking for its corresponding y-value in the desired CDF. This is the same as replacing the dark pixel values with the x-values from the desired CDF. Figure 5-7 show the three modified images using three different desired CDF. The histogram and the CDF of the modified images are shown in Figure 8-10. The CDFs shown are equivalent to those previously desired.

(a)

(b)

Figure 8. The (a) histogram and (b) CDF of the image in Figure 5.

(a)

(b)

Figure 9. The (a) histogram and (b) CDF of the image in Figure 6.

(a)

(b)

Figure 10. The (a) histogram and (b) CDF of the image in Figure 7.

I also tried using the histogram manipulation that can be found in advanced image processing software such as Photoshop or GIMP. I opened the same image in the GIMP and converted it to a grayscale image. Using the Curves function, the histogram CDF was manipulated by simply dragging the diagonal line to a desired position. Figure 11 shows the original image and its histogram. Figure 12 shows the images modified by manipulating the histogram.

Figure 11. Original image with its histogram opened in GIMP.

(a)

(b)

(c)

Figure 12. Modified Images by manipulating the histogram.

Reference:

M.N. Soriano, "Applied Physics 186 Enhancement by the Histogram Manipulation", 2010.

Credits:

Self Evaluation: