Activity 19 - Restoration of Blurred image

This activity will show how to restore blurred and noisy images from an a priori knowledge of degradation function and noise. The main idea of image degradation and restoration is presented using the model below:

General model in spatial domain:

where g(x,y) is the degraded image, f(x,y) is the original image, h(x,y) is the degradation function and n(x,y) is the noise model.

General model in frequency domain:

The variables presented are the Fourier transforms of the above equation.

From the above equations, we are to find the degrading function and the noise model. The noise model is described in the previous activity, so we are left with the degrading function, or the motion blurring effect. It is given by:

The image restoration procedure is summarized in the following way:

One filter that is shown in this activity is Weiner filter which is also known as the minimum mean-square error (MMSE) filter by modeling the error through statistical characterization of the noise. It minimizes the average error giving the Weiner filter:

or

where , , and K is a specified constant.

Results

The result of filtering for different types of images are presented. The first image is the original, the second is the blurred and noisy image, the last two are the two types of Weiner filter the last being specified by K = 1.

As seen, the Weiner filter restores the blurred and noisy image but it can't recover the original image itself. The phrase can still be discerned from the restoration but not fully. The last image is not a good restoration using the constant K. Even if K value is varied, no visible change in the restored image is observed.

I give myself 10 pts for understanding and finishing the activity. :)

References
M. Soriano. A19 – Restoration of blurred image. pdf file
Digital Image Processing. Gonzales and Woods.

Scilab Code

// Degradation

chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act19\');
img = gray_imread('k.bmp');
//img = gray_imread('phrase_ei.bmp');
//img = gray_imread('people1.bmp');
FTimg = fft2(img);
FT = fftshift(abs(FTimg));

a = .05;
b = .05;
T = 1;

// noise
y = grand(size(img,1), size(img,2), 'nor', 0.1, 0.1); // gauss

FTy = fft2(y);

for i = 1:size(FT,1)
for j = 1:size(FT,2)
H(i,j) = T/(%pi*(i*a + j*b)) * sin(%pi*(i*a + j*b)) * exp(-%i*(%pi*(i*a + j*b)));
end
end

G = H.*FTimg + FTy;
Sn = FTy.*conj(FTy);
Sf = FTimg.*conj(FTimg);

a = (H.^(-1));
b = ((H.*conj(H))./((H.*conj(H)) + (Sn./Sf)));
Fh1 = a.*b.*G;

K = 10000000000;
b = ((H.*conj(H))./((H.*conj(H)) + K));
Fh2 = a.*b.*G;

//scf(), imshow(FTimg, [])
scf(), imshow(abs(ifft(G)), [])
scf(), imshow(abs(ifft(Fh1)), [])
scf(), imshow(abs(ifft(Fh2)), [])
//scf(), imshow(img, [])

Activity 18 - Noise Models and basic Image Restoration

This activity allows one to digitally restore an image or by improving its appearance through a procedure of restoration and a priori knowledge or known mathematical model of degradation applied to the image. In this case, different types of noise are added on the original image.

There are different types of noise presented in this activity. These noise are experimentally obtained from the variations of temperature in the experimental area, light source flickering and also during the transmission of data through wireless connections where atmospheric disturbances may happen. Random variables comprise noise and are characterized by probability density function (PDF as discussed in Activity 4). Most common noise include:

• Gaussian
• Rayleigh
• Gamma
• Exponential
• Uniform
• Impulse (Salt and Pepper)

To restore images with these added noise, a very basic method to use is to filter in spatial domain. Filtering methods employed in this activity are as follows:

• Arithmetic
• Geometric
• Harmonic
• ContraHarmonic

where S(x,y) represent the set of coordinates in a rectangular subimage window of size m x n with center at x,y. Thus, the window S(x,y) moves through the image and proper restoration process occurs. For the ContraHarmonic Filter, Q is the order of filter with Q>0 eliminating pepper noise and Q<0 eliminating salt noise;however, it can't do both. Below is the test pattern:

RESULTS
Noise is being added to the original image. For a kind of noise added, all forms of filtering are imposed.

It can be observed that filtering the noisy images is a difficult task to do. The restored images does not contain sharp boundaries between the grayscale values. Also as seen from the PDF's of the restored images, there is really a distortion that has happened. Among the filters, the harmonic filter provides the best result except for the pepper noise. It best eliminates salt noise. The Q value used for the contraharmonic filter is 1. It better eliminates pepper than salt noise. As observed, geometric filter also does not work best for images with pepper noise.

From left to right:

original image+noise, arithmetic filtering, geometric, harmonic, contraharmonic

Gaussian noise

Rayleigh Noise



Gamma Noise



Exponential Noise



Uniform Noise



Salt Noise

Pepper noise

OTHER images

I also applied the filtering to other types of images.

Gaussian

Rayleigh

Gamma

Exponential

Uniform

Salt

Pepper

I give myself 10 pts for finishing and understanding the activity. :)

References
M. Soriano. A18 – Noise models and basic image restoration.

Scilab Code

chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act18')

function [mat] = res(arr);
mat = matrix(arr,1,-1);
endfunction

function [hist] = pdf(arr);
mx = max(arr);
hist = [];
for i = 1:mx
f = find(arr == i);
hist(i) = length(f);
end
endfunction


// constant pattern of three levels of grayscale increment
//sz = 127;
//x = -127:127;
//[X,Y] = meshgrid(x);
//s1 = (X-32)*Y;
//s2 = (X+32)*Y;
//szx = size(X,1); szy = size(Y,1); sz = szx*szy;
//s = zeros(szx, szy) + 1;
//w = 100; s(127 - w:127+w,127-w:127+w) = 128;
//w = 50; s(127 - w:127+w,127-w:127+w) = 256;
//scf(), imshow(s, [])

//
//s = gray_imread('pattern.bmp')*255 + 1;
s = gray_imread('people1.bmp')*255 + 1;
//s = gray_imread('wiki1.bmp')*255 + 1;

szx = size(s,1); szy = size(s,1); sz = szx*szy;

// noise

//Y = grand(szx, szy, 'nor', 0.1, 0.1)*255 + 1; // gauss
//Y = grand(szx, szy, 'gam', 20, 100)*255 + 1; // gamma
//Y = grand(szx, szy, 'exp', 0.1)*255 + 1; // exponential
//Y = matrix(genrayl(0.05,sz), szx, szy)*255 + 1; // rayleigh
//Y = grand(szx, szy, 'unf', 0, 0.1)*255 + 1;

//
if min(Y) < 0
Y = abs(min(Y)) + Y + 1;
end

S = s + Y;
S = round((S/max(S)) * 256);


// salt and pepper
d = 0.009;
//S = round(imnoise(s, 'salt & pepper', d));

// salt only
solt = s;
sz = size(s,1)*size(s,2);
bits = round(d*sz);
solt(round(rand(bits,1)*(sz-1)+1)) = 256;
//S = round(solt);

// pepper only
peper = s;
peper(round(rand(bits,1)*(sz-1)+1)) = 1;
S = round(peper);


// histogram
//scf(), plot(pdf(s))
//scf(), plot(pdf(S))
scf(), imshow(S, [])


// arithmetic mean filtering
filtered_ari = zeros(szx,szy);
filtered_geo = zeros(szx,szy);
filtered_har = zeros(szx,szy);
filtered_con = zeros(szx,szy);
Q = 1;
win = 3;

for i = 1:szx
for j = 1:szy
iwm = i - win;
iwp = i + win;
jwm = j - win;
jwp = j + win;
sxy = 0;
su = 0;
pr = 0;
hr = 0;
cr = 0;
fac = 1E-9;
//[i j]
if iwm <= 0 & jwm <= 0 // upper left corner
sxy = cat(2, res(S(1:i, 1:j)), res(S(i:iwp, 1:j)), res(S(1:i, j:jwp)), res(S(i:iwp,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);
//su = mean( mean(S(1:i, 1:j)) + mean(S(i:iwp, 1:j)) + mean(S(1:i, j:jwp)) + mean(S(i:iwp,j:jwp)) );

elseif iwm > 0 & iwp <= szx & jwm <= 0 // left column
sxy = cat(2, res(S(iwm:i,1:j)) , res(S(i:iwp,1:j)) , res(S(iwm:i,j:jwp)) , res(S(i:iwp,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif iwm > 0 & iwp <= szx & jwp >= szx // right column
jj = abs(j-szx);
sxy = cat(2, res(S(iwm:i,jwm:j)) , res(S(i:iwp,j:j+jj)) , res(S(iwm:i,jwm:j)) , res(S(i:iwp,j:j+jj))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif iwp > szy & jwm <= 0 // lower left corner
ii = abs(i-szy);
sxy = cat(2, res(S(iwm:i,1:j)), res(S(i:i+ii,1:j)) , res(S(iwm:i,j:jwp)) , res(S(i:i+ii,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif jwm > 0 & jwp <= szy & iwm <=0 // upper row
sxy = cat(2, res(S(1:i,jwm:j)) , res(S(i:iwp,jwm:j)) , res(S(1:i,j:jwp)) , res(S(i:iwp,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif jwm > 0 & jwp <= szy & iwp >= szy // bottom row
ii = abs(i-szy);
sxy = cat(2, res(S(iwm:i,jwm:j)) , res(S(i:i+ii,jwm:j)) , res(S(iwm:i,j:jwp)) , res(S(i:i+ii,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif iwm <= 0 & jwp > szx // upper right corner
jj = abs(j-szx);
sxy = cat(2, res(S(1:i,jwm:j)) , res(S(i:iwp,j:j+jj)) , res(S(1:i,jwm:j)) , res(S(i:iwp,j:j+jj))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);


elseif iwp > szy & jwp > szx // lower right corner
ii = abs(i-szy);
jj = abs(j-szx);
sxy = cat(2, res(S(iwm:i,jwm:j)) , res(S(i:i+ii,jwm:j)) , res(S(iwm:i,jwm:j+jj)) , res(S(i:i+ii,j:j+jj))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

elseif i > win & i <> win & j < szx-win // no problem region
sxy = cat(2, res(S(iwm:i,jwm:j)) , res(S(i:iwp,jwm:j)) , res(S(iwm:i,j:jwp)) , res(S(i:iwp,j:jwp))) + fac;
su = mean( sxy );
pr = max(cumprod(sxy))^(1/size(sxy,2));
hr = size(sxy,2) / sum(sxy.^(-1));
cr = sum(sxy).^(Q+1) ./ sum(sxy).^(Q);

end
filtered_ari(i,j) = su;
filtered_geo(i,j) = pr;
filtered_har(i,j) = hr;
filtered_con(i,j) = cr;
end
end

//scf(), imshow(s, [])
//scf(), imshow(S, [])
scf(), imshow(filtered_ari, [])
scf(), imshow(real(filtered_geo), [])
scf(), imshow(real(filtered_har), [])
scf(), imshow(real(filtered_con), [])

//scf(), plot(pdf(round(filtered_ari)))
//scf(), plot(pdf(round(filtered_geo)))
//scf(), plot(pdf(round(filtered_har)))
//scf(), plot(pdf(round(filtered_con)))
//

Activity 17 - Photometric Stereo

Photometric stereo is a technique for 3D modeling which estimates the surface contour with the cues given by the shadows of the object, taken at one viewpoint under different illumination directions. Technically, method computes for the surface normals of the object and from there develops the height variations of the object.

In this activity, synthetic photometric stereo images are used as shown below.

Using these images, we can calculate for the surface information as encoded in the shadings defined by some matrix:

with each column corresponding to x,y,z component of the source. Having multiple images allow us to compute for g given by:

After that, we take the normal vector by dividing g by its length and relating the surface normal to height variations give us:

.

The V's are given as:

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

Doing all the rigorous mathematical derivations and implementing algorithm gives us a 3D reconstruction of the synthetic object.

Note that the reconstruction is not a perfectly hemispherical surface. This is maybe due to the black and white boundary between quadrants of the hemisphere.

I give myself 10 pts for understanding and doing the activity correctly.

Reference
A17 - Photometric Stereo handout. M. Soriano, 2008.

Code
chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act17');

loadmatfile('Photos.mat')
//who // list of variables
//scf(), imshow(I1, [])
//scf(), imshow(I2, [])
//scf(), imshow(I3, [])
//scf(), imshow(I4, [])


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

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

g = inv(V'*V)*V'*I;
lg = sqrt(g(1,:).^2 + g(2,:).^2 + g(3,:).^2);
for i = 1:size(g, 1), n(i,:) = g(i,:)./(lg + 1E-9), end

dx = -n(1,:) ./ (n(3,:) + 1E-9);
dy = -n(2,:) ./ (n(3,:) + 1E-9);
sx = size(I1,1); sy = size(I1,2);

z = cumsum(matrix(dx,sx,sy),2) + cumsum(matrix(dy,sx,sy),1);

scf(), plot3d(1:sx,1:sy,z)

Activity 16 - Artificial Neural Network (ANN)

Artificial Neural Network employs the ability of the brain to learn and classify objects that is fed on to the system. The idea is that the algorithm learns by itself from the examples or training sets that are given to it. Once the network learns from the examples, it can easily identify classes that are presented to it with higher processing speed.

The input layer is where the features are fed and the hidden layer process these features by having interconnected neurons after the preceding input layer. Then the output layer integrates the learning process and gives the learned classification.

Detailed explanation and examples are presented to Cole's Blog.

In my previous activities, I classified the 25c and P1 coins using Minimum Distance Classification and Linear Discriminant Analysis and found out that LDA presents better classification accuracy compared to MDC. In this activity, I'll try to compare ANN classification with the previous algorithm.

Again, the rg values of my samples are the features used to classify the classes. Here are some sample images of my objects:

To analyze the results, I preferred to present them visually. Here are the results of my ANN:

Black = 0 and white = 1. The left image is the input training class. Black means that it is of class P1 coin and white of class 25c coin. The black boundary is meant for visualizing the white part. The right image is the output classification. Grayscaling is observed since the output is not an integer. The outputs are close to 0 and 1. This set presents an arranged training classes for the algorithm to learn. That is, ANN learns the P1 coin features first an next are the 25c coin features. This training style exhibits 100% classification accuracy. Next is we determine if the arrangement of training classes affect the accuracy of the ANN.

The left image is the input training style and the output is at the right. This presents that even if we scramble the training class style, a 100% accuracy classification is still obtained. Thus, ANN presents a more flexible way of classification compared to LDA and MDC. Also, no stricter assumptions are carried into the algorithm unlike LDA which assumes that objects are linearly separable. Also, this algorithm trains and learns itself an accurate way of classification.

I give myself 10 pts for understanding and doing the activity correctly.

Reference:
A16 - Neural Networks handout. M. Soriano, 2008.

Acknowledgement:
I acknowledge M. Sison for the gift I received from him in debugging the problem for ANN toolbox. ;)

Code:
chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act14\' + 'data');

x1 = fscanfMat('piso_t.txt');
x2 = fscanfMat('ben_t.txt');
fname = 'piso';
tst1 = fscanfMat(fname + '_t.txt');
fname = 'ben';
tst2 = fscanfMat(fname + '_t.txt');

rand('seed',0);
N = [2,2,1];

//x = cat(1, x1(:,1:2), x2(:,1:2))';

//set 1
x = cat(1, tst1(:,1:2), tst2(:,1:2))';
t = [0 0 0 0 0 1 1 1 1 1];

//set 2
x=cat(1,tst1(1:2,1:2),tst2(1:2,1:2),tst1(3:4,1:2),tst2(3:4,1:2),tst1(5,1:2),tst2(5,1:2))';
t = [0 0 1 1 0 0 1 1 0 1];

lp = [4, 0];
W = ann_FF_init(N);
T = 1000;
W = ann_FF_Std_online(x,t,N,W,lp,T);

//results
c = ann_FF_run(x,N,W)';
t = mtlb_repmat(t', [1,50]); //input
c = mtlb_repmat(c, [1 50]); //output
//boundary
c(:,1:2) = 0; c(:,49:50) = 0; t(:,1:2) = 0; t(:,49:50) = 0;
scf(), imshow(t, [])
scf(), imshow(c, [])

Activity 15 - Probabilistic Classification (Linear Discriminant Analysis)

LDA 2-feature
This is another type of classification in which the probability of class belonging is calculated and is determined through early observations of trained data sets. The assumption in this technique is that the classes are linearly separable through linear combinations of features that best distinguish the classes. In this activity, two features are used to separate the P1 and 25c coins in the previous activity. The features used are the rg values, similar to Activity 14.

For more detailed explanation of LDA visit this site.

Shown above is the original plot of the features in rg space. We can see that the data is linearly separable but we need to find the line to separate them and rotate them to see clearly the separability of the samples.

Applying the LDA algorithm, we obtained the following result:

As seen, the LDA separates the P1 and 25c coins very closely. The algorithm obtained 100% or 10/10 classification based on the two features.

This makes me conclude that LDA is better than Minimum distance classification since it has higher classification accuracy.

I give myself 10 pts for understanding and doing the activity correctly.

Reference:
A15 - Probabilistic Classification handout. M. Soriano. 2008.

Acknowledgement
I acknowledge Jaya for the help in understanding LDA algorithm. ;)

Code:
LDA
chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act14\' + 'data');
fname = 'ben';
x1 = fscanfMat('piso_t.txt');
x2 = fscanfMat('ben_t.txt');
tst = fscanfMat(fname + '_s.txt');

x = [cat(1,x1(:,1),x2(:,1)) cat(1,x1(:,2),x2(:,2)) cat(1,x1(:,5),x2(:,5))];

u1 = mean([x1(:,1) x1(:,2) x1(:,5)], 'r');
u2 = mean([x2(:,1) x2(:,2) x2(:,5)], 'r');

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

x1c = [x1(:,1) - u(1) x1(:,2) - u(2) x1(:,5) - u(3)];
x2c = [x2(:,1) - u(1) x2(:,2) - u(2) x2(:,5) - u(3)];

n1 = size(x1c,1);
n2 = size(x2c,1);
n = size(x,1);

c1 = (x1c'*x1c) / n1;
c2 = (x2c'*x2c) / n2;

C = [(n1/n)*c1(1)+(n2/n)*c2(1) (n1/n)*c1(2)+(n2/n)*c2(2) (n1/n)*c1(3)+(n2/n)*c2(3);...
(n1/n)*c1(4)+(n2/n)*c2(4) (n1/n)*c1(5)+(n2/n)*c2(5) (n1/n)*c1(6)+(n2/n)*c2(6);...
(n1/n)*c1(7)+(n2/n)*c2(7) (n1/n)*c1(8)+(n2/n)*c2(8) (n1/n)*c1(9)+(n2/n)*c2(9)];

h = [];
hp = [];
for tst_num = 1:size(tst, 1)
xk = cat(2, tst(tst_num,1:2), tst(tst_num,5));
f1 = u1*inv(C)*xk' - 0.5*u1*inv(C)*u1' + log(n1/n);
f2 = u2*inv(C)*xk' - 0.5*u2*inv(C)*u2' + log(n2/n);
hp = cat(1, hp, [f1 f2]);
h = cat(1, h, [1*(f1>f2) 1*(f1
end

chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act15');
fprintfMat(fname + '_res3' + '.txt', hp, '%0.5f');

PLOTTING
chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act14\data');
ben_t = fscanfMat('ben_t.txt');
piso_t = fscanfMat('piso_t.txt');
scf(), plot(ben_t(:,1), ben_t(:,2), 'sb', 'markerface', 'b')
plot(piso_t(:,1), piso_t(:,2), 'sg', 'markerface', 'g')
xlabel(' r')
ylabel(' g')
legend('25c', 'P1',2)

chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act15');
ben_res2 = fscanfMat('ben_res2.txt');
ben_res3 = fscanfMat('ben_res3.txt');
ben_tres2 = fscanfMat('ben_tres2.txt');

piso_res2 = fscanfMat('piso_res2.txt');
piso_res3 = fscanfMat('piso_res3.txt');
piso_tres2 = fscanfMat('piso_tres2.txt');

x = 200:0.001:300;
scf(), plot(ben_tres2(:,1), ben_tres2(:,2), 'ob', 'markerface', 'b')
plot(ben_res2(:,1), ben_res2(:,2), 'sy', 'markerface', 'y')
plot(piso_tres2(:,1), piso_tres2(:,2), 'og', 'markerface', 'g')
plot(piso_res2(:,1), piso_res2(:,2), 'sm', 'markerface', 'm')
plot(x,x,'-r')

legend('train 25c', 'sample 25c','train P1', 'sample P1', 2)
xlabel('f1')
ylabel('f2')
//mtlb_axis([100 400 100 400])

Activity 14 - Pattern Recognition

In this activity, we are to classify a set of images into its corresponding classes. A class is a type of an object or group of objects (eg. P1 coin, 25 c coin, flower, etc.). To recognize a sample object with unknown class, we need to do pattern recognition. Minimum distance classification is used in this activity to classify the objects. Basically, a training set of objects is fed into the algorithm and patterns are extracted from these training set. To specify the patterns, features are selected. These features are the patterns used to determine whether a sample object belongs to a kind of class using minimum distance classification.

Three features are used: normalized r, g and per^2 / area to identify the class belonging of a certain object. If the distance of an object's feature vector from the mean of a certain class feature space is at minimum, then that object belongs to that class. This is governed by the Euclidean distance:

.

Equating Dj = y, we shall compute the distance using the equation below. The distance to which the class is minimum, that sample object belongs to that class.

Basically, the closeness of a sample feature vector to the mean of a certain class feature vector is determined by obtaining the minimum distance. To apply this, the original image taken is presented to the left. Then, white balancing is applied to the image (right).

The objects are named as piso, ben, leafe, leafs and santan with 5 classes in total. Half of the set of classes are used as training and half are used as samples. Again, the features r and g are obtained by taking a part of the object and the normalized red and green computed (A12). Then, thresholding is used to segment the ROI and determine the perimeter using follow of scilab and the total number of pixels are calculated. The area is calculated using the morphological operations (A9). I tried to use only as third feature the area or perimeter, but there has been no consistency in the results and the feature is shape and size invariant, that's why I choose to use the ratio of per^2 / area since it is a dimensionless unit.

I also tried to use the eccentricity of the region to as feature since it is invariant in the capture angle of images but I've got a problem to non-regular polygons because I don't know where to obtain the major and minor axes. That's why I opt to use other features such as perimeter and area ratio.

Shown below is the plot of the results in rg axes (left). It can be observed that the graph follows that of the rg Normalized chromaticity space as shown to its right. The 3D graph of the distribution of mean of training and samples feature vectors in the [r g per^2/area] coordinate system is presented right below of the two plots. It is interesting to note that the P1 and 25c coins have almost equal ratio of per^2/area since for a circle that ratio is equal to 4pi. But for the sake of non-circular objects, this feature is very distinct. Thus, the two coins present errors in the object classification. As seen also, the cluster of the 25c and P1 sample data are very close to each other. This is because of the dilapidated texture of the coins. The rg plot is presented to exemplify that these features describe the normalized rg values fo each classes correctly.

The algorithm presents 84% accuracy in object classification. The errors obtained are incurred from the P1 and 25c samples classification missing 2 correct classification for both P1 and 25c samples. In the next activity using Linear Discriminant Analysis, the 25 c and P1 data is expected to be classified at higher accuracy.

I give myself 9 points for finishing and understanding the activity but failed to finish early.

Acknowledgments:
I acknowledge Jaya and Orly for bringing up intellectual discussions about the algo.

CAVEAT
Be careful to put '.' during element division as shown in red texts in the code below.

Code:

Scilab (image processing)

function [img] = imopen(img_bw, se);
img = dilate(erode(img_bw, se), se);
endfunction

function [img] = imclose(img_bw, se);
img = erode(dilate(img_bw, se), se);
endfunction

// intialize
folder = string('piso sample');
chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act14\' + folder);
data = [];
fname = string('piso_s');
n = 5;

// rg feature
for i = 1:n
imgc = imread(fname + string(i) + 'c.jpg');
//imgc = imread('piso_t1c.jpg');
R = imgc(:,:,1);
G = imgc(:,:,2);
B = imgc(:,:,3);
I = R + G + B;
r(i) = mean(R ./ I);
g(i) = mean(G ./ I);
end


// area feature
for i = 1:n
img = imread(fname + string(i) + '.jpg');
//img = imread('piso_t3.jpg');
img = img(:,:,3); // blue channel
bw = abs(im2bw(img, 0.65) - 1); // 0.65 // 0.55
se1 = [1 1 1; 1 1 1;1 1 1];
imgb_m = imclose(imopen(bw, se1), se1);
scf(), imshow(imgb_m, [])
ma(i) = sum(imgb_m);
[t, y] = follow(imgb_m);
//scf(), plot2d(t,y)
per(i) = size(t, 1);
// eccentricity
//[row, col] = find(imgb_m == 1);
// vert = max(row) - min(row);
// horz = max(col) - min(col);
// if vert > horz, e = sqrt(vert^2 - horz^2) / vert^2;
// else e = sqrt(horz^2 - vert^2 ) / horz^2;
// end
// ecc(i) = e + 1;
end

data(:,1) = r;
data(:,2) = g;
data(:,3) = per;
data = cat(2, data, ma);
data = cat(2, data, per.^2 ./ ma);

chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act14\' + 'data');
fprintfMat(fname + '.txt', data, '%0.5f');

Matlab (plotting and classification)
%%
load ben_s.txt
load ben_t.txt
load leafe_s.txt
load leafe_t.txt
load leafs_s.txt
load leafs_t.txt
load piso_s.txt
load piso_t.txt
load santan_s.txt
load santan_t.txt

%%
x = 1; y = 2; z = 5;
hold on
plot3(mean(ben_t(:,x)), mean(ben_t(:,y)), mean(ben_t(:,z)), 'bo', 'MarkerFaceColor', 'k')
plot3(mean(leafe_t(:,x)), mean(leafe_t(:,y)), mean(leafe_t(:,z)), 'go', 'MarkerFaceColor', 'k')
plot3(mean(leafs_t(:,x)), mean(leafs_t(:,y)), mean(leafs_t(:,z)), 'yo', 'MarkerFaceColor', 'k')
plot3(mean(santan_t(:,x)), mean(santan_t(:,y)), mean(santan_t(:,z)), 'ro', 'MarkerFaceColor', 'k')
plot3(mean(piso_t(:,x)), mean(piso_t(:,y)), mean(piso_t(:,z)), 'mo', 'MarkerFaceColor', 'k')
xlabel('r'); ylabel('g'); zlabel('per^2 / area');

%%
plot3(ben_s(:,x), ben_s(:,y), ben_s(:,z), 'bs', 'MarkerFaceColor', 'b'), hold on
plot3(leafe_s(:,x), leafe_s(:,y), leafe_s(:,z), 'gs', 'MarkerFaceColor', 'g')
plot3(leafs_s(:,x), leafs_s(:,y), leafs_s(:,z), 'ys', 'MarkerFaceColor', 'y')
plot3(santan_s(:,x), santan_s(:,y), santan_s(:,z), 'rs', 'MarkerFaceColor', 'r')
plot3(piso_s(:,x), piso_s(:,y), piso_s(:,z), 'ms', 'MarkerFaceColor', 'm')
grid on;
% xlabel('r'); ylabel('g'); zlabel('per^2 / area');

%%
m1 = mean(cat(2, ben_t(:,1:2), ben_t(:,5)), 1);
m2 = mean(cat(2, leafe_t(:,1:2), leafe_t(:,5)), 1);
m3 = mean(cat(2, leafs_t(:,1:2), leafs_t(:,5)), 1);
m4 = mean(cat(2, santan_t(:,1:2), santan_t(:,5)), 1);
m5 = mean(cat(2, piso_t(:,1:2), piso_t(:,5)), 1);

%%
x1 = cat(2, ben_s(:,1:2), ben_s(:,5));
x2 = cat(2, leafe_s(:,1:2), leafe_s(:,5));
x3 = cat(2, leafs_s(:,1:2), leafs_s(:,5));
x4 = cat(2, santan_s(:,1:2), santan_s(:,5));
x5 = cat(2, piso_s(:,1:2), piso_s(:,5));

%% minimum distance classification
y = [];
h = [];
for i = 1:5 % classes
for j = 1:5 % samples
for k = 1:5 % mean classes
mm = eval(['m' int2str(k)]);
xx = eval(['x' int2str(i)]);
y = abs(xx(j,:) - mm);
z(k) = abs(sqrt(y*y'));
end
h = cat(1, h, [find(z == min(z)) i]);
end
end

Activity 13 - Correcting Geometric Distortion

In this activity, we correct an image with geometric distortion. We used an image that has grid patterns in the image. This is to determine the kind of distortion present in the image and correct the image by making curve lines make straight. There are two kinds of distortions: 1. pincushion and 2. barrel distortion. Barrel distortion introduces bulging effect in the images while pincushion distortion introduces inward bending of lines. To correct the image, a regularly occurring pattern is required to be part of the image to be able to visualize the distortion effect.

In the correcting process, a transformation of the pixel coordinates and interpolation of gray level values are required to complete the correction. Pixel coordinate transformation allows one to straighten the coordinates and gray level interpolation to visualize the transformation. It can be expressed in matrix notation,

and solving for the coefficients,

.

After which, we can now compute for the undistorted pixel coordinates using:

.

Then, we project the grayscale values from the original image to the undistorted using:

where the coefficients can be solved similar to the matrix notation above.

By generating an ideal grid of known vertex coordinates, we can proceed to pixel coordinates transformation.

Shown to the left is the original image with barrel distortion. At the right is the ideal rectangle vertices and below them is the corrected image. Note that there are some pixels that are not interpolated in the corrected image.

I give myself 9 pts for understanding the activity but blogging it late.

Reference:
A13 - Correcting Geometric distortion. M. Soriano. 2008.

Acknowledgment:
I acknowledge the help of Jaya in understanding the notations in the procedure.

Code:
chdir('E:\Documents and Settings\vergara\My Documents\acads\1st sem 09-10\186\act13');
//165, 92 //182, 92 //165, 109 //182, 109
stacksize(50000000)
img = gray_imread('barrel1.png');

// for visualization of ideal rect
sx = size(img,1); sy = size(img,2);
z = zeros(sx, sy);
X = 1:110:sx; Y = 1:110:sy; //17
//z(X, Y) = 1;
//imwrite(z, 'ideal_rect.jpg');
//scf(), imshow(z, [])


// locate real vertices
pts = 25; //25, 49, 81
scf(), imshow(img, [])
l = locate(pts); l = l';
l = [l(pts:-1:1,2) l(pts:-1:1,1)];

fprintfMat('barrel.txt', l, '%0.2f');
l = fscanfMat('barrel.txt');

// for visualization of real rect
z1 = zeros(sx, sy);
for i = 1:pts
z1(l(i,1), l(i,2)) = 1;
end
scf(), imshow(z1, [])
//imwrite(z1, 'vertices_rect.jpg');

// reshaping of data
ideal = [X' mtlb_repmat(Y(size(Y,2))',size(Y,2),1)];
rreal = l;
for i = size(X,2) - 1:-1:1
ideal = cat(1,ideal,[X' mtlb_repmat(Y(i)',size(X,2),1)])
end

// test
for j = 1:4// row
img(ix(j,1), iy(j,1)) = 1;
end

//reshape again
rx = matrix(rreal(:,1),size(X,2),size(Y,2)); rx = rx(:,size(Y,2):-1:1);
ry = matrix(rreal(:,2),size(X,2),size(Y,2)); ry = ry(:,size(Y,2):-1:1);

ix = matrix(ideal(:,1),size(X,2),size(Y,2)); ix = ix(:,size(Y,2):-1:1);
iy = matrix(ideal(:,2),size(X,2),size(Y,2)); iy = iy(:,size(Y,2):-1:1);

//
undist = zeros(sx, sy);
for i = 1:size(X,2) - 1
for j = 1:size(Y,2) - 1
T = [ix(i,j) iy(i,j) ix(i,j)*iy(i,j) 1;...
ix(i,j+1) iy(i,j+1) ix(i,j+1)*iy(i,j+1) 1;...
ix(i+1,j) iy(i+1,j) ix(i+1,j)*iy(i+1,j) 1;...
ix(i+1,j+1) iy(i+1,j+1) ix(i+1,j+1)*iy(i+1,j+1) 1];
Xx = [rx(i,j); rx(i,j+1); rx(i+1,j); rx(i+1,j+1)];
Yy = [ry(i,j); ry(i,j+1); ry(i+1,j); ry(i+1,j+1)];
C14 = inv(T)*Xx;
C58 = inv(T)*Yy;

for kx = ix(i,j):ix(i+1,j) - 1
for ky = iy(i,j):iy(i,j+1) - 1
xh = C14(1)*kx + C14(2)*ky + C14(3)*kx*ky + C14(4);
yh = C58(1)*kx + C58(2)*ky + C58(3)*kx*ky + C58(4);
// not this easy, solve for the 4 unknowns and 4 equations
if modulo(xh, round(xh)) == 0 & modulo(yh, round(yh)) == 0 // integer coordinates
undist(round(xh), round(yh)) = img(round(xh), round(yh));
else // non-integer coordinates
if xh == 1 | xh == sx | yh == 1 | yh == sy // at boundaries
undist(xh,yh) = 0;
else // not at boundaries
U = [xh-1 yh-1 (xh-1)*(yh-1) 1;...
xh-1 yh+1 (xh-1)*(yh+1) 1;...
xh+1 yh-1 (xh+1)*(yh-1) 1;...
xh+1 yh+1 (xh+1)*(yh+1) 1];
V = [img(xh-1, yh-1) img(xh-1, yh+1) img(xh+1, yh-1) img(xh+1, yh+1)]';
A = inv(U)*V;
undist(xh,yh) = A(1)*xh + A(2)*yh + A(3)*xh*yh + A(4);
end
end
end
end
end
end

scf(), imshow(undist, []);