Cole's AP186

Activity 20: Neural Networks

2008-09-26T15:55:00.000-07:00

In this activity, we use neural networks as another method to classify our objects into groups depending on the values of the extracted features. Similar to activity 19, we use 2 classes: kwekkwek and pillows. However, only 2 features are used to differentiate the 2 groups for simplicity. These are the length to height ratio and the red component of color.

A neural network is a computational model of how neurons in the brain work. It is a preferred alternative to linear discriminant analysis in the pattern recognition because one does not need heuristics and recognition rules to perform classification. Instead neural networks “learn” the rules of a mapping by example. Although it may take a long time for a network to train, its recognition processing speed is faster once it has learned.

source: M. Soriano, A20 - Neural Networks.pdf
____________________________________________________________________

A code was already prepared by J. Tugaff and it only needs to be modified depending on the number and values of the input features. For our purposes, we have:

2 input features
4 training objects from each class (kwekkwek and pillows)
4 test objects from each class

We first train our neural network using the training set.

Length/Height, Red component --(kwek kwek; pillows)
0.97, 0.74
1.00, 0.77
0.96, 0.73
0.91, 0.74
0.96, 0.30
1.00, 0.30
1.02, 0.33
1.02, 0.33
____________________________________________________________________

//original code c/o Jeric Tugaff

//ensure the same starting point each time
rand('seed',0);

//network def.
//neurons per layer, including input
//2 neurons in the input layer, 2 in the hidden layer and 1 in the output layer
N = [2,2,1];

//inputs, 1st column = length/height, 2nd column = red component
//training set
x = [0.97, 0.74;
1.00, 0.77;
0.96, 0.73;
0.91, 0.74;
0.96, 0.30;
1.00, 0.30;
1.02, 0.33;
1.02, 0.33]';

//targets, 0 if kwekkwek and 1 if pillows
t = [0 0 0 0 1 1 1 1];
//learning rate is 0.1 and 0 is the threshold for the error tolerated by the network
lp = [0.1,0];

W = ann_FF_init(N);

//training cycles
T = 1000;
W = ann_FF_Std_online(x,t,N,W,lp,T);
//x is the training, t is the output, W is the initialized weights,
//N is the NN architecture, lp is the learning rate and T is the number of iterations

//full run
ann_FF_run(x,N,W)
//the network N was tested using x as the test set, and W as the weights of the connections
____________________________________________________________________

We expect that this code will output values close to (0, 0, 0, 0, 1, 1, 1, 1), meaning that the first 4 objects will be classified as members of the 'kwekkwek' class while the last 4 objects will be classified as members of the 'pillows' class. The actual results are (0.1448030, 0.1291128, 0.1506565, 0.1342528, 0.8876302, 0.8956291, 0.8738883, 0.8738883). We can now say that our neural network has "learned" the features of each class and has sorted the objects into the correct groups. We now use the test set to see whether the network can classify the other objects accurately as well.

Length/Height, Red component --(kwek kwek; pillows)
0.67, 0.77
0.81, 0.79
1.08, 0.78
1.00, 0.81
1.09, 0.35
1.12, 0.35
1.05, 0.35
0.94, 0.32

Using the training parameters, the resulting values are (0.0900170, 0.0954682, 0.1355222, 0.1071037, 0.8711712, 0.8782258, 0.8610966, 0.8642343). A 100% correct classification was achieved. Next, we try to tweak the parameters (learning rate and training cycles) to see whether the values will go closer to 0 and 1 or whether the clasification will still remain 100% correct.

learning rate = from 0.1 --> change to 1.0
training cycles = from 1000 --> change to 400
result: (0.0201587, 0.0220962, 0.0384002, 0.0265042, 0.9631881, 0.9668789, 0.9575001, 0.9582095)

The values are closer to 0 and 1 using the new training parameters and the classification is still 100% correct. Now, let us try to change the order of the test set to make sure that the 100% correct classification does not depend on the order.

Length/Height, Red component --(kwek kwek; pillows)
0.67, 0.77
1.05, 0.35
1.08, 0.78
1.09, 0.35
0.81, 0.79
1.12, 0.35
0.94, 0.32
1.00, 0.81

learning rate = 1.0
training cycles = 400
result: (0.0201587, 0.9575001, 0.0384002, 0.9631881, 0.0220962, 0.9668789, 0.9582095, 0.0265042)

The classification is still 100% correct and we have shown that the order of the objects in the input matrix is not important.
____________________________________________________________________

I give myself 10 points for this activity because of the 100% correct classification of the objects into their respective classes. Thanks to Jeric Tugaff for the ANN toolbox and the code and Mark Leo Bejemino for helping me to apply the code to my objects and extracted features.

Activity 19: Probabilistic Classification

2008-09-18T11:13:00.000-07:00

In this activity, we apply a different method to classify objects into groups. In the previous activity, mean distance classification was used. Here, we employ probabilistic classification via Linear Discriminant Analysis (LDA). A detailed discussion on the method and equations used can be viewed from Pattern_Recognition_2.pdf by Dr. S. Marcos. We use two classes of objects from the previous activity and see whether this method classifies them correctly or not.

Kwekkwek and Pillows

Also shown below are the values for length/height, R, G, B for each object.Again, we separate the 8 objects per class into a test set and a training set. The training set is given below. From here the method is implemented and the code and results are also shown below.//act 19 code
//this code just follows the procedure described in Pattern_Recognition_2.pdf
x = [0.97 0.74 0.37 0.03;
1.00 0.77 0.42 0.03;
0.96 0.73 0.38 0.03;
0.91 0.74 0.38 0.06;
0.96 0.30 0.16 0.07;
1.00 0.30 0.16 0.07;
1.02 0.33 0.18 0.08;
1.02 0.33 0.19 0.09];

y = [1 1 1 1 2 2 2 2]';
g = 2;
x1 = x(1:4,:);
x2 = x(5:8,:);

u1 = [mean(x1(:,1)) mean(x1(:,2)) mean(x1(:,3)) mean(x1(:,4))];
u2 = [mean(x2(:,1)) mean(x2(:,2)) mean(x2(:,3)) mean(x2(:,4))];
u = [mean(x(:,1)) mean(x(:,2)) mean(x(:,3)) mean(x(:,4))];

x1_o = [x1(1,:)-u; x1(2,:)-u; x1(3,:)-u; x1(4,:)-u];
x2_o = [x2(1,:)-u; x2(2,:)-u; x2(3,:)-u; x2(4,:)-u];

n1 = 4; n2 = n1;
c1 = (x1_o'*x1_o)/n1;
c2 = (x2_o'*x2_o)/n2;

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

Cinv = inv(C);
p = [4/8; 4/8];

for k = 1:8
f1(k) = u1*Cinv*x(k,:)' - 0.5*u1*Cinv*u1' + log(p(1));
f2(k) = u2*Cinv*x(k,:)' - 0.5*u2*Cinv*u2' + log(p(2));
end
____________________________________________________________________

RESULTS:
From the table above, we can see that this method resulted in a 100% correct classification for the training set. The higher value of f will determine which group the object is classified under. We implement this on the test set and see what happens.Again, a 100% correct classification is obtained. We have illustrated the success of this method in classifying objects into groups.
____________________________________________________________________

I give myself 10 points for this activity since the method was implemented correctly and there is 0% error.

Activity 18: Pattern Recognition

2008-09-16T23:25:00.000-07:00

In this activity, we try to test if a random object can be classified using pattern recognition. First, we obtain 8-10 pieces each of similar objects like leaves, fishballs, potato chips, etc. One category of objects is called a class. We assemble 8-10 objects that can be classified into 2-5 classes. Half of this will be used as a test set and the other as training set. In this case, I have assembled 8 pieces each of:

1. pillows (a chocolate-coated snack)
2. piatos potato chips
3. squidballs
4. kwek kwek (orange, flour-coated quail egg)

We extract certain features of these objects such as size, shape, and color and arrange these numbers into an organized set called a feature vector. For each individual object, the feature vector that I will use will have four dimensions namely:

1. ratio of length to height of object (x-axis/y-axis)
2. red component of color
3. green component of color
4. blue component of color

After finding these values for our 32 objects in total, we tabulate them in order to try to visualize clustering in 4-D feature space. We expect clustering to happen if objects are similar to each other. Thus, piatos chips should be clustered together with other piatos chips since they are similar in both size and color.

To calculate for the size, I employed the following algorithm. Assuming that the each object under one class is depicted in one image, we open it in Scilab, convert to grayscale, look at its histogram to find the proper threshold value in order to convert it to B&W. After that, we do morphological operations to clean the image (close holes and connect/remove stray pixels). Then we use the follow function to measure the contour of the object, get the maximum and minimum values in the x and y axes, and find their ratio. To calculate for the color, we need to first separate the region of interest from the background. Once that is finished, we get the mean for each color channel. These numbers constitute our feature vector for each object. There are 8 objects each for 4 classes used so there are 32 feature vectors.
____________________________________________________________________
Shown below are the values for the feature vectors. For each class, half will be used as a test set and the other half as a training set.____________________________________________________________________

If we have any other random object, we can say whether it belongs to the class of squidballs or pillows through a method called minimum distance classification.source: M. Soriano, A18- Pattern Recognition.pdf

We use the training feature vectors to find the class representatives mj. Then, we classify the test feature vectors using this method.
____________________________________________________________________
class representatives or mean feature vector, mj:
1. kwek kwek - [0.96, 0.75, 0.39, 0.04]
2. squidballs - [1.03, 0.55, 0.35, 0.09]
3. pillows - [1.00, 0.32, 0.17, 0.08]
4. piatos - [1.35, 0.59, 0.45, 0.16]From this table it can be seen that only 2 out of the 16 objects were classified incorrectly. Thus, the success rate of this method is 87.5%. Errors may have come from the calculation of the length to height ratio because the images might not have been thresholded and/or morphologically enhanced correctly.
____________________________________________________________________

I give myself 10 points for this activity since the method was implemented resulting in only a few errors. The success rate of this method in pattern recognition was found to be 87.5%, which is a fairly acceptable number. Thanks to JC Nadora, Benj Palmares, Raf Jaculbia, Mer Camba, Ed David, Billy Narag, and Jorge Presto for buying the objects/food tested (kwek kwek, squidballs, piatos, pillows).

Activity 17: Basic Video Processing

2008-09-09T11:00:00.000-07:00

Video, whether analog or digital, is a string of still images shown in rapid succession such that an impression of movement is perceived by an observer. Therefore, all image processing techniques we've learned so far may be applied on image frames of a video. The advantage of video is that it allows us access to another physical dimension, time. Thus, the dynamics and kinematics of a system can be extracted from video.

In this activity, we aim to extract some kinematic constants and variables like gravitational constant, position, and velocity upon taking a video of an event like dropping a ball or watching a spring reach equilibrium when initially stretched. These extracted variables will be plotted against time and the kinematic constants will be calculated. The results will be compared against analytical values.

In order to do this with Scilab, we need to install a video processing toolbox. However, the SIVP toolbox cannot read more than two frames from an avi video so we have to use other software to achieve our goal. These are Windows Movie Maker - for capturing the video and saving it directly on the PC, Stoik Video Converter - for converting the .wmv file generated by Movie Maker into an .avi file, and Virtual Dub - to parse the .avi file into image sequences in either .jpg or .bmp formats. Once stored as individual images, the position, shape or area can be extracted using image processing techniques. We loop through the image sequence to get the time variation of the variables of interest.

source: M. Soriano, A17 - Basic Video Processing.pdf

The kinematic event that I used was using a steel ruler as a physical pendulum. From this, I can calculate for some quantities like the center of oscillation, its position after some time, the period of oscillation, etc.
____________________________________________________________________

Physical Pendulum
http://en.wikipedia.org/wiki/Pendulum#Physical_pendulum

The simple pendulum assumes that the rod is massless and that the bob is small, so has negligible angular momentum in itself. A physical pendulum has significant size and mass, and hence a significant moment of inertia.

A physical pendulum behaves like a simple pendulum but the expression for the period is modified.

where:

I

is the moment of inertia of the pendulum about the pivot point

L

is the distance from the center of mass to the pivot point

m

is the mass of the pendulum

For a physical pendulum, since this has the same form as the simple pendulum, it is possible to define the center of oscillation. This is the distance from the pivot point, at which, if all the mass of the pendulum was concentrated (and thus forming a simple pendulum), would give the pendulum the same period.

The moment of inertia in this case is given by

with height h, width w, and mass m (Axis of rotation at the end of the plate)
http://en.wikipedia.org/wiki/List_of_moments_of_inertia
____________________________________________________________________

Before calculating the desired quantities, we first need to process our images. This activity utilizes most of the techniques we have learned so far. Step by step, these are:

1. Separate the video per frame using other software.
2. White balance the images.
3. Normalize the color.
4. Do color segmentation on the object (in our case, the ruler that acted as a physical pendulum).
5. Convert to a B&W image with the proper threshold value.
6. Do morphological operations to close holes and enhance the image.

After these, we can calculate for the centroid proceed in solving for the desired quantities numerically.
____________________________________________________________________

RESULTS
ORIGINAL, WHITE BALANCED, NORMALIZED, SEGMENTED, B&W, MORPHED IMAGES

measured quantities: m = 175 g, length = 0.63 m, width = 0.03 m

Analytic solutions:
Moment of inertia = 175g(0.63m*0.63m)/3 + 175g(0.03m*0.03m)/12 = 23.165625 g m^2
Center of oscillation = I/mL = 23.165625 m/(175*0.315) = 0.4202381 m
Period of oscillation = 2*pi*sqrt(23.165625/(175*9.8*0.315)) = 1.3011116 s

We now proceed in solving these quantities numerically from the data.

Length of ruler = 105 pixels --> 0.63 m
Width of ruler = 6 pixels
Moment of inertia = 175g(0.63m*0.63m)/3 + 175g((6*0.63m/105)*(6*0.63m/105))/12
= 23.1714 g m^2

%error = 0.02493

Center of oscillation = 23.1714 m/(175*(52.5*0.63/105)) = 0.4203429 m

%error = 0.02494

Period of oscillation:
Our video was captured with a fps of 15. The graph of the position vs time of the centroid of the ruler is shown below.
The wave is not a perfect sinusoid due to some air resistance present. The peaks and troughs are relatively equally spaced between each other so the period can be calculated. Peak to peak, the heights for one frame is given from the data:
94, 93, 90, 87, 83, 78, 73.5, 68.5, 65, 64.5, 64.5, 65, 68.5, 73, 77.5, 83.5, 87.5, 90, 92.5, 92.5
Thus, one period was captured in 20 frames. If 15 frames = 1 second, then using ratio and proportion, 20 frames = 1.3333333 s.

%error = 2.47647
____________________________________________________________________

chdir('act17images\');
img = 216; //total no. of frames in video

//open cropped image for color segmentation
im2 = imread('act17crop1.jpg');
wid = zeros(img,1); leng = zeros(1,1);

for i = 1:img
im = imread('images'+string(i)+'.jpg'); //open images
white = im(1:120,120:160,:);

//white balance
r = mean(white(:,:,1));
g = mean(white(:,:,2));
b = mean(white(:,:,3));
im(:,:,1) = im(:,:,1)/r;
im(:,:,2) = im(:,:,2)/g;
im(:,:,3) = im(:,:,3)/b;
im = im/max(im);
m = find(im>1);
im(m) = 1;

//color normalization
I = im(:,:,1)+im(:,:,2)+im(:,:,3);
im(:,:,1) = im(:,:,1)./I;
im(:,:,2) = im(:,:,2)./I;
im(:,:,3) = im(:,:,3)./I;

//histogram backprojection
R = linspace(0,1,32); G = R;
P = zeros(32,32);
[x2,y2] = size(im2);
for j = 1:x2
for k = 1:y2
xr = find(R <= im2(j,k,1));
xg = find(G <= im2(j,k,2));
P(xr(length(xr)),xg(length(xg))) = P(xr(length(xr)),xg(length(xg))) + 1;
end
end

P = P/sum(P);
[x1,y1] = size(im);
bpj = zeros(x1,y1);
for j = 1:x1
for k = 1:y1
xr = find(R <= im(j,k,1));
xg = find(G <= im(j,k,2));
bpj(j,k) = P(xr(length(xr)), xg(length(xg)));
end
end
bpj = bpj/max(bpj);

//convert to B&W
bpj = im2bw(bpj, 0.02);

//morphological operations
strel = ones(4,1);
bpj = erode(bpj,strel); //opening operation
bpj = dilate(bpj,strel);
bpj = dilate(bpj,strel); //closing operation
bpj = erode(bpj,strel);

//get centroid
[p,q] = find(bpj == 1);
wid(i) = (max(q)+min(q))/2;
leng(i) = (max(p)+min(p))/2;
end

//plot position vs time
plot([1:img]/15,wid); //consider fps of video
title('physical pendulum');
xlabel('time (seconds)');
ylabel('relative amplitude');
____________________________________________________________________

I give myself 10 points for this activity since the calculated quantities were close to their analytic values. My collaborators were Jeric Tugaff and Mark Leo Bejemino.

Activity 16: Color Image Segmentation

2008-09-02T18:14:00.000-07:00

In image segmentation, a region of interest (ROI) is picked out from the rest of the image such that further processing can be done on it. Selection rules are based on features unique to the ROI. One such feature is color. Color has been used to segment images of skin regions in face and hand recognition, land cover in remote sensing, and cells in microscopy.

In this activity, we aim to do image segmentation using 2 different methods: parametric and non-parametric probability distribution estimation. First, we convert image RGB to normalized chromaticity coordinates (NCC) and then apply the 2 techniques.

We can transform the image RGB values to NCC by,
Since r + b + g = 1, then b = 1 - r - g and it is enough to represent chromaticity by just two coordinates, r and g. When segmenting 3D objects it is better to first transform RGB into rgI. The figure below shows the r-g color space.When r and g are both zero, b = 1. Therefore, the origin corresponds to blue while points corresponding to r=1 and g=1 are points of pure red and green, respectively. Any color therefore can be represented as a point in NCC space. Thus the pixel chromaticities of a red ball will only occupy a small blob in the red region regardless of its shading variation.

source: M.Soriano, A16 - Color Image Segmentation.pdf

The histogram of a green chair is shown below.
____________________________________________________________________

Parametric method:
Original & Probability image

Here, the original object is a green chair. We cropped a region of interest (center part of the backrest) as our color model and tried to see where it matched that of the original image. The pixel values of the resulting image is the probability that a pixel will have the same color as the model. This method segmented the image well and we can infer that the bright parts of the seat are of the same color as that of the backrest.
____________________________________________________________________

Histogram backprojection (non-parametric) method:
Original & Probability image

Histogram backprojection is one technique where based on the color histogram, a pixel location is given a value equal to its histogram value in chromaticity space. We implemented this and the resulting image was generated. Upon comparison with the result from the parametric method, there are more white pixels that correspond to the cropped color model. We can really see where the similar colors lie. The segmentation is better using the 2nd method. This also has the advantage of faster processing because no more computations are needed, just a look-up of histogram values. In the parametric method, the mean and standard deviation of r and g must be found in order to calculate a Gaussian distribution that tests the probability of pixel membership to the ROI.
____________________________________________________________________

stacksize(4e7);

im1 = imread('green chair.jpg');
im2 = imread('green chair crop.jpg');
//imshow(im1);

R1 = im1(:,:,1); G1 = im1(:,:,2); B1 = im1(:,:,3);
I1 = R1 + G1 + B1;
R2 = im2(:,:,1); G2 = im2(:,:,2); B2 = im2(:,:,3);
I2 = R2 + G2 + B2;

r1 = R1./I1; g1 = G1./I1; b1 = B1./I1;
r2 = R2./I2; g2 = G2./I2; b2 = B2./I2;

//parametric estimation
ur = mean(r2); sr = stdev(r2);
ug = mean(g2); sg = stdev(g2);

pr = 1.0*exp(-((r1-ur).^2)/(2*sr^2))/(sr*sqrt(2*%pi));
pg = 1.0*exp(-((g1-ug).^2)/(2*sg^2))/(sg*sqrt(2*%pi));

prob = pr.*pg;
prob = prob/max(prob);
//imshow(prob,[]);

//2d histogram
r = linspace(0,1,32); r = g;
P = zeros(32,32);
[x,y] = size(im2);
for i = 1:x
for j = 1:y
xr = find(r <= im2(i,j,1));
xg = find(g <= im2(i,j,2));
P(xr(length(xr)), xg(length(xg))) = P(xr(length(xr)), xg(length(xg)))+1;
end
end
P = P/sum(P);
//surf(P);

//backprojection
[x,y] = size(im1);
recons = zeros(x,y);
for i = 1:x
for j = 1:y
xr = find(r <= im1(i,j,1));
xg = find(g <= im1(i,j,2));
recons(i,j) = P(xr(length(xr)), xg(length(xg)));
end
end
imshow(recons, []);
____________________________________________________________________

I give myself 10 points for this activity since the colored image was properly segmented using 2 different methods. Histogram backprojection is more efficient since less calculations are needed, only a look-up table. It also produced a better probability image. Thanks to Jeric Tugaff for his help in the histogram backprojection code.

Activity 15: Color Image Processing

2008-08-28T10:46:00.000-07:00

A colored digital image is an array of pixels each having red, green and blue light overlaid in various proportions. Per pixel, the color captured by a digital color camera is an integral of the product of the spectral power distribution of the incident light source S(l), the surface reflectance r(l), and the spectral sensitivity of the camera h(l). That is, if the color camera has spectral sensitivity hR(l), hG(l) , hB(l) for the red, green, and blue channels respectively, then the color of the pixel is given by

is a balancing constant equal to the inverse of the camera output when shown a white object (r(l) = 1.0). Equations (1) to (4) explain why sometimes the colors captured by a digital camera appear unsatisfactory.

In this activity, we aim to observe the differences between images obtained using several white balancing settings in a camera and compare two white balancing algorithms: Reference White and Gray World. We capture an obviously wrongly white balanced image, apply the two algorithms, and comment on their white balancing quality.

In the Reference White Algorithm, you capture an image using an unbalanced camera and use the RGB values of a known white object as the divisor. In the Gray World Algorithm, it is assumed that the average color of the world is gray. Gray is part of the family of white, as is black. Therefore, if you know the RGB of a gray object, it is essentially the RGB of white up to a constant factor. Thus, to get the balancing constants, you take the average red, green and blue value of the captured image and utilize them as the balancing constants.

source: M. Soriano, A15 - Color Image Processing.pdf
____________________________________________________________________

The following images were taken using an Olympus Stylus 770SW digital camera with an exposure of -2 under a fluorescent light source. I utilized a Macbeth chart since the major hues were represented and white was also present.

(order of images from left to right, top to bottom: auto WB, daylight, cloudy, tungsten, fluorescent1, fluorescent2, fluorescent3)

Upon comparison, we can observe that the colors are represented in high quality in image 1 and 5, which are the auto white balance and fluorescent1 settings, respectively. The sunlight and cloudy settings made the colors look dimmer, flourescent2 and 3 had a larger blue component, and the tungsten setting was just plain bad. What it did was since a tungsten light source has a low blue component upon looking at its spectra, the camera increased the blue to compensate. I will use this obviously wrongly white balanced image in the next part of the activity.
____________________________________________________________________

Results:
ORIGINAL, REFERENCE WHITE, GRAY WORLD

We can see from the images that the Reference White Algorithm produced a better result. The colors are more vibrant. However, some patches became white, like yellow and light orange. But in comparison with the Gray World algorithm, this is still acceptable. Majority of the colors are now green or white in the last image, showing poor quality. Next, we investigate the effects of these two methods for objects with similar hues.

Again, I use the tungsten white balance setting because it is really not suited for the illumination condition (fluorescent). We implement the two algorithms and the results are shown below.

Results:
ORIGINAL, REFERENCE WHITE, GRAY WORLD

Again, the Reference White Algorithm produced the better result. This is because the red component of white is really red and we use this as the divisor in equations 1-3. In comparison with the Gray World algorithm, the red component of the whole image is not necessarily also red and thus the error is produced. From both the Macbeth chart and the blue objects, we observe that the blue parts become green and those with some red component turns white.
____________________________________________________________________

I give myself 10 points for this activity since both algorithms were implemented and the results were very clear. Thanks to Marge Maallo for her help in this activity.

Activity 14: Stereometry

2008-08-26T16:44:00.000-07:00

Multiple 2D views allow the computation of depth information. In this activity, we use two identical cameras to capture and reconstruct a 3D object. Images are essentially 3D objects projected onto a 2D frame. From object points at (x,y,z), an image is reduced to (x,y) with z projected as a function of x and y and camera object geometry. In applications such as microscopy or terrain imaging, it is sometimes more informative to preserve the depth information z. In this way, the 3D image may be inspected at different view angles. Stereo imaging is the technique we will use in this experiment and it is inspired by how our two eyes allow us to discriminate depth.

From similar triangles, the following equations are generated. If done on several points on the image, we can reconstruct the object's 3D shape using the z's calculated.

source: M. Soriano, A14 - Stereometry.pdf

Shown below are two images of a simple 3D object that I captured using an Olympus Stylus 770SW digital camera. The vertices of the Rubik's cube will serve as reference points. From figure 1, b was measured to be 30mm using a ruler.
____________________________________________________________________

Next, we calibrate our camera using the algorithm in activity 11 to find the focal length f and the distances from the camera center, x1 and x2. The matrix A was calculated to be,
A = [
-12.892329
8.19365
-1.2825517
121.45776
-5.9037883
-6.0874536
12.008296
80.296042
-0.0214988
-0.0197136
-0.0140380
1]

We now get A(1:3,1:3) and factorize using RQ factorization in order to convert the matrix into the upper diagonal matrix K. We also need to ensure that the factorized A(3,3) element is 1 by dividing the whole matrix by K(3,3).RQ factorization code by Michael Overton
(source: http://www.cs.nyu.edu/faculty/overton/software/uncontrol/rq.m)

function [R,Q]= rq(A)

m = size(A,1);
n = size(A,2);
if n < m
error('RQ requires m<=n')
end

p = mtlb_fliplr(mtlb_eye(m));
Atp = A’*P;
[Q2,R2] = qr(Atp);
bigperm = [P zeros(m,n-m; zeros(n-m,m) mtlb_eye(n-m)];
q = (Q2*bigperm)’;
r = (bigperm*R2*P)’;

K =
[-14.756525, 0.0554842, 121.46413;
0, 14.493566, 80.290264;
0, 0, 1]

We can see that the camera center is located at xo = 121.46413 and yo = 80.290264. The actual size of the image is 256x217. Thus, we can say that the camera center is relatively close to the image center (considering only x since we only need x1 and x2, i.e. 121.5*121.5 = 243, which is close to 256). From this value of xo, we can now solve for x1 and x2 for different points on the cube. Then the height z will follow from the formula. We plot our results in 3D. (see Scilab help on splin2d, example 2)
____________________________________________________________________
Reconstructed images of the Rubik's cube:
Here we have 3 views of the reconstructed image. Although the edges are curved, the shape of the cube is clearly observed. A better reconstruction may be possible by locating more points. For the resulting images, only the 7 visible corner points in the cube were used in finding the correspondence between the 1st and 2nd images.

____________________________________________________________________

I give myself 9 points for this activity since the reconstruction was properly done but the result was not very convincing. The result looked like a cube with smoothened corners. Thanks to JC Nadora, Julie Dado, and Jeric Tugaff for their help in the use of splin2d and interp2d.

Activity 13: Photometric Stereo

2008-08-07T11:13:00.000-07:00

In this activity, we will use the concept of photometric stereo to reconstruct an object in 3D from image using shadow cues. The matlab file Photos.mat contains the 4 images that we will work with and the elements of matrix V is also given. These are the locations in space of the light source illuminating the object.

We can estimate the shape of the surface by capturing multiple images of the surface with the sources at different locations. The information about the surface will be coded in the shadings obtained from the images. We can define a matrix V such that,
where each row is a source and each column corresponds to the x,y,z component of the source's location in space. If we take N images of the surface using each of these N sources, then for each point (x,y) on the surface we have,
We can now solve for g since I and V are known and we normalize g by its length in order to get the normal vector.After this, we calculate the object's shape from the normal vectors using the following equations.
source: M. Soriano, A13 - Photometric Stereo.pdf
____________________________________________________________________
//load images
loadmatfile('Photos.mat');

//enter V values
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];

//make matrix I
I1 = I1(:)';
I2 = I2(:)';
I3 = I3(:)';
I4 = I4(:)';
I = [I1;I2;I3;I4];

a = 1e-6; //additive factor
g = inv(V'*V)*V'*I;
mod = sqrt((g(1,:).*g(1,:))+(g(2,:).*g(2,:))+(g(3,:).*g(3,:))); //length of g
mod = mod+a;

for i = 1:3
n(i,:) = g(i,:)./mod; //calculate normal vectors
end

nz = n(3,:)+a;
dfx = -n(1,:)./nz; //partial differentiation
dfy = -n(2,:)./nz;
z1 = matrix(dfx,128,128); //reshape the matrix to 128x128
z2 = matrix(dfy,128,128);
int1 = cumsum(z1,2); //a simpler way to integrate: cumulative sum
int2 = cumsum(z2,1);
z = int1+int2;
scf(0); plot3d(1:128, 1:128, z); //show reconstruction
____________________________________________________________________

RESULTS:

4 sample images

reconstructed image
This 3D image clearly looks like a sphere. The cross mark at the center is also observed. Thus we can conclude that our image was properly reconstructed in 3d using the photometric stereo technique.

____________________________________________________________________

I give myself 10 points for this activity since the sphere was properly reconstructed. Thanks to Jeric Tugaff, Benj Palmares, and Mark Bejemino for their help in this activity.

Activity 12: Correcting Geometric Distortions

2008-08-03T12:20:00.000-07:00

In this activity, we aim to correct geomertic distortions produced by the camera. A sample distorted image is given and is shown below. We need to make the straight lines look straight.

Procedure:
1. Using the given image, we count the number of pixels down and across one box where the grid is most undistorted. Or we can capture another image of a regular grid. Just make sure that the image plane is parallel to the grid plane.
2. We generate the ideal grid vertex coordinates and compute coefficients c1 to c8.

3. Pixel per pixel in the ideal image, use equations 8 and 9 to determine the location of a point in the distorted image.

4. If the computed distorted location is integer-valued, copy the grayscale value from the distorted image onto the blank pixel. If the location is not integer-valued, compute the interpolated grayscale value using equation 10. In pixel locations where there are no transferred pixels, interpolate the graylevel values from the neighboring filled pixels.

source: M. Soriano, A12 - Correcting Geometric Distortions.pdf
____________________________________________________________________

The resulting enhanced image is shown below. There were still some pixels that were not replaced to the proper location but the distortions were corrected. I give myself 9 points for this activity since the enhanced image was not perfectly done. My collaborators for this activity were Mark Leo Bejemino and Jeric Tugaff.

Activity 11: Camera Calibration

2008-07-29T16:27:00.000-07:00

In this activity, we will try to model the physical processes involved in the geometric aspects of image formation. We will try to calibrate the camera by looking at corresponding world points and image points.

An image is a projection of brightness values from surfaces existing in 3D world space to 2D sensor space. One whole dimension is lost in the projection. Nevertheless, under certain constraints it is possible to recover the lost information. The process involves carefully calibrating the camera by taking several views of a scene.

First, we take a picture of a checkerboard pattern, also known as a Tsai grid, using an ordinary digital camera. We pick an origin and 20-25 more points in the image while taking note of its real world coordinates. The squares are 1x1 inches thick. Using the locate function in Scilab, we get the image coordinates of the same 20-25 points we picked. Next, we set up the matrix Q shown by equation 13 for values of i from 1 to 25. We then solve for a using equation 15 and then plug its elements into equations 11 and 12. For the last two equations, we set a_34 to be equal to 1.

RESULTS:

Shown below is the Tsai grid together with the points (red) we wish to use in the calibration. The origin is set at (0,0,0). The left side of the board is x, the right side is y, and the vertical axis is z.

(xo,yo,zo)
1. (0,0,2)
2. (0,0,5)
3. (0,0,9)
4. (0,0,12)
5. (2,0,1)
6. (5,0,2)
7. (3,0,4)
8. (2,0,7)
9. (7,0,6)
10. (4,0,6)

11. (6,0,9)
12. (2,0,10)
13. (8,0,11)
14. (1,0,3)
15. (0,0,0)
16. (0,6,12)
17. (0,8,10)
18. (0,3,10)
19. (0,2,1)
20. (0,4,3)
21. (0,7,6)
22. (0,5,7)
23. (0,6,2)
24. (0,4,6)
25. (0,2,8)

Using locate, their corresponding 2D image coordinates are as follows:
(yi, zi)
(130.53181, 67.545109)
(130.53181, 114.07882)
(131.00665, 178.18139)
(130.53181, 229.46344)
(106.31529, 45.702754)
(65.954416, 51.40076)
(92.545109, 90.811966)
(105.36562, 143.04368)
(33.665717, 115.50332)
(78.300095, 122.151)
(47.910731, 172.95821)
(104.89079, 194.32574)
(15.147198, 209.52042)
(118.18613, 79.890788)
(130.05698, 37.155745)
(211.25356, 228.51377)
(241.16809, 188.62773)
(168.04368, 192.90123)
(152.849, 45.702754)
(178.96486, 70.868946)
(223.12441, 115.50332)
(194.63438, 136.87085)
(206.03039, 47.127255)
(179.91453, 122.151)
(155.22317, 159.18803)

We set up matrix Q given by equation 13 and solving for a using equation 15 will give us the elements of matrix a.
a =
- 13.945461
8.8516947
- 0.6362905
130.1377
- 4.0454015
- 4.1159857
14.796702
37.040714
- 0.0178636
- 0.0176363
- 0.0052799

Using these values, we plug them into equations 11 and 12 and see whether these new yi's and zi's are close to the values previously obtained.

(new_yi, new_zi)
(130.24043, 67.345268)
(130.39869, 114.03467)
(130.61791, 178.7028)
(130.78881, 229.11771)
(105.9554, 45.617241)
(65.69975, 51.556447)
(92.680325, 90.881065)
(105.45811, 142.91477)
(34.035999, 115.62425)
(78.64952, 122.2471)
(48.193941, 172.64721)
(105.19649, 194.0998)
(14.486398, 209.56016)
(118.26943, 80.084525)
(130.1377, 37.040714)
(211.37148, 228.57471)
(241.39165, 188.6588)
(168.09928, 193.06864)
(153.42664, 45.448497)
(179.10783, 71.109685)
(222.85396, 114.82174)
(194.25085, 137.20802)
(205.9424, 47.461688)
(180.14164, 121.80884)
(154.74534, 159.5493)

We can compute for the accuracy of the values we obtained by finding the difference between the original and new coordinates.

diff =
- 0.2913883, - 0.1998413
- 0.1331196, - 0.0441513
- 0.3887400, 0.5214104
0.2569913, - 0.3457294
- 0.3598862, - 0.0855131
- 0.2546660, 0.1556874
0.1352158, 0.0690990
0.0924834, - 0.1289110
0.3702816, 0.1209236
0.3494252, 0.0961056
0.2832094, - 0.3110003
0.3057051, - 0.2259333
- 0.6608001, 0.0397375
0.0832921, 0.1937365
0.0807186, - 0.1150312
0.1179217, 0.0609362
0.2235590, 0.0310715
0.0556003, 0.1674038
0.5776353, - 0.2542572
0.1429685, 0.2407392
- 0.2704486, - 0.6815880
- 0.3835267, 0.3371792
- 0.0879896, 0.3344329
0.2271150, - 0.3421524
- 0.4778311, 0.3612629

We get a mean of 0.26442074 difference in yi and a mean of 0.218553368 difference in zi. Thus, we can say that our calibration is very accurate since there is only a difference of less than a pixel in the values gathered.
____________________________________________________________________

im = gray_imread('tsai grid.jpg');
imshow(im);
x = locate(25,1);
yi = x(1,:);
zi = x(2,:);

xo = [0 0 0 0 2 5 3 2 7 4 6 2 8 1 0 0 0 0 0 0 0 0 0 0 0];
yo = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 8 3 2 4 7 5 6 4 2];
zo = [2 5 9 12 1 2 4 7 6 6 9 10 11 3 0 12 10 10 1 3 6 7 2 6 8];

for i = 1:length(xo)
Q((2*i)-1,:) = [xo(i) yo(i) zo(i) 1 0 0 0 0 -(yi(i)*xo(i)) -(yi(i)*yo(i)) -(yi(i)*zo(i))];
Q(2*i,:) = [0 0 0 0 xo(i) yo(i) zo(i) 1 -(zi(i)*xo(i)) -(zi(i)*yo(i)) -(zi(i)*zo(i))];
d((2*i)-1,:) = yi(i);
d(2*i,:) = zi(i);
end

a = inv(Q'*Q)*Q'*d;

for j = 1:length(xo)
yi_new(j,:) = ((a(1)*xo(j))+(a(2)*yo(j))+(a(3)*zo(j))+a(4))/((a(9)*xo(j))+(a(10)*yo(j))+(a(11)*zo(j))+1);
zi_new(j,:) = ((a(5)*xo(j))+(a(6)*yo(j))+(a(7)*zo(j))+a(8))/((a(9)*xo(j))+(a(10)*yo(j))+(a(11)*zo(j))+1);
end
____________________________________________________________________

To verify if the calibration is correct, we predict the image coordinates of some cornerpoints of the checkerboard which were not used in the calibration (yellow). The number of points should be greater than 11. These points are, in world coordinates:

(3,0,8)
(0,1,5)
(0,4,4)
(0,6,5)
(0,3,9)
(0,7,11)
(4,0,12)
(6,0,1)
(8,0,5)
(2,0,3)
(5,0,5)
(7,0,2)

We predict that these points will be transformed into the following by using equations 11 and 12.
(92.030252, 158.4637)
(142.06384, 111.83289)
(179.44842, 87.812934)
(207.50173, 99.481483)
(167.82, 175.48687)
(226.15495, 208.91806)
(77.116722, 229.33722)
(51.635668, 31.057816)
(18.52981, 94.693389)
(105.79334, 77.32756)
(64.717941, 102.67897)
(36.148858, 44.327313)

//prediction code
xo = [3 0 0 0 0 0 4 6 8 2 5 7];
yo = [0 1 4 6 3 7 0 0 0 0 0 0];
zo = [8 5 4 5 9 11 12 1 5 3 5 2];

a = [-13.945461 8.8516947 -0.6362905 130.1377 -4.0454015 -4.1159857 14.796702 37.040714 -0.0178636 -0.0176363 -0.0052799];

for j = 1:length(xo)
yi(j,:) = ((a(1)*xo(j))+(a(2)*yo(j))+(a(3)*zo(j))+a(4))/((a(9)*xo(j))+(a(10)*yo(j))+(a(11)*zo(j))+1);
zi(j,:) = ((a(5)*xo(j))+(a(6)*yo(j))+(a(7)*zo(j))+a(8))/((a(9)*xo(j))+(a(10)*yo(j))+(a(11)*zo(j))+1);
end

After running the original code, these points were converted into:
(91.919444, 158.17825)
(141.77925, 111.69473)
(179.29545, 87.811292)
(207.46461, 99.490062)
(167.2938, 175.12682)
(225.77969, 208.89718)
(76.901337, 228.70553)
(51.993604, 31.260393)
(18.750844, 94.88024)
(105.83348, 77.396374)
(64.879166, 102.73954)
(36.482614, 44.520168)

Again, we get the difference between the two sets of values and find the mean difference for yi and zi.
mean difference for yi = 0.234702
mean difference for zi = 0.179843167
____________________________________________________________________

I give myself 10 points for this activity since the results of the camera calibration was good. The values differed by a value that is less than 1 pixel. Good predictions resulted because of this. My collaborators for this activity were Jeric Tugaff and Julie Dado.

Activity 10: Preprocessing Handritten Text

2008-07-22T11:16:00.000-07:00

In this activity, we will try to extract handwritten text from an imaged document with lines. In handwriting recognition, individual examples of letters must be extracted. This is another real-world problem like Activity 9 wherein we will need to use all the image processing techniques that we have learned so far.

We are given two images to choose from. Both have handwritten text within lines. I will be using the 2nd image. It is shown below together with a cropped portion that we will be processing.

<-- middle left portion First, we have to remove the horizontal lines and isolate the text. We can do this by looking at its FT and masking the vertical lines in Fourier space since these correspond to the horizontal lines in image space. After inversion, the text should be isolated from the lines.

FT image
From this FT image, we need to mask the prominent frequencies in the vertical direction in order to isolate the text from the horizontal lines.

MASKED FT image and ENHANCED image

im = gray_imread('im.jpg'); //size of image is 89 rows x 125 columns
IM = fftshift(fft2(im));
IM(1:44,62:63) = 0;
IM(46:89,63:64) = 0;
F = real((IM).*conj(IM));
scf(0); imshow(log(F+1),[]); xset('colormap',jetcolormap(256));
newim = abs(fft2(IM));
scf(1); imshow(newim,[]);

The line removal is not perfect but the words are still comprehensible. Next, we binarize the image and clean it of stray pixels, holes, and the like. However, the text is in black, so we have to invert first and make the text white and the background black before binarizing.

BINARIZED image and CLEANED image

The thickness of most of the letters are already at 1 pixel after binarization. Thus, we only need to close the gaps made by the removal of the horizontal lines. We use the closing operation using a 4x1 rectangle structuring element and the result is shown in the right image above. Unfortunately, this is the best we can do without making the words indistinguishable. The words, "remote", "cable", and a hint of "control" can still be recognized. By using the bwlabel function, we assign a label to each letter (or in this case, to each cluster of white pixels).

No. of clusters = 84
No. of actual letters in text = 46
Error = 82.6%

im = gray_imread('enhanced.jpg');
im = 1 - im;
im = im2bw(im,0.7);
scf(0); imshow(im);

SE = ones(4,1);
im = dilate(im,SE);
im = erode(im,SE);
scf(1); imshow(im);

[L,n] = bwlabel(im);
____________________________________________________________________

I give myself 10 points for this activity even if the result of the cleaning part is not so nice. The horizontal lines were removed and the words are still comprehensible after doing the morphological operations. My collaborator for this activity is Julie Dado.

Activity 9: Binary operations

2008-07-17T20:50:00.000-07:00

In this activity, we aim to use what we have learned so far in finding the area of a piece of circular punched paper. This may seem easy but these many circles are scattered: some are touching or even on top of each other. We will therefore analyze several regions of interest in one image.

PROCEDURE:
1. If the image is large, divide it into several 256x256 sub-images so that less memory space is used. Another way is to just resize the image. Open the image in grayscale using im2gray or gray_imread.
2. Plot the histogram of the sub-image and determine its thresholding value. Binarize the image using im2bw.
3. Do morphological operations to remove isolated pixels, separate connected blobs or fill out holes. The opening and closing operations can do this "image cleaning." The opening operation is just an erosion followed by a dilation. On the other hand, the closing operation is a dilation followed by an erosion.
4. Use bwlabel to label each blob on the image.
5. Find the areas of each disk of punched paper. Repeat for the other sub-images and give the best estimate of the area.
____________________________________________________________________

The image used in this activity and one of its sub-images are shown below.

After binarization and morphology (closing first then opening operations),

HISTOGRAM OF AREAS
From the histogram, it seems as if the correct value of the area lies between 500 and 600 pixels since the histogram bars are highest around these values. The values smaller than 500 are those disks that got cropped off when the sub-images were created or those that do not look like a circle anymore after the morphological operations. Values larger than 600 are those overlapping disks or those too close together that they got combined by the morphological operations. After removing these values, we get the mean and standard deviation. These tell us the best estimate for the area of an individual punched paper.

Mean area = 538 pixels
Standard deviation = 24.83

We need to verify if this value for the area is acceptable. We can do this by finding non-overlapping disks and then calculating the areas.

Mean area = 549 pixels
Standard deviation = 13.05

error for area = 2%

Also, the standard deviations are quite close. Smaller standard deviations indicate that the other areas are close to the mean value.

____________________________________________________________________

stacksize(4e7);

subimg = 11; //no. of sub images
c = 1;
for i = 1:subimg
im = gray_imread('c2_0'+string(i)+'.jpg');
im = im2bw(im, 0.85); //convert to BW with corresponding threshold value
//scf(0); imshow(im);

image = ones(5,5);
SE = mkfftfilter(image,'binary',2); //structuring element
im = dilate(im,SE); //closing operation
im = erode(im,SE);
im = erode(im,SE); //opening operation
im = dilate(im,SE);
//scf(i); imshow(im);

[L,n] = bwlabel(im); //label the individual blobs

//area calculation of individual blobs
for j = 1:n
roi = (L==j);
A(c) = sum(roi);
c = c + 1;
end
end

scf(12); histplot(length(A),A);
title('Histogram of disk areas');
xlabel('bins'); ylabel('frequency');

//remove too small and too large areas
range_area = find(A>500 & A<600);
mean_area = sum(A(range_area))/length(range_area); //mean area
std_area = stdev(A(range_area)); //standard deviation
____________________________________________________________________

I give myself 10 points for this activity since a good estimate of the area of the punched paper was given. Thanks to Jeric Tugaff for help in the histogram and in removing the unwanted small and big area values.

Activity 8: Morphological Operations

2008-07-15T16:05:00.000-07:00

Morphology refers to shape or structure. In image processing, classical morphological operations are treatments done on binary images, particularly aggregates of 1's that form a particular shape, to improve the image for further processing or to extract information from it. All morphological operations affect the shape of the image in some way, for example, the shapes may be expanded, thinned, internal holes could be closed, disconnected blobs can be joined. In a binary image, all pixels which are “OFF” or have value equal to zero are considered background and all pixels which are “ON” or 1 are foreground.

Source: M. Soriano, A8 - Morphological Operations.pdf

In this activity, we will investigate two kinds of morphological operations: dilation and erosion. These have functions in Scilab and we will apply them to simple geometric shapes for now.

(1) Dilation

In set theory,
(2) Erosion

The erosion operator is defined by

We apply this to 5 shapes: 50x50 pixel square, triangle (50 pixel base, 30 pixel height), circle with radius 25, hollow square (60x60 pixels with 4 pixel wide edges), and plus sign (8 pixels wide and 50 pixels long) using 4 different structuring elements, namely, 4x4 square, 2x4 rectangle, 4x2 rectangle, and a 5 pixel long, 1 pixel thick cross. We have already predicted the outcome of this operation on these images and we will now verify them using Scilab. The predictions are submitted separately on a sheet of yellow pad paper.

ORIGINAL SHAPES

i) SE = 4x4 square
DILATE

ERODE

ii) SE = 2x4 rectangle
DILATE

ERODE

iii) SE = 4x2 rectangle
DILATE

ERODE

iv) SE = 5 pixels long, 1 pixel wide cross
DILATE

ERODE

(3) THIN

Function thin performs thinning of binary objects. The resulting image, the skeleton, is not always connected and is very sensible to noise. It is also slower than a good skeletonization algorithm (see skel). For thin shapes, it should work faster and provide better quality.

(4) SKEL

Function skel performs skeletonization (thinning) of a binary object. The resulting medial axis is multi-scale, meaning that it can be progressively pruned to eliminate detail. This pruning is done by thresholding the output skeleton image. The algorithm computes skeletons that are guaranteed to be connected over all scales of simplification. The skeletons are computed using the euclidean metric. This has the advantage to produce high-quality, isotropic and well-centered skeletons in the shape. However the exact algorithm is computationally intensive.

____________________________________________________________________

I give myself 9 points for this activity because some of the predictions were wrong especially in the erosion part. Thanks to Julie Dado, Angel Lim, Lei Uy, JC Nadora, Jeric Tugaff, and Mark Leo Bejemino for their help in the predictions in this activity.

Activity 7: Enhancement in the Frequency Domain

2008-07-12T18:07:00.000-07:00

This activity aims to investigate on 3 things, namely: (1) Anamorphic property of the FT, (2) Fingerprints: Ridge enhancement, and (3) Lunar landing scanned pictures: Line removal.

(1) Anamorphic property of the FT

The Fourier Transform (FT) of a signal gives the spatial frequency distribution of that signal. Essentially, the Fourier Theorem states that any image or signal can be expressed as a superposition of sinusoids. Unique to the 2D FT (as compared to the 1D FT) is the fact that rotation of the sinusoids result to the rotation of their FT's.

Here, we look into what the FT of a sinusoid looks like and see what happens when its frequency is changed, when it is rotated, and when it is superimposed with another sinusoid. We create a sinusoid using Scilab, take its FT, and expect that the rotation of the sinusoid results into the rotation of its FT.

nx = 100; ny = 100;
x = linspace(-1,1,nx);
y = linspace(-1,1,ny);
[X,Y] = ndgrid(x,y);
f = 4; //frequency
theta = 30; //angle of rotation
z = sin(2*%pi*f*X); //sinusoid
z1 = sin(2*%pi*f*(Y*sin(theta) + X*cos(theta)));
scf(1);
imshow(z,[]);
Fz = fft2(z); //or change to z1 for rotated sinusoid
scf(2); imshow(abs(fftshift(Fz)),[]);

SINUSOID (f=4) and FT
Here we observe the FT of the sine wave. The sinusoid is horizontal while the two dots are oriented vertically. This is because doing a Fourier transform maps an object into the Fourier space (or k-space or inverse space). We can say that the frequency axis in the FT image is the vertical one and the frequency value corresponds to how far apart a dot is from the origin. We only need to consider one dot since the other is just a mirror image. Thus, we expect that increasing the frequency of the sinusoid will result in making the dots in the FT image farther apart.

SINUSOID (f=20) and FT
Now that the frequency is increased, the dots are farther away as expected. Next, we have a rotated sinusoid which is oriented 30 degrees to the left with respect to the vertical axis. Since the FT maps this wave into inverse space, we expect that the resulting image will have dots also rotated 30 degrees but now with respect to the horizontal axis. Indeed, this is what happens. Finally, we use a combination of 2 sine waves and look at its resulting FT image. The form is of z = sin(2*%pi*4*X).*sin(2*%pi*4*Y);

ROTATED SINUSOID (f =4) and FT; COMBINATION and FT

The product of the 2 sinusoids generates the 3rd figure above. A checkerboard pattern results and the FT image consists of 4 dots forming the corners of a square. We can think of this as the product of 2 horizontal dots and 2 vertical dots. The axes being inverted are not evident here since the pattern is square. One can also show that when we add the two sinusoids, its FT looks like a cross rather than a square. We can conclude that the FT of the combination of two signals, whether you add, subtract, multiply, or divide them, will result in the combination of the FT's of the individual signals.
____________________________________________________________________

(2) Fingerprints: Ridge Enhancement

We obtained a scanned image of a fingerprint for this part of the activity. If it is not yet in grayscale, it should be converted. After this, we must subtract the mean from the image in order to remove the dark current (DC). Doing this will make all the other spatial frequencies present in the image visible. (source of fingerprint image: A7 - Enhancement in the Frequency Domain.pdf)

ORIGINAL & MEAN-CENTERED IMAGES

We generate the FT image for the fingerprint and is shown at the left.

We want to enchance the appearance of the ridges while removing the blotches. We can do this by applying a filter using the mkfftfilter function in Scilab. We will use this instead of mkfilter because it is easier to work in the Fourier domain than in the image domain. We then apply the filter to the image in Fourier space in order to enhance certain features. To select the appropriate filter, we need to take a look at the image's FT first.

Our objective is to enhance the ridges or make the image sharper. There are two types of image enhancement in the Fourier domain in general, low-pass filtering and high-pass filtering. Applying a low-pass filter corresponds to reducing the high frequency content, which in turn results into blurring or smoothing. High-pass filtering, on the other hand, increases the magnitude of the high-frequency components relative to the low-frequency components. This results in sharpening. Thus, we need to apply a high-pass filter to the FT image in order to enhance the fingerprint ridges. Such a filter is shown at the left.

What this does is that only those frequencies in the FT image that fall in the white parts of the filter is allowed to pass and all others are filtered out. Note that the frequencies near the center are low frequencies. This way, we block out the low frequencies and the image obtained after doing an inverse FT is sharpened (and of course, also inverted). The code for generating these images is as follows.

//Original code c/o Scilab Help - mkfftfilter
//Increase the stack size because images are memory consuming
stacksize(4e7);

image=gray_imread('fingerprint.png');
image2=image-mean(image); //removes DC
scf(1); imshow(image2);

IM=fft2(image2);

F=real((IM).*conj(IM));
scf(2); imshow(fftshift(log(F+1)),[]); //the use of log(F+1) allows to observe great amplitude variations
xset('colormap',jetcolormap(256));

//Transfer function
h=mkfftfilter(image2,'exp',20); //exponential with "radius" = 20
scf(3); imshow(h);

//High-pass filter
IM2=IM.*fftshift(1-h);
im=real(fft2(IM2));
scf(4); imshow(im);

ORIGINAL & ENHANCED IMAGES

Clearly, the fingerprint ridges are enhanced.

____________________________________________________________________

(3) Lunar Landing scanned pictures: Line Removal

Here we aim to remove the vertical lines from the given picture by filtering in the Fourier domain. We should make sure that the image is in grayscale first. Its FT is also shown below.

ORIGINAL & FT IMAGES

From part 1 of this activity, we saw that the spatial frequencies of an image are mapped into the Fourier plane at perpendicular directions. Hence, vertical lines in the image will correspond to dots in the horizontal axis of the FT image. To remove the vertical lines from the satellite moon image, we need to mask those prominent dots at the horizontal axis while still retaining the center spot so that no information is lost. The mask that was used is shown here.

ORIGINAL & SMOOTHENED IMAGE

The vertical lines are now removed and the resulting image still retains its information, meaning it still looks like a photograph of the moon.
____________________________________________________________________

I give myself 10 points for this activity since the objectives were met and the images were enhanced. Special thanks to Julie Dado and Jeric Tugaff for their help in this activity.

Activity 6: Fourier Transform Model of Image Formation

2008-07-08T11:28:00.000-07:00

This activity is divided into 4 parts namely: (1) Familiarization with the discrete Fast Fourier Transform (FFT), (2) Simulation of an imaging device, (3) Template matching using correlation, and (4) Edge detection using the convolution integral.

(1) Familiarization with the discrete FFT

This part is straightforward in the sense that we try to understand the properties of the FFT and see what happens to an image after doing the operation. As a start, we generate a 128x128 image of a circle in MSPaint and apply the FFT.

I = imread('circle.bmp');
Igray = im2gray(I);
FIgray = fft2(Igray);
imshow(abs(FIgray),[]); //intensity image
//imshow(abs(fftshift(FIgray)),[]); //fft shifted image
//imshow(abs(fft2(FIgray))); //fft twice
xset("colormap",hotcolormap(32)); //add color to your images

CIRCLE: ORIGINAL image, FFT INTENSITY image, FFT SHIFTED image, DOING FFT TWICE

For the 1st FFT image (2nd from left), we observe spots at the corners of the image. These are actually the 4 quarters of a circle and we predict that once this image is FFT shifted, we will observe the correct analytical Fourier Transform of a circle, which is an Airy disk: a bright spot centered about concentric rings outside of it. Indeed, this is what we observe in our image (3rd from left). Upon applying the FFT twice on our original image, what we observe is that we return to our original image. However, this is actually an inverted image of the original and it is just not obvious. We can verify this by using an "A" as our image.

"A": ORIGINAL image, FFT INTENSITY image, FFT SHIFTED image, DOING FFT TWICE

The "A" is flipped with respect to the horizontal and is now inverted with respect to the original image. Thus, doing FFT twice will result in an inverted image.

____________________________________________________________________

(2) Simulation of an imaging device

Here we try to simulate what happens in a camera. We vary the camera lens radii and see what the FFT images look like.

c = gray_imread('circle.png'); //medium aperture
//c = gray_imread('bigcircle.png'); //big aperture
//c = gray_imread('smallcircle.png'); //small aperture
vip = gray_imread('vip.png');
Fc = fftshift(c); //aperture is already in the Fourier Plane and need not be FFT'ed
Fvip = fft2(vip);
Fcvip = Fc.*(Fvip);
Icvip = fft2(Fcvip); //inverse FFT
FImage = abs(Icvip);
imshow(FImage, []);

ORIGINAL, SMALL APERTURE, MEDIUM APERTURE, BIG APERTURE

A finite lens radius in a camera means that the lens can only gather a limited number of rays reflected off an object, therefore reconstruction of the object is never perfect. A smaller lens aperture will gather fewer rays reflected off an object and this results in a blurred image. Comparing this with the bigger aperture, the image is sharper and most closely resembles the original.
____________________________________________________________________

(3) Template matching using correlation

Template Matching is a pattern recognition technique suitable for finding exactly identical patterns in a scene such as in the case of finding a certain word in a document. Here, we apply this technique in finding where the letter A's are in a sample text.

rain = gray_imread('rain.png'); //image of the text
a = gray_imread('a2.png'); //image of the letter to correlate with text, same font type and size
Frain = fft2(rain);
Fa = fft2(a);
p = Fa.*conj(Frain);
im = ifft(p); //inverse FFT
imshow(abs(fftshift(im)),[]);

ORIGINAL & RESULTING IMAGE

In the resulting image, we see that bright dots are located where there are letter A's in the text. These words are "rain", "Spain", "stays", "mainly", and "plain." There are 5 bright dots that correspond with 5 A's in the text. We have demonstrated here how the shape of the letter A was matched with the text using correlation.

(4) Edge detection using the convolution integral

Edge detection can be seen as template matching of an edge pattern with an image. We can detect the edges of the previous "VIP" image, for example, by applying a pattern oriented in a particular detection. The code and resulting images are shown below.

hor = [-1 -1 -1; 2 2 2; -1 -1 -1]; //horizontal pattern
vert = [-1 2 -1; -1 2 -1; -1 2 -1]; //vertical pattern
spot = [-1 -1 -1; -1 8 -1; -1 -1 -1]; //spot pattern
pattern = hor;
vip = gray_imread('vip.png');
im = imcorrcoef(vip, pattern);
imshow(im);

ORIGINAL, HORIZONTAL EDGE, VERTICAL EDGE, SPOTS

As we can see, the edges are more prominent depending on what kind of pattern the image was convolved with. If a horizontal pattern was used, only the horizontal portions of the "VIP" image is shown (2nd to the left). The same is true for the vertical pattern. The spot pattern only shows portions containing 1 white pixel and the other pixels are merged with the background. That is why there are some holes in the letter V and P in particular upon close inspection.
____________________________________________________________________

I give myself 10 points for this activity since the all of the 4 objectives were met. I found this activity to be interesting and the techniques discussed here useful either for research and even for leisure. Thanks to Jeric Tugaff for help in the code for parts 3 and 4.

Activity 5: Physical Measurements from Discrete Fourier Transforms

2008-07-03T11:20:00.000-07:00

In this activity, we will be applying the Discrete Fourier Transform (DFT) in image processing. First, we will try to understand what the DFT is, how it works, what it does to signals/images, and what information we can get out of it. Scilab already has a built-in fft function and we will use this so that we no longer have to implement the math involved and make another code for it.

As an example, we generate a 1-D sinusoid of the form sin (2*pi*f*t), compute its Fourier Transform (FT) and display it.

//Original code c/o Dr. Soriano
//Generate 1-D sinusoid specifying dt, N, and f
T = 2;
N = 256;
dt = T/N;
t = [0:dt:(N-1)*dt];
f = 5;
y = sin(2*%pi*f*t);
subplot(1,2,1);
plot(t,y);
title('1-D sinusoid: f=5');
xlabel('time (seconds)');
ylabel('amplitude');

//Get FT of signal and compute the frequency scale
FY = fft(y);
Fmax = 1/(2*dt);
df = 2*Fmax/N;
f = [-Fmax:df:df*(N/2-1)];

//Display the fftshifted output with frequency axis
subplot(1,2,2);
plot(f,fftshift(abs(FY)));
title('Fourier transform of signal');
xlabel('frequency (Hz)');
ylabel('amplitude');

____________________________________________________________________

The Fourier Transform is an important image processing tool which is used to decompose an image into its sine and cosine components. The output of the transformation represents the image in the Fourier or frequency domain, while the input image is in the spatial domain equivalent. In the Fourier domain image, each point represents a particular frequency contained in the spatial domain image.

The Fourier Transform is used in a wide range of applications, such as image analysis, image filtering, image reconstruction and image compression. (Source: http://homepages.inf.ed.ac.uk/rbf/HIPR2/fourier.htm)

So for images, we work in the spatial and not the temporal domain. We first find the FT 1-D at a time, i.e. rows first then columns or vice versa. From there we can manipulate the FT depending on the particular application.

Answers to questions:

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

From the Nyquist Theorem, Fmax = 1/2Δt. One can show that the value of Δt is equal to 0.004167 seconds.

2. What is the effect of increasing the number of samples N in the FT?

Working equations: Δt = T/N, Fmax = 1/2Δt, Δf = 2Fmax/N

Increasing N will mean that for the same Δt, T has to increase as well. When this happens, you have the same Fmax but Δf, or the discrete frequency steps in the Fourier domain, will decrease since it has an inverse relation with N as seen from the 3rd equation. The frequency measurement is more precise due this decrease in Δf, thus the lines become sharper (better resolution). We verify this by graphing.

3. What is the effect of decreasing the sampling interval Δt in the FT?

Decreasing Δt will mean that for the same N, T has to decrease as well. Fmax will increase, N stays the same, so Δf will also increase. Δf increasing will result in not so precise frequency measurement and frequency lines that are not sharp (poor resolution). We verify this by graphing.

4. What is the effect of fixing the total time interval T but increasing the number of samples N?

Fixing T but increasing N will mean that Δt decreases accordingly. When this happens, Fmax increases since it has an inverse relation with Δt. Also, Fmax has a coefficient of 2, thus it increases faster than N and from the 3rd equation, Δf also increases. We verify this by graphing. The distinction between sharp frequency lines and those that are not sharp are not so evident here.

I give myself 10 points for this activity since the answers to the questions were verified to be correct by observing the graphs produced in Scilab. Acknowledgment goes to Jeric Tugaff and Julie Dado for help in the Fourier Transform concepts.

Activity 4: Enhancement by Histogram Manipulation

2008-07-01T11:20:00.000-07:00

In this activity, we are asked to enhance a poor contrast image by histogram manipulation and backprojection. First, a poor contrast image was obtained from the internet and its histogram or probability distribution function (PDF) was calculated using Scilab. These are shown below. An image is considered to be of poor contrast if its histogram is not well spread out across the gray level range.

source:http://homepages.inf.ed.ac.uk/rbf/HIPR2crimmins.htm

Next, we solve for the image's cumulative distribution function (CDF) by using the cumsum function in Scilab. For example, if the area of the image is 256x256 and H is the histogram function, then cdf = cumsum([H/(256*256)]); and we plot it. The CDF is just the sum of the histogram values as the gray level increases. In essence, the PDF and CDF give a complete description of the probability distribution of a random variable. The original image's CDF here is crooked but the minimum and maximum values are at the same gray levels as the PDF.

We now proceed in backprojecting the image pixel values, pixel per pixel, by finding its corresponding y-value in the CDF and assigning these values to the enhanced image.

im = imread('gray.jpg');

//Plot histogram
val = [];
num = [];
c = 1;
for i = 0:1:255
[x,y] = find(im == i);
val(c) = i;
num(c) = length(x);
c = c + 1;
end
//plot(val,num/(length(im)));
//title('Normalized histogram of image');
//xlabel('grayscale values');

//Plot CDF
cdf = cumsum([num/(length(im))]);
//plot(cdf);
//title('Normalized CDF');
//xlabel('grayscale values');

//Backprojection
newim = [];
for ii = 1:1:size(im,1)
for jj = 1:1:size(im,2)
newim(ii,jj) = cdf(im(ii,jj));
end
end
imwrite(newim, 'newim.jpg');

Here is the resulting enhanced image.
ORIGINAL & ENHANCED

In theory, the histogram of the enhanced image should be more spread out across the gray levels, indicating good contrast. Due to this, the resulting CDF of the enhanced image should resemble a line with increasing slope. These graphs are shown below.

As expected, the histogram is well spread and the CDF looks fairly like a line. However, the image is observed to be highly pixelated and this is due to low resolution of the images (111 x 111). More computational power will be required for images with higher resolution.

____________________________________________________________________

Another method of image enhancement using histogram manipulation is by mapping the original image's gray levels into a new set of gray levels such that a desired CDF is obtained. We can do this by (1) finding the CDF value from pixel grayscale, (2) trace this value in desired CDF, (3) replace pixel value by grayscale value having this CDF value in desired CDF.
The desired CDF can be any function so long as the y-values are from 0 to 1. I will use the function G(z) = z^2 where z ranges from 0 to 255. The code for enhancing the image is also shown below.

//function of desired CDF
z = [0:255];
G = z^2;
G = G./max(G);

//method implementation
newim = [];
for ii = 1:1:size(im,1)
for jj = 1:1:size(im,2)
E(ii,jj) = cdf(im(ii,jj));
end
end

//interpolation
newim = interp1(G,z,E);
imshow(newim, [0 255]);
newim = round(newim);

//plot histogram
value = [];
number = [];
c1 = 1;
for i = 0:1:255
[x,y] = find(newim == i);
value(c1) = i;
number(c1) = length(x);
c1 = c1 + 1;
end
//plot(value, number/(length(newim)));
//title('Histogram of enhanced image');
//xlabel('grayscale values');

//plot CDF
cdf = cumsum([number/(length(newim))]);
plot(cdf);
title('Enhanced CDF');
xlabel('grayscale values');

And here is the result:
ORIGINAL, 1ST METHOD ENHANCEMENT, 2ND METHOD ENHANCEMENT

Shown below are the histogram and CDF of the enhanced image in the 2nd method. The PDF is also well spread across the graylevel range while the CDF followed the desired CDF, which is a parabola (G(z) = z^2).

The enhanced image from the 1st method is darker than in the 2nd method. In terms of contrast, the 1st method is better. Any of the 2 methods can be used for any other poor contrast image in attempting to enhance it. I cannot decide as of the moment which of the 2 is "better" per se.

I give myself 10 points for this activity since the images were enhanced significantly compared to the original image. This can be shown not only in the images themselves but also from the PDF and CDF of the enhanced images. Thanks to Jeric Tugaff for help in implementing the 2nd method in this activity.

Activity 3: Image Types and Basic Image Enhancement

2008-06-24T11:18:00.000-07:00

In this activity, we are asked to:

1. Collect images of different types and show their properties using the imfinfo function of Scilab.

2. Open the previously acquired scanned image (I will be using a leaf) and convert it into grayscale if it is still colored.

3. Examine the histogram of the images. Convert them into black and white using im2bw. Choose the best threshold from the histogram.

4. Finally, apply the area calculation from last week to determine the area of the image.

____________________________________________________________________

1. Binary Image
FileName: smith.gif

FileSize: 1104
Format: GIF
Width: 244
Height: 153
Depth: 8
StorageType: indexed
NumberOfColors: 2
ResolutionUnit: centimeter
XResolution: 72.000000
YResolution: 72.000000

source: http://www.sharewareplaza.com/Signature-Creator-download_31667.html

2. Grayscale Image
FileName: xray.jpg
FileSize: 4424
Format: JPEG
Width: 183
Height: 193
Depth: 8
StorageType: indexed
NumberOfColors: 256
ResolutionUnit: inch
XResolution: 72.000000
YResolution: 72.000000
source: http://www.answers.com/topic/x-ray?cat=health

3. Truecolor Image

FileName: butterfly.jpg
FileSize: 14300
Format: JPEG
Width: 267
Height: 200
Depth: 8
StorageType: truecolor
NumberOfColors: 0
ResolutionUnit: centimeter
XResolution: 72.000000
YResolution: 72.000000
source: http://monarch-butterfly.info

______________________________________________________________

Here is the scanned image of a leaf I mentioned in last week's report and its properties.

FileName: leaf2.jpg
FileSize: 2851
Format: JPEG
Width: 170
Height: 181
Depth: 8
StorageType: indexed
NumberOfColors: 256
ResolutionUnit: inch
XResolution: 72.000000
YResolution: 72.000000

It was converted into grayscale by using the following lines of code.

image = imread('leaf.jpg');
grayimage = im2gray(image);
imwrite(grayimage, 'leaf2.jpg');

Next, we create a histogram plot of the number of pixels falling under a specific grayscale value. A value of 0 corresponds to black while 255 is for white. The code for plotting the histogram and the resulting graph are shown below.

//original code c/o Jeric Tugaff
image = imread('leaf2.jpg');
val=[];
num=[];
counter=1;
for i=0:1:255
[x,y]=find(image==i); //finds where image==i
val(counter)=i;
num(counter)=length(x); //find how many pixels of leaf have a value of i
counter=counter+1;
end
plot(val, num); //plot

From this histogram, we can see that there is a dip in the 100-180 region. We can now convert our grayscaled leaf image into black and white by using im2bw. Our threshold value can range from 0.4 to 0.7. I used 0.5 in this case. However, note that the parameters for the follow function (the function that follows the contour of the image) require that the object has a value of 1(white) and the background has a value of 0(black). This means that the colors of our current B&W image have to be inverted before we can calculate for its area. We can do this in Scilab by simple matrix manipulation: BWleaf = 1 - BWleaf; and again we write the resulting image into a file.

BEFORE, AFTER, & AFTER STILL...

We can do the same for any other image. The procedure would be as follows:

open the image
convert to grayscale
plot histogram
convert to B&W using a suitable threshold value
invert colors

Whew! After all of that, we are now ready to calculate for the area of the leaf using Green's theorem.

[x,y] = follow(BWleaf);
x_1(1) = x(length(x));
x_1(2:length(x)) = x(1:(length(x) - 1));
y_1(1) = y(length(y));
y_1(2:length(y)) = y(1:(length(y) - 1));
A = abs(1/2*(sum((x.*y_1) - (y.*x_1))));
theo = sum(BWleaf);

Upon comparing the values, we get theo = 7202 and A = 7000.5. This yields an error of only 2.8%. This error comes from the image of the leaf having some concavities. According to the other proof of Green's theorem (from A2 - Area estimation of images with defined edges.pdf),

so long as a closed curve is convex (no concavities) its area may be computed by

summing the area of its “pie slices.”

Thus, Green's theorem for estimating areas will not work all of the time. This may be why the error between calculated and theoretical area is large compared to the regular polygons' areas from last time.

____________________________________________________________________

I give myself 10 points for this activity since like last week, the error between the theoretical and calculated areas are less than 5%. Special thanks goes to Jeric Tugaff and Jorge Presto for the help in plotting the histogram of the grayscaled image.

Activity 2: Area Estimation of Images with Defined Edges

2008-06-17T11:20:00.000-07:00

In this activity, we are asked to estimate areas of different images using Scilab. First, I created several simple images using MS Paint. My first image was a rectangle. It is important that the background is black while the shape is white. The image can be saved under any extension (.jpg, .bmp, etc). Here is the Scilab code that generates the area of this white rectangle.

rect = imread('rect.jpg'); //load image
rect = im2bw(rect,0.5); //convert to binary
[x,y] = follow(rect); //function that follows the contour
x_1(1) = x(length(x)); //create x_i+1 array
x_1(2:length(x)) = x(1:(length(x) - 1));
y_1(1) = y(length(y)); //create y_i+1 array
y_1(2:length(y)) = y(1:(length(y) - 1));
theo = (max(y)-min(y))*(max(x)-min(y)); //Theoretical Area of rectangle
A = abs(1/2*(sum((x.*y_1) - (y.*x_1)))); //Green's function for area estimation

After running the code, the theoretical area for the rectangle is theo = 13152, while the area calculated from Green's theorem is A = 13056. There is an error of 0.73%.

We can also try this method out using various shapes. For a circle, the code is similar. Just replace the image, image name, and the theoretical formula of its area.

theo = 4536
A = 4532
error = 0.09%

Again, for another circle with a smaller radius:

theo = 346
A = 341
error = 1.45%

It can be observed that there is a higher error for a smaller circle than a bigger one. This may be because a smaller circle generated in MS Paint looks less like a "real circle" due to the pixels being square. A bigger circle will look more like a circle since more pixels are being used to smoothen the circle's boundary.

In summary, the Green's method for calculating areas of images with definite edges is nice and efficient since the computed areas are within 5% error of their theoretical values.

I give myself 10 points for this activity since the computed areas were within 5% of their theoretical values. Thanks to Jeric Tugaff, Julie Dado, and Rica Mercado for their help in the Scilab code syntax and in loading the SIP toolbox.

Next time, I will try to see if this method works on a scanned image of a relatively flat object (i.e. leaf, ID, cards, etc...).

Activity 1: Digital Scanning

2008-06-12T10:40:00.000-07:00

We obtained a scanned graph from the Canadian Journal of Physics (1951) and tried to reproduce it in MS Excel using some basic image processing. Here are the steps I used in doing this activity:

1. Get pixel location of the graph's origin, points in the graph and pixel distances of axes tick marks.
2. Subtract the offset due to orientation of the graph's origin and MS Paint's origin.
3. Use ratio and proportion to find corresponding pixels and points on the graph.
4. Graph the data using MS Excel and overlay the original graph for comparison.

And here is the result:The red points are the reconstructed data. The original points on the graph are in black.

I give myself 10 points for this activity since there is a good correspondence between the original and reconstructed data. (The black points were overlapped by the red points, which show good correspondence.)