Scilab Playground

Sayonara AP 186!

2011-10-14T03:41:00.000+08:00

Created using Wordle.

Basic Video Processing

2011-10-14T03:21:00.000+08:00

In this activity, we integrate all the techniques we have learned so far in image processing in order to perform video processing. A video can be considered as a series of images shown in rapid succession. The number of frames captured within a second depends on the frame rate of the camera used. For a camera with frame rate of 30fps (frames per second), 0.03 s is the interval between successive frames. Thus, image processing can be applied on a video in order to extract necessary information from it. What is interesting in such investigation is the appearance of a variable t, i.e., time, which is important in the dynamics of a certain system.

Shown below is a video of a Christmas ball swinging in pendulum-like motion. We aim to track the position (x, y) of the ball (in pixel coordinates) in the whole duration of motion. This video contains 30 frames per second (30 fps).

The image sequences are obtained using Virtual Dub. Each image is then processed in Scilab. Some of the results after processing the images are shown below.

The ball is segmented from the background by finding the value in the histogram of each image that separates the background from region of interest. It is found by looking for the dip between two peaks which correspond to the region of interest and background. In this case, it is 0.2. The results of thresholding are shown above for frame nos. 5, 20, 50 and 35 in clockwise manner. What is done next is to find for the exact coordinates (x, y) of the center of the ball in each frame. In finding this, the center of mass of all pixels belonging to the ball is calculated. The image below shows the location of the center of mass throughout the duration of motion. It somehow tells us about the trajectory of the ball.

Take note that the positions are in pixel coordinates wherein the origin is located at the uppermost left hand corner. Our world can be compared to the 4th quadrant in Cartesian coordinates system. Symmetry with respect to the y-axis can be observed in the figure above. The right hand of the parabolic trajectory being longer is due to the ball being released at a high location in the beginning of its motion.

In systems like this, it is important or interesting to look on the evolution of a certain variable with time t. In the video shown, damped oscillation occurs. To characterize such behavior, we look on the plot of x/y versus frame no. (30 frames in a second) wherein 1 frame = 0.03s.

Plot of x versus frame no. (t = frame no. X 0.03s)

Plot of y versus frame no. (t = frame no. X 0.03s)

The figures show that indeed the system investigated displays damped oscillation. We know that after sometime t, the oscillation of the ball will stop due to some factors (e.g. frictional forces).

For this activity, I am giving myself a grade of 10. This is an easy yet interesting activity. Most of the time, information available are stored in videos.It is then a good thing to learn how to perform video processing.

Sayonara, AP 186!

Neural Networks

2011-09-28T14:05:00.001+08:00

In this activity, we try using the artificial neural network or ANN toolbox for Scilab in object classification. A neural network is a computational model of how the neurons in our brain work. This is an alternative to linear discriminant analysis or LDA in pattern recognition. In neural network, a pattern is learned through example. The learning process may take some time but once a pattern is learned, a faster recognition process is expected.

Below is an example of an artificial neuron modeled by McCulloch and Pitts in 1943. A neuron receives weighted inputs w_ix_i (where i = 1, 2, 3 for the example shown below) from other neurons. These weighted inputs are then added to yield a and then acts on an activation function g. The output z is then fired to other neurons.

An artificial neuron.

A neural network is formed by connecting many neurons. An example of this is shown below. A typical neural network consists of 3 layers: input layer, hidden layer and output layer. The input layer can be a set of features extracted from the objects to be classified. The hidden layer then acts on such set and then passes the result to the output layer.

A neural network.

We are to utilize the ANN toolbox in object classification. There are two classes of objects, A and B, to be classified based on two features extracted. The code used is shown below. Take note that the input training matrix x consists of 12 rows with 2 columns since there are 12 objects to be classified based on 2 features. The desired output t represents the class of the object classified such that a value equal to 1 means class A and a value equal to 0 means class B. Take note that the values for both sets of features are normalized.

The parameters learning rate and training cycles can be tweaked. The results for different combinations of values of learning rate and training cycles.

learning rate = [0.1, 0], T = 400

learning rate = [2.5, 0], T = 400

learning rate = [0.1, 0], T = 1000

learning rate = [2.5, 0], T = 1000

It is obvious from the results shown that as the learning rate is increased, the more closer the obtained values for the output are. Also, increasing the number of training cycles results to values for the output that are closer to the desired ones. It turns out that the predictions neural network made are 100% accurate. Using ANN toolbox in classification is easier compared to using LDA for many reasons.

This is again an easy activity. Thanks to J. Tugaff for the sample code which made the activity much easier. I am giving myself a grade of 10 for this one.

Probabilistic Classification

2011-09-20T21:58:00.001+08:00

This is a continuation of the previous activity. Instead of using minimum distance classification in pattern recognition, another technique called linear discriminant analysis or LDA is used. Given a pattern x, the probability that it belongs to class ω_k is given by p(ω_k| x). In the case wherein a classifier wrongly identifies x to be part of class ω_j but actually belongs to ω_k, a penalty loss λ_jk is acquired for such mistake. A loss matrix L stores all these penalty terms λ_jk at their corresponding indices (k, j). The aim here is to minimize the risk or loss in classification given by

In LDA, the main assumption is that classes or groups of objects are linearly separable, i.e., classes are separable by linear combination of the features of these objects. Also, what we are interested in is to find the class where a test object belongs to which is represented by a so-called dependent variable y. On the other hand, the independent variable x pertains to any quantifiable or measurable features of the object.

In order to minimize the total error of classification or TEC, an object is assigned to a class with highest conditional probability. To do this, there is a need to find for the conditional probability P(i|x) that an object having a set of features x belongs to i. But such conditional probability is difficult to find. What is readily available instead is the conditional probability P(x|i) that describes the chance of obtaining a set of features x given that an object belongs to i. Luckily, Bayes theorem relates these two conditional probabilities and is given below.

P(i) is the probability about the group i. In our case, P(i) is based on the number of samples for each group. However, it is not advisable to use Bayes rule directly for some reasons. But if we assume that each group has multivariate normal distribution and all groups have the same covariance matrix, a so-called linear discriminant analysis formula shown below can be used in the classification.

We assign an object with feature vector x_kto class i with maximum f_i. The table below is the result after applying LDA on two classes of objects.

It is clear from the table above that the results of classification using LDA is 100% accurate, i.e., it was able to classify objects correctly.The training sets used are classes of Piattos chips and 1-peso coins.

This activity is an easy one since no more image processing were to be done. Just plain computations. And I enjoyed doing this. I learned a lot from these activities on pattern recognition. And I am planning right now to have my project based on the concepts on pattern recognition. I am, therefore, giving myself a grade of 10 for a job well done. :-)

Pattern Recognition

2011-09-18T22:20:00.001+08:00

Humans have the ability to sort and label things according to features based on previously observed examples. We learned from experience that a circle is different from a square and that red is not equal to green. Computer vision aims to mimic such ability of the human brain in pattern recognition. A set of features defines a pattern. This set of features contains information such as size, color, shape, etc. It can be represented as a vector and is called a feature vector. A certain feature vector is unique for every class of objects. The tasks in pattern recognition is to determine the features for all classes to be classified and to come up with a classifier based on these features collected.

I assembled 3 classes of objects. Each class is divided into two sets namely training set and test set. The objects are shown below.

Left: Training sets. Right: Test sets.

The features are extracted from these classes, i.e., both training sets and test sets, to form unique feature vectors for the 3 classes. Each class has a set of training set feature vector and test set feature vector. The training set feature vectors are then used to compute for the class representatives m_j. In minimum distance classification, the class representative m_j of class ω_j is calculated from the ensemble of training set feature vectors using the equation below and is also called the mean feature vector. This is done for each class of objects. N_j is the total number of objects contained in the training set of class ω_j. The calculated class representatives are shown in the table below as highlighted values wherein each color represents a certain class of objects. Note that the mean feature vector contains only two elements, i.e., mean area and mean perimeter. Also shown below are the calculated mean test set feature vector for each class.

The figure below is a plot of Perimeter versus Area in pixels. 3 clusters are evident in the plot below which correspond to the 3 classes of objects investigated, i.e., 1 peso coin, Piattos Chips and 5 peso coin.

We are now ready to classify the calculated test set feature vectors. This is done by finding the class wherein the test set feature vector x is closest to. The equation below is utilized for this task. The membership of a test set is determined by finding the class wherein d is maximum.

The table below shows the calculated d for different combinations of test sets and training sets. The values in red are where d is maximum. The classifier outputs the correct membership for each test set as shown below.

After applying minimum distance classification on all test sets, all objects for each test set are classified correctly which gives us a 100% accuracy for this case. Such accuracy obtained is I think due to the well clustered features we see earlier.

There are a lot of interesting features that were not extracted and investigated in this activity. And color is one these features. I am giving myself a grade of 10 for this activity since I was able to output the necessary results.

Color Image Segmentation

2011-09-09T01:05:00.000+08:00

The main goal of this activity is to segment a region of interest ROI from an image based on the color of the ROI. The easiest way to segment a region of interest is, of course, by thresholding the grayscale image. This technique, unfortunately, does not work at all times. In the case wherein the ROI cannot be segmented through thresholding alone, color image segmentation is the solution. Here, the color of the region to be segmented is used as a basis in the segmentation. But note that an image of a 3D object can have shading variations. It is then important to represent the color space in such a way that it separates brightness and chromaticity (pure color). A very useful color space is the normalized chromaticity coordinates or NCC. Per pixel, let I = R + G + B wherein R, G and B represent the relative intensity of the red, green and blue channels, respectively. R, G and B are normalized such that

Take note that r + g + b = 1 which then allows us to express b in terms of r and g, i.e., b = 1 - r - g. This gives us the normalized chromaticity space shown below which is a 2D color space. In this representation, the r-g plane gives us the chromatic information. A pixel in the image can be represented as a point in the NCC space. A region of interest ROI, for example, when represented in the NCC space, occupies a small blob positioned somewhere in space depending on the color of the ROI.

Normalized chromaticity space with the x and y axes representing r and g, respectively.

In color segmentation, determining whether a pixel belongs to the ROI can be done by looking on the probability of the pixel to belong on the color distribution of the region of interest. This probability distribution function PDF can be obtained by finding the histogram of a portion in the ROI and normalizing it by the total number of pixels of the subregion selected. There are two ways to segment ROI: parametric segmentation and non-parametric segmentation. Shown below is the image to be segmented together with the cropped ROI.

a) The image to be segmented with the yellow object as the region of interest. b) A cropped portion from the ROI.

Parametric Segmentation

In parametric segmentation, a Gaussian distribution is assumed independently along r and g and a joint probability p(r)p(g) function is used to determine whether a pixel belongs to the region of interest.

The histogram of the subregion together with the normalized chromaticity space is shown below.

Left: The histogram of the subregion. Right: The normalized chromaticity space

Note that the histogram displays a peak which corresponds to the color yellow in the normalized chromaticity space. Shown below is the result using parametric segmentation.

The results using parametric segmentation with the left image in binary.

Non-parametric segmentation
In non-parametric segmentation, the 2D histogram of the subregion chosen is used in tagging whether a pixel belongs to the ROI. Moreover, histogram back projection is done in assigning a pixel a value in chromaticity space based on its value on the color histogram.

Left: The histogram of the subregion selected. Right: The normalized chromaticity space.

It is obvious that the peak observed in the histogram of the subregion chosen corresponds to the yellow region in the normalized chromaticity space. Having found the histogram of the subregion allows us now to perform histogram backprojection in order to segment the ROI. And the result for this method is shown below.

The results using non-parametric segmentation with the left image in binary.

The results obtained for both techniques suggest that parametric segmentation is better than non-parametric segmentation. The use of the Gaussian fit in the parametric segmentation makes the quality of the output better than that of non-parametric. The output obtained for this case reflects the smoothness of the Gaussian curve used.

This was another difficult activity for me. But I still enjoyed doing this activity. Everything just paid off when I saw the results of the segmentation. I am giving myself a grade of 10 for this since I was able to show results for both segmentation techniques.

Image Compression

2011-09-02T00:13:00.000+08:00

Principal component analysis or PCA is utilized in this activity as image compression technique. This activity has the goal to demonstrate how image compression works. Nowadays, storing images of large file sizes is no longer a problem with the emergence and availability of cheap flash disks, and probably hard disks. But before we experienced such comfort, there was a time that storage was really a problem. This is where image compression became useful. Several tutorials on PCA are available over the internet. The technique itself deals with a lot of linear algebra but basically PCA dwells on the idea that a given image can be represented as a summation of weighted basis images. PCA, as a compression technique, uses a given number of eigenfunctions (eigenimages) which depends mainly on the degree of accuracy of the reconstruction. A lesser number of eigenimages used suggests a smaller file size and that quality is sacrificed.

The image below is the grayscale image to be compressed. It has a file size of 735kB. Using the built-in PCA function in Scilab, the shown image is then compressed.

A sample grayscale image to be compressed with file size of 735kB.

But before PCA takes the job, the image to be compressed must be cut into 10x10 blocks! Take note that our image is 2440 x 3260! And that is a lot of square blocks to make. Each subblock is then converted into a column matrix having 100 elements. These are then arranged to form an n x p matrix x where n is the number of blocks and p is the total number of elements in each subblock, i.e., 100. PCA is then applied on the matrix x. The PCA outputs a set of eigenvalues (lambda), together with their corresponding eigenvectors, that are arranged in decreasing order. Each lambda represents the weight of each eigenvector. Given 1000 subblocks, for example, with size 10x10, 100 eigenvectors can be formed with length equal to 1000. Each eigenvector contains the intensity value of the same pixel location for the 1000 subblocks. The eigenvalue (lambda) dictates the number of eigenimages to be used in the reconstruction of the image depending on some threshold to be set [1]. The decrease in the number of eigenvectors (eigenimages) used leads to the decrease in the file size of the image. Below is an example of the compressed image with the file size reduced to 567kB.

Compressed image with the file size reduced to 567kB.

There is an observed difference between the original image and the compressed image in terms of resolution. It appears that image compression using PCA sacrifices the quality of the image to be reconstructed since fewer number of eigenvectors are used to represent the given image. Increasing the number of eigenimages to be used leads to a more accurate reconstruction of the image but does not help in reducing the file size of the image. Below shows the trend observed as the number of eigenimages is increased from 1 to 99.

The total number of eigenimages is 100. And the size of the original image is 735kB.

Figuring out what the outputs of the function PCA mean just consumed much of my time. And do not forget about the complicated PCA routine. Working with matrices is such a mess. But despite these problems, I was able to fulfill the objective of this activity which is to compress a sample image. And thank God for that. And yes, I am giving myself a grade of 10.

References:

[1] Principal component analysis. Retrieved Sepetember 1, 2011 from http://en.wikipedia.org/wiki/Principal_component_analysis.

Preprocessing Text

2011-08-31T22:51:00.001+08:00

This activity aims to integrate and apply all the things we have learned in image processing so far. The first time I had a glimpse on the manual for this activity, I thought it was much easier than the previous activity (click here). I still don't know why it took me more than two meetings to finish it. Compared to the previous activity, this one requires a lot of subjectivity. Yes, it is clear that the main objective is to extract handwritten text from an imaged document. But one's judgement regarding the quality of the output is so subjective that determining when to stop from repeating the process again and again is not definite.

The imaged document used for this activity.

The whole image was not used. Instead, we cropped a small portion containing handwritten text which is the one shown below.

Left: Fourier transform FT of the cropped image on the right. Right: the cropped image from the imaged document which contains handwritten text.

From the imaged document itself, it is very clear that the image is tilted. Looking on its FT is helpful in finding the tilt angle. It was suggested that we use the mogrify function in Scilab. The documentation in Scilab is not helpful so I searched over the net for the right arguments for the function. For some weird reason, Scilab kept on returning error messages for each attempt I made just to rotate the image using mogrify. This just consumed much of my time which made me rotate the image using MATLAB. I know this is a mortal sin but... mogrigy is not really cooperating.

Finding the tilt angle was done manually. Approximately, a tilt angle of 1degree was chosen. The rotated image and its Fourier transform are shown below.

Left: the rotated cropped image. Right: The Fourier transform of the image on the left.

To remove the horizontal lines, a filter must be created such that it covers the vertical line that appears on the FT of the image. Recall that a sinusoid along the x direction produces a Fourier transform that is composed of two dots that is symmetric along the y-axis. The spacing between these observed symmetric peaks is dependent mainly on the frequency of the signal.

Left: the resulting image upon the addition of the filter. Right: the FT of the rotated image with a filter added.

What we have done so far is to remove the horizontal lines in the selected portion of the imaged document. We now proceed to converting the image to binary. The histogram of the said image is shown below together with the resulting binary image using a threshold of 0.51.

Left: histogram of the image. Right: the resulting binary image using a threshold of 0.51.

It is evident on the histogram that a single peak exists. In this case, it is very difficult to determine the perfect threshold value to use in order to separate the background and the region of interest, i.e., the handwritten text. I decided to choose a threshold value that is near 0.5, which I think is a safe value to use.

Now, our goal is to clean the image using the techniques we have learned so far. The figure below summarizes the things I have done on the binary image.

Top row: dilate, erode. Bottom row: close, open.

I tried different morphological operations on the binary image. I am not quite satisfied with the quality of the results I obtained for this part of the activity. Cleaning the image seems to be impossible. First, not only handwritten text was present in the chosen portion but also typewritten ones. Some of these components after the operation appear to be dirt on the processed image. Also, the effect of filtering the vertical line on the FT is that the handwritten letters appear to be sliced horizontally. This hinders us from labeling each handwritten letter in the image. I find it hard to decipher the text written by manual inspection. The letters were written so close together that preprocessing the text was too hard to accomplish. Yes!

Another thing I made was to use the imconv function on the resulting image after a filter was applied. This is helpful in the case of edge detection. Somehow, what I wanted to happen on this is for imconv to detect the edges of the text present on the image, Two types of mask were used, i.e., a cross pattern and a dot pattern.

Left: the result using a cross mask. Right: the result using a dot mask.

Apparently, the results when these images were converted to binary are not okay, i.e., the results obtained earlier using the morphological operations are better.

Below are my results for the part wherein the other instances of the word "DESCRIPTION" are to be found. The very first technique that crossed my mind is template matching. Given a template containing the word "DESCRIPTION" which must at least appear to have a font style and size similar to that of the typewritten text, its correlation with the imaged document will be found using the imcorrcoef function. The positions where the other instances occur result to a higher correlation values. The result is shown below.

Left: the imaged document with the horizontal lines removed by adding a filter. Right: the output after using imcorrcoef.

The output after obtaining the correlation of the template with the image appears to display the maximum correlation values at the locations of the instances of the word "DESCRIPTION". However, there are observed locations where the resulting correlation values appear to be maximum but do not correspond to the positions of the instances to be found. I know that this technique is too qualitative to be used.

I am giving myself a grade of 9 for the lousy results I have. I am not completely contented with the things I have done for this activity. I think I have not fully applied the things that I have learned.

Playing Notes by image processing

2011-08-19T04:10:00.001+08:00

This activity is, I think, the most challenging so far. First, I need to learn how to read a musical score sheet, i.e., from notes to their corresponding tones. I got so frustrated during the 1st day that I only continued doing this a few hours before deadline (that explains why this blog post is late). Okay let us just forget about these rants and proceed to our main concern.

It is only when I totally concentrated on doing this that I realized how cool this activity is. Yes! This is true especially when everything is working already, i.e., from the code to the saved sound file. I chose a song which has a musical score sheet that is composed of only two kinds of note: half note and quarter note. The song is the famous 'Twinkle Twinkle Little Star'. Its score sheet is shown below.

This musical score sheet is downloaded from mamalisa.com.

The very first thing I performed on the digitized score sheet was to cut it into 4 parallel rows. Each strip containing the necessary information would then undergo image processing. Since we aim to play the song, detecting what types of note are present in the image is needed. My idea was to use template matching. There only exist two types of note in the whole score sheet. This is advantageous since only a single template will be used. In my case, I opted to use the quarter note as the template. In template matching, the size of the two images being matched must be the same. The figure shown below is the template used for this activity.

The template, in binary, used for this activity.

The first row in the musical score sheet.

Shown below are the results when the 1st row in the musical score sheet is image processed. The first row image shows the correlation of the binarized template and score sheet. This is obtained by multiplying the FT of the template to the conjugate of the FT of the score sheet. The second row presents the result when the obtained correlation of the two function is converted to binary using a good threshold value which is 0.55, in our case.

Notice that difference between the forms of a quarter note and a half note. Our problem regarding how to detect the right note is solved! To further quantify things, the histogram of the areas is investigated and from that a threshold value or criterion in determining whether the connected pixels represent a quarter note or a half note is set.

My next concern is to determine the corresponding tone for a given note. The third row in the figure shows the exact coordinates where maximum correlation in each region exists. This is better than computing for the center of mass of such regions. The basis for such returned coordinates is the computed correlation of the two patterns earlier. Knowing the coordinates/locations of the notes will help us in determining what specific tone is to be played. My concern at this moment is to calculate the vertical distance of a given note in the score sheet with respect to the note having the lowest vertical position. Such note will then serve as our reference point. It happens that all notes having the lowest vertical position in the 4 rows in the sheet have C as their corresponding tone. A histogram of the spacing of the notes with respect to the reference note is investigated and hints us the existence of constant spacing between tones, e.g., C-D-E-F and so on. A spacing of 4, on average is observed.

C = 261.63
D = 293.66
E = 329.63
F = 349.23
G = 392
A = 440

The image processed musical score sheet.

The histogram of the areas of the regions observed. It is safe to choose a threshold value of 10 in this case.

Having the necessary tools on how to determine the correct note and tone, we are now ready to let Scilab sing the song for us. Below is the final result of everything.

I am giving myself a grade of 10 for this since I was able to accomplish the main goal of this activity which is to make Scilab sing the song written on a musical score sheet.

Binary Operations

2011-08-04T20:29:00.002+08:00

In this activity, we will apply everything we have learned so far in image processing. We are asked to estimate the area of the cells found in the figure shown below. Note that the image contains cells whose area distribution is uniform. We then hope to proceed to the case wherein some cells differ in sizes with respect to the typical cell size in the figure. The activity is synonymous to what is being done in cancer cell detection. Knowing the size of a typical cell, one can easily detect any possible cancer cell in a given image using image processing techniques.

A digitized image of punched paper using a flatbed scanner.

The first thing that needs to be done is to separate the region of interest ROI from the background. A threshold value can be obtained from the PDF of the grayscale image. The PDF of the given image is shown below. The peaks observed in the figure below represents the background and ROI. A good threshold value lies somewhere between these two peaks and is a minimum, i.e., 0.825.

Histogram of the given image. The threshold value used is 0.825.

Cutting the image to produce 12 sub-images and converting these to binary give the following results. Notice that the resulting binary images are a little bit noisy as a result of the threshold process.

The result of converting the gray scale image to binary using a threshold value of 0.825.

To remove such noise, some morphological operations are needed to be done on these binary images. Previously we have investigated how dilation and erosion work (see previous post). We have the option to utilize some available morphological operations such as opening and closing. These operators are basically a combination of dilation and erosion.

The opening of A by B is obtained by first eroding A by B and then dilating the resulting image by B.

The closing of A by B is obtained by first dilating A by B and then eroding the resulting image by B.

After comparing the results using the possible combinations of the opening and closing operators, I resorted back to using the opening operator only. The results it produce appear to be the best among the many options I have. And the result of using such operator is shown below. Also, note that the structuring element B used here is a circle and its origin is set to be at the origin.

The result using the opening operator.

Using bwlabel function in Scilab to label all connected parts (or blobs) in the sub-images and compute for the area of each blob, the histogram of the area computed is shown below.

The histogram of the area computed using bwlabel. The peak is at around 531.25. And the range of values clustering around this peak is from 387.5 to 571.875.

To find the best estimate for the cell size of the given image, we find the mean of the area values that fall on the specified range, as well as its standard deviation using the mean and stdev functions in Scilab, respectively. The best estimate for the cell size is 490.71 ± 46.97. The actual area of the cell is found by manually measuring its pixel dimensions in Paint. The value obtained is 530.93. This differs from the mean obtained but still falls on the range found.

Now we move on to a more complicated task: that is to identify cancer cells in a given image containing normal cells and those that are cancerous, i.e., cancer cells. The following figure below shows the image used in this part of the activity.

A digital image containing punched papers of two different sizes which correspond to normal cells and cancer cells.

We are to apply binary operations on the image in order to identify the 5 cancer cells that are present. The first thing that needs to be done is to view the histogram of the grayscale version of the image. This histogram is shown below. Notice that there exist two regions where the pixel values tend to cluster. Between these lie the threshold value that we want and it is obviously at around 0.843.

Histogram of the image containing cancer cells. A good threshold value is somewhere around 0.843 as seen on the image.

We are now ready to convert the image to binary using the threshold value that we found. The result is shown below.

The result upon converting the image to binary.

The binary image shown above is too noisy to process. We then apply morphological operations to it in order to remove the unnecessary pixels in the image as a result of the threshold process. The same technique utilized earlier on sub-images containing punched papers of equal sizes will be used here, i.e., the opening operation. The result of this operation on our binary image is shown below.

The result after using the opening operator.

These time the noise found on the previous image is no longer present. In order to tag a cancer cell found on the image, we need to have an idea of its size. This can be done by looking on the histogram of area of the cells in the figure presented previously. This histogram is displayed below. Notice that there exist three regions on the histogram wherein cell sizes tend to cluster. We may assume that these correspond to the normal cells, cancer cells and finally overlapping normal cells. Previously we have determined where the normal cells in the histogram lie and it is found to be at around 500. It is okay to assume that the range of areas of the cancer cells in the image lies somewhere between 870 to 1020. The remaining region belongs to those normal cells that are overlapping. Averaging the areas that lie between the predicted range and computing for the standard deviation give us a best estimate of 940.83333 ± 58.020111 pixels for the cell size of a cancer cell.

Histogram for areas of the cells contained in the processed image.

For tagging the cancer cells, we are to evaluate each labeled part in the image whether its area is within the range of values for area of a cancel cell or not. If it falls in this range, there is a high chance that the said labeled part is a cancer cell. The result is shown below.

The final result! The 5 cancer cells are tagged but also included are some areas corresponding to overlapping cells in the image given.

Our result suggests that the process done is not 100% accurate, i.e., only the cancer cells are tagged. There are a lot of factors why such result is obtained. The regions of overlapping normal cells may have an area which falls in the given range. If the image used contains cells that are completely isolated from each other, this screening process will no longer be a problem at all. But remember, in real images of cells, for example, it is very rare to see cells that are completely isolated from each other. The solution is, therefore, to improve the implementation process, i.e., from obtaining a good threshold value to the morphological operations to be used and predicting the region in the area histogram corresponding to cancer cells. There are a lot of improvements to be done. But for a beginner like me, these results are already satisfying. And for that I am giving myself a grade of 10.

Morphological Operations

2011-07-27T00:29:00.000+08:00

Morphology in image processing refers to shape or structure. Morphological operations are done on binary images in image processing in order to extract more information or to enhance the quality of the image. These said operations affect the shape of the image depending on the output of the a given operation. Using these operations, disconnected blobs, for example, can be joined, shapes found on the image may either be expanded or thinned, and loops can also be closed. These are just some of the many things that can be done using morphological operations.

In a binary image, all pixels having a value of 0 are considered to be part of the background. On the other hand, all pixels having a value of 1 are part of the foreground. It is important to recall some basics in Set Theory since morphological operations have its backbone on this.

Let A be a set in 2-D integer space Z².The elements of Z² are the (x, y) locations of pixels in the image.

If a = (x, y) is an element of A and b is not an element of A, then

If A is a subset of another set B, then

The union of two sets A and B is the set of all elements that belong to either A or B.

The intersection of two sets A and B is the set of all elements that are both in A and B only.

Two sets are said to be mutually exclusive if they have no common elements, i.e.,

The complement of a set A is the set of elements not contained in A.

The difference of two sets A and B denoted by A-B is defined as all elements of A excluding all elements of B.

The reflection of set A is defined and shown below.

The translation of set A by point z = (z₁, z₂) denoted by (A)_z is defined and illustrated below.

Note that the NOT operator is similar to complement and AND is equal to intersection. Moreover, the OR operator is the same as union and XOR can be written as a combination of NOT, AND and OR which is expressed as

We now proceed to some morphological operations such as dilation and erosion.

Dilation
The dilation of A by B is defined as

and thus involves all z's which are translations of the reflection of B that when intersected with A the result is not an empty set. Note that the translation of B is relative to its origin. In an image, the origin of a set is located at its left topmost corner. Moreover, the origin can also be the geometric center of the set. Here, B is the structuring element. Therefore, the effect of dilation is to expand or elongate A in the shape of B. The figure shown below illustrates how dilation works.

The dilation of A by B results to an expanded version of A.

Erosion
The erosion of A by B is defined as

and it is the set of all points z such that B translated by z is contained in A. Unlike dilation, the effect of erosion is to reduce the image such that its shape is similar to B. The figure below illustrates how erosion works.

The erosion of A by B leads to a reduced version of A that inherits the shape of B.

We now proceed on the application of these operations on different image shapes. The structuring elements shown below will be used for both dilation and erosion of different image shapes. Note that the origin is set to be at the top leftmost corner of the structuring element.

The structuring elements used in the dilation and erosion of different image shapes: 1) 2x2 ones; 2) 2x1 ones; 3) 1x2 ones; 4) cross, 3 pixels long and 1 pixel thick; 5) a diagonal line, 2 boxes long.

We first find the dilation and erosion of a 5x5 square by the different structuring elements shown above.

5x5 square.

My predictions for the dilation of the 5x5 square shown on the left by the different structuring elements.

The shaded part (both green and gray) on the images above corresponds to the white part of the images below.

The result using the function dilate() in Scilab.

For the erosion of the given shape, my predictions are shown below.

My predictions for the erosion of the 5x5 square shown on the left by the different structuring elements.

And the results using Scilab are presented below.

The output using the erode() function in Scilab.

Time to move on now to the dilation and erosion of a triangle having a base of length 4 pixels and height of length 3 pixels.


A triangle having a height of 3 pixels and base of 4 pixels.

My predictions for the dilation of the 3x4 triangle shown on the left by the different structuring elements.

The result using dilate() in Scilab.

My predictions for the erosion of the 3x4 triangle shown on the left by the different structuring elements.

The result using erode() in Scilab.

For a hollow 10x10 square with thickness of 2 pixels on all sides, the predicted dilation and erosion of the said image are shown below together with the results using Scilab.

A 10x10 hollow square with a thickness of 2 pixels on all sides.

My predictions for the dilation of the given image by the different structuring elements.

The results using the dilate() function in Scilab.

My predictions for the erosion of the given image by the different structuring elements.

The output using the erode() function in Scilab.

Finally, we investigate the dilation and erosion of a plus sign of length 5 pixels along each line.

A plus sign having a length of 5 pixels along each line.

Its dilation and erosion are tabulated below together with some predictions I made for both operations.

The predictions for the dilation of the given object by the different structuring elements.

The result using the dilate() function in Scilab.

My predictions for the erosion of the said image.

The output using the erode() in Scilab.

There are a lot of morphological operations available in Scilab like thin() and skel(). The function skel(), for example, performs skeletonization of a binary object. This is similar to the function thin() which performs thinning of binary objects. The figure below is created in Paint will serve as a test object for these operations. The results are also shown below.

Sample image.

Using thin(), the following result is obtained. Notice that it returns the similar figure but this time it has thinner boundary lines.

The output of thin().

Moreover, using skel(), the results are shown tabulated below for different key words in the argument of skel().

The "skeleton" of the sample image (left: exterior skeleton, middle: interior skeleton, right: both exterior and interior)

The function skel() is capable of computing both the background and foreground skeleton. The default that it outputs is the interior skeleton which is almost similar to the result using thin().

Understanding how both operations work is a must-have when predicting their corresponding output for each case. And believe it or not, this had consumed much of my time. Nonetheless, I enjoyed doing this activity despite working 12 hours straight in front of my computer. I will give myself a grade of 10 for this since my predictions match with the results from Scilab.

Enhancement in the Frequency Domain

2011-07-21T23:47:00.001+08:00

Any unnecessary pattern, noise for example, in an image can be eliminated by filtering its corresponding frequency in the Fourier domain. If this can be done, enhancement of the desired frequencies corresponding to some patterns in the image is also possible. These techniques sound so simple but some basics are to be mastered. And here are some of the things that must be considered and remembered when utilizing such methods.

The Fourier transform of a convolution of two functions f and g in space is the product of the Fourier transform of the functions.
The convolution of a dirac delta and a function f(t) outputs a replication of f(t) in the location of the dirac delta.

Convolution of a dirac delta and and function f(t).

We now investigate these properties and see how they are applied in image enhancement in the frequency (or Fourier) domain. First, we will show how the convolution theorem applies. Each image shown below contains two identical objects, i.e., circles. The size of the circles is varied increasingly.

The size of the circles are increased with the two dots as the reference.

We now take the Fourier transform of the images shown above. And the results are shown below. Notice that as the size of the circles increases, the size of the fringes in the Fourier plane decreases.


The Fourier transforms of the images shown above.

We have seen the effect of the increase in the size of circular apertures in the Fourier transform. Let us now see the effect on images containing squares.

Increasing aperture size.

The resulting Fourier transform for each aperture size is shown below. Similar trend is observed here, i.e., as the size increases, the size of the interference patterns (fringes) decreases.

The resulting FT's for different aperture sizes.

Now, we try this on images containing a pair of Gaussian functions.

Images containing Gaussian function for different variance values. The variance is varied increasingly from left to right.

The Fourier transform of the above images are shown below. Notice again that the same trend that we observe earlier also applies in here. The size of the spot size that we see in the Fourier plane decreases as the variance is increased.

The Fourier transform of the images containing Gaussian function for different variance values.

Now, we try to convolve a pattern with an image having 10 dirac delta functions at random locations.

The pattern chosen: a cross.

The result of the convolution is shown on the right side of the figure below. Notice that the pattern appears at the 10 locations of the dirac delta functions.

Left: 10 dirac deltas are distributed at different locations randomly. Right: The result of the convolution of the image on the left and the cross pattern.

This shows that indeed when a dirac delta is convolved with a function f, i.e., a given pattern, the result will be an image having the pattern at the location of the dirac delta.

We now see how the Fourier transform of a lattice looks like.

Left: spacing = 50. Right: spacing = 25.

Their Fourier transforms are shown below. Notice that the spacing in the spatial domain is inversely proportional to the spacing in the frequency domain or Fourier space.

The corresponding Fourier transform of the images shown above.

Having some idea of how things work, we now proceed to some applications of image enhancement in frequency domain. Below is a fingerprint image and we want to enhance the ridges that appear on the given image.


Source: http://www.topnews.in/fingerprint-technology-beats-worlds-toughest-tests-2229048

Shown below are the images for the mask used, the Fourier transform of the image and the resulting image when the mask is applied on the Fourier transform.

Left: the mask used in eliminating unwanted frequency values. Middle: the Fourier transform of the fingerprint image. Right: the product of the mask and the FT.

Basically, we create a mask that will cover the unnecessary frequencies on the Fourier domain. The resulting image is then shown below.

The resulting fingerprint image.

Now we move to line removal. The image shown below is taken from http://www.lpi.usra.edu/lunar/missions/apollo/apollo_11/images. Our goal this time is to remove the vertical lines that appear on the given image.

The image we want to enhance, i.e., removing the vertical lines.

The Fourier transform of the above image is shown below. The other image is the result when a mask is applied on the Fourier transform. Since we want to remove the vertical lines, we need to mask the frequencies along the horizontal line, leaving the frequencies at the center and along the vertical line. Recall that the Fourier transform of a grating along the x-axis produced a series of frequency values arranged in the horizontal manner. This justifies why only the horizontal component is masked.

Left: the Fourier transform of the given image. Right: the result when a mask is applied on the image on the left.

The result of the masking is shown below. The vertical lines are no longer prominent in the image shown below.

The vertical lines are successfully eliminated!

We move on to the next part of this activity, i.e., canvas weave modeling and removal. The image that will be used is shown below.

This is a part of an oil painting from the UP Vargas Museum Collection.

Shown below are images for the Fourier transform and the result when a mask is applied to filter the frequencies corresponding to the canvas weave. The mask used here is manually generated unlike the previous ones used in this activity. The choice of what to filter relies on the physical appearance of the canvas weave which is periodic and arranged in a way that it has some kind of symmetry. After taking the Fourier transform, I noticed that six dots forming a hexagon when connected appear on the resulting image. I am having a hint that perhaps these correspond to the canvas weave since these dots also display some kind of symmetry similar to that of the canvas weave.

Left: Fourier transform of the image shown above. Right: the result when a mask is applied on the image shown on the left.

The resulting image is shown below. Observe that the canvas weave is now removed.

Finally! The output after filtering the frequencies corresponding to the canvas weave.

What if we want to recover the canvas weave from the given image? This can be done by inverting the filter used previously so it will appear like the one shown below on the left. The figure on the right is the recovered canvas weave.

The canvas weave on the right is not noticeable.

This activity is perhaps the most tedious and difficult one I have encountered. First, it requires me to remember and apply the things I learned from the previous activities, especially on the part wherein I need to create masks or filters for my images. I think much of my time was wasted on the fingerprint image enhancement. I am not really satisfied with my results for that part. I don't know whether the image I got from the internet is too perfect already to be enhanced. I just realized that image processing is a bit subjective. And sometimes, it is very hard to determine whether the result I am obtaining is already good enough. Overall, I learned a lot from this activity! I am giving myself a grade of 9 for this.

Properties of the 2D Fourier Transform

2011-07-14T14:32:00.000+08:00

The spatial frequency distribution of a signal can be found by finding its Fourier Transform FT. The Fourier theorem states that any given signal or image is just a superposition of sinusoids having different frequency values. In this activity, we investigate how a rotation in the orientation of a sinusoid gives rise to a unique and different FT that is also rotated.

Some of the sample images investigated are shown below.

First row: square, square annulus. Second row: circle, annulus. Third: double-slit along x-axis, two dots along x-axis.

The Fourier transform of the above images are then obtained and presented below in an order similar to the one above.

The Fourier transform of the sample images.

Observe the difference in the interference pattern produced by the sample images.

Having an idea of the FT of different image patterns, we now move to the investigation of the anamorphic property of Fourier transform. The image below shows a sinuoid captured and generated at different frequency values, i.e., f = 2, 4, 16, 32.

First row: f = 2, f = 4. Second row: f = 16, f = 32.

Their corresponding Fourier transforms are shown below. Notice that as f increases, the two dominant frequencies seen move away from each other. The output generated is consistent with our definition of what Fourier transform is. For a function having a single-valued frequency f, we expect that a single peak will appear in its Fourier transform. But since the function investigated is odd, another peak will be observed which corresponds to the negative frequency f.

The Fourier transforms for different f values.

But since in digital images there are no negative values, we eliminate these by adding a bias to the sinusoid. The peak that appears in the origin is the result of the addition of bias since the frequency f of any constant is equal to zero. For the addition of a non-constant bias, i.e., a sinusoid with frequency f = 1, to the original sinusoid having a frequency f = 4, the Fourier transform is observed to have two additional symmetric peaks. These correspond to both positive and negative frequencies f of the non-constant bias.

Top: Sinusoid with f = 4 having a constant bias = 10. Bottom: Sinusoid with f = 4 with a non-constant bias in the form a sinusoid having a frequency f = 1.

We try to rotate the sinusoid and observe what happens to its Fourier transform.

The sinusoid rotated at different angles theta. First row: theta = -10, 10. Second row: theta = -45, 45.

The Fourier transforms of the above images are shown below. Notice the change in the orientation.

The FT of the sinusoid rotated at different angles.

We generate a sinusoid that is both dependent on x and y. We expect that in its Fourier transform, four peaks will be observed, i.e., 2 peaks for x and another 2 peaks for y (both positive and negative frequency values).

Observe the four peaks in the right image.

We can also try different combinations of sinusoids. Below are images generated by combining the product of sinusoids along x and y with sinusoid/s rotated at a given angle. Also, shown are their corresponding Fourier transforms.

First row: Product of two sinuoids along x and y having a frequency f =4 added to a sinusoid inclined at an angle -10. Second row: Product of two sinuoids along x and y having a frequency f =4 added to sinusoids inclined at an angles -10 and 45.

The resulting FT's are predictable since the signal being investigated is made up of components that we have just tackled. The output is just then the superposition of all the FTs of the components.

This activity is easy but tedious. But I still enjoyed doing this. It is amazing to see patterns that are familiar in the field of optics. Doing this activity made me further understand how Fourier transform works. And for that, I am giving myself a grade of 10.

Fourier Transform Model of Image Formation

2011-07-12T22:23:00.000+08:00

Lens as Fourier Transformer

Fourier Transform (FT), as its name implies, recasts information from one domain to another. In particular, it decomposes a signal, which is in time domain, to its constituent frequencies (frequency domain). Here, we see the transformation from the time domain to the frequency domain. FT is well utilized in the field of signal processing. If one's goal is to eliminate the noise in a given signal, working on the frequency domain is advantageous compared to the case in time domain. Looking on the FT of the said signal gives one the hint of what frequency value gives rise to the noise observed in the signal. Thus, eliminating the noise is easier with the help of FT. To find the Fourier transform of one signal, the equation shown below is used.

Fourier transform equation

An implementation of this, that is fast and efficient, had been created by Cooley and Tukey. And it is called the Fast Fourier Transform algorithm (FFT). This is so widely used today that software like Scilab uses FFT in the implementation of Fourier transform. In Scilab, it is available as fft() for 1D and fft2() for 2D. The list below enumerates the properties of FFT for both 1D and 2D.

FFT2 is most efficient if the size of the image is a power of 2 (e.g. 64x64, 128x128, 256x256, etc.)
The output of FFT2 has quadrants along the diagonal interchanged.

An example of what happens when FFT2 is applied to a given input. To return the quadrants to their original position, fftshift() is utilized.

Take note that FFT outputs a complex number array. The abs() function can be used to view the intensity of the output.
The inverse FFT of an image is the same as the forward FFT with the image inverted.

To see how these things work, we apply the said functions on images. The first is a 128x128 image containing a white circle in black background.

Test subject no. 1

Its Fourier transform is shown below. Notice the corners of the square.

The result when fft2() is applied on test object no. 1.

Remember that when using FFT, the quadrants in the image are interchanged. Thus, in order to see the actual result, fftshift() is used.

The result when fftshift() is used.

Performing another fft2() on the test object results to the image shown below.

The result of implementing fft2() twice on the image.

It appears that after two consecutive fft2() implementation on the test object the result is similar to the original image, i.e., the test object. But we do not know whether the result is just the original image itself. If ever this is not the case, the symmetry of the circle hinders us to any possible observation on the output image. Thus, it is perfect now to have a second test object that has no perfect symmetry.

Test object no. 2

Applying fft2() on the above test object, the output is shown below.

Always remember how fft2() works. Notice the corners once again.

In order to see the real output, i.e., similar quadrant locations, fftshift() is used.

The output image is better to analyze when fftshift() is utilized.

Again, another fft2() is applied on the test object resulting to the image shown below.

The result is similar to the test object but this time inverted.

Remember that the inverse Fourier transform of an image is just the Fourier transform of the image inverted. Thus, if we want to obtain the original image, i.e., inverse FT of the FT of the image, the FT of the image must first be inverted before another FT is applied. This explains why the recovered image is inverted.

Convolution
The convolution between two 2-D functions f and g is given by

Convolution integral

or is simply denoted as

h = f *g

Take note that it is a linear operation. This implies that if both f and g are recast by linear transformation (e.g. Laplace, Fourier transform), they obey the convolution theorem given by

H = FG

where H, F and G are the transforms of h, f and g, respectively. Convolution is basically the smearing of one function against another function. This gives us an idea that the resulting output h appears to be almost a combination of f and g, i.e., it looks like both f and g.

Convolution has its application in simulating the linear regime of instruments or detector devices used in imaging. The functions to be convolved, f and g, can be the object and the transfer function of the system being investigated, respectively. Thus, when convolution is done, the result represents the output of the system. It is expected that the output is not perfectly similar to the original image.

Now, we try to simulate how imaging happens using two images: one representing the lens, the other which is the test object.

The test object.

The size of the aperture of the lens is varied and the resulting image is observed.

The lens used in varying aperture sizes.

When the test object is convolved with the lens having different aperture sizes, the results are as follows.

The resulting images are shown for different aperture sizes. First row images correspond to the first two aperture sizes and the second row is for the last two aperture sizes presented previously.

As the size of the aperture of the lens increases, the quality of the image produced also increases. The quality of the reconstruction is dependent on the amount of light rays reflected from the object entering the lens. For the reconstruction to be of good quality, the size of the aperture must be wide enough to gather all the rays that reflected from the object. But in reality, this is not the case because most lens used in imaging has finite aperture size, i.e., the quality of the image produced is also limited.

Correlation
The correlation between two functions f and g, which are both 2-dimensional functions, is given by

Correlation integral

or is written simply as

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

Observe the similarity between the correlation integral and convolution integral. The difference is on the signs inside the function g. We expect that there also exists a correlation theorem similar to that of the convolution theorem. This then holds for linear transforms of both f and g and is given by

P = F*G

where P, F and G are the Fourier transforms of p, f and g, respectively. The asterisk above F indicates complex conjugation.The correlation p measures the degree of similarity between two functions f and g. High correlation value at a certain point signifies that the two functions are most likely identical at that point. This has applications on template matching or pattern recognition.Something nice happens when the function g is even. Both integrals for correlation and convolution are the same.

To see how correlation works in template matching, we correlate the image shown below with another image containing only the letter A having the same font style and size with that of the other image.

Shown is the test object which we want to correlate with an image containing the letter A. Note that the letter A has the same font style and size similar to that of the A's found on the test object.

Our goal is to locate all the A's in the test object using a template containing the letter A. We utilize the correlation theorem to do this. The result is shown below.

The result is pretty amazing. Take a look on the locations of the A's in the previous image and locate it on the image. What do you see?

The locations of A have the highest correlation values as marked by the observed twinkles on the sites. This suggests that indeed correlation can be applied on template matching such as this one.

Edge detection can also be done using the convolution integral. We want to convolve the VIP image shown above with different patterns. Some of the patterns used are shown below.

The patterns used. Take note that the sum of all the elements of the matrix must be zero.

The results are shown below for the convolution of the VIP image with the different patterns.

First row: horizontal, verticals. Second row: spot, diagonal patterns.

Since convolution took place, we expect that the resulting image appears like the convolved images, i.e., the VIP image and the pattern. We have seen that the above results contain the word VIP or parts of it. It is observed that edge detection is best when the image is convolved with a spot pattern. Convolution with a horizontal patterns outputs most edges that are horizontal. On the other hand, a vertical pattern is good in displaying the vertical edges. The diagonal pattern has the least accuracy in terms of edge detection.

For this activity, I am giving myself a grade of 10. This activity is a good starting point for beginners in FT model of image formation. I have encountered the terms convolution and correlation, as well as Fourier transform, so many times that I have memorized them. It feels good to know how these things work, i.e., their applications, the things that can be done using these concepts. There was this one activity in 185 wherein I encountered the Scilab function fftshift(). It is quite embarrassing to admit right now that through this activity I am able to discover what this function really does. When trying to attack something that appears to be complex, starting on the basics is really and should be the first step to be done.

Enhancement by Histogram Manipulation

2011-07-06T02:02:00.002+08:00

The graylevel probability distribution function (PDF) of a grayscale image is equal to its histogram normalized to the total number of pixels found in the image. The quality of an image can be improved by manipulating its so-called histogram. By doing so, certain features of an image can be improved. Also, the response of different imaging systems can also be mimicked by histogram manipulation. The grayscale PDF of an image can then be changed by backprojecting the grayscale values using the cumulative distribution function (CDF) of a desired PDF.

The following statements further explain how the technique works.

Given the probability distribution function p1(r) of the graylevels r of an image and the corresponding cumulative distribution function T(r)

Equation 1

where g is a dummy variable here, we plan to map all the graylevels r to a different set of graylevels z. This then gives us a new cumulative distribution function G(z) given by

Equation 2

wherein we have a new PDF p2(z) corresponding to the transformed image. Take note also that t is a dummy variable. To summarize all, the equation below shows how the graylevels r are map to a different set of graylevels z.

Equation 3

1-2. For a grayscale level r, find its corresponding CDF value T(r).
3. Find the value of T(r) in the y-axis of the CDF G(z).
4. Finally, find the corresponding z value of the obtained G(z).

To apply all the concepts that we just learned, we are to improve the quality of the image shown below, i.e., through histogram manipulation.

A photo of the sunset at Roxas boulevard courtesy of my friend Desiree Badua.

The image is then converted to grayscale and the result is as follows. Also, the histogram of the grayscale image is shown below. The histogram shows that most of the graylevel values are concentrated near the area corresponding to black. This is evident in the image shown below. In order to make the image appear brighter, we want a histogram wherein the graylevel values are well spread from 0 to 255.

The grayscale image.

The histogram of the grayscale image.

In order to find for the probability distribution function of the grayscale image, we need to normalized the histogram shown above to the total number of pixels in the image, i.e., 1024x768 pixels.

The PDF or the normalized histogram of the grayscale image. Take note that the range of x is from 0 to 255.

The CDF of the above PDF can be found by using equation 1. When implemented using Scilab, the CDF can be found by performing a cumulative sum over all values of the PDF. The result is shown below.

The values of x range from 0 to 255.

We then create our desired CDF to be able to start the image enhancement. Based on this CDF, the image will now be transformed. We first try a linear CDF shown below.

A desired linear CDF given by the equation G(z) = 1/255 z.

The main idea is to map the original set of graylevels r (0-255) to a different set of graylevels z. This is what the following code snippet performs. Having obtained the new set of graylevels z, we are to set the graylevel value ri of each pixel in the image to its corresponding graylevel zi.

Remember that each element of vector z corresponds to each element of r. Also, the key here is equation 3.

Instead of getting analytically the inverse of the CDF G(z), I tried using the function [yp] = interp1(x, y, xp) in Scilab 5.3.2. Apparently, I am given a set vectors for x, y and yp, i.e., z, G(z) and T(r), respectively. When using this function, xp and yp should not be used interchangeably. Unfortunately, we can not utilize this one. This made me proceed to my last resort: finding the inverse of the desired CDF G(z).

The modified image using a linear CDF.

Also, its corresponding PDF and CDF are shown below. It is clear that this time the graylevels are at least well distributed from 0 (black) to 255(white).

The PDF of the enhanced image.

The corresponding CDF of the enhanced image.

Let us then try using a nonlinear CDF and see the effect on the image. A parabolic CDF given by the equation of the form G(z) = Az² is shown below.

A nonlinear CDF G(z) = 3z²/195075.

The code used to generate the modified image for this CDF is similar with that of the linear CDF.

The code for image enhancement using the nonlinear CDF.

The resulting image is shown below.

Notice that this appears to be whiter in color compared to the previously enhanced image using a linear CDF.

For the above modified image, its PDF and CDF are shown below.

The PDF of the enhanced image using a nonlinear CDF.

It appears that the recovered CDF for this enhanced image is no longer smooth unlike the original CDF that we wanted.

Histogram manipulation can be found in some image processing software (e.g. GIMP). Shown below are snapshots showing the GIMP environment and how histogram manipulation is easily done.

Notice that the desired CDF in this case is not actually a function since there is no one-to-one correspondence that exists between the x and y-axes.

In this case, I tried to mimick a step function for the desired CDF.

What I have learned in this activity is that through histogram manipulation, information stored in an image that appears to be invisible to the eye can be recovered or seen. We are lucky that most of the software used in image processing today has functions that deal with histogram manipulation. Doing this activity made me appreciate a lot of things. Slowly, those stereotypes I had for image processing are vanishing. I did enjoy doing this activity. And for that, I am giving myself a grade of 10.

References:

Enhancement by Histogram Manipulation 2010, Soriano M., Applied Physics 186.

Area Estimation for Images with Defined Edges

2011-06-29T01:04:00.008+08:00

Area measurement in images play a significant role in research. As much as possible, people nowadays want everything to be automated. With this, estimating an area based on images is possible. There are a lot of things that one can do if equipped with the necessary skills in area estimation using images. In research, for example, this technique is used to automate the detection of cancerous cells. Also, given an image, the area of a forest can be obtained using the said technique (remote sensing).

There are two ways in area estimation. The first technique uses Green's Theorem. The theorem relates the line integral around a closed curve C, for example, to the double integral over the plane region D bounded by C, and vice versa.

The closed curve C, in the counterclockwise direction, that bounds the region D.

Since the theorem relates a closed integral to a double integral, the equation below summarizes it with the need to consider two continuous functions, F1 and F2, which have continuous partial derivatives in the region containing D.

Equation 1

Note that if the contour C is in the counterclockwise direction, the region of interest D is located to the left of the boundary. C is then a collection of points (x, y) that make up the contour or the boundary of region D. The succeeding lines demonstrate how the area of region D can be obtained.
If we let F1 = 0, F2 = x and F1 = -y, F2 = 0, equation 1, for each case, becomes

Equation 2

Equation 3

Both equations 2 and 3 represent the area of region D. Averaging the two expressions for area, we have

Equation 4

The discrete form of equation 4 is given by

Equation 5

wherein N is the total number of pixels in the boundary or contour C. If the contour is made up of N pixels, i.e., the total number of pixel coordinates (x, y) is equal to N, then there exists a problem when we try to access the pixel coordinates (xN+1, yN+1). This is solved by letting (x1, y1) = (xN+1, yN+1). This is done because we want to close our loop.

The second technique in area estimation is through morphological operations. This considers connected pixels as blob. It estimates the area in the image by counting the pixels contained inside the blob.

To test the accuracy (validity?) of equation 5, I generated some figures whose area can be easily computed, i.e., regular geometric shapes.

In this activity, I used the 4.1.2 version of Scilab together with the SIP toolbox since SIVP does not contain the follow command. Basically, the follow command is a contour follower. It returns the coordinates of the points that make up the contour or boundary of the region of interest. Knowing the coordinates simplifies the job. The only thing to be done then is to insert it into equation 5.

The area of any regular geometric shape can be obtained using the following snippets.

The code used to calculate the area of any regular geometric shape. It is recommended that the image to be used is saved as a Bitmap image.

And here are computed values for the geometric shapes I have chosen, together with the corresponding actual area measurement and the percentage error for each case.

The percentage error values suggest that exploiting equation 5 in the calculation of the area of given regular shape is quite reasonable up to some degree. The higher the resolution of the image, the greater the possibility that the actual and calculated area coincide.

The method I used in measuring the actual area of the given shape is to measure its dimensions in pixels. This method applies whenever the shape to be measured is regular, i.e., the known formulas for area calculation can be applied directly.

A problem arises when the area to be measured has an irregular shape. This is the case in the real world. It is then important to extend our knowledge of area measurement to objects that are irregular in shape. As an example, I want to measure the area bounded by Roces St., Roxas and Osmeña Avenue.

A satellite image obtained using Google Maps showing a portion of land located in the University of the Philippines Diliman.

In order to begin with the area measurement, the following steps are to be done. I opted for the use of Paint in tracing the boundary of the parcel of land shown. The result is shown below.


The boundary of the parcel of land now traced using Paint.

Since Green's theorem can only be implemented when the region of interest is white in color and the background is black, I planned to recover the black outlines I drew using the technique I learned from Activity 3.

The recovered black outlines from the map shown previously. This appears to be smaller compared to the images shown above but they all have the same area. This is because of the white spaces surrounding the black outline. Trust me. I have measured and compared the dimensions of the three images.

The code I used in converting the image to binary, together with the histogram plot (which is the basis for the threshold value used), is shown below.

The number of bins, which is 1000 in this case, may vary. Also, we want the output image to have the same size as the input image. I used imshow (im_bw) instead of the regular imshow (im_bw, []) that rescales the image.

Notice the very sharp peak at around 0 which corresponds to the color black. We then have an idea that if our only interest is the black outline that we saw earlier, the safe threshold value to be used is around 0, which in my case 0.001.

The resulting image is then processed using Paint and Scilab.

And this image is ready to undergo area measurement!

The same code used in calculating the area of regular objects is used. And it gives as the following result.

The calculated area is in pixels. We need to devise a way on how to convert it to a real quantity.

It is then important to note the scale used in the map taken using Google Maps. I plan to count the number of pixels in between tick marks, i.e., the number of pixels that make up a 50m.

The scale corresponding to the image I obtained using Google Maps.

Based on pixel measurement of the above image, there are 87 pixels in between the two tick marks in meters. This means that 50m = 87pixels. Then using dimensional analysis,

The calculated area using Green's Theorem.

Okay. I need then to compare the value I obtained with the actual area measurement of the parcel of land. Hmmm. I tried to Google for the area of the lot I selected. Unfortunately, I found no results. What to do now? I am lucky that somehow the lot I selected has a shape that is almost a combination of regular geometric shapes, i.e., a triangle and a rectangle. My plan now is to compare my measurement with the area of the triangle-plus-rectangle combination. I am very aware that this method is not very accurate. Thus, I will just be comparing in terms of the order of magnitude between the two measurements.

Notice the red triangle and the green rectangle.

After the tedious measurement in pixel, the following table tabulates the area I obtained using my clever method.

Hmmm. It seems that the actual measurement is almost equal to the calculated area using Green's theorem with a percentage error of 0.51%. Not bad at all.

Having presented everything, it is time for me to rate myself for this activity. I am giving myself a grade of 10 for this. I am learning a lot from this course. And the good thing is I am able to apply what I learned from the previous activities to the current one. I am not quite satisfied with how I obtained the actual measurement of the lot. I believe that there exists a better way of figuring it out. There is, I believe, some data available concerning parcel of lots located in UP Diliman. Doing this activity reminds me of my father's work. He is a geodetic engineer. And at home, I am quite exposed to what he is doing, i.e., land titling, surveying. This also reminds me of the AutoCAD that he is using. What I remember is that you input coordinates, and after some clicking, it returns the area of a selected parcel of land in the plot.

And I hope to learn more in the next activities.

Image Types and Formats

2011-06-25T23:51:00.008+08:00

When saving an image, I usually use the JPEG/JPG format. It was only last year that I came to ponder on the existence of other formats, such as .png, .tif, .gif and .bmp. There must be a difference among these file formats. And of course, each must have its odd traits. But if my only aim at a given moment is to save the picture, then saving it using the very well-known .jpeg/.jpg file extension is an accomplishment. Apparently, there exist huge differences among available image formats. And in image processing, the format of the image does matter. I jot down some of the common image formats and the benefits of using such.

JPEG - Its name stands for Joint Photographic Experts Group. This group of people are responsible for the development of .jpeg. It is well known for its lossy compression. When you want your file to be compact, i.e., small file size, this image format is the best to use. As its name suggests, not all of the data found on the original file will be saved if a .jpeg extension is used. One needs to choose small storage size over image quality.
GIF - It stands for Graphics Interchange Format. It was introduced by CompuServe in 1987. It supports animation but has some problem with color limitation. It only supports a maximum of 256 colors. Also, unlike JPEG, it uses lossless data compression technique. Examples (from http://www.anikaos.com/0001-gifs-hamtaro1.html) are shown below. It is best to use GIF when creating logos. It is not advisable for people who are into photography and publishing industry to use this since it only holds up to 256 colors.

PNG - (pronounced Ping) It is short for Portable Network Graphics. It is a better alternative for the GIF format. In terms of image compression, PNG is better than JPEG.
TIFF - The name stands for Tagged Image File Formal. It was developed by Aldus and as of 2009, it is under Adobe Systems. It is popular among people who are into photography, publishing industry and the like. As its name suggests, developers can include tags containing proprietary information inside a TIFF file.
BMP - It is more commonly known as Bitmap Image File or Device Independent Bitmap (DIB) file format. This file format is used in storing bitmap images. This is a popular image format used on Windows and OS/2 operating systems. Data compression is optional in this case.

Now that we have an idea of the available file formats, it is good to have some basic knowledge of the different image types. The basic types are enumerated below, together with some examples.

Truecolor Images - Examples of this image type are those captured by the digital cameras that we have today. The colors that we see are the results of the overlapping of its 3 channels that contain the primary colors red, blue and green.

A photo of a cute baby taken from: http://www.flickr.com/photos/asam/432194779/.

The image shown above has the following properties:

Note that dpi means dots per inch. It expresses the resolution of the image. The higher the resolution of an image, the more you can zoom in on the image.

Binary Images - The first things that come into your mind when you hear the word binary are 1's and 0's. The image is basically made up of 1's and 0's which in color are represented by the colors white and black, respectively.

Photo taken from: http://free-extras.com/images/black_and_white_cats-1541.htm

Here are the image properties of the binary image shown above:

Grayscale images - These are images whose colors ranges between 0 (black) and 1(white). Images of this type are gray in color.

Photo of a baby and a dog taken from: http://www.likecool.com/Cute_Baby--Pic--Gear.html.

The above image has the properties shown below:

Indexed images - These are basically images whose colors are represented by numbers indicating the corresponding index of colors in a color map. Images of this type contains two sets of information: the original image and its color map.

Photo taken from: http://annaleeblysse.blogspot.com/2010/09/lava-fountain.html.

The indexed image shown above has the properties displayed below.

The properties of an image can be found by right clicking the file then choosing Properties and then Details. This applies when the operating system you are using is Windows 7. The disadvantage of using this is that it does not give you the image type. In Scilab 5.3.2, these properties can be seen by simply typing

info = imfinfo(filename)

This will output the following set of properties:

And here are some of the advanced image types available.

High dynamic range images

Photo taken from: http://www.modny73.com/artwork/photography/30-hdr-images-of-famous-cities-of-the-world

3D images

Photo taken from: http://www.fotocommunity.com/pc/pc/display/16115020.

Multi or hyperspectral images

A satellite image taken from http://www.gearthblog.com/blog/archives/2008/10/new_satellite_imagery_galore_coming.html.

Temporal images or videos - Presented below is a sample video I found in my laptop c/o Windows 7.

One can convert an image from one type to another. The image shown below is an example of a truecolor type.

Photo taken from: http://www.desktoprating.com/wallpapers/4-cute-puppies-wallpaper-640x480.htm.

This image can be converted to a binary type using the snippet of code below.

The code I used to convert the image from truecolor to binary.

And the result is shown below. Note that in converting an image to binary, a threshold value which ranges from 0 to 1 needs to be set. In this case I chose a value of 0.5.

The binary image with a threshold value of 0.5.

To have a hint of what is happening when we increase or decrease the threshold value, look at the figure below. It shows how the threshold value plays a significant role in the process of converting the image to binary.

From left to right: threshold value starting from zero with an increment of 0.1.

When using Scilab 5.3.2 together with the SIVP toolbox, the corresponding function for im2gray in Scilab 4.1.2 plus SIP toolbox is rgb2gray. In converting a truecolor image to grayscale, the following set of lines was used.

From truecolor to grayscale.

The resulting image is this.

The cute puppies in grayscale.

The next thing I needed to do is to convert the scanned image I used in Activity 1 to grayscale. The common error encountered when doing this is that the stacksize is smaller than the required value.


When you encounter problem like this, try increasing the stacksize you are using.

To increase the stacksize,

stacksize(x); //where x = [180000,268435454]

Earlier, we converted a truecolor image to binary. Converting an image to binary needs a threshold value. Looking on the histogram of an image gives you clue on what threshold value needs to be used in order to separate a region of interest from the background. This is possible when in the histogram of an image, the pixel intensities are crowded around two values. In this case, the threshold value to be used to separate the region of interest from the background is somewhere between the two peaks. The image shown below is the gray scale version of the scanned image I used in Activity 1.

The scanned image converted to grayscale.

And for the histogram of the above image we have


Notice the two peaks in the histogram shown below. Each peak either represents the background or the region of interest, i.e., the lines. A good threshold value is somewhere between the two peaks (145.71).

Using the value selected as the basis for thresholding, the result is shown below.

The binary image as a result of a threshold value of ~0.57.

The code used to generate the binary image shown above.

From grayscale to binary. Plus displaying the histogram of the gray scale image.

To compare how good the thresholding process, the original scanned image is shown below. Notice that due to the quality of original image, i.e., the scanning process or how it was scanned, there exists on the right side a column that is black. The quality of the thresholding process is therefore affected. Nevertheless, the result is quite okay upon ignoring the limitations due to the quality of the image.

The original image.

I have learned a lot from this activity especially on when to use a specific image format. I was introduced to the vast variety of image formats available, including the image types. I would give my self a grade of 10 for this activity.

Installing SIVP 0.5.2 on Scilab 5.3.2

2011-06-21T10:33:00.000+08:00

I finally was able to install SIVP on Scilab 5.3.2! Apparently the Install file that I read was misleading. Installing SIVP is easy. There are two simple ways to install it. Click here to find out how.

Scilab Basics

2011-06-19T16:46:00.042+08:00

I had a hard time figuring out how to install the SIVP toolbox for Scilab. The version of Scilab that I have is 5.3.0 which for compatibility reason needs an SIVP toolbox, not an SIP.

The error I encountered. Hmm, what to do now?

The common suggestion I got from those I asked was to use Scilab 5.3.2 instead of 5.3.0. I installed the 5.3.2 together with the SIVP toolbox 0.5.2. And guess what the same error appeared. A C compiler is needed in the installation process. But I don't have any idea how to get rid of this error. Instead, I went to the older version of Scilab which is 4.1.2 and installed SIP toolbox 0.4.0 and ImageMagick 6.2.5-5. I again encountered a problem but this time a solution is available.If you are having a difficult time porting SIP toolbox into Scilab, read this link.

I am a newbie in Scilab but had already used Matlab for my research. Scilab, according to many, is almost similar to Matlab. The good thing about Scilab is that it is a free software. Anyone can use it.

Of course, whenever a programming language is new to you, practice is needed. The first thing I generated, given the snippets of code as an example, was a circular aperture.


The circular aperture I generated using the snippets shown below.

The free code found in the activity manual.

And now it is my time to generate my own figures. The easiest one I think to make is the annular. Here, I used the technique of masking. The code is shown below.


The code I used to generate an annular. Here, a mask B was created.

I recycled the code I used in generating a circular aperture and added some few lines since an annular is almost similar to a circular aperture with a smaller circle inside it. I created a matrix B having the same size as A but this time it is composed of 1's. I set to zero all elements of B which were less than 0.3 (r < 0.3). I then add A and B and set it to C.

The annular I generated.

The second one I generated was a square aperture. To generate this is simple. First, you need to know the center of the square, i.e., the matrix. In this case, it is matrix A with all its elements set to zero. The dimension of the square can be any number as long as it is less than that of the matrix.


The code for a square aperture.

And ta-da! Here is the square aperture, with an area of 40x40, that I generated.

The square aperture with an area of 40x40.

The next thing I generated was a sinusoid along the x-direction, a.k.a. corrugated roof. Since we wanted a sinusoid along the x-direction, for sure, the X matrix would be the argument of sine. I had some time figuring out the right argument for sine. For our synthetic image to be a corrugated roof, a spatial frequency factor, together with a pi, must be inserted as an argument of sine. And here it is!

The corrugated roof with a spatial frequency of 10!

And for the code... The factors f and pi in the argument of sine played significant role in making the synthetic image looked like a corrugated roof. A plain A = sin(X) would not output the right figure.


The code for a corrugated roof. Notice the factors f and pi in the argument of sine!

Okay. For the next synthetic image, we have a grating along the x direction. I took advantage of the sine function in this case since I realized that a grating has a repetitive property. Together with the functions round and abs, I was able to generate a grating. Notice that the argument of sine contained a factor of 10. Varying this, i.e., inserting other number, changes the output image. This consumed much of my time figuring out how to generalize it. That means controlling the number of gratings.

This still needs to be generalized.

And the result...


Not bad at all.

And finally, a circular aperture with graded transparency was created. Depending on the value of sigma used, the spread of the Gaussian function varies. It is good to know that the variance measures the spread of the function. Adjusting the value of sigma in the snippet of code below gives you a different output.

Try different values for sigma. And compare the results. Also, you may try adding a factor for the Gaussian function.

And the result is this!

A circular aperture with a Gaussian transparency (sigma = 0.5).

The double-slit experiment or Young's experiment is very known in the field of optics. Since most of the synthetic images I created are being used in optics, I decided to come up with patterns that are at least familiar and useful in optics.

The well-known double- slit!

The code I used to generate such image is shown below. Notice that the width of each slit a and the separation distance d of the slits may vary. In this case, I chose a = 5 and d = 20.

You may try different values for a and d too!

I remember some of the apertures I encountered in Fresnel diffraction and I want to duplicate them. An example is shown below.

A man-with-a-hat aperture.

Half of the bigger circle needs to be cropped. Masking technique was used to remove the lower half of the circle. Also, a smaller circle was created inside the bigger one. See the code below.

The code for the man-with-a-hat aperture.

There are a lot of interesting synthetic images one can create using Scilab or any numerical software. Imagination and algorithm are the key! I enjoyed doing this activity and I would not be hesitant to give my self a grade of 13 for this. I hope to be able to enhance my Scilab skills in the future, i.e., in the succeeding weeks!

Hello AP186!

2011-06-16T12:48:00.000+08:00

Digital scanning is the first activity that we performed in my Applied Physics 186 class. During the time when computers were not yet used in generating figures, people had the last resort (or only option?) to draw their figures and graph manually. Such a tedious job for them! The french curve is one among the many useful tools used in drawing smooth curves. Thanks to everyone who contributed in making things easier for us these days.

graphicsdirect.co.uk

I am not really into editing pictures or anything that has something to do with Paint, Adobe Photoshop and the like. I became a bit scared when I learned that our activity would be about counting pixels and doing ratio and proportion. Contrary to what I expected, everything went well.

The figure that I used in this activity was from Tracy. It came from a 1973 thesis. The first problem I encountered was how to rotate the figure since it was kind of tilted.

The original image I got from Tracy.

It was Rusty's idea to rotate it using Word. For about 15 minutes, I tried figuring out how to do it in Word. And finally after doing multiple right clicks in order to find the proper command, the image was now ready to be digitally scanned!

The rotated image ready to be digitally scanned!

With the image now ready, I removed all the excess areas on the graph leaving everything enclosed by the intersecting x and y axes.

And now let's begin!

I located the origin and recorded its pixel coordinates. I also measured the length of each axis and the corresponding location of the tick marks in the pixel world. The first technique I used was to know how many pixels make up an interval in both axes. Then I selected few points in the graph and obtained their coordinates in pixel. With this, I could go directly to ratio and proportion in finding the exact physical quantities.

The final result! The original image superimposed with the reconstructed graph.

And now for the second technique, the following plots were made.

A plot of the y-axis tick marks versus their corresponding y-coordinate in pixel.

A plot of the x-axis tick marks versus their corresponding x-coordinate in pixel.

This technique consumed much of my time because the image I used in determining the x and y-coordinates in pixel of the tick marks for both axes was the one that was not yet cropped. So I needed to start all over again. I was kind of disappointed to see that the result of this is not the same as that of the first one.

Not a good one. But better than the first one I got using this technique.

It seems that the reconstructed graph is shifted some units away from the original one.

With all the results shown, a grade of 8 will be okay for me. The first technique did work but a lot of approximations were done. The second technique presented here is better than the first but asks you to be accurate and precise in determining the coordinates in pixel. It seems to me now that my measurements were not that accurate and precise.