certified dormer

Activity 19 - Restoration of Blurred image

2009-10-07T09:07:00.001-07:00

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:

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

2009-10-06T08:03:00.000-07:00

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

2009-09-13T09:19:00.000-07:00

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)

2009-09-13T06:25:00.000-07:00

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)

2009-09-12T08:10:00.001-07:00

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

2009-09-09T09:27:00.000-07:00

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

2009-09-02T17:53:00.000-07:00

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, []);

Activity 12 – Color Image Segmentation

2009-08-06T23:21:00.000-07:00

During the previous activities, we have been separating the regions of interest using histogram manipulation. In this activity, ROI separation is obtained using color. Grayscale value is not enough since it only provides a two-dimension space to select unique features from the scene. Using color, it provides three representations to provide unique features of a scene. However, colored objects will have variations in brightness or in shading that is RGB is not enough to represent color distinguish ROI. In this activity, the Normalized Chromaticity Coordinates (NCC) is being introduced to separate color and brightness information.

With this kind of representation, we can now present 3D to 2D coordinate system only.

There are two ways of obtaining the ROI using NCC. One is the Parametric probability distribution where this method calculate the probability that a certain color belongs to the ROI while in the Non-parametric probability distribution, it is mainly concerned on the histogram of the ROI and using the idea back in Activity 4, pixel per pixel values are backprojected from the rg histogram to the image. The first technique is purely computational while the latter involves ‘manual’ pixel per pixel replacement.

Above is sample ROI segmentation using parametric estimation. Note that the flag image doesn’t have a shading variation in its color so it is expected that this image will be segmented accurately. A yellow patch is obtained for sample segmentation. Even if non-parametric estimation is used, there will be no change in the quality of segmentation produced.

Displayed above are ROI segmentation results using parametric and non-parametric estimates. Beside the original image is the patch that is used to select the ROI. Parametric estimate is the one that is displayed to the left and for the non-parametric is the one displayed to the right. It can be observed that parametric estimates have a poor segmentation quality as compared to the non-parametric. This is because only a probability of the color information is obtained. From the formula, the method computes for the mean and the standard deviation of a color from the image and it is compared to the probability that it has the same statistics as that of the ROI. The non-parametric estimate then provide better segmentation results because it follows a look-up table process and is one-to-one pixel per pixel replacement.

Reply to comment:

Sorry mam, I failed to rate my work. I give myself 10 points for finishing and understanding this activity.

Activity 11 – Color Image Processing

2009-08-06T22:51:00.000-07:00

This activity is mainly about proper white balancing of images under certain lighting conditions and also about understanding how incorrect white balanced images be corrected. Color in images is represented by the RGB channels. Mathematically, the color perceived by the camera is the sum of the product of the object reflectance, source illumination and the camera spectral sensitivity channels for RGB. It can be presented using:

where

K is the white balancing constant and all of these parameters constitute to the image taken by the camera.

Shown above are four images with the two outmost images as white patches corresponding to the near images. Each white patch is used to do automatic white balancing in the two images. These two images are different because they are taken at different NA and shutter speed values that is why the image at the left looks darker.

Two most popular ways of obtaining auto white balancing is used. The first is the White Patch algorithm where the mean of the RGB channel values of the white patch is used to divide the corresponding RGB channel of the incorrectly-white balanced object. For the Gray Patch method we need not use the white patch instead we use the mean of the RGB channels of the image itself and use this to divide the RGB channel of the image.

Presented above are the results of the White Patch (left) and Gray World (right) method. Both for the two sample images, by the way image is wrongly white balanced because the illumination is flourescent light while the white balancing setting of the camera is sunny, the White Patch algorithm produces a better white balancing estimate than the Gray World algorithm. For the White patch recon, one can observe good rendering of colors. Green appears green and so with other colors. In the Gray World recon, there is an observed saturation of colors with white light. The brightness in the latter method is greater than the other.

These images are taken at the same lighting condition as above but the white balancing setting is incandescent with corresponding white patch taken.

Similar to the result for the sunny white balancing setting, White Patch recon is visually better than Gray World. It is observed however that in the White Patch recons, there is an increased saturation of the color. The maroon table for example becomes a striking color maroon and the red hue becomes intense red especially to the image taken at high shutter speed.

The images above are taken at incandescent white balancing setting and the illumination is still flourescent lamp. Again, these two images differ in the NA and shutter speed settings. Objects of interest are the different shades of blue as represented by those different objects. What we want is to obtain correct white balancing using the two methods.

Left images correspond to White Patch rendering algorithm while right images for Gray World. It can be clearly seen that shades of blue color are better rendered for White rather than for Gray. This is due to the incorrect blue shade rendering for Gray. This is because the Gray World algorithm uses the mean of the RGB which in turns averages the whole channel. White patch algo refers to the white patch itself to properly white balance an image. This provides better scaling value because it is based on the mean of the white patch. White Patch method then assures that white will look white when white balanced.

I give myself 9 points for finishing and understanding the activity since I got images of saturated intensities.

Activity 10 – Preprocessing Text

2009-08-06T17:48:00.000-07:00

This activity includes usage of techniques to pre-process an image. Specifically, this is a real-world problem to extract handwritten text from a sample that presents challenges in cleaning object of interest.

The image to the left is the image to be processed. It can be noticed that the image is tilted based from the horizontal lines. Before we process this image, we must rotate it to make horizontal lines horizontal and ease filtering. We used the function mogrify in scilab. The image is tilted at 1.2 degrees. The tilted image is shown at the right.

Next is to crop desired object. I cropped the part where there are horizontal lines, handwritten text and printed text. After cropping, I filtered the horizontal lines in the image using the mask below from its Fourier transform. Then, I binarized the image and apply morphological operations such as opening, closing, erode,dilate and thin. Below are the processed images.

It can be seen from above that after binarizing the image, one can still read the handwritten text VGA Cable and printed text DESCRIPTION. However after applying closing, then thin operations of SE = [1 1]', the one-pixel thick text can't be read. I then chose other handwritten texts to process. Below is another processed images with the same method as above.

Also, I chose a printed text to convince myself that I can digitally extract texts from challenging samples.

I can say that I have extracted the handwritten text DEMO III from above although there is really no rigorous cleaning of the image involved. That is, if I try to clean the image using FT filtering, some parts of my handwritten text can be removed which makes me difficult to read the text after morphological operations. But the printed text, after FT filtering, still can be read. Although some parts or letters are distorted due to filtering of the horizontal lines. Maybe if I can just obtain a good mask and avoid filtering information from the text itself, I may have better results.

Also, correlation is presented in this activity. I used the word DESCRIPTION and try to see if it will produce good correlations in the whole image. The left image is the binarized pattern and to the left is the correlation map. To zoom in, cllick the image. There are white spots that can be noticed in the correlation map. These are the words corresponding to the pattern.

I give myself 9 points for finshing and understanding the activity although I did not get accurately the pre-processing of one pixel handwritten text.

Activity 9 – Binary Operations

2009-08-06T16:46:00.001-07:00

This activity uses histogram manipulation to separate ROIs in binary images and to estimate the areas of cell-like regions or blobs using closing and opening functions.

Shown below is a grayscale of the test image. This test image is then cropped into 12 parts and each subimage is binarized with proper threshold. The threshold value can be obtained using histogram function in GIMP.

Then, each binarized subimage is processed using closing and opening functions. Closing function is the erosion of a dilated image using the same structuring element and opening function is the dilation of an eroded image using the same SE. The SE used in this activity is se1 = [0 1 1 0 ;1 1 1 1;1 1 1 1;0 1 1 0]. The morphological process used is first to open the image then closed it. Opening function allows groups of blobs to be separated from other blobs and to remove smaller blobs (dots). Closing functions allows the filling of holes within blobs. Shown below is the image before morphological operation.

The next of images are obtained after morphological operations.

It can be noticed that the noise (dots) are removed and holes within blobs are filled. Also the edges of the blobs become well recognized. From these set of images, each blob is labeled using BWLABEL. This function labels every connected pixels from 0 - n blobs. Thus, to estimate the mean area of a blob, the areas are tallied and its histogram is plotted. The point where there will be the peak value corresponds to the number of individual cell with that area. This area of the cell will be the highest frequency area because it is dominant in the image.

From the histogram, I took the top two bins of highest frequency. From its range of areas, I calculated the mean and the standard deviation. I chose these top two bins because I don't have the idea what the ideal cell area is so it would be safe to choose these two bins. The range of these two bins are from 475 - 585. Thus the mean is 532.81 square px and the stddev is 18.96. Also I calculated the area of an individual cell and I obtained 536 sq px. The mean estimate is still within comparison range noting the value of the standard deviation.

I give myself 10 points for finishing and understanding the activity.

Activity 8 - Morphological Operations

2009-07-27T07:04:00.000-07:00

In this activity, basic morphological operation is introduced. The first operation is by erosion. In image erosion, images are being thinned according to the structuring element that is used. The second is image dilation. Dilation allows the image to cover more area based on the structuring element.

Before simulation, these to operations are solved manually. Shown below are the simulations for these two operations. First is erosion.

The first column is the original image, the second is for an SE = 4 x 4 ones, third is for 2 x 4 ones, fourth is for 4 x 2 ones and lastly is for a cross that is 5 pixels long and one pixel thick.

My predictions matched for the square, triangle and circle patterns. But for the hollow square and cross, my prediction are not as accurate. This is because I thought that there will no side lengths in the 2nd to last SE because i thought that the sides will be completely removed. For the cross, I didn't expect the last pattern for the cross SE.

Almost all of my predictions are correct in the Dilation of images. The hard one is on the pattern triangle where there will be no pointed corner.

I give myself 10 points for finishing and understanding the activity.

Activity 4 – Enhancement by Histogram Manipulation

2009-07-22T09:15:00.001-07:00

This activity involves the manipulation of image histograms to enhance the quality of low contrast images. The histogram of an image is equivalent to the Probability Distribution Function taking into account that the whole image is normalized with the total number of pixels. Taking the cumulative sum of this PDF, a CDF or Cumulative Distribution function can be obtained. Using this CDF, one can backproject the gray values of a poor-contrast image to this CDF and in turn enhancing the image.

Presented in this activity are the different types of CDF to enhance poor-contrast images. Below are the plots of these desired CDFs.

The blue plot represents the CDF of an original image, yellow corresponds to a linear CDF, green is from a non-linear CDF specifically an arctan function and red plot, also a non-linear CDF is a parabolic CDF.

Displayed below are the corresponding images of the CDFs above. The first image is the original image, 2nd is the enhanced image via backprojection of a linear CDF, 3rd is by using the arctan function and the last set of images are for the parabolic CDF.

It can be observed that there is a significant visual difference among the set of processed images. There is a good visual contrast enhancement for both the linear and arctan CDF application. It is observed that compared from the original image, there is good shading and is pleasing to the eye the reconstructions using these CDFs. However for the parabolic CDF, it seems that the image was saturated because this CDF increases the value of each pixel being backprojected as observed.

Shown above are the histograms of sample- set of images. From the orginal histogram, there is a concentration of low value pixels whereas when compared with the linear and arctan - enhanced images, there is a good distribution of pixel values both in the high and low pixel values. For the parabolic CDF, there is a lump of pixel value at the high value pixel range making the image look saturated.

When the CDF of histograms of these images are obtained, the result is the applied CDF in to the image.

Non-linear CDFs mimic that non-linear human eye response. The results are shown above where there are good enhancements using the arctan and linear CDFs rather than the parabolic CDF.

I give myself 9 points for finishing and understanding the activity but posted late.

References:
images are taken from

http://z.about.com/d/geology/1/0/u/E/1/sutbutaerialorig.jpg
http://static.photo.net/attachments/bboard/00F/00FzkU-29347284.jpg
http://myweb.lsbu.ac.uk/dirt/museum/margaret/08--452-1000120.jpg
http://www.dphclub.com/learn-photoshop/images/war-time-1.jpg

Activity 7 - Enhancement in the Frequency Domain

2009-07-20T10:25:00.001-07:00

In this activity, the idea of Fourier transform is being evaluated.

2 dots (left) modulus of FT(right)
Thus the FT of two dots is a sine pattern.

Increasing the size of the circle decreases the area where the sine pattern is observed. Thus, the smaller the size of the circles, the more the sine pattern will be observed.

The same is true with a 2 square aperture. The smaller the 2 square apertures are, the more the sine pattern (Bessel function) becomes observed.

In terms of the Gaussian aperture, there is a distinct circular area where the sine pattern is observed. Increasing the variance increases the area of the observed sine pattern.

The image to the left are the gaussian dots, center are the real parts and to the right are the imaginary parts. It can be observed that both the two Gaussian dots (inverted and not) have the same imaginary parts but different real parts. For this case, the apertures have the same odd signals since they have the same imaginary components.

Next is to enhance images using FT. This part is to enhance the fingerprint by removing blotches and enhancing the ridges. A mask is used to allow only the desired high frequency patterns for the fingerprints to be enhanced.

In this part, instead of allowing certain high frequency patterns, these high frequency patterns are blocked. This is to limit the repetition of the vertical lines in the image.

This is also true for the canvass weave filtering. This is a different mask used to filter out the canvasweave.

To visualize the canvas weave, just invert the mask pattern and the resulting image will be the canvasweave patterns.

In this activity, I learned how to filter and enhance images in the FT space. Also, the formation of certain masks to enhance the images is a skill to develop so that important features of the images will not be sacrificed and only the desired patterns will be filtered.

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

Activity 6 – Properties of the 2D Fourier Transform

2009-07-08T18:58:00.000-07:00

Activity 6.A Familiarization with FT of different 2D patterns

Below are the 2D patterns used in this activity:

SquareCircular annulussquare annulustwo slitstwo dots

Below are the respective FT of the 2D patterns shown above:FT of Square
The FT of a square aperture is a Bessel function intensity symmetric with x and y axes.FT of circular annulus
The FT of a circular annulus is a Bessel function intensity that is radially symmetric with origin.
FT of square annulus
This is similar with that of the square pattern only that there are other far-side intensities and smaller amplitude intensity at the origin.FT of two slits along y direction
The FT is at the opposite direction of the aperture. That is in this example, the slit is along the y direction but the FT plot is along the x-direction. Also, this FT is similar to the square annulus aperture missing the slit along the x-direction.FT of two dots
The FT of the two dots along the x-axis is just like a french fries that is oriented along the y-axis. Also, it can be observed that the vertical lines are not purely black and white stripes. It looks like the edges are grayscaling. I suspect that these are somewhat sinusoids in y-axis direction.

Sinusoidal waves pattern

sine patternFT of sine pattern f = 4
The FT of the sine pattern oriented horizontally with freq = 4 is represented by two dots along the y axis.FT of sine pattern f = 6FT of sine pattern f = 8
As the frequency of the pattern is increased, the separation between the two dots also increased. We can say that the FT is symmetric with the x-axis since our sine pattern is oriented along the x direction. Also, we can say that both x and y axes are frequency units.sine pattern with bias y + 1FT of sine pattern with bias y + 1
Observing the FT of the sine pattern above with see that the original sine wave has frequency = 8 which is represented by the previous result. There is an additional dot at the center of the plot representing 0 frequency. Thus we say that the FT has a property of linear combination since the FT plot shows the FT of sine with f = 8 and FT of the bias with obviously with f = 0 since it does not repeat.rotated sine pattern theta = 30FT of rotated sine pattern theta = 30
When the vertically oriented sine pattern is tilted by theta = 30, the spot concentration in the FT plot also tilts with respect to the x axis and the spot spreads along the y direction.sine product f = 4
FT of the above pattern
We can see that the FT follows the superposition since the FT displays that of two sine patterns oriented in opposite direction with same f = 4.
sine product pattern f = 4, f = 6
FT of the above pattern
This is just the same as the previous pattern the difference is just there are two distinct frequencies comprising the mat patternsine product plus two rotated sines theta = 30 & 60 f = 4FT of the above pattern
The above example is a combination of added sine patterns of the same freq along x and y direction and two rotated sines with theta = 30 and 60.

I give myself 10 points for finishing this activity. Its fun!

Activity 5 - Fourier Transform Model of Image Enhancement

2009-07-06T20:17:00.000-07:00

Image Enhancement through Fourier Transform

Below is an image of a circle of size 128 x 128 px. It is created using scilab and saved as an image to process.
The Fourier transform of a circle aperture or a lens is a bessel function that is symmetric at the origin. Although in the image it is quite small but if you zoom in the image you see that is just not a spot.Taking twice fft2 of the grayscale image before absolute value results in obtaining a circular image again. However, when this procedure is applied to a letter A aperture, it can be noticed that there is flipping from the original aperture.
In the next activity, convolution of two functions is investigated. A circular aperture image is to be convolved with a VIP pattern (as shown below) of the same size. That is, the fftshifted circular image pattern is element wise multiplied to the fft2 of the VIP pattern image. After getting the inverse fft, the resulting pattern is a flipped image along the diagonal and is distorted a little by the circular aperture.
Decreasing the size of the circular pattern results to a more blurring of the VIP pattern. Physically, as you see through smaller holes your vision is impaired. In the next activity, we investigated on how correlation works for in pattern recognition. First we have an image phrase shown below. The objective is to allow a search pattern of the letter A on the phrase image. This is done using correlation. The Fourier transform of A pattern is taken and it is multiplied element wise to the conjugate of the fft of phrase image.

The result is that the letters A in the phrase is highlighted. Thus, the highlight means that there is high correlation between the pattern A and to some of the attributes of phrase image.

For the final activity, we do template matching using correlation. Different templates are used like a vertical template = [-1 2 -1; -1 2 -1; -1 2 -1], horizontal and diagonal. The respective templates used are vertical, horizontal, left diagonal, right diagonal and spot for the following processed images. It can be seen that the kind of template used highlights the direction that is in the same direction at.

I give myself 10 points for finishing the activity. I would like to acknowledge Mam Jing and Earl for helping me find loopholes in the code. :)

Activity 3 - Image types

2009-06-22T18:46:00.000-07:00

In this activity, different types of images are being discussed namely binary, grayscale, indexed and true color images. Their properties are being discussed and sample types of images are presented.

Below is a binary image.

FileName: 6_m_airpistol.bmp
FileSize: 237882
Format: BMP
Width: 1481
Height: 1265
Depth: 8
StorageType: indexed
NumberOfColors: 2
ResolutionUnit: centimeter
XResolution: 62.990000
YResolution: 62.990000

http://www.louiscandell.com/pdf/targets/6_m_airpistol.BMP
Below is a grayscale image.

FileName: gray parrot.png
FileSize: 29342
Format: PNG
Width: 150
Height: 200
Depth: 8
StorageType: indexed
NumberOfColors: 256
ResolutionUnit: centimeter
XResolution: 72.000000
YResolution: 72.000000

http://en.wikipedia.org/wiki/File:Grayscale_8bits_palette_sample_image.png

Below is a true color image.

FileName: pahiyas.jpg
FileSize: 27871
Format: JPEG
Width: 400
Height: 300
Depth: 8
StorageType: truecolor
NumberOfColors: 0
ResolutionUnit: inch
XResolution: 160.000000
YResolution: 160.000000

Below is an indexed image

FileName: indexed.gif
FileSize: 16029
Format: GIF
Width: 155
Height: 148
Depth: 8
StorageType: indexed
NumberOfColors: 256
ResolutionUnit: centimeter
XResolution: 72.000000
YResolution: 72.000000

http://www.devx.com/projectcool/article/19997/0/page/8

CODE:

imfinfo('6_m_airpistol.bmp', 'verbose');
imfinfo('pahiyas.jpg', 'verbose');
imfinfo('gray parrot.png', 'verbose')
imfinfo('indexed.gif', 'verbose');

bw = imread('circle.bmp');
gr = imread('gray parrot.png');
in = imread('indexed.gif');
tc = imread('pahiyas.jpg');

Binary images are easy to find in the net since they are images with pixel values of 1 or 0. They are commonly known as black and white images. Also when presented in scilab, they are those with only 2 number of colors. Grayscale images are those images with 256 number of colors and of one layer only unlike truecolor images where they have 3 color channels for red, green and blue. That's why truecolor images are referred to as RGB images. I am not so sure with regards to indexed images since they are somewhat similar with grayscale images although they have additional stored information -- colormap value. Pixel values are based (ranges) only to some colormap values.

I give myself 9 points for taking all images of different types and also obtaining important information although i am not so sure with the indexed image sample.

References:
http://docs.gimp.org/2.2/en/gimp-images-in.html
http://www.mathworks.com/access/helpdesk/help/toolbox/images/index.html?
A3 – Image types and basic image enhancement -- manual

ROI extension

FileName: circles.bmp

FileSize: 179874

Format: BMP

Width: 324

Height: 185

Depth: 8

StorageType: truecolor

NumberOfColors: 0

ResolutionUnit: centimeter

XResolution: 37.800000

YResolution: 37.800000

FileName: circles_gray.bmp

FileSize: 61018

Format: BMP

Width: 324

Height: 185

Depth: 8

StorageType: indexed

NumberOfColors: 256

ResolutionUnit: centimeter

XResolution: 37.800000

YResolution: 37.800000

//histogram

x = imread('circles_gray.bmp');

mx = max(x);

sz = size(x,1)*size(x,2);

h = [];

for i = 1:mx

f = find(x == i);

h(i) = length(f);

end

scf(); plot2d(h/sz) // histogram plot

Looking at the histogram, it can be noticed that each colored circle has a distinct gray value (ROI’s) and the white background can be easily separated from the ROI’s through proper thresholding. We then convert the image to binary image and select first the ROI with the lowest gray value up to the highest gray value ROI. The following code and set of images demonstrates how I did it. Separating the objects of interest can be done using proper histogram thresholding.

//blue

xb = im2bw(x,0.12);

scf(); imshow(abs(xb-1), 2);

//blue and red

xbr = im2bw(x, 0.22);

scf(); imshow(abs(xbr-1), 2);

//red

scf(); imshow(xr-xbr, 2);

//blue, red, violet

xbrv = im2bw(x, 0.3);

scf(); imshow(abs(xbrv-1), 2);

//violet

scf(); imshow(xbr-xbrv, 2);

//blue, red, violet, green

xbrvg = im2bw(x, 0.8);

scf(); imshow(abs(xbrvg-1), 2);

//green

scf(); imshow(xbrv-xbrvg, 2);

//blue, red, violet, green, yellow

xbrvgy = im2bw(x, 0.94);

scf(); imshow(abs(xbrvgy-1), 2);

//yellow

scf(); imshow(xbrvg-xbrvgy, 2);

Since the sum of the individual color’s analytic area matches with all colors, then it must be correct. Also, the obtained total individual color’s area using Green’s is comparable with all color’s area. Thus, the validity of the code is exemplified with less than 1% error.

Activity 2 - Area Estimation for images with defined edges edges

2009-06-22T18:39:00.001-07:00

Using the Green's Theorem, areas of some polygons are obtained in square pixels unit. This analytically estimates the area of the polygon by summing all areas of smaller triangles constituting the polygon.

Images used in the activity:

Area Comparison:

For the analytical way of obtaining the area of the white region, I just summed up all the pixels since this is a binary image pixel will only have values of 1 or 0 with 1 representing the white area and 0 as black. Thus I can obtain the area analytically and compare the results with that of Green's theorem. The results are posted below:
Code:
[pol,map] = imread('I.bmp');
polbw = im2bw(pol,map,0.1);

[x, y] = follow(polbw);
plot2d(x,y)
sz = size(x,1);
xi = x(1:sz -1);
yi = y(1:sz -1);
xi1 = x(2:sz);
yi1 = y(2:sz);

area = 0.5 * sum(xi.*yi1 - yi.*xi1) + (sz/2)
areat = sum(polbw)
err = 100 - ((areat - area) / areat) * 100

Above is the code used to implement the Green's Theorem. It is then vectorized and limited the use of for loops to improve the computing time of the code. It is just a simple code to solve Green's.

The obtained accuracy is very significant since it is at 98 - 100 % accurate in estimating the area. Also, the algorithm has a correction factor added to the computed area of (sz/2) because it covers the contour plot from the follow function. This contour is still part of the original area and a small correction has to be added.

The code is limited to regions only with one filled area. The code does not apply to polygons like a "doughnut" since the follow function of scilab considers only the first connected contour like. In this case, the inner contour for this region will be neglected and a wrong calculated area will be obtained.

I give myself 10 points for satisfying the objectives of the activity.

Extension:

Estimating UP-Diliman's total land cover

I obtained a UP Diliman map from Wikimapia. To make it easier, I mouse-over the UP area and its boundary higlighted in yellow color. I then have the page print-screened. Using Gimp, I manipulated the image so that only the UP area will be left. I used the select by color tool (subtract), added an alpha channel and copied the selected part. Visit http://www.gimp.org/tutorials/Changing_Background_Color_2 for detailed steps. Since I am not skilled in gimp, i used morphological operations in Matlab to enhance more the image to fill in the holes. Specifically I used imclose function to "fill" the map. I finally have a filled-in UP Diliman map and have to invert the color so that white is inside the contour. The following images will illustrate more of the steps.

Then, the final image can now be processed for area estimation in square pixel units. To convert this in square meters, a conversion factor must be obtained. From the print-screened image, I determined how many pixels there are in 500 m after loading it in Paint. This can be seen at the lower left part of the image (scale).

I obtained a conversion ratio of 106 px:500 m and computed UP Diliman area summary:
area = 215179.5 sq px
areat = 215469. sq px
err = 0.1343581 %
up_area_m2 = 4787724.7 sq m
up_area_h = 478.77247 hectares
accuracy = 97.114092 %.
I can say that the algorithm is very good in the area estimation comparing with the reported 493 hectares UP Diliman land cover.

Activity 1 - Digital Scanning

2009-06-17T19:11:00.000-07:00

Objective: To reproduce a scanned plot using ratio and proportion of pixel values and its corresponding physical values

Softwares: OpenOffice.org Calc, Paint and Nero Photoshop

Procedure:

We looked for a one-plot graph in the library, photocopied it and scanned it.
We opened the image using paint and check if the axes of the graph are in a perpendicular position with respect to the image itself. If not, we used Nero Photoshop to rotate the graph in such a way that the axes are in a right angle position.
After that, we obtained the pixel location of the origin. From my image, it is at (2465, 1637) with the corresponding physical value of (0, 15).
Next is to obtain how many pixels correspond to a physical value both in the x and y tick marks. To do this, I obtained several data and get the average of the pixel values both in x and y axes. On the average, I obtained a 91.7 pixel values for 60 x-physical values and 69.83 pixel values for 5 y-physical values.
Then, we collected pixel locations of points along the graph and tabulated them in Spreadsheet. Pixel locations are seen at the lower right part of the Paint window. There are at least 20 data points collected.
We compute the x and y-factor multipliers using (physical value/pixel value) or (60/91.17) for x and (5/69.83) for y. This will be used in the conversion process.
Conversion of the obtained pixel locations involves shifting and flipping of graph. For the x-axis, we use the equation:

x_new = ((x-xorigin) * xfactor) + xphysorigin

while for the y-axis, we follow:

y_new = (((-1 * y) + yorigin) * yfactor) + yphysorigin

where

x and y - pixel points along the graph
xorigin and yorigin - x and y pixel location of origin (2465, 1637)
xfactor and y factor - pre-computed from above ((60/91.17), (5/69.83))
xphsyorigin and yphysorigin - x and y physical values of origin (0,15).

From there, we can now reproduce the scanned plot by plotting the x_new and y_new. We want also to check if we got a good new plot so we compare it with the original. We overlay the original graph on the plot chart. We do the following steps:

Crop only the chart area of the original graph excluding tick mark labels and titles. Only the x and y axes box should be included. Save it as another image. This will be used to overlay in the plot area.
Insert any image in the Spreadsheet, right click it and select Area-Bitmaps-Import. Import the image that should be overlayed on the graph. Save the current bitmap file and overwrite the existing file. Close the spreadsheet application and try filling the chart (reconstructed plot) with the saved bitmap image. It should now then appear on the list of bitmaps.
Uncheck the tile button and check autofit. That should do the overlaying of image on the graph.

The graph above is the superimposed plots of the reconstructed (green) and the original plot (black). It can be seen that the reconstructed plot is very similar with the scanned graph. I give myself 10 points.

url name 'plagiarized'

2009-06-15T19:46:00.000-07:00

ang url kong i-am-not-jonnocon.blogger.com ay hango sa i-am-yonin ni yonin dionisio na kaibigan ko sa yakal!

mabuhay blogger na ako!