Liknayang Inaplay 186

Activity 19 - Restoration of Blurred Image

2009-09-24T08:49:00.070+08:00

This last activity is about demonstration of restoration of a corrupted grayscale image with a known degradation function which in this case is motion blur and additive noise using Weiner filtering.

The following schematic diagram is a model of the image degradation and restoration process.

From the diagram, f_hat(x, y) is the approximate original image after the degraded image passes through one or more restoration filters, then the degraded image in the spatial domain can then be written in an equation as follows.

Here, h(x, y) represents the degradation function, f(x, y) is the original image and n(x, y) embodies the added noise. Meanwhile, its frequency representation in Fourier space is shown below.

The transfer function of image degradation implemented here is the following.

The variable T corresponds to the duration of exposure while a and b are the blurring displacements in the horizontal or x and vertical or y direction respectively.
In restoring the original image, the estimate image f_hat and the corrupted image f must have minimum mean square error between them. This condition is satisfied through the expression as follows called the minimum mean square error filter or the Weiner filter.

Note that the power spectrum of the noise and the original image are included in the Weiner filter which are shown below.

The power spectrum is similar to the modulus of the degradation transfer function as in the following equation.

As usual the first factor in the right hand side of the expression above denotes the complex conjugate of H(u, v).
In most cases, the power spectrum of the noise and undergraded image are not known, hence the Weiner filter can be approximated in the expression as follows.

Instead of the power spectra, it is replaced by K which is a constant that can be varied.

At first, the following grayscale image is obtained from the internet.

Source: http://tvtropes.org/pmwiki/pmwiki.php/Main/Ptitlefjcz80qe

A Gaussian noise with mean and standard deviation of 0.02 is then applied to the image. Next, the degradation transfer function parameters a, b, and T are varied resulting to the set of images as shown below (Arranged from top to bottom, left to right).

For T = 1: (1) a = 0.1, b = 0.1; (2) a = 0.01, b = 0.01; (3) a = 0.001, b = 0.001

For a = 0.01, b = 0.01: (1) T = 0.01; (2) T = 0.1; (3) T = 10; (4) T = 100

For constant exposure time, it can be observed that the degraded image becomes more discernible as the blurring displacements a and b decrease. Meanwhile, for constant blurring displacements, the added Gaussian noise is masked by the blurred image as the exposure time T increases.
Shown as follows are the restored images using Weiner filtering with the spectrum of the noise and undergraded image for different degradation transfer function parameters (Arranged from top to bottom, left to right).

For T = 1: (1) a = 0.1, b = 0.1; (2) a = 0.01, b = 0.01; (3) a = 0.001, b = 0.001

For a = 0.01, b = 0.01: (1) T = 0.01; (2) T = 0.1; (3) T = 10; (4) T = 100

Apparently, the restoration approaches the original image as the blurring displacements decrease as well as the exposure time increases.
The following last set of images is composed of restored images using Weiner filter with varying K values (Arranged from top to bottom, left to right).

For a = 0.01, b = 0.01, T = 1: (1) K = 0; (2) K = 0.000001; (3) K = 0.00001; (4) K = 0.0001; (5) K = 0.001; (6) K = 0.01; (7) K = 0.1; (8) K = 1

It can be noticed that since a constant, which is based from the ratio of the power spectrum of the noise to the undergraded image, is added in all the elements of the matrix in the process Weiner filtering, the added Gaussian noise is eliminated at the same time magnifying the blurring effect as K increases.
Weiner filtering with the power spectrum of the noise and the undergraded image is then recommended to be used when the power spectrum of the added noise and the undergraded image is known and also with minimal degradation while Weiner filtering with the constant K is more applicable for various corrupted image even without the power spectra as long as K is carefully chosen.

Since I am able to demonstrate a corrupted image with motion blur and additive noise, and restore it using Weiner filtering successfully, I grade myself 10/10.

I have worked individually in this activity however I have shared my insights to my classmates.

Appendix
The Scilab code below is utilized in this activity.

stacksize(4e7);

image = gray_imread('dissidia.jpg');
//scf(0);
//imshow(image, []);
//imwrite(normal(image), 'bwdissidia.bmp');

s = size(image);

noise_gauss = grand(s(1), s(2), 'nor', 0.02, 0.02);
noisy_image = image + noise_gauss;

N = fft2(noise_gauss);
F = fft2(image);

a = 0.01;
b = 0.01;
T = 1;

H = [];
for i = 1:s(1)
for j = 1:s(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.*F + N;
noisy_blurred = abs(ifft(G));

//scf(1);
//imshow(noisy_blurred, []);
//imwrite(normal(noisy_blurred), 'noisy blurred_a001_b001_T1.bmp');

SN = N.*conj(N);
SF = F.*conj(F);
restore1 = (((1)./H).*((H.*conj(H))./((H.*conj(H)) + (SN./SF)))).*G;
restore1 = abs(ifft(restore1));

//scf(2);
//imshow(restore1, []);
//imwrite(normal(restore1), 'restore_a001_b001_T1.bmp');

K = 0.01;
restore2 = (((1)./H).*((H.*conj(H))./((H.*conj(H)) + K))).*G;
restore2 = abs(ifft(restore2));

//scf(3);
//imshow(restore2, []);
//imwrite(normal(restore2), 'restoreK001.bmp');

Activity 18 - Noise Models and Basic Image Restoration

2009-09-15T09:06:00.150+08:00

In this activity, various noise models are applied in grayscale images, then different basic image restoration techniques are tested in recovering the original images.

Like in a grayscale image, noise is characterized by a probability distribution function (PDF). However, its PDF consists of random variables with parameters such as mean and standard deviation that may be altered.
The following are the most common noise models encountered in image processing. First below is the PDF of Gaussian noise which obviously follows a Gaussian curve with mu as its mean and sigma as its standard deviation.

The random variable z represents the grayscale value in the noise.
The PDF of Rayleigh noise is given as follows.

Its mean and variance are the following.

Next below is the PDF of Erlang or Gamma noise.

Whereas its mean and variance are shown as follows.

PDF of exponential noise is described by the following equations.

Its mean and variance are shown below.

Uniform noise has the following PDF.

Then its mean and variance are as follows.

Lastly, impulse or better known as salt-and-pepper noise is considered here. This has PDF expressed below.

Meanwhile, there are four basic image restoration techniques implemented in this activity. These are spatial filters which reduce added noise in the image.
Arithmetic mean filter operates through the following equation.

S_xy is a set of coordinates within a chosen rectangular sub-image window with size m x n in the noisy image centered at (x, y). Hence, it filters the noise by averaging the grayscale value at (x, y).
A noisy image under geometric mean filtering works with the equation as follows.

The harmonic mean filter which is shown below works well for salt noise or noisy pixels with grayscale values close to white but not with pepper noise or noisy pixels with grayscale values near to black.

Lastly, contraharmonic mean filter has expression given by the following.

The variable Q is called the order of the filter which may be positive that is good for pepper noise or negative that eliminates salt noise. Hence, it cannot remove both noise simultaneously.

Initially, a grayscale image as follows is created from Paint which only consists of three grayscale values.

Its PDF with grayscale values ranging from 0 to 255 is plotted in the following figure.

Apparently, the three peaks correspond to the three grayscale values of the image. Then, using the built-in function grand of Scilab, majority of the noise models are generated with the same size as the image (150 x 150). The image degradation is then simulated by simply adding the generated noise to the test image. Note that the parameters such as mean and standard deviation of the noise models in grand are carefully set to prevent the merging of the three grayscale levels in the original or test image and so the noisy image is still recognizable.
After adding Gaussian noise with mean 0.5 and standard deviation of 0.075 to the image, the image and its PDF becomes the following.

As can be observed, the Gaussian profiles have maximum near the three grayscale values implying that the noisy image is still discernable. The next set of figures are the filtered images using the arithmetic, geometric, harmonic, and contraharmonic mean filter arranged from left to right, top to bottom.

The sub-image window used for all the filters has size 3 x 3. However, the algorithm does not apply at the corners and edges of the noisy image, thus necessary boundary conditions are considered at these locations. For the contraharmonic mean filter, Q is equal to 1. It can be observed that the arithmetic and the contraharmonic mean filters are better compared in the geometric and harmonic mean filter for the added Gaussian noise.
Next, the Rayleigh noise is generated by first downloading a module called modnum. Rayleigh noise generator can then be used when the loader file of the module is loaded. The function genrayl in the module is then used with noise generator parameter set to 0.1. The resulting noisy image and its PDF after adding the Rayleigh noise are shown below.

The filtered images are as follows.

In this type of noise, all filters work well in restoring the test image.
Gamma noise is then added to the test image again using grand. Its parameters shape and scale are set to 1 and 10 respectively. The following are the resulting noisy image and its PDF.

The figures with Gamma noise after using the four filters are shown below.

The contraharmonic mean filter somehow fails in restoring the test image for the Gamma noise while the harmonic mean filter is the best filter for this case.
The test image with added exponential noise and its PDF are as follows. The mean of the exponential noise is set to 0.1.

Notice that the PDF is similar when Gamma noise is applied. The following set of figures is the result after using the four filters in the test image with Gamma noise.

And so the filtered images are similar to the added Gamma noise.
Again from grand, another noise model which is uniform noise is generated here where a uniform distribution is set from 0 to 0.2. Shown below are the noisy image and its PDF with this uniform noise.

The PDF indeed illustrate a uniform random grayscale values within the range of the three graylevels. Then, the restored noisy images with uniform noise through the four filters are as follows.

Noting that uniform noise is present, the four filters are effective in enhancing the noisy image.
The impulse or salt-and-pepper noise is simulated through the function imnoise where the noise density or probability of making salt or pepper noise is set to 0.01. The resulting noisy image and its PDF are the following.

The PDF is similar to the PDF of the test image without any noise because impulse noise only adds white (salt) or black (pepper) pixels while the test image is composed of only white, black, and gray pixels. Thus, what is only different in this PDF is the frequency of the three graylevels. Illustrated below are the filtered noisy image with impulse noise through the four filters.

It can be observed that the arithmetic mean filter is best in the test image with impulse noise. It is proven here that the harmonic mean filter eliminates salt noise but fails in pepper noise. Meanwhile, the contraharmonic mean filter with Q = 1 removes pepper noise but some salt noise still remain. Different Q values of the contraharmonic mean filter are also tested here. The filtered noisy image with Q equals to -1, 2, -2, and 0 arranged from left to right, top to bottom are as follows.

From the results above, positive Q values in the contraharmonic mean filter indeed eliminate pepper noise and not salt noise while negative Q values are truly good to salt noise but not to pepper noise. However, with Q equal to 0, both salt and pepper noise are minimized.

The whole process is also tested for a grayscale image which is not restricted to only three discrete grayscale values. The grayscale image of a beach obtained from the internet is shown below.

Source: http://www.advimglib.com/samples.html

The PDF of the image above is shown as follows.

Apparently, the PDF of this test image does not contain only three graylevels. The test image with Gaussian noise and its PDF are the following.

Below are the filtered test image with Gaussian noise using the four filters.

The test image with Rayleigh noise and its PDF are as follows.

The results of the filtered test image with Rayleigh noise through the four filters are the following.

Shown below are the test image with Gamma noise and its PDF.

The filtered test image with Gamma noise using the four filters are as follows.

Next are the test image with exponential noise and its PDF.

The following set of figures is the filtered test image with exponential noise.

Then, the test image with uniform noise and its PDF are shown below.

The filtered test image with uniform noise after applying the four filters are as follows.

Lastly, the following are the test image with impulse noise and its PDF.

Below are the filtered test image with impulse noise through the four filters.

The Q values of the contraharmonic mean filter are also varied to -1, 2, -2, and 0 resulting to the figures as follows.

Since the test image now has distribution of grayscale values, the four filters resulted to almost the same filtered image except in the test image with impulse noise unlike in the test image with only three grayscale values where a certain filter is preferred over the others. Overall, it can be deduced that the arithmetic mean filter is consistent in restoring both test image with three discrete graylevels and with distributed grayscale values to every type of noise added.

I grade myself 10/10 because I am able to simulate all the noise models introduced in this activity and I have successfully implemented the four basic image restoration techniques in the test images with different the noise models applied to them.

I thank Neil in providing me the test image with only three grayscale values.

Appendix
Below is the whole Scilab code utilized in this activity.

stacksize(4e7);

image = gray_imread('beach.jpg');
s = size(image);

//scf(0);
//imshow(image, []);

a = 1;
grayscale = [];
pixels = [];
for i = 0:255
[x, y] = find(round(255*normal(image)) == i);
grayscale(a) = i;
pixels(a) = length(x);
a = a + 1;
end

//scf(1);
//plot(grayscale, normal(pixels));
//title('PDF of Image without Noise');
//xlabel('Grayscale Values');

noise_gauss = grand(s(1), s(2), 'nor', 0.5, 0.075);
noise_ray = matrix(genrayl(0.1, s(1)*s(2)), s(1), s(2));
noise_gamma = grand(s(1), s(2), 'gam', 1, 10);
noise_exp = grand(s(1), s(2), 'exp', 0.1);
noise_uni = grand(s(1), s(2), 'unf', 0, 0.2);
noisy_image_salt = imnoise(image, 'salt & pepper', 0.01);

//noisy_image = round(255*normal(noisy_image_salt));
noisy_image = round(255*normal(image + noise_gauss));

//scf(2);
//imshow(noisy_image, []);
//imwrite(normal(noisy_image), 'image gauss noise.bmp');

a = 1;
grayscale = [];
pixels = [];
for i = 0:255
[x, y] = find(round(255*normal(noisy_image)) == i);
grayscale(a) = i;
pixels(a) = length(x);
a = a + 1;
end

//scf(3);
//plot(grayscale, normal(pixels));
//title('PDF of Image with Gaussian Noise');
//xlabel('Grayscale Values');

// Filtering Noisy Image

// Arithmetic Mean Filter

// Size of Rectangle Subimage Window (3 x 3)
m = 3;
n = 3;

for a = 1:s(1)
for b = 1:s(2)
if a == 1 & b == 1
g1 = [noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered1(a, b) = (1/(2*2))*sum(g1);
elseif a == 1 & b == s(2)
g1 = [noisy_image(a, b - 1) noisy_image(a, b); noisy_image(a + 1, b - 1) noisy_image(a + 1, b)];
filtered1(a, b) = (1/(2*2))*sum(g1);
elseif a == s(1) & b == 1
g1 = [noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b) noisy_image(a, b + 1)];
filtered1(a, b) = (1/(2*2))*sum(g1);
elseif a == s(1) & b == s(2)
g1 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b); noisy_image(a, b - 1) noisy_image(a, b)];
filtered1(a, b) = (1/(2*2))*sum(g1);
elseif a == 1 & b ~= 1 & b ~= s(2)
g1 = [noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b - 1) noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered1(a, b) = (1/(2*3))*sum(g1);
elseif b == 1 & a ~= 1 & a ~= s(1)
g1 = [noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered1(a, b) = (1/(3*2))*sum(g1);
elseif b == s(2) & a ~= 1 & a ~= s(1)
g1 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b); noisy_image(a, b - 1) noisy_image(a, b); noisy_image(a + 1, b - 1) noisy_image(a + 1, b)];
filtered1(a, b) = (1/(3*2))*sum(g1);
elseif a == s(1) & b ~= 1 & b ~= s(2)
g1 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1)];
filtered1(a, b) = (1/(2*3))*sum(g1);
else
g1 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b - 1) noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered1(a, b) = (1/(m*n))*sum(g1);
end
end
end

//scf(4);
//imshow(filtered1, []);
//imwrite(normal(filtered1), 'gauss filtered1.bmp');

// Geometric Mean Filter

for a = 1:s(1)
for b = 1:s(2)
if a == 1 & b == 1
g2 = [noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered2(a, b) = prod(g2)^(1/(2*2));
elseif a == 1 & b == s(2)
g2 = [noisy_image(a, b - 1) noisy_image(a, b); noisy_image(a + 1, b - 1) noisy_image(a + 1, b)];
filtered2(a, b) = prod(g2)^(1/(2*2));
elseif a == s(1) & b == 1
g2 = [noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b) noisy_image(a, b + 1)];
filtered2(a, b) = prod(g2)^(1/(2*2));
elseif a == s(1) & b == s(2)
g2 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b); noisy_image(a, b - 1) noisy_image(a, b)];
filtered2(a, b) = prod(g2)^(1/(2*2));
elseif a == 1 & b ~= 1 & b ~= s(2)
g2 = [noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b - 1) noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered2(a, b) = prod(g2)^(1/(2*3));
elseif b == 1 & a ~= 1 & a ~= s(1)
g2 = [noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered2(a, b) = prod(g2)^(1/(3*2));
elseif b == s(2) & a ~= 1 & a ~= s(1)
g2 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b); noisy_image(a, b - 1) noisy_image(a, b); noisy_image(a + 1, b - 1) noisy_image(a + 1, b)];
filtered2(a, b) = prod(g2)^(1/(3*2));
elseif a == s(1) & b ~= 1 & b ~= s(2)
g2 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1)];
filtered2(a, b) = prod(g2)^(1/(2*3));
else
g2 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b - 1) noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered2(a, b) = prod(g2)^(1/(m*n));
end
end
end

//scf(5);
//imshow(filtered2, []);
//imwrite(normal(filtered2), gauss filtered2.bmp');

// Harmonic Mean Filter

small_number = 0.000001;
for a = 1:s(1)
for b = 1:s(2)
if a == 1 & b == 1
g3 = [noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered3(a, b) = (2*2)/(sum((1)./(g3 + small_number)));
elseif a == 1 & b == s(2)
g3 = [noisy_image(a, b - 1) noisy_image(a, b); noisy_image(a + 1, b - 1) noisy_image(a + 1, b)];
filtered3(a, b) = (2*2)/(sum((1)./(g3 + small_number)));
elseif a == s(1) & b == 1
g3 = [noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b) noisy_image(a, b + 1)];
filtered3(a, b) = (2*2)/(sum((1)./(g3 + small_number)));
elseif a == s(1) & b == s(2)
g3 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b); noisy_image(a, b - 1) noisy_image(a, b)];
filtered3(a, b) = (2*2)/(sum((1)./(g3 + small_number)));
elseif a == 1 & b ~= 1 & b ~= s(2)
g3 = [noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b - 1) noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered3(a, b) = (2*3)/(sum((1)./(g3 + small_number)));
elseif b == 1 & a ~= 1 & a ~= s(1)
g3 = [noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered3(a, b) = (3*2)/(sum((1)./(g3 + small_number)));
elseif b == s(2) & a ~= 1 & a ~= s(1)
g3 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b); noisy_image(a, b - 1) noisy_image(a, b); noisy_image(a + 1, b - 1) noisy_image(a + 1, b)];
filtered3(a, b) = (3*2)/(sum((1)./(g3 + small_number)));
elseif a == s(1) & b ~= 1 & b ~= s(2)
g3 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1)];
filtered3(a, b) = (2*3)/(sum((1)./(g3 + small_number)));
else
g3 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b - 1) noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered3(a, b) = (m*n)/(sum((1)./(g3 + small_number)));
end
end
end

//scf(6);
//imshow(filtered3, []);
//imwrite(normal(filtered3), 'gauss filtered3.bmp');

// Contraharmonic Mean Filter

Q = 0;
for a = 1:s(1)
for b = 1:s(2)
if a == 1 & b == 1
g4 = [noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
elseif a == 1 & b == s(2)
g4 = [noisy_image(a, b - 1) noisy_image(a, b); noisy_image(a + 1, b - 1) noisy_image(a + 1, b)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
elseif a == s(1) & b == 1
g4 = [noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b) noisy_image(a, b + 1)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
elseif a == s(1) & b == s(2)
g4 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b); noisy_image(a, b - 1) noisy_image(a, b)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
elseif a == 1 & b ~= 1 & b ~= s(2)
g4 = [noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b - 1) noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
elseif b == 1 & a ~= 1 & a ~= s(1)
g4 = [noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
elseif b == s(2) & a ~= 1 & a ~= s(1)
g4 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b); noisy_image(a, b - 1) noisy_image(a, b); noisy_image(a + 1, b - 1) noisy_image(a + 1, b)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
elseif a == s(1) & b ~= 1 & b ~= s(2)
g4 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
else
g4 = [noisy_image(a - 1, b - 1) noisy_image(a - 1, b) noisy_image(a - 1, b + 1); noisy_image(a, b - 1) noisy_image(a, b) noisy_image(a, b + 1); noisy_image(a + 1, b - 1) noisy_image(a + 1, b) noisy_image(a + 1, b + 1)];
filtered4(a, b) = (sum((g4 + small_number).^(Q + 1)))/(sum((g4 + small_number).^Q));
end
end
end

//scf(7);
//imshow(filtered4, []);
//imwrite(normal(filtered4), 'gauss filtered4.bmp');

Activity 17 - Photometric Stereo

2009-09-08T09:03:00.041+08:00

The shape of a surface is estimated and then plotted 3-dimensionally in this activity using photometric stereo.

In photometric stereo, the 3-dimensional form of an object can be determined through captured intensity profiles of that object. The intensity profile is varied through different locations of an illuminating source. A point source is considered here. The location of N point sources is defined in a matrix as follows.

A row represents a point source and each row corresponds to its x, y, and z component. Now, for N number of images of the object, the intensity at (x, y) on the surface is given by the following.

This can be written in a matrix operation shown below.

Note that g is a function of the reflectance of the surface multiplied to the normal vector at (x, y) as shown in the following equation.

Since g is initially unknown, this can be obtained using the matrix equation as follows.

Then, the normal vector of the surface is calculated through the expression below.

Note that the denominator from the equation above is the length of g computed employing Pythagorean theorem in 3-dimension. From the normal vector, the shape of the object can be determined in the following definition of its surface elevation at (x, y).

The expression above is then rewritten as follows.

The gradient of c(x, y, z) shown below is next obtained in order to relate the estimated surface normals n_x, n_y, and n_z using photometric stereo and the partial derivative of f.

Meanwhile, the relation between the surface normals and the partial derivative of f is the following.

Finally, the object's surface elevation at point (u, v) is evaluated through the line integral as follows.

A Matlab file containing four images of a surface is loaded in Scilab using the function loadmatfile. These images are shown below.

Each image above is illuminated by a point source located respectively at
V_1 = [0.085832, 0.17365, 0.98106]
V_2 = [0.085832, -0.17365, 0.98106]
V_3 = [0.17365, 0, 0.98481]
V_4 = [0.16318, -0.34202, 0.92542]

These four point source locations are then placed in matrix V with size 4 x 3 as well as the four intensity profiles in matrix I with size 4 x (128*128) or 4 x 16384. Note that each loaded intensity profile has size 128 x 128. The size of matrix g and surface normal n obtained is thus 3 x 16384. To prevent division by zero, a very small value 0.000001 is added to every element of the length or magnitude of g and to all elements of n_z. The calculated partial derivatives of f with respect to x and y are recasted into a matrix with size 128 x 128 for evaluation of the surface elevation at (u, v). The line integrals are operated using cumsum of Scilab where the limit of integration from 0 to u means the cumulative summation along the column (from 1st column to 128th column) of the partial derivative of f with respect to x. Conversely, the limit of integration from 0 to v is the cumulative summation along the row (from 1st row to 128th row) of the partial derivative of f with respect to y.
The following figure is then the plotted surface of the object in 3-dimension.

Apparently, the shape of the surface is successfully reconstructed 3-dimensionally. From the 3-D plot, it is revealed that the object is hemispherical with four divided quadrants as observed from the intensity profiles.

Since I am able to show the shape of the given surface successfully using photometric stereo, I grade myself 10/10.

I worked with Gary in completing this activity.

Appendix
The whole Scilab code utilized in this activity is as follows.

loadmatfile('Photos.mat');

//scf(0);
//imshow(I1, []);
//scf(1);
//imshow(I2, []);
//scf(2);
//imshow(I3, []);
//scf(3);
//imshow(I4, []);

I1 = matrix(I1, 1, size(I1, 1)*size(I1, 2));
I2 = matrix(I2, 1, size(I2, 1)*size(I2, 2));
I3 = matrix(I3, 1, size(I3, 1)*size(I3, 2));
I4 = matrix(I4, 1, size(I4, 1)*size(I4, 2));

I = [I1; I2; I3; I4];

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

small_number = 0.000001;
g = inv(V'*V)*V'*I;
magnitude = sqrt((g(1,:).^2) + (g(2,:).^2) + (g(3,:).^2)) + small_number;

for i = 1:3
n(i,:) = g(i,:)./magnitude;
end

n_x = n(1,:);
n_y = n(2,:);
n_z = n(3,:) + small_number;
partial_f_to_x = -n_x./n_z;
partial_f_to_y = -n_y./n_z;
partial_f_to_x = matrix(partial_f_to_x, 128, 128);
partial_f_to_y = matrix(partial_f_to_y, 128, 128);
z = cumsum(partial_f_to_x, 2) + cumsum(partial_f_to_y, 1);

//scf(4);
//plot3d(1:128, 1:128, z);

Activity 16 - Neural Networks

2009-09-03T08:51:00.046+08:00

For the third time after Activity 14 and 15, grouping objects belonging to different classes is done in this acitvity. The artificial neural network (ANN) algorithm is implemented here to classify these sets of objects.

Again, the features of a sample which are the pixel area, normalized chromaticity values in red, green, and blue obtained from Activity 14 are used in this activity. This time, 3 classes from Activity 14 which are Snacku rice crackers, Superthin biscuits, and Nova multigrain snack with 10 samples each are accounted here unlike from Activity 15 which considers only 2 classes. These are shown below.

ANN recognizes patterns through learning by example. This is similar to how the human brain works in classifying objects. Once it has learned after some training cycles, ANN can quickly and accurately operate pattern processing.
The model of an artificial neuron which is the basic unit of ANN is illustrated as follows.

From the model above, the weighted inputs x_1, x_2, ..., x_i multiplied to their corresponding synaptic strengths or weights w_1, w_2, ..., w_i of a neuron are summed to produce a which then acts to the activation function g. Then, the output z is transmitted to other neurons as part of the ANN's learning period.
ANN has 3 layers namely the input layer, hidden layer, and the output layer. The input layer consists the features that are then transmitted to every neuron of the hidden layer noting that all neurons of the preceding layer are interconnected to every neurons of the succeeding layer. The results are then passed on to the output layer which determines the class to where the input features belong. The following is an illustration of a simple ANN.

Using the ANN_Toolbox_0.4.2 loaded in Scilab, built-in functions ann_FF_init, ann_FF_Std_online, and ann_FF_run are employed to simulate an ANN to classify 15 test samples. The first function ann_FF_init initializes the network which has l number of neurons in the input layer, m number of neurons in the hidden layer, and n number of neurons in the output layer. This then becomes the weight of ann_FF_Std_online training period. The features of the training set, their corresponding class matrix, the network, its learning rate, and the number of training cycles are other arguments of the second function. The argument learning rate is together with the threshold of the error tolerated by the network which is usually set to 0. The last built-in function ann_FF_run with weight as the result of the second function is finally applied to the features of the 15 test samples. This also has the argument of the network same in the first two functions. It returns the classes to where each sample belongs.
Note that the 30 samples are again divided into two where the first serves as the training set with 5 samples in each class for the ANN to learn and the second represents the test set to be classified, again with 5 samples in each class. The samples of the training set have class matrix [0, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1, 1, 1, 1] because it is observed that the built-in functions in the ANN toolbox of Scilab operate only between 0 and 1. However, the output classes are made to return rounded integer values of 0 for Snacku, 1 for Superthin, and 2 for Nova. Furthermore, it is ascertained that these functions cannot recognize features greater than 1. Since the pixel area feature is obviously greater than 1, it is normalized for the training and test samples. Meanwhile, the network used is composed of 4 neurons in the input layers where each neuron represents a feature, 8 neurons in the hidden layer which may be altered to multiples of 4, and 1 neuron in the output layer corresponding to the class assigned to a test sample.
The tables below are the summary of the test samples' classification where the training parameters are varied to observe the ANN's behavior (click the tables for larger view).

The first table above has 500 training cycles while the table below it has 1000 training cycles. Both are tested with learning rates of 1.0 and 5.0. The 15 test samples are classified correctly for all different training parameters used. Moreover, it can be observed that as the number of training cycles and the learning rate are increased, the outputs approach 0, 0.5 and 1 and thus approach 0 for Snacku, 1 for Superthin, and 2 for Nova.

The order of the test samples are also randomized to see if they are still to be grouped correctly. Here, 1000 training cycles and learning rate of 5.0 are employed. The results are summarized in the table as follows (click the table for larger view).

Apparently, all test samples are still grouped appropriately to their corresponding classes even these inputs are in different order. Thus, ANN's class recognition does not require inputs in arranged manner.

I grade myself 10/10 in this activity because I have successfully implemented ANN for class recognition and all test samples are grouped correctly.

Thanks to Miguel who sent a working version of the ANN toolbox for Scilab 4.1.2 for the whole class.

Appendix
The following Scilab code is utilized in this activity. This is based from the code of Mr. Cole Fabros in his blog (http://cole-ap186.blogspot.com) where it is mentioned that the code is originally made by Mr. Jeric Tugaff (Both are 2008 AP 186 students).

snacku_training = fscanfMat('snacku_training.txt');
superthin_training = fscanfMat('superthin_training.txt');
nova_training = fscanfMat('nova_training.txt');
snacku_test = fscanfMat('snacku_test.txt');
superthin_test = fscanfMat('superthin_test.txt');
nova_test = fscanfMat('nova_test.txt');

snacks_training = [snacku_training; superthin_training; nova_training];
snacks_training(:,1) = snacks_training(:,1)/max(snacks_training(:,1));
snacks_training = snacks_training';
snacks_test = [snacku_test; superthin_test; nova_test];
snacks_test(:,1) = snacks_test(:,1)/max(snacks_test(:,1));
snacks_test = snacks_test';

rand('seed', 0);

network = [4, 8, 1];
groupings = [0 0 0 0 0 0.5 0.5 0.5 0.5 0.5 1 1 1 1 1];
learning_rate = [5, 0];
training_cycle = 1000;

index = [1:15]';
random_index = grand(1, 'prm', index);

random_snacks = [];
for i = 1:length(index)
random_snacks(:,i) = snacks_test(:,random_index(i));
end

training_weight = ann_FF_init(network);
weight = ann_FF_Std_online(snacks_training, groupings, network, training_weight, learning_rate, training_cycle);
class = ann_FF_run(random_snacks, network, weight);

for j = 1:length(index)
if class(j) < 0.2
class(j) = class(j);
elseif 0.4 < class(j) & class(j) < 0.6
class(j) = 0.5 + class(j);
else
class(j) = 1.0 + class(j);
end
end

rclass = round(class);

Activity 15 - Probabilistic Classification

2009-08-27T09:16:00.024+08:00

This activity is similar to the previous one except that another method is implemented here. Objects are grouped here using probabilistic classification.

Instead of 3, only 2 classes are considered here which are all from Activity 14, Snacku rice crackers and Superthin biscuits. As assembled previously, 10 samples represent each class. These are shown below.

Also, the 4 features which are the sample's pixel area, normalized chromaticity values in red, green, and blue and their corresponding data from Activity 14 are used for probabilistic classification. Each class is again divided into a training set and a test set where each set contains 5 samples.
The probabilistic classification is rooted from the algorithm of linear discriminant analysis (LDA). From the following discriminant function, a test object k is assigned to class i that has maximum f_i.

Note that the transpose matrix operation is represented by T and -1 is the inverse matrix operation. The last term contains the prior probability of class i shown below.

From the equation above, n_i is the number of samples in class i and N is the total number of samples in the training set. Meanwhile, matrix C is called the pooled within class covariance matrix and this obtained from matrix c which is called the covariance matrix. These are as follows.

Note that there are g total number of classes. Covariance matrix is determined from the mean corrected data and the features of the training set which are represented by the following expressions from top to bottom respectively.

The global mean vector mu consists the average of each feature w_i,j,m coming from all classes of the training set and M is the total number of features. Going back to the discriminant function f_i, mu_i is composed of the mean of each feature just from class i in the training set. Thus, using LDA to classify an unknown or test sample, a training set is still required in order to establish the discriminant function which categorize the test sample with the same features possessed by the training set to a class.

The table below summarizes the data gathered from the training set (click the table for larger view).

Finally, the results of the LDA to the test set are summarized in the following table (click the table for larger view).

It can be observed that the values between the two discriminant functions are very small for all samples of the test set compared to the large deviations between the minimum distances acquired in Activity 14. The large values of discriminant function are because of the quantified feature pixel area which is relatively greater than the normalized chromaticity values. Nevertheless, all 10 samples of the test set are correctly grouped to their respective classes.

I grade myself 10/10 in this activity because I have implemented LDA properly and grouped all test samples to their expected classes.

I have successfully finished the activity individually. However, I have helped my classmates on how to handle this activity.

Appendix
The whole Scilab code below is utilized in this activity.

snacku_training = fscanfMat('snacku_training.txt');
superthin_training = fscanfMat('superthin_training.txt');
snacku_test = fscanfMat('snacku_test.txt');
superthin_test = fscanfMat('superthin_test.txt');

x_training = [snacku_training; superthin_training];
x_test = [snacku_test; superthin_test];

myu1 = [mean(snacku_training(:,1)), mean(snacku_training(:,2)), mean(snacku_training(:,3)), mean(snacku_training(:,4))];
myu2 = [mean(superthin_training(:,1)), mean(superthin_training(:,2)), mean(superthin_training(:,3)), mean(superthin_training(:,4))];

myu = [mean(x_training(:,1)), mean(x_training(:,2)), mean(x_training(:,3)), mean(x_training(:,4))];

x01 = [snacku_training(1,:) - myu; snacku_training(2,:) - myu; snacku_training(3,:) - myu; snacku_training(4,:) - myu; snacku_training(5,:) - myu];
x02 = [superthin_training(1,:) - myu; superthin_training(2,:) - myu; superthin_training(3,:) - myu; superthin_training(4,:) - myu; superthin_training(5,:) - myu];

n1 = 5;
n2 = 5;

c1 = (x01'*x01)/n1;
c2 = (x02'*x02)/n2;

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

invC = inv(C);

P1 = n1/(n1 + n2);
P2 = n2/(n1 + n2);
P = [P1; P2];

for k = 1:(n1 + n2)
f1(k) = myu1*invC*x_test(k,:)' - (1/2)*myu1*invC*myu1' + log(P(1));
f2(k) = myu2*invC*x_test(k,:)' - (1/2)*myu2*invC*myu2' + log(P(2));
end

f = [f1, f2];
//fprintfMat('f_test.txt', f);

class = [];
for i = 1:(n1 + n2)
[a, b] = find(f(i,:) == max(f(i,:)));
class = [class; b];
end

Activity 14 - Pattern Recognition

2009-08-25T09:20:00.031+08:00

The main objective of this activity is to classify various objects with different characteristics through pattern recognition.

In pattern recognition, the characteristics of an object belonging to a type or class are called features. These are quantified which may be the object's area, color, and many others. Pattern is then a set of features and a class is a set of patterns which share common features. The task is thus finding the features of an unknown object belonging to one of the several classes.

There are 3 classes assembled in this activity. They are all edible which include Snacku rice crackers, Superthin biscuits, and Nova multigrain snack. Each class is composed of 10 samples. The setup is shown below.

The 30 samples are grouped according to their respective classes for simplified image processing. Each class is divided into two, 5 samples serve as the training set and the other 5 correspond to the test set. Features that have been decided to use for pattern recognition are the sample's pixel area, average normalized chromaticity values in red, green, and blue for a total of 4.
The pixel area is attained by binarizing and inverting the cropped image of a class using the most appropriate predetermined threshold value. The funciton bwlabel of Scilab is then implemented to sum up the white pixels of each cluster equivalent to the pixel area of a sample and this is similar to the algorithm in Activity 9.
Meanwhile, parametric segmentation in Activity 12 is employed to get the normalized chromaticity values in red, green, and blue of each sample. This is operated through cropped individual sample images of each class. Note that the normalized chromaticity values in red, green, and blue are ensured that these just came from the samples (not with the background) by extracting the pixel coordinates of the segmented images.
Next, quantified features are ordered into a matrix called feature vector. Hence, in this case, a sample's feature vector consists 4 elements. Mean feature vectors representing each class are then obtained using the following equation.

Here, N_j is the total number number of samples in a class w_j which is here, equal to 5 while the 1 x 4 matrix x_j is composed of all 4 feature vectors acquired from image processing and then this is summed for N_j samples. Through minimum distance classification stated by the following expression, the features of each sample in the test set is then compared to the matrix of mean features m_j in the training set.

From the expressions above, x is the feature vector of a test sample, W is the total number of class and T is the transpose matrix operation. The class index j of the minimum distance D_j(x) is then assigned to the Nth sample from the test set of a class.

The results from the training set are shown in the table below (click the table for larger view).

Finally, the resulting classification of the test set is summarized in the following table (click the table for larger view).

It can be observed that the calculated distance D_j(x) is relatively large when a test sample does not match the expected class. Notice that the distances are big numbers because recall that the quantified feature pixel area is relatively larger than the normalized chromaticity values. But most importantly, all 15 samples in the test set are classified correctly.

Since I have understood the process of pattern recognition using the minimum distance classification and at the same time, successfully classifying all the samples in the test set, I grade myself 10/10 in this activity.

I have successfully completed this activity through discussions with Gary.

Appendix
The whole Scilab code below is utilized in this activity.

stacksize(4e7);

//snacku = gray_imread('snacku.jpg');

//scf(0);
//imshow(snacku);
//imwrite(snacku, 'gray snacku.bmp');

snack_A = [];
for a = 1:3
snack = gray_imread('gray snack'+string(a)+'.bmp');

if a == 1
snack = im2bw(snack, 175/255);
elseif a == 2
snack = im2bw(snack, 160/255);
else
snack = im2bw(snack, 155/255);
end

snack = 1 - snack;

//scf(a);
//imshow(snack);

[L, n] = bwlabel(snack);

A = [];
for i = 1:n
b = (L == i);
A(i) = sum(b);
snack_A(i, a) = A(i);
end

end

snacku_A = snack_A(:,1);
superthin_A = snack_A(:,2);
nova_A = snack_A(:,3);

training_snack_A = [snacku_A(1:5), superthin_A(1:5), nova_A(1:5)];
test_snack_A = [snacku_A(6:10), superthin_A(6:10), nova_A(6:10)];

// Average Normalized RGB Values

mean_snacku_r = [];
mean_snacku_g = [];
mean_snacku_b = [];
for j = 1:10
image = imread('snacku'+string(j)+'.jpg');
I = imread('snacku ROI.jpg');

R = I(:,:,1);
G = I(:,:,2);
B = I(:,:,3);

RGB = R + G + B;

r = R./RGB;
b = B./RGB;
g = G./RGB;

Ri = image(:,:,1);
Gi = image(:,:,2);
Bi = image(:,:,3);

RGBi = Ri + Gi + Bi;

ri = Ri./RGBi;
bi = Bi./RGBi;
gi = Gi./RGBi;

// Parametric Segmentation

rmean = mean(r);
gmean = mean(g);
rsigma = stdev(r);
gsigma = stdev(g);

rp = (1/(rsigma*sqrt(2*%pi)))*exp(-((ri - rmean).^2)/(2*rsigma^2));
gp = (1/(gsigma*sqrt(2*%pi)))*exp(-((gi - gmean).^2)/(2*gsigma^2));

rgp = round(rp.*gp);

//scf(j + 3);
//imshow(rgp);

[x, y] = find(rgp ~= 0);

snacku_r = [];
snacku_g = [];
snacku_b = [];
for k = 1:length(x)
snacku_r = [snacku_r, ri(x(k), y(k))];
snacku_g = [snacku_g, gi(x(k), y(k))];
snacku_b = [snacku_b, bi(x(k), y(k))];
end

mean_snacku_r(j,:) = mean(snacku_r);
mean_snacku_g(j,:) = mean(snacku_g);
mean_snacku_b(j,:) = mean(snacku_b);

end

training_mean_snacku_r = mean_snacku_r(1:5);
training_mean_snacku_g = mean_snacku_g(1:5);
training_mean_snacku_b = mean_snacku_b(1:5);
training_color_snacku = [mean(training_mean_snacku_r), mean(training_mean_snacku_g), mean(training_mean_snacku_b)];
//fprintfMat('training_color_snack1.txt', training_color_snacku);

test_mean_snacku_r = mean_snacku_r(6:10);
test_mean_snacku_g = mean_snacku_g(6:10);
test_mean_snacku_b = mean_snacku_b(6:10);
test_color_snacku = [test_mean_snacku_r, test_mean_snacku_g, test_mean_snacku_b];
//fprintfMat('test_color_snack1.txt', test_color_snacku);

// Test Object Classification

D_test_snack = [];
for l = 1:3
for m = 1:3
training_color_snack = fscanfMat('training_color_snack'+string(l)+'.txt');
m_snack = [mean(training_snack_A(:,l)), training_color_snack];

test_color_snack = fscanfMat('test_color_snack'+string(m)+'.txt');

for n = 1:5
x_snack = [test_snack_A(n, m), test_color_snack(n, 1), test_color_snack(n, 2), test_color_snack(n, 3)];
D_test = [x_snack - m_snack];
D_test_snack = [D_test_snack; sqrt(D_test*D_test')];
end

end

end

D_test_snacku = D_test_snack(1:15);
D_test_superthin = D_test_snack(16:30);
D_test_nova = D_test_snack(31:45);

D = [D_test_snacku, D_test_superthin, D_test_nova];

class = [];
for o = 1:15
[a, b] = find(D(o,:) == min(D(o,:)));
class = [class; b];
end

Activity 13 - Correcting Geometric Distortion

2009-08-18T08:57:00.093+08:00

The barrel geometric distortion in an image is removed in this activity based on the image's pixel coordinates and graylevel values.

The following distorted image of a grid is obtained from the internet.

Source: http://www.jmg-galleries.com/blog/2007/08/01/photo-term-series-post-13-barrel-distortion

The correction starts with the image's pixel coordinates. A subgrid from a distorted image can be thought of transformation of this subgrid from the ideal image through functions in the horizontal and vertical coordinates as shown below.

The transformation of a vertex from a subgrid is then illustrated as follows.

The following transformation functions are thus bilinear if the vertex from the ideal image is mapped simultaneously in both directions.

Eight equations are enough to determine the unknown eight coefficients c_1 to c_8 and these come from the two coordinates of the four vertices of a closed region in the grid both in the distorted and ideal image. Hence, a different set of coefficients is yielded for other closed region.

Initially from the distorted image obtained from the internet, a vertex from the most undistorted region is located. In this case, it is approximately near the image's center. Then, the length in pixels of the ideal rectangle is measured. The size of the ideal rectangle in the grid is found to be 28 x 28 pixels. Starting from the vertex, the pixel coordinates by counting 28 pixels to its left, right, top and bottom are determined and so the vertex points of the ideal grid are generated which is shown as follows.

Notice that the pixel coordinates at the edges are also obtained in order to cover the whole distorted image in correcting the barrel distortion. Meanwhile, the pixel coordinates of the vertices in the distorted image corresponding to the ideal image are gathered using the function locate of Scilab.
The sets of coefficients for each rectangle are then determined from the following equations.

The matrix expressions can be rewritten as follows.

So the coefficients are computed as shown below.

Now, per pixel inside a rectangle, the coordinates from the distorted image are attained using the corresponding set of coefficients in that rectangle through the following equations.

The grayscale values of each pixel are the next to be handled. At first, a matrix of zeros with the same size of the distorted image (150 x 151) is created. Then, interpolated graylevel values of each computed locations in the distorted image are copied to this blank matrix.
Two types of graylevel interpolation are implemented here. The first one is the nearest neighbor interpolation technique where the graylevel value of the nearest integer-valued coordinate is used and this just comes from the rounded pixel coordinates in the distorted image calculated using the eight coefficients.
Bilinear graylevel interpolation is the second technique employed to correct the barrel distortion. Since many of the computed pixel coordinates in the distorted image are not integer-valued, this technique is expected to result to a better corrected image compared to the nearest neighbor graylevel interpolation. Otherwise, the graylevel value at integer-valued locations is automatically copied to the blank matrix. An illustration of bilinear graylevel interpolation is shown below.

The pixel location pointed above is clearly not integer-valued. Through the equation below, the graylevel value is obtained at this location.

Notice that the equation above has form exactly the same in calculating the pixel coordinates in the distorted image. Thus, the same procedure applies in the bilinear graylevel interpolation which is based from the matrix expression shown as follows.

The coefficients a, b, c, and d must first be determined using four equations. These are then coming from the four nearest pixels of each computed pixel coordinate in the distorted image. The integer coordinates are assumed to be at the center of each pixel hence the four nearest pixels are located at the rounded values of 0.5 pixel away from calculated pixel location in the distorted image for both directions. The 4 x 4 matrix from the expression above are then assembled for each computed pixel coordinate in the distorted image which contains the location of the four nearest pixels. The following is the simplified version of the matrix expression above.

Then, the coefficients a, b, c, and d for each computed pixel location in the distorted image are acquired as shown below.

The interpolated graylevel value of each calculated pixel location in the distorted image is then transferred to the blank matrix for each pixel inside a rectangle in the ideal grid. The whole algorithm is repeated for the other rectangles in the ideal grid.
Recall that vertices at the perimeter of the ideal grid and distorted image are also considered. Thus to prevent incorrect indexing, the graylevel value is directly copied from the distorted image onto the black matrix at the pixel coordinate inside a rectangle in the ideal grid, if the calculated pixel location is less than 1 or greater than the length of the image which is 150 for the vertical direction and 151 for the horizontal direction.

The next figure is the result of the nearest neighbor graylevel interpolation technique applied to the barrel distorted grid image.

It can be observed that the lines are pixelated and jagged even though the curved grid lines are straightened.
Meanwhile, the figure below is the result of the bilinear graylevel interpolation.

Apparently, bilinear graylevel interpolation is better than the nearest neighbor interpolation technique. The barrel distortion is completely eliminated at the same time producing smooth grid lines.

I say that this activity has the most complicated code so far but since I have successfully corrected the barrel distortion of a grid image and understand how the whole process is implemented, I grade myself 10/10.

Earl has helped me in this activity and I thank Mandy for explaining how she has implemented the bilinear graylevel interpolation technique.

Appendix
The following Scilab code is utilized in this activity.

stacksize(4e7);

distort = gray_imread('barrel_distortion.jpg');

pixelx = 28;
pixely = 28;
tran = 20;

s = size(distort);
ideal = zeros(s(1), s(2));

xit = [];
yit = [];
for i = 1:5
for j = 1:5
xit = [xit, tran + pixelx*(j - 1)];
yit = [yit, tran + pixely*(i - 1)];
ideal(tran + pixelx*(i - 1), tran + pixely*(j - 1)) = 1;
end
end

xit = matrix(xit, 5, 5);
yit = matrix(yit, 5, 5);

xi = zeros(7, 7);
yi = zeros(7, 7);
xi(2:6, 2:6) = xit;
yi(2:6, 2:6) = yit;

xi(1,:) = 1;
xi(7,:) = s(1);
yi(:,1) = 1;
yi(:,7) = s(2);
ideal(1, 1) = 1;
ideal(1, s(2)) = 1;
ideal(s(1), 1) = 1;
ideal(s(1), s(2)) = 1;

for k = 1:5
xi(k + 1, 1) = xi(k + 1, 2);
xi(k + 1, 7) = xi(k + 1, 2);
yi(1, k + 1) = yi(2, k + 1);
yi(7, k + 1) = yi(2, k + 1);
ideal(xi(k + 1, 2), 1) = 1;
ideal(1, yi(2, k + 1)) = 1;
ideal(xi(k + 1, 2), s(2)) = 1;
ideal(s(1), yi(2, k + 1)) = 1;
end

//scf(0);
//imshow(ideal);
//imwrite(ideal, 'idealgrid.bmp');

//scf(1);
//imshow(distort, []);
//xdyd = (locate(49, 1))';
//fprintfMat('xdyd.txt', xdyd);

xdyd = fscanfMat('xdyd.txt');

xd = s(1) - xdyd(:,2);
yd = xdyd(:,1);

xd = (matrix(xd', 7, 7))';
yd = (matrix(yd', 7, 7))';

newgrid_bilinear = zeros(s(1), s(2));
newgrid_nearest_neighbor = zeros(s(1), s(2));
for a = 1:6
for b = 1:6

vertex_upperleft = [xd(a, b), yd(a, b)];
vertex_upperright = [xd(a, b + 1), yd(a, b + 1)];
vertex_lowerleft = [xd(a + 1, b), yd(a + 1, b)];
vertex_lowerright = [xd(a + 1, b + 1), yd(a + 1, b + 1)];

xd4 = [vertex_upperleft(1); vertex_upperright(1); vertex_lowerleft(1); vertex_lowerright(1)];
yd4 = [vertex_upperleft(2); vertex_upperright(2); vertex_lowerleft(2); vertex_lowerright(2)];

T = [xi(a, b) yi(a, b) xi(a, b)*yi(a, b) 1; xi(a, b + 1) yi(a, b + 1) xi(a, b + 1)*yi(a, b + 1) 1; xi(a + 1, b) yi(a + 1, b) xi(a + 1, b)*yi(a + 1, b) 1; xi(a + 1, b + 1) yi(a + 1, b + 1) xi(a + 1, b + 1)*yi(a + 1, b + 1) 1];

C1to4 = inv(T)*xd4;
C5to8 = inv(T)*yd4;

for m = (xi(a, 1)):(xi((a + 1), 1))
for n = (yi(1, b)):(yi(1, (b + 1)))

points_in_one_box = [m n m*n 1];

xnew = points_in_one_box*C1to4;
ynew = points_in_one_box*C5to8;
rxnew = round(xnew);
rynew = round(ynew);

if xnew < 1 | xnew > s(1) | ynew < 1 | ynew > s(2)
newgrid_bilinear(m, n) = distort(m, n);
newgrid_nearest_neighbor(m, n) = distort(m, n);
else

if xnew - rxnew == 0 & ynew - rynew == 0
newgrid_bilinear(m, n) = distort(xnew, ynew);
else

nearest_lowerleft = [round(xnew - 0.5), round(ynew - 0.5)];
nearest_upperleft = [round(xnew - 0.5), round(ynew + 0.5)];
nearest_lowerright = [round(xnew + 0.5), round(ynew - 0.5)];
nearest_upperright = [round(xnew + 0.5), round(ynew + 0.5)];

nearest_pixels = [nearest_lowerleft(1) nearest_lowerleft(2) nearest_lowerleft(1)*nearest_lowerleft(2) 1; nearest_upperleft(1) nearest_upperleft(2) nearest_upperleft(1)*nearest_upperleft(2) 1; nearest_lowerright(1) nearest_lowerright(2) nearest_lowerright(1)*nearest_lowerright(2) 1; nearest_upperright(1) nearest_upperright(2) nearest_upperright(1)*nearest_upperright(2) 1];
distortgray_lowerleft = distort(nearest_lowerleft(1), nearest_lowerleft(2));
distortgray_upperleft = distort(nearest_upperleft(1), nearest_upperleft(2));
distortgray_lowerright = distort(nearest_lowerright(1), nearest_lowerright(2));
distortgray_upperright = distort(nearest_upperright(1), nearest_upperright(2));
distortgray_nearest = [distortgray_lowerleft; distortgray_upperleft; distortgray_lowerright; distortgray_upperright];

abcd_coefficients = inv(nearest_pixels)*distortgray_nearest;
interpolated_gray = [xnew ynew xnew*ynew 1]*abcd_coefficients;

newgrid_bilinear(m, n) = interpolated_gray;
newgrid_nearest_neighbor(m, n) = distort(rxnew, rynew);

end

end

end
end

end
end

//scf(2);
//imshow(newgrid_bilinear, []);
//imwrite(newgrid_bilinear/max(newgrid_bilinear), 'newgrid_bilinear.bmp');

//scf(3);
//imshow(newgrid_nearest_neighbor, []);
//imwrite(newgrid_nearest_neighbor/max(newgrid_nearest_neighbor), 'newgrid_nearest_neighbor.bmp');

Activity 12 - Color Image Segmentation

2009-08-04T09:01:00.028+08:00

A region of interest (ROI) from a colored image is segmented in this activity using two types of segmentation - Parametric and Non-parametric.

The following image is searched from the internet.

Source: http://www.justfurfun.org/doghtml/yorikpage.htm

The blue ball which the dog holds is the ROI used for segmentation. A monochromatic patch is then cropped from the ROI and this is shown below.

This monochromatic patch is the basis of extracting the ROI from the whole image. Its normalized primary color channels red (R), green (G), and blue (B) separately with the whole image's are obtained next. This is done by dividing each color channel element-by-element to the sum of the RGB arrays resulting to the normalized chromaticity coordinates for red (r), green (g), and blue (b). The normalized chromaticity space is as follows.

Here, the 3 dimensional color space is reduced to only 2 having red as the horizontal axis and green as the vertical axis. Of course, blue has its region in the normalized chromaticity space since the relation between the normalized chromaticity coordinates is r + g + b = 1. Normalizing the color channels of the whole image and the patch of the ROI is very important because the shade or brightness of the colors of the whole image and the ROI is separated from their pure colors or chromaticities and so the image is segmented correctly.

Parametric segmentation is the first type implemented to the image. In this process, the probability of each pixel from the whole image belonging to the chromaticity of the ROI is determined through a Gaussian distribution having the form as shown below.

From the expressions above, sigma is the standard deviation of the normalized chromaticities of the patch of the ROI while mu is the mean. Then, r and g come from the normalized chromaticities of each pixel of the whole image. The segmented image is then generated from the multiplied or joint probabilities in the red and blue normalized channels. This is illustrated as follows.

Apparently, the white region represents the ROI that is the blue ball. Comparing from the original image with the blue ball, the image is well segmented since various shades of the ROI are indeed ignored and thus the whole blue ball is almost separated from the image.

Non-parametric segmentation operates in the probability distribution of the normalized chromaticity of the ROI alone without depending to a function such as the Gaussian distribution employed in the parametric segmentation. A 2 dimensional histogram is acquired here as shown below which is then verified from the normalized chromaticity space.

Both histograms above are comparable to the normalized chromaticity space, their bins are just varied for the sake of visualization (left histogram with 32 bins) and application of the histogram backprojection process (right histogram with 256 bins). The histogram agrees with the normalized chromaticity space since the peaks are located within the blue region knowing that the ROI has chromaticity within this area.
The histogram backprojection is where each pixel in the image is assigned to a value corresponding to the histogram value in the normalized chromaticity space. Hence, it is essential to increase the bins of the histogram to 256 in order to accommodate the 256 rounded integer values in the normalized chromaticity space. Finally, the result of the non-parametric segmentation is the following.

It can be observed that the image is not that well segmented compared to the parametric approach since majority of the blue ball particularly at the bottom part where its brightness is minimal is not extracted. The good thing about non-parametric segmentation is that no function is needed, instead, just the histogram of a monochromatic patch from the ROI is required. Therefore, for this case, a Gaussian function works great as basis of probability distribution for image segmentation.

Since I have successfully done and understand parametric and non-parametric segmentation of a colored image, I give myself 10/10 in this activity.

Throughout the activity, I have worked independently in its completion. However, I shared useful ideas that have helped my classmates.

Appendix
Below is the whole Scilab code utilized in this Activity.

stacksize(4e7);

image = imread('big_blue_ball.jpg');
I = imread('ROI.jpg');

R = I(:,:,1);
G = I(:,:,2);
B = I(:,:,3);

RGB = R + G + B;

r = R./RGB;
b = B./RGB;
g = G./RGB;

Ri = image(:,:,1);
Gi = image(:,:,2);
Bi = image(:,:,3);

RGBi = Ri + Gi + Bi;

ri = Ri./RGBi;
bi = Bi./RGBi;
gi = Gi./RGBi;

// Parametric Segmentation

rmean = mean(r);
gmean = mean(g);
rsigma = stdev(r);
gsigma = stdev(g);

rp = (1/(rsigma*sqrt(2*%pi)))*exp(-((ri - rmean).^2)/(2*rsigma^2));
gp = (1/(gsigma*sqrt(2*%pi)))*exp(-((gi - gmean).^2)/(2*gsigma^2));

rgp = round(rp.*gp);

//scf(0);
//imshow(rgp);
//imwrite(rgp, 'parametric segment.bmp');

// Non-parametric Segmentation
// 2D Histogram Plot

BINS = 256;
rint = round(r*(BINS - 1) + 1);
gint = round(g*(BINS - 1) + 1);
colors = gint(:) + (rint(:) - 1)*BINS;
hist = zeros(BINS, BINS);

for row = 1:BINS
for col = 1:(BINS - row + 1)
hist(row, col) = length(find(colors == (((col + (row - 1)*BINS)))));
end
end

histo = hist/max(hist);

//scf(1);
//mesh(histo);

// Backprojection

rib = ri*255;
gib = gi*255;

[a, b] = size(image);
segment = zeros(a, b);

for i = 1:a
for j = 1:b
c = round(rib(i, j)) + 1;
d = round(gib(i, j)) + 1;
segment(i, j) = hist(c, d);
end
end

//scf(2);
//imshow(segment);
//imwrite(segment, 'nonparametric segment.bmp');

Activity 11 - Color Image Processing

2009-07-30T09:33:00.053+08:00

After previous activities of processing grayscale and binary-type images, this time, color images are handled here specifically true color images which are taken by a digital color camera using its various white balancing (WB) settings.

In a colored digital image, each pixel from it is composed of the three primary spectral colors Red (R), Green (G), and Blue (B) which are overlaid in different proportions. These RGB values are given by the following equations.

From the expressions above, rho corresponds to the spectral reflectance of a surface, S is the spectral power distribution (SPD) of a Planckian light source illuminating the surface and nu is the spectral sensitivity of a color camera in the RGB channels. Meanwhile, K is the normalizing constant or WB constant dependent from the SPD of the light source and spectral sensitivity of the camera, again different in each color channel. The wavelength lambda ranges within the visible spectrum usually from 380 nm to 780 nm.

Initially, an image which includes the three major hues RGB and a known white object is taken by a digital color camera applying its different WB settings. The digital camera used is Fujifilm Finepix with maximum image resolution of 6.3 mega pixels while the WB settings include automatic, fine, shade, fluorescent 1, fluorescent 2, fluorescent 3, and incandescent. The exposure value (EV) of the camera is set to -2 ensuring that the RGB values of each pixel do not exceed the maximum pixel value of 1.
The following set of images are the pictures taken under constant fluorescent light source illumination implementing the 7 WB settings in order of automatic, fine, shade, fluorescent 1, fluorescent 2, fluorescent 3, and incandescent from top to bottom and left to right.

As expected, each WB setting of the camera results to different color rendering of the picture. Hence, the normalizing constant K varies with each WB setting of the camera and thus producing different RGB values in pixels. Of course, the automatic WB setting generates the best rendering of the picture since it automatically selects the most appropriate K values which adjust according to the light source illumination. The fine and shade WB settings are considerably good enough for the rendering while the three fluorescent WB settings approach the blue hue as the number increases from 1 to 3.
The last one which uses incandescent WB is observed to be the most obvious incorrectly balanced image. This is then enhanced using two automatic white balancing (AWB) algorithms.

The first one is the White Patch Algorithm where a patch from a known white object of the unbalanced image is taken and its averaged RGB values are determined. This is why the picture must have a known white object. Then, the RGB color channels of the unbalanced image are divided to the corresponding averaged RGB values of the white patch. The resulting RGB arrays are finally overlaid for the white patch balanced rendering of the wrongly balanced image. Note that the final RGB arrays in saving the rendered image are made sure not to exceed the value of 1 by dividing them to their maximum value even though initially there are no RGB values greater than 1 because of the EV of -2. This is because some RGB values less than 1 are divided to the averaged RGB values from the white patch which are less than them. Therefore, not doing this makes the rendering still incorrectly white balanced.
Shown below is the white patch selected from the picture taken with incandescent WB setting.

Then the following are the original picture and its rendering using the White Patch Algorithm.

The incorrectly balanced image using the incandescent WB setting is now appropriately white balanced using the White Patch Algorithm. The quality is as good as in using the automatic WB setting of the camera.

The second AWB algorithm is called the Gray World Algorithm. Here, the average pixel values of each RGB color channel are obtained. The RGB arrays of the unbalanced image are then divided to the corresponding average RGB values. Hence, this algorithm works even the picture does not have a known white object. Again, in saving the rendered image, the values are clipped to 1 like in the White Patch Algorithm.
The original picture taken with incandescent WB setting and its rendering using the Gray World Algorithm are as follows.

The name of the Gray World Algorithm itself implies the effect when applied to an image. The average values taken from each color channel here represent the RGB values of white however these are already weighted. Hence, when the RGB arrays are divided by these average values, gray which is between white and black becomes the basis of balancing, and thus the image enters the gray world.

As another illustration of the two AWB algortihms, the picture taken with fluorescent 3 WB setting is rendered using these algorithms. First below is the white patch chosen.

The next set of images consists the original picture, its rendering using the White Patch Algorithm, and using the Gray World Algorithm.

It can be observed that the rendered picture taken with fluorescent 3 WB setting using the White Patch Algorithm is almost the same as in the rendered picture taken with the incandescent WB setting. However, the rendering using the Gray World Algorithm in this WB setting is better than in the incandescent type.

For a more conclusive observation, a picture is taken again but this time the objects have colors belonging to a primary hue. Incandescent WB setting is employed here since this is the most inappropriate WB setting for the illumination condition. Some red objects are gathered together and of course a known white object is included in the group.
The white patch selected is shown below.

The following are the original picture, its rendering using the White Patch Algorithm, and using the Gray World Algorithm.

Apparently, the White Patch Algortihm is better than the Gray World Algorithm. This is due to the fact that the picture is indeed white balanced since the average RGB values where the RGB color channels are normalized, come from a known white object. While the Gray World Algorithm just uses the average RGB values of each color channel of the unbalanced image in effect just balancing it to gray. However, the Gray World Algorithm generates good rendering if the WB setting used to take the picture blends well with the illumination condition.

Since I successfully processed some true color images taken with different WB settings of a digital color camera using the two AWB algorithms , I give myself 10/10 in this activity.

I am able to complete this activity successfully through discussions with Gary and Raffy.

Appendix
Below is the Scilab code utilized in this activity.

stacksize(4e7);

imageRGB = imread('imageRGBincan.jpg');

//scf(0);
//imshow(imageRGB);

R = imageRGB(:,:,1);
G = imageRGB(:,:,2);
B = imageRGB(:,:,3);

// White Patch Algorithm

white_patch = imread('white patch RGBincan.jpg');

Rw = sum(white_patch(:,:,1))/length(white_patch(:,:,1));
Gw = sum(white_patch(:,:,2))/length(white_patch(:,:,2));
Bw = sum(white_patch(:,:,3))/length(white_patch(:,:,3));

nR = R/Rw;
nG = G/Gw;
nB = B/Bw;

newimageRGB = [];
newimageRGB(:,:,1) = nR;
newimageRGB(:,:,2) = nG;
newimageRGB(:,:,3) = nB;

//scf(1);
//imshow(newimageRGB);
//imwrite(newimageRGB/max(newimageRGB), 'newimageRGBincan.jpg');

// Gray World Algorithm

Rgray = sum(R)/length(R);
Ggray = sum(G)/length(G);
Bgray = sum(B)/length(B);

ngrayR = R/Rgray;
ngrayG = G/Rgray;
ngrayB = B/Rgray;

newimagegrayRGB = [];
newimagegrayRGB(:,:,1) = ngrayR;
newimagegrayRGB(:,:,2) = ngrayG;
newimagegrayRGB(:,:,3) = ngrayB;

//scf(2);
//imshow(newimagegrayRGB);
//imwrite(newimagegrayRGB/max(newimagegrayRGB), 'newimagegrayRGBincan.jpg');

Activity 10 - Preprocessing Text

2009-07-28T09:46:00.021+08:00

Another activity applying various image processing techniques is done here. This is divided into two major parts. First, a handwritten text from a scanned image is preprocessed. In the second part, incidences of a typewritten word are found in the same scanned image.

The image of a scanned demo checklist form is initially downloaded given by the activity. This is shown below.

A portion of handwritten text along the horizontal lines is then cropped from this image. Since the text image is tilted, its Fourier transform (FT) is determined in order to obtain the angle of rotation. Either way, this angle can just be set through trial and error. In this case, the text image is rotated by 1 degree in the clockwise direction using the function mogrify in Scilab. The following are the cropped text image and its rotation in grayscale.

The text image arranged from top reads: VGA Cable, Power Cord, Remote Control, RCA Cable, USB Cable. The primary reason of rotation is for the ease enhancement of the text image since the horizontal lines must be eliminated and hence a filter that removes this pattern can easily be created. Next, the FT of the rotated text image is attained and this serves as the template for making the filter. Henceforth, frequencies of unwanted pattern can be blocked here. The FT of the rotated text image and the filter created are shown as follows.

Since the horizontal lines are to be diminished, the central vertical line frequencies must be blocked. However, its center must not be obstructed because it contains the primary information of the rotated text image.
Below is the filtered rotated text image.

The text can now be extracted from the background after filtering. From the histogram plot of the grayscale values of the filtered rotated text image using GIMP, the threshold value for its conversion to binary type is obtained. The following are the binarized filtered rotated text image and its inversion.

The threshold value used for binarizing is 0.275. The inversion of the binarized filtered rotated text image is essential since the final step in preprocessing the text is the utilization of image morphological operations. During these operations, the region of interest, which is the text, must be the foreground with pixel value 1 while the background must be 0. Notice that some traces of the horizontal lines reside in the inverted binarized filtered rotated text image. Therefore, a closing operation is needed for these to be removed. Recall that this morphological operation is equivalent to the dilation of the eroded image by a structuring element to that same structuring element. Its implementation to the inverted filtered rotated text image is applied with a 3 x 1 sized matrix of ones as the structuring element. Its effect is shown below.

The image above is now the preprocessed handwritten text from the scanned image of a demo checklist form. It can be noticed that the remnants of horizontal lines are eliminated through the closing operation. This is because the structuring element is enough to reconnect two clusters along the vertical direction. Hence, increasing the number of rows in the matrix of the structuring element makes the letters that are few pixels adjacent vertically to each other cluster and this is not good for preprocessing the text.
For further analysis of the preprocessed text, the built-in function bwlabel in Scilab is again used here for labeling the clusters which are the letters reconstructed. There are 54 clusters detected and originally there are only 46 letters in the text image. This result is fair good enough because in the first place, some handwritten letters in the original text image are not distinctive. Some also coincide with adjacent letters and are not written in the same magnitude.

The incidences of the word DESCRIPTION from the same scanned image of a demo checklist form are found in the second part of this activity by correlation.
This time the whole image is needed since it is required to locate all incidences of the word from the scanned image. Now, this is rotated so that a sample image of the word can easily be extracted that is to be used in correlation. The original scanned image and its rotation in grayscale are as follows.

Again, the whole image is rotated by 1 degree in the clockwise direction using the function mogrify in Scilab like the text image in the first part. Then, the rotated image is binarized base to the threshold value from the histogram of its grayscale values obtained using GIMP. This is illustrated below.

The threshold value employed in binarizing the rotated image is 0.49.
Due to the fact that it is customary in a binary image where the region of interest must be white value and the background to be black, then the binarized rotated image is inverted. Shown in the following are the inverted binarized rotated image and the sample image of the word DESCRIPTION extracted from it.

In getting a sample image of the word DESCRIPTION from the inverted binarized rotated image, it is important that it is placed in an image with black background and this sample image must have the same size as the inverted binarized rotated image because correlation occurs at the frequency domain of their FTs.
Recall that correlation works here when the FT of the sample image is multiplied element-by-element to the conjugate of the FT of the inverted binarized rotated image. Then, the modulus of the shifted inverse FT of their product displays the correlated image. Now, the following is the result of correlation of the sample image of the word DESCRIPTION to the inverted binarized rotated image.

Apparently from the correlated image above, there are three incidences of the word DESCRIPTION in the inverted binarized rotated image which indeed agrees comparing from the scanned image of a demo checklist form. The incidences are marked by the bright spots as observed in the correlated image. One is found at the upper left and two are located below it and are arranged in a same line horizontally.

Although the preprocessing of a handwritten text in the first part is not perfectly done, I have applied various techniques of image processing successfully in this part also with the correlation of the word DESCRIPTION to the scanned image in the second part, hence I grade myself 10/10 in this activity.

This activity is successfully completed through my collaboration with Gary and Ed.

Appendix
The following is the source code for this activity.

stacksize(4e7);

// Preprocessing Handwritten Text

text = gray_imread('htext.bmp');
//scf(0);
//imshow(text);

textrot = mogrify(text, ['-rotate', '1']);
//scf(1);
//imshow(textrot);
//imwrite(textrot, 'htextrot.bmp');

Ftext = log(abs(fftshift(fft2(textrot))));
//scf(2);
//imshow(Ftext, []);
//imwrite(Ftext/max(Ftext), 'Fhtextrot.bmp');

filter = gray_imread('filterhtextrot.bmp');
filtext = abs(ifft(fftshift(filter).*fft2(textrot)));
//invfiltext = max(filtext) - filtext;
//scf(3);
//imshow(filtext, []);
//imwrite(filtext/max(filtext), 'filteredhtextrot.bmp');

bintext = im2bw(filtext, 0.275);
//imwrite(bintext, 'bintext.bmp');
invbintext = 1 - bintext;
//scf(4);
//imshow(invbintext);
//imwrite(invbintext, 'invbintext.bmp');

SE = ones(3, 1);

morphtext = dilate(invbintext, SE);
morphtext = erode(morphtext, SE);

//scf(5);
//imshow(morphtext);
//imwrite(morphtext, 'morphtext.bmp');

[L, n] = bwlabel(morphtext);

// Typewritten Text Correlation

checklist = gray_imread('Untitled_0001.jpg');
checklistrot = mogrify(checklist, ['-rotate', '1']);
//scf(6);
//imshow(checklistrot);
//imwrite(checklistrot, 'checklistrot.bmp');

binchecklist = im2bw(checklistrot, 0.49);
//scf(7);
//imshow(binchecklist);
//imwrite(binchecklist, 'binchecklist.bmp');

invbinchecklist = 1 - binchecklist;
//scf(8);
//imshow(invbinchecklist);
//imwrite(invbinchecklist, 'invbinchecklist.bmp');

description = imread('description.bmp');
cor = abs(fftshift(fft2(fft2(description).*(conj(fft2(invbinchecklist))))));
//scf(9);
//imshow(cor, []);
//imwrite(cor/max(cor), 'correlated text.bmp');

Activity 9 - Binary Operations

2009-07-23T10:00:00.034+08:00

This activity requires to give the best estimate of the area of a circular punched paper which is called here as cell. Some techniques of image processing that have been worked on in the past activities are employed here in order to calculate the area of this cell.

First, a scanned image of scattered punched papers is selected from the two images given by the activity. This is shown below.

The image above has size 748 x 618 and it is subdivided into 9 subimages each of which has size 256 x 256 pixels. Thus, some subimages overlap with the other which is better since it results to more area to be calculated and thus the result to be more accurate. Paint is used in cropping the scanned image.
A representative subimage that is the parcel labeled 'circles_2' located at the middle of the one-third topmost part of the whole original image is illustrated as follows.

The histogram of the grayscale values of each subimage is then analyzed using GIMP. This is where the threshold value is determined for the binarizing of the 9 subimages, that is, they only possess two values which are 0 for black corresponding to the background and 1 for white corresponding to the foreground or called the Region of Interest (ROI).
A threshold value of 0.85 is used in making the 9 subimages be of Binary-type. The binarized subimage 'circles_2' is the following.

Notice that binarizing alone is not enough to get the ROI for area calculation. Thus some image morphological operations are utilized in order to remove unnecessary noise and hence to connect overlapping cells which are coined blobs.
Closing and opening operators are then implemented here. The closing operator is defined to be the erosion of the dilated image A by a structuring element (SE) B to that same SE B. Meanwhile, the opening operator is equivalent to the dilation of the eroded image A by a SE B to that same SE B. Therefore, a series of dilation, erosion, erosion, and dilation operation is applied to each subimage for its further enhancement.
The SE used takes the form of a small circle as the matrix below.

The effect to the representative binarized subimage 'circles_2' after applying the closing and opening operations is shown as follows.

It can be observed that noise are eliminated, blobs are created due to the overlapping of some cells and in the same time single-cell units are more defined.
After applying the methods presented above to all the subimages, the blobs and the single-cell units are then isolated and labeled in preparation for area calculation using the built-in function in Scilab namely bwlabel. Thus, the area, which is the total number of white pixels, of blobs and single-cell units to all subimages are independently determined and gathered.

Now, the histogram of the calculated areas is shown below.

The occurrence of the area of a single-cell unit is found to be between 500 and 550 square pixels. Of course, the blobs and other cells that are not fully reconstructed give very large and very small areas respectively. Hence, these are eliminated in the process of determining the mean area and the standard deviation of the distribution.

The calculated mean area of the cell is 519.79 square pixels. This is derived from the sum of the elements within the interval of 500 to 550 square pixels divided by the number of elements in this range. Using the built-in function stdev in Scilab, the standard deviation attained is 12.12.

For validity of the result, the mean area is compared to the area of a single-cell unit which is not overlapping to other cell(s) and this is cropped from a subimage. This is also binarized but the closing and opening operations are not already implemented. The following images represent the basis of comparison to the calculated mean area.

Pixels of the cell are then summed for its area and the result is 517 square pixels. Finally, its percent difference to the calculated mean area is 0.54 %.

Since I have applied different image processing techniques for this activity to be successful, and obtained very small or almost insignificant percent difference for the area of a cell, I grade myself 10/10.

Before starting, Mimie has shared some insights to me on how this activity is handled. However, I have worked independently throughout its completion at the same time sharing my ideas to my classmates after I have finished this.

Appendix
The following Scilab code below is utilized in this activity.

circles = [];
binarycircles = [];
N = 9;
for i=1:N
circles1 = gray_imread('circles_'+string(i)+'.jpg');
binarycircles1 = im2bw(circles1, 0.85);
circles = [circles, circles1];
binarycircles = [binarycircles, binarycircles1];
end

M = 2;
graysubimage = circles(:,((M*256) - 255):(M*256));
binarysubimage = binarycircles(:,((M*256) - 255):(M*256));

//scf(0);
//imshow(graysubimage);

//scf(1);
//imshow(binarysubimage);
//imwrite(binarysubimage, 'binarysubimage2.bmp');

SE = [0 1 0; 1 1 1; 0 1 0];
newbinarysubimage = dilate(binarysubimage, SE);
newbinarysubimage = erode(newbinarysubimage, SE);
newbinarysubimage = erode(newbinarysubimage, SE);
newbinarysubimage = dilate(newbinarysubimage, SE);

//scf(2);
//imshow(newbinarysubimage);
//imwrite(newbinarysubimage, 'newbinarysubimage2.bmp');

[L, n] = bwlabel(newbinarysubimage);

A = [];
for j = 1:n
b = (L == j);
A(j) = sum(b);
end

Area = fscanfMat('area data.txt');

//scf(3);
//histplot(length(Area), Area);
//title('Histogram of Calculated Cell Areas');
//xlabel('Area');
//ylabel('Frequency');

index_goodArea = find(Area>500 & Area<550);
mean_Area = sum(Area(index_goodArea))/length(index_goodArea);
stddev_Area = stdev(Area(index_goodArea));
gray1circle = gray_imread('1circle.bmp');

//scf(4);
//imshow(gray1circle);

binary1circle = im2bw(gray1circle, 0.85);

//scf(5);
//imshow(binary1circle);
//imwrite(binary1circle, 'binary1circle.bmp');

theo_Area = sum(binary1circle);
percent_error = abs((theo_Area - mean_Area)/(theo_Area))*100;

Activity 8 - Morphological Operations

2009-07-21T09:30:00.058+08:00

We have applied the basic image morphological operations dilation and erosion to some basic shapes in this activity. In Scilab, there are already built-in image morphological functions namely dilate and erode which are used here. Also, different structural elements are implemented in the dilation and erosion of images.

Dilation is defined by the expression

It states that the intersection of the translation of the reflection of B translated at z for all z's to A is not an empty set. Here, A is the image to be dilated by the structuring element B and z is the coordinate of image A.
Image dilation operates in the following algorithm. A certain location of B, which in this activity is its geometric center, is to be superimposed to every z in A. If at least one element of B coincides with the foreground of A, then the location in A being superimposed by the geometric center of B is to be replaced by one. Thus, if no elements of B coincide with any element in the foreground of A, then it is retained to zero. Recall that the background of A is equal to zero and the object or its foreground under morphological operation is equal to one.
Meanwhile, erosion is defined as the expression below

This means that all elements of B translated at z for all z's are contained in A.
Image erosion then has the following algorithm. The geometric center of B is to be superimposed to all z in A. If all elements of B are contained in A then the location in A being superimposed by the geometric center of B is to be retained to one and otherwise, if at least one element of B coincides with the background of A, it is replaced by zero.

The shapes that have been created and morphologically operated are a square with size 50 x 50, a right triangle with base 50 pixels long and height 30 pixels long, a circle with 25 pixels radius, a hollow square with size 60 x 60 and its edges are 4 pixels thick, and a plus sign with 50 pixels long for each line having thickness of 8 pixels. They are all made in Paint.
Meanwhile, the following structuring elements are used in dilating and eroding the images. A 4 x 4 ones, 2 x 4 ones, 4 x 2 ones, and a cross with 5 pixels long and 1 pixel thick. These are generated in Scilab as matrices shown below respectively.

Before simulating the effects of the morphological operations in the images, we have first predicted the outcome and draw this in yellow paper which is submitted to our professor Dr. Maricor Soriano.
Now, the following sets of figures are the results of the simulated morphological operations in the five images.

Square (Original Image)

Dilation

Structuring Element: 4 x 4, 2 x 4, 4 x 2, cross