miguel AP 186 blog

Activity 19 Restoration of blurred image

2009-10-05T06:16:00.003-07:00

Our final activity is again about image restoration. However unlike the previous activity, where we tackled images with additive noise, this time we will be restoring blurred images. Basically a blurred image is due to the effect of a degradation function to a given image. Let us also remember that in our earlier activities we have already encountered blurred images as a result of convolving apertures or filtering in the Fourier space. In this activity we will only be focused on blurred images caused by uniform linear motion with an additive gaussian noise.

Whatever the degradation process maybe, it can always be modeled as a convolution of the image and the degradation function such that in inverse space we obtained a blurred image given by

Equation 1

where H is the transfer function of the blurring process and N is the additive noise. For this activity we will be using the known degradation function of a uniform linear motion blur (equation 2). In equation 2, parameters a and b are the total displacements along x and y directions respectively and T is the total exposure time.

Equation 2

Just like in our previous activity the restoration process involves applying a filter or operation to obtain an estimate of the image from the degraded one. This time we will be using the Weiner (minimum mean square error) filter which incorporates both the degradation function and the statistics of the additive noise. The application of this filter operates in inverse space and is given by

Equation 3

where Sn and Sf are the power spectra of the noise and the ideal (undegraded) image. However, in real life situations it is highly improbable for us to have information on both the additive noise and the ideal (undegraded) image. For such cases equation 3 will be approximated into

Equation 4

where Sn/Sf is replaced by a constant K.

In this activity we use equations 1 and 2 on a gray scale image (figure 1) to produce a blurred image. Then we restore this blurred image using equation 3 or equation 4 with different values of K. We also do all of these for different gaussian noise strengths.

Figure 1

Figure 2. Blurred and restored image for different values of K, T, a, and b. We set the additive noise to be equal to zero

Figures 2 shows our results if we only have an image that is motioned blurred and no additive noise is present. From these images we clearly see that decreasing K improves the quality of our restoration. Furthermore, at the extreme case when K = 0 we arrive at a perfect restoration. These results are actually expected if we analyze equations 1 and 4. Doing simple math, at K = 0, equation 4 reduces to equation 5.

Equation 5

And since in this case we set the additive noise to zero, the last term in equation 5 vanishes and we are left with a perfect restoration.

(click to enlarge)
Figure 3. Blurred and restored image for different values of K, T, a, and b. Noise strength is sigma=0.001

In figure 3 we show the results of restoration when we consider the blurred image with an additive gaussian noise. Again similar to the case when there is no noise having a small value of K gives a very good quality of restoration. However, in this case if we further decrease K we then start to obtain dark noisy images. And, unlike the case of no noise, at K = 0 we obtain a meaningless restored image (almost pure noise).

One more notable result shown here is that even by using the exact power spectra of the noise and the original image, as in equation 3, we still do not obtain a perfect reconstruction. The restored image is indeed a significant improvement compared the blurred image but it is still not perfect.

I give myself a grade of 9 in this activity. The results we have obtained were very convincing but still left some questions unanswered. I would like to thank Master, Jaya, Orly, and Thirdy for helping me realize that I was wrong.

Main reference

A19 – Restoration of blurred image, AP186 2009

Activity 18 Noise models and basic image restoration

2009-10-05T06:16:00.001-07:00

Just like any other measuring device or data gathering techniques noise is present in almost all forms of digital image acquisition. Noise, in general, is any unwanted aberration, random signal, or artifacts present in the image. Some of our previous activities already involved dealing with poor quality images using image enhancements. Previously we have used filtering in Fourier space, morphological operations, and even white balancing to improve the images. In this activity we will be using different image restoration methods to improve images with additive noise. In contrast to image enhancement, image restoration involves reversing the degradation process that produced the noise in the image. This means that a priori knowledge of the noise is important for a successful image restoration.

For this activity, we are tackling additive noise in the form of random variables following a specific probability distribution functions (PDF). We will be using six different noise models specifically: gaussian, rayleigh, gamma, exponential, uniform, and impulse noise.

Figure 1

There are many different types and methods of image restoration but in our case we will limit ourselves with spatial filters since they are most suited for additive noise. We will be using four types of filters namely, arithmetic, geometric, harmonic, and contraharmonic mean filters. With this filters we are basically changing the pixel values using the information provided by the image. All of this filters work by considering a subimage window centered at pixel(x,y) and using all the values within this window to calculate the new value of pixel(x,y) (figure 1). Equation 1 shows the formula followed by the arithmetic mean filter. This filter replaces the pixel value with the average of all the values in the window. On the other hand the geometric mean filter uses the product of all the values and takes its Ath root where A is the total area of the window (equation 2). The harmonic and contraharmonic mean filters follows equation 3 and 4 respectively.

Equation 1

Equation 2

Equation 3

Equation 4

In this activity we first create a noisy image by adding noise to a grayscale image then apply a restoration filter. We do this for each of the noise types and each of the restoration filters. Furthermore, to aide us in assessing our results we also take the histograms of all our images before and after restoration.

Figure 2 to 7 shows the result of applying the different restoration filters for each of the six noise types. It is clearly seen that even though the result is slightly more blurred it is actually much closer to the original image. This is further highlighted by looking at the image histograms. The histograms of the restored images are much closer that of the original image. The broadening or shifting introduced by the noise was greatly reduced by all the restoration filters. The minor blurring is simply due to the windowing; that is a smaller window results in less blurring but sacrifices the efficacy of the restoration. Here we used a 5x5 pixel window for a 256x256 pixel image.

Figure 2 Figure 3

Figure 4 Figure 5

Figure 6 Figure 7

However, looking at figure 7, we see that the restoration processes don't seem to negate the effect of the impulse (salt and pepper) noise. comparing the images we see very little to no improvement after implementing the arithmetic and geometric mean filters. Furthermore, looking at the histograms we even see that the restoration process degraded the image. The only semblance of improvement we can see is that the harmonic filtering and contraharmonic filtering with negative Q were both able to completely remove salt noise but not pepper. On the other hand using the contraharmonic filter with positive Q removes pepper noise but not salt.

We also did this whole procedure for a real grayscale image obtained from the internet (figure 8).

Figure 8

Overall this activity has been very enlightening and the results were very convincing. We were able to show that the different spatial filters are indeed useful for image restoration except for salt and pepper noise. I give myself a grade of 10 in this activity.

Main Reference:
AP186 "A18 – Noise models and basic image restoration", 2009

Activity 17 Photometric Stereo

2009-09-09T18:29:00.002-07:00

Many of our previous activities involved using an objects reflectance, its interaction with a light source, to do color segmentation, white balancing, and even classification. In this activity we will again use the this phenomenon to reconstruct a 3D surface. The method we will be implementing is called photometric stereo.

To give us a better understanding of our method it is best to backtrack a little. Let us consider a camera, such as those we have benn using in our previous activities. Esentially, the image obtained by a camerais simply the light intensity, I, reflected by its subject which is dependent on the materials reflectance. However aside form its reflectance this intensity pattern is also dependent on the geometry of the object which results in light and dark (shadow) areas. This is much more understood by looking at figure 1, where we can see that the reflected intensity at point P is affected by the angle at which the light hits the surface, resulting to either a bright area and/or a shadow. Hence an intesity profile of an object also contains information on its geometry. Equation 1 gives an expression for this intensity where rho is the reflectance and k is a scaling parameter.

Figure 1. Visualization of Reflection on a 3D surface

Equation 1

Looking at our model on figure 1 and equation 1, it becomes apparent that we should also take into account the type or our light source. That is, point sources, line sources, and area sources would result in different reflected intensity patterns since the light coming from them would arrive at the surface differently (different angles and intensities as a function of distance). However in this activity we limit ourselves to an approximation that our point source is at infinity which results to equation 1 having the form

Equation 2.

Since we know the intensity profile from our image, if we have information on the location of the light source we can obtain the geometry of the surface by solving for the normal vector n.

In this method, photometric stereo, we will be using multiple images of the same object but with light sources at different locations to obtain an estimate of an objects 3 dimensional geometry. Since we have images for different light source locations we can rewrite equation 2 in matrix form where each row corresponds to a image (intensity profile) light source pair. Solving for the normal vector n we obtain

Equation 3

Equation 4

We also know that this normal vector n is related to the surface c by the taking the gradient of c such that

Equation 5

Finally we obtain the surface f(x,y) by integrating

Equation 6

We implement this algorithm in Scilab and make use of matrix properties and operation to reduce the computation time.

Figure 2. Images of the surface taken with 4 different location of the light source. Their respective locations were previously provided.

Figure 2 shows the raw images we used to calculate the 3D surface. We used 4 images with sources at different locations from the matlab file photos.mat. It is easily seen from this images that indeed there is a noticeable difference with different light source position.

Figure 3. 3D plot of the calculated surface of the 3D object. Errors are seen as spikes and dips over the surface

Using the algorithm we discussed and the raw images shown in figure 2 we reconstruct a hemispherical surface shown in figure 3. The resulting surface we obtained is in good agreement with our expectation. Notice that the shape is indeed that of a sphere and the background remains flat. However we also see that the surface is very rough and there is even a gentle dip in a the shape of a cross over the hemishpere. These imperfections is inherent to our method, since it is still only an approximation. A much better result may be obtained by having more images (hence more information resulting to a better approximation). Still overall we can say that the method is a success and is fairly simple to implement.

I give myself a grade of 10 in this activity. I thank Dr. Gay Pereze for her patience and useful adice.

Main Reference

A17– Photometric Stereo by Dr. Maricor Soriano, 2008

Activity 16 Neural Networks

2009-09-09T18:29:00.001-07:00

Do you remember when I said machines are stupid? Well, I am not the first or the only one to believe so. Genius scientists and programmers that came before us also believed that machines and computers are sorely lacking when compared to the human brain. Just like in our pattern recognition activities (14 and 15), after creating our feature vectors we still needed to implement a formula so that the computer would classify the objects. So the brilliant scientists and programmers of the 19th century developed the “Neural Network”. The “Neural Network” works by mimicking how our own brain works; that is each neuron takes in multiple weighted inputs, processes the information, and then passes its output to a succeeding neuron. A collection of many interconnected neurons forms a “Neural Network.” The main advantage of using “Neural Network” over LDA and minimum distance classification is that we do not need to apply various classification rules and formulas.

Although I said that we do not use formulas with “Neural Network”, strictly speaking we still do. However this various calculations happen or exists within the structure of the network and the user won’t ever encounter it. We use “Neural Network” by simply giving it examples of a class and telling it that that particular example is a member of a specific class. For a given amount of time or iterations the “Neural Network” will adjusts the input weights of each neuron to minimize its error (identifying an object as a member of the wrong class). After all the iterations the “Neural Network” has already “learned” and will identify an object as a member of any class it has learned.

In this activity we will use Neural Networks from the ANN toolbox of SciLab to do pattern recognition. Our input to the neural network is the feature vectors we have from activities 14 and 15. We train the neural network by giving it 9 examples of each of the 4 classes and telling it which class each object belongs to (Activity 14 Table 1). Also we test the accuracy of the “Neural Network” network by letting it classify 40 more objects, one for each class, and counting how many it classified correctly. We test the accuracy of classification for different learning rates of the “Neural Network” from 0.01 to 0.99.

Table 1. Classification result and percent accuracy at different learning rates.
(click to enlarge)

Table 1 shows the results and accuracy of classification of the “Neural Network”. We see that as expected as we increase learning rate the accuracy of classification also increases and reaches a maximum of 100% at a learning rate of 0.23 to 0.55. The red box over table 1 highlights the area with 100% accuracy. However, as we increase the learning rate beyond 0.55 the accuracy slightly decrease to 97.5%. This is not expected since a higher learning rate should result in better classification. The classification accuracy versus learning rate is plotted in figure 1 and this trend of increase and slight decrease is observed much better.

Figure 1. Classification accuracy versus Neural Network learning rate. Red line highlights range with 100% accuracy

It would have been better to do multiple trials at different seed values to determine the standard deviation of the classification accuracy for each learning rate. If we have this information we might find out that this decrease in accuracy is not significant and is within the standard deviation. It would have also been better if we explored the effect of number of learning iterations but this was not explored due to time constraints.I would like to thank Kaye for reminding me of what to do in this activity. I give myself a grade of 10.

Main Reference:
Maricor Soriano, A16 – Neural Networks, AP186 2008
Cole’s AP186 Blog

Activity 15 Probabilistic Classification

2009-09-09T18:28:00.000-07:00

In the previous activity we used minimum distance classification to assign the classes of our objects. However, even with our best choice of features, relying on the distances amongst feature vectors may still be ambiguous. This was not a problem in our previous activity but such a scenario is indeed possible. In activity 14 if we had just opted to use the normalized Red value and the Area/Perimeter^2 it would have been very difficult to differentiate between heart and diamond. In this activity we implement another method of classification called Linear Discriminant Analysis (LDA). We will still be using the same feature vectors that we described in activity 14 but this time we will use a different method for calculating an object’s class.

In pattern recognition we may think of the all the events as having a given probability. Let us consider two classes and a single object. Assuming that the object and the classes have similar features we may consider that the error or accuracy of our classification is given by the probability of assigning the object to the wrong or correct class. Moreover, the probability that an object came from a given class not only depends on its features but also the knowledge of how many objects belong to a given class. LDA is simply a method of calculating the probability of an object to belong to a given class using their features. The membership of an object is decided by the class that has the highest probability of containing it. The LDA method is based on the assumption that each class has a multivariate Normal distribution and that all the classes share the same Covariance matrix.

LDA calculates the probability, fij, of an object to belong in a class according to the equation

Equation 1

In this equation wi is the feature vector of object i and mj is the mean feature vector of class j. The value pj is the probability or ratio of the number of elements in class j over the total population of the classes.

Equation 2

Matrix C is the shared covariance matrix of all the classes and is calculated as the weighted mean of all the covariance matrices, cj, of the classes. We calculate the covariance matrix of each class following equation 3, where m is the mean feature vector of all the objects in all classes and Wj is the matrix containing all the feature vectors of all the members of class j. Finally nj is the number of membership in a given class.

Equation 3

Using the feature vectors we obtained in activity 14 we use LDA to classify the objects into the four classes of heart, diamond, spade, and club. We used the nine training objects in activity 14 to determine the covariance matrices and mean feature vector of the classes. All the raw data we used are shown in activity 14 Table 1 and Table 2.

Table 1. Resulting probability of each object per class and classification results.
(click to enlarge)

The result of our classification is shown on Table 1. Again we obtain a perfect 100% classification using this set of feature vectors. We however did not implement this method on the other 140 objects we used in activity 14. Still we believe that a very high classification accuracy will result if done so.

I give myself a grade of 10 in this activity. I want to thank Orly, Kaye, and Jaya for answering my questions.

Main Reference
Teknomo, Kardi. “Discriminant Analysis Tutorial”. http://people.revoledu.com/kardi/ tutorial/LDAKardi

Activity 14 Pattern Recognition

2009-09-09T18:27:00.001-07:00

Us humans already have the inborn instinct or skill that allows us to recognize or differentiate one thing from the other. And most of the time, our own observations and experiences are enough for us to learn which is which. However machines and computers do not posses this instinct (they don't have instincts); much like how we teach babies and kids, we must also exert effort to teach machines how to recognize various things it may encounter. The difficulty with machines is that, although their language is based on math, they are still much harder to teach compared to kids (machines are stupid). In spite of this, the growing demand to increase productivity through automation and the physical limitations of the human body (there is just so much our body can do) pushes us to rely on machines and computers for various tasks. This demand and need for automated machine recognition resulted to a field known as machine vision. One sub domain of machine vision is pattern recognition which, from the name implies, basically means having a machine recognize and classify/segregate objects of patterns. In this activity we will be using Scilab to do pattern recognition of objects from an image. Our method for classification we will be using is the minimum distance classification.

All pattern recognition problems start with having a set of patterns or objects with a finite number of classification groups or classes. Each group or class is differentiated from each other based on a set of identifiable features and parameters. Hence in all pattern recognition tasks the first step is selecting the parameters or features we base our differentiation from. Since we are dealing with images it is logical that we use visually obtained information as our features. These features may be, color, shape, size, area, etc. and we may use as many as we see fit. We just have to remember that we need to select features that are distinctive to our objects or patterns otherwise we wont be able to classify them. In our method we place these features into vector form. That for each object there corresponds a vector that contains all of its features (equation 1).

Equation 1

The next step is to train the computer to recognize the different patterns. In our method this simply means that for each class we select a few samples and we take the average of their feature vectors (equation 2). This mean vector serves as the representation of the whole class; that is, each class should have a unique mean feature vector.

Equation 2

The final step is to test whether our feature selection and training was sufficient. We do this by letting the computer classify the other objects (excluding those used in training). The accuracy of classification depends on how much of the objects were classified into their proper classes. The minimum distance classification, which we use in this activity works by following equation 3.

Equation 3

As in equations 1 and 2 wi is the feature vector of object i and mj is the mean feature vector of class j. Using this equation we simply calculate d for each object for each class. We then classify object i into the class that resulted in the highest value of d.

In this activity I chose the suites of playing cards as my object class. Since I don't have a camera, I scanned the cards and this served as my raw data. From this images I classified whether the pattern is heart, diamond, spade, or club. Take note that I only consider the figures as my object and not the actual playing card (in other words the 10 of spades has 10 spade objects).

Figure 1

I used three features to construct the vectors describing each object. The first feature is their normalized R value (see activity 12). As expected the R value of diamond and heart is the same and so is the value of spade and club. Diamond and heart had a very high R value whereas club and spade had a low one. This made it easy for me to segment these objects from the background (figure 1) and calculate the other features I used. However, the segmentation was still not perfect so I used closing and opening (activity 8 and 9) operations to further improve the quality.

The next 2 features I used both characterizes the shape of each object. I calculated the ratio of the total area of each object with the area of the minimum rectangle that can contain the object which served as the second feature. The area of the minimum rectangle is simply the product of the maximum dimensions of the object ([max(x)-min(x)]*[max(y)-min(y)]) The final feature was simply the ratio of the object's total area with the square of its perimeter.

Table 1. Features and Group of Training set
(click to enlarge)

Table 1 contains all the features and group classification of the 9 samples per class I used to construct the training set and the mean feature vector of each class. I also labeled each group such that grp 1 are hearts, grp 2 diamonds, grp 3 spades, and grp 4 clubs. Now using the feature vectors of the other 40 objects in figure 1 (10 per suite/class) I calculated the corresponding d for each class. Table 2 shows the result of classifying these 40 objects and we clearly see a 100% accuracy of classification. This is further highlighted in figure 2 where I plotted the feature vectors of each object. Although the plot of Area/maxDim vs R could not distinguish between heart and diamond, it was compensated by the other parameters as seen in the 3D plot (figure 2d).

Table 2. Parameters and Classification Results of the 40 Test Objects
(click to enlarge)

(click to enlarge)
Figure 2. a)Area/max Dimensions vs R, b) Area/Perimeter^2 vs R, c) Area/max Dimensions vs Area/Perimeter^2, d) 3D plot R (x axis), Area/max Dimensions (y axis), Area/Perimeter^2 (z axis). The classes heart, diamond, spade, club are represented by the red, green, blue, black markers respectively.

Having tested the accuracy of our method I further tested its flexibility using larger images with mixed objects (figure 3) with a total of 140 objects. The results shown in Table 3 and Figure 4 still shows a 100% accurate classification.

Figure 3.

Table 2. Classification Results of the 140 Test Objects in figure 3
(click to enlarge)

(click to enlarge)
Figure 4. a)Area/max Dimensions vs R, b) Area/Perimeter^2 vs R, c) Area/max Dimensions vs Area/Perimeter^2, d) 3D plot R (x axis), Area/max Dimensions (y axis), Area/Perimeter^2 (z axis).

I give myself a grade of 10 in this activity for a job well done. I would like to thank Ms. Cindy Liza Esporlas for helping me scan the cards and sending me the images. I also would like to thank Mr. Jay Samuel Combinido for teaching me to do 3D plots.

Main Reference

Maricor Soriano, A14 - Pattern Recognition, AP186 2008

Activity 13 Correcting Geometric Distortion

2009-09-09T18:26:00.000-07:00

Have you ever tried looking through one of those door peepholes, the ones that allow you to see who is knocking on your door? I, myself, have only tried it ones when I visited my cousin who lived in a hotel and the door of their room had one. Anyway, for those of you who have tried using one, I am sure that you have also noticed that the image you see through it is highly distorted. This geometric distortion is similar to that found in some cameras (in a much lower degree) and other instruments that use spherical lenses. Basically geometric distortion results from the varying magnification as the distance from the optical axis increases. The two most common forms of this optical abberation is the barrel (magnification decreases farther from the optical axis) and pincushion distortion (magnification increases farther from the optical axis). In photography barrel distortion is commonly encountered when taking wide angle photos. In this activity we will be correcting geometrically distorted images numerically.

Figure 1. Example of distorted images

Our method of correcting geometric distortion involves two main steps. The first step is correcting the pixel coordinates, which refers to mapping the coordinates of the distorted image back to the ideal undistorted image. Our experience with image processing allows us to understand that the process of distortion simply involves the ideal undistorted image and a distorting transfer function such that

Equation 1

Where T is the transfer function and f and g are the undistorted and distorted images respectively. So given the distorted image we just need the transfer function to recover our ideal image. However the process to determine the transfer function requires us to have information on the undistorted image. Our method allows us to calculate the transfer function using only vertices of polygons or regularly occurring patterns present in our image.

Figure 2

Consider figure 2 above where we show an undistorted and distorted polygon/grid. assuming that the area of the polygon is not to large we can approximate the transformation as

Equation 2

The 8 unknown coefficients in this transformation can be solved using the for vertices of the polygon such that in matrix form

Equation 3

This resulting transfer function is valid for all points with in the polygon as long as it is not to large.

The ideal grid required in this method can easily be generated using the least distorted area in the distorted image. The dimensions of the least distorted "box" can be used to construct the ideal grid. We would then just need to determine the transfer function for each pair of distorted and undistorted "boxes." Using these transfer function we can now solve equation 2 to determine the distorted coordinates of each point in the ideal image. That is we can now determine the graylevel values for each ideal undistorted image coordinate by copying the graylevel value of its corresponding distorted coordinate. This is the second step of our correction process.

However, most of the time our calculated distorted coordinates are not integers. For such cases we need to interpolate the graylevel value for our ideal undistorted image. Three commonly used interpolation techniques are the nearest neighbor, bilinear interpolation, and cubic convolution. In this activity we choose to implement bilinear interpolation which has a much better quality than nearest neighbor and is much faster compared to cubic convolution.

Equation 4

Equation 5

Bilinear interpolation works as in equation 4. That is, the greylevel v is assigned to the coordinate (xhat,yhat). so by solving equation 4 we can obtain the gray values even for non-integer coordinates. The unknown coefficients of equation 4 can likewise be determined similar to equation 3 but know we use the four nearest pixels around (xhat,yhat) (figure 3, equation 5).

Figure 3. Visualization of Bilinear Interpolation

In this activity I used an image with barrel distortion downloaded from the internet. Figure 4 shows the original image and the ideal grid I generated using the least distorted area in the image. After implementing the method we have discussed I was able to remove the obvious barrel distortion of the image.

Figure 4. (left) Original distorted image with indication of corrected area (red box) and area of least distortion (yellow box). (center) Ideal grid generated from area of least distortion. (right and bottom) Resulting corrected image.

The result presented in figure 4 shows that our method is very effective in removing barrel distortion. It is clearly seen that the previously curved lines are now vertically and horizontally straight. The bricks that previously had different sizes are now more uniform. The only inaccuracy we can see is that some of the bricks are slightly tilted. It is hard to determine whether they should be tilted or not but I believe it is safe to assume that most of them are squarely laid-out. Still, even with this slight error, we see very promising and convincing results.

Figure 5. Corrected image with color. The red circles identify areas with erroneous color.

Figure 5 shows a colored version of our result. This was generated by implementing the correction process for the 3 color channels and combining them into one image. This colored result shows one more error. We clearly see that at some points have colors that are obviously not accurate. This is because at these points one or two of the three color components are not properly represented. Using the same reference points in all the color channels may result in some points having the wrong color component values. In other words although the gray scale result looks almost perfect it cannot be directly implemented to a colored image and expect the same level of accuracy.

I thank Ate Cherry for helping me understand the correction procedure. I give my self a grade of 10!!

Image source:
http://www.nicovandijk.net/lenses.htm

Main Reference:
A13 – Correcting Geometric Distortions by Maricor Soriano, AP186 2009

Activity 12 Color Image Segmentation

2009-08-06T09:11:00.002-07:00

So far we have been doing image segmentation and morphology on binary and gray-scale images. In the past activities simple thresholding is enough to seperate our region of interest (ROI) but in some cases this wont be sufficient. In image processing there instances where binarizing an image or converting it into gray-scale wont separate it from the background, in such cases we need to utilize its color information. In this activity we use normalized chromaticity coordinates to do segmentation via color information.

Figure 1. Normalize Chromaticity Coordinates

The normalized chromaticity coordinates (NCC), shown in figure 1, can be thought of as a way of expressing 3 dimensional RGB information to a simpler 2 Dimensional r-g color space. The basis for this transform is that r+g+b=1 (done by normalizing RGB), hence the third value is redundant and it is enough to express color by only two values, in this case r and g.

After expressing an image in NCC we can now do the segmentation by using the probability distribution function (PDF) of the ROI in NCC. By mapping the image to the PDF of the ROI we are able to separate the ROI from the rest of the image. In this activity we calculate the PDF of the ROI using two methods, one of these is the parametric method which assumes that the PDF is a gaussian. In the parametric method we simply obtain the mean and standard deviation of the r and g of the ROI and use this to calculate its 2 dimensional PDF. The other method is the non-parametric method, here we do not assume any form of the PDF and we simply use the r and g (2 dimensional) histogram of the image and use this as the PDF. In both methods we obtain the PDF of the ROI by only taking a patch of it.

I used an image of multicolored candies to test these two methods for color image segmentation. Excluding the background, there are only three main colors in this image; red, orange, and green.

Figure 2 below shows the results of the segmentation using both parametric and non-parametric methods. For both methods separating the green region was accomplished almost perfectly. The pixelated result for the parametric method is caused by the actual texture of the candy and the parametric method wont be able to compensate for such a feature. On the otherhand the holes on the result of non-parametric method are most likely due to the patch I used to sample the ROI. This simply means that there are other colors which are part of the ROI that were not included in the sample patch. Using a better sampled patch would improve the result in this case.

Compared to the non-parametric the parametric method did not fare as well when segmenting the red or orange regions of the image. Since the parametric method assumes a gaussian PDF it is higly probable that the PDF of the red or orange ROI includes each other. Hence the parametric method fails when we try to separate colors that are very close to each other.

The non-parametric method does not suffer from the same problem as the parametric method since it uses the actual histogram of the ROI as the PDF for mapping. This means that even if the two colors are very close to each other there wont be any confusion as long as the other color does not appear frequently inside the ROI. As we can see from the histograms the orange and red are very close to each other but their histograms are almost mutually exclusive.

Figure 2. Original unsegmented image with the results of both parametric and non-parametric color segmentation. The patch of the ROI is also included as well as its histogram.

Overall, by looking at the results, I have to say that the non-parametric method is a much better color segmentation technique. Compared to the parametric method, the non-parametric is much more versatile, accurate, and the quality is also better.

I give myself a grade of 10 in this activity.

Activity 11 Color Image Processing

2009-08-06T09:11:00.001-07:00

We have learned in our Applied Physics 187 course that an image captured by a color digital camera is composed of an array of pixels each having R, G, and B information. These three color information are each a product of the spectral profile of the light source, objects reflectivity and the cameras spectral reflectivity for R, G, and B. To obtain the actual appearance of the object we must divide these R, G, and B values with a scaling/balancing factor equivalent to the R, G, and B values of purely white imaged captured by the camera. This scaling factor is just the product of the cameras sensitivity and the illuminator's spectral profile. If the wrong scaling factor is used or it is miscalculated then the resulting image would not appear as good. The process of applying these scaling factors to a digital image is called White Balancing. In this activity we try out the White Patch and Gray World algorithms in doing White Balancing.

To test the effectiveness of these algorithms; we take pictures of objects containing the major hues which are obviously not properly balanced. We do so by not using the "auto white balance" function of the camera and just selecting the different balancing conditions. In the camera we used these conditions are called "daylight", "cloudy", "fluorescent", and "incandescent".

The White Patch algorithm makes use of information "obviously white" regions in the image. That is we calculate the balancing factor by taking the average R, G, and B values within a white patch of the image. On the other hand the Gray World algorithm assumes that the average color of the world is gray. Since the ratio of the R, G, and B values for gray is the same as white or black, this method takes the RGB of gray and uses it for White Balancing. In the Gray World algorithm the RGB values of gray is taken as the average R, G, and B of the whole image (whole world).

We must remember though that these White Balancing Methods wont work if the captured image is saturated, even for just a single color channel. A saturated image actually contains no color information and would not be useful in calculating the balancing factor.

Figure 1 below shows samples of captured images and their corresponding RGB channels. The blacked out regions on the RGB channels indicate the areas of saturation and wont be useful for our purposes. Therefore, in selecting the white patch we create a composite image (figure 1 last column) using the RGB channels that blacks out all the areas that are not useful. For the Gray World algorithm we simply take out the saturated regions in calculating the white balancing factors.

Figure 1. Sample captured images with their RGB channels and a composite image
indicating (blacked out) the areas that saturated.

Taking into account the saturated regions the results of applying the White Patch and Gray world algorithms on an image containing the major hues (RGB) are shown in figure 2 . First let us note that the results for both White Patch and Gray world algorithms for the incandescent and fluorescent settings are obviously wrong. The resulting image takes on an extreme shade of blue (or red) which is clearly not the true color of the objects. This error is mainly due to the fact that the original images taken with these settings were mostly saturated on the blue channel. On the other hand the images taken with Daylight and Cloudy settings were seen to improve by using the White Patch algorithm. In the daylight setting the original image is slightly bluish which after processing took on more accurate colors of the object.s This was also observed in the cloudy setting which initially had a reddish or brownish that was cured after processing. As for the Gray world algorithm, it was only seen to improve the images taken with cloudy setting resulting to a final image that displayed white as "white" even better than the White Patch algorithm.

Figure 2. Original image taken with different white balancing settings and
their respective results after applying White Patch and Gray World.

Next, figure 3 shows the results of applying White Patch and Gray World on images that had a dominant color, which in this case is blue. Again it is observed that applying a white balance on the images taken with fluorescent and incandescent settings resulted in a dark image. But even though these processed images are darker they are stilled considered improvements on the original since the original pictures are actually very bluish are actually worse representations of the true color of the image. As for the results of the Daylight settings we first note that the original image actually seems to be properly white balanced. We also see very little difference in comparing the original with the result after the White patch algorithm. But looking at the results of the Gray world algorithm we see that the image lost its yellowish color after it was processed. It is though to determine which is the more correct result. Finally the image taken with the cloudy setting has a very yellowish color which is obviously not correct. After doing the White patch we see that the yellowish color is removed which is a significant improvement. But after doing the Gray world algorithm the final image completely lost its yellowish shade but also became darker.

Figure 3. Original image taken with different white balancing settings and their respective results after applying White Patch and Gray World. The images have blue as the dominant color

Overall the results suggest that in general the White Patch algorithm results in a much better white balancing. The results of the gray world algorithm is seen to be highly dependent on the dominant color of the image.

I thank ate Cherry for lending us her camera and Irene for taking the pictures. I give my self a grade of 10 in this activity.

Activity 10 Preprocessing Text

2009-08-06T09:08:00.001-07:00

In our past few activities we learned various ways of manipulating images both in the frequency domain and the morphologically. This time we will use all that we have learned to extract handwriting and other texts from a scanned document.

Figure 1. (Left) Original document and (Right) the schematic
on how the tilt angle was calculated.

Looking at the given document we first notice that it is slightly rotated so our first task is to straighten the document. We do this by using the 'rotate' argument of the mogrify() function in Scilab. This is a function that edits or distorts an image in numerous ways, in this case we do rotation. But before we do the rotation, we must first determine the angle of tilt of the image. We measure this angle by creating a horizontal line across the lines of the document and we calculate the angle by simple trigonometry (figure 1). And doing this method the angle we calculated was 1.05 degrees. The results of the rotation are shown on figure 2.

Figure 2. The result of rectifying the whole image and also a selected area.

We also note that the horizontal lines of the image are actually not needed and that it interferes with our objective. To remove these lines we must first invert our image (white over black) and convert it to gray-scale then apply a filter or mask in inverse space. This is similar to what we did in activity 7 and since we already know that the FT of horizontal lines would be along the vertical we simply create a mask that removes this portion of the images FT. After we apply the mask we then binarize the image to remove the background. The mask and resulting images are shown on figure 3, Fourier transform of the image is also included.

Figure 3. (top) The inverted image with its Fourier transform and (below) to the right is the mask used to filter the lines. To the left is the resulting binarized image after filtering.

After removing the lines we now go on to extracting the handwritings. We only cut a select few regions in the image where we will extract the handwritings. Our first step is to binarize the selected regions using a threshold just like in the past activities. Then we did some morphological operations hoping that it would result to a clear, readable, and distinct letter characters separated from each other. Then we label each character using bwlabel(). Unfortunately this task is very very difficult to accomplish. The handwriting through out the document overlaps with unwanted structures. And since the handwriting also has a lower gray value compared to the unwanted structures separating them becomes even more difficult. Further more if we are not careful the morphological operations destroy the letters instead of making them more distinct. Figure 4 shows the results after morphology, here I used the skel() function of Scilab to reduce the characters to single pixel width. The results, however, show that the final product does not actually look like letters. Still it is also true that these are actually similar to the original handwriting and to some degree each character was separated from each other.

Figure 4. (Left) The slected area of the document with handwriting and (Right) the result by extracting the handwriting using different morphological operations.

Our last task is to find instances of the word description with in the document. We did this by taking a template of the word description and correlating it with the image. The correlation is done similar to how we did it in activity 5. We simply multiply the FT of the template and the conjugate of the FT of the image and then take the FT of this product. If the correlation is correct it would result to an image that have the highest value at the instances of the word description. Figure 5 shows the template and the result of locating the instances of the word description. By comparing the image with the result of the correlation we see that as expected the correlation peaks at the locations of the word description. This result means that our method for locating a template works even if the quality of the image is not ideal.

Figure 5. The binarized image and mask used to locate the instances of the word description. The result of the correlation indicating the instances is also shown at the bottom.

I give myself a grade of 9 in this activity.

Activity 9 Binary Operations

2009-08-06T08:33:00.001-07:00

In image processing it is often more efficient if we are able to separate our region of interest (ROI) from the rest of the image. The simplest way of doing this is by binarizing the image using a threshold. But even in binarized images having multiple ROIs offers some problems in the analysis. Overlaps among ROIs or noise introduced by insufficient thresholding would hinder the analysis of the image. Such problems can be addressed using different morphological operators, such us those we explored in the previous activity. The opening operator can be used to separate ROIs that overlap while the closing operator is used to fill in holes within the ROI. In this activity we make us of different morphological operations in improving an image of multiple disks/circles for area estimation.

To do the opening operation we first erode the image then dilate it. This would have an effect of separating objects that were initially adjacent or slightly overlapping. On the other hand the closing operation first erodes the image then dilates its. The effect of this operation would be to fill in holes of solid objects.

Figure 1. Original gray-scale image of disks.

Our image is a gray-scale ensemble of circular disks with some that are overlapping (figure 1). In this activity we aim to estimate the area of a single disk by taking the mean and standard deviation of the measured areas. Our first step is to cut the image into smaller images to make it more manageable for Scilab. We cut the image in way such that the slices overlap, this increases our sampling and prevents loss of information. Next we binarize each slice to exclude the background and highlight the disks. Then we do opening and closing operations to separate overlapping disks and also remove the noise left after binarizing. The structuring element I used is also a circle but is much smaller compared to those in the image. Finally for each slice we label each white area using the function bwlabel() and then measure its area and store the information. The results for each step we did are shown in figure 2.

Figure 2. From left to right: a slice of the image, a binarized version of the slice,
image after morphological operations, and the gray-scale labeled slice.

After we have completed the area measurements we now express our results in the form of a histogram (figure 3 left). Looking at the histogram we see that there are area measurements that are extremely large or small. Intuitively we know that this are not correct measurements and are mainly due to disks that were not properly separated or those that were cut in when we sliced the image. These measurements are outliers and we can discard them when we calculate the mean and standard deviation of the area. Another histogram is shown on the right side of figure 3 excluding these outliers.

Figure 3. (Left) Histogram of all the measured areas and
(Right) histogram of the measured areas excluding the outliers.

We know calculate the mean and the standard deviation of the measured areas excluding the outliers. Doing so we obtain a value of 501.8585 for the mean and 70.24116 for the standard deviation. These values are well within the expected value since the disks from the raw image have a radius ranging from 12 to 14 pixels. The mean area is the best estimate we obtained for the area of the disks whereas the standard deviation provides the range at which we can consider the areas equal or the objects the same. That is the area of the disks is 501.8585 ± 70.24116.

I will give myself a 10 in this activity. It was difficult but I think the results I obtained were accurate enough.

Activity 8 Morphological Operations

2009-08-06T03:45:00.000-07:00

Our previous activities dwelled on operations, enhancements, and other image processing techniques done in the Fourier space. This time our activity involves morphological operations to augment or distort binary images. In this activity we will be 'eroding' and 'dilating' binary images of a square, triangle, circle, hollow square, and a cross using different structuring elements. But before we do this operations in Scilab we will first do some manual, 'pen and paper', predictions.

A structuring element (SE) can be thought of as the distorting pattern used in morphological operations. The dilate operation results in an increase in the size of the dilated object whereas the erode operation decreases the its size. Both operations follow the shape of the structuring element. The result of a dilation can be thought of as the collection of points that the center of the SE can occupy such that the SE and the dilated object intersects. On the other hand the result of an erosion is the collection of points the center of the SE can occupy such that the whole SE is inside the object.

The structuring elements we use in this activity are a 4x4 square, 2x4 and 4x2 rectangles, and a cross 1pixel thick and 5 pixels long per arm. We apply the morphological operations to five shapes namely: 50x50 square, 50x30 triangle (base x height), radius 25 circle, 60x60 hollow square with 4 pixel thick arms, and a large cross 8 pixels thick and 50 pixels long per arm.

Figure 1 shows the result of the dilation of these shapes. My predictions for the square, circle, hollow square, and cross were all correct. I was able to predict both the shape and the actual size (pixel counts) perfectly. The square SE increased the sizes of these four shapes along the horizontal and vertical uniformly. Hence all retained their general shape except for the circle which became slightly squarish. As for the rectangle SE's they distorted the shapes more along their dominant axis. Finally, the cross SE caused the shapes to have a hole on the corners similar to the cross whereas the circle slightly had a greater increase along the horizontal and vertical.

I only made a mistake in predicting the result of the triangle. I got the proper shape right but the actual dimensions were off. I was confused on how the edges and the corners of the triangle would result. Compared to the results my predictions were slightly smaller.

Figure 1. Result of the dilation of five shapes (left column)
with four different structuring elements (top row).

Next, figure 2 shows the result when the shapes are eroded. Again my predictions were mostly correct but this time I was only wrong for the circle eroded with the cross SE. The results for erosion using the square and rectangle SE were simply the opposite of their dilation. The interesting result here is for the cross SE on the hollow square and cross. As I predicted only the corners of the hollow square can fit the cross SE, hence its edges would vanish upon erosion. As for the cross SE on cross it is interesting to see that near the center the resulting shape is thicker. This is because near the intersection of the cross there are more points where the cross SE can fit.

Figure 2. Result of the erosion of five shapes (left column)
with four different structuring elements (top row).

The last part of this activity involves exploring the thin and skel functions of Scilab. Figure 3 shows an image included in the SIP package which we chose to apply thin and skel on. Figure 4 shows the result of thin, just like in its description thin removes the edges of a figure until it becomes a skeleton of the original image. But the result seen in figure 4 is not actually a good skeleton of figure 3, which is part of the limitations of thin. Finally figure 5 shows the result of using the skel function on figure 3. skel is a better function that thin because (1) you can controll what type of algorithm is used and (2) both exterior and interior skeletons can be obtained. The output of skel is also in gray-scale which means we can controll the sensitivity of the resulting skeleton using im2bw() with variable threshold. Further more skel is an actual function designed to take the skeleton of an image not like thin.

Figure 3. Image of Birds (Escher)

Figure 4. Result of applying thin on figure 4.

Figure 5. Result of applying skel on figure 4, taking both
exterior and interior skeletons and using exact euclidian algorithm.

Information on the functions skel and thin is included in the help and documentation of Scilab.

I give my self a grade of 10 in this activity.

Activity 7 Enhancement in the Frequency Domain

2009-07-20T18:29:00.001-07:00

In the past few activities we familiarized ourselves with the use and applications of the 2-dimensional Fourier Transform. Along with the basic properties of the 2-dimensional FT,we also learned and observed the effects of imaging with convolutions and correlations. When we think about it, the past few activities were exercises geared to build on our skills for manipulation in the inverse space or frequency domain. In this activity we utilize the skills we learned to do image enhancements in inverse space.

Most of the time, when we do image processing and analysis, it is more appropriate to do manipulations in the frequency domain. For example, unwanted repetitive patterns over an image would be very difficult to remove without loss of important information. Such tasks can be easily handled when by going into the frequency domain. As we have seen in the previous activities, repetitive patterns are just represented by a few points or lines in inverse space and removing them becomes much simpler. In this activity we will be doing similar tasks in enhancing images of fingerprints, paintings, etc.

The first part of this activity is again aimed to further familiarize ourselves with inverse space. Here we compare the FT of different structures, namely: single points, circles, squares, and gaussians. We can see on Figure 1 below shows the Fourier transform of two single pixel dots compared with the Fourier transforms of two circles of different radius. As expected the Ft of the single pixel dots are multiple parallel lines just like a sinusoid in activity 6. On the other hand the FT of the two circles are very similar to the multiplication of the FT of the two dots and that of a circle (airy disk). When we think about it the image of these two circles is just like the convolution of the two dots and a circle of the same radius. Hence in inverse space the result should be the product of their FT, as seen in figure 1. Increasing the radius of the circle simply results to a decrease in the size of the FT pattern, similar to our results in previous activities.

Figure 1. 1st column: Two single pixel dots and its FT.
2nd-5th column: Two circles of varying radius with FT

Figures 2 and 3 below are just like figure 1 except that the objects are squares and gaussians respectively. Again as expected the FT patterns are simply the product of the FT of the two dots and the objects. That is the FT of a square is two perpendicular sinc patterns and for a gaussian the FT is also a gaussian. The behavior of the FT with varying object size is still the same (inversely proportional).

Figure 2. 1st column: Two single pixel dots and its FT.
2nd-5th column: Two squares of varying side length with FT

Figure 3. 1st column: Two single pixel dots and its FT. 2nd-5th column:
Two gaussians of varying spread (standard deviation/sigma) with FT

Next we explore the difference in the FT of an inverted gaussian. Figure 4 shows the gaussian and inverted gaussian patterns with different standard deviations. When we look at their normalized FT we see no difference at all so instead we compare their Real and Imaginary parts.

Figure 4. Double Gaussian pattern of different spreads and their inverted counterparts.

Figure 4 and 5 shows the comparison of the real and imaginary parts of the gaussian and inverted gaussian FT's. First let me say that though the images look like they are just flat gray images there are actually fine detail around the center which can be seen when the images are enlarged. Still even with very close inspection there are hardly any difference with the gaussian and inverted gaussian. But this is actually expected since inverting the gaussian simply involves multiplying by -1 which can be factored out of the transform. This simply results to a switching of the high and low values of the FT Real and Imaginary parts which is, in this case, hardly distinguishable.

Figure 5. Real part of the Fourier transform of the gaussian and inverted gaussian. The large image corresponds to the first gaussian and inverse gaussian pattern (figure 4) and the other three pairs are for the succeeding gaussian patterns.

Figure 6. Imaginary part of the Fourier transform of the gaussian and inverted gaussian. The large image corresponds to the first gaussian and inverse gaussian pattern (figure 4) and the other three pairs are for the succeeding gaussian patterns.

The second part of this activity involves enhancing a fingerprint image using inverse space filtering. We do this by applying a mask to remove the 'unwanted' regions in the FT of the finger print. Then we take the inverse Fourier transform of the masked FT and hope to obtain clearer and sharper fingerprint pattern.

Figure 7 shows the original fingerprint image and its Fourier transform. Looking at the original image we see that aside from the fingerprint a background pattern is present. This background pattern is probably the grain of the paper which only interferes with the fingerprint pattern. We will remove this background by applying a mask on the FT of the image. To create the mask we must first observe the FT pattern and determine which regions most likely contribute to the fingerprint. Once we have determined this region we will create a mask that only retains this part of the FT.

Figure 7. Original fingerprint image and its FT.

We can see from the FT in figure 7 that the center, having a much higher magnitude, contains most of the information of the image. Removing the center would result to losing not only the background but virtually the whole image. We also note that there are also two curved regions in the FT that are brighter than the other areas and most likely contain information on the fingerprint. Using these observations we now create a mask that enhances this regions hoping that it would result to an enhanced fingerprint.

Figure 8 below shows three different filtering masks I used, hoping to enhance the fingerprint. The first mask was made by simply blacking out all unwanted regions in the FT. The second mask was generated by creating two gaussians and rotating them about the center. Finally the last mask is a convolution of 2 gaussians with the first mask. Each of these masks highlights the regions I mentioned above and should produce an enhanced fingerprint. But as I found out using just one of these masks did not enhance the fingerprint enough to be considered a successful enhancement. So after numerous trials and even more errors I ended up using a combination of these three patterns, multiplying them with different scaling factors, to produce a good enough enhancement of the fingerprint (figure 9).

Figure 8. Three different filtering mask used together to
obtain an enhanced fingerprint pattern.

Figure 9 is a comparison of the original fingerprint and the enhanced one. The left image is the original fingerprint which was inverted and cropped. The right image is the resulting enhancement after the left image was filtered with the masks in figure 8. Comparing the two images we first note that in the enhanced fingerprint the background paper grain patter is not visible anymore. Moreover the text on the background have been significantly reduced. The most notable improvement though is that the fingerprint itself is much clearer. The lines are thinner and sharper, the overall fingerprint is also more uniform in color, and some lines are now distinguishable. But looking at the right edge of the fingerprint we can see some loss. A small part of the lower right edge of the fingerprint can not be seen clearly. Still this is a considerable improvement from the original.

Figure 9. Inverted original fingerprint image and the enhanced fingerprint image using frequency domain filtering.

The third part of the activity involves removing vertical lines found in a composite image of stitched photographs of the moon. This task is considerably much simpler compared to the enhancing a fingerprint. For this task it is clear what structure of the image we do not want and we simply need to create a mask that removes this in the frequency domain.

Figure 10 shows the original and gray-scale image of the moon and in this two images the unwanted vertical lines are clearly seen. But from our previous activities we know that this array of vertical lines is represented by areas along the horizontal of the FT so we create a mask to filter this region out (figure 10).

Figure 10. Original moon image and below it is the gray-scale image.
To the right is the FT and to the lower right is the filtering mask used.

Using the mask in figure 10 we removed the vertical lines and the result is seen on figure 11. The left image in figure 11 shows the resulting image after the mask was applied. It is clearly seen that there are no more vertical lines present in the image. But the final image is also seen to be darker than the original and the lighting is also different. This is probably due to some losses caused by the mask but was not enough to distort the detail in the photograph. Still the objective was achieved and there was very little lost information.

Figure 11. Gray-scale images of the original moon surface and
the resulting image after the mask was applied.

The last part of this activity is aimed to enhance an image of a section of a painting. The canvass weave is clearly seen in this image and distorts the brush strokes and the actual painting. Again we use a mask to remove this canvass weave and recover a clearer image. Figure 12 shows the original and gray-scale image of the painting as well as the FT. We know that since the canvass weave is repetitive then its image in the inverse space are probably symmetric points. So from the FT we identify areas that might be part of the weave and then create a mask to remove this areas. The mask I used is also shown in figure 12. We can see that the mask removes only selected symmetric bright points in the FT.

Figure 12. Original and Gray-scale image of the painting, the Fourier transform
and the mask used to remove the canvass weave.

Figure 13 shows the resulting image after the mask was applied. It is clearly seen that most of the canvass weave has been removed and the painting has much more distinguishable detail. But when we look closely there are still some traces of the canvass weave. The enhancement may not be perfect but the canvass weaved was sufficiently removed, enough to highlight the brushstrokes and details of the painting.

Figure 13. Gray-scale image of the painting and the
resulting image after applying the filtering mask.

Finally we look at the Fourier transform of the mask with the idea that it should result to the canvass weave, but as seen in figure 14 this is not exactly the case. The FT of the mask does resemble the canvass weave but the intensity is not uniform. An elliptic figure seems to be distort the canvass weave pattern. And also along the edges there are lines that cut across the uniformity of the pattern.

Figure 14. Inverted image of the mask and the FT. The FT mimics
the canvass weave pattern with some obvious distortions.

I got the fingerprint image using google image search but I forgot where (sorry). I give my self a grade of 9 in this activity since the results were not perfect.

Activity 6 Properties of the 2D Fourier Transform

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

In Physics, the Fourier Transform is widely used in various applications, including signal processing and analysis, modeling, and many more. The 2D Fourier Transform, in particular, is very useful in optics; calculations in the inverse space utilizes the Fourier Transform. In this activity we further improve our familiarity with the Fast Fourier Transform functions of Scilab and also the properties of the 2D Fourier Transform.

For the first part of this activity we generate some common shapes with familiar Fourier Transforms; figure 1 shows these shapes together with their 2D Fourier Transform. Note that the FT of a circle is an airy disk, while the FT of a double slit look very much like Young's double slit experiment. As for the square, we know that the Ft of a rectangular wave is the sinc function so the Ft of a square will resemble sinc functions along the x axis and y axis. Finally we can see that the FT of two dots is very much like a sinusoid (or cosine) pattern which is expected since the FT of a cosine is 2 peaks.

Figure 1. Different shapes or patterns generated using Scilab
and their corresponding Fourier Transform

The second part of the activity is aimed to familiarize ourselves with the properties of the 2D Fourier Transform using 2D sinusoid patterns. The sine function was chosen since, as I have said, the absolute value of the FT of a sinusoid (or cosine) is known to be 2 peaks.

Figure 2 below shows the 2D Fourier Transform of sinusoid patterns with different frequencies. As expected, the FT images contain 2 peaks vertically symmetric about the center. Also as the frequency of the sinusoid is increased we observe that the peaks move away from the center. This is also the expected result since low frequencies are concentrated on the center. Finally we see that with higher frequency the FT deviates slightly from a spot. For the higher frequencies the spots of the FT is less distinct and is slightly spread, which is probably due to aliasing.

Figure 2. Sinusoid patterns of different frequencies with their Fourier Transforms.
The frequencies are 1, 4, 40, from top to bottom.

For the next part we explore the effect of combining different functions, specifically a sine with an added constant bias and 2 sine functions with different frequencies. The base sinusoid we used has a frequency of 4 (figure 2 mid row) and we see that by adding a constant bias (figure 3 top row) that the FT only differs by the appearance of a central spot. This central spot corresponds to a 0 frequency signal which is exactly the constant bias. As for the 2 sinusoids, we simply added the sinusoid with frequencies 4 and 1. And again we see that the resulting FT is simply just the addition or superposition of their FT seen on figure 2. There for in general we can say that in 2D Fourier Transform the FT of an addition of 2 functions is just the addition of their FT.

Real life signals obtained from actual sensors usually contain a constant bias. For example, an interferogram as detected by a photodiode most likely rides above a bias. In trying to recover the frequencies of such a signal by taking its Fourier Transform, you would sure obtain a high valued zero order or zero frequency term just like in the top row of figure 3. If we want to recover the signal we can just ignore the zero order. But if it is too strong that it floods out the signal it would be better to just center your signal at zero.

Figure 3. (top) A sinusoid of frequency 4 with a constant bias and
(bottom) addition of sinusoids with frequencies 1 and 4. Beside them are their FT.

In figure 4 we see the effect of rotating the sinusoid pattern. Looking at the images we see that the FT is still 2 peaks except that it is also rotated. And just like in figure 2 the location of the peaks are symmetric about the center, also the line joining the 2 peaks is perpendicular to the pattern (same with figure 2).

Figure 4. Rotated sinusoid and its Fourier Transform.

We know look at the effect of combining a sinusoid pattern along the horizontal and the vertical. The first result is a multiplication of the two patterns, which resulted in a checkered pattern. The resulting FT of this pattern have 4 peaks, still symmetric about the center, but looks to be located at the corners of a smaller square. This makes sense since the pattern actually looks like an addition of 2 perpendicular sinusoids (horizontal and vertical) rotated by 45 degrees. In other words we see a pattern with a periodicity that is angled by 45 degrees. The next set of images is just an addition of a vertical and a horizontal sinusoid patterns.

Figure 5. (top) Multiplication and (bottom) addition of a vertical
and horizontal sinusoid with their Fourier Transforms.

Finally we make a pattern by adding multiple sinusoid patterns each rotated by an arbitrary angle. The resulting pattern doesn't seem to mimic the previous pattern we have presented but still looks to have some periodicities and symmetries. Looking at the FT of this pattern we can say that it is just like an addition or superposition of the FTs of the different rotated sinusoids. This is observation is confirmed by the circular appearance of the FT. So over all we can say this final image just reaffirm what we have already demonstrated in figures 3 and 4.

Figure 6. Addition of arbitrarily rotated sinusoid patterns and its Fourier Transform.

I give my self a grade of 10 in this activity. I also want to thank master, jica, and jaya for staying with me in the lab and for reminding me to do the blog.

Activity 5 Fourier Transform Model of Image Formation

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

This is an activity which explores the functionality of the built in two-dimensional Fourier Transform of Scilab 4.1.2, fft2(). We will also be showing how this function, along with the convolution and correlation theroems, can be used in various image processing and analysis tasks.

The first task is to familiarize a\ourselves with the fft2()function. Here we first constructed an image of a circle and the letter A (first row). Then using the fft2() and the fftshift() functions we are able to display both of their Fourier Transforms. We know that the results are correct by examining the images. We already know that the fourier transform of a circle or a point is that of an Airy disk and that a larger circle will produce a smaller Airy disk and vice versa. Such results are seen in the images bellow (middle row). Finally we take the fourier transform of the fourier transforms (last row). We know see that the final images are inverted with respect to the original. This is also a common result when using the fourier transform functions of other programming languages. Shifting and scaling constants should be used to be able to recover the original image.

Next we the fft2() function and convolution theorem to model an imaging device with a circular aperture. Here we see that when we use an aperture to image an object, the final image is slightly blurred or deformed. Also we only multiplied the FT of the image to the circular aperture and not its FT since the aperture is already at the focus or the inverse space. We can also see that a larger aperture produces a sharper image which is consistent with our previous results; that is a larger aperture will have a smaller Airy disk FT.

The third task is that we use correlation to locate a pattern with in a certain image. Here we have an image with texts "THE RAIN IN SPAIN STAYS MAINLY IN THE PLAIN" and another image with just "A" of the same text size. We used the fft2() function and implemented the correlation theorem to locate the occurence of the letter A in the texts. After doing correlation and taking its FT again we recover the same image but with locations of "A" highlighted.

The final task is to do edge detection using convolution. Our object is an image of "VIP" and we use different patterns in edge detection. The patterns we used are 3x3 matrices of horizontal, vertical and a point. The resulting images clearly show that the method most accurately detected figures in the image that matched the patterns.

Thank you Martin for talking.

I give myself a grade of 10 for a job well done.

Activity 3 Image Types and Basic Image Enhancements

2009-07-04T21:12:00.000-07:00

Isn't it surprising to discover how much we don't know about things we commonly use or encounter; just like in this activity for example. I am sure that most of us, if not all of us, are already aware of the different image or picture formats that are recognized by our personal computers. But to be honest I don't really know the difference between these formats or even how they are interpreted by our machines. In this activity I learned that there are just a few basic types of digital images. The simplest is that of a binary images which is suitable for black and white images since it can only store a single bit of information per pixel; ones or zeros, black or white. Next would be a grayscale image which commonly has 8bits of information, meaning it has 256 different gray values between black and white. The third are RGB images or true color images that stores about 3 times the information than a grayscale image, a value each for Red, Green, Blue. This image format imitates how are eyes interpret color as different combinations of Red, Green, and Blue. Finally there are the indexed images which are also colored images like the RGB. But instead of storing informations for the RGB channels an indexed image stores the image as index in a colormap. Hence it stores only 2 information; the image and its color map. As I gather the various image formats that we are more familiar with such as jpg, gif, png, tiff, etc. are just different ways of handling the these image types. Examples of these images are presented at the end of this post.

Also in this activity we employ the skills we obtained in Activities 1 and 2 to calculate the area of a real scanned image. For this activity I opted to create my own picture by writing my name on a piece of paper and scanning it with a ruler to give an idea of the actual size of the image. The resulting image is shown bellow together with its information obtained using the function imfinfo(). I drew a green rectangle over the scanned image to indicate my designated region of interest.

Original Scanned Image with Indication of the Region of Interest

FileName: miguel.bmp

FileSize: 372262
Format: BMP
Width: 834
Height: 444
Depth: 8
StorageType: indexed
NumberOfColors: 256
ResolutionUnit: centimeter
XResolution: 59.060000
YResolution: 59.060000
Size: 444 rows X 834 columns
Indexed Image

Just like in activity 2 I needed to convert this image to a binary black and white image in order for me to utilize the same algorithm in calculating the are. So to determine what the adequate threshold should be I made an algorithm to obtain the histogram or Probability Distribution Function of gray values for this image. As seen below, the resulting histogram is dominated by

//create histogram or PDF
hist = zeros([1:256]);
for i = 1:256,
hist(i)=sum(name==i);
end
hist = hist/total; //noramlize

the black background, which is expected by simply looking at the picture. The only light colored regions in the image are the ruler and my name. These areas are clearly indicated in the histogram, encircled by the green highlight, and are completely separated from the background. Using the information provided by the histogram I selected a suitable threshold indicated by the dashed red line.

To be able to calculate the area accurately I separated the ruler from my region of interest before I converted it to a binary image. This was simply done by selecting the range with in the image matrix that you are interested with. After this we are now ready to calculate the area of my name.

ROI = name((10:444),(85:820));
ROI = im2bw(ROI, 0.5);

Binary Image of the Region of Interest

Boundary of the Region for Area Calculation

By using the same algorithm used in Activity 2, I calculated the area of this image to be 56,140 pixels. Since we also scanned a ruler beside the image we can use the skills we employed in Activity 1 to estimate the actual physical area occupied by my name. Doing the same method of calibration as in Activity 1, I calculated that the calibration is 1mm/5.9 pixels or 1pixel = 0.169492 mm. From these calibration we can say that 1x1 pixel has an are of 0.028727 square mm. So the total area of my name is 56,140pixel x 0.028727 square mm/pixel =1612.755 square mm or 16.13 square cm. This is probably not a 100% calculation of area since we can see from the image that there are some black pixels inside my name. Also there can be an error introduced when the image was scanned which would affect the calibration. Still I believe that the combination of these 2 errors is very minimal.

Overall I think this was a fruitful activity. I was able to combine the things I previously learned in Activities 1 and 2. I think I deserve a pretty high grade for this one, maybe a 9 or a 10 plus a bonus since I was able to make my own histogram. Hahaha!

Binary Image

Grayscale Image

Indexed Image

True Color Images

Activity 4 Enhancement by Histogram Manipulation

2009-07-02T07:44:00.000-07:00

During my senior year in high school I chose this visual communication class as my elective, and one of our activities were photography. At that time we were required to use a classic SLR camera which was shared by 3 people per group. It was a really fun activity but there were also some very frustrating moments. I clearly remember times when our shots were either to bright or too dark. Since we were newbies, we had a very hard time getting the exposure time, flash, and aperture size correctly. To add to our difficulties, we could only see the pictures ones we had them developed. This resulted to us spending money on buying rolls and rolls of film which really burned our pockets. At that time I would really have appreciated anything that would give me the ability to improve my shots which did not include taking more pictures and buying more film. Here I am now, five years later, doing the exact activity that could have saved me money five years ago.

In this activity we are tasked to improve or enhance the contrast of grayscale images by histogram manipulation. We do this by obtaining the images histogram or probability distribution function(PDF) and, from here, calculating its cumulative distribution function (CDF). Then to improve the image we map its current CDF to the CDF of our desired PDF which we believe would improve the image. This mapping of the CDF is done in 3 steps:

We take the grayscale value of a pixel i and locate its value in the CDF, cdf(i).
We then match this CDF value, cdf(i), to the desired CDF and locate its corresponding grayscale value.
Finally we replace the existing grayscale value of the pixel with the new one obtained in step 2.

The final result of this mapping would be a new image which would have a CDF similar in shape to the mapping CDF.

There are also other simpler methods of enhancing an image, such us simple contrast stretching. I used this method before, during my old research topic, to fully utilize all 8bits. In this method we simply take the minimum and maximum values of the image and implement the algorithm:

Image = (Image - min(Image))/(max(Image) - min(Image)).

This method is considerably much simpler both theoretically and computationally but its use is very limited. It works by taking your range of values, from minimum to maximum, and stretching it so that you will use all 8bits or 256 grayscale values. This being the case you can only use this if the image doesn't already have a value of 256 and 0, otherwise you would have to do some thresholding and truncate the higher and lower values. This would then cause some loss of possibly useful information.

The images below are the results of the enhancements using CDF mapping. The first column shows the original image and its normalized PDF and CDF. Whereas the second and third column shows the result after CDF mappings on a linear CDF and logarithmic CDF respectively.

The first CDF mapping enhancement made use of a linear CDF which corresponds to a flat PDF. A flat PDF means that all gray values have equal probability, in other words a detector having such a PDF would have no bias on any grayscale value. According to our Proffesor, Dr. Soriano, research grade cameras would ideally have a linear CDF. But she also said that such a camera would usually show a dark image which would not appeal to most people. Still this makes sense since in research accuracy and calibration is important, so having a favored range of grayscale values would make it difficult for the researcher. The other CDF mapping made use of a logarithmic CDF because according to John C. Russ the human eye has a logarithmic response [1].

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Before commenting on the results of the enhancements I would first like to note that by looking at the plots of the CDF we are assured that the implemented algorithm works correctly since the resulting CDFs match the theoretically expected.

By examining figures 1-5 we can see that in general the linear CDF works very well as an image enhancing tool. As expected all of the 5 images had a better contrast compared to the original.; this observation is supported by the resulting PDF of grayscale values. A broader PDF means that the image makes use of a wider range of grayscale values which is could for higher contrast and detail. Such a result was achieved in all of the pictures and specially for figures 4 and 5. Looking at the images we can easily see that most of them did not only had a better contrast but also had more distinguishable detail. Only figure 2 did not fair as well as the others. Though the image had a much better contrast, some of the detail previously seen was lost after enhancement.

As for the logarithmic CDF mapping, it is clearly seen that it did not fair as well as the linear one. It is true though that the CDF indeed looked match a logarithmic curve, still the image contrast actually suffered after the enhancement. All of these 5 images became much darker than the original and most of their detial became indistinguishable. Also when we look at their resulting PDF we see that it is concentrated on the dark values just like the image. This result should have been expected since this logarithmic curve has a most of its values in the lower x axis, or the dark region. Probably if we stretch the logarithmic curve a little it might result to a much better enhancement.

I tried various types of nonlinear CDFs, such as parabolic, exponential, cubed, and cube root, hoping to find one that would enhance all of the grayscale images that I tried. Soon enough I realized that such a CDF does not exist, or at least it is beyond me. What I did learn though was that a CDF would enhance your image if it stretched your histogram in the range of your region of interest. That is if you have little experience, such as myself, it would be better to look at the histogram and determine which range corresponds to your ROI. Then look for a CDF that would stretch this region and this would amount to a very good ehnacement.

I would like to thank Dr. Soriano for giving me very useful insights. I also thank Dr. Perez for helping me understand what to do. Finally I thank Jaya, Master, and Jica for the discussions and the laughs while doing this activity.

I also thank and award myself with a grade of 10. Thank you very much... HAHAHA!

Reference:
[1] Russ, John C. "The image processing handbook" 5th Edition, Taylor and Francis Group, CRC Press 2006-12

Activity 2 Area Estimation for Images with Defined Edges

2009-06-28T01:29:00.000-07:00

When I first heard or thought about it, area estimation did not strike me at all as something complicated or even useful. I was already aware of how area estimation played a role in agriculture, agrarian reform, and other endeavors involving large land areas (such as land grabbing and proffesional squating... haha). But as Ma'am Jing narrated and enumerated all the instances where estimating an area was the key task in the creater scheme of things then I started to realize what a powerful tool it can be. It was quite an eye opener when she shared how the study and treatment of skin cancer heavily relied on area estimation which was done manually. I also remembered that my cousin (a marine science student) once told me that in marine science measuring the total area of leaving corals is a very important task. This made me realize that area estimation is not as simple as it once seemed. Moreover measuring irregularly shaped areas is difficult enough but teaching a computer to do it is even more challenging. In this activity we utilize scilab and its SIP toolbox to write an area estimation code using Green's theorem.

Green's theorem has been and is still used in various area estimation algorithms. There is even an old area estimating instrument called the planimeter (http://en.wikipedia.org/wiki/Planimeter), used in irregularly shaped areas, that is justified by Green's theorem. Named after British mathematician George Green, it is a mathematical relationship of a double integral and the line integral along its boundary (http://en.wikipedia.org/wiki/Green%27s_theorem).
Given that the partial derivatives of F1 and F2 exists and are continuous within the region R we can take F1 = 0, F2 =x and F1 = y F2 = 0 and obtain two line integrals which both gives the area of R. The average of this two, equation 2, gives a better measure of R and equation 3 is the discreet form if this integral, where N is the number of points along the boundary contour.
Based from this concept I wrote a program in Scilab making use of the SIP toolbox to access the image and the required information from it. In this basic implementation we make use of the function follow() to obtain the boundary contour of the area of interest. But follow() can only handle binary images so the function im2bw() is utilized to convert non-binary images. Also, as we can see from equation 3, the descreet algorithm for Green's theorem makes use of one more point over the length of the contour. The simple adjustment I made was to simply add the starting point/coordinate as the end point of the contour. This is a valid adjsutment since there is no extraneous information added and in fact this last point necessarily closes the contour which was initially open. Also I used pixel counting as the basis for checking; that is, I was to compare the area calculated by the algorithm against the total number of white pixels. This made sense since in this case all of the images were created using paint and none were real photographs. As it is the program was fit enough to run and implement Greens theorem.

The trial run of the program aimed to calculate the area of a square with one side measuring 178 pixels but, contrary to what I expected, it returned the wrong value. At first I suspected that the difference was that the algorithm did not account for the boundary and Dr. Soriano shared this thought, but I already saw from the values that it was not the case. Then I immediatley noticed that the calculated area was that of a 177 pixel sided square. After some serious thinking and contemplating I had an indea of adding half of the boundary contour. For the square area this was obviously sufficient since the only difference of 177pixel and 178pixel on a side is half of the perimeter, still I had a feeling that this would be valid for all the other shapes. After modifying the code and testing it on sqaures, circles, triangles, and various irregular shapes I was amazed that the program was 100% accurate each time (table 1). Even very crooked curves were not a problem.
Then I tried the program for more irregular shapes and saw a limitation. once there is a discontinuity to the black background the program could not handle it. This means that if there are two seperate white areas then the program would only see one of them. Also if the shape was a donut or any other that had a whole or if the white areas were only connected by a white line of 1 pixel width then the program would return the wrong value. I traced this limitation back to the function follow(). This function can only follow a single contour that did not intersect itself, faced with any other case and it would not work for our purposes.

Overall I still consider this as a huge success and I have to give myself a pat on the back for such good work. For this I'm giving myself a grade of 10!!! woooohhooooo!!!!

I would first like to thank Dr. Soriano for all the information and insights she shared. I would also like to thank Thirdy for helping me in discovering the limitations of the follow command.

square
triangle

circle

tshirt 1

tshirt 2

irregular

irregular 2

irregular 3

Im = imread('circle.jpg');
Im = im2bw(Im, 0.5);

[x,y] = follow(Im);
len=length(x);
x(len+1)=x(1);
y(len+1)=y(1);

//theoretical area via pixel counting
A1 = sum(Im);
//the additive term 'round((len+1)/2)' accounts for the sections of the boundary that is not included in the sumation
A = round(sum((x(1:len).*y(2:len+1)-y(1:len).*x(2:len+1)))/2) + round((len+1)/2);

//area checking
A==A1

Activity 2 and 3 pictures

2009-06-22T19:04:00.000-07:00

Activity 1

2009-06-17T20:31:00.001-07:00

Once we have the required data or information, making graphs, charts, or plots is never a problem for us. With the aide of our modern computers and the numerous calculating and graphing softwares such tasks are usually done efficiently and hassle free. But before the age of even the earliest computers fairly accurate graphs and plots already existed; which is, when you think of it, an amazing feat. I don't know how this plots were made back then, but by simply using the powers of wikipedia we can somewhat get an idea of how scientists back then made use of different types of graphing papers (http://en.wikipedia.org/wiki/Graphing_paper) and odd looking implements (http://en.wikipedia.org/wiki/French_curve) to generate their graphs. In this first activity we recreated a very old plot by using MS Paint and MS Excel.

We took very old handwritten plots from books and journals as far back as the 1920's ,or even older, and scanned them. Now using the select tool of MS paint we took note of the pixel coordinates of all the significant points of the plot. Specifically we took the pixel values of the axis of the plot along with their corresponding physical values indicated on the plot. From these data we obtained the conversion from the pixel values to their physical values. For each tick mark on both axes, we took the difference in their pixel values and also their physical values. From the ratio of these two values we would already have a conversion to use on the curve:

But since the plots were hand drawn these difference in pixel values were not constant, still the variations were not that large so I simply took the average. The more significant concern is that this particular doesn't have an origin with pixel coordinates (0,0) and the plot itself did not even start at 0,0 of the physical values. Also in the image file the direction of increasing pixel value on the vertical was opposite of that with the vertical axis of the plot. Still these problems can be easily addressed by using some basic algebra, point-slope or two-point form of a line, along with some concentration and patience. Doing this I arrived at a modified equation for conversion of the form:

Using the pixel coordinates of the points along the curve in this equation I tried to recreate this plot in MS Excel and compared it with the original. In comparing the two plots we simply use the original scanned plot as the background of the plot in Excel. But upon comparing the recreated curve with the original it was obvious that the two did not match. Then I found out that the x-axis of the plot I chose was not the same through out the plot.

It appears that the origin of the second segment coincided with the 1 value of the first. Also, I found out that the ratio of the pixel values and the physical values of the axes was not the same. I tried to add value of 1 to offset the second segment of the curve and I also used a different conversion for the two parts to improve the appearance of the recreated curve. But as you can see from the figure bellow, though the general shape and trend of the plot was replicated the corresponding points and tick marks did not coincide with the original plot.

After doing a lot of tweaking with the math and even more lucky trial-and-errors I found out that it is not possible to fully reconstruct this plot using this method. So with the permission of my professors I decided to break the curve into two plots and then do the reconstruction.

Finally the recreation is a success...!!! The recreated curve fits very well with the original and even the points and the tick marks are spot on. I would like to give myself a 10 for this but since I was late I guess a 9 would do... hahahahahahhah!!!!!

And since this was a success I have to thank all those who helped me... first my professors Dr. Soriano and Dr. Perez, Thirdy for helping me find this plot, Jica for sending me my files, and Winsome for telling me that I need to sign in before I can do anything with my blog.

Source of Plot: A. Wilmen Duff, et al, "Physics for Students of Science and Engineers", Maple Press Co., York PA, 1926