Hough Transform

Hough transform is used to detect lines in an image. It uses the parametric representation of a line:


where rho is the distance of the line from the origin and theta is the angle between the x-axis and the vector perpendicular to the line.

Fortunately, there are built-in functions in MATLAB that are useful for this activity.

Let us start with this image, taken in Cibo. The design of the place uses a lot of lines.




First, convert the image above to grayscale and use the edge function to detect the edges in the image. I tried the Canny algorithm.






So now we have the edges in the image. Let us now implement the hough function. The hough function returns this plot:



The x- and y- axis are theta and rho, respectively. The bright spots correspond to the most prominent lines in the image.

To detect the brightest spots or the peaks, let us use the function houghpeaks. The function returns a list of theta and rho coordinates that correspond to the peaks in the hough transform plot.




Shown in the figure above are the three peaks. To plot the lines that correspond to these peaks, we use the function houghlines. The figure below shows the superposition of the plotted lines on the grayscale of the original image. We can see that the lines coincide with the lines in the ceiling and the table. The lines from the table are somewhat unexpected because visually, the ceiling lines are more prominent. I expected that before the lines from the table will be detected, all the ceiling lines will be detected first.



I also noticed that the detected lines depend on the edge detection algorithm. Shown below are the different results for the different edge detection methods.



I give myself 10 points for this activity. The hough transform was used to identify and draw the lines detected in the image. I also investigated the effect of using different edge detection algorithms, although I cannot explain the differences yet.

Stereometry

Stereometry

We can see in three dimensions because we have two eyes, with one eye seeing a slightly different scene than the other. The two images of the scene are taken with a horizontal offset that is the distance between our eyes. In this activity, we will simulate how the eyes see three dimension with two images, or stereo vision.

There are two ways to do stereo imaging. We can use two identical cameras separated by a distance b from each other. Alternatively, we can use just one camera and move it by a distance b. For practicality, we went with the second option.

Two images of the object are taken. The second image is taken by moving the camera horizontally by a distance b apart from the location where the first image is taken. The vertical position of the camera was maintained at a constant height. The same camera settings were used for both images.

Figure 1 shows the geometry of the setup.

Figure 1. Geometry of the setup

Using similar triangles, we can see that:
and that:


Solving for z gives us:


If we do this for many points, we can reconstruct the 3D surface of the object.

The calibration of the camera was done in a previous activity, and the focus was determined to be f = 1063.6 pixels. The distance b is b = 10cm.

The 3D object used here is a box, covered with graphing paper such that it will be easy to correspond the points on the surface. Figure 2 shows the two images captured with the camera displaced a distance b.

Figure 2. The two images taken with the camera displaced a distance b.

Several (x,y) coordinates of corresponding vertices on the two images were selected. Only the x-coordinate is important because the y values are kept constant.

The z values are calculated for each point following the equations above.

Figure 3 shows the reconstruction of the surface.

Figure 3. Reconstruction of the 3D surface of the cube.

The 3D surface of the box is rendered, shown in Figure 3. The reconstruction was successful, the 3D shape of the box can be recognized. I give myself 9 points for this activity.

Thanks to Kaye and Miguel Sison, my groupmates in a previous class where this activity was previously done.

3D Reconstruction using Gray Code Illumination

Introduction

Gray code illumination (GCI) is one of the many techniques in rendering a 3D reconstruction of an object. It is a non-invasive technique, where only images of object are taken, under different pattern projection conditions. Another advantage is that it is insensitive to ambient illumination.
GCI works on the idea that patterns illuminated against a flat background will be distorted by objects placed against the background. These distortions are the cues to the shape of the object.

Gray code patterns

The patterns projected are a series of black and white stripes in the Gray code pattern. Each pattern has a bit weight, which is multiplied to the pattern itself. When these weighted patterns are summed, we now have what is called the bit plane stack (BPS). Each column in the BPS has a unique number. The distortion of this BPS due to the object is what we will keep track. Without the object, we have a background bit plane stack (BBPS). With the object, we have the OBPS.

Figure 1. Making the bit plane stack by summing the gray code patterns multiplied by the corresponding bit weights.

Methodology

Gray code patterns (Figure 2) are projected against a white background. Their images are taken, which will serve as the references for the distortion. The 3D object to be reconstructed is placed against the same white background. Several images of black and white patterns in the Gray Code configuration is projected onto the 3D object. For each gray code pattern, two images are taken. One is for when the object is present, and one is for background, in the absence of the object.

Figure 2. Object and background illuminated by different gray code patterns.

The images are then binarized. Each binarized image (from the background and object) is multipled by its corresponding bit weight, and added to make the background bit plane stack and object bit plane stack (Figure 3).

Figure 3. Left: Object bit plane stack. Right: Background bit plane stack

In the bit plane stacks, each stripe or color is a unique number. The shape of the object manifests in the distortion of the object bit plane stack. For each color, the displacement of the stripe is taken by subtracting the right-side edges. The differences are now the height of the object. The resulting 3D surface is plotted below (before and after median filtering).

Results

Reconstruction of the 3D surface 3 levels, which correspond to the levels of the pyramid. However, noise is still present, seen as sharp peaks in the height reconstruction (Figure 4a). This is remedied by applying median filtering. The resulting surface reconstruction is now smoother (Figure 4b).

The heights of each level are not the same, with the lowest level being the thickest.

(a) (b)

Figure 4. Left: Reconstruction before median filtering. Right: Reconstruction after median filtering.

To summarize, different patterns in Gray code configurations were used to successfully reconstruct a 3D surface. Several noise filtering methods were applied in the intermediate steps of the whole process. I give myself 10 points for this activity.

Many many thanks to Kaye for guidance in implementing the algorithms for this activity.

Camera Calibration using MATLAB toolbox

Camera Calibration

The camera calibration toolbox of MATLAB downloaded from http://www.vision.caltech.edu/bouguetj/calib_doc/ is an easier was of calibrating the camera.

It has an interface, which makes things easier.

But first several images of the grid oriented in different directions were taken. Shown below are the 16 images of the grid. In each image, the camera is fixed and the grid is moved.

Figure 1. Grid images.

The images in the folder must be named in sequence, for example DSC1.jpg, DSC2.jpg, DSC3.jpg, and so on. To input the images, click the button [Image names].

Figure 2. Calibration toolbox interface

The images are read by the program. The corners of the grid are selected one by one. A window showing the grid is displayed. The inner corners of the grid are selected, because the outer corners might be damaged.

Figure 3. Corner selection window.

The corners are selected starting from the upper left, clockwise. The GUI will ask for the dimension of the boxes in the grid. For this grid, each box is 1 inch wide and 1 inch long.

The corners are then extracted by the toolbox (use [Extract grid corners] button), and are labeled in the image (Figure 4).

Figure 4. Corners are extracted.

This is done for the rest of the images. You must be careful not to make any mistake. Otherwise, you will have to do this all over again, starting from the image input process. (This is most frustrating when you are in the last few images)

After the corners are extracted for all the images (in this case 16), the camera can now be calibrated. Click the button [Calibration]. The intrinsic parameters of the camera is now calculated by the program. Shown in Figure 5 is the output of the program.

Figure 5. Intrinsic camera parameters.

We have now calculated for the intrinsic parameters of the camera, such as focal length, principal point, skew, etc.

The extrinsic parameters can also be measured (click [Show extrinsic] button). The camera centered view or fixed-camera view is displayed (Figure 6). Here, the different orientations of the grid as captured by the camera is projected.

Figure 6. Extrinsic parameter (camera-centered)

Figure 7 shows the extrinsic parameters (world-centered). Here, the different locations of the camera is projected, as seen by the grid. This is basically the same as in Figure 6, because location is relative.

Figure 7. Extrinsic parameter (world-centered).

The camera was calibrated using the Camera Calibration Toolbox in Matlab. For more detailed instructions, refer to this link: http://www.vision.caltech.edu/bouguetj/calib_doc/htmls/example.html

I give myself 10 points for this activity. I was able to get the intrinsic parameters the camera using the toolbox.

Curve Spread

A paper by Dr. Soriano, Ms. Araullo, and Dr. Saloma on front-view biometrics published in Pattern Recognition Letters uses features from the silhouette of the object. In this activity, we will learn the technique used to extract the curve spread of an object.

I used an image of a key because it has different contour types. It has straight, concave, and convex edges. This image is not yet an application of this method, just an aid to help me program the algorithm.

First, the image of the key is thresholded so that the region of interest, which is the key, is well separated from the background. Morphological cleaning operations were performed to remove floating pixels and close holes in the object. Now that we have a cleaned image, the edge pixels are taken. The command edge in Matlab returns an array the size of the original image, with the edge pixels highlighted.


Now we have a one-pixel thick outline of the original image. What we want next is to get the coordinates of these pixels. Initially, I used Matlab to program this activity but I encountered a problem. I tried to use regionprops then PixelList to get the coordinates of the edge pixels but the coordinates it returns are arrange per column. What we need here is a list of pixels arranged as it was followed. The command follow in Scilab does exactly this.

Now that we have a list of pixels, we convert it to Freeman chain code by starting at an arbitrary pixel then following the line in a certain direction. For the image I used, I followed the contour in a clockwise direction. In making a Freeman chain code to store the relative directions of one pixel with respect to a previous pixel, we assign a number based on the following:



For example, we start at the certain pixel, and following the contour clockwise the next pixel happens to be in the upper left corner of the current pixel, the next pixel will be labeled '3'. If after that pixel labeled '3', the next pixel is directly below, this next pixel will be labeled '6', and so on.

Now we have the Freeman code chain for the contour of our image. To get its curvature information, we take the first order discrete difference of this array, or the difference of the current pixel and the previous pixel. Then for each pixel, we sum the values of the previous, current, and next pixel. This would produce an array of zeros, positive and negative numbers, which now gives us clues to the curvature of the portions it represents.

RESULTS
Let us look at the result of the curvature information for a portion of the key:


Straight segments are represented by consecutive zeros, concave regions are represented by positive numbers, and convex regions are represented by negative numbers.

One must be careful in using the follow command in Scilab. The outline of the object must be closed, otherwise the inner and outer edge of the outline will be followed, doubling the length of the coordinate list.

I give myself 10 points for successfully using the Freeman vector code in identifying the curvatures in the image.

Adaptive Skin Color Modeling (based on Dr. Soriano's paper)

This activity is on adaptive skin color tracking. The objective is to be able to track the skin under changing illimination conditions. This means that the skin must be tracked even if the subject in the video changing location, walking from a fluorescent lit room to a sunny outdoors.

SKIN LOCUS

First, all possible colors of skin under different illumination conditions must be determined. Its values in the normalized chromaticity coordinates is called the skin locus. Ideally, this is calculated by taking images with all possible permutation of illumination condition and white balance setting of the camera and adding up all the histograms of the skins.

In our case, we were not able to take images under all possible illumination settings that the camera can white balance. Instead, we took the histograms of the skins under different illumination as captured in our video. Shown below are the cropped portions of the face where the histogram was taken.

The combined histograms of the cropped portions of skin is known as the skin locus, shown below.

ADAPTIVE HISTOGRAM BACKPROJECTION

First, a portion of the face or skin is cropped. This will give the initial histogram of which will be used to locate regions in the image corresponding to the skin. Filter its histogram with the skin locus histogram to remove values which do not correspond to skin.

This filtered histogram is used to segment the current image. The segmented image is then binarized and cleaned using morphological operations. The box enclosing the largest cluster is used to select a portion in the current image. The histogram of this portion is used to segment the next image or frame. But in the next frame, the search region or the region where the histogram backprojection will be applied is limited only to a box slightly larger than the box enclosing the largest cluster. This is done under the assumption that the object or face has not moved too swiftly, moving only a few pixels between frames. This is applied to all the images in the video.


RESULTS

Presented below are slices of the video where the subject is under different illumination conditions. The success of the technique is shown by the correct tracking even if the light illuminating the subject changes.


SOME COMMENTS

Sometimes the largest cluster unusually jumps a large distance, probably because the clusters so disconnected or maybe the wall and the subject becomes similar in color. To address this solution (as suggested by Ma'am Jing), the distance between the centroid of the current largest cluster and the centroid of the previous largest cluster is checked. If it moved more than a set threshold distance, the current largest cluster is rejected. The previous largest cluster is used in the process.

In this activity, I worked with Kaye and Thirdy. Thanks to Kirby for sharing his face for the video hehe :D

I give myself 10 points for this activity. Face/skin tracking was successfully done even with changing illumination. Also, I had to code some functions myself because they were not available in Scilab. :)

Reference:
Soriano, Martinkauppi, Huovinen, Laaksonen, Adaptive skin color modeling modeling using the skin locus for selecting training pixels, Pattern Recognition, Vol. 36, No. 3, March 2003.

Activity 19 - Restoration of blurred image

In this activity, the objective is to demonstrate restoration of an image corrupted with a known degradation function (in this case, motion blur) and additive noise. [1]

A grayscale image containing text was copied from the web and corrupted by applying motion blur and noise, following the equation:
where f(x,y) is the original image, h(x,y) is the blurring (degradation) function, n(x,y) is the added noise. Or in Fourier space, the corrupted image is given by:

where G,H,F, and N are the Fourier transforms of the corresponding terms in the previous equation [1].

The degradation function is modeled by the equation:

where a and b are the total distance for which the image has been displaced in the x- and y-direction and T is the total time of displacement [1].

The image is restored by Weiner filtering or minimum mean square error filtering or least square error filtering by applying the following equation, which gives the frequency-domain estimate of the restored image:
where:
H(u,v) is the degradation or blurring function
|H(u,v)|^2 = H(u,v)* H(u,v)
Sn(u,v) = |N(u,v)|^2 (power spectrum of the noise)
Sf(u,v)| = |F(u,v)|^2 (power spectrum of original image)

Most of the time, the power spectrum of the original image is unknown, in which case we apply the other form of the Weiner filter given by:
where K is a specified constant.

The original image (left of the divider) is blurred by applying the degradation function (where a=b=0.1 and T was varied) and Gaussian noise was added. The restored images are shown below (right side of divider).

Original image | T=1, T=2, T=3.

All the images above are restored assuming the power spectrum of the undegraded image is known. We can see that as long as the power spectrum of the original image is known, the image can be restored quite well no matter how the image was degraded.

Now we look at the effects of varying a and b to the restoration. Shown below are the restored images blurred with varying a and b.

Original image | a=b=0.01, 0.03, 0.05, 0.07, 0.09

We can see that varying a and b has minimal effect on the restoration, as long as the power spectrum of the undegraded image is known.

What if we do not know the power spectrum of the original image? We try to predict the value for K. Shown below are the restored images using different values of K.

Original image | (top row) K = 1, 0.1, 0.01 (bottom row) K = 0.001, 0.0001, 0.00001

Restoring the corrupted image without knowledge of the power spectrum of the original image proved to be more difficult. Different K values were tried out, and it was observed that starting from K=0.001 below, the image is restored quite well. But for values greater than K=0.001, the restored image is still blurred.

Problems encountered:
I, as well as other classmates, had trouble displaying the blurred image. The resulting image which was supposedly corrupted with motion blur did not display the expected blurred image with streaks but instead showed a superposition of the "initial and final frames" of the supposedly blurred image. However, the information contained in the array in the program seems correct, because reconstruction was done successfully. The problem is just in displaying the blurred image.

I give myself 9 points for this activity. I would like to thank Thirdy Buno for useful discussions.



References:
[1] A19 - Restoration of blurred image - Dr. Maricor Soriano

Activity 18 - Noise Models and Basic Image Restoration

In this activity, we investigate how different noise reduction filters restore images with different kinds of added noise:

Noise models: Gaussian, Rayleigh, Gamma, Exponential, Uniform, and Impulse (or salt-and-pepper) Noise
Noise filters: Arithmetic Mean, Geometric Mean, Harmonic Mean, Contra-harmonic Mean

We used each of the above noise models in conjunction with each of the noise filters. Below are the results:


The arithmetic mean filter gave the cleanest restoration for all the noise models. The harmonic mean filter comes close, except for the salt and pepper noise. The harmonic mean is also better than arithmetic mean for the Gaussian noise. The geometric mean filter was the least effective filter. Except for the gamma noise, the geom mean filter converted the gray pixels to black.

I give myself 10 points for this activity. I was able to perform noise addition for different noise models, and restore the images using different noise filters.

Thanks to Kaye and Miguel for help in some of the noise models.

Activity 17 - Photometric Stereo

Photometric stereo is a method that gets information from shadow to reconstruct a shape. Images of an object under different illumination sources is taken, and from these images which contain shadow information, the shape can be retrieved. We were given the file Photos.mat which contains four images of a hemisphere (shown below), with illumination coming from different directions. The x,y,z components of each of the four light sources stored as a vector are given below:

V1 = [0.085832, 0.17365, 0.98106];
V2 = [0.085832, -0.17365, 0.98106];
V3 = [0.17365, 0, 0.98481];
V4 = [0.16318, -0.34202, 0.92542];

Above are the images of a hemisphere illuminated from different directions. The directions stored in an array shown below, where V1 and V2 are the directions for the top row image, respectively, and V3 and V4 are the directions for the bottom row, respectively. The light source matrix is defined as

where each row is a light source and each column is x,y,z component of the source [1]. The image is represented in matrix form as I=Vg. I here is the matrix of the values of each pixel in each image, given by:

Since we know I and V, we can solve for g using:

The normal vector is given by:

The surface elevation f(u,v) at a point (u,v) is given by:

where

The line integral above was implemented by first computing all partial derivatives for each pixel such that we have two 2D matrices (with the same size as the images) containing the partial derivatives. The first term is computed by summing the uth row from columns 1 to v, and the second term is computed by summing the vth column from rows 1 to u. The code is shown below:

for u=1:x
for v=1:y
t1 = dfdx(u,1:v);
t2 = dfdy(1:u,v);
f(u,v) = sum(t1)+sum(t2);
end
end

Plotting f(u,v) in 3D gives us:

I give myself 10 points for this activity. The hemisphere is successfully reconstructed. However, I cannot explain the visible grooves in the reconstruction, which seems to form a cross in the middle.

Thanks to Martin for clarifying the integral limits.