Activity 6: Properties of the 2D Fourier Transform

Fourier Theorem states that any function can be expressed as infinite sum of sinusoids of different frequencies. In imaging, this translates to finding the spatial frequencies of an image. In 2D FT, rotation of the image results to rotation of the resulting FT, we will investigate this property in a while but let us first familiarize ourselves with FT of different patterns [1]

Activity 6.A: Familiarization with FT of different 2D patterns

Below is the FT of different patterns. A quick way to check if the result is correct is to imagine placing a mask(aperture) in front of a light source and observe the image in a screen far away. The resulting FT should have the somewhat the same symmetry as the object.

Observe the difference with the square above, i.e. similar to subtracting terms in your diffraction. In Fresnel diffraction, this amounts to removing zones in your image.

Textbook example of two slits along the x-axis (see Hecht, section diffraction on Fraunhoffer diffraction)

The FT of two dots is a sinousoid. This is textbook example of Young's double slit experiment.

Activity 6.B: Anamorphic Property of the Fourier Transform









Sinusoid of frequency 4 generated using sin(2*pi*f*x).

Sinusoid of frequency 8. The dots are more widely separated. Imagine placing a frequency axis on the y-axis above. The higher the frequency of the sinusoid, the farther apart the dots.

In digital imaging, there is no negative values, adding a bias to an image results to the presence of zero order. This is exactly what happens when we add bias to the image of the sinusoids above.

Sinusoid of frequency 4 with bias. The center dot is the zero-order frequency. In electronics, we can imagine a DC bias or in imaging, the image is offset by the background.

Suppose then that we took a picture of an interferogram, we can obtain the actual frequencies of the object by adding a constant bias to the image (i.e., min value). Or we just obtain the FT and neglect the zero order frequency. However, if the resulting interferogram has non-constant bias then we can add a filter to obtain the desired frequencies. For example, suppose that the image was added with very low frequency sinusoids. We can obtain the frequencies by obtaining the FT of the image and multiply the resulting FT with a high pass filter (i.e, to remove the low frequencies) and then obtain the inverse FT to get the final image of the object. Depending on the bias, we can adjust the filter to obtain the desired frequencies.

In 2D FT, rotating the sinusoids result to a rotated FT.

Rotated sinusoid with its corresponding FT. Observe that the resulting FT is also rotated in the same direction as the object.

Fourier transform of combination of sinusoids. Left: multiplication of two corrugated roofs, in the X and Y direction and right: FT of the object, 4 dots of the same distance signifying same frequency in the X and Y direction.

We added a rotated sinusoid to the corrugated image above using different frequency (f=16 and rotated for 90 degrees). Observe the presence of the 4 dots similar above and the addition of two dots.

We added a different rotated sinusoid using different frequency (f=12 and rotated for 30 degrees).

Added rotated sinusoid of f=8 and rotated for 60 degrees.

Combination of all the above rotated sinusoids. Here we can see that the result is just the addition of all the FT of the different rotated sinusoids. This is because FT is a linear transform, i.e.,

let X=A+B+C then F{X} = F{A} + F{B} + F{C}

Here, we investigated the different properties of the 2D Fourier Transform, in particular

• rotation of the object results to a rotated FT
• FT is a linear transform
• FT obtains the spatial frequencies of the image

For this activity, I give myself a grade of 10 for doing the required objectives.

I would like to acknowledge miguel and martin for useful discussions and master for knowing the spelling of Hecht =).

References:
[1] Activity 6 Manual

Activity 5: Fourier Transform Model of Image Formation

Lens as a Fourier Transform

Activity 5.A Familiarization with discrete FFT
Using the built-in FT of scilab, fft2, we are going to obtain the FFT of different patterns.

Left: image of circle, right: Fourier transform of circle. As expected, the FT is a bright spot most intense in the center, imagine imaging uniform disk of same intensity.







Left: image of letter 'A' and right: Fourier tranform of circle. The resulting FT has higher frequency components as expected from the shape of letter A.

Activity 5.B Convolution: Simulation of an imaging device
The representation of an image using an imaging device is not a perfect reconstruction but is a "smeared/convolve" image of the original image. Consider the equation below in frequency space

H=FG

Where F is the transfer function of the imaging device and G is the original image. The convolution of f and g results to the image we actually see. In frequency space, by the Convolution theorem, this is just the multiplication of the linear transform (i.e, Laplace of Fourier) of f and g. We simulate an imaging device using an aperture as our transfer function (i.e, this can be modeled as a lens) and convolve to an image of our choice. Below is the image of letters "VIP" which we are going to image using different sets of apertures.

Original image that we will image.

Left: Convolution of the original image using the aperture on the right. Observe that the resulting image is blurry (i.e., diffraction effect). We can imagine in our real life model, i.e lens, that the aperture did not collect sufficient amount of photons resulting to blurry image or the other way to put it is that the NA (numerical aperture) of the lens is very low.

Left: Convolution of the FT of the aperture in the right with the original image. Observe the significant change in resolution when we increase the size of the aperture. (increase NA)








Left: Convolution of the FT of the aperture in the right with the original image. Further increase in the aperture results to better resolution.








Left: Convolution of the FT of the aperture in the right with the original image. The largest apertuer that we can construct. Although observe that even with the highest NA, we cannot replicate the original image (i.e, presence of blurring at the edges) this is because in the mathematics of Fourier Analysis, in order to perfectly reconstruct we must have infinite coefficients to represent a function, but in our case, we are limited to discrete frequencies.

Activity 5.C: Template Matching using correlation
Correlation is basically finding the overlap of two functions. Thus, it can be used for pattern matching.








"The rain...": Image that we will correlate with the image with letter 'A'. The resulting image results with the all the letter 'A' in the image to have high intensity (appear white in the result).


Activity 6.D: Edge Detection
Using similar concept as above, we can correlate an image with different patterns and thus highlighting those parts that we need depending on the pattern.









Edge detection using the pattern in the right. As observed, the pattern can be deduced to detect all point in the edges.









Edge detection using a vertical pattern, as expected, all the vertical parts of the letter VIP are highlighted.









Edge detection using horizontal pattern. The horizontal edges are highlighted.
It must be noted that the sum of the matrix in the pattern is zero.

For this activity, I give myself a grade of 10 for performing all the requirements needed.

Activity 4: Enhancement by Histogram Manipulation

By looking at the histogram of an image, one can can determine wether the image is too bright or too dark, has poor or good contrast, i.e, if it has good dynamic range. Usually, we want an image to span all the gray levels resulting to a high contrast image. By manipulating the histogram of an image we can improve the quality of the image highlighting or enhancing certain features, or depending on the user, can mimick other imaging system[1].

The histogram of an image when normalized is equal to graylevel probability distribution function (PDF). We can alter the PDF of our image to the desired PDF by backprojecting the grayscale level values of the desired cumulative distribution function (CDF) to the grayscale values of the original CDF.

The steps in overview are as follows:
1. Per pixel in the image having grayscale r, find its CDF value T(r).
2. Find the value of T(r) in the y-axis of the desired CDF G(z).
3. Find the z value of this G(z).
4. Replace the image with grayscale value r to grayscale value z.
This process is called HISTOGRAM EQUALIZATION.

I have provided two methods in solving the problem, the difference is essentially in obtaining the PDF and CDF of the image. In the first solution, we use binning (middle of bin, i.e, 0.5, 1.5...254.5) in obtaining the histogram of the image. In the second solution, we use the function tabul of scilab which tabulates the frequency of distinct grayscale level values.

1st solution:

directory = 'C:\Documents and Settings\Orly\Desktop\AP186\Activity 4';
chdir(directory);
//load image in grayscale
img1=gray_imread('img2.png');
img2=img1; //holder
//obtain the histogram
bits=255;
bins = 0.5:bits; bins = bins./bits; //bins for [0,1] image.
for i=1:bits,
frequency(i) = sum(img1>=(bins(i)-1/(2*bits))&img1<(bins(i)+1/(2*bits))); end; frequency = frequency./max(cumsum(frequency)); f1=scf(1); plot(bins.*bits, frequency); //normalized histogram f2=scf(2); cdf = cumsum(frequency); //cumulative distribution function plot(bins.*bits,cdf); //desired cdf is a straight line, can be change descdf = (0:1.0/(length(cdf)-1):1); f3=scf(3); plot(bins.*bits, descdf); //mapping for i=1:bits-1, index = find(descdf>=cdf(i) & descdf=1 then //satisfied condtion above
//get first index
imgval = bins(index(1)); //desired pixel value
imgcur = bins(i); //current pixel value of image
img2(img2>=(imgcur-1/(2*bits))&img2<(imgcur+1/(2*bits))) = imgval; end, end; //for comparison of images f4=scf(4); imshow(img1); f5=scf(5); imshow(img2); //check the resulting CDF and PDF of the equalized image for i=1:bits, frequency2(i) = sum(img2>=(bins(i)-1/(2*bits))&img2<(bins(i)+1/(2*bits))); end; f6=scf(6); plot(bins.*bits, frequency2./max(cumsum(frequency2))); f7=scf(7); plot(bins.*bits, cumsum(frequency2)./max(cumsum(frequency2)))

The image below is the original image which we will use to demonstrate contrast enhancement.

Poor contrast imageHistogram of the original image, min = 0.5, max =181.5

CDF of the original image.
Desired CDF which is just the CDF of a uniform distribution

Using the above code, the resulting image and histograms is as follows.

Resulting image: Note that in our process, we span all the gray levels and redistribute them in the image in a uniform process. The midtones are enhanced but the problem is the image got saturated on parts with high gray level and too much shadows on parts with low gray level value. Also, the result is a bit noisy compared to the original.

Histogram of the resulting image, min = 0.5, max =253.5. The contrast of the resulting image base on the minimum and maximum gray values increased dramatically.

The CDF of the resulting image follows the shape of the desired CDF above.

Main problem encountered with the above solution is that for those images with limited dynamic range (i.e, 100-150), the resulting image would be too saturated for those pixels with higher pixel value (130-150) and too dark for those pixels with lower pixel value (100-120) because essentially, we are uniformly stretching the image into all the gray levels (0-255). The resulting image might not be too pleasing but it would have better dynamic range.

Another solution for this is to uniformly distribute the image (i.e, linear CDF) in the range of the original histogram. (i.e., from 100-150). In essence, we LOCALIZED the histogram equalization.

2nd Solution
directory = 'C:\Documents and Settings\UP DILIMAN\Desktop\ORLY\Acads\AP186\Activity 4';
chdir(directory);

//load image in grayscale
img1=gray_imread('img2.png');
info = imfinfo('img2.png');
img2=img1; //holder
bits=(2^info.Depth)-1;
//find frequency of distinct values
hist=tabul(img1,'i'); //arrange in increasing order
frequency = hist(:,2)./max(cumsum(hist(:,2))); //normalized frequency
bins=hist(:,1);

//displaying histogram
histogram=fhistogram(img1, info.Depth);
frequency2 = histogram(:,2)./max(cumsum(histogram(:,2)));
f1=scf(1);
plot(histogram(:,1).*bits, frequency2);

//find the minimum and maximum grayscale values
index=find(frequency~=0);
minB = bins(index(1))*255;
maxB = bins(index(length(index)))*255;

f2=scf(2);
cdf = cumsum(frequency);
plot(bins.*bits,cdf);
//desired cdf is a straight line, can be change
descdf = (min(cdf):(max(cdf)-min(cdf))/(length(cdf)-1):max(cdf));
desbins = (min(bins):(max(bins)-min(bins))/(length(bins)-1):max(bins));
f3=scf(3);
plot(desbins.*255, descdf);

//mapping
for i=1:length(bins)-1,
//returns the index in descdf where cdf(i) belongs
index = find(descdf>=cdf(i)); //returns a list of indices
if sum(index)>=1 then
//get first index
imgval = desbins(index(1)); //desired pixel value
imgcur = bins(i); //current pixel value
img2(img2==imgcur) = imgval;
else
x_message('cdf index is' + string(i));
end,
end;

f4=scf(4);
imshow(img1);
f5=scf(5);
imshow(img2);
f6=scf(6);
histogram2=fhistogram(img2, info.Depth);
frequency3 = histogram2(:,2)./max(cumsum(histogram2(:,2)));
plot(histogram2(:,1).*bits, frequency3);
f7=scf(7);
plot(histogram2(:,1).*bits, cumsum(frequency3));

code of filename: fhistogram.sce

function histogram = fhistogram(img, depth)
bits=(2^depth)-1;
histogram = zeros(bits,2);
bins = 0.5:bits; bins = bins./bits; //bins for [0,1] image.
for i=1:bits,
frequency(i) = sum(img>=(bins(i)-1/(2*bits))&img<(bins(i)+1/(2*bits))); end; histogram(:,1)=bins'; histogram(:,2)=frequency; endfunction
Using tabul, the bins are the distinct grayscale values and the frequency is just the occurrence of that distinct grayscale value.

Resulting image after mapping. Note that the image has better contrast compared to the reconstruction using the first method. The advantage of this reconstruction is that, the image would not be too saturated on some areas and too dark on other areas because we uniformly distributed the image in the range of the min and max grayscale values.

Resulting histogram, min=0, max=182

Resulting CDF, note that it is almost similar to the desired CDF.

By defining another set of bins ('desbins') that are also uniform in the min and max grayscale level values, we have essentially equally redistributed the histogram of the image.

Using other images:

Original: left, Contrast enhanced: right. The reconstructed image has better contrast compared to the original image. Image obtained from source [2].
Original: left, Reconstructed: right. Note that the effect of the histogram equalization is not quite visible. This is because the dynamic range of the image is limited as can be seen in the histogram and CDF of the reconstructed image. Image was obtained from source [3].
Note that in this histogram and CDF, we can see clearly how the method works. The histogram was equally distributed in the range of the histogram (min - max) hence the resulting CDF is linear in that range only.

USING OTHER CDF:
As a proof of principle, we use an exponential CDF to demonstrate what happens on the image when using nonlinear CDF. Depending on the user, we can opt to highlight the darker areas (low grayscale values), the whiter areas (high grayscale values) or the midtones (middle grayscale values). When highlighting the darker areas, we can use logarithmic CDF, when highlighting whiter areas, we can use exponential CDF and when highlighting the midtones we can use S-curve CDF.

The code above can be altered by defining the appropriate desired CDF. We just inserted the following code:
desbins = (min(bins):(max(bins)-min(bins))/(length(bins)-1):max(bins)); sigma =2; //strength of exponential descdf = exp(sigma*desbins);
descdf = (descdf - descdf(1))./max(descdf-descdf(1));

Resulting image using exponential CDF. Note that the resulting image highlighted the higher grayscale values (whiter areas) as expected from the shape of an exponential curve.

Left: Desired CDF, Right: CDF of the resulting image.

We have provided a solution for contrast enhancing an image by histogram manipulation. Our solution involves locally equalizing the histogram of an image using the available min and max pixel values of the image.

It must be noted that there is no single CDF that can enhance the contrast of every image. Each CDF has its own limitations and use. Care must be taken when contrast enhancing an image and knowledge of how each CDF alters the image must be understood.

Please email me (orly.tarun@gmail.com) if you plan to use the code above! Please leave your comments and suggestions. =)

For this activity, I give myself a grade of 10 for performing acceptable contrast enhancement algorithms.

References:
[1] Activity 4 Manual.
[2] http://myweb.lsbu.ac.uk/dirt/museum/margaret/08--452-1000120.jpg
[3] http://www.andrew.cmu.edu/user/timothyz/hw3/low_contrast.jpg