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.

ACTIVITY 12 - Color Image Segmentation

In this activity, we are going to pick out a region of interest from the rest of the image through image segmentation.

A digital image of 3D objects, shown below, was utilized in this activity. Different brightly colored candies are shown in the image.

We cropped a monochromatic region of interest in the scene. In this case, we wanted the blue candies to be our region of interest.

We notice that these objects found in the image, as well as the cropped image have shading variations which is inherent in 3D objects. To separate brightness and chromaticity (pure color) information, it is better to represent the color space by the normalized chromaticity coordinates or NCC. The normalized chromaticity coordinates can be obtained by dividing each pixel from each channel by the summation of the pixel corresponding on each channel.

This was implemented by the code below.

image =imread('candy.jpg');

patch = imread('patch.jpg');

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);

Parametric Segmentation

The Gaussian PDF in the r values was derived using the equation

where sigma is the standard deviation and mu is the average of the normalized chromaticity coordinate in red of the patch image. The same equation was used for g values.

Rmu = mean(Rp);

Gmu = mean(Gp);

Rsigma = stdev(Rp);

Gsigma = stdev(Gp);

pr = (1/(Rsigma*sqrt(2*%pi)))*exp(-((r-Rmu).^2)/(2*(Rsigma)^2));

pg = (1/(Gsigma*sqrt(2*%pi)))*exp(-((g-Gmu).^2)/(2*(Gsigma)^2));

Then, the joint probability was taken as the product of rho(r) and rho(g).

product = round(pr.*pg);

The resulting image was shown below.

Non-Parametric Segmentation

The 2D histogram of the ROI was obtained using the give code below.

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;

scf(1);

mesh(hist);

To test the correctness of the histogram, we compare the location of the peaks with the rg chromaticity diagram (at the right). Clearly, we can see that the histogram of the ROI peaks at the are corresponds to the blue color in the rg chromaticity diagram. This histogram willl be important in segmenting the image using histogram bakcprojection.

This histogram backprojection is similar to what has done in Activity 4 excepth that in this case, the histogram is in 2D. Backprojection was done by replacing the value of the pixel location by its histogram value in chromaticity space. This was implemented by the code below.

rib = R*255;

gib = G*255;

[a, b] = size(image);

np = zeros(a, b);

for i = 1:a

for j = 1:b

c = round(rib(i, j)) + 1;

d = round(gib(i, j)) + 1;

np(i, j) = hist(c, d);

end

end

scf(2);

imshow(np);

The resulting image was shown below.

By implementing two different segmentation techniques, different parts of the region of interest was segmented. Parametric segmentation enables us to segment the outer part of the blue candies while non-parametric segmentation enables us to segment the inner body of the blue candies. The non-parametric segmentation favors in segmenting the region of interest although not all ROI were segmented.

I give myself a grade of 10 for finishing this activity. I was able to implement the parametric and non-parametric segmentation on an image composed of an ensemble of colored objects.

I would not have finished my blogs without the overwhelming help of Gilbert and Rommel. Gilbert was the one who patiently discussed every activity with me. He help me a lot in understanding previous activities. Rommel, after finishing his blog was able to lend me his laptop. It was a big help considering that I did not work on my blog in a computer shop for almost 4 hours.

:-) I acknowledged all my classmates who made me feel their concerns for me. I want to extend my deepest gratitude for all of you guys.

ACTIVITY 11 - Color Image Processing

Primarily, the objective of this activity is to be able to enhanced an image that is wrongly white-balanced. Images of an ensemble of colorful objects were captured using a digital camera with different white balance settings.

White Balancing Setting: Fire

White Balancing Setting: Shade

White Balancing Setting: Fluorescent-1

White Balancing Setting: Fluorescent-2

White Balancing Setting: Fluorescent-3

White Balancing Setting: Incandescent

Different white balancing settings resulted to different qualities of images. The Fire and Shade are almost the same resulting to a slightly darker image than what was taken. Fluorescent-1 settings resulted to a brownish image. Fluorescent-2 and Fluorescent-3 produced shades of blue images. The Incandescent setting produced a more bluish image. Obviously, this image is wrongly balanced.

White Patch Algorithm

Using the image captured by the camera in Incandescent setting, the White Patch algorithm was implemented. First, a patch of an white object in real world which appears differently in the image was cropped.
The RGB values of pixels belonging to this white object: Rw, Gw, and Bw was taken. And then, all pixels in the red layer, green layer and blue layer of the image was divided by their corresponding Rw, Gw and Bw. This was implemented with the code below.

image = imread('filename.jpg');
w = imread('wpatch.jpg');

rw = sum(w(:,:,1))/length(w(:,:,1));
gw = sum(w(:,:,2))/length(w(:,:,2));
bw = sum(w(:,:,3))/length(w(:,:,3));

nr = image(:,:,1)/rw;
ng = image(:,:,2)/gw;
nb = image(:,:,3)/bw;

The resulting image was then displayed as shown below.

Compared with the wrongly white balanced image, the resulting image is less brighter as the bluish color decreases. White hues that appear bluish in the wrongly balanced image now appear more white. The White Patch algorithm assumes that the maximum response of the image is caused by the white patch. Thus, by dividing each channels of the image with the constants determined from the white patch we are essentially bringing back true colors of the image.

Gray World Algorithm

The Gray World algorithm was also implemented on the image captured using Incandescent setting. For this algorithm, the R, G and B layers of the unbalanced image were averaged and let to be Rw, Gw, and Bw, respectively. Then the original R, G and B layer values were divided by their corresponding constant.

image = imread('filename.jpg');

rw = sum(image(:,:,1))/length(image(:,:,1));
gw = sum(image(:,:,2))/length(image(:,:,2));
bw = sum(image(:,:,3))/length(image(:,:,3));

nr = image(:,:,1)/rw;
ng = image(:,:,2)/gw;
nb = image(:,:,3)/bw;

The resulting image was displayed below.
The image looks old after the implementation of Gray World algorithm. The blue hues also lessen and the white hues look grayish. The Gray World algorithm assumes that the average of the colors on each channels of the image taken in normal light is gray. Using the Incandescent setting in capturing the image, we have disrupted the Gray World assumptions. We force the Gray World assumption again by dividing each channel by their average. Thus, we are reaquiring the true color of the image in real world.

Doing both algorithm with the other images will result to the same image as shown above for white patch and gray world algorithm. Implementing both algorithms is reacquiring the true colors of the image based on certain assumptions. However, it turns out that the White Patch is better than the Gray World algorithm.
------------------------------------------------------------------------------------------------------------

In this part of the activity, the white patch and gray world algorithm were implemented on an image of objects with the same hue (in this case, green).
The image above was obviously wrongly white balanced. After implementing the two algorithm as in above procedure, the resulting images were shown below.

White Patch Algorithm

Gray World Algorithm

The white patch algorithm uses the white wall at the background of the image. Clearly, the original color was reacquired using white patch algorithm. Gray world algorithm allows the removal of the bluish hue in the image. However, the image becomes slightly reddish. This is because the unbalanced image mainly composed of green hues.

For this activity, I will get 10/10. By producing images from each algorithm, I was able to understand how the white patch and gray world algorithm assumes the color of the real world. Applying these algorithm on an image captured with colored lighting (i.e. Incandescent white balancing setting ) will reacquire the true colors of the image. However, after finishing the activity, it turns out that for the images used, the White Patch algorithm assumes a better color of the real world. I acknowledged Gilbert for discussing the activity with my and for lending me pictures he captured.