Pixel of Knowledge

finally found the meaning of ‘I Like it on the Floor, I like it on my desk..I like it anywhere’ Facebook updates

2010-10-07T17:53:00.000+08:00

Ok. I'll just go direct to the point. Read this article http://thenewsoftoday.com/about-the-i-like-it-on-the-floor-i-like-it-anywhere-facebook-status-updates/3239/

October is Breast cancer month and women posting these updates want to show their support. That's it.

lessons from drinking sessions part 2..

2010-10-07T17:41:00.000+08:00

Second round of drinking session, this time with my good friend wert, we talked to our friends dad and he shared with us several key points. Below are the things that I remember :).

Key thing that I remember is "never be complacent, strive for greatness and avoid the trap of mediocrity".
Other things he also mentioned are the following:

1. too much analysis leads to paralysis - deal with problems, don't think and analyze too much.
2. marijuana is good. smells good - hehe.. :)
3. be frank, use your position, directly say what you want.
4. experience life, everything is a dose of experience.
5. make mistakes while young...
6. when talking to a crowd, start from where you came from - your foundation and branch your ideas from it. In this way, you are confident that you know what you're talking.

7. In correspondence to number 6, you must establish a solid foundation. whatever that maybe..
8.

..there are still couple of things. I just don't remember, :) hehe

lessons from drinking sessions..part 1.

2010-10-03T23:31:00.002+08:00

Nakausap ko si Fidel Rillo habang ako ay nakipag-inuman sa bahay ng aming kaibigan (tawagin natin Mambo) bilang despidida nya papuntang Germany. Nais kong ibahagi ang mga bagay na napagusapan namin.

Tinanong ako ng "What is God?"

napakahirap na tanong, isang tanong na kahit alam mong iba ang nagtanong sayo ay isa iyong tanong sa sarili mo. Bilang hindi ko alam kung ano ang tamang sagot, sinagot ko na lang ay..uhh, hindi ko pa pinagiisipan yan. hehe
Isang sagot ang binigay sa akin ng tatay ng aking kaibigan. "it is the search for truth"
Ako ay isang estudyante ng agham, at kung gagamitan mo ng rason, mahirap mapatunayan kung ano ang tamang sagot. "its the search for truth" -- sa gantong paraan, maaring pagdugtungin ang dalawang bagay na tila magkasalungat.
It is the search for truth..science and God...

isa itong palaisipan sa akin, isang bagay na dapat kong sagutin...

abangan ang iba pang mga lessons from drinking sessions :)

Simple way to fight Dengue

2010-10-03T22:49:00.002+08:00

I got this email from my colleague and it is worth trying. I recommend doing this at home and provide feedback if it will work.

I don't know where this came from but below is an exact copy of the email. Thank you to whoever made this. I haven't proven its worth but it should be worth a try.. here it goes...

Its just a mix of water, brown sugar and yeast.

Cut a plastic bottle in half, keep both parts. Can be soft drink bottle.
Take the lower portion of the bottle. Dissolve the brown sugar in hot water. Let it cool down to ~70 deg F.
Add the yeast. Carbon dioxide will form (This will attract the mosquitoes)
Cover the bottle with a dark wrap and insert in the top portion upside down like a funnel. Place it in a corner in your house.
In 2 weeks you will be surprised by the number of mosquitoes killed.

So there you go. It sounds simple and I hope through this it will help some families out there who are in danger of mosquito diseases..

Activity 19: Restoration of blurred image

2009-10-10T22:36:00.000+08:00

In the previous activity, we only considered adding noise but without the application of the degradation function. In this activity, we perform restoration for an image degraded by a known function h(x,y) and added with noise n(x,y) (see FIG 1 of activity 18 for the image degradation and restoration model).

FIG 1: Degradation process of an image f(x,y) by convolution of the degradation function h(x,y) with the image f(x,y) and addition of nosie n(x,y). By the convolution theorem this is represented as multiplication in the frequency domain represented in Eq (2) where capital letters denote the Fourier transform of the smaller letters in Eq (1). The '*' operator denotes convolution.

In this activity we consider uniform linear motion blur and obtain an estimate of the analytical form of this kind of blurring. Thus, we find the transfer function H(u,v) of the motion blur.

Estimate of the degradation function, h
Assuming that the opening and closing of the shutter of the imaging device occurs simultaneously and considering perfect optical imaging process, the degraded image g(x,y) in the frequency domain can be represented by Eq. (2) where the transfer function H follows the form below.

FIG 2: Estimate of the transfer function for uniform linear motion assuming perfect optical imaging and instantaneous opening and closing of the shutter of the imaging device. Where T is the duration of the exposure, a and b are the total distance for which the image was displaced in the x and y direction respectively.

Restoration of the degraded image g(x,y) to f^(x.y)
In this activity, we use the Wiener filter which minimizes the mean square error between the restored image (f^) and the original image f. i.e., min{|f^-f|} The Wiener filter assumes that the image and noise are random processes and ncorporates the statistical properties of the noise and degradation function in the restoration process. The minimum of the mean square error is given below in the frequency domain as:

FIG 3: Wiener filter. H(u,v) is the degradation function from Eq. (3), |H(u,v)| = H*(u,v)H(u,v) where H*(u,v) above denotes complex conjugation, Sn(u,v) = |N(u,v)|^2 and Sf(u,v) = |F(u,v)| ^2 are the power spectrum of the noise and the undergraded image.

FIG 4: Transfer function H (frequency space, left) and the equivalent function h in real space, (right) for T=1, and a=b=0.01.

The image below shows the results for an image degraded by H and added with a Gaussian noise of mean = 0.01 and var = 0.01 generated using scilab's grand function.

FIG 5: Test image degraded by H following the form of Eq. (3) with T=1, a=0.01, b=0.01 and added by a Gaussian noise of mean=0.01 and variance = 0.01. The restoration using Wiener filter of Eq. (4) is visually successful.

Usually when we are dealing with real images, the power spectrum Sf(u,v) of the undergraded image is unknown. We can estimate Eq. (4) by assuming a white noise with a constant power spectrum reducing Eq. (4) to the equation below.

FIG 6: Estimate of the Wiener filter.

For different values of K, we reconstruct the degraded image as shown in the images below.

FIG 7: Reconstruction for varying values of K = [1e-4 1e-3 1e-2 1e-1 1] for degraded image with T=1, a=b=0.01 added with Gaussian noise of mean = 0.01 and variance=0.01.

Comparing the Wiener filter reconstruction with the 'estimate of the wiener filter' , the Wiener filter yields better reconstruction as compared to the 'estimate wiener filter'. The reconstruction for the 'estimate of the wiener filter' has good quality when the value of K is close to zero.

For other combinations of a, b and T, the following reconstruction results.

FIG 8. Reconstruction for different combinations of T and a=b. Generally, increasing T increases the quality of the image while increasing a=b decreases the quality of the image.

From the above image, increase in T leads to better reconstruction. This result is most logical since increasing the exposure time of the imaging device increases the number of photons detected and thus better image quality. Since a and b indicates the amount that the image is displaced, increasing a and b will lead to further blurring leading to poor image quality.

In this activity, I would give myself a grade of 9 for not knowing if the result is correct. I have doubts with the appearance of my transfer function and from the picture of the blurred image, it does not appear to be resembling that of a uniform linear blur. Instead, blurring seems to occur in both the x,y directions as oppose to single direction. The image cannot possibly move in both axis simultaneously. It is recommended that the form of the transfer function be checked for future research.

Acknowledgements
I would like to acknowledge Mark Jayson Villangca, Luis Buno and Miguel Sison for useful discussions.

References:
[1] Applied Physics 186 Activity 19 manual.

[2] image source

Activity 18: Noise model and basic image restoration

2009-09-22T10:47:00.000+08:00

In image restoration, the objective is to recover the original image f(x,y) from the observed image g(x,y) with an a priori knowledge of the degradation process h(x,y) and the noise n(x,y) added on the image. This kind of modeling is shown in the figure below.

FIG 1: Image degradation and restoration model. The original image f(x,y) undergoes degradation by a known process h(x,y) and added with noise n(x,y) leading th the observed image g(x,y). With the use of restoration filters, image restoration aims to recover f(x,y). [1]

In this activity, we only cover half of the image restoration modelling by adding a known noise n(x,y) to the image and performing restoration algorithm to recover the original image. The strength and properties of the noise n(x,y) can be characterized from its probability distribution function PDF shown below for different noise models.

FIG 2: PDF of the different noise models used characterized by its mean and variance. (click to zoom)

The different noise models were generated using scilab's builtin grand and imnoise function. The Rayleigh noise was generated using a module modnum downloaded from [2] while the other noise were generated using grand and imnoise (see scilabs documentation for more detail on these functions). Below is the resulting image after adding different noise to a test image.

FIG 3: Image histograms after addition of different noise models.

Restoring in the presence of noise

Below are the different noise-reduction filters used when only noise is added to the image [1].

FIG 4: Mean filtering methods for a window of size m x n.

Consider a window, i.e., 3x3 window, we perform the 'mean filtering' on this window. For each pixel in the noisy image g(x,y), we choose a window and perform mean filtering for that pixel. This was done for all the pixels in the images. In order for the algorithm to be effective around the boundaries, the image was tiled in the manner below:

FIG 5: Tiled image to wrap the boundary conditions. The red window represents the area where the spatial noise filtering was applied.

The following pictures below show the result of the restoration process using the different mean filtering methods in FIG 4 for each noise model. The window used is 5x5 and the restoration for the contraharmoic was divided into two depending on the value of Q, {Q=2 & Q=-2}.

FIG 6: Image restoration results for an image added with a Gaussian noise. Looking at the histogram, the restoration was generally successful for all methods except for some error on the shift of the graylevel value of the three peaks in the original image. In particular, geometric, harmonic and contra-harmonic (Q=-2) appear to shift the graylevel value to the right (towards white) relative to arithmetic and contra-harmonic (Q=2). This means that the image shifted to the left will appear darker relative to those shifted to the right.

FIG 7: Image restoration results for an image added with a Rayleigh noise. The restoration of the arithmetic and geometric filtering appears visually correct while for the contra-harmonic filtering, the image is either shifted to the left (Q=2) or right (Q=-2) depending on the value of Q used. Unlike geometric and arithmetic filtering, the histogram of harmonic and contra-harmonic filtering is wider along the peaks of the original image.

FIG 8: Image restoration results for an image added with a Gamma noise. In general, the restoration appears successful except for the contraharmoic (Q=2) restoration where the histogram of the image is more wider.

FIG 9: Image restoration results for an image added with an exponential noise. The contraharmonic restoration (Q=2 & Q=-2) failed in restoring the image due to wide histogram and overshifting to the left of the grayvalues. This degree of shifting is also manifested in the harmonic filtering.

FIG 10: Image restoration results for an image added with a uniform noise. In general, the restoration is quite successful for all methods.

FIG 11: Image restoration results for an image added with a salt & pepper noise. The restoration seems to have artifacts - presence of extra graylevel values. This could be attributed to the strength of the noise used which might be too strong and filtering on a 5 x 5 window is not sufficient to recover the image.

FIG 12. Applying the restoration process to a grayscale image shown above obtained from my brothers collection of images.

Applying this restoration to a grayscale image (obtained from my brother's picture collection), we perform the same method, adding a known noise term to the image and restoring the image using the mean filters of FIG 4.

FIG 13: Grayscale image added with Gaussian noise with the corresponding restorations.

FIG 14: Grayscale image added with Rayleigh noise with the corresponding restorations.

FIG 15: Grayscale image added with Gamma noise with the corresponding restorations.

FIG 16: Grayscale image added with an exponential noise with the corresponding restorations.

FIG 17: Grayscale image added with a uniform noise with the corresponding restorations.

FIG 18: Grayscale image added with salt & pepper noise with the corresponding restorations.

In general, geometric and arithemetic filtering work for all kinds of noise models presented. In particular, the author recommends geometric mean as it seems to be flexible and suitable enough for all noise models studied above.

In this activity, I give myself a grade of 10 for performing the restoration and doing all the required activities.

Acknowledgement

I would like to acknowledge Gilbert Gubatan for providing the downloaded modnum module for generating Rayleigh noise and my brother for lending one of his collection pictures.

References
[1] Image Restoration. fromhttp://www.ph.tn.tudelft.nl/~lucas/education/et2720in/2002/ImageRestoration.pdf
[2] http://www.scicos.org/ScicosModNum/modnum_web/src/modnum_421/interf/scicoslab/help/eng/htm/rayleigh_sim.htm

Activity 17: Photometric Stereo

2009-09-17T19:24:00.000+08:00

Photometric stereo is a technique that extracts surface elevation by computing the surface normals assuming the intensity is directly proportional to brightness. Using N light sources (V) , the surface normal of each point in the image can be computed following the steps in the image below. (click to zoom for clearer picture)

applying these steps to four light sources with images below, we reconstruct the corresponding shape elevation.

The result is shown below:

The reconstruction is relatively correct but with 'horns' on top and on the edges. Also the image is not smooth and some filtering is recommended.

In this activity, I give myself a grade of 10 for successfully reconstructing the image.

I would like to thank Martin Tensuan for useful discussions.

References:
[1] Applied Physisc 187 Manual (c)2008 Maricor Soriano

Activity 16: Neural Networks

2009-09-14T03:19:00.001+08:00

In activity 15, we used Linear discriminant analysis (LDA) as predictor to classify chronic lymphocytic leukemia (CLLs) from lymphocytes using two different sets of features, {area perimeter} and {%red %green}. We already know that using perimeter and area as features would yield 100% classification for LDA and minimum distance algorithm while using %red and %green as features yield 90% success rate for LDA. In this activity we used neural network as predictor to classify CLLs and lymphocytes using the %red and %green as features. Similar to activity 15, we used 5 images of both CLLs and lymphocytes as training set and 5 images from both classes as test set. (see fig 2 and fig 3 of activity 14).

Neural network is a computational model that mimics the behavior of biological neurons, the brain for example[1]. Essentially, the idea of neural network is to assign weights on the input (say features), sums them and pass it to an activation function that processes the weighted input into an output. The process of determining the weights of the input is an iterative process with an a priori knowledge of the output. What neural network does is to compare its output with the desired output and adjust the weights of the input depending on the deviation from the desired output. This process is called the learning process of the network. After the network learns how to classify objects based from their features, classification of the test set is straightforward by feeding the network with the objects' features. For detailed discussion on the subject matter, refer to article in [2]

In order to understand more on how neural networks are modeled, fig 1 shows an artificial neuron as modeled by McCulloch and Pits (1943). A connection of many such neurons comprises a network.

Figure 1: Artificial neuron (left) and many connections of an artificial neuron leading to a neural network (right). It typically consists of an input layer followed by a hidden layer and an output layer.

In this activity, we used the ANN (artificial neural network) toolbox of scilab in performing artificial neural network computation. The process is straightforward and the parameter of interest adjusted are the learning rate and the training cycles. Table 1 shows the result of the computation for learning rate of 0.1, i.e., l=(0.1 0) and training cycle of 10000

Table 1: Neural network classification yielding 100% success rate for learning rate=0.1 and training cycle = 10000

We investigate the effect of learning rate and training cycle on the output of the neural network. We used the lymphocytes as test object, if the network's output is close to 1 the classification is correct otherwise it is wrong. For a constant learning rate (i.e 0.01) and increasing training cycle (1000-10000) we expect the output of the network to approach a limiting value as shown in fig 2. This implies that for increasing training cycle, the neural network is improving and is able to classify the object better.

Figure 2: Effect of different training cycle on the neural network's output for a constant learning rate(0.01). As shown, increasing the training cycle increases the output approaching the value 1 indicating that the network is learning better.

For a constant training cycle (1000) , and increasing learning rate (0.1-1), the result also approaches the limiting value of 1 as shown in figure 3.

Figure 3: Effect of different learning rate on the neural network's output for a constant training cycle (1000). As shown, increasing the learning rate increases the output approaching the value 1 indicating that the network is learning better.

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

I would like to acknowledge Luis Buno for very useful discussions and for mentioning the idea of plotting the output for different training cycle and learning rate.

References
[1] http://en.wikipedia.org/wiki/Neural_Network#Learning_paradigms
[2] C. Bishop. Neural Networks and Their Applications. Rev. Sci. Instrum. 65(6),
June 1994.
[3] Applied Physics 186 Activity 16 (c)2008 Maricor Soriano

Activity 15: Probabilistic Classification

2009-09-14T01:40:00.000+08:00

In the previous activity we used minimum distance classification as our predictor in classifying objects. Essentially, what we did is to calculate the deviation of an object from the mean features of different classes. However such classification fails when used with features that overlaps in both classes. In this activity we used Linear Discriminant analysis as predictor in classifying objects. We limit the number of classes to two and use the same features as the previous activity.

Consider the plot below of the features in phase space plot.

Figure 1: Phase space plot using area and perimeter. The regions are well segmented and we expect that simple minimum distance classification should work well as demonstrated in activity 14.

Figure 2: Phase space plot using red and green composition of the image. As can be seen, there are overlaps in both coordinates for both classes such that when the features (independent variables) are used alone, the classification would fail.

LDA transforms a class as a linear superposition of features and assumes therefore that features are linearly independent from each other. The aim of LDA is to minimize the error of classification that an object belongs to a certain class. The classification rule is that an object belongs to a class if it has the highest conditional probability in that class. The probability of "belongingness" can be computed using the formula below. It is worthwhile to state that the equation below was derived with the assumption that each class has normally distributed conditional probability[1-3].

where

ui: mean features in group i.
C: pooled covariance defined as

where ci is defined as the covariance matrix of group i, n is the total number of samples and ni is the number of samples in group i

Pi: prior probability usually assumed to be the number of samples in the group over the total number of samples, ni/n
xi: features data of group i.
xk: feature data of row k for a certain group.

Using the above formulas we compute for fi and determine the class with the highest f for each test object using the features in figure 1 and 2.

Table 1: Calculated values of f1 and f2 for different test objects using area and perimeter as features.

The classification is 100% as expected since the phase space plot of area and perimeter are already well segmented as shown in figure 1.

Table 2:Calculated values of f1 and f2 for different test objects using red and green composition as features.

The classfication failed 1/10 resulting to 90% success classification rate. This error is most probably due to the limited number of training sets (5) and the assumption of normal distribution conditional probability might not yet be satisfied.

In this activity, I give myself a grade of 10 for understanding a bit on LDA and performing the classifications well.

I would like to acknowledge Mark Jayson Villangca and Jay Samuel Combinido for useful discussions.

References:
[1] http://people.revoledu.com/kardi/tutorial/LDA/Numerical%20Example.html
[2] Applied Physics 186 Activity manual (c)2007 Maricor Soriano
[3] http://en.wikipedia.org/wiki/LDA

Activity 14: Pattern Recognition

2009-09-10T11:03:00.000+08:00

In pattern recognition, the idea is to mimic the human brain in classifying objects into class. Similar to object oriented programming, an object has methods or attributes called features which are common to objects belonging to a certain class. Pattern recognition has two tasks, (a) finding the appropriate features such as color, area or perimeter for efficient separation of classes and (b) finding a decision function that can classify objects into class. We can define the features as quantifiable properties and combine them together to form a feature vector. From the feature vector of objects, we use a method/function to classify them.

In this activity we employ the minimum distance classification to classify objects. We use five features: area, perimeter, % red, % blue and % green in classifying chronic lymphocytic leukemia (CLL) cells from lymphocytes (T cells or B cells) [1]. In order to classify objects, an a priori knowledge of a class and its set of features must first be establish. We can think of it as training the classifier for further classification later on. Below is the image of a collection of CLL cells and lymphocytes which we will use in our classification algorithm.

Figure 1: Chronic lymphocytic leukemia cells (i.e., those pointed by the arrows, blue) and lymphoytes (light violet).

From the above image, five cropped images of lymphocytes and CLL cells were used as training sets and five more images of both classes were used as test sets. In classifying the objects, the minimum distance classification is used which can be represented using the equation below:
where x is an array where each row represents an object and each column represents a feature, m is the mean feature for a certain class W. For each class W, we can obtain d and an object can be classified to a class if it has the largest value of d.

Below are the images used for the training sets and test sets for both classes (CLLs and lymphocytes). In order to extract the area and perimeter features, thresholding was used to segment the images (b) and bwlabel was used to clean the images from isolated spots (c) (see fig 2). Pixel count and the scilab's builtin command follow was used to extract the area and perimeter respectively.

Figure 2: Training set of images for both CLLs and lymphocytes with the corresponding image processing operations used (a) rgb (b) bw/thresholding and (c) bwlabel.

Figure 3: Test set of images for both CLLs and lymphocytes with the corresponding image processing operations used (a) rgb (b) bw/thresholding and (c) bwlabel.

The table below shows the mean features obtained from the training set. It can be seen that the area and perimeter features segments the classes pretty well as shown in figure 4.

Table 1: Mean features of CLLs and lymphocytes.

Figure 4: Phase space plot of area and perimeter as major classifiers for lymphocytes and CLLs.

The table below shows the result of the classification after computing d for both classes (CLL and lymphocytes) applying on the test sets of figure 3. The algorithm correctly identifies 10/10 objects giving 100% correct classification.

Table 2: Computed d of lymphocytes and CLLs.

In this activity, I give myself a grade of 10 for attaining the objectives of the activity, correctly classifying objects.

I would like to acknowledge Jay Samuel Combinido for very useful discussions.

References
[1] http://www.vet.uga.edu/vpp/clerk/waikart/index.php
[2] Applied Physics 186 activity 14 manual, (c)2008 Maricor Soriano

Activity 13: Correcting Geometric Distortions

2009-09-10T10:39:00.000+08:00

Inherent in most imaging device is the presence of diffraction and aberration due to the finite collection aperture and geometry of the lens. Because of the shape and composition of the lens (i.e spherical), magnification and focusing is not isotropic resulting to unwanted distortions. Two of the most common distortions are pincushion and barrel distortion (fig 1).

Figure 1: Two kinds of radial distortion. (a) Barrel distortion appears expanded at the middle and pushed at the sides while (b) pincushion distortion is of the opposite, pushed at the middle and expanded at the sides. Images are taken from [1].

Distortions are usually undesirable in image processing especially in machine vision where slight distortions have major drawbacks [2]. Intuitively, a distorted image is the result of a transformation from an ideal image. Essentially, the process in correcting image distortion is to find the transformation from the ideal image to the distorted image.

Using a reference grid, we can correct image distortion as follows:

Figure 2: Transformation from ideal to distorted image of a grid. (a) Left distorted image can be mapped to the (b) ideal image using the transformation r and s assuming linear transformation in both directions simultaneously. This scheme is valid as long as the grid is not large relative to the image.

Consider fig 2, given the coordinates (x,y) of an ideal image (i.e., reference grid) we can find the distorted image coordinate x',y' using the transformation r, s by finding the coefficients c1-c8. The problem now reduces to solving systems of equation with eight unknown, c1-c8 (see fig 3). The eight equations can be obtained from the coordinates of the vertices (four corners) of the grid. It must be noted that the transformation is only valid inside the grid and the coefficients must be recomputed for a different grid.

Figure 3: Finding the coefficients c1-c8. From the matrix transformation T (x1-x4 and y1-y4 can be obtained from the four corners of the grid) we can obtain C's by matrix inversion of T and multiplying them with the x's and y's coordinate of the ideal grid.

Per pixel in the ideal grid, we can find the corresponding distorted coordinate r,s using the computed coefficients, c1-c8 (fig 2).

After finding the corresponding distorted coordinate, we have to find the graylevel value of the ideal image. If the computed distorted coordinate is integer valued, we just have to carry the graylevel value of the distorted image to the ideal image otherwise we have to interpolate. We use two kinds of interpolation, nearest neighbor and bilinear interpolation.

The nearest neighbor interpolation uses the graylevel value of the nearest integer-valued coordinate in the distorted image. This kind of reconstruction is very fast but the quality is poor and is ideal for thumbnails and images with relatively not varying graylevel value. On the other hand bilinear interpolation assumes linear interpolation in both directions simultaneously similar to what we did above. The graylevel value v can be obtained as follows (fig 4).

Figure 4: Finding the graylevel value v of the distorted image using the four neareast neighbors encompassing the distorted coordinate. The coefficients abcd can be found in a similar manner using fig 3.

Below is the result of the reconstruction using checkerboard and a distorted image of a window for nearest neighbor interpolation. To find the reference grid coordinates of the distorted image, scilab's locate was used to manually obtain the data points. Although this is very taxing, locating reference grid coordinates embedded with images is sometimes a formidable task. It is highly recommended though to obtain the coordinates using image processing techniques.

Figure 5: Nearest neighbor interpolation of (top) checkerboard and (bottom) distorted window. Typical of nearest neighbor interpolation is the presence of jagged regions due to the nature of the algorithm (picking the nearest neighbor). click to zoom.

For the checkerboard, nearest neighbor interpolation is already acceptable since the image is only black and white but for the window, the details of the image was slightly diminished and we need better interpolation scheme.

In bilinear interpolation, the quality is better compared to the nearest neighbor scheme but the speed is very much slower, figure 6 shows this result.

Figure 6: Bilinear interpolation of (top) checkerboard and (bottom) distorted window. The quality of the image is better compared to the nearest neighbor interpolation. click to zoom.

In this activity, I give myself a grade of 10 for correcting the distortion thereby achieving the goal of the activity.

I would like to acknowledge Mark Jayson Villangca, Jay Samuel Combinido and Luis Buno for all the help and jokes in doing this activity.

References
[1] http://en.wikipedia.org/wiki/Image_distortion
[2] Applied Physics 186 manula (c)2008 Maricor Soriano

Activity 12: Color Image Segmentation

2009-08-06T16:15:00.000+08:00

Segmentation by thresholding will not always work especially when the area of interest has the same grayscale value as the background. In such cases, segmentation can be done using the color of the image of interest. In this activity, we're going to demonstrate parametric and non-parametric color segmentation.

Normalize Chromaticity Coordinates (NCC)

Before we proceed with color segmentation, it is convenient to convert the RGB space into NCC space. NCC has the advantage of separating the chromaticity and intensity of the image which is ideal for representing 3D objects. Essentially, a 3D object obtains its 3D effect through shading variations.

To convert from an RGB space to NCC space we do the following: for each pixel,

Note that the chromaticity of the image is essentially reduced to two dimensions since b can be obtained from r and g usng b=1-(r+g). The image below represents the normalize chromaticity coordinates where the x-axis represents the r and y-axis represents the g.

Normalize chromaticity coordinate.

Color Segmentation: Parametric vs Non-parametric

To segment the image using its color, we must determine the probability that a pixel belongs into the color of interest. To do this, we need to extract the histogram of the color of interest and normalize it to the number of pixels to obtain the normalized PDF. In the NCC space, chromaticity is reduced to two dimensions, r and g such that we need to obtain the PDF for each space. If we assume that r and g are independent, we can obtain the PDF for each space separately and multiply them to obtain the overall PDF of the image.

Parametric Segmentation

Parametric segmentation segments the color by fitting a PDF in the image. Assuming a Gaussian distribution independently along r and g axis, we can obtain the PDF by using the equation below: mu is the mean value and sigma is the standard deviation of the color of interest, r is the NCC r-space of the image.Essentially, what this PDF does is to highlight all those pixels having pixel value near the mean and darkens those pixel values that deviates from the mean. Similarly we can obtain a similar PDF for the NCC-g space. The joint probability is obtained by multiplying p(r) and p(g). This PDF then contains all the highlighted pixel values near the desired color of interest.

Nonparametric Segmentation

The drawback of parametric segmentation is that it assumes a PDF independent of the pixel values of the image. In non-parametric segmentation, we base the segmentation base on the histogram of the color of interest itself and use this to backproject it to the image. Essentially, we obtain the 2D histogram of the color of interest by first converting it to NCC space. The result of this histogram then is a 2D matrix with the r-axis (representing the x-axis), g-axis (representing the y-axis) and the intensity (representing the count, frequency). For each pixel in the image, we backproject the obtained histogram. The steps are as follows:

obtain histogram of the color of interest, convert to NCC space first
for each pixel in the image, obtain its NCC r and g values
using the r and g values, find the pixel value from the histogram
replaec the pixel value of the image with the value obtained in the histogram

Below are reconstructed images of parametric and non-parametric segmentation.

Parametric reconstruction, observe that the segmentation is not perfect and contains two or more colors.

Non-parametric segmentation. The images are segmented quite well except for the brown and yellow hounds. Image was taken from source [2]

Parametric reconstruction. Note that for the gren shoe, artifacts of blue shoe are present in the segmentation. Image was taken from source [3].

Non-parametric segmentation. Note that the segmentation is clearer as compared to the parametric segmentation.

In general, non-parametric segmentation has better results in that it was able to segment the colored objects solely and cleanly as compared to the parametric segmentation. This is expected since in the paramtric segmentation, we assumed a Gaussian PDF whereas in the non-parametric segmentation, we base our segmentation solely from the image itself. It must be noted however that non-parametric segmentation highly depends on the cropped color of interest. For this reason, it is probably convenient to apply parametric segmetation for those that require automatic color segmentation.

In this activity, I give myself a grade of 10 for segmenting the images properly.

Acknowledgement
I would like to acknowledge jaya for useful discussions regarding the non-parametric segmentation

References
[1] App Physics 186 activity 12 manual
[2] www.greytsoaps.com/greys.htm
[3] http://www.whatsalltheracquet.com/archives/pictures/renelacoste.jpg

Activity 11: Color Image Processing

2009-08-06T09:55:00.000+08:00

When imaging using a detector (i.e., a camera), the amount of light detected depends on the sensitivity of the detector, the intensity of the illuminating light and reflectance properties of the object. Similarly in digital imaging, a pixel is compose of different proportions of red, green and blue values overlaid in different proportions. The equation below describes how the RGB values are obtained. It is an integral product of the reflectance spectra of the object, spectral power distribution of the light source and camera sensitivity.

The presence of the K values above normalizes the effect of camera sensitivity to obtain what we call white balancing. Note that if K is remove, cameras sensitive to certain spectral range will affect the overall color of the image. Also, the kind of illumination will affect the color of the object, i.e., incandescent light appears yellowish as compared to fluorescent light.

In this activity, we're going to apply two algorithms to achieve white balancing. They are the White Patch Algorithm and Gray World Algorithm.

White Patch Algorithm

Observe from equation above that the K's are essentially the inverse of the camera output when shown with a white object. Thus white balancing here is just dividing the raw camera output by the image of the white object. This is precisely the White Patch Algorithm.

Gray World Algorithm

In the gray world algorithm, it assumes that the average color of the world is gray. Taking the average of the red, green and blue channels can then serve as white balancing constants for the red, green and blue channels respectively. For a certain channel, when you divide the pixel values by the average value of that channel, it is probable that many pixels would contain values greater than 1 (i.e., if pixel value > ave. pixel value). By normalization, (divide by the max) we can lower this value to 1. However, this has the overall effect of darkening the image. As compared to the white patch which divides each channel by the pixel value of the 'white', most probably contains the highest pixel value of the image.

In summary, the algorithm is as follows

read the image
obtain the constants for white balancing

white patch algorithm - obtain a white patch in the image and average the pixel values of this patch for each channel. This will serve as the white balancing constants

for the gray world - obtain the average of each channel and use as the white balancing constants
note that before we do the averaging, we remove all pixels that are saturated.

divide each channel by the corresponding white balance constants
normalize the image by dividing the image by its maximum. Note that we avoided the use of clipping the image for in the gray algorithm, it is probable that for each channel many pixel values are greater than the average pixel value making the image look saturated when clipped.

Below are images obtained from different lighting conditions and corrected using the white patch algorithm and the gray world algorithm. It is arrange such that each column represents the white balancing used by the camera while for each row is the corresponding algorithm used. Raw means without the use of the white and gray algorithm.

The picture below was taken around 9AM in the CSRC garden. (Cloudy condition)

From the above image, the raw, white patch and gray algorithm has generally the same color hue but the gray patch has darker image compared to the white patch and raw. Generally, the raw image is still better since this image was taken using the correct white balancing of the camera.

The above picture has wrong white balancing set by the camera. The gray world algorithm was able to correct the image, adjusting the hue of the red and blue channels. The same can be said of the white patch algorithm. However, the image of the gray world algorithm appears darker as compared to the white patch. Again we attribute this to the constant we used in dividing the pixel values of the image and by normalization (see discussion above).

We tried the algorithm for different lighting conditions and the results are seen below.

Image of a table in the CSRC lobby (daylight)

Note that the white and gray algorithm was able to correct the "light blue hue" background of the raw image.

Note that the raw image has yellowish hue which was corrected by the gray and white algorithm.

Another image taken inside the room (i.e., fluorescent lighting condition)

Note that the result of the gray world algorithm and white patch algorithm brings out the color of the table better compared to the raw image.

In this picture, we can see that gray world algorithm is better compared to the white patch algorithm. In the raw picture, the image has yellowish hue background indicating that some channel is dominant over the other channels. Note that using the gray world algorithm, we are obtaining the average of each channel reducing the effect of dominance of other channels.

Another image inside the room of different hues of red. (fluorescent lighting conditions)
Observe the corrected gray and white algorithm using incandescent white balancing of the camera.

Generally, I think that the gray world algorithm is much convenient and sometimes better compared to the white patch algorithm. Note that unless, we can obtain an algorithm for finding the white patch image in the white patch algorithm, our reconstruction would always depend on the patch that we're going to use. As compared to the gray world algorithm which makes use of averaging. The drawback of gray world algortithm is that it makes the image look darker. However, this can be remedied by using contrast enhancement, i.e., use an exponential cumulative distribution function to lighten up the image.

In this activity, I give myself a grade of 10, for obtaining the reconstructions and explaining the results

Acknowledgement

I would like to acknowledge Kaye for lending me the McDo stuff toy and Carmen for the blue notebook.

References
[1] App Physics 186 Activity 11 Manual

Activity 10: Preprocesing Text

2009-08-06T09:46:00.000+08:00

When handwriting documents, usually the papers we are using are embedded with lines. In this activity, we are going to extract handwritten text from these kinds of papers. We're going to apply different image processing techniques to remove the lines and extract the text.

Scanned image of a piece of paper with handwritten texts.

We crop a portion of the above image and try to segment the handwritten texts individually.

The image was rotated using mogrify at an angle ~1.219. Click to enlarge the image

The lines can be removed by filtering in Fourier space by noting that the horizontal lines produce frequencies along the vertical axis in frequency space.

Filtering in Fourier space. Click to enlarge the image

The filtered image was further polish using scilab's sharpen and enhance (i.e., mogrify(img, ['-sharpen','2'])to highlight the letters better. Enhance is applied to remove the noise which will be useful when thresholing the imge.

Sharpening and enhancing the image for thresholding. Click to enlarge the image

Note that after thresholing, the image is noisy. We apply closing followed by opening operator to clean the image and segment the letters as much as possible.

From leftmost: (a) structuring element (b) result of closing operator and (c) result of opening operator. Click to enlarge the image

We use scilab's built-in thin function to segment the image into 1 pixel width. Scilab's bwlabel was applied to separate each continuous blobs which represent the segmented letters.

Segmented letters using bwlabel and thin. Click to enlarge the image

Note that although we are able to separate the letters into 1 pixel width, the segmentation is not perfect and needs lot of improvement.

Recommendation:
For the case of the image above, the handwritten letters were written using blue colored pen. Using color segmentation, we can separate the handwritten texts without having to worry about filtering. Morphological operations can then be used to segment the image better. It must be noted that most of the error in the above segmentation is due to thresholding after filtering the image. Segmentation by color might improve the reconstruction

For the next part of this activity, we find other instances of the word description in the original image. We simply correlate the original image with the image of the word description and find the locations where the pixel value is highest.
The images that we're going to correlate. (a) binarize image and (b) image of the word DESCRIPTION.

Result of the correlation. The highlighted pixels correspond to the word DESCRIPTION.The bright dots correspond to high correlation. Click to zoom in

In this activity, I give myself a grade of 8 for below average reconstruction. However, I recommended an alternate solution to the problem and for this reason I give myself a grade of 8.5 (bonus of 0.5) =).

References
[1] App Physics 186 activity 10 manual

Activity 9: Binary Operations

2009-08-01T16:49:00.000+08:00

In this activity, we're going to obtain the best estimate of cell area (pixel count) by using all the morphological and binary operations we have learned. Below is the image of the "simulated" cells that we're going to measure, punched papers imaged using a flatbed scanner.

Image of simulated cells, punched paper digitize using a scanner.

A quick look at the image tells us that we need to analyze several ROI (region of interest) in order to obtain the best estimate. Ofcourse, we could opt to obtain an isolated cell and perform area measurement but this will not always apply especially when the cells have areas that deviate slightly from one another. Also, the cells shape may not always be uniform. In order to obtain our region of interest, we could do thresholding and binarized our image for easier area measurements later on.

Below is the histogram of our image:

Histogram of the simulated cells.

We can see that most of the information in our image is roughly between 0.5 and 0.85. Base on this value, we could binarize our image using a treshold of around 0.85.

Statistically, if we can obtain several measurements of different cells, we can plot the histogram of our measurement and base from this histogram, we can obtain our best estimate for area measurement. This is precisely what we're going to do. We obtain several subimages (256 x256) from the original image and perform area measurement on each image. We do this by thresholding the image and perform morphological operations like opening and closing to separate nearly touching cells and to remove isolated spots. A literature on opening and closing operation is available in source [2].

In this activity, we obtain 20 subimages of sizes (256x256) randomly on the image.

Subimages of size 256 x 256.

The above image was thresholded at around 0.820 resulting to the image below:

Thresholded subimages.

Observe the presence of nearly touching cells and isolated spots on the subimages. We can futher clean this image by performing morphological operations. In this activity, we perform opening operation on the above image using a circle as structuring element having diameter = 8 pixels. Note that we have to be careful when choosing our structuring element for this could affect the sizes of the cells drastically.

Below is the resulting image after the opening operation.

Comparing to the "unopened image" above, we can see that the images are relatively clean and we are able to separate several nearly touching cells. Although, this separation is not 100% successful, we are somehow able to enhance our image for area measurement.

After performing morphological operations on our image, we are now in the position to obtain area measurement by labeling all contiguous blobs and couting their areas (pixel count). We can do this in scilab using the function bwlabel which labels all contiguous blobs and tabulates them in an array. Below is the resulting images after applying bwlabel.

Labeled subimages using scilabs bwlabel.

By looping on all the subimages and counting the area of each contiguous blobs, we can obtain a frequency distribution of the area of the cell. In our measurement, we disregarded area measurements greater than 800, since by visual investigation, the cells' pixel area is clearly below this value. Below is the histogram of the obtained cell areas.

Histogram of cell measurement

From the above histogram, we obtain an estimate of area = 530 pixels plus or minus 10 pixels since we use bins of spacing 20 pixels in our histogram. In the above histogram the number of cells having areas between 520-540 is 87, significantly higher than the other measured area.

Since our approach is statistical, the higher the number of subimages, the better. Although, we could perform morphological operations on nearly touching cells, it must be stressed that the approximation is much easier if we can obtain images of cells that are relatively sparse. We must think of this when we obtain images of cells in our experiment.

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

Acknowledgement
I would like to thank Irene for very useful conversations.

References

[1] AP 186 Activity 9 manual
[2] http://en.wikipedia.org/wiki/Mathematical_morphology#Opening_and_Closing

Activity 8: Morphological Operations

2009-07-23T19:17:00.000+08:00

An image contain informations that can be extracted using image processing. Usually before image processing, we try to enhance the input through morphological operations. Morphology refers to shape and structure manipulation, thinning, expanding, closing of holes and the like. Because of its versatility,morphological operations find many applications in area measurement, counting, and tracking to name a few [1]. Morphological operations make use of set theory. A brief rundown of set theory is available in source[2].

In this activity, we're going to demonstrate two basic morphological operations, dilation and erosion. We proceed by defining dilation and erosion using set theory and provide examples using different images. To futher understand dilation and erosion better, the reader is recommended to predict the outcome of dilation and erosion on a piece of paper before programming it on a computer.

DILATION
From our activity manual, "dilation of A by B denoted by A dilation B is defined as".

It is important to note that dilation of A by B is the set of all Z's such that B translated by Z union A is not an empty set. B here is called the structuring element. The figure above demonstrates this effect of dilation. In general, we can say that dilation expands the image by the shape of B.

In this activity we use images of square (50×50) , a triangle (base = 50 , height= 30), a circle radius 25), a hollow square (60×60, edges are 4 pixels thick), and a plus sign (8 pixels thick and 50 pixels long for each line). The structuring elements used are 4 x 4 ones, 2 x 4 ones, 4 x 2 ones and a cross. Below are the result of my predictions and the corresponding dilated images.

my predictions were mostly correct, except for the dilation of a cross using a cross as structuring element. I failed to account for the corners of the cross.
Results below can be summarize in general as follows:

4x4 ones just increase the size uniformly adding 3 pixels

2x4 makes the image wider (horizontally) by about 3 pixels

4x2 makes the image longer (vertically) by about 3 pixels

cross adds 2 pixels on all sides except on corners

Or as a rule of thumb, dilation expands the object similar to the shape of the structuring element (see definition above).

Dilation of different images (leftmost) by different structuring elements (columns 2 to 5)

EROSION
From our activity manual, "The erosion of A by B is the set of all z's such that B translated by z is contained in A."

Similar to dilation, the effect of erosion is to reduce the image by the shape of B.

Summarize below are my predictions and the result of my simulation.

I failed to predict erosion of the hollow figure using a cross (i.e., presence of 4 dots) and the erosion of the cross using a cross as structuring elements.
similar to the dilation, the effect of erosion can be thought of as eroding the shape of the element by the structuring element.

Erosion of different images (leftmost) by different structuring elements (columns 2 to 5)

I give myself a grade of 10 for performing the activity on time and for predicting the results quite well.

I want to acknowledge master, jaya, thirdy for the laughs while programming the code and predicting the results of the erosion and dilation. Miguel for the gimp and the incident about his computer :-p.

References
[1] AP 186 Activity 8 manual
[2] http://en.wikipedia.org/wiki/Set_(mathematics)

Activity 7 Enhancement in the Frequency Domain

2009-07-17T14:21:00.000+08:00

Filtering in Fourier Space

Signal processing as the name implies is an area that deals with analyzing and processing of signals. A major concern in signal processing is the signal to noise ratio. By definition, anything that you do not want is considered noise [1]. One of the methods use in signal processing is filtering in Fourier space. In essence, if we know the characteristic frequency of a signal we can use a filter to remove unwanted frequencies and enhance the signal thereby increasing the signal to noise ratio. In this activity, we demonstrate basic use of filtering in Fourier space to remove unwanted signals.

Activity 7.A: Convolution Theorem

Before we dwell on filtering, we must first familiarize ourselves with the convolution theorem. This is a crucial step since the process of filtering in Fourier space is basically just a convolution. For example, the convolution of a Dirac delta is just equal to the function itself on the location of the Dirac delta. Below are Fourier Transform images of different patterns.

FT of two dots is a sinusoid with a distinct frequency(Dirac delta at f).

FT of circles of different radius. As we increase the radius of the circle the FT decreases in radius. This is expected. i.e., similar to the rect function since FT is anamorphic.

FT of squares of different sides. Again as we increase the length of the side, the width of the FT decreases. In fact, the limit (side -> inf) will result to a dirac delta (point).

FT of gaussian of varying variance. The unique property of a gaussian is that its FT is also a gaussian. Again similar to square and cirlces, as you increase the variance of a gaussian its FT decreases in variance.

Comparing the imaginary and real part of gaussian and inverse gaussian.

click on the image to ZOOM

Gaussian of varying variance with the (a)[leftmost] image of gaussian (b) its real part (c) imag part (d) and FT (e) [bottom] inverse gaussian image (f) real part of inv gaussian (g) imag part of inv gaussian and (h) FFT of inv gaussian.

Observe that the real part of the gaussian and inverse gaussian are the same.

Activity 7.B Fingerprints: Ridge Enhancement
One of the use of filtering in Fourier space is through ridge enhancement of fingerprints. The process can be basically summarize as follows

the FT of the image (i.e, fingerprint - input) was obtained. In order to enhance the visibility of the frequencies, the FT was obtained in log scale.
The obtained FT was analyzed and was used as a template to create a mask that will enhance the frequencies of ridge while decreasing other unwanted frequencies.
It must be noted that use of binary mask is not too adviseable as this will pick only the desired frequencies but not enhance them relative to the other frequencies. Use of gradient mask or grayscale mask is more adviseable. If we know exactly the frequency of the object we want, we can maximize the pixel value of the mask at that location while blurring or decreasing the pixel value of other frequencies that we do not want.
In this example, the author was having a hard time separating the background from the frequency of the ridges. As as a solution, I use a gradient mask on the frequencies of the ridges while blurring it as you go away from the signal. Similarly for the zero order (background), a gradient mask was used which is darker at the center and increasingly bright moving away from the center. This was done because complete removal of the zero-order will not result to a very good reconstruction.

For our purposes, this reconstruction is not the best reconstruction that could be achieve but its a start. The ridges are well defined and darkening of the background provided additional contrast. It is recommended however that further refining be done for better reconstruction.

Activity 7.C Lunar Landing Scanned Pictures: Line Removal

From our previous discussion of FT, we know that the FT of vertical lines are two points symmetric about the center along the horizontal axis. The picture of the lunar above was filled with vertical lines resulting from stitching of several digital "framelets" to obtain a composite image. We can remove this in Fourier space by filtering the image using a mask by filtering the frequencies along the horizontal direction. (see mask used)

Activity 7.D Canvass Weave Modelling and Removal

Similar to ridge enhancement we can remove the pattern created by the canvass weave using filtering in Fourier space. Observe that the canvass weave looks a lot like corrugated sinusoids, thus, we must remove frequencies along the vertical and horizontal side. Also by virtue of rotation, other frequencies might also be evident around the center. The presence of dots in the masked pattern was due to thresholding in scilab (i.e, img(img~=0)=1).

For this activity, I give myself a grade of 10 for obtaining the requirements needed to complete the activity. However, I admit that the reconstructions I presented above are not the best reconstructions. Further research can be devoted to this area to provide better algorithms for filtering and pattern recognition.

I would like to acknowledge Irene and Miguel for all the help they've given while doing this activity. (jokes, company and intellectual discussions)

References
[1] http://en.wikipedia.org/wiki/Noise
[2] Applied Physics 186 Activity 7 Manual.

Activity 6: Properties of the 2D Fourier Transform

2009-07-09T10:47:00.000+08:00

Fourier Theorem states that any function can be expressed as infinite sum of sinusoids of different frequencies. In imaging, this translates to finding the spatial frequencies of an image. In 2D FT, rotation of the image results to rotation of the resulting FT, we will investigate this property in a while but let us first familiarize ourselves with FT of different patterns [1]

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

Below is the FT of different patterns. A quick way to check if the result is correct is to imagine placing a mask(aperture) in front of a light source and observe the image in a screen far away. The resulting FT should have the somewhat the same symmetry as the object.

Observe the difference with the square above, i.e. similar to subtracting terms in your diffraction. In Fresnel diffraction, this amounts to removing zones in your image.

Textbook example of two slits along the x-axis (see Hecht, section diffraction on Fraunhoffer diffraction)

The FT of two dots is a sinousoid. This is textbook example of Young's double slit experiment.

Activity 6.B: Anamorphic Property of the Fourier Transform

Sinusoid of frequency 4 generated using sin(2*pi*f*x).

Sinusoid of frequency 8. The dots are more widely separated. Imagine placing a frequency axis on the y-axis above. The higher the frequency of the sinusoid, the farther apart the dots.

In digital imaging, there is no negative values, adding a bias to an image results to the presence of zero order. This is exactly what happens when we add bias to the image of the sinusoids above.

Sinusoid of frequency 4 with bias. The center dot is the zero-order frequency. In electronics, we can imagine a DC bias or in imaging, the image is offset by the background.

Suppose then that we took a picture of an interferogram, we can obtain the actual frequencies of the object by adding a constant bias to the image (i.e., min value). Or we just obtain the FT and neglect the zero order frequency. However, if the resulting interferogram has non-constant bias then we can add a filter to obtain the desired frequencies. For example, suppose that the image was added with very low frequency sinusoids. We can obtain the frequencies by obtaining the FT of the image and multiply the resulting FT with a high pass filter (i.e, to remove the low frequencies) and then obtain the inverse FT to get the final image of the object. Depending on the bias, we can adjust the filter to obtain the desired frequencies.

In 2D FT, rotating the sinusoids result to a rotated FT.

Rotated sinusoid with its corresponding FT. Observe that the resulting FT is also rotated in the same direction as the object.

Fourier transform of combination of sinusoids. Left: multiplication of two corrugated roofs, in the X and Y direction and right: FT of the object, 4 dots of the same distance signifying same frequency in the X and Y direction.

We added a rotated sinusoid to the corrugated image above using different frequency (f=16 and rotated for 90 degrees). Observe the presence of the 4 dots similar above and the addition of two dots.

We added a different rotated sinusoid using different frequency (f=12 and rotated for 30 degrees).

Added rotated sinusoid of f=8 and rotated for 60 degrees.

Combination of all the above rotated sinusoids. Here we can see that the result is just the addition of all the FT of the different rotated sinusoids. This is because FT is a linear transform, i.e.,

let X=A+B+C then F{X} = F{A} + F{B} + F{C}

Here, we investigated the different properties of the 2D Fourier Transform, in particular

rotation of the object results to a rotated FT
FT is a linear transform
FT obtains the spatial frequencies of the image

For this activity, I give myself a grade of 10 for doing the required objectives.

I would like to acknowledge miguel and martin for useful discussions and master for knowing the spelling of Hecht =).

References:
[1] Activity 6 Manual

Activity 5: Fourier Transform Model of Image Formation

2009-07-07T10:29:00.000+08:00

Lens as a Fourier Transform

Activity 5.A Familiarization with discrete FFT
Using the built-in FT of scilab, fft2, we are going to obtain the FFT of different patterns.

Left: image of circle, right: Fourier transform of circle. As expected, the FT is a bright spot most intense in the center, imagine imaging uniform disk of same intensity.

Left: image of letter 'A' and right: Fourier tranform of circle. The resulting FT has higher frequency components as expected from the shape of letter A.

Activity 5.B Convolution: Simulation of an imaging device
The representation of an image using an imaging device is not a perfect reconstruction but is a "smeared/convolve" image of the original image. Consider the equation below in frequency space

H=FG

Where F is the transfer function of the imaging device and G is the original image. The convolution of f and g results to the image we actually see. In frequency space, by the Convolution theorem, this is just the multiplication of the linear transform (i.e, Laplace of Fourier) of f and g. We simulate an imaging device using an aperture as our transfer function (i.e, this can be modeled as a lens) and convolve to an image of our choice. Below is the image of letters "VIP" which we are going to image using different sets of apertures.

Original image that we will image.

Left: Convolution of the original image using the aperture on the right. Observe that the resulting image is blurry (i.e., diffraction effect). We can imagine in our real life model, i.e lens, that the aperture did not collect sufficient amount of photons resulting to blurry image or the other way to put it is that the NA (numerical aperture) of the lens is very low.

Left: Convolution of the FT of the aperture in the right with the original image. Observe the significant change in resolution when we increase the size of the aperture. (increase NA)

Left: Convolution of the FT of the aperture in the right with the original image. Further increase in the aperture results to better resolution.

Left: Convolution of the FT of the aperture in the right with the original image. The largest apertuer that we can construct. Although observe that even with the highest NA, we cannot replicate the original image (i.e, presence of blurring at the edges) this is because in the mathematics of Fourier Analysis, in order to perfectly reconstruct we must have infinite coefficients to represent a function, but in our case, we are limited to discrete frequencies.

Activity 5.C: Template Matching using correlation
Correlation is basically finding the overlap of two functions. Thus, it can be used for pattern matching.

"The rain...": Image that we will correlate with the image with letter 'A'. The resulting image results with the all the letter 'A' in the image to have high intensity (appear white in the result).

Activity 6.D: Edge Detection
Using similar concept as above, we can correlate an image with different patterns and thus highlighting those parts that we need depending on the pattern.

Edge detection using the pattern in the right. As observed, the pattern can be deduced to detect all point in the edges.

Edge detection using a vertical pattern, as expected, all the vertical parts of the letter VIP are highlighted.

Edge detection using horizontal pattern. The horizontal edges are highlighted.
It must be noted that the sum of the matrix in the pattern is zero.

For this activity, I give myself a grade of 10 for performing all the requirements needed.

Activity 4: Enhancement by Histogram Manipulation

2009-07-01T21:08:00.001+08:00

By looking at the histogram of an image, one can can determine wether the image is too bright or too dark, has poor or good contrast, i.e, if it has good dynamic range. Usually, we want an image to span all the gray levels resulting to a high contrast image. By manipulating the histogram of an image we can improve the quality of the image highlighting or enhancing certain features, or depending on the user, can mimick other imaging system[1].

The histogram of an image when normalized is equal to graylevel probability distribution function (PDF). We can alter the PDF of our image to the desired PDF by backprojecting the grayscale level values of the desired cumulative distribution function (CDF) to the grayscale values of the original CDF.

The steps in overview are as follows:
1. Per pixel in the image having grayscale r, find its CDF value T(r).
2. Find the value of T(r) in the y-axis of the desired CDF G(z).
3. Find the z value of this G(z).
4. Replace the image with grayscale value r to grayscale value z.
This process is called HISTOGRAM EQUALIZATION.

I have provided two methods in solving the problem, the difference is essentially in obtaining the PDF and CDF of the image. In the first solution, we use binning (middle of bin, i.e, 0.5, 1.5...254.5) in obtaining the histogram of the image. In the second solution, we use the function tabul of scilab which tabulates the frequency of distinct grayscale level values.

1st solution:

directory = 'C:\Documents and Settings\Orly\Desktop\AP186\Activity 4';
chdir(directory);
//load image in grayscale
img1=gray_imread('img2.png');
img2=img1; //holder
//obtain the histogram
bits=255;
bins = 0.5:bits; bins = bins./bits; //bins for [0,1] image.
for i=1:bits,
frequency(i) = sum(img1>=(bins(i)-1/(2*bits))&img1<(bins(i)+1/(2*bits))); end; frequency = frequency./max(cumsum(frequency)); f1=scf(1); plot(bins.*bits, frequency); //normalized histogram f2=scf(2); cdf = cumsum(frequency); //cumulative distribution function plot(bins.*bits,cdf); //desired cdf is a straight line, can be change descdf = (0:1.0/(length(cdf)-1):1); f3=scf(3); plot(bins.*bits, descdf); //mapping for i=1:bits-1, index = find(descdf>=cdf(i) & descdf=1 then //satisfied condtion above
//get first index
imgval = bins(index(1)); //desired pixel value
imgcur = bins(i); //current pixel value of image
img2(img2>=(imgcur-1/(2*bits))&img2<(imgcur+1/(2*bits))) = imgval; end, end; //for comparison of images f4=scf(4); imshow(img1); f5=scf(5); imshow(img2); //check the resulting CDF and PDF of the equalized image for i=1:bits, frequency2(i) = sum(img2>=(bins(i)-1/(2*bits))&img2<(bins(i)+1/(2*bits))); end; f6=scf(6); plot(bins.*bits, frequency2./max(cumsum(frequency2))); f7=scf(7); plot(bins.*bits, cumsum(frequency2)./max(cumsum(frequency2)))

The image below is the original image which we will use to demonstrate contrast enhancement.

Poor contrast imageHistogram of the original image, min = 0.5, max =181.5

CDF of the original image.
Desired CDF which is just the CDF of a uniform distribution

Using the above code, the resulting image and histograms is as follows.

Resulting image: Note that in our process, we span all the gray levels and redistribute them in the image in a uniform process. The midtones are enhanced but the problem is the image got saturated on parts with high gray level and too much shadows on parts with low gray level value. Also, the result is a bit noisy compared to the original.

Histogram of the resulting image, min = 0.5, max =253.5. The contrast of the resulting image base on the minimum and maximum gray values increased dramatically.

The CDF of the resulting image follows the shape of the desired CDF above.

Main problem encountered with the above solution is that for those images with limited dynamic range (i.e, 100-150), the resulting image would be too saturated for those pixels with higher pixel value (130-150) and too dark for those pixels with lower pixel value (100-120) because essentially, we are uniformly stretching the image into all the gray levels (0-255). The resulting image might not be too pleasing but it would have better dynamic range.

Another solution for this is to uniformly distribute the image (i.e, linear CDF) in the range of the original histogram. (i.e., from 100-150). In essence, we LOCALIZED the histogram equalization.

2nd Solution
directory = 'C:\Documents and Settings\UP DILIMAN\Desktop\ORLY\Acads\AP186\Activity 4';
chdir(directory);

//load image in grayscale
img1=gray_imread('img2.png');
info = imfinfo('img2.png');
img2=img1; //holder
bits=(2^info.Depth)-1;
//find frequency of distinct values
hist=tabul(img1,'i'); //arrange in increasing order
frequency = hist(:,2)./max(cumsum(hist(:,2))); //normalized frequency
bins=hist(:,1);

//displaying histogram
histogram=fhistogram(img1, info.Depth);
frequency2 = histogram(:,2)./max(cumsum(histogram(:,2)));
f1=scf(1);
plot(histogram(:,1).*bits, frequency2);

//find the minimum and maximum grayscale values
index=find(frequency~=0);
minB = bins(index(1))*255;
maxB = bins(index(length(index)))*255;

f2=scf(2);
cdf = cumsum(frequency);
plot(bins.*bits,cdf);
//desired cdf is a straight line, can be change
descdf = (min(cdf):(max(cdf)-min(cdf))/(length(cdf)-1):max(cdf));
desbins = (min(bins):(max(bins)-min(bins))/(length(bins)-1):max(bins));
f3=scf(3);
plot(desbins.*255, descdf);

//mapping
for i=1:length(bins)-1,
//returns the index in descdf where cdf(i) belongs
index = find(descdf>=cdf(i)); //returns a list of indices
if sum(index)>=1 then
//get first index
imgval = desbins(index(1)); //desired pixel value
imgcur = bins(i); //current pixel value
img2(img2==imgcur) = imgval;
else
x_message('cdf index is' + string(i));
end,
end;

f4=scf(4);
imshow(img1);
f5=scf(5);
imshow(img2);
f6=scf(6);
histogram2=fhistogram(img2, info.Depth);
frequency3 = histogram2(:,2)./max(cumsum(histogram2(:,2)));
plot(histogram2(:,1).*bits, frequency3);
f7=scf(7);
plot(histogram2(:,1).*bits, cumsum(frequency3));

code of filename: fhistogram.sce

function histogram = fhistogram(img, depth)
bits=(2^depth)-1;
histogram = zeros(bits,2);
bins = 0.5:bits; bins = bins./bits; //bins for [0,1] image.
for i=1:bits,
frequency(i) = sum(img>=(bins(i)-1/(2*bits))&img<(bins(i)+1/(2*bits))); end; histogram(:,1)=bins'; histogram(:,2)=frequency; endfunction
Using tabul, the bins are the distinct grayscale values and the frequency is just the occurrence of that distinct grayscale value.

Resulting image after mapping. Note that the image has better contrast compared to the reconstruction using the first method. The advantage of this reconstruction is that, the image would not be too saturated on some areas and too dark on other areas because we uniformly distributed the image in the range of the min and max grayscale values.

Resulting histogram, min=0, max=182

Resulting CDF, note that it is almost similar to the desired CDF.

By defining another set of bins ('desbins') that are also uniform in the min and max grayscale level values, we have essentially equally redistributed the histogram of the image.

Using other images:

Original: left, Contrast enhanced: right. The reconstructed image has better contrast compared to the original image. Image obtained from source [2].
Original: left, Reconstructed: right. Note that the effect of the histogram equalization is not quite visible. This is because the dynamic range of the image is limited as can be seen in the histogram and CDF of the reconstructed image. Image was obtained from source [3].
Note that in this histogram and CDF, we can see clearly how the method works. The histogram was equally distributed in the range of the histogram (min - max) hence the resulting CDF is linear in that range only.

USING OTHER CDF:
As a proof of principle, we use an exponential CDF to demonstrate what happens on the image when using nonlinear CDF. Depending on the user, we can opt to highlight the darker areas (low grayscale values), the whiter areas (high grayscale values) or the midtones (middle grayscale values). When highlighting the darker areas, we can use logarithmic CDF, when highlighting whiter areas, we can use exponential CDF and when highlighting the midtones we can use S-curve CDF.

The code above can be altered by defining the appropriate desired CDF. We just inserted the following code:
desbins = (min(bins):(max(bins)-min(bins))/(length(bins)-1):max(bins)); sigma =2; //strength of exponential descdf = exp(sigma*desbins);
descdf = (descdf - descdf(1))./max(descdf-descdf(1));

Resulting image using exponential CDF. Note that the resulting image highlighted the higher grayscale values (whiter areas) as expected from the shape of an exponential curve.

Left: Desired CDF, Right: CDF of the resulting image.

We have provided a solution for contrast enhancing an image by histogram manipulation. Our solution involves locally equalizing the histogram of an image using the available min and max pixel values of the image.

It must be noted that there is no single CDF that can enhance the contrast of every image. Each CDF has its own limitations and use. Care must be taken when contrast enhancing an image and knowledge of how each CDF alters the image must be understood.

Please email me (orly.tarun@gmail.com) if you plan to use the code above! Please leave your comments and suggestions. =)

For this activity, I give myself a grade of 10 for performing acceptable contrast enhancement algorithms.

References:
[1] Activity 4 Manual.
[2] http://myweb.lsbu.ac.uk/dirt/museum/margaret/08--452-1000120.jpg
[3] http://www.andrew.cmu.edu/user/timothyz/hw3/low_contrast.jpg

Activity 3: Image types and basic image enhancement

2009-06-23T09:27:00.001+08:00

Digitize image comes in four different types. They are,

Binary image which are either 1 or 0, black or white. For an image of size 128 x 128 pixels, each pixel occupies 1 bit, so 128 x 128 = 2048bytes or 2KB.
Grayscale image which contain 8 bits of information. That is, their pixel value can span the range between 0-255. Since it contains 8 bits of information, each pixel occupies 1 byte, for an image of 128x128 pixels = 16384 bytes or 16KB.
Truecolor image is equivalent to 3 grayscale images called channels. Truecolor images are also know as RGB because each R, G, B represent the 3 different channels of the image (red, green and blue). Per pixel, there are 256^3 possible combination of information. Most colored photos, i.e, taken using digital cameras are RGB.
Indexed images are color images that contain index numbers mapped to a color map. Thus, they contain two data, the colormap and image. They are generally smaller in file size compared to RGB but less color information because the maximum colors it can contain are indicated in the colormap.

Below are examples of different types of images.
JPEG images of Lena obtained from this source [1]. Note of the difference in file size between the different image file types: Truecolor (RgB), grayscale and binary (black and white).

File Properties: RGB (72.4KB) Grayscale (26KB) Black and white (5KB) Index(31.6KB)
*(http://www.ee.columbia.edu/~cylin/course/mss/MSS_hw1.html)

In the next part, we are going to apply what we learned in activity 2 by obtaining the area of a scanned image. I would like to acknowledge Mark Jayson 'master' Villangca for lending me his scanned image for this activity. In order to check if our area calculation is correct, a ruler was placed side by side with the image. We can then obtain a pixel:physical unit calibration to obtain the true area of the image.

FileName: master.png

FileSize: 208681
Format: PNG
Width: 905
Height: 446
Depth: 8
StorageType: truecolor
NumberOfColors: 0
ResolutionUnit: centimeter
XResolution: 59.050000
YResolution: 59.050000

In order to obtain the area of the image above, we need to binarize the image using im2bw. In doing this, we need to select the correct threshold using the histogram of the image.
Histogram of the above image using my own code for histogram below,
function histogram = fhistogram(img, depth)
bits=(2^depth)-1;
histogram = zeros(bits,2);
bins = 0.5:bits; bins = bins./bits; //bins for [0,1] image.
for i=1:bits,
frequency(i) = sum(img>=(bins(i)-1/(2*bits))&img<(bins(i)+1/(2*bits))); end; histogram(:,1)=bins'; histogram(:,2)=frequency; endfunction

From the histogram we can choose a suitable threshold (i.e, ~ pixel value of 100). Using im2bw, we can then set the threshold (since 0-1, divide 100 by 255) and specify the ROI in the image by defining another matrix. In the above image, we opted to remove the ruler in the thresholding and focused on the image.

Resulting image after converting to binary and applying a region of interest
After using the contour follow, we obtain the contour of the image.

Using the algorithm in activity 2, we obtain the following results
Area of master = 16576 pixels.
Pixel count = 16894 pixels.
Observe the converted binary image has extra white pixels outside the contour 'master'.

In order to calibrate this to physical units we can obtain the pixel:physical ratio using paint. (i.e., we obtained 1cm:30pixels ~ ). Thus the area of the contour in physical units is
Area = 552.533 cm^2.

For this activity I give myself a grade of 9 for finishing the activity although a bit late and additional bonus for providing my own code for the histogram.

I would like to acknowledge Mark Jayson Villangca for lending me the scanned image. Jica and Miguel for useful discussions.

References:
[1] Activity 3 manual.

[2] http://www.ee.columbia.edu/~cylin/course/mss/MSS_hw1.html

Activity 2: Area Estimation of Images with Defined Edges

2009-06-23T09:23:00.000+08:00

Area estimation finds many applications in the field of material science, biology, geology and the natural science in general . Biologists are especially interested in the difference in size between cancer cells and healthy cells. In material science, area measurement of leads and semiconductor materials among others is an effective tool for processing and failure analysis. In geology, remote sensing is an active area of research used widely in environmental and marine science researches [1].

The method below discusses image area measurement using Green's Theorem. It must be noted that before using this method, there are several limitations that must be considered. We will discuss these limitations at the end of the article.

We start with the following.
Below are sample images (256 x 256 ) of triangle, square, circle and several irregular polygons which we will use in our area estimation method.

Circle of diameter=125px, Square of side=130px, triangle of base=height=150px. These measurements are obtained from MS Paint while creating the images. It must be noted that these measurements are just guidelines for checking our results and must not be used as the theoretical area of the image. In creating these images, the border already occupied several pixels, so the area of the shape (white) does not necessarily reflect the analytically expected result. (i.e, Area_square = 130*130)

Pixel Count: 16129
Computed Area: 16129

Pixel Count: 11644
Computed Area: 11644

Pixel Count: 10731
Computed Area: 10731

Pixel Count: 17522
Computed Area: 17522

Pixel Count: 14058
Computed Area: 14056

Pixel Count: 14913
Computed Area: 7414

In order to compute the area of the above images, we use the SIP toolbox of scilab-4.1.2 using the command follow. Excerpt from the description of follow command from scilab:

"follow - A contour follower. Function follow extracts parametric contours of binary objects. This is useful for further extracting object features such as curvature and bending energy."

Note that before we use the command follow, our images must be first binarize. This can be done in scilab using the command im2bw. The command follow returns the [x,y] coordinates of the contour. We can then use this contour coordinates to compute the area using Green's Theorem found below. For more detailed discussion of area measurement using Green's Theorem please follow this source [2] (http://www.attewode.com/Calculus/AreaMeasurement/area.htm).

Since the image is binarized, we can compute the thoretical area by summing the pixels of the image (white~1 and black~0). We designated the theoretical area as "pixel count" in the above images and "computed area" when using Green's Theorem.

Below is the code in scilab used to obtain the pixel count and the computed area:

//load image
directory='C:\Documents and Settings\2004-29578\Desktop\Activities\Activity 2\'
filename='quadrilateral_130side.jpg';
img1= imread(directory+filename);

//binarize image
img2=im2bw(img1, 0.5);

//get the contour using follow
[x,y]=follow(img2);
//get the dimensions of x and y
[M,N]=size(x)
//complete the contour setting the last pixel=first pixel
x(M+1)=x(1);
y(M+1)=y(1);

plot(x,y);

//get the area using Greens Theorem. we use round since the area (number of pixels (1 or 0)) can only take //integer values.
Area = round(0.5*sum((x(1:M).*y(2:(M+1)))-(y(1:M).*x(2:(M+1)))))+ round((M+1)/2);
//pixel count
TArea = sum(img2);

//imshow to scale to grascale, i.e., highest pixel becomes white, lowest pixel becomes black
//binarize image using im2bw
//increase memory for scilab use, stacksize.

Note that in obtaining the computed area, we added a correcting factor "round((M+1)/2)" in our code where M is the length of the x or y vector.

Results and Limitations
(1) For the regular polygons and the first irregular polygon, the pixel count and the computed area are the same, meaning that the algorithm successfully estimated the area of the images.
(2) For the next irregular polygon (figure batman) there is a discrepancy in the pixel count and the computed area.

In the figure above, note the extra 2 pixels in 'batman's" horn which accounted for the discrepancy. This error is due to the command "follow" of scilab which assumes only ONE contour per image. In order to check this, we made two contours in the next image and as expected there is a huge discrepancy between the pixel count and the computed area.

(3) The algorithm and theory for computing the area is pretty solid (Green's Theorem) hence most of our errors should come from the input and processing of the input. That is, our limitations will most probably due to the command follow and im2bw of scilab.

Overlayed image of the binared square and its contour using the command follow. Note that there is an offset in the contour compared to the image of ~2 pixels. Similar observation is found in the next image.
Note the offset between the contour and the image.

(4) The command im2bw is similar to thresholding and thus depending on your threshold, the computed area will change. This effect would most probably depend on the user setting the threshold and making measurements.

(5) It is expected that in large scale images (2048x2048) the computed area will approach the analytical area. This would be very useful in applications like remote sensing and geological survey.

For this activity, I give myself a grade of 10 for exact area measurement as proven by the sets of images used.

References:

[1] Activity 2 Handout
[2] http://www.attewode.com/Calculus/AreaMeasurement/area.htm

Activity 1: Digital Scanning

2009-06-16T10:38:00.000+08:00

excerpt from activity 1 handout...
"In this activity, we use ratio and proportion to find the numerical values of a digitally scanned hand-drawn plot. For now, you will need the following software, Paint for finding the pixel locations on the graph, and Excel of Open Office Spreadsheet to tabulate the values of the graph."

What I did.
1. Photocopied a plot from the College of Science Libray from the following journal.
The Renormalization Group and Cricitical Phenomena I: Renormalization Group and the Kardanoff Scaling Picture by Kenneth Wilson from the Laboratory of Nuclear Studies, Cornell University, Ithaca New York, 14850 on 1974.

The photocopied image was then digitally scanned. (see picture below)
For convenience, I cropped the image at the origin resulting to the image below.
The cropped image was then opened in MS Paint.
Producing the Plot in Open Office:
Using MS Paint, the pixel location of the origin was recorded. Note that in paint, the (0,0) mark is at the top-left of the document. In this activity, the pixel location of the origin of the cropped image is at (220,1053).

Using Open Office, I successively obtained the pixel locations of the graph at 50 pixel interval in the x-axis amounting to 33 data points. Using the pixel locations of the origin, I calibrated the pixel location obtained previously by noting that the top-left part of Paint is the (0,0) mark using the following formula:
Xexcel=Xpixel-220,
Yexcel=1053-Ypixel (see picture below for sample data points)

Next, I obtained the number of pixels/unit in the graph and obtained the following:
Pixels/unit in the x-axis = 648
Pixels/unit in the y-axis = 1276
Using this calibration I divided all the values of the Xexcel by 648 and all values of Yexcel by 1276. This will give the calibrated values of the graph in physical units.

In order to check if the reconstructed image and the original image are perfectly matched, the two images were overlayed. This can be easily done in open office by importing the scanned image as graphic. (Format->graphic->bitmap->import the image). After importing the image, click on the reconstructed graph and set the Area to the imported image (click on the graph->object properties->area->select bitmap->select your imported image).

At first, the images will not be perfectly matched. In order to overlay it perfectly, you must obtain the maximum locations of the grid in the x and y axis. Use this and modify the axes of the graph and set it to the obtained maximum.

The image below is the reconstructed image overlayed with the scanned image.
Observe that the image grids in the scanned image and the reconstructed image are perfectly matched. Furthermore, the graph is well-matched with no visible offset. For this activity, I give myself a grade of 10 for perfectly reconstructing the scanned image.