A17-Color Image Segmentation

We gathered a video of freely-falling objects, two pieces of paper with the same dimensions. One of them was crumpled. My goal is to track the distance traversed by the object from its original height as time increases. I then compared this with what should be expected theoretically. Above is the first frame of the video we gathered.

Since it is hard to segment the paper from the hand, and since in the last frames of the video, the crumpled paper had already reached the floor, I only used frames 3 to 12. For each frame, I got two portions and binarized them to track the falling papers. I then got the centroid for each portion. From the original height located in frame 1, and the y-components of the centroids, I obtained the position versus time graph for the two objects. I compared them by the theoretical values d = gt^2 where d stands for distance, g for the acceleration due to gravity and t is time. The white markings in the frames served as scale and their real world interval was 10 cm. In pixels, the interval was 33 pixels.

The final distance for the uncrumpled paper is 0.23 m, for the crumpled paper 0.60 m, and for the theoretical distance 0.78 m. We expect the final distance for the uncrumpled paper to be less than that for the crumpled paper because of drag force. The discrepancy between the measured distances versus time for the crumpled paper and the theoretical must be due to a little drag force on the crumpled paper. But at least they are of the same order of magnitude.

SELF--GRADE: 9/10

Thank you to Lei, Aiyin, Beth, Rica for the video and sliced frames.

Project - Image Deblurring

I'm still on the process of searching references:
[1] http://www.picturesolve.com,
[2] http://ieeexplore.ieee.org/iel3/4760/13241/00600843.pdf
[3] How to Sharpen an Image in Photoshop, http://www.photoshopsupport.com/tutorials/sharpen-an-image/photo-sharpening.html
[4] Sharpen Effect, http://www.websupergoo.com/helpie/source/2-effects/sharpen.htm

Hehe, we changed topic. Sharpening of images taken under a microscope. Due to the limited depth of field of the microscope, and the sample's varying topography, some parts of the image become blurred. The sharpening filter to be used is a high-boost filter given by
H = (1/9)*[-1 -1 -1; -1 w -1; -1 -1 -1]; (1)

where w = 9A-1 with A>1. The effect of the above mask is an ouput image related to the original by this equation
High-boost Output = (A - 1) Original + Highpass [2].

A20 - Neural Networks

This is another way of classifying piattos and pillows chips. First, we train our system. We have an input feature set. The columns are the features (red-and-green contrast, and area normalized to the largest area in the set) and the rows are the individuals.

x =[0.313023 0.7962447;
0.2596636 1.;
0.2721711 0.7661728;
0.3666842 0.842306;
0.8614653 0.8345313;
0.8959559 0.7132170;
0.9718898 0.5795805;
0.9224472 0.5188499]';

The first four rows pertain to piattos chips while the last four refer to pillows chips. We designate piattos as 0 and pillows as 1.Then our target set is
t = [0 0 0 0 1 1 1 1];
We input this into the neural network program and choose a test set whose member classifications we do not know.

testset = [0.3322144 0.8215565;
1.0195367 0.3562358;
0.3121461 0.9116466;
1.043473 0.3339846;
0.4000078 1.;
1.0175605 0.3679583;
0.3930316 1.;
0.9543794 0.2963204 ];

The test result is
0.1190987
0.9510804
0.0952159
0.9528344
0.1154394
0.9507229
0.1122649
0.9477875

which when binarized becomes
0
1
0
1
0
1
0
1.

Our system has classified the test samples with perfect accuracy. Indeed the test set comprises of alternating piattos and pillows photos.

SELF-GRADE: 10/10

Acknowledgment: Jeric Tugaff for the code and Cole Fabros for explaining it to me.

Mr. Tugaff's code modified to process my data:
// Simple NN that learns 'and' logic
// ensure the same starting point each timerand('seed',0);
// network def.// - neurons per layer, including input//2 neurons in the input layer, 2 in the hidden layer and 1 in the ouput layerN = [2,2,1];
// inputsx = [0.313023 0.7962447; 0.2596636 1.; 0.2721711 0.7661728; 0.3666842 0.842306; 0.8614653 0.8345313; 0.8959559 0.7132170; 0.9718898 0.5795805; 0.9224472 0.5188499]'; x2 = [0.3322144 0.8215565; 1.0195367 0.3562358; 0.3121461 0.9116466; 1.043473 0.3339846; 0.4000078 1.; 1.0175605 0.3679583; 0.3930316 1.; 0.9543794 0.2963204]';
// targets, 0 if there is at least one 0 and 1 if both inputs are 1t = [0 0 0 0 1 1 1 1];// learning rate is 2.5 and 0 is the threshold for the error tolerated by the networklp = [0.1,0];
W = ann_FF_init(N);
// 400 training cylesT = 1000;W = ann_FF_Std_online(x,t,N,W,lp,T);

A19-Probabilistic Classification

I used data from the previous activity. There were again two classes (piattos and pillows), each with four training images . There was only one test image chosen for analysis. The features set was the same as in the previous activity.

----RG----Area--Classification----

0.3130230 5428. piattos
0.2596636 6817. piattos
0.2721711 5223. piattos
0.3666842 5742. piattos
0.8614653 5689. pillows
0.8959559 4862. pillows
0.9718898 3951. pillows
0.9224472 3537. pillows

We have our global feature set

x = [0.3130230 5428.;

0.2596636 6817.;

0.2721711 5223.;

0.3666842 5742.;

0.8614653 5689.;

0.8959559 4862.;

0.9718898 3951.;

0.9224472 3537.];

and its classification vector

y = [1;1;1;1;2;2;2;2];

We separate the global feature set to class feature sets

x1 = [0.3130230 5428.; 0.2596636 6817.; 0.2721711 5223.; 0.3666842 5742.];
x2 = [0.8614653 5689.; 0.8959559 4862.; 0.9718898 3951.; 0.9224472 3537.];

We get the mean for each feature in each group

u1 = [0.3028855 5802.5];

u2 = [0.9129396 4509.75];

and the global mean vector

u = [0.6079125 5156.125];

Then we get the mean-corrected data (xnaught1 and xnaught2) and the covariance matrix for each group (c1 and c2)

xnaught1 = [- 0.2948895 271.875;

- 0.3482489 1660.875;

- 0.3357414 66.875;

- 0.2412283 585.875];

xnaught2 = [0.2535528 532.875;

0.2880434 - 294.125;

0.3639773 - 1205.125;

0.3145347 - 1619.125 ];

c1 = [0.0947876 - 205.58834;

- 205.58834 795035.89];

c2 = [0.0946674 - 224.37948;

- 224.37948 1111089.3];

The pooled within group matrix C is then

C = [0.0947275 - 214.98391;

- 214.98391 953062.61];

Its inverse is

inv(C) = [21.629463 0.0048790;

0.0048790 0.0000021];

The prior probability vector is p = [4/8; 4/8];

The linear discriminant functions yield

----f1----- ----f2----- Class

40.193166 38.215655 1.

57.712379 55.641361 1.

35.908829 33.609495 1.

46.444826 44.89887 1.

62.954231 64.805784 2.

52.618282 54.544249 2.

42.555135 44.824398 2.

35.055327 36.902364 2.

Let our test image be that of a member of the piattos group (suppose we do not know this fact). It has a feature vector of

test = [0.3322144 0.8215565];

Its f1 and f2 values are 43.269644 and 41.458361, respectively.

We show the linear discriminant functions plot and use it to classify the test image. Through the graph, we can say that the test sample is a piattos chip.

Code:

x = [0.313023 0.7962447; 0.2596636 1.; 0.2721711 0.7661728; 0.3666842 0.842306; 0.8614653 0.8345313; 0.8959559 0.7132170; 0.9718898 0.5795805; 0.9224472 0.5188499];

y = [1;1;1;1;2;2;2;2];

x1 = [0.313023 0.7962447; 0.2596636 1.; 0.2721711 0.7661728; 0.3666842 0.842306];

x2 = [0.8614653 0.8345313; 0.8959559 0.7132170; 0.9718898 0.5795805; 0.9224472 0.5188499];

u1 = mean(x1,'r');

u2 = mean(x2,'r');

u = mean(x,'r');

xnaught1 = [];

for i = 1:size(x1,1);

xnaught1(i,:) = x1(i,:) -u;

end

xnaught2 = [];

for i = 1:size(x2,1);

xnaught2(i,:) = x2(i,:) -u;

end

c1 = xnaught1'*xnaught1/size(xnaught1,1);

c2 = xnaught2'*xnaught2/size(xnaught2,1);

C = (4*c1 + 4*c2)/8;
Cinv = inv(C);

P = [4/8; 4/8];

f1 = u1*Cinv*x' - (u1*Cinv*u1')/2 + log(P(1));

f2 = u2*Cinv*x' - (u2*Cinv*u2')/2 + log(P(2));

test = [0.3322144 0.8215565];

f1test = u1*Cinv*test' - (u1*Cinv*u1')/2 + log(P(1));

f2test = u2*Cinv*test' - (u2*Cinv*u2')/2 + log(P(2));

SELF-GRADE: 9/10 . I enjoyed the activity but it took me long to do this.

A18 - Pattern Recognition

For this activity, the aim is to automatically classify images of piattos and pillows chips using features extracted from the training images. There were four training images for the piattos and another four images for the pillows. The features used were red-and-green contrast and area.
The area was obtained by simply binarizing the images and summing the pixels. red-and-green contrast is defined by the following.
We let r be the contrast in the red channel: r = (max(r) - min(r))/(max(r) + min(r))
We let g be the contrast in the green channel: g = (max(g) - min(g))/(max(g) + min(g))
The red-and-green contrast is given by the rg = sqrt(r*r + g*g).

I designated piattos as 1 and pillows as 2.


Piattos features:
---RG---- ----Area----
0.3130230 5428.
0.2596636 6817.
0.2721711 5223.
0.3666842 5742.

Pillows features:

---RG---- ----Area----

0.8614653 5689.

0.8959559 4862.

0.9718898 3951.

0.9224472 3537.

Then I input a series of images consisting of four consecutive piattos and four consecutive pillows.

Test features:

---RG---- ----Area----

0.3322144 7569.

0.3121461 8399.

0.4000078 9213.

0.3930316 9213.

1.0195367 3282.

1.043473 3077.

1.0175605 3390.

0.9543794 2730.


The output of the program below, which made use of minimum distance classification, was

1 1 1 1 2 2 2 2

The chips were correctly classified with 100% accuracy.

If I add another class, the vcut class, using the same features, the accuracy will drop to 50%. This is because piattos and vcut nearly have the same color and the same area.

//for piattos
piattos = [];
for i = 1:4
M = imread(strcat('D:\ap186\september17\piatos'+string(i)+'.jpg'));
r = M(:,:,1);
g = M(:,:,2);
contrast_r = (max(r)-min(r))/(max(r)+min(r));
contrast_g = (max(g)-min(g))/(max(g)+min(g));
piattos(1,i) = sqrt(contrast_r*contrast_r + contrast_g*contrast_g);
M_gray = im2gray(M);
M_bw = im2bw(abs(1-M_gray),0.78);
piattos(2,i) = sum(M_bw);
end

//for pillows
pillow = [];
for i = 1:4;
M = imread(strcat('D:\ap186\september17\pillow'+string(i)+'.jpg'));
r = M(:,:,1);
g = M(:,:,2);
contrast_r = (max(r)-min(r))/(max(r)+min(r));
contrast_g = (max(g)-min(g))/(max(g)+min(g));
pillow(1,i) = sqrt(contrast_r*contrast_r + contrast_g*contrast_g);
M_gray = im2gray(M);
M_bw = im2bw(abs(1-M_gray),0.6);
pillow(2,i) = sum(M_bw);
end

m = [];
m(:,1) = sum(piattos,'c')/size(piattos,2);
m(:,2) = sum(pillow,'c')/size(pillow,2);

//for test samples
test = [];
for i = 1:8
M = imread(strcat('D:\ap186\september17\test\'+string(i)+'.jpg'));
r = M(:,:,1);
g = M(:,:,2);
contrast_r = (max(r)-min(r))/(max(r)+min(r));
contrast_g = (max(g)-min(g))/(max(g)+min(g));
test(1,i) = sqrt(contrast_r*contrast_r + contrast_g*contrast_g);
M_gray = im2gray(M);
M_bw = im2bw(abs(1-M_gray),0.63);
test(2,i) = sum(M_bw);
end

//Minimum distance classification
d = [];
for i = 1:8
for j = 1:2
d(j,i) = test(:,i)'*m(:,j) - m(:,j)'*m(:,j);
end
end

d = abs(d);
x =[];
for i = 1:8
x(i) = find(d==min(d(:,i)));
end
for i = 1:length(x);
x(i) = x(i) - 2*(i-1);
end

x

SELF-GRADE: 10/10 because I have 100% accuracy

A16 - Image Color Segmentation

September 2, 2008

Suppose I have a reference image Mpatch and an image I want to segment M, they have NCC (normalized chromaticidty coordinates) components rpatch, gpatch and bpatch, and r, g, and b

Code for non-parametric image segmentation by color:


delta_pix = 0.01;
val=[];

num_=[];

num=[];

counter=1;

for i=0:delta_pix:1;

[x,y]=find((rpatch>=i) & (rpatch<(i+delta_pix)));

val(counter)=i;

num_(counter)=length(x);

counter=counter+1;

end

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

plot(val, num_);

B = zeros(size(M,1),size(M,2),3);
delta_pix = 0.01;

val=[];

num_=[];

num=[];

counter=1;

for i=0:delta_pix:1;

[x,y]=find((gpatch>=i) & (gpatch<(i+delta_pix)));

val(counter)=i;

num_(counter)=length(x);

counter=counter+1;

end

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

plot(val, num_);

for i = 0:delta_pix:1

[x,y] = find((g>=i) & (g<(i+delta_pix)));

for k=1:length(x) for l=1:length(y)

B(x(k),y(l),2) = num_((i+delta_pix)/delta_pix);

end

end

i

end
B(:,:,3) = 1-B(:,:,1)-B(:,:,2);

Code for parametric segmentation:

delta_pix = 0.01;
val=[];

num_=[];

PDF = [];

counter=1;

for i=0:delta_pix:1;

[x,y]=find((rpatch>=i) & (rpatch<(i+delta_pix)));

val(counter)=i;

num_(counter)=length(x);

counter=counter+1;

end

PDF = num_./((size(rpatch,1))*size(rpatch,2));

plot(val, PDF);

mn = find(PDF==max(PDF));

mu_red = val(mn);

sigma_red = stdev(PDF);
val=[];

num_=[];

PDF = [];

counter=1;

for i=0:delta_pix:1;

[x,y]=find((gpatch>=i) & (gpatch<(i+delta_pix)));

val(counter)=i;

num_(counter)=length(x);

counter=counter+1;

end

PDF = num_./((size(gpatch,1))*size(gpatch,2));

plot(val,PDF);

mn = find(PDF==max(PDF));

mu_green = val(mn);

sigma_green = stdev(PDF);
x = [0:delta_pix:1];

pr = (exp(-((x-mu_red)^2)/(2*sigma_red)))/(sigma_red*sqrt(2*%pi));

pg = (exp(-((x-mu_green)^2)/(2*sigma_green)))/(sigma_green*sqrt(2*%pi));

K = zeros(size(M,1),size(M,2));

for i = 0:delta_pix:1

[x,y] = find((r>=i) & (r<(i+delta_pix)));

for k=1:length(x)

for l=1:length(y)

green = round(g(x(k),y(k)));

probability = pr((i+delta_pix)/delta_pix)*pg((green+delta_pix)/delta_pix);

K(x(k),y(l)) = probability;

end

end

i

end

Display

G = im2gray(B);

L = im2gray(M);

subplot(311)

imshow(L,[]);

subplot(312)

imshow(G,[]);

Figure 1. The raw image.

Figure 2. The reference patch which belongs to the hand.

Figure 3. The raw image (top) and the maps showing the probability that a pixel belongs to the region of interest which in this case is the skin. The center image is the result of non-parametric image segmentation while the bottom image is the result of parametric image-segmentation.

Self-grade: 10/10 because I was able to finish the activity on time and my parametric and non-parametric image segmentations result looks reasonable

Thank you, Lei for helping me with histogramming.

A15 - Color Camera Processing

August 28, 2008


We performed white balancing on images under different white balancing settings. In figure 1, row 1, we see the original images under white balancing settings inside-cloudy, inside-incandescent and auto-white balance, respectively. We see that the colors they contain per point are different and the white of the first two settings do not appear white. We try to process the first two images using Reference White Algorithm (results in second column) and Gray World Algorithm (results in third column). Now their white looks whiter, with the processed image for the inside-cloudy setting closer to the AWB setting. We may say that processed images from the Gray World Algorithm are slightly darker than those from the Reference White Algorithm.Figure 1. Myriad-colored objects processed using Reference White Algorithm and Gray World Algorithm. The white patches used for processing are in the lower-right corner of the collage figure.

Next we, try processing objects with different shades of the same hue. We try this for a high contrast image (with dark and bright background), and low contrast image (bright background only). We now see that Gray World Algorithm results are brighter. And white appears whiter for the high contrast image.Figure 2. Monochromatic Objects Processed. 1st row: raw images; 2nd row: Reference White Algorithm outcomes; 3rd row: Gray World Algorithm results. Right most image: the white patch

Code for performing the Reference White and Gray World Algorithms

//1. Defining the stacksize and calling of image
n = 100000000;
stacksize(n);
M = imread('C:\Documents and Settings\AP186user20\Desktop\act15\i - incandescent.jpg');

//2. Reference White Algo
white = imread('C:\Documents and Settings\AP186user20\Desktop\act15\whitepatch_i-incandescent.jpg');
Rw = mean(white(:,:,1));
Gw = mean(white(:,:,2));
Bw = mean(white(:,:,3));

//3. Gray World Algo
K(:,:,1) = M(:,:,1)/Rw;
K(:,:,2) = M(:,:,2)/Gw;
K(:,:,3) = M(:,:,3)/Bw;
K = K/max(max(max(K)));
imwrite(K,'C:\Documents and Settings\AP186user20\Desktop\act15\ReferenceWhite_i-incandescent.jpg');

SELF-GRADE:
10/10 because I was able to finish the activity on time

CREDITS:
Rica, Aiyin, Beth for helping me acquire images and explaining the algorithms

A14 -Stereometry

I have acquired the images. Still processing...

Above are two images of the same object at slightly different perspectives (The object was moved to the right by 15cm).

My code:

f = 6; //focus

b = 150;

M1 = imread('C:\Documents and Settings\AP186user20\Desktop\august12_2008\Cube1.jpg');

M2 = imread('C:\Documents and Settings\AP186user20\Desktop\august12_2008\Cube2.jpg');

M1 = im2gray(M1);

M2 = im2gray(M2);

imshow(M1,[]);

M2 = im2gray(M2);

imshow(M1,[]);

x1 = locate(10,1);

imshow(M2,[]);

x2 = locate(10,1);

z = b*f/(x2(1,:) - x1(1,:));
nx = 40;

ny = 40;

x = linspace(1,1,nx);

y = linspace(1,1,ny);

[XP,YP] = ndgrid(x,y);
ZP = interp2d(XP, YP, x1(1,:), x1(2,:), splin2d(x1(1,:), x1(2,:), z, "natural"));
//plot(x2(1,:),x2(:,1), z);

mesh(ZP)

SELF-GRADE: 6/10, hindi ko pa yari

Thanks to Beth for pictures, Rica and Aiyin for storage

A12 - Geometric Distortion Correction

August 19, 2008

We have an image of a Capiz window, an analog to a grid, whose left side is slightly tilted. We aim to correct this distortion by creating an ideal grid, then mapping into it the grayscale values of the corresponding pixels in the distorted image.


First, we choose an eminently distorted area in the image. This area is covered by a rectangle in figure 1. We designate the rectangle to be the ideal polygon, how the shape of the set of small squares in the window appear. The ideal polygon vertices are labeled as [xi, yi] , while the vertices of the distorted area are tagged as [x^i, y^i].

The transformation of [x, y] to [x^, y^] can be expressed by a function r and s.

x^ = r(x,y)
y ^ = s(x,y) (1)

If the polygon is not too large, we could assume that the transformations for both directions happen linearly. Thus, we guess.

x^ = c1x + c2y + c3xy + c4
y^ = c5x + c6y + c7xy + c8 (2)

From equation (2), we could obtain the transformation vectors C14 = [c1, c2, c3, c4]t and
C58 = [c5, c6, c7, c8] t, where t stands for transposed via

C14 = (inv(T) ) X
C58 = (inv(T) ) Y

where inv means inv(T) means the inverse of T, T is a 4 by 4 matrix with row i given by [xi, yi, xiyi, 1], X = [[x^1,x^2, x^3, x^4] t and Y = [y^1, y^2, y^3, y^4] t.

In the ideal rectangle, we find the corresponding pixel in the distorted image. If the corresponding coordinates are integer-valued, we simply duplicate the point intensity value into the ideal space. Otherwise, we employ the bilinear interpolation wherein if x^ or y^ is not an element of the set of integers, the value at (x^, y^) is given by

v(x^,y^) = ax^ + by^ + cx^y^ + d (3)

To find the unknowns a, b, c and d, we use the nearest neighbors of (x^,y^).

We, see that the image has been quite corrected. The sides of the squares without using bilinear interpolation are crooked. Meanwhile, the sides of the squares when bilinear interpolation was employed were straight, though the resulting image had less contrast,

//1 Call the image
M = imread('D:\ap187\distorted_capiz_window.jpg');
M = im2gray(M);
imshow(M,[]);
[m,k] = size(M);

//2 Select large subgrid vertices
g = locate(4,1);
g = [46.723164 238.70056;
99.265537 239.83051;
100.96045 41.525424;
50.112994 40.960452 ];
gprime(:,1) = g(:,2);
gprime(:,2) = g(:,1);
gprime(:,1) = abs(gprime(:,1) - m - 1);
//3 Generate ideal grid
x_1 = gprime(1,1);
x_2 = gprime(1,1);
x_3 = gprime(3,1);
x_4 = x_3;
y_1 = gprime(1,2);
y_2 = gprime(2,2);
y_3 = gprime(2,2);
y_4 = gprime(1,2);

//4 Obtain transformation matrix C
T = [x_1 y_1 x_1*y_1 1;
x_2 y_2 x_2*y_2 1;
x_3 y_3 x_3*y_3 1;
x_4 y_4 x_4*y_4 1];

X = gprime(:,1);
Y = gprime(:,2);

C14 = (inv(T))*X;
C58 = (inv(T))*Y;

//5 Map the distorted image into the ideal space (w/ bilinear interpolation)
v = zeros(m,k);

for x = 5:m-5;
for y = 5:k-5;
t = [x y x*y 1];
xhat = t*C14;
yhat = t*C58;
xhat_integer = int(xhat);
yhat_integer = int(yhat);

if xhat_integer == xhat & yhat_integer == yhat then
if xhat_integer == 0 then
xhat_integer = 1;
end
if yhat_integer == 0 then
yhat_integer = 1;
end

v(x,y) = M(xhat_integer, yhat_integer);

else
xplus = xhat_integer + 1;
xminus = xhat_integer;
yplus = yhat_integer + 1;
yminus = yhat_integer;

nearestneighbors = [xplus yplus xplus*yplus 1;
xplus yminus xplus*yminus 1;
xminus yminus xminus*yminus 1;
xminus yplus xminus*yminus 1];

vhat = [M(xplus,yplus); M(xplus,yminus); M(xminus,yminus); M(xminus,yplus)];
a_b_c_d = inv(nearestneighbors)*vhat;
nu = [x y y*x 1];
v(x,y) = nu*a_b_c_d;
end
end
end

//6 Mapping without bilinear interpolation
v2 = zeros(m,k);

for x = 5:m-5;
for y = 5:k-5;
t = [x y x*y 1];
xhat = t*C14;
yhat = t*C58;
xhat_integer = int(xhat);
yhat_integer = int(yhat);
v2(x,y) = M(xhat_integer, yhat_integer);
end
end

//7 Showing the images, M - original, v - corrected with interpolation, v2 - corrected without interpolation
subplot(221)
imshow(M,[]);
subplot(222)
imshow(v2,[]);
subplot(223)
imshow(v3,[]);

I give myself 7/10 because it took me this long to solve the problem.

Thank you to Mark Leo for some explanations.

A13 - Photometric Stereo

August 7, 2008
We have 4 synthetic images, namely I1, I2, I3 and I4, of a sphere illuminated from 4 different directions.
The illumination matrix is a composite of 4 vectors given by

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].

The illumination is from an infinitely far away source. The relationship between the reflected intensity I, the illumination vector L and the object surface normal N is given by

I = L*N (1).

Since we know I and L, we could get the surface normals via

g = [(inv(V'V))V'] I (2).

We normalize g by dividing by its magnitude to obtain the normalized surface normal

n^ = g/|g| (3).

This components of the normal vector can be transformed to the partial derivatives by using the relation

df/dx = -nx/nz and df/dy = -ny/nz.

The derivatives are then integrated via line integration to obtain the surface profile of the hemisphere (See Figure )2. The drawback of using line integration is that it produces small spikes in the reconstruction.

code:
loadmatfile('C:\Documents and Settings\AP186user20\Desktop\august7_2008\photos.mat',['I1','I2','I3','I4']);
subplot(221)
imshow(I1,[]);
subplot(222)
imshow(I2,[]);
subplot(223)
imshow(I3,[]);
subplot(224)
imshow(I4,[]);
[m,k] = size(I1)
I = [I1(:)'; I2(:)'; I3(:)';I4(:)'];
V = [0.085832 0.17365 0.98106;
0.085832 -0.17365 0.98106
0.17365 0 0.98481
0.16318 -0.34202 0.92542];
MoorePennrose = (inv(V'*V))*V';
g = MoorePennrose*I;
[M,K] = size(g);
for i = 1:M;
for j = 1:K;
absolutevalue = sqrt( g(1,j)*g(1,j) + g(2,j)*g(2,j) + g(3,j)*g(3,j) );
n(i,j) = g(i,j)./(absolutevalue + 1.0e-015);
end
end
//derivative along the x-axis
p = -n(1,:)./(n(3,:)+1.0e-20);
//derivative along the y-axis
q = -n(2,:)./(n(3,:)+1.0e-20);

Imagexderivative = matrix(p,m,k);
Imageyderivative = matrix(q,m,k);
x=cumsum(Imagexderivative,2);
y= cumsum(Imageyderivative,1);
Image = x+y;
mesh(Image);
Reference
[1]. Dr. Soriano's lecture notes

Self-grade: 10/10

A11 - Camera Calibration

July 31, 2008



In this activity, we try to align the world coordinate system, which pertains to the real three-dimensional world of an object, with the two-dimensional image coordinate system captured by a camera.


We get the transformation matrix that would translate all world coordinates to the image coordinates. This transformation matrix depends on the angle the optic axis makes with respect to the object to be captured and on the intrinsic properties of the camera which are mainly affected by its focus. Denoted by A, it can be found by first realizing that













(Equation taken from Dr. Maricor Soriano's lecture notes.)


where xo, yo, zo are real world coordinates and yi, zi are image coordinates. Thus, A is given by

A = [ (inv(Q'Q)) Q'] d

where Q is the first matrix given in equation 1 and d is the column vector at the the right hand side of the equation. When there are N sample points, Q is a 2N x 11 matrix while d has a dimension of 2N x 1.

After we have determined A, we tried to predict the image coordinates of world points that were not previously sampled. The whole process is summarized by figure 1. A celfone camera was used to capture a checkerboard with two arms positioned perpendicular to each other. Each cell in the checkerboard is 1 inch x 1 inch. Twenty-five random points and their 2D projection on the image plane were used to determine the transformation matrix. Afterwards, the position in the image of five other points were interpolated using the transformation matrix. With the z plane in the real world chosen to be vertical, the left arm was designated to be the xz plane. The right arm is then the yz plane.

We could see that the predictions are accurate since some points even blocked the labeling numbers. Prediction #2, which missed the theoretical position, is not too far from where it should be.
























Self-grade: 10/10
Credits: Benj for some clarifications on the later part of the problem


//Camera Calibration
//1 Image space
d = [135 138;230 165;177 188;229 259;178 91;279 141;261 216;110 213;244 89;159 260;179 233;135 89;315 137;36 137;63 266;64 58;315 311;63 322;317 21;39 32;38 242;180 47;336 227;88 189;297 53];
//2 Object space
s = [3.5 0 1.5;0 1.5 0.5;1.5 0 -0.5;0 1.5 -3.5;1.5 0 3.5;0 4.5 1.5;0 3.5 -1.5;4.5 0 -1.5;0 2.5 3.5;2.5 0 -3.5;1.5 0 -2.5;3.5 0 3.5;0 6.5 1.5;7.5 0 1.5;6.5 0 -3.5;6.5 0 4.5;0 6.5 -4.5;6.5 0 -5.5;0 6.5 5.5;7.5 0 5.5;7.5 0 -2.5;1.5 0 5.5;0 7.5 -1.5;5.5 0 -0.5;0 5.5 4.5];
//3 Finding the transformation matrix a
Q = zeros(50, 11);
i = 1;
j = 2*i-1;
while i<>
Q(j,1:3) = s(i,1:3);
Q(j,4) = 1;
Q(j,5:8) = 0;
Q(j,9:11) = -d(i,1)*s(i,1:3);
Q(j+1,1:4) = 0;
Q(j+1,5:7) = s(i,1:3);
Q(j+1,8) = 1;
Q(j+1,9:11) = -d(i,2)*s(i,1:3);
i = i+1;
j = 2*i-1;
end;
d = matrix(d',50,1);
a = inv(Q'*Q)*Q'*d;
//Predicting the image coordinates of unsampled world points
test = [3 0 0 1;0 5 -3 1;0 3 2 1;1 0 2 1;3 0 -2 1];
newd = a*test';
newd = matrix(newd,5,3);
newd = newd';
[m,k] = size(M);
overlay = zeros(m,k);
for i = 1:size(newd,2);
overlay(round(newd(2,i)),round(newd(1,i))) = 1;
end
M = imread('D:\ap186\july29_2008\cropped.jpg');
M = im2gray(M);
subplot(121)
imshow(M,[]);
subplot(122)
imshow(overlay);

A11 - Preprocessing Handwritten Text

First, we get chose a portion of the grayscaled image. Then we, removed the vertical lines by knocking-out the frequencies in the Fourier transform of the image that correspond to vertical lines. These are those frequencies that lie along the vertical axis of the shifted FFT. So the mask used was a black vertical strip with a hole in the middle so as not to remove the DC component of the image.

Next, I binarized the image. I used thresholds = 80, 90 and 100 intensity values and found out that threshold = 80 is the best. So I, tried applying the operations closing of holes and opening of holes for the image and also opening the holes for holes-closed image. I think that the image became cleaner after performing closing of holes only. Next, I thinned the characters/regions into one pixel. The final image is shown in figure 6. While not all letters are recognizable, the letters U,P and S have been, in my opinion, enhanced.

Enjoy mag-label ng continuous regions.

Code:

//1 Calls the image and converts it to grayscale

M = imread('D:\ap186\july20_2008\toEnhance3.jpg');

M = im2gray(M);

M = M-min(min(M));

M = 255*M/max(max(M));

//3 Shows the FFT of the best binarized image,K2

FM = fft2(M);

FMshift = fftshift(abs(FM));

//4 Creation of a mask to filter the horizontal and vertical lines in the image

[m,n] = size(M);mask = ones(m,n);

mask(:,n/2-5:n/2+5) = 0;mask(m/2-3:m/2+3,n/2-3:n/2+3) = 1;

//5 Shows the image FFT and the mask to be used

subplot(121)

imshow(log(FMshift),[]);

subplot(122)

imshow(mask);

//6 Filtering of the lines; the new image is called H

FHM = FM.*fftshift(mask);

H = abs(ifft(FHM));

//7 Shows the original image, the filtered image, and the latter's FFT

subplot(131)

imshow(M,[]);

subplot(132)

imshow(H,[]);

subplot(133);

imshow(fftshift(abs(FHM)),[]);
H = H-min(min(H));

H = 255*H/max(max(H));

subplot(221)

histplot([0:255],H);

subplot(222)

L1 = im2bw(H,90/255);

L1 = 1-abs(L1);

imshow(L1);

subplot(223);

L2 = im2bw(H,100/255);

L2 = 1-abs(L2);

imshow(L2);

subplot(224);

L3 = im2bw(H,110/255);

L3 = 1-abs(L3);

imshow(L3);

//Close or Open holes for L1
L = ones(2,2);

zero = zeros(2,2);

negL = abs(zero-L);

D = dilate(L1,negL);

Close = erode(D,negL);
E = erode(Close,L);

Open = dilate(E,L);

subplot(221)

imshow(L2);

subplot(222)

imshow(Close);

subplot(223)

imshow(Open);

subplot(224)

thn = thin(Close);

imshow(thn);
lbl = bwlabel(thn);

imshow(lbl,jetcolormap(max(max(lbl))));

Grade: 8/10

processed image is not so readable

A10 - Morphological Operations

We tried to get the best estimate of the area of a punched paper by statistical analysis of a set of punched papers. We cut the image into six sub-images such as the one in figure 1 and converted them into grayscale. We looked at the histogram of their intensity values and determined that the threshold value should be set at intensity = 200.

We then performed the Morphological operations of opening and closing of holes to see what would happen. Opening of holes means to dilate with a structuring element L the image eroded with L. Closing of holes means to erode with the reflection of L the image dilated by the reflection of L. It seems that closing of holes had a better effect of producing a cleaner binary image. So I chose to continue processing this.

As shown in figure 3, the closed regions were labeled with a specific number. The area of each labeled region is determined by pixel counting.




























The histogram of the above sample is merged with the histograms of other subimages. The result is the histogram shown above. It seemed that the histogram needed no further noise eradication. The best estimate for the area was 549 px.

To check, we get the number of included pixels in a binarized image of a single punched paper. It was 519. The deviation of the best estimate from this value was around 6%.














Code:
for i = 1:7;
//1 Call image, transform it to black-and-white
M = imread(strcat(string(i)+'.jpg'));
M = im2gray(M);M = M*255;
//subplot(121)
//imshow(M,[]);
//subplot(122)
//histplot([0:255],M);
threshold = 200;
K = im2bw(M,threshold/255);
//2 Close-openL = ones(5,5);
zero = zeros(5,5);
negL = abs(zero-L);
E = erode(K,L);
D = dilate(K,negL);
Open = dilate(E,L);
Close = erode(D,negL);
subplot(131)
imshow(K);
subplot(132)
imshow(Close);
subplot(133)
imshow(Open);
//3 Label Opened image because it looks good
lbl = bwlabel(Open);
N = max(max(lbl));
//imshow(lbl,jetcolormap(N));
//4 Count area of each similarly labeled region
Area = zeros(7,N);
for j = 1:N; [x,y] = find(lbl == i);
Area(i,j) = length(x);
end
end
histplot([1:max(max(Area))],Area);
mn = mean(mean(Area))
std = stdev(Area)
//5 Eradicate noiseupperbound = mn + std + 1;
lowerbound = mn - std - 1;
Area = Area.*(Arealowerbound);
histplot([1:max(area)],Area);

Credits:

Lei Uy and Beth Prieto for helpful discussions, Benj Palmares for help in histogramming.

Self-evaluation: 9/10

I fell happy when labelling the regions because the resulting image is colorful

A9 - Binary Operations

July 17, 2008

The goal of this activity is to perform area estimation of cells represented by punched paper.

In figure 1, part of the entire image is binarized using a threshold pixel value that is equal to 200 where the valley in the histogram occurs.
















Next, we performed higher order morphological operations on the binarized image: opening of holes and closing of holes. We let M be the original image, K the binary image and L the structuring element. The two operations are implemented in SciLab as

M = imread('C:\Documents and Settings\AP186user20\Desktop\Circles1.jpg');
M = M - 1;
histplot([0:255,M);
threshold = 200;
K = im2bw(M,threshold/255);
L = ones(5,5);
zero = zeros(5,5);
negL = abs(zero-L);
E = erode(K,L);
D = dilate(K,negL);
Open = dilate(E,L);
Close = erode(D,negL);





















Reference:
[1] http://www.ph.tn.tudelft.nl/Courses/FIP/noframes/fip-Morpholo.html#Heading98

A8 Morphological Operations

July 15, 2008


In this activity, we tried to investigate the effects of dilation and erosion of an image by a structuring element. We predicted first what would happen and then simulated them in SciLab.
The structuring elements used were:

• 4x4 ones
• 2x4 ones
• 4x2 ones
• cross, 1pixel thick, 5 pixels long

We dilated and eroded images of a square, a triangle, a hollow square (60 x 60, 4px thick), a circle and a cross.

So far, most of my predictions were wrong. I was fairly good in guessing for the square, triangle and hollow square images and slightly for the cross. I was terrible with the circle.

I also realized that the image size I should should be larger than the shape of the actual square or triangle, etc. This is because if the input image Img is m x n, the dilated image is also m x n. Some parts were cut in the dilated image. But since in erosion the result is smaller than the image, the effect of this operation was more imminent.

Remarks:

Beth Prieto, Ayeen Laganapan, Rica Mercado for helping me grapple with concepts.

Julie Dado for assisting me in publishing a video.

Grade: 9/10

A. Square

B. Triangle

C. Hollow square

D. Circle

E. Cross

A7 - Enhancement in the Frequency Domain

July 10, 2008


A. Anamorphic Property of the Fourier Transform

We tried to get the Fourier transform of a sinusoidal object similar to a corrugated roof. The first sinusoid has a frequency of 2. Its Fourier transform is shown as two dots above and below the center of the FFT and symmetric about the origin. We explored what would happen if we increase the frequency to 4. The two spots move further away from each other. We rotated the sinusoid. The FFT was rotated in the opposite direction. Next, we superpose a horizontal and a vertical sinusoid, both with frequency = 2, by simply adding them. The resulting FFT is the superposition of the FFTs of each sinusoid, since we already know that the FFT of the vertical sinusoid is the same as that of the horizontal sinusoid only that it is rotated 90 degrees.





















































B. Fingerprints: Ridge Enhancement

Here, we try to enhance the picture of a fingerprint in such a way as to make the ridges more visible. We created an aperture shown below which we multiplied by the Fourier transform of the fingerprint image.











We made sure that the position of the mask corresponds to that of the secondary bright spots in the fingerprint image's FFT.

Code:
mask = imread('square6.bmp');
finger = imread('fingerprint2.jpg');
maskgray = im2gray(mask);fingergray = im2gray(finger);
FFT_finger_mask = fft2(fingergray).*fftshift(maskgray);
Image_filtered = abs(ifft(fftshift(FFT_finger_mask)));
Image_filtered = Image_filtered-mean(mean(Image_filtered));
scf,imshow(Image_filtered,[]);
scf, imshow(fingergray,[]);
scf,imshow(log(fftshift(abs(fft2(fingergray)))),[]);
scf,imshow(log(fftshift(abs(FFT_finger_mask))),[]);



C. Lunar Landing Scanned Pictures: Line Removal

Here, we try to remove the vertical lines regarded as artifacts in the image. The original photo is taken from http://www.lpi.usra.edu/lunar/missions/apollo/apollo_11/images. We get the FFT transform of the image. If it fringe-free, we expect that it should only possess a DC component or a zero-order frequency indicated by a very bright central spot. However, if this is not the case, we expect the spots that are less bright than the DC would appear. For the original image with vertical lines, the FFT had horizontal slits to the left and right of the DC component. We therefore, make a mask that would cover these slits and multiply them to the FFT of the original image. We see in the final image, that the vertical lines reduced in visibility.

Code:

M = imread('hi_res_vertical_lg.gif');

K = imread('slit2.bmp');

K = im2gray(K);

FFT_M = fft2(M);

FFT_MK = FFT_M.*fftshift(K);

imshow(log(fftshift(abs(FFT_MK))),[]);

MK = (ifft(FFT_MK));

MKimage = abs(MKimage)/max(max(abs(MKimage)));

scf,imshow(abs(MK),[]);

scf,imshow(log(fftshift(abs(FFT_M))),[]);

scf,imshow(M,[]);

SELF-EVALUATION: 10/10.

Fun activity, practical but tedious.

Credits:

Ma'am Jing, Andrew Banas, Gerold Pedemonte and Beth Prieto for answering some of my questions.