ACTIVITY 19 - Restoration of Blurred Image

This activity is a demonstration on how to restore an image which is corrupted with a known degradation function such as motion blur and additive noise.

The image degradation and restoration process can be modeled by the diagram shown below.

A degradation function H along with an additive noise terms, n(x,y) operates on an input image f(x,y) producing a degraded image g(x,y). With the given g(x,y) and some knowledge about the degradation function H and additive noise terms n(x,y), we can obtain an estimate restoration f'(x,y) of the original image.

The degraded image is given in the spatial domain

Thus, it can be written in an equivalent frequency domain by getting the Fourier transforms of the corresponding terms.

So we need to have our original image to be transformed in the frequency domain. Now, we will start the degradation of an original image by having our degradation function that will cause the original image to blur. H(u,v) is the Fourier transform of the degradation function given by

where a and b is the total distance for which the image has been displaced in the x- and y- direction, correspondingly and T is the duration of the exposure from the opening and closing of the shutter in the imaging process. We will use a = b = 0.1 and T = 1 and investigate also on other values for these parameters.

The additive noise terms will also be transformed in frequency domain.

To restore these corrupted images, we will apply the Weiner filter expressed by

This expression is also commonly referred to as the minimum mean square error filter or the least square error filter. The terms on the expression are as follows:

This expression is very useful when the power spectrum of the noise and the original image are known. Things can be more simplified when we are dealing with spectral white noise (the power spectrum of the noise is constant). However, the power spectrum of the original image is not often known wherein another approach can be used with the expression

where K is a specified constant. This expression yields the frequency-domain estimate. The restored image must be in the spatial domain which will be obtain by taking the inverse Fourier transform of estimation.

The image utilized for this activity is shown below

**Taken from: http://www.background-wallpapers.com/games-wallpapers/final-fantasy/final-fantasy-vii.html

The following images show the degradation of the original image with the corresponding parameters used. The first set of images has T = 1 and a = b = 0.001, 0.01, and 0.1.

Their corresponding restored images are also obtained and shown below.

We observed that by having small values of a and b, we are putting less blur on the image. Thus, its corresponding Weiner filtered image yields better image.

By keeping a and b constant, T was varied. (T = 0.01, 0.1, 10, and 100 )

Their respective restored image was shown below.

As observed, lesser exposure time makes the noise more obvious with constant values of a and b. On the contrary, greater exposure time will make the blur more visible. The Weiner filter, which works on the blur applied on the original image restores the image with higher exposure time much better.

For cases with unknown power spectrum of the original image, the constant K was varied with a = b = 0.01 and T = 1 . (K = 0.001, 0.01, 0.1 and 1, respectively)

Using a constant value of K shows a significant deviation on the restoration of the degraded image. It shows that knowing the power spectrum of the ungraded image as well as the power spectrum of the added noise is crucial in filtering blurred images. However, in instances that the two parameters are unknown, choosing a good value of K will yield a more enhanced image.

I will grade myself 9/10 for this activity. I was able to understand the activity and obtain the needed output given a short working time. However, I know that there are still so much to learn from this activity. I acknowledged Gilbert for helping me finish this activity.

The code below was use in implementing the Weiner filter on the degraded images.

image = gray_imread('FF7.bmp');
noise = grand(size(image,1), size(image,2), 'nor', 0.02, 0.02);

a = 0.01;
b = 0.01;
T = 1;

H = [];
for i = 1:size(image,1)
for j = 1:size(image,2)
H(i, j) = (T/(%pi*(i*a + j*b)))*(sin(%pi*(i*a + j*b)))*exp(-%i*%pi*(i*a + j*b));
end
end
F = fft2(image);
N = fft2(noise);

G = H.*F + N;
noisyblurredimage = abs(ifft(G));

scf(1);
imshow(noisyblurredimage, []);
imwrite(normal(noisyblurredimage), 'filename.bmp');

// Weiner filtering

Sn = N.*conj(N);
Sf = F.*conj(F);
K = Sn./Sf;
W = H.*conj(H);

Fres = ((1)./H).*((W)./(W+K)).*G;
Fres = abs(ifft(Fres));
scf(2);
imshow(Fres, []);
imwrite(normal(Fres), 'filename.bmp');

ACTIVITY 18 - Noise Models and Basic Image Restoration

This activity aims to be able to be familiar with different noise models by applying them on an image and then restore the degraded image by implementing various spatial filters.

Noise are random variables that are characterized by a probability density function or PDF. In this activity, different noise models were applied on an ungraded image.

First is the Gaussian noise which is also known as normal noise model. The PDF of a Gaussian random variable, z, is given by
where z represents gray level, myu is the average value of z and sigma is its standard deviation.

Next is the Erlang or Gamma noise which is given by
The mean and variance of this density are given by
, respectively.

The Exponential noise has the PDF given by
and has the mean and variance given by
Then we have the Uniform noise with PDF given by
and has mean and variance given by
The Impulse or Salt-and-Pepper noise has PDF given by
Finally, we have the Rayleigh noise which has the PDF given by
The mean and variance of the Rayleigh noise is given by
These noises were created using the built-in function grand in Scilab and then applied on an ungraded image.

The corrupted images were then restored using different filters. First is the Arithmetic mean filter. represent the set of coordinates in a rectangular subimage window of size m x n, which has center at point (x,y). The arithmetic mean filtering process computes the average value of the corrupted image g(x,y) in the are defined by S sub xy. The value of the restored image f at point (x,y) is simply the mean computed using the pixels int he region defined by S sub xy.


The Geometric mean filter will restore an image with the use of the equation below.
The Harmonic mean filtering operations is given by the expression
This filter works better for salt noise than the pepper noise.

The Contraharmonic mean filtering operations yields a restored image based on the expression
where Q is the order of the filter. This filter is suited in treating salt-and-pepper noise. Positive values of Q will eliminate pepper noise while negative values omits salt noise. Thus, it cannot remove both noise simultaneously.


I was able to finish this activity and obtained MANY images. All files were deleted after the CSRC reset our computers. Unfortunately, I didn't have backup of my files. I guess, 5/10 is enough for me since I don't have pictures to show that I have finished this activity. :-(
Anyway, I want to thank Gilbert for guiding me in this activity.

ACTIVITY 17 - Photometric Stereo

In this activity, we estimated and then extracted the shape of an object from shadow with the use of different sources and shading models.

First, we utilized the images of of synthetic spherical surfaces that are illuminated by a far away point source.

It looks like there are no significant difference among the pictures above but the shading of the images tells much information about the surface of the object. It gives the intensity captured by the camera at point (x,y). These images are captured from the surface of the object with the sources located respectively at
This numbers were put into matrix form where each row is a source and each column is the x,y, and z component of the source.
Now we are given with I and V which is related by
We can solve for the surface normal vector by first getting g with the use of the equation
This is the reflectance of the of the object at the point normal to the surface. The surface normal vector is obtained from g divided by its magnitude
From the surface normals, we computed the elevation z = f(u,v) and the 3D plot of the shape of the object was displayed.

The surface normals (nx, ny, nz) obtained using photometric stereo are related to the partial derivative of f(x,y) as
Then, the surface elevation z at point (u,v) which is given by f(u,v) is evaluated by a line integral
Finally, the 3D plot of the shape of the object was displayed.
The shape of the object with spherical surfaces was successfully extracted and displayed.

The Scilab code used in this activity was shown below. I want to acknowledge the help of Gilbert in working on this activity.

loadmatfile('Photos.mat');

// intensity captured by camera at point (x,y)
I1 = matrix(I1, 1, size(I1, 1)*size(I1, 2));
I2 = matrix(I2, 1, size(I2, 1)*size(I2, 2));
I3 = matrix(I3, 1, size(I3, 1)*size(I3, 2));
I4 = matrix(I4, 1, size(I4, 1)*size(I4, 2));
I = [I1; I2; I3; I4];

// location of the point sources
V1 = [0.085832, 0.17365, 0.98106];
V2 = [0.085832, -0.17365, 0.98106];
V3 = [0.17365, 0, 0.98481];
V4 = [0.16318, -0.34202, 0.92542];
V = [V1; V2; V3; V4];

// calculation of surface normal vector
g = inv(V'*V)*V'*I;
magnitude = sqrt((g(1,:).^2) + (g(2,:).^2) + (g(3,:).^2))+.0001;
n = [];
for i = 1:3
n(i,:) = g(i,:)./magnitude;
end

// computation for the elevation z = f(x,y)
nx = n(1,:);
ny = n(2,:);
nz = n(3,:) + 0.0001;
dfx = -nx./nz;
dfy = -ny./nz;
z1 = matrix(dfx,128,128);
z2 = matrix(dfy,128,128);
Z1 = cumsum(z1,2); // integration from 0 to u
Z2 = cumsum(z2,1); // integration from 0 to v
z = Z1 + Z2;
scf(0);
plot3d(1:128, 1:128, z);

ACTIVITY 16 - Neural Networks

This activity is again related from the two previous activities, Activity 14 and 15. The purpose of this activity is to classify objects from their corresponding class using neural networks. The features of of the samples from the two classes used in Activity 15 was also used in this activity. The Clorets mint candy was tagged with the value of 0 while the Pillows chocolate snack was tagged with the value of 1.

The code below was made to implement the Artificial Neural Network algorithm.

clorets_train = fscanfMat('clorets_train.txt');
pillows_train = fscanfMat('pillows_train.txt');
clorets_test = fscanfMat('clorets_test.txt');
pillows_test = fscanfMat('pillows_test.txt');

cp_train = [clorets_train; pillows_train];
cp_train(:,1) = cp_train(:,1)/max(cp_train(:,1));
cp_train = cp_train';
cp_test = [clorets_test; pillows_test];
cp_test(:,1) = cp_test(:,1)/max(cp_test(:,1));
cp_test = cp_test';

rand('seed', 0);

network = [4, 4, 1];
groupings = [0 0 0 0 0 1 1 1 1 1];
learning_rate = [1, 0];
training_cycle = 1000;

training_weight = ann_FF_init(network);
weight = ann_FF_Std_online(cp_train, groupings, network, training_weight, learning_rate, training_cycle);
class = ann_FF_run(cp_test, network, weight);

** the source code was from Cole Fabro's work

The training parameters, learning rate and training cycle were tuned. It was observe that for a given training cycle, the recognition is more accurate with large learning rate. On the other hand, with the learning rate being constant, the recognition is also more accurate.

I will give myself a grade of 10/10 for this activity. Although the code was already given, I fully understand the effect of tuning the training parameters on the accuracy of the recognition. I thank Gilbert for helping me on this activity.

ACTIVITY 15 - Probabilistic Classification

This activity is related from the previous activity, Activity 14. Two classes from the pattern recognition activity was chose to be used in this activity. The Clorets and the Pillows classes will be the classified by applying the Linear Discriminant Analysis or LDA.

Like the Pattern Recognition, the purpose of LDA is to classify objects into groups according on a set of features that describe the object. Again, we need a training set with their features which will represent the predetermined groups.

The LDA process was comprehensively explained on http://people.revoledu.com/kardi/tutorial/LDA . The procedure can be explained also by the code made for this activity.

clorets_train = fscanfMat('clorets_train.txt');
clorets_test = fscanfMat('clorets_test.txt');
pillows_train = fscanfMat('pillows_train.txt');
pillows_test = fscanfMat('pillows_test.txt');

train = [clorets_train; pillows_train];
test = [clorets_test; pillows_test];

u1 = [mean(clorets_train(:,1)), mean(clorets_train(:,2)), mean(clorets_train(:,3)), mean(clorets_train(:,4))];
u2 = [mean(pillows_train(:,1)), mean(pillows_train(:,2)), mean(pillows_train(:,3)), mean(pillows_train(:,4))];

u = [mean(train(:,1)), mean(train(:,2)), mean(train(:,3)), mean(train(:,4))];

x01 = [clorets_train(1,:) - u; clorets_train(2,:) - u; clorets_train(3,:) - u; clorets_train(4,:) - u; clorets_train(5,:) - u];
x02 = [pillows_train(1,:) - u; pillows_train(2,:) - u; pillows_train(3,:) - u; pillows_train(4,:) - u; pillows_train(5,:) - u];

n1 = 5;
n2 = 5;

c1 = (x01'*x01)/n1;
c2 = (x02'*x02)/n2;

for r = 1:4
for s = 1:4
C(r, s) = (n1/(n1 + n2))*(c1(r, s) + c2(r, s));
end
end

invC = inv(C);

P1 = n1/(n1 + n2);
P2 = n2/(n1 + n2);
P = [P1; P2];

for k = 1:(n1 + n2)
f1(k) = u1*invC*test(k,:)' - (1/2)*u1*invC*u1' + log(P(1));
f2(k) = u2*invC*test(k,:)' - (1/2)*u2*invC*u2' + log(P(2));
end

f = [f1, f2];

The object with maximum f1 will be assigned to class of Clorets while the object with maximum f2 will be assigned to class Pillows.

The table below shows the corresponding feature vector of each sample from each class. The calculated mean feature for each class and the global mean vector was shown in the table.

The table below shows the results.

The values of the discriminant function obtained are very large. This is due to the large discrepancy between the pixel area and the normalized chromaticity values which were used as featured vector. However, it can be observed that the difference between the two discriminant function were so small.

I will grade myself 10/10 for completing this activity. I was able to implement the Linear Discriminant Analysis thus I was able to classify objects to their expected class. I want to acknowledged Gilbert for helping me finish this activity.

ACTIVITY 14 - Pattern Recognition

This activity was done to be able to qualify different objects through their various characteristics such as shape, size and color by pattern recognition. A set of features of the objects will be used as pattern. We aim to identify from where class an unknown object belongs.

First, three sets of objects were assembled. Each set represent a class composing of 10 samples. This activity utilized Clorets mint candy, Pillows chocolate snack and Potchi gummy candy. The picture of the objects obtained was shown below.

Each of the three class was divided into two sets. The first five samples will be the training set while the remaining five will serve as the test set. The pixel area, the average normalized chromaticity values in red, green and blue will be used as features for pattern recognition. These features were obtained using the techniques that have been learned in previous activities. The pixel area was determined by thresholding, inverting and counting the number of pixels of the black and white image of the set objects. On the other hand, the chromaticity values were obtained by using the non-parametric segmentation technique. After obtaining each feature, they are arranged into matrix form producing the feature vector. The training feature vectors were separated from the test feature vectors. From the training feature vectors, the mean feature vector was determined from the expression
Here, N is the total number of samples in the class wj and xj is feature vector in the training set of the corresponding class. Basically, this is just the mean of the pixel area and the normalized chromaticity values for red, green and blue from 5 samples comprising the training set. Each sample from the test set will be compared to the mean feature vector, which is a 1 x 4 matrix. This will be done through minimum distance classification where the Euclidean distance is applicable given by
where x is the feature vector of the test set and j is the number of classes which in this case is 3.

The feature vector of every sample on each set and their corresponding mean feature vector is shown in the table below.
The resulting mean distance and the classification of the samples from the test set was shown in another table below.
As have been observed, the calculated mean distance, D is relatively larger when the test sample does not belong to the class. Relatively small mean distance tells that the test sample belongs to the class. However, there is a deviation between the Clorets and Potchi sample. Some Clorets test samples were classified as Potchi while some Potchi samples were classified as Clorets. We can say that the deviation was brought by the shape of both sample. Overall, the technique was useful in classifying objects of different characteristics.

I will give myself a grade of 10/10 in this activity. I was able to understand how to use the characteristics of objects as pattern in classifying them through pattern recognition. I had a hard time capturing pictures of the samples that can be threshold correctly in order to use the bwlabel function of Scilab in getting the pixel area. After obtaining the picture, acquiring the data was made easy with the help of Gilbert.

The code below was utilized in this activity.


image1 = 1-gray_imread('bwclorets1.bmp');

// Area computation for training set
[L, n] = bwlabel(image1);
area = [];
for i = 1:n
area(i) = sum(i==L);
end
mean_area = mean(area(1:5));

for i=1:10
image =imread('clorets'+string(i)+'.bmp');
patch = imread('patchclorets.bmp');

r = image(:,:,1);
g = image(:,:,2);
b = image(:,:,3);

rp = patch(:,:,1);
gp = patch(:,:,2);
bp = patch(:,:,3);

R = r./(r+g+b);
G = g./(r+g+b);
B = b./(r+g+b);

Rp = rp./(rp+gp+bp);
Gp = gp./(rp+gp+bp);
Bp = bp./(rp+gp+bp);

Rmu = mean(Rp);
Gmu = mean(Gp);
Rsigma = stdev(Rp);
Gsigma = stdev(Gp);

// Non-parametric Segmentation
BINS = 256;
rint = round( Rp*(BINS-1) + 1);
gint = round (Gp*(BINS-1) + 1);
colors = gint(:) + (rint(:)-1)*BINS;
hist = zeros(BINS,BINS);
for row = 1:BINS
for col = 1:(BINS-row+1)
hist(row,col) = length( find(colors==( ((col + (row-1)*BINS)))));
end;
end;

rib = R*255;
gib = G*255;

[a, b] = size(image);
np = zeros(a, b);

for r = 1:a
for s = 1:b
c = round(rib(r, s)) + 1;
d = round(gib(r, s)) + 1;
np(r,s) = hist(c, d);
end
end

// Mean color per channel
imageR = [];
imageG = [];
imageB = [];

[x , y] = find(np ~= 0);
for j=1:length(y)
imageR = [imageR, R(x(j),y(j))];
imageG = [imageG, G(x(j),y(j))];
imageB = [imageB, B(x(j),y(j))];
end
mean_imageR(i,:) = mean(imageR);
mean_imageG(i,:) = mean(imageG);
mean_imageB(i,:) = mean(imageB);
end

mean_imageR_training = mean(mean_imageR(1:5));
mean_imageG_training = mean(mean_imageG(1:5));
mean_imageB_training = mean(mean_imageB(1:5));

M=[mean_area,mean_imageR_training,mean_imageG_training,mean_imageB_training];

// Area computation for test set
image2 = 1-gray_imread('bwclorets2.bmp');
[L, m] = bwlabel(image2);
area2 = [];
for i = 1:m
area2(i) = sum(i==L);
end

// Mean Distance
D = [];
X = [];
for i=6:10
X=[area(i-5),mean_imageR(i),mean_imageG(i),mean_imageB(i)];

d = X-M;
D = [D, sqrt(d*d')];

end

ACTIVITY 4 - Enhancement by Histogram Manipulation

In this activity, grayscale images will be enhanced by manipulating their histogram which is the graylevel probabilit function or PDF when normalized. Histogram can be modified by having the cumulative distribution function (CDF) of the image and a desired CDF. The image was modified by backprojecting the grayscale values of the desired CDF.

A grayscale image with poor contrast was downloaded from the internet. Then the grayscale histogram of the image was obtained using Scilab.

grays = gray_ imread('gray_image.jpg'); //read from grayscale file

//compute and plot the histogram
j=1;
pix = [];
for i = 0:255
[x,y]=find(grays==i);
pix(j)=length(x);
j=j+1;
end





















The grayscale image was show above along with its grayscale histogram. We then determine the CDF of the image from its PDF. The CDF is given by
where p1(r) is the PDF of the image. The CDF of the grayscale image above was shown at the left side below.

We are going to remap the grayscales of the image in a way that the resulting CDF of the reconstructed image will look like our desired CDF. The CDF of a uniform distribution is a straight increasing line. If we wanted our image to have a uniform distribution of gray values, we will use a straight line, shown above right as our desired CDF.


The process of backprojection was illustrated below.


Each pixel value has a CDF value. This value was traced in the desired CDF. The pixel value will then be replaced by the new pixel value in the desired CDF.

b = pix;
c = cumsum(b);
f = c/max(c);
x = 0:255;

imsize = size(grays);
dCDF = x;

for i = 1:imsize(1)
for j = 1:imsize(2)
ind= find(x == grays(i, j));
grays(i, j) = f(index);

The resulting image was shown below.
The transformed image shows a significant increase in contrast. This image is now called a histogram equalized image.To have more idea about the image, we displayed its grayscale histogram and its CDF.
















Compared with the original histogram, the histogram of the transformed image showed a more uniform distribution of grayvalues. As with our goal, the CDF of transformed image looked like just the same as the desired CDF.
----------------------------------------------------------------------------------------------------------------
We have manipulated the histogram of a grayscale image with the use of a linear desired CDF. However, the human eye has a nonlinear reponse. So we can mimic the human eye response by using a nonlinear CDF. We can use a tanh function as our desired CDF which is nonlinear. Its plot was shown below.

x = 0:255;
tanhf = tanh(16.*(x-128)./255);
tanhf = (tanhf - min(tanh))./2;

Again, the backprojection was done. The CDF value of the pixel value was traced along the CDF value of the desired CDF. The pixel value in the desired CDF will be the new pixel value.

b = pix;
c = cumsum(b);
f = c/max(c);
x = 0:255;

imsize = size(grays);

for i = 1:imsize(1)
for j = 1:imsize(2)
ind = find(x == grays(i, j));
pixel1 = f(index);
pixel2 = find(tanhf<= pixel1); grays(i, j) = (pixel2(max(pixel2)));

The transformed image using a non-linear CDF was shown below.
Clearly, the transformed image is not as good as the image transformed using linear CDF. The image becomes darker with very poor contrast. We display its grayscale histogram and its CDF.















We can see from the grayscale histogram the distribution of gray values at a certain region at the center. The CDF shows that we were able to transformed the grayscale image base from the desired CDF. However, we were not able to enhance the contrast of the gray image.

I give myself a grade of 9 for this activity. I was able to easily understand the backprojection of pixels however, I had a hard time implementing it in a code. Also, I was able to show the effect of having a linear and non-linear CDF desirable for transforming a grayscale image. In this activity, I learned how to qualify an image with just looking on its grayscale histogram. The more uniform the distribution of the grayscale values, the better the contrast of the grayscale image. I acknowledged the help of Gilbert, who tirelessly helped me debug my code and waited for me until I gathered all my results.