Pure awesomeness..

Activity 19 - Restoration of blurred image

2009-10-09T10:42:00.000-07:00

In this activity, we are tasked to perform restoration algorithm to motion blurred images with additive noise.

Model

Degradation process is modeled as a degradation function (see figure below), together with an additive noise, applied to an input image.

The blurring process is taken as the convolution of the degradation function h(x,y) with the input image f(x,y). With the addition of noise, the output image becomes:

g(x,y) = h(x,y)*f(x,y) + η(x,y)

We take the Fourier transform (FT) of the output function which results to a simple algebraic form:

G(u,v) = H(u,v)F(u,v) + N(u,v)

Now, the degradation transfer function H(u,v) is given by:

H(u,v) =[^T/_π(ua+vb)] sin [π(ua+vb)]e^-jπ(ua+vb)

where a,b are the displacement of the image along the x and y direction, T is the exposure time of the capture (shutter speed).

Below is an example of an blurred image (a = 0.1, b = 0.1, T = 1) with an additional Gaussian noise (μ = 0.01,σ = 0.001).

Figure 1. (a) Image of Hilary Duff[3] and (b) the blurred image using a = b = 0.1, T = 1

In restoring blurred images, we use the Minimum square error (Weiner) filtering. This technique incorporates the degradation function (blurring process) and the noise to reconstruct a corrupted image (f) and finding an estimate f^{^} such that, as the name implies, the mean square error between them is minimum. The best estimate of the image in the frequency domain is given by:

where S_η = |N(u,v)|² (power spectrum of the noise) and S_f = |F(u,v)|² (power spectrum of the undergraded image).
Using the above equation, we demonstrate the following restorations at various exposure times (T) and motion blur displacements (a,b).

Figure 2. Restored images (down) increasing motion displacements a = b = 0.01, 0.1, 1 and (from left) increasing T (0.1, 1, 5, 10 and 50).

From the figure, we observed that increasing T results to a better restoration of the image. As previously mentioned, T is the exposure time of capture. It measures the amount of light coming into the sensor of an imaging device. This is the reason why we notice that images using small T tend to be relatively darker than when using larger T. Moreover, T also affects how movements appear in a picture[2]. Less Ts are intended for freeze fast moving objects, such as in sports. High Ts, on the other hand, are intended for intentional blurring of images for artistic purposes[2]. Our results suggest that higher T renders better image restoration. At these values, we gain more information of the blurring process; thus we are able to have improved estimate of the original image.
Also, as seen from the figure, the quality of the restored image decreases with increasing motion blur displacements. As expected, increasing the values of a and b enhance the blur resulting to an inevitably poor restoration of the image.

Most of the time, given an image regardless as being blurred or not, we do not have an idea of the distribution of the noise contained in that image. Hence, the equation above cannot be solved if we do not know the power spectrum of the noise. We thus approximate the equation by replacing the last term of the denominator with a specified constant K (shown below).

We explored on the effect of the value of K to the quality of the restoration. The figure below shows the restored images as K varies. As seen, the restoration quality decreases with K. From the equation above, K replaces the ratio S_η/S_f (the ration of the spectrum of noise and the spectrum of the undergraded image). K determines the relative degree of noise with respect to the image. The choice of higher values of K comes with the assumption that the image contains a high degree of noise. Hence, increased K (K higher than the actual S_η/S_f) implies lesser possibility of rendering a good result.

Figure 3. Restored images (down) increasing motion displacements a = b = 0.01, 0.1, 1 and (from left) increasing K (0.00001, 0.001, 0.1, and 10).

Note: Most of the results presented were qualitatively described i.e, the qualification of a good or bad quality of the images is subjective. A better and meaningful presentation of the results is to take the Linfoot's criterion (FCQ) as a qualifier and plot it with the parameters a,b, K and T.

In this activity, I give myself a grade of 10 for performing the activity well.

I thank Kaye, Irene and Mark Jayson for useful discussions.

References:
[1] Applied Physics 186 Activity 19 manual
[2] http://en.wikipedia.org/wiki/Shutter_speed
[3] http://www.smashingapps.com/wp-content/uploads/2009/04/grayscale-technique.jpg

Activity 18 - Noise model and basic image restoration

2009-09-24T07:40:00.000-07:00

Several factors affect the quality of a digital image. For one, the lenses of the camera, especially for single-lens cameras, distort the captured image. In our previous activity, we performed image processing techniques to correct such distortions. Aside from distortion, noise which is random variation of brightness and color information of images[2], also cause undesirable degradation to the captured image that usually arises during acquisition. For instance, temperature which alters the response of a camera sensor affects the amount of noise in an image[1]. In restoring such images, we use our knowledge of the degradation phenomenon, model it and perform an inverse process in order to recover the image.

Noise are composed of random variables characterized by certain probability distribution function (PDF)[1]. In this activity, we explore on most common noise models found in image processing applications as well as some basic noise-reduction filters used in image restoration.

Noise Models

Gaussian noise:

p(z) = (¹/√_2πσ )e^{^-(z-μ)²/_2σ²}

Rayleigh noise:

p(z) = (²/_b)(z - a)e^{^-(z-a)²/_b} for z > a else p(z) = 0

Erlang (Gamma) Noise:

p(z) = [^{a^bz^b-¹}/_(b-1)!] e^-az for z > a else p(z) = 0

Exponential Noise:

p(z) = ae^-az for z > a else p(z) = 0

Uniform Noise:

p(z) = ¹/_(b-a) if a < z < b else p(z) = 0

Impulse (salt an pepper) Noise:

p(z) = P_a for z = a
p(z) = P_b for z = b
p(z) = 0 otherwise

The figure below demonstrates the effect of each noise PDF to an image. Here we initially observe a grainy effect of the noise.

Figure 1. (left) Original image and (right) the image added with (clockwise from upper left) Gaussian, gamma, exponential, Rayleigh, uniform and, salt&pepper noise.

From the original grayscale value histogram of the image (shown in Figure 2), upon application of the noise to the image, the histogram was transformed as in Figure 3.

Figure 2. Grayscale value histogram of the image

Figure 3. Histogram of grayscale values from the generated noisy images: (clockwise from upper left) Gaussian, gamma, exponential, Rayleigh, uniform and, salt&pepper noise.

Basic noise-reduction filters

These filters takes a window S_xy in an image with center x,y and perform averaging.
Arithmetic mean filter (AMF)

Geometric mean filter (GMF)

Harmonic mean filter (HMF)

Contraharmonic mean filter (CMF)

CMF is good in removing salt&pepper noise. If Q>0, it removes pepper noise. If Q<0, it removes salt noise.

Restoration

We have tested the effect of the window size to the quality of the reconstructed image. The figure below shows the reconstructed images using the four filters for sizes of S_xy = {3x3, 5x5, 9x9}. We observe that smaller window size resulted to a relatively grainy restoration than with that of larger window size. On the other hand, large window size resulted to a more blurry image. In the following reconstructions, we used a window size of 5x5.

Figure 4. Reconstructed image of Gaussian-noised image using (from top to bottom) AMF, GMF, HMF and CMF (Q= 1) at different window sizes (left to right) 3x3, 5x5, 9x9.

Here we show the effects of the filters to the noised images.

Figure 5. (right)Reconstruction of (left)noisy images: (top to bottom) Gaussian, gamma, exponential, Rayleigh, uniform and, salt&pepper noise using (from 2^nd column to the right) AMF, GMF, HMF and CMF (Q=1).

From Fig. 5, we see that for all noised images, the AMF produced a consistent results although may not be the best result. This technique can thus be use generically for noisy images. GMF also produced a very good restoration for all noised images except for the salt&pepper-noised image. HMF is also resulted to a good reconstruction except for Guassian and salt&pepper-noised images. CMF (Q= 1) produced a good result except for exponential and salt&pepper-noised image. Salt&pepper-noised image was harder to restore compared to others. Moreover, we see that CMF (Q=1) was not able to correctly restore a salt&pepper-noised image, as previously said. Now we demonstrate the effect of CMF to salt-only and pepper-only-noised image.

Figure 6. Image reconstruction of Q = 1 and Q= -1 CMF for (top to bottom) the image with salt and pepper, salt-only and pepper-only noise. CMF cannot properly reconstruct a salt& pepper-noised image. Salt is removed when Q is negative. Pepper is removed when Q is positive.

We indeed see that CMF can perfectly restore images with salt noise and pepper noise alone.

Figure 7. Test image taken from [3]

Below are reconstructions of noised images using (from upperleft clockwise) AMF, GMF, CMF, and HMF.

Figure8. Reconstruction of a Gaussian-noised image.

Figure 9. Reconstruction of a gamma-noised image.

Figure 10. Reconstruction of a exponential-noised image.

Figure 11. Reconstruction of a Rayleigh-noised image.

Figure 12. Reconstruction of a uniform-noised image.

Figure 13. Reconstruction of a salt&pepper-noised image.

In this activity, I give myself a grade of 10 for finishing it well.
I thank Irene, Janno, Martin, Master, Cherry for the help.

References:
[1] AP186 Activity 18 manual
[2] http://en.wikipedia.org/wiki/Image_noise
[3] http://bluecandy.typepad.com/.a/

Activity 17 - Photometric Stereo

2009-09-08T01:00:00.002-07:00

Objects appear different depending on where the light source illuminating it is located. Shadows form at different locations. In this activity, we demonstrate a shadow model and extract shape from the shadows. Models depend on the type light source (point source, line source, etc.) used and where the light source is located. One of the more simple model is when the light source is placed at infinity. In general, when the light source is nearby, the brightness of an object illuminated by it is given by:

B(P) = [ρ(P)n(P)•S(P)]/R(P)²

where B(P) is the brightness of the object, S(P) is a vector from the light source and R(P) is the distance of the light source from the object.
In an infinite light source, 1/r dependence of B(P) disappears and the vector S(P) becomes constant for every point in the object. Thus the equation becomes:

B(P) = ρ(P)n(P)•S_o

In capturing an image of an object, we assume that it is proportional to the brightness of the object:

I(x,y) = kB(x,y) = kρ(x,y)n(x,y)•S_o = g(x,y)•V_o

where g(x,y) = ρ(x,y)n(x,y) and V_o=kS_o.

When an object is illuminated by light N light sources, Vo takes the form:

Each row represents the position (x,y,z) of the source in 3D space. Given an image with its intensity equal to its grayscale value, we can determine g(x,y) using the equation above. If we took a N captures of the surface of object illuminated by N light sources, then for each point (x,y):

I₁(x,y) = V₁₁g₁ + V₁₂g₂ + V₁₃g₃
I₂(x,y) = V₂₁g₁ + V₂₂g₂ + V₂₃g₃
I_N(x,y) = V_N1g₁ + V_N2g₂ + V_N3g₃

or in matrix form:

I = gV

Since we know I and V, we can calculate g using:

g = (V^TV)^-1VI

The normal vector can be solved by:

To get the surface elevation of point (u,v) on the surface, we use:

where

In this activity, we captured four images of a hemisphere illuminated by a far away light source located at four positions V:

(0.085832, 0.17365, 0.98106)
(0.085832, -0.17365, 0.98106)
(0.173650, 0.0, 0.98481)
(0.163180, -0.34202, 0.92542)

Below are the images taken for each light source locations.

Figure 1. Images of a hemispheres illuminated by a light source located at different positions

Using the procedure above, we obtained a 3D image of the surface shown in Fig.2.

Figure 2. 3D reconstruction of the hemisphere

Here, we have successfully reconstructed the image. The reconstruction, however, is not smooth. Some parts do not follow the contour of a hemisphere.

In this activity, I give myself a grade of 10 for a successful reconstruction of the surface.

I thank Irene, Alva, Orly for very very useful help. :)

Activity 16 - Neural Networks

2009-09-08T01:00:00.001-07:00

In the previous activities, we have used Euclidean distance and LDA as classifiers for pattern recognition. In this activity, we tried the neural networks. A neural network is a computational model of how neurons in the brain works. Basically, just as the brain does pattern recognition, the neural network is fed with features of known patterns; in which case, it "learns" the rules for classification from the input it receives.
Just as the neurons pass electrical signals to other neurons, this model (McCoulloch, Pitts 1943) essentially works the same way. With input signals x_i, (see Fig.1) a neuron receives weighted inputs(w_i*x_i) and lets the sum of all the signal it receives act on an activation function g. The resulting signal z is then passed to other neurons[1].

Figure 1. Artificial neuron

The same process takes place as the signal is passed from neuron to neuron. This network of neurons is called the neural network.

Figure 2. Neural Network

The neural network is composed of an input layer, a hidden layer, and the output layer. The input layer receives the signal (features) which is passed onto the hidden layer. The hidden layer then passes it to the output layer[1].

Results

Using a code originally written by J. Tugaff [2], we use neural network to classify objects as 1-peso coin, 25-cents coin, long leaf, short leaf or flower using r-g chromaticity and dimension ratio (length/width) as features. (Code needs modification depending on the number of classes and features used.)

Learning process
As training , we used the features of 5 samples/class as inputs to the neural network. It is trained such that it outputs: [1 0 0 0 0] for 1-peso coin, [0 1 0 0 0] for 25-cents coin and so on.

Figure 3. Plot of the distinguishing capability of the neural network with (left) learning rate and (right) training cycles T.

Before the classification, we checked the distinguishing capability (the ability to properly separate classification from other classes)of the neural network as a function of the learning rate and training cycles T. We checked it by using the training set as our test set. We took the average value of the "supposedly 1" outputs. Fig.3(left) shows the average against learning rate. An increasing trend is observed. As learning rate increases, the output approaches to 1 thus the network is able to classify the objects properly. We also observed an increasing trend when T increases. High T improves learning process resulting to better classification.
Based on these results, it is quite a good idea to choose learning rate and T as high as possible. However, this is very time-consuming to do. This will slow down the learning process which makes it impractical.¹ For our case, we use T = 5000 and learning rate = 1 for training.

After training, we tested the training set to verify that the network has "learned". The result of classification is summarized on the table below.

Table 1. Classification of the training set used in the learning process of the neural network.

Here we see that the network is able to identify properly identified the training set hence the network has already learned, we can now use it to classify our test objects.

Classification of test objects
In the table below, we show the result of classification of the test objects. Here, we have achieved a 96% correct classification.

Table 2. Classification of test objects using the neural network

¹ If one, on the other hand, is willing to wait for better results, he might as well choose higher learning rate and T.

I thank Master, Kaye and Thirdy and Irene for very useful help.

I give myself a grade of 10 for doing a good job.

References:
[1] Applied Physics 186 Activity 16 manual
[2] http://cole-ap186.blogspot.com/

Activity 15 - Probabilistic Classification

2009-09-07T07:13:00.000-07:00

In the previous activity, we tried to classify objects into classes based on a set of features that describe the objects. Particularly, we used the Euclidean distance as a classifier for our dataset. In this activity, on the other hand, we use another classifier called the Linear Discriminant Analysis (LDA). Discriminant analysis is a statistical approach in pattern recognition. LDA assumes that the classes of objects are linearly separable (as the term linear applies)-- meaning, the classes of objects can be separated by linear combination of features that describe the object[1]. The LDA formula is given by:

f_i = u_iC^-1x_k^T - 1/2u_iC^-1u_i^T + ln(p_i) .

where u_i is the mean features of the class i, C is the covariance matrix obtained from the training set, x_k is the features of test k^th test object, and p_i is the prior probability of a class i.

Mean Features of a Class

Given features of a set of test objects of the classes 1 and 2, (in this case, we have 3 test objects and 2 features for class 1)

x₁ = [a_o b_o; a₁ b₁; a₂ b₂] and x₂ = [c_o d_o; c₁ d₁]

mean feature of class 1 is calculated as:

u₁ = [mean(a_o,a₁,a₂) mean(b_o,b₁,b₂)]

Covariance Matrix
From the features of the whole training set (contains all training samples of all classes), one can calculate the global mean vector. This is essentially the mean features of the whole training set.
For example, given a data set

x = [a_o b_o; a₁ b₁; a₂ b₂; c_o d_o; c₁ d₂] ,

the global mean vector is given by

u = [mean(a_o,a₁,a₂,c_o,c₁) mean(b_o,b₁,b₂,d_o,d₁)]

From the u, we can solve the mean corrected data x_i^o for each class i

x_i^o= x_i -u .

The covariance matrix of a class i is:

c_i = [(x_i^o)^Tx_i^o]/n_i

where n_i is the number of test samples used for class i.
The covariance matrix of the whole test set is then solved using:

C = 1/n Σ n_i⋅c_i

where n is the number of samples in the whole data set.

Prior Probability

The prior probability of class i is given by:

p_i = n_i/n .

Results

We used the LDA formula in classifying objects in our data set. The object k is classified as belonging to class i if it resulted the highest f at i. Below is the result of our classification. Highlighted values are the maximum f_i obtained per sample. Samples are classified according to i where highlighted values occur.

Table 1. LDA results of 5 test samples of 1-peso coins,25-cent coins, long leaves, short leaves and flowers.

Based on the table, we were able to get a 92% classification of objects.

In this activity, I give myself a grade of 10 for understanding ang implementing properly LDA.

I thank Mark Jayson and Luis for happy conversations while doing the activity in the IPL lab.

References:
[1]http://people.revoledu.com/kardi/tutorial/LDA/Numerical%20Example.html

Activity 14 - Pattern recognition

2009-08-31T06:46:00.000-07:00

The human brain has the capacity to discriminate from one object to another. From the stigma it receives from the eyes (-or from other sensory organs), one of the things that human does is to find significant features unique to the objects and draw,as well, overlapping features between different objects that will allow to recognize things. The aim of computer vision is to approximate the pattern recognition capabilities of the human brain. Various features can be found in a certain object or pattern -- color, size, shape, area, etc. We use such features to identify an object as a member of a set that has the same properties (this set is called a class).
In this activity, we are tasked to classify objects according to their classes using their physical features. We should be able to (1) find the correct features that will effectively separate the classes and (2) create a classifier (or an algorithm) that is suitable for the task.

In Task 1, we determine the useful features and train a representative set of objects of known class and calculate the feature values (for example, area). An average value per feature is taken from the whole class.

where m_j is the average for all feature x_j of class w_j. It is a 1xn matrix for all features n. These average values will be used as the basis for classification. In our case, we used the r-g chromaticity (discussed in Activity 12)and dimension ratio (length/width) as features for the set of classes we gathered (image is shown below).

Figure 1. Dataset of classes used that contains a set of short and long leaves, flowers, and 25-cent and 1-peso coins.

Once the m_j for all the classes are done, we are now ready to classify a test object based on its features. We extract all the needed features x and determine what class it belongs to. For Task 2, we used the Euclidean distance as a classifier for our test objects. It is given by:

d_j(x) = x^Tm_j - 1/2 m_j^Tm_j

This determines how close the features x of an object to the features of class j. The object will be classified as belonging to a class which yields the largest d_j.

Results

Below is 3D plot of features extracted from the objects in our dataset. We readily see an obvious separation of class clusters of the leaves (short and long) and of the flowers. One reason is that the color of the flower is different from the color of the leaves. The difference in the dimension ratio of the short leaves and long leaves has also separated the two leaf types. However, class clusters of the coins were not well separated. Based on the criteria used, the dimension ratio cannot possibly separate the two since both are circular (dimention ratio ~ 1). Only the color determined by the r-g chromaticities will be able to separate them.

Figure 2. Separation of class clusters using the r, g and dimension ratio as features

We take a look at the distribution of points using only the r-g chromaticities -- Fig.1 as viewed from the above. Here, we see that there is an overlap between the clusters of 1-peso and 25-cent coins. In r-g chromaticity diagram (see Fig. 3 in Activity 12), the color yellow and white, which are the colors of the 25-cents and 1-peso, respectively, lie close to each other. It is thus expected that the clusters be consequently close. The overlap is caused by the highly reflective surface of the coins. In our case, a 25-cent coin is misidentified as a 1-peso coin. Highly reflective 25-cents (the new ones) will appear as though it is a 1-peso coin.

Figure 3. r-g chromaticities of the samples. A single data point overlap of the between two clouds of data points -- cloud of peso coins (inverted triangle) and 25 centavo coins (filled diamonds) is seen indicating that the r-g chromaticity wrongly classifies a 25-cent as a 1-peso coin.

From the plots, we can already observe whether the features we used is capable of classifying objects. However, this can only be handy when dealing with at most 3 features. More than this will be very hard to visualize. The advantage of a classifier such as the Euclidean distance is that it essentially decomposes all the features into one number, thereby resulting to easier classification of objects. Below is a table of test objects and there classification using the algorithm.

Table 1. List of Euclidean distances of the samples from the training set. Sample x - i are samples taken from a set of samples i where i = {peso, 25 cents, long leaf, short leaf, flower}.

Values highlighted with yellow belong to the column where the samples are classified. For example, Sample 1-1 resulted the highest value at d_peso hence classified as a 1-peso coin. All the classifications were correct except for one sample. This is a 25-cent coin wrongly classified as a 1-peso coin. The features used are 96% accurate in classification.
Note that the features used are not enough to completely separate the classes of coins. One can use the area as a feature to be able to do such.

In this activity, I give myself a grade of 10 for doing the job well.

I thank Luis, Mark Jayson and Irene for useful discussions. Especial thanks to Mark Jayson, Irene and Kaye for lending the dataset.

Activity 13 - Correcting Geometric Distortion

2009-08-31T04:15:00.000-07:00

Digital images captured by cameras, most of the time, are different from the one we see with our eyes -- although some are unrecongnizable. The color for example is altered by the spectral sensitivity and white balance settings of the camera. Another factor that contributes to these differences is the distortion it adds to the captured image due to its lens configuration. Pincushion and barrel distortions are observed in spherical lenses. In pincushion distortion, the image seems pinched in the middle and expanded in the boundaries. In barrel distortion, the image seems bloated in the middle and pinched at the sides *[1]. Below are images that demonstrates the two distortions.

Figure 1. (left) Pincushion and (right) barrel distortion of a square grid

Cameras using compound lenses have considerably minimized these effects. Single lens cameras render a much more obvious distortion.
In this activity, we are tasked to correct distorted image. We download an image from the net with distorted grid patterns. The idea is to perform an image processing technique to make straight lines look straight.

Steps in correcting a distorted image:
1. Find the image coordinates the four corners of a square where the effect of the distortion is small. Usually, this can be found at the center of the image or at the optical axis of the camera.

2. Create an ideal grid from the obtained coordinates. (The grid only contains corner points.)

3. Calculate c's from the corner points using the equation below to be able to transfer distorted image coordinate to the ideal image coordinates.(c's are unique to a specific square)

(x,y) is from the ideal image coordinates, (x,y) is from the distorted image coordinates. The c's allow us to determine the location of a pixel from the ideal image to the distorted image.

4. Transfer grayscale values.

Nearest neighbor technique: In determining the distorted image coordinates that correspond to the ideal image coordinates (using the c's), cases may occur that the calculations would yield non-integer coordinates. In the simplistic case, one can round the coordinates, then the grayscale value of the that pixel will be transferred to the ideal image. The advantage of this is that it is fast, although the reconstruction is pixelated and jagged.

Bilinear interpolation: If the location of the pixel from the ideal image corresponds to integer distorted image coordinates, pixel values from the distorted image can be directly transferred to the ideal image. If not, we use the equation below:

a, b, c and d can be calculated using the surrounding pixels of the pixel with non-integer coordinates.

Results

In this work, we demonstrate the two techniques used in the correcting distorted images as well as the other application of the techniques.

Figure 2. Distorted grid (above) and the corrected image (below)

The above figure consist of the distorted grid and the corrected grid. In this particular example, the nearest neighbor technique was used. We can see that this technique can fairly reconstruct the image. The 'jagged' effect of the technique is less manifested because of the thick lines of the grids.

Note: In this work, the distorted image coordinates of the corner points were determined using the locate function in Scilab. Given a very large grid, this procedure is very tedious and time-consuming. Alternatively, one can perform image processing techniques to determine the coordinates. However, when given a real distorted image as in Fig.3, locate is much more convenient to use since it would be harder to separate the grid lines.

The procedure can also be applied to distort an image. Here's how it's done.

We have a very good, undistorted image. We overlay a distorted grid to the image. Then we perform the same procedure as done previously. In this particular case, we use the bother technique, bilinear interpolation to reconstruct the image.Below is the image and its distortion.

Figure 3. (above) image and (below) its distortion

Here, we observe the success of the procedure. We have straightened the straight lines and were able to distort the image. Image distortion will be highly appreciated if the pixels from the ROI is contained in obviously distorted grids.

I give myself a grade of 10 fro exploring to the two techniques and the possible application of the procedure.

I thank Mark Jayson, Luis and Orly for useful discussions.

References:
[1] Activity 13 manual
[2]http://upload.wikimedia.org/wikipedia/en/9/97/Pincushion-distortion.jpg
[3]http://upload.wikimedia.org/wikipedia/commons/thumb/6/63/Barrel_distortion.svg/600px-Barrel_distortion.svg.png
*Sentences are stated vervatim.

Activity 12 - Color Image Segmentation

2009-08-05T11:46:00.001-07:00

From the previous activities, we find that thresholding is indeed a powerful way in separating our region of interest with the rest of the image. This technique, however, has its limitations. Thresholding can only be applied to images with single-value pixels such as grayscale images. When dealing with colored images (where each pixel has 3 values -- RGB), one has to reduce it first to grayscale before thresholding. Upon reduction, different colors may appear the same in grayscale. As such, separation of the ROI may not be possible just by thresholding. In this activity, we perform another technique for ROI separation in colored images -- color image segmentation.

Colored images may have different shades (or brightness) of colors inherent to it. Thus, it is better represent the image into brightness and chromaticity information. This is called the normalized chromaticity coordinates (NCC).

Transforming RGB to NCC
Per pixel, we add up the RGB values I = R+G+B. The NCC is then:

r = R/I
g = G/I
b = B/I

Note that b is also equal to 1-r-g hence only the r and g values are needed to represent the chromaticities r,g and b. We essentially have transformed RGB to rgI values where I accounts for the brightness of the pixel.

Two techniques will be used in this activity - parametric and non-parametric segmentation.

Parametric segmentation

In parametric segmentaion, we crop a portion of the ROI and take the mean and standard deviation for the chromaticities r and g. To check whether a pixel belongs to the ROI, we check its probability of belonging to the ROI. Since pixel has two chromaticities, we calculate the effective probability

p = p(r)*p(g)

where p(r) and p(g) are the probabilities of a pixel with chromaticities r and g, respectively belongs to the ROI. We now assume a Gaussian function using the calculated mean and standard deviation for each chromaticity such that the probabilities p(r) and p(g) takes the form:

where x is the r-g chromaticity of the pixel.

The result of the whole process is a martix of p values the same size as the image. High values of p will only occur at ROI and thus we effectively separated the ROI from the rest of the image. Below are examples of the results using this method.

Fig.1 (upper left) the raw image and (rest) the segmented images of red, yellow and blue patches.

Fig.2. (left) the raw image ang (right) the segmented image using parametric segmentation

Here we see a complete segmentation of the ROI.

Non-parametric segmentation

This technique does not assume any form of function unlike the previous one. What is need is just the 2D chromaticity histogram of a portion of the ROI. (Note that the axes are the r and g chromaticities in different levels). In creating the histogram, binning is important because it will determine the quality of the segmentation. If bin sizes are small the details of the ROI is preserved but the toleration of the technique in terms of chromaticity would be very low which results to dark regions in the ROI. If on the other hand, the bin sizes are large, the toleration becomes high such that most of the pixels in the ROI is bright but we compromise its details. Our choice of bins depends on what we want. In this particular case, we used a bin size equal to 256, 32 and 10 to demonstrate their reconstuction differences. Once the histogram is created, we now back-project the chromaticity values of each pixel in the image to the histogram. Steps of backprojection is shown below.

Figure 3. Verification of the histogram

Before the backprojection, we first verify if the peak of the histogram corresponds to its color in the chromaticity space. The patch used is yellow. We observe that the peak coincides to the yellow region.

Histogram backprojection
- Obtain the chromaticity (r,g) values of each pixel in the image.
- For each pixel, find its position (r,g) on the histogram and find its value.
- The value at the pixel location is changed to the histogram value found.

In this technique, we can infer that high values will be obtained when the (r,g) values of the pixel in the image corresponds to the (r,g) where high values occurred in the histogram. Shown in Fig.3 is the result of this technique.

Figure 4. (clockwise from upperleft): the raw image and the sementated images using (1) 256 bins, (2) 10 bins and (3) 32 bins in the histogram. (Note that large bins correspond to small bin sizes).

In Fig.4, we indeed see that the segmentation quality differs at different bin sizes. This technique is much better than the previous one since it does not any assumptions. The only challenge is that one should be able to find the proper bin sizes depending on what he wants.

In this activity, I give myself a grade of 9/10 for doing a good job.

I thank Ma'am Jing for helping me get through some problems encountered in this activity.

Activity 11 - Color Image Processing

2009-08-05T11:39:00.000-07:00

Objects have their innate color reflectance. The amount and color of light, however, does not depend on reflectance alone. Light reflected also depends on the spectrum of the light source that illuminates the object. Moreover, the device used to capture effective reflectance also comes into play. Cameras are one of these devices. Cameras have their own spectral sensitivities like the eye and thus render colors they detect depending on such sensitivities.

Images are rendered in different modes -one of which is in RGB mode represented in the following equations.

ρ(λ) and S(λ) are the reflectance and light source, respectively and, η(λ) are the spectral sensitivities of the camera (subscripts indicate the channel sensitivity - red, green and blue). K_i are called the white balancing constants of the form:

These constants are important so that the camera renders the correct reflectance of an object. Digital cameras have this option to change their white balance depending on the current source that illuminates on the object i.e. cloudy, sunny, fluorescent, etc. When the white balance of the camera does not match the light source, what happens is that the reflectance of the object is rendered incorrectly by the camera. In such case, we say that the camera is improperly white balanced. We can readily see this if, for example, we see a known white object as a non-white (bluish probably). When faced with such images, image processing techniques has to be done to correct them.

In this activity, we explore two of the white balancing algorithms used for such images -- the White Patch algorithm and the Gray World algorithm.

Algorithms
White Patch (WP) algorithm is implemented by first finding a patch of a known "white" object in the image. The average RGB values (Rw Gw and Bw) of the entire patch is calculated. All 3 channels in the image is normalized using those average values. Essentially, we force the "white patch" to be white. As such, all the other pixels are affected also.
Gray World (GW) algorithm, on the other hand, has the same implementation as the WP algorithm only that the image is normalized by the average RGB values of the entire image. The algorithm , given an image with a sufficient amount of color variations [1], assumes that the average color of the world is gray - a family of white. In this case, we force the image to have a common gray value. As a result, the effect of the illuminant on the image is removed. The image then becomes white-balanced.

Below are images of a scene under constant lighting environment captured by a Motorola W156 camera phone with different white balancing (WB) settings. Both algorithms were performed in each image.

Fig1. (left column)Images containing different colors under different WB settings of the camera and, white-balanced images using (middle) WP and (right) GW algorithms.

Comparing the raw images from the reconstructed images using either of the algorithms, we see that we have successfully implemented the techniques. The background (which is in actual white) in the raw images initially does not appear white. For example, the background of the image using indoor WB camera settings appears blue but now appears white upon white balancing.
Also, we observe that the both algorithms produces different results - which is expected since they have different implementation. It can be seen that the images produces using GW algorithm is much brighter than that of the WP algorithm. In this case, we can infer that after the GW implementation, a lot of pixel values saturated which means that these pixel values exceed the supposedly maximum value of 1. This is not what we want. Here we see the disadvantages of the the GW algorithm. One is that, the image must satisfy a condition (that it must contain a sufficient amount of color variations) for the GW algorithm to work. Another is that since the algorithm normalizes the image using the average RGB values of the image it is expected that saturation of RGB values will occur.

In this part we now use the two algorithms onto monochromatic images.

Figure 2. Monochromatic images capture using different WB camera settings and the reconstructed images using the two algorithms.

From the figure above, we see that the reconstruction of the images using the WP algorithm is way better than the GW algorithm. Again, we go back to the first condition an image must satisfy for the GW to work. Since the images are monochromatic; that is, they do not contain a variety of colors, GW is not appropriate to use.

Strengths and Weaknesses of Each Algorithm.

From the results obtained, we can pretty much gauge the strengths and weaknesses of each algorithm. The GW algorithm does not need any white patch to reconstruct the image hence it is useful when one does not know anything about the image.
Based on the results, we can infer that WP algorithm can work whether the images are monochromatic or colorful. And algorithm can be generically used for white balancing images. However, this will only work if one knows a "white" object in the image.

In general, GW algorithm is more convenient to use than WP. It is easier to implement since we do not need to find white patches in the image and very handy when dealing with colorful images.

I give myself a grade of 9/10 in this activity since I was able to explain the strengths and limitations of the algorithms.

I thank Orly, Thirdy and Master for useful discussions.
References:
[1]http://scien.stanford.edu/class/psych221/projects/00/trek/GWimages.html

Activity 10: Pre-processing Text

2009-08-01T08:24:00.000-07:00

In an image, most often than not we only need a certain part or a certain element -- the region of interest(ROI). It may be an object or a text. Either of the two, we need to perform image processing techniques to extract what we need. In this activity, we exploit these techniques particularly in text retrieval. One of its applications is in handwriting recognition.
In this case, we used a scanned image with texts shown in Fig.1(left) as our raw image.

Figure 1. Raw image used in pre-processing text.

We cropped a portion from Fig.1.(left) to be cleaned. Initially the image is rotated so we have to rotate it back so that image processing will be easier. This was done by taking the Fourier transform(FT) of the cropped image in Fig.1(upper right). The frequencies of rotated horizontal lines manifest as a tilted vertical line in the Fourier space. The angle of the tilt was measured using MSPaint. This angle was used to rotate the image back using mogrify command in Scilab. Alternatively, one can measure the tilt angle using Fig.1(upper left) immediately without performing FT.

Figure 2. Removal of horizontal lines on the selected portion of Fig.1(left).

Now that the image is rotated back, the next step in cleaning the image is to remove the horizontal lines. To do this, we get the FT of the image then remove the frequencies of the horizontal lines and take FT again to obtain the reconstructed image.

Removal of frequencies: Horizontal lines as explained earlier will manifest its frequencies along the vertical axis of the Fourier domain. Blocking out these frequencies will effectively remove the horizontal lines. We do this by multiplying the FT (Fig.2 upper right) to a mask (Fig. 2 lower left) that will completely block out the frequencies along the vertical axis (excluding the center so as not darken the image).

After the multiplication, we now take its FT to obtain a partially cleaned image (without the horizontal lines). The image is shown in Fig.2(lower right). Using this image, we are now ready to separate the region of interest (the text) and the background.

We use a simple technique for the separation process -- thresholding. Below is a comparison between Fig.2(lower right) and the thresholded image.

Figure 3. (left) the pre-cleaned image and (right) a thresholded image of the left.

Here we see a separation of the text and the background. Notice, however, we observe a black line dashing through most of the text due to the horizontal line removal previously done. What we need to do next is to remove this black lines. We then preform a closing of the image using a vertical structuring element. In this case, since the texts in the image are tilted, a 4x2 diagonal was used as a structuring element.

Figure 4. The structuring element used and the resulting image after closing

Here we see that that most of the black lines were removed. This is very promising to see. We have used a proper structuring element for the morphological operation. Also, we see that some letters are distinguishable. In Fig.4, we have labeled the blobs found in the image (indicated by a different gray level) and we indeed see some single-letter blobs. It is also interesting to note that the word'"cable" (encircled) is fairly readable although some words are not like the word "power" which we can blame to the quality of the handwriting i.e. some pen strokes are faint that it disappeared upon thresholding.

One of the problems encountered is that it is hard to reduce the letters into 1 pixel thick. To do this, a erosion process must be done using a 1x2 structuring element. However, some strokes on the letters are thin enough such that they are removed upon erosion which we do not want. Hence this operation was not done. Also, it can be observed that the previous "DESCRIPTION" is now almost unreadable. This has been the trade-off of the operation. Since the "DESCRIPTION" did not have the same problem as the rest of the text (dashed by black lines), performing of the close operation made it worse.

Now, we explore if we can find instances of the word "DESCRIPTION" from the whole image. We do this using the following steps:

Step1: Rotate the image
Step2: Binarize the image
Step3: Create a "DESCRIPTION" pattern with the same font and size as in Fig.1(bottom right).
Step4: Correlate the FT of the binarized image with the FT of the produced pattern.
Step5: FT the result in Step 4.

Below are images as a result of Steps 1-3.

Figure 5. (left) binarized rotated image of Fig.1(left); and (right) the pattern used for word
detection. The pattern has font italized Arial Black with fontsize 9.

We then proceed with Step 4 and 5. Patterns that are the same with Fig.5(right) will have a produce a peak at their location - an indication of high value correlation. Below is the image as a result of Step 4 and 5.

Figure 6. Image as a result of correlation

We observed three peaks in the image which correspond to the word "DESCRIPTION". This can be readily verified in Fig.1 (left).

In this activity, I give myself of 10 for doing the job well.

Activity 9: Binary Operations

2009-07-29T08:17:00.000-07:00

In image processing, especially if one is dealing with a region of interest (ROI) in an image, it is necessary to separate the ROI from the background. Thresholding is one way to effectively separate the two. This will reduce the image into a binary one. However, this requires an optimum threshold but this optimum threshold does not guarantee us a perfect separation of the ROI and the background. Instances that would sacrifice some of the ROI to zero out the background and vise-versa may be encountered. Hence, a threshold somewhere in between must be chosen that would leave some white pixels in the background and holes on the ROI.

In such cases (which happens more often), morphological operations must be performed for two basic reasons: (1) to remove the small white pixels on the background and holes on the ROI and (2) to separate touching blobs in the image. Blob separation is very important for accurate determination of cell sizes, etc. Here,we use the closing and opening operators to serve our purpose. Closing is just a dilation of the image followed by an erosion. Opening is the opposite.

Using the image below, we segment this into 256x256 sub-images. Threshold is applied to binarize each sub-image. Morphological operations is then performed per binarized sub-image. The area of each blob (in pixels) is calculated per sub-image. A histogram of blob areas is then obtained from all the sub-images. The thresholded is shown in Fig.2.

Figure 1. Grayscale image of ''cells''

Figure2. Thresholded sub-images in Fig.1

Before the process of opening and closing, choice of the structuring element (SE) is needed. In choosing, we have to take into account that the SE should not be large enough that upon dilation will not cause originally separated blobs to connect. As such, performing an erosion operation will not be able to separate them again. If, on the other hand, the SE is small enough, closure of relatively large holes will not be effective. Even repetitive closing would be futile. A balance between the two must be chosen.
Here, a 3x3 cross was chosen as the SE, since there are no large holes in the binary images obtained. One of the advantage of this SE is that this will not join blobs upon closing thus we can accurately determine the cell sizes.

Figure 3. Structuring element used for morphological operations

Using thus SE, we perform three opening to effectively remove the white pixels in the background and one closing to close the holes in the ROI. Then we use the bwlabel to label the blobs in each sub-image and calculate their areas in pixels. A figure below shows all the blobs obtained from all sub-images. Each blob are shaded differently.

Figure 4. Cleaned images from Fig.2.

From the areas computed, we plot the its histogram to approximate the sizes of the cells. In Fig.4, we can see that most blobs are composed only of one cell. Judging from this, we can pretty much predict that the peak of the histogram would occur at area of single-cell blobs. Below is the histogram obtained.

Figure 5. Histogram of areas in pixels

The peak corresponds to the approximate size of the cells. Examining Fig.5, we can infer that the cell size is approximately 550 pixels.

In this activity, I give myself a grade of 10 for doing the job well. I acknowledge Orly, Mark Jayson and Thirdy for useful discussions and laughter.

Activity 8 - Morphological Operations

2009-07-22T20:57:00.000-07:00

In this activity, we were tasked to study dilation and erosion -- two of the morpholigical operation techniques applied in image processing. These operations (usually done in binary images) are used to enhance the image for better extraction of information. Morphological techniques have different effects on the image.
Let's take a look a dilation first. Using a structuring element; the dilation of an image is all the translation of the center of the element in the image such that their intersection is non-zero. Dilation, in effect, enlarges the image. The figure below demonstrates dilation of various images using different structuring elements.

Figure 1. Dilation using different structuring elements

Dilation of images can be done by hand. We just have to know the structuring element that would reshape the image. Before the simulation was done, the images as a result of dilation were first predicted. All my predictions were correct in comparison to the images obtained in Fig. 1.

Erosion, on the other hand, shrinks the image. One of the application of erosion is to remove unnecessary pixel on the image. Erosion takes all the translations of a structuring element such that all its pixels are on the image. Erosion can also be predicted by hand. The figure below demonstrates erosion of different images using different structuring elements.

Figure 2. Erosion of images using different structuring elements.

All my predictions were correct compared to images above.

In this activity, I give myself a grade of 10 for predicting properly the images as a result of dilation and erosion and also for properly creating the images above.

I thank Orly Tarun for uselful help.

Activity 7: Enhancement in the Frequency Domain

2009-07-20T07:48:00.000-07:00

7A. Convolution Theorem

In this part, we were tasked to observe the changes in the modulus of the FT of the images as the size of the apertures are increased.

Figure1. 2-dot images and their FTs

Using two one-pixel dots, Fig.1(upper left), a grating-like FT was observed. As the radii of the dots increase, the FT becomes circle with decreasing radius located at the center of the two dots. this is expected since the FT is in inverse space.

Figure 2. square images and their FTs

Same was observed using squares. As the width of the squares increases, the width of their FTs decreases. Again, this is expected.

Also, for Gaussian dots, same was observed in Fig.1. As the variance of the Gaussian is increased - which results to an increase on the width of the Gaussian dots, their corresponding FTs is decreased.

Figure3. Guassian dot images at different variances and their FTs

We also take a look on how the FT of these Gaussian dot images differ from that of the inverted images (max becomes min; min becomes max). Here we plot in Fig.4 the images of the real part and imaginary part of the original and inverse Gaussian dots for comparison.

Figure 4. (upper row) images for the real part of original and inverse of 2 Gaussian dots; (lower row) images for the imaginary part of the original and inverse of 2 Gaussian dots.

We observe from the figure that only the real part changed upon inversion. The imaginary part remained the same.

7B. Fingerprints: Ridge Enhancements

In this part of the activity, we are to remove the blotches present in the image. Here we use the FT, we use a mask to remove the frequencies that result to the blotches and retain the information needed. We then take the inverse FT, we obtain a better image.

Figure 5. Removal of blotches present in the image in top left. A mask in bottom left was multiplied to the FT of the image (top right) to obtain a reconstructed image in bottom right.

As observed in the reconstructed image in Fig. 5 (bottom right), we essentially remove the blotches from the original image. Although the area surrounding the blotches darkened. This can be further improved if a better mask is used.

7C. Lunar Landing scanned Picture: Line removal

In this part of the activity, we were tasked to remove the vertical lines present in image below. Here, we use the Fourier transform to take information of the image in the frequency domain. We remove the frequencies that cause the vertical lines. From the previous exercise, we have learned that the frequencies of vertical lines are contained on the horizontal axis of the FT. Therefore, removal of frequencies lying on the horizontal axis removes the vertical lines.

Figure 6. Steps done (from upper left to lower right) to remove the vertical lines in the image.

From the FT in Fig.6(upper right), a mask was created Fig.6lower left) that will remove certain frequencies on the FT. The mask was then multiplied to the FT. Its FT was taken to recover the better image. In Fig.5, we see that the vertical lines were removed by removing the corresponding frequencies. Using FT, we were able to obtain a better image.

7D. Canvas Weave Modeling and Removal

In this part, we use again the technique demonstrated in 7C to remove the weave cross-hatches on the image in Fig.6(upper left). Since, we are dealing with cross-hatches (intertwined vertical and horizontal lines), we must take the frequencies on the horizontal and vertical axis of the FT domain. We might also consider that the cross-hatches are a result of the product of sinusoids (which has also been observed in Activity 3 Fig.5) and hence we also remove that constitute these sinusoids.

Figure 7. Steps of weave cross-hatches removal

Using the mask in Fig.7(lower left), a reconstructed image in Fig.7(lower right) was obtained. Here we see a decrease of cross-hatch visibility.

Figure 8. Weave pattern extracted from the image.

Furthermore, we can see in Fig. 8 a weave pattern extracted from the image. This was done by inverting the mask and multiplying it to the FT of the image and taking the inverse FT.

Although the technique produced promising results, our limit on reconstructing the image (or in other words our limit in creating better reconstructions) is dependent on the visibility of the peaks seen in the frequency domain. We can only block the frequencies that we see.

In this activity, I give myself a grade of 10 for doing a good job.

I thank Cindy Esporlas very much for helping me to create a Gaussian dot. I thank also Master and Thirdy for their help.

Activity 6 - Properties of the 2D Fourier Transform

2009-07-08T19:56:00.000-07:00

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

In this part of the activity, Fourier transforms (FT) of different apertures were investigated. Here, we familiarize on the patterns formed by this images upon applying FT.

Figure 1. (lower row) FTs of the images on the upper row.

The figure above shows the patterns obtained as a result of taking the FT of images of different apertures (Fig.1 upper row). We see here interesting patterns. For example, the FT of a square is a sinc function along the x and y-axis. Many other shapes can be applied with FT and may observe other patterns as well.

6.B. Anamorphic property of the Fourier transform

Figure 2. (below) FT of (above)sinusoid images of different frequencies (left) 1Hz (middle) 4 Hz and (right) 10 Hz)

The two points in the FT of the images are the frequencies (f and -f) of the sinusoid used (0 Hz is located between the two points). As the frequency of the sinusoid increases, the distance between the two points also increases. The points get farther away from the center as the frequency is increased.
Usual images don't have negative pixel values -- in other words the image is biased. The problem now arises, since the image is biased, is how to determine correctly the frequencies contained in the image. This is one of the powerful uses of Fourier transformation.

Figure 3. FT of 10 Hz-sinusoid images with (upper left)constant bias and (upper right)sinusoidal bias with 1 Hz frequency

The figure above shows the a biased sinusoid image with their FTs. Fig.3(upper left) is constant-biased and the Fig.3(upper right) is sinusoid-biased. As seen, the FT of the constant-biased image has three points -- the peripheral points are the frequencies (f and -f)of the image and the point at the center (at 0 Hz) accounts for the bias. The frequency of the image can now be correctly determined ignoring the bias.
What if the image has a sinusoidal bias? These biases include noise, etc. Usually, in signals and images, non-constant biases are of low frequencies. If the frequencies contained in the image is sufficiently high, we can separate easily separate the low from the high frequencies using the FT. An example is demonstrated in Fig. 3 (right). Fig.3(upper right) is a 10 Hz-sinusoid with a 1 Hz-sinusoid bias. using FT, 2 distinct frequencies were observed: (1) the one closer to the center and (2) farther away the center. The frequency closer to the center (meaning of low frequency) accounts for the bias and the farther accounts for the frequency of the image. The low-frequency can be disregard and we are able to find the frequencies contained in the image.
The obtained FT can also be multiplied with a function that to remove low frequencies. We can take its inverse FT to retrieve a better image. Note, however, that this technique cannot always serve its purpose right. Once the frequency contained in the image is close to the frequency of the bias, doing such might destroy the image.

In this part, we rotate the image the image and look into how the FTs change with the rotation.

Figure 4. (below) FTs of the (above) rotated images at : (left) 30^o (mid) 40^o and (right) 50^o

We can observe in Fig. 2 and 3, that the axis of the dots of the FT is perpendicular to the bands of the image. This observation is the same even when the image is rotated. The FTs rotate as the image is rotated. This is another property of the Fourier transformation. The axis of the FTs are just rotated and same frequencies were obtained.

Images may contain different frequencies (assuming the image has no biases). It can be a combination of sinusoids of with these frequencies (i.e summation or multiplication of sinusoids). Below are images of sinusoidal combinations and their FTs.

Figure 5. (upper row) Images formed from sinusoidal combinations and (lower row) their FTs. (upper left) z = sin(2*pi*f*X)+ sin(2*pi*f*Y) ; (upper right) z = sin(2*pi*f*X)sin(2*pi*f*Y)

The FT of the sum of 2 sinusoids (Fig. 5(lower left)) is just the sum of the FT of each sinusoid:

F_t =F{sin(2*pi*f*X)} + F{sin(2*pi*f*Y)}

This is another property of FT -- linearity. The FT of sum is the sum of FT.

This property holds regardless of the number of additions to the orignal pattern. This is illustrated in the figure below.
For the product combination,

z = sin(aX)sin(bY),

the expression can be expanded into sum of cosines given below:

z = [cos(aX - bY) - cos(aX + bY)]/2 or

z = [sin(aX - bY + pi/2) - sin(aX + bY + pi/2)]/2

This is just sum of rotated sines at opposite directions (due difference in signs of bY). This results to a pattern in Fig.5(lower right).

Figure 6. Product combination of 2 sinusoids as in Fig. 5(upper right) with additional two rotated sinusoids with frequencies 10 Hz and 25 Hz. The 10 Hz signal is rotated at 30^o. The 25 Hz-sinusoid is rotated from 30^o, 60^o and 90^o.

Using the pattern in Fig.5 (upper right) plus additional rotated patterns, since FT is a linear transform, an expected pattern is:

F_overall = F{pattern1} + F{pattern2} + F{pattern3}

Indeed, the prediction is right. The FT of the product combination is retained, the FT of the rotated 10 Hz sinusoid (as in Fig.4(lower left)), and the FT of the 25 Hz sinusoid at different rotation angles.

In this activity, I give myself a grade of 10.

I thank Jica Monsanto, Miguel Sison, Orly Tarun, and Jayson Villangca for useful discussions. Special thanks to Sison for making me realize that some of my images were wrong.

Activity 5 - Fourier Transform Model of Image Formation

2009-07-06T19:40:00.000-07:00

Activity 5.A - Familiarization with discrete FFT

In this part of the activity, we are to familiarize how the FT works. We apply the FFT algorithm to images. We also the effect of fftshift() in the FT of the image and the results of double FFT.

Figure 1. (middle)Fourier transform and (right) double Fourier transform of the (left)image.

In 3D space, the circle (left) can be imagined as a cylinder. It can further be imagined as rectangles rotated about its center. A rectangular signal has a sinc transform. Now, a sinc rotated at its center produces an airy disk as observed in Fig.1 (middle). Also, the diameter of the center of the airy disk is small which is the inverse of the diameter of the circle which is a property of FT - changing to inverse space.
Fftshift() was already used on Fig.1(middle). Without using fftshift(), the FT looks like below:

Figure 2. FT (left)with and (wihtout) using fftshift

As observed, the center of the airy disc are separated to the four corners of the image. This is an effect of the FT. The diagonal quadrants of the FT are interchanged. fftshift() interchanges the quadrants back resulting to Fig. 1(right).

Using a different image, as in Fig. 3(left), images alongside with it were obtained.

Figure 3. (middle) FT and (right) double FT of the (left) image.

Notice that Fig.3(right) is an inverted image of Fig.3(left). The same happened in Fig.1 - although not obvious since a circle is symmetric.
This should be expected since:

F{F{f(x)}} = f(-x)

leading to an inversion of the image.

Activity 5.B - Simulation of an Imaging Device

In this part of the activity, images formed in an imaging system were simulated. Particularly, we have an object f (Fig.4) and an aperture g in Fig.5(upper) and simulate how image produced looks like using Convolution theorem. The algorithm is that the FT of f and g are multiplied and then inverse FT is done. This results to an image that looks like f and g. This is the image formed by g.

Figure 4. Object used in the simulation

Figure 5. Images formed using different sizes of circular apertures.

It can be observed from the figure that the readability of the image decreases (the image becomes blurry).

Activity 5C. Template matching using correlation

Correlation is a measure of how the a certain function is similar to another function, how similar an image to another image when talking about 2D arrays. One of its application is in pattern detection. Correlation values of highly similar functions is high as well. In this activity, we were tasked to find the letter A's contained in a text "The rain in Spain stays mainly in the plain."

Correlation is done using the equation below:

p = f⊙g = f(-x,-y)* g(x,y)

where (*) represents the complex conjugate of the function.

Figure 6. (left) the text used; (middle) the letter to be detected and (right) the image as a result of correlation

Upon applying correlation, Fig.6(right) was obtained. Bright spots on the image are high correlation values indicating an "A" pattern.

Activity 5D. Edge detection using convolution integral

In this part, using the image in Fig.4, we again exploit convolution to detect edges in the image. The idea is to convolve a desired pattern to find with the image being edge-detected. In other words, we find the certain pattern on the image.

Figure 7. Images upon application of edge-detection using (upper left) horizontal pattern; (upper right) vertical pattern; (lower left) diagonal pattern; and (lower right) a point.

We can see from the images that image detection algorithm allowed us to detect patterns on the image teh same as the desired pattern used(i.e, horizontal lines appear on Fig.7(upper right) since the image was convolved with a horizontal pattern). it is also interesting to note that using a point pattern allows us to detect all edges of the image.

In this activity I give myself a grade of 10 for doing the activity right. =)

Activity 4: Enhancement by Histogram Manipulation

2009-06-30T06:38:00.000-07:00

In this activity, we were tasked to get a poor-contrast grayscale image from the web. Using the histogram manipulation technique, the goal is to make a high-contrast image out of the original image.
In histogram manipulation, the first step is to take a histogram of the pixel values in the image (histogram is plotted with pixel values on the x-axis and frequency of occurrence of the pixel value on the y-axis. The histogram gives an idea of how many of a pixel value is contained in the image. Once the histogram is made, a cumulative distribution function (CDF) can be plotted. CDF is just the cumulative sum of the values on the y-axis plotted with pixel values. (It is better to work with CDF than histogram as we would see later.)

To start with, below is the graysacle image taken from the web.

Figure 1. Grayscale image taken from http://homepages.inf.ed.ac.uk/rbf/HIPR2/images/wom1.gif. The figure is then converted to a jpeg file.

Using Scilab, the histogram, as well as the CDF as shown below, of the image in Fig.1 was obtained.

Figure 2. (above) Histogram and (below) CDF of the image in Fig. 1.

As observed from the plots in Fig.2, the CDF is cleaner to see than the histogram. Having the advantage, CDF is easier to manipulate than the histogram.
In histogram manipulation, a desired CDF is needed. This will be used to alter the pixel values contained in the original image. The method is illustrated below.

Figure 3. Data manipulation procedure

The desired CDF is of a uniform distribution.

Data manipulation procedure done in the following:
(Step 1) Take a pixel value from the CDF of the image in Fig.1
(Step 2) Obtain the value on the y-axis corresponding to the pixel value
(Step 3) Project the y-axis value onto the desired CDF
(Step 4) Take the pixel value from the desired CDF corresponding to the y-axis value.
The pixel value in the image is then replaced by the pixel value from the desired CDF.

The code implementation for a linear desired CDF is written below:

h1 = zeros(256,256);
for i=0:255
new = m(i+1)*255;
h1 = h1 + (im == i)*new;

Performing this procedure for all pixel values (0-255), a reconstructed image shown below was obtained.

Figure 4. Reconstructed image

As seen, the image highly contrasted compared to the previous image as a result of the procedure. We also characterize the image by taking its pixel histogram and CDF.

Figure 5. (above) Histogram and (below) CDF of the reconstructed image.

The histogram of the reconstructed image is different from the original image. It is also interesting to note that the spread of the histogram of the reconstructed image is larger than that of the original image indicative of the high contrast of the reconstructed image. Moreover, looking at the CDF of the reconstructed image, it resembles a straight line closer to the desired CDF.
To further explore the technique, a desired CDF that mimics the nonlinear response of the human eye was used. Below is the desired CDF used in reconstructing the image.

Figure 6. Exponential CDF used in reconstructing the image

The same procedure was done as before and a reconstructed image shown below was obtained.

The code implementation for a nonlinear desired CDF is written below:

h1 = zeros(256,256);
for i=0:255
new = int(50*log(m(i+1)*(exp(255/50)-1)+1));
h1 = h1 + (im == i)*new;

Figure 7. Reconstructed image using a nonlinear desired CDF

Compared to the previous reconstruction, the image in Fig.1 is brighter (many of the pixel have high values). we can readily attribute this to the desired CDF used. The latter favors high-valued pixels. As a result, a bright reconstructed image was obtained.

Figure 8. (above) Histogram and (below) CDF of the reconstructed image using nonlinear desired CDF

It can also be seen from Fig. 8 that the histogram of the constructed image has high amplitudes on high-value pixels which indicates a bright image. Also the CDF of the reconstructed image followed the desired CDF.

In this work, i will give myself a grade of 10 for doing it right. Good job for myself.

Activity 2: Area Estimation for Images with Defined Edges

2009-06-23T09:36:00.001-07:00

In this activity, we were tasked to draw shapes on MSPaint. The images are binary. The image themselves have pixel value of '1'. The background having '0'. The images drawn are shown below.

Figure 1. Circle (Area = )

Figure 2. Rectangle (Area = 252x140 sq.px)

Figure 3. Triangle

Figure 4. Irregular shapes

The main goal of the task is to take the area of the images comparing the measured area by just counting the 1's and by using Green's theorem given below.

The summation is over the contour of the image.
The results are shown below.

Table 1. Area calculations using pixel counting and Green's formula

Here, we can see that the Green's formula is accurately able to determine the actual area of the figure. Even when the shape is irregular as in Fig 4. as long as the contour is well-defined, the Green's theorem can correctly calculate the area enclosed by the contour.

In this activity, i will give myself a grade of 10 because i was able to correctly determine the area using Green's theorem.

Activity 3: Basic Types and Basic image enhancements

2009-06-22T20:14:00.000-07:00

In first part of this activity, we are tasked to find different modes of images,

Figure 1. True color image

Dimensions: 504 x 390
Width: 504
Height: 390
Horizontal Resolution: 72 dpi
Vertical Resolution: 72 dpi
Bit Depth: 24

Figure 2. Grayscale image

Dimensions: 600 x 400
Width: 600
Height: 400
Horizontal Resolution: 230 dpi
Vertical Resolution: 230 dpi
Bit Depth: 8
Resolution Unit: 2

Figure 3. Indexed image

Dimensions: 600 x 410
Width: 600
Height: 410
Bit Depth: 8

Figure 4. Binary image

Dimensions: 600 x 400
Width: 600
Height: 400
Bit Depth: 1

The images were taken from the following sources:

1(left): http://www.thesunblog.com/gourmetgal/fruit.jpg
1(right): http://www.alexmorais.com/Blog%20images/cat-grayscale.jpg
2(left):http://www.gingerbreadartgallery.com/images/archi/large/Gingerbread-house.gif
2(right): http://www.codeproject.com/KB/GDI-plus/AnnotatedImageV2/Mask.png

The binary image was really hard to find.

For the second part of the activity, a scanned image was used. Here we are tasked to apply thresholding to the image to:
(1) enhance the contrast of the image and
(2) to remove unnecessary pixel values on the region of interest (ROI).

In the figure below, the ROI is the name "JAYA".

Figure 5. (left) Scanned image and (right)pixel histogram of the image

Scanned images are usually in RGB mode. In this figure, the mode of the image was changed to grayscale using GIMP. The histogram in Fig.5(left) shows the pixel value distribution on the image (from 0 -255 pixel values). Here we can see the distinct separation between the high value pixels (white) and the low value pixels (black). As such, our ROI is sufficiently separated from the background.

One can also use Scilab to plot the histogram of the image using this code:

g = imread('(path of the image)');
g = int(255*im2gray(g));
[x,y] = size(g);
b = 0:255;
a =[];
for i = 0:255,
a(i+1) = sum((g == i)*1);
end;

plot(b,a/(x*y));

Using this we obtain:

Figure 6. Pixel histogram of image in Fig.5(left) using Scilab

Now, we apply GIMP to threshold the image. In thresholding, we set a certain value (threshold) such that any value less than the threshold will be discarded (pixel value is set to zero when dealing with images) and anything above will be retained or set to a certain value depending on the purpose.

Figure 7. (left) Raw image and (right)the thresholded image

Using a threshold, we are now able to completely separate the ROI from the background.
We can area calculations learned in Activity 2 to Fig.7(right).

We have calculated the area of the ROI in pixels using Green's theorem and pixel counting. Below is the same image with an attached ruler alongside with it. Using this ruler, we can calibrate the pixel-mm ratio which allows us to determine the actual area of the ROI.

Figure 8. Original image with attached ruler

Using MS Paint, the pixel-mm ratio was calculated.

pixel-cm ratio = 313 pixels:(262 mm - 156mm)
= 1 pixel : 0.339 mm

Hence, the length of a pixel is equal to 0.35 cm. The area of a pixel is then:

A = (0.35 mm)² = 0.115 mm²

Using the pixel counts earlier, we obtain these actual areas.

In this activity, I will give myself a grade of 8 for finishing this blog late.

Activity 1 - Digital Scanning

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

Figure 1. Graph taken from Introduction to Physiological Chemistry by Meyer Bodansky, 1930.

In this activity, from the graph (Fig. 1) we found last meeting, we are to replicate it using Paint and Spreadsheet. From the pixel locations (x,y) of the points on the plot, we record it to the spreadsheet and perform scaling operations to the values to replicate the physical values found in the graph.

To scale the axis:
1. Take the x-pixel locations (x- if scaling the x axis, y- if scaling the y-axis) of the smallest value and largest value on the respective axis of the graph (i.e. 4.7 and 12.2).
2. Solve for the the average number of pixels in between grids:
pbg = (pixel location of max value on the axis - pixel location of min value on the axis)/[(max value -min value)/value increments]

For scaling on the x- axis:
physical value = min x value + x-pixel location* 0.5/pbg

For scaling on the y-axis:
physical value = max y value - y-pixel value*1/pbg
(Scaling on the y- axis is different from the x-axis since the origin (0,0) of the image is different from the origin of the Cartesian coordinate system.)

The numerical values on the above equations accounts for the physical value spacing between grids. Also note that pbg on x and y axis are not necessarily the same. The graph below is the superimpose plot of the original figure and the digitally scanned figure.

Figure 2. Superimposed plot of the original figure and the digitally scanned figure. It can be observed that the scanned figure fits the original one.

Th colored dots represent the digitally scanned points from the original graph. We can see that the points coincide with original one. I give myself 10 for doing a good job!

The search for the oldest plot...

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

Today, we were tasked to find the oldest plot we can get. Preferably a hand written one. I got a year-1930 plot. I found it in a chemistry book entitled Introduction to Physiological Chemistry by Meyer Bodansky PhD.