Fourier Transform Model of Image Formation

July 8, 2008

In this activity, we try to get the Fourier transform of an image.

A. Familiarization

The images are created in Paint. Programs are run in Scilab. First, we investigate a circle such as the one below. The input image size is 128 x 128.


Performing the built-in Fast Fourier transform (FFT) algorithm in Scilab on it, we get quadrants of concentric circles positioned in the corner of the FFT image.



We shift the FFT so that it is symmetric about the origin.
We obtain the Airy disk.





Next, we render the Airy disk in Scilab's hotcolormap, which is composed of shades of red and yellow.
B. Simulation of an Imaging Device

Now, suppose we have a camera which has a lens aperture. Upon passing through the aperture, the signals from this object are distorted or convoluted. Convolution of two functions f and g, representing the object and the aperture respectively, is equivalent to multiplying element-per-element their Fourier transforms, and getting the Fourier transform of their product. The resulting is what we see when the image forms at the back of the lens.

For example, we have this object below.


Its Fourier transform is,




The aperture is already in Fourier space, so we multiply it with VIP's FFT. Then, we get the FFT again. The result is a blurred mirror reflection, with visible fringes, of the object.










When we replace the lens with a larger one, the object is better resolved. Meanwhile, shifting to a smaller lens distorts the image more. This must be because the zero-order of the Airy disk broadens as the aperture gets smaller, thus causing more distortion.



































C. Template Matching Using Correlation


We also tried the application of the Fourier transform in Template Matching.
We have on the left side the input image, and below it the pattern we want to trace. Using the following code, to track down A, an inverted and blurred image of the input is retrieved. Notice the locations of the bright spots. They correspond to the positions of A in the input.




m = imread('C:\Documents and Settings\instru\Desktop\ap186(2)\fig11.bmp');
k = imread('C:\Documents and Settings\instru\Desktop\ap186(2)\fig12.bmp');
mgray = im2gray(m);
kgray = im2gray(k);
M = fft2(mgray);



K = fft2(kgray);
MK = M.*conj(K);
invMK2 = fft2(MK);

invMK2 = fftshift(invMK2);
imshow(abs(invMK2),[]);












D. Edge Detection

Let's go back to VIP. Now our aim is to detect the horizontal (top image) and vertical edges (bottom image).



m = imread('C:\Documents and Settings\instru\Desktop\ap186(2)\fig5.bmp');
k = [-1 2 -1; -1 2 -1; -1 2 -1];
mgray = im2gray(m);
l = imcorrcoef(mgray,k);
imshow(abs(l),[]);




















SELF -EVALUATION: 10/10.
FFT is amazing, especially in edge detection.



Reference:
[1] Dr. Soriano's lecture notes

A5 - Physical Measurements from Discrete Fourier Transforms

To obtain the Fourier transform of images, we apply the two-dimensional discrete Fourier transform (DFT) algorithm.

The input image is the spatial domain equivalent of the sinusoid which is in temporal domain. The variable used to represent an input image pixel in the above formula is x n,m.

In Scilab, two-dimensional Fourier transform is executed through the command fft2(img) where img is the image variable.

In measuring physical signals, one need to know the signal frequency, the sampling rate must be twice the frequency of the signal. The sampling interval must be equal or less than the inverse of the sampling rate.

a) Light from a fluorescent lamp is known to flicker at 120Hz. What should be the threshold sampling interval for accurate FT analysis?

The threshold sampling interval dt should be dt = 1/(2*120Hz) or 0.00417 s.

b) What is the effect of increasing the number of samples N in the FT?

Guess: Increasing the number of samples N would also change the time interval dt. This would make the FT more accurate as it is based on the assumption that the signal is infinite and periodic.

Right answer: increasing the number of samples N means decreasing the sampling interval. This makes the FFT resolution better so that the peaks become narrower.

c) What is the effect of decreasing the sampling interval Δt in the FT?
Decreasing the sampling interval Δt in the FT would increase the number of samples N for the same time period. This would make the FT more accurate.

d) What is the effect of fixing the total time interval T but increasing the
number of samples N?
Fixing the total time interval T but increasing the number of samples is also good because it fixes the number

SELF-EVALUATION: 9/10. The guesses were not always correct but the simulations were carried out correctly.

References:
[1] Dr. Maricor Soriano's lecture notes
[2] http://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm.
[3] Nyquist Theorem, http://searchcio-midmarket.techtarget.com/sDefinition/0,,sid183_gci812005,00.html

Image Enhancement

July 1, 2008
10:28 pm







This is a continuation of the image enhancement via cumulative distribution function (CDF) manipulation. This time the desired CDF is of the form G(z) = k tanh (z-128)/125 where z is the gray level pixel value and k is a positive constant here chosen to be 16. The original image is shown on the left.















Its enhanced image by forcing the CDF to be equal to G is shown on the right. Notice the changes in intensity. While there is more contrast in the original image, the modified image is more pleasing to look at.



















What happens to their histograms and CDFs. The histogram of the altered image become
s bell-shaped. Its CDF takes the shape of a tanh which is our desired CDF.

SELF-EVALUATION RATING: 10/10 because I got the tanh shape for the output image CDF


Code:
If A is the original image with pixel values ranging from 0 to 255, and if its Cumulative Distributive Function is represented as cdf and B is the new image with a tanh shape CDF, then to obtain B:
z = [0:255];
G = (1+tanh(16*(z-128)/255))/2;



dimensions = (size(A,1)*size(A,2));
B = zeros(size(A,1),size(A,2));
for i = 0:255
m = cdf(i+1);
desired = interp1(G,z,m);

[x,y] = find(A==i);
for k = 1:length(x)
B(x(k),y(k)) = desired;

end
end

Image Enhancement

June 30, 2008

11:51 am

The goal for this activity is to enhance the contrast of an image by manipulating its cumulative distribution function (CDF). The standard way of doing this is by applying a transformation to the image such that its (CDF) is a straight line. We do this by replacing the gray level values of each pixel by the CDF value of its gray level. This is called histogram equalization. This is to distribute the intensity over the whole range of gray level values, i.e., 0-255. We may also create our own image enhancement protocols. I did it by replacing the gray level values of each pixelin the original image by the squared CDF value of its gray level.

Images and Their Corresponding Histograms

a) original image

b) histogram-equalized image

c) custom enhanced by replacing the original pixel values by the CDF value squared of the original

The CDF for (c) takes the shape of a quadratic function that opens to the right. It is fitted to the function y = sqrt(x/255).

Codes:

This code. It aims to perform histogram equalization and get the probability distribution function and cumulative distribution function of both the original and the enhanced. Its histogram-getting part is taken from Mr. Tugaff's code.

A=imread('sunflower3.jpg');

A = A-1; //because the image that would be read would have pixel values from 1 to 256.

val=[];

num=[];

counter=1;
for i=0:1:255

[x,y]=find(A==i);

val(counter)=i;

num_(counter)=length(x);

counter=counter+1;

end

num = num_./((size(A,1))*size(A,2));

plot(val, num_);

cdf = [];c

cdf(1) = num(1);

for i = 2:256

cdf(i) = num(i) + cdf(i-1);

end

plot(val,cdf)

dimensions = (size(A,1)*size(A,2));

B = zeros(size(A,1),size(A,2));

for i = 0:1:255

[x,y] = find(A==i);

for k=1:length(x)

B(x(k),y(k)) = (cdf(i+1) - min(cdf))/(dimensions -min(cdf));

end

end

B = B/max(max(B));

B = round(B*255);

counter = 1;

val2=[];

num = [];

num2_ = [];

for i=0:1:255

[x,y]=find(B==i);

val2(counter)=i;

num2_(counter)=length(x);

counter=counter+1;

end

num2 = num2_./((size(B,1))*size(B,2));

plot(val2, num2_);

cdf2 = [];

cdf2(1) = num2(1);

for i = 2:256

cdf2(i) = num2(i) + cdf2(i-1);

end

plot(val2,cdf2)

To get the images for section (c), custom enhancement, we simply this part of the above code

B = zeros(size(A,1),size(A,2));
for i = 0:1:255
[x,y] = find(A==i);
for k=1:length(x)
B(x(k),y(k)) = (cdf(i+1) - min(cdf))/(dimensions -min(cdf));
end

by

B = zeros(size(A,1),size(A,2));
for i = 0:1:255
[x,y] = find(A==i);
for k=1:length(x)
B(x(k),y(k)) = (cdf(i+1) - min(cdf))*(cdf(i+1) - min(cdf))/(dimensions -min(cdf));
end

SELF-EVALUATION RATING: 10/10 because I got a straight line for the CDF of the histogram-equalized image.

Credits:

Orignal sunflower photo was taken from: Christopher Maxwell, PhotographerWedding Photography: Black and White Info, http://www.christophermaxwell.com/black-white.htm

The part of the above code that get the histogram of an image is lifted from a code by Mr. Jeric Tugaff.

Digital Scanning part 2

June 26, 2008

Thursday

The goal for today is to get the area of a real object via Green's Theorem. My chosen object is a leaf.

The leaf image I used for processing was scaled to half of the above image. As can be seen a length of 1 centimeter in the vertical direction ate the same number of pixels as that in the horizontal direction. So no further rescaling was required. The image was saved with a bit depth of 8 using Paint software. Its histogram was then used to find out which gray level should be the threshold. It was found to be at gray level = 100 or 0.4. When the image had been converted to binary, its colors were afterwards inverted so that the leaf is bright and the background is dark. Then Green's Theorem was applied. (See preceding post for my code for executing the theorem.)

Scilab code for converting the leaf image to binary:

A = imread('leaves2.jpg');

histplot([0:255],A); //to determine the appropriate threshold level

bw = im2bw(A, 0.4); //threshold level= 0.4

imwrite(bw,'bw2.jpg');

Via the pixel counting method, getting sum(sum(bw)), the resulting area was 8491 pixels.

Meanwhile, Green's Theorem yielded an area of 8245 pixels (%difference from the pixel counting result = 2.9%). This is approximately 21 cm squared.

............

I rate myself: 10/10! because I think 21 cm squared is reasonable. The leaf's length was about 10 cm and its width, measured across the middle, was nearly 2 cm.

Credits: Benj and Ed for helping with syntax.

Different Image Types according to Bit Depth

June 24, 2008
An image of each type is opened in Scilab to get their pixel size and bit count - one for binary, one for grayscale, one for true color, and one for indexed image. Since most of the internet files are saved as JPEG. They were converted to other image types using Paint.

A. Binary
Original photo from www.unm.edu
File size : 2922 bytes
Format: BMP
Size: 148 x 143
Depth: 8
Number of Colors: 2
ResolutionUnit: centimeter
XResolution: 37.800000
YResolution: 37.800000

B. Grayscale
Original photo from purplezoe.blogspot.com
FileName: C:\Documents and Settings\AP186user20\Desktop\ap186_june24\grayscale.bmp
FileSize: 18254
Format: BMP
Width: 150
Height: 113
Depth: 8
StorageType: indexed
NumberOfColors: 256
ResolutionUnit: centimeter
XResolution: 37.800000
YResolution: 37.800000

C. True Color
photo from twisted physics.typepad.com
FileName: C:\Documents and Settings\AP186user20\Desktop\ap186_june24\truecolor.jpg
FileSize: 2202
Format: JPEG
Width: 135
Height: 90
Depth: 8
StorageType: truecolor
NumberOfColors: 0
ResolutionUnit: centimeter
XResolution: 72.000000
YResolution: 72.000000

D. Indexed
photo from chamorrobible.org
FileName: C:\Documents and Settings\AP186user20\Desktop\ap186_june24\indexed.png
FileSize: 9610
Format: PNG
Width: 150
Height: 107
Depth: 8
StorageType: truecolor
NumberOfColors: 0
ResolutionUnit: centimeter
XResolution: 72.000000
YResolution: 72.000000

Credits:
Beth Prieto for helping me with scilab

Area Measurement I

June 19, 2008 11:21 am
Grade: 10/10

The goal for today is to determine the area of a regular shape saved as a binary image, which is created in Paint.

Area is computed for using Green's Theorem, an approximation for contour integration.
My code, which has been written for scilab is as follows.
A = imread('C:\Documents and Settings\AP186User20\Desktop\circ4.bmp');
[x,y] = follow(A);
length_x = size(x,1);
length_x_minus1 = length_x-1;
x_(2:length_x) = x(1:length_x_minus1);
length_y = size(y,1);
length_y_minus1 = length_y-1;
y_(2:length_y) = y(1:length_y_minus1);
Area =(sum( x.*y_ - y.*x_))/2
theo = sum(sum(A))
error = 100*((theo - abs(Area))/theo)

The pixel count of the binary image is taken to be the theoretical area.

There were three shapes investigated: triangle, rectangle and circle. For the circle, the effect of radius on the error is investigated. Error has been found to be decreasing with increasing radius.

Triangle
Area calculated: 3162
Theoretical area = 3326
%error = 4.93%

Rectangle
Area calculated: 2517.5
Theoretical area: 2652
%error: 5.07%

Circle (radius = 15)
Area calculated: 600.5
Theoretical area: 649
%error: 7.47%

Circle (radius = 25)
Area calculated: 1707
Theoretical area: 1788
%error : 4.53%

Circle (radius = 35)
Area calculated :3509.5
Theoretical area: 3625
%error: 3.19%

Credits:
Thanks to Beth Prieto for helping me with scilab; Rica for loading ImageMagick; Aiyin for a pdf copy of the emailed attachment; Paul & JC for letting me use their internet; and to Ma'am Jing for pinpointing that scilab can have memory problems.