The correct segmentation of the iris and the pupil may fail for different reasons (See Fig 5). Concerning the pupil, the algorithm may fail in the detection of its center if in the image are present dark areas, like in case of non uniform illumination and if the subject is wearing glasses (a-b). One other source of problems is if the reflexes are in the pupil area and are not properly removed, because they can be detected as pupil edges, leading to a smaller pupil (c). Moreover the pupil edge can be detected erroneously if a part of its edge is occluded by the eyelid or by the eyelashes (d-f). Concerning the iris, since its fitting is constrained by the pupil shape, if the pupil detection is wrong consequently the iris can be completely wrong or deformed. Even if the pupil is detected correctly, when the edge between the iris and the sclera come to have low contrast, for reasons like not uniform illumination, not correct camera focus, or bright color of the eye, the algorithm may fail in finding the limbus (g-h).

Over the whole set of images, the correct segmentation rate is 94, 5%, attesting a good efficacy of the algorithm. In fact it is able to segment properly the iris area (See Fig 4) with changhing size of the puil, in presence of glasses and heavy reflexes (a-c), bushy eyelashes (d-e), iris and pupil partially occluded (f).

S5001−R−06 S5017−L−03 S5048−R−06 S5044−R−05

Fig. 4. Examples of correct segmentation on the CASIA Iris Database, on the Thousand subset. The magenta ellipse defines the pupil contour, while the green one is the limbus. The red dots represent the points used to compute the limbus ellipse equation, while the white ones are those removed for their possible belonging to the eyelashes or to wrong estimation of the edge.

**5.4.1 Hardware implementation**

stream of images.

of the image elaboration.

**5.4.2 Image acquisition and segmentation**

(see for example Fig. 6, center-right).

**6. Eye-tracking**

**6.1 Tracking**

produces that projection.

properly the pupil and the iris in the 100% of the cases.

The head-mounted device is endowed with two cheap USB web cams (Hercules Deluxe Optical Glass) that provide images at a resolution of 800 × 600 pixels, with a frame rate of 30 fps. The cameras were mounted in the inner part of a chin strap, at a distance of ≈ 60 mm from the respective eye. At this distance, the field of view provided by the cameras, is [36◦, 26◦], that is more than enough to have a complete view of the eyes. To make them work in infra-red light, the IR-cut filter were removed from the optics, and substituted with a IR-pass filter, with cut frequency of 850 nm. To have a constant illumination of the images, both in daylight or in indoor environments and during night time, the system were endowed with three IR illuminators, that help to keep constant the contrast and the illumination of the

The Perspective Geometry of the Eye: Toward Image-Based Eye-Tracking 23

In fact, the illuminators produce visible corneal reflexes, that are used as reference feature in other kinds of eye trackers (Eizenman et al., 1984; Morimoto et al., 2000). In our case, since we are seeking to use the limbus position to track the eye, if the reflex, depending on the position of the eye falls in its correspondence, it can lead to the detection of a wrong edge, thus to a wrong gaze estimation. To prevent this case, and considering that the points affected by the reflexes are few respect to the entire limbus edge, these points are removed at the beginning

The developed algorithm was tested on three sets of images, taken from different subjects. In each set, the subjects were asked to fixate a grid of points, in order to have the gaze ranging from −30◦ and −30◦ of azimuth, and from −20◦ and 20◦ of elevation, with a step of 5◦. In this way each set is composed by 117 images where the gaze direction is known. The azimuth and elevation angles were defined following a Helmholtz reference frame (Haslwanter, 1995). The use of a transformation of the image from a Cartesian to an elliptic domain allows the algorithm to work properly on the segmentation of the pupil and consequently, as explained in Sec 3, on the segmentation of the iris, even in the cases where the iris is drastically occluded

Considering that the images are captured in an optimal condition, *i.e.* in an indoor environment where the only sources of IR light are the illuminators and the subjects do not wear glasses, and with the eye correctly centered in the image, the algorithm is able to segment

Once the iris circle is detected steadily on the image plane, and its edge is fitted with an ellipse, knowing the coefficient matrix **Z** of the quadratic form, it remains to estimate the gaze direction. This can be obtained computing which is the orientation in space of the circle that

**Σ** = diag(*λ*1, *λ*2, *λ*3)

Because **Z** is a symmetric matrix, it can be diagonalized, leading to a matrix:

Fig. 5. Examples of wrong segmentation on the CASIA Iris Database, on the Thousand subset. The magenta ellipse defines the pupil contour, while the green one is the limbus. The red dots represent the points used to compute the limbus ellipse equation, while the white ones are those removed for their possible belonging to the eyelashes or to wrong estimation of the edge.

#### **5.4 The proposed system**

The images available in the CASIA-Iris-Thousand are characterized by a gaze direction that is close to the primary, and are taken by a camera positioned directly in front of the eye. The images provided by such a configuration are characterized by a pupil and an iris whose edge is close to a perfect circle. In fact, since the normal to the plane that contains the limbus is parallel to the optical axis of the camera, the projected ellipse has eccentricity close to zero. While the feature to be searched in the image has a circular shape, the technique of re-sampling the image with a transformation from Cartesian to polar coordinates is an effective technique (Ferreira et al., 2009; Lin et al., 2005; Rahib & Koray, 2009). In such domain, starting from the assumption that its origin is in the center of the iris, the feature to be searched, *i.e.* the circle of limbus, is transformed to a straight line, and thus it is easier to be individuated than in the Cartesian domain.

On the other side, considering not only the purpose of biometrics but also the eye tracking, the eye can be rotated by large angles respect to the primary position. Moreover, in our system, the camera is positioned some centimeters lower than the eye center in order to prevent as much as possible occlusions in the gaze direction. The images obtained by such a configuration are characterized by the pupil and the iris with an eccentricity higher than zero, that increases the more the gaze direction differs from the optical axis of the camera (see for example Fig. 6, top-right).

Since the presented eye-tracking system is based on the perspective geometry of the pupil and iris circle on the image plane, it is important for the system that the relative position between the eye and the camera stay fixed. On the oder side, for a good easiness of use and freedom of movement of the subject it is important the the system allows free head movement. For this purpose, we developed an head mounted device, in order to guarantee both these features.

### **5.4.1 Hardware implementation**

22 Will-be-set-by-IN-TECH S5018−R−02 S5111−R−07 S5052−R−09 S5086−L−08

> (a) (b) (c) (d) S5082−L−03 S5065−L−01 S5012−R−03 S5012−R−05

> (e) (f) (g) (h)

The images available in the CASIA-Iris-Thousand are characterized by a gaze direction that is close to the primary, and are taken by a camera positioned directly in front of the eye. The images provided by such a configuration are characterized by a pupil and an iris whose edge is close to a perfect circle. In fact, since the normal to the plane that contains the limbus is parallel to the optical axis of the camera, the projected ellipse has eccentricity close to zero. While the feature to be searched in the image has a circular shape, the technique of re-sampling the image with a transformation from Cartesian to polar coordinates is an effective technique (Ferreira et al., 2009; Lin et al., 2005; Rahib & Koray, 2009). In such domain, starting from the assumption that its origin is in the center of the iris, the feature to be searched, *i.e.* the circle of limbus, is transformed to a straight line, and thus it is easier to be individuated than in the

On the other side, considering not only the purpose of biometrics but also the eye tracking, the eye can be rotated by large angles respect to the primary position. Moreover, in our system, the camera is positioned some centimeters lower than the eye center in order to prevent as much as possible occlusions in the gaze direction. The images obtained by such a configuration are characterized by the pupil and the iris with an eccentricity higher than zero, that increases the more the gaze direction differs from the optical axis of the camera (see for example Fig. 6,

Since the presented eye-tracking system is based on the perspective geometry of the pupil and iris circle on the image plane, it is important for the system that the relative position between the eye and the camera stay fixed. On the oder side, for a good easiness of use and freedom of movement of the subject it is important the the system allows free head movement. For this purpose, we developed an head mounted device, in order to guarantee both these features.

Fig. 5. Examples of wrong segmentation on the CASIA Iris Database, on the Thousand subset. The magenta ellipse defines the pupil contour, while the green one is the limbus. The red dots represent the points used to compute the limbus ellipse equation, while the white ones are those removed for their possible belonging to the eyelashes or to wrong estimation

of the edge.

**5.4 The proposed system**

Cartesian domain.

top-right).

The head-mounted device is endowed with two cheap USB web cams (Hercules Deluxe Optical Glass) that provide images at a resolution of 800 × 600 pixels, with a frame rate of 30 fps. The cameras were mounted in the inner part of a chin strap, at a distance of ≈ 60 mm from the respective eye. At this distance, the field of view provided by the cameras, is [36◦, 26◦], that is more than enough to have a complete view of the eyes. To make them work in infra-red light, the IR-cut filter were removed from the optics, and substituted with a IR-pass filter, with cut frequency of 850 nm. To have a constant illumination of the images, both in daylight or in indoor environments and during night time, the system were endowed with three IR illuminators, that help to keep constant the contrast and the illumination of the stream of images.

In fact, the illuminators produce visible corneal reflexes, that are used as reference feature in other kinds of eye trackers (Eizenman et al., 1984; Morimoto et al., 2000). In our case, since we are seeking to use the limbus position to track the eye, if the reflex, depending on the position of the eye falls in its correspondence, it can lead to the detection of a wrong edge, thus to a wrong gaze estimation. To prevent this case, and considering that the points affected by the reflexes are few respect to the entire limbus edge, these points are removed at the beginning of the image elaboration.

#### **5.4.2 Image acquisition and segmentation**

The developed algorithm was tested on three sets of images, taken from different subjects. In each set, the subjects were asked to fixate a grid of points, in order to have the gaze ranging from −30◦ and −30◦ of azimuth, and from −20◦ and 20◦ of elevation, with a step of 5◦. In this way each set is composed by 117 images where the gaze direction is known. The azimuth and elevation angles were defined following a Helmholtz reference frame (Haslwanter, 1995). The use of a transformation of the image from a Cartesian to an elliptic domain allows the algorithm to work properly on the segmentation of the pupil and consequently, as explained in Sec 3, on the segmentation of the iris, even in the cases where the iris is drastically occluded (see for example Fig. 6, center-right).

Considering that the images are captured in an optimal condition, *i.e.* in an indoor environment where the only sources of IR light are the illuminators and the subjects do not wear glasses, and with the eye correctly centered in the image, the algorithm is able to segment properly the pupil and the iris in the 100% of the cases.

### **6. Eye-tracking**

Once the iris circle is detected steadily on the image plane, and its edge is fitted with an ellipse, knowing the coefficient matrix **Z** of the quadratic form, it remains to estimate the gaze direction. This can be obtained computing which is the orientation in space of the circle that produces that projection.

### **6.1 Tracking**

Because **Z** is a symmetric matrix, it can be diagonalized, leading to a matrix:

$$\boldsymbol{\Sigma} = \text{diag}(\lambda\_1, \lambda\_2, \lambda\_3)$$

and consider **x**� = **R***T***x**, the equation of the projective cone in the new basis is:

at the intersection of the cone with the plane *w* = 1

reference frame.

**6.2 Calibration**

(4) the radius of the eye.

**6.3 Estimation of the gaze direction**

which is a cone expressed in canonical form, whose axis is parallel to **e3**. Now, look for a while

The Perspective Geometry of the Eye: Toward Image-Based Eye-Tracking 25

*λ*1*x*�2 + *λ*2*y*�2 = 1 whose axes are parallel to **e**1 and **e**2, and semiaxis length are √*λ*1 and √*λ*2. If we consider to cut the cone in Eq.12 with a plane tilted along **e**1 there exist a particular angle *θ* which makes tha plane to intersect the cone in a circle: this circle will be the limbus and *θ* its tilt angle in the basis described by the rotation matrix **R**. As suggested in Forsyth et al. (1991), to find *θ* it is possible to exploit the properties of circle to have equal semiaxes or, equivalently, to have equal coefficient for the *x*�2 and *y*�2 terms in the quadratic form. Equality of the *x*�2 and *y*�2 coefficients is achieved by a rotation along the *x*� axis by an angle *θ* = ± arctan ��*λ*2−*λ*1

which set both the coefficients equal to *λ*1. The normal to the plane that intersects the cone in

10 0

and **R***cam* is the matrix describing the rotation of the camera with respect the fixed world

Since it is not possible to measure with the desired precision the relative position between the center of the eye and nodal point of the camera, and particularly the relative orientation, we developed a calibration procedure with the purpose of estimating the implicit parameters. The subject is asked to fixate a calibration grid of 25 points in a known position. The grid of points is designed to make the subject fixate with an azimuth angle between −30◦ and 30◦ with steps of 15◦, and with an elevation angle between −20◦ and 20◦ with steps of 10◦.

The calibration procedure is based on a functional whose minimization provides: (1) the eye position respect to the camera, (2) the camera orientation respect to a fixed reference frame,

In order to have a validation of the algorithm, the estimation of the fixation angle is computed over a different set of points respect to the calibration grid. The grid of point used for the test is designed to make the subject fixate with an azimuth angle between −20◦ and 20◦ with steps

of 5◦, and with an elevation angle between −10◦ and 10◦ with steps of 5◦.

0 cos *θ* − sin *θ* 0 sin *θ* cos *θ* ⎤ ⎥ ⎥ ⎦

a circle, expressed in the camera coordinate system, is **n** = **R***cam***RR***θ* [0; 0; −1], where:

⎡ ⎢ ⎢ ⎣

**R***θ* =

*λ*2 3

*λ*1*x*�2 + *λ*2*y*�2 + *λ*3*w*�2 = 0 (12)

. This is the ellipse:

*λ*1−*λ*3 � ,

Fig. 6. Subset of the images taken from a subject with the IR camera. The subject is fixating to a grid of nine points at the widest angles of the whole set. The magenta ellipse defines the pupil contour, while the green one is the limbus. The red dots represent the points used to compute the limbus ellipse equation, while the white ones are those removed for their possible belonging to the eyelashes or to wrong estimation of the edge. The blue lines represent the parabolas used to remove the possible eyelashes points.

with *λ*1, *λ*2, *λ*3 ∈ **R**. This transformation from **Z** to **Σ** is just a change of basis, and thus **Σ** may be expressed as **Σ** = **R**−**1ZR**, where **R** is the matrix changing between the actual orthonormal basis to a new one, formed by three eigenvectors of **Z**. The columns of **R** are the components of the eigenvectors **e1**, **e2**, **e3**, and the elements *λ*1, *λ*2, *λ*3 of the diagonal of **Σ** are the associated eigenvalues. For the Sylvester's Law of inertia the signature of **Σ** is equal to that of **Z**, (−, +, +); thus only one eigenvalue is negative, and the other are positive. We assume *λ*3 < 0 and *λ*2 > *λ*1 > 0. If we apply the transformation matrix **R** to **QZ**(**x**), we obtain:

$$\begin{aligned} \mathbf{Q} \mathbf{z} (\mathbf{R} \mathbf{x}) &= (\mathbf{R} \mathbf{x})^T \mathbf{Z} (\mathbf{R} \mathbf{x}) \\ &= \mathbf{x}^T \mathbf{R}^T \mathbf{Z} \mathbf{R} \mathbf{x} \\ &= \mathbf{x}^T \mathbf{E} \mathbf{x} \end{aligned}$$

and consider **x**� = **R***T***x**, the equation of the projective cone in the new basis is:

$$
\lambda\_1 \mathbf{x}'^2 + \lambda\_2 \mathbf{y}'^2 + \lambda\_3 \mathbf{w}'^2 = \mathbf{0} \tag{12}
$$

which is a cone expressed in canonical form, whose axis is parallel to **e3**. Now, look for a while at the intersection of the cone with the plane *w* = 1 *λ*2 . This is the ellipse:

3

$$
\lambda\_1 \mathbf{x}'^2 + \lambda\_2 \mathbf{y}'^2 = 1
$$

whose axes are parallel to **e**1 and **e**2, and semiaxis length are √*λ*1 and √*λ*2. If we consider to cut the cone in Eq.12 with a plane tilted along **e**1 there exist a particular angle *θ* which makes tha plane to intersect the cone in a circle: this circle will be the limbus and *θ* its tilt angle in the basis described by the rotation matrix **R**. As suggested in Forsyth et al. (1991), to find *θ* it is possible to exploit the properties of circle to have equal semiaxes or, equivalently, to have equal coefficient for the *x*�2 and *y*�2 terms in the quadratic form. Equality of the *x*�2 and *y*�2 coefficients is achieved by a rotation along the *x*� axis by an angle *θ* = ± arctan ��*λ*2−*λ*1 *λ*1−*λ*3 � , which set both the coefficients equal to *λ*1. The normal to the plane that intersects the cone in a circle, expressed in the camera coordinate system, is **n** = **R***cam***RR***θ* [0; 0; −1], where:

$$\mathbf{R}\_{\theta} = \begin{bmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & -\sin \theta \\ 0 & \sin \theta & \cos \theta \end{bmatrix}$$

and **R***cam* is the matrix describing the rotation of the camera with respect the fixed world reference frame.

#### **6.2 Calibration**

24 Will-be-set-by-IN-TECH

*H* = −30◦ *V* = 20◦ *H* = 0◦ *V* = 20◦ *H* = 30◦ *V* = 20◦

*H* = −30◦ *V* = 0◦ *H* = 0◦ *V* = 0◦ *H* = 30◦ *V* = 0◦

*H* = −30◦ *V* = −20◦ *H* = 0◦ *V* = −20◦ *H* = 30◦ *V* = −20◦

Fig. 6. Subset of the images taken from a subject with the IR camera. The subject is fixating to a grid of nine points at the widest angles of the whole set. The magenta ellipse defines the pupil contour, while the green one is the limbus. The red dots represent the points used to compute the limbus ellipse equation, while the white ones are those removed for their possible belonging to the eyelashes or to wrong estimation of the edge. The blue lines

with *λ*1, *λ*2, *λ*3 ∈ **R**. This transformation from **Z** to **Σ** is just a change of basis, and thus **Σ** may be expressed as **Σ** = **R**−**1ZR**, where **R** is the matrix changing between the actual orthonormal basis to a new one, formed by three eigenvectors of **Z**. The columns of **R** are the components of the eigenvectors **e1**, **e2**, **e3**, and the elements *λ*1, *λ*2, *λ*3 of the diagonal of **Σ** are the associated eigenvalues. For the Sylvester's Law of inertia the signature of **Σ** is equal to that of **Z**, (−, +, +); thus only one eigenvalue is negative, and the other are positive. We assume *λ*3 < 0 and *λ*2 > *λ*1 > 0. If we apply the transformation matrix **R** to **QZ**(**x**), we obtain:

**QZ**(**Rx**)=(**Rx**)*T***Z**(**Rx**)

= **x***T***R***T***ZRx** = **x***T***Σx**

represent the parabolas used to remove the possible eyelashes points.

Since it is not possible to measure with the desired precision the relative position between the center of the eye and nodal point of the camera, and particularly the relative orientation, we developed a calibration procedure with the purpose of estimating the implicit parameters. The subject is asked to fixate a calibration grid of 25 points in a known position. The grid of points is designed to make the subject fixate with an azimuth angle between −30◦ and 30◦ with steps of 15◦, and with an elevation angle between −20◦ and 20◦ with steps of 10◦.

The calibration procedure is based on a functional whose minimization provides: (1) the eye position respect to the camera, (2) the camera orientation respect to a fixed reference frame, (4) the radius of the eye.

#### **6.3 Estimation of the gaze direction**

In order to have a validation of the algorithm, the estimation of the fixation angle is computed over a different set of points respect to the calibration grid. The grid of point used for the test is designed to make the subject fixate with an azimuth angle between −20◦ and 20◦ with steps of 5◦, and with an elevation angle between −10◦ and 10◦ with steps of 5◦.

procedure. In fact, considering to calibrate the algorithm over the range of interest for the tracking of the eye, the pupil is searched in an area where it is likely to be, preventing to detect the initial point on the glasses frame or on other dark regions of the image. Moreover, since the system is not endowed with a frontal scene camera, it comes to be more difficult both to calibrate correctly the algorithm and to test it. Currently for the calibration, the subject is posed manually in the desired position respect to the grid, without any chin rest, and she/he is asked to remain steady all along the procedure. Without any visual feedback from where the subject is fixating, any movement between the subject and the grid (due to undesired rotations and translations of the head, or to physiologic nystagmus) becomes an unpredictable and meaningful source of error. The next steps of our research are to implement of a more comfortable and precise calibration procedure, as through a chin rest or a scene camera, and

The Perspective Geometry of the Eye: Toward Image-Based Eye-Tracking 27

In conclusion, the proposed method, resorting on visible and salient features, like pupil and limbus, and exploiting the known geometry of the structure of the eye, is able to provide a reliable segmentation of the iris that can be in principle used both for non-invasive and

Portions of the research in this paper use the CASIA-IrisV4 collected by the Chinese Academy

This work has been partially supported by the Italian MIUR (PRIN 2008) project "Bio-inspired models for the control of robot ocular movements during active vision and 3D exploration".

Al-Sharadqah, A. & Chernov, N. (2009). Error analysis for circle fitting algorithms, *Electron. J.*

Brolly, X. & Mulligan, J. (2004). Implicit calibration of a remote gaze tracker, *IEEE Conference*

Chernov, N. & Lesort, C. (2004). Statistical efficiency of curve fitting algorithms, *Comput.*

Chessa, M., Garibotti, M., Canessa, A., Gibaldi, A., Sabatini, S. & Solari, F. (2012). A

Chojnacki, W., Brooks, M., van den Hengel, A. & Gawley, D. (2000). On the fitting of surfaces to data with covariances, *IEEE Trans. Patt. Anal. Mach. Intell.* 22(11): 1294–1303. Cornsweet, T. & Crane, H. (1973). Accurate two-dimensional eye tracker using first and fourth

Crane, H. & Steele, C. (1978). Accurate three-dimensional eyetracker, *J. Opt. Soc. Am.*

Daugman, J. (1993). High confidence visual recognition of persons by a test of statistical

independence, *Pattern Analysis and Machine Intelligence, IEEE Transactions on*

stereoscopic augmented reality system for the veridical perception of the 3D scene layout, *International Conference on Computer Vision Theory and Applications (VISAPP*

*on CVPR Workshop on Object Tracking Beyond the Visible Spectrum*.

to extend the system from monocular to binocular tracking.

low-cost eye tracking and for iris recognition applications.

of Sciences' Institute of Automation (CASIA).

CASIA-IrisV1 (2010). *http://biometrics.idealtest.org*.

*Statist. Data Anal.* 47: 713-728.

purkinje images, *J. Opt. Soc. Am.* 63(8): 921–928.

*Stat.* 3: 886–911.

*2012)*.

17(5): 691–705.

15(11): 1148–1161.

**7.1 Acknowledgment**

**8. References**

The error is measured as the angle between the estimated gaze direction and the actual direction of the calibration points. Over the whole set of 45 points, the algorithm is able to provide a mean error of ≈ 0.6◦.

Fig. 7. Azimuthal (horizontal) and elevation (vertical) angles of the grid of fixational points (blue squares), with respect to the angles estimated by the proposed algorithm (red crosses).

#### **7. Discussion and conclusion**

We developed a novel approach for iris segmentation and eye tracking that resorts on the geometrical characteristics of the projection of the eye on the image plane.

Once that the pupil center is roughly located and the ellipse that describes the pupil is fitted, the parameters of the pupil ellipse can be exploited to improve the search of the limbus. We developed a transformation of the image to an elliptical domain, that is shaped by the pupil, in order to transform the limbus in a straight line, thus easier to be detected. The points that do not belong to the limbus are removed considering that the border of the superior and inferior eyelids is well described by two parabolas intersecting at the eyelids intersections. The similarity of the projections of iris and pupil allows a proper segmentation even if large parts of the iris are occluded by eyelids. We developed a method that takes into account the orientation and eccentricity of the pupils ellipse in order to fit the limbus ellipse. The iris segmentation algorithm is able to work both on an iris image database and on the images acquired by our system. Since the limbus can be considered a perfect circle oriented in 3D with respect to the image plane, its imaged ellipse is used to compute the gaze direction finding the orientation in space of the circle that projects the fitted ellipse.

Even though the iris segmentation demonstrates a good effectiveness in a large variety of cases and a good robustness to perturbations due to reflections and glasses, the gaze tracking part is in a preliminary implementation, and many improvements can be implemented in the current algorithm. In order to restrain the wrong matching of the pupil center, the pupil search area can be constrained to a circle defined by the pupil points found during the calibration procedure. In fact, considering to calibrate the algorithm over the range of interest for the tracking of the eye, the pupil is searched in an area where it is likely to be, preventing to detect the initial point on the glasses frame or on other dark regions of the image. Moreover, since the system is not endowed with a frontal scene camera, it comes to be more difficult both to calibrate correctly the algorithm and to test it. Currently for the calibration, the subject is posed manually in the desired position respect to the grid, without any chin rest, and she/he is asked to remain steady all along the procedure. Without any visual feedback from where the subject is fixating, any movement between the subject and the grid (due to undesired rotations and translations of the head, or to physiologic nystagmus) becomes an unpredictable and meaningful source of error. The next steps of our research are to implement of a more comfortable and precise calibration procedure, as through a chin rest or a scene camera, and to extend the system from monocular to binocular tracking.

In conclusion, the proposed method, resorting on visible and salient features, like pupil and limbus, and exploiting the known geometry of the structure of the eye, is able to provide a reliable segmentation of the iris that can be in principle used both for non-invasive and low-cost eye tracking and for iris recognition applications.

#### **7.1 Acknowledgment**

26 Will-be-set-by-IN-TECH

The error is measured as the angle between the estimated gaze direction and the actual direction of the calibration points. Over the whole set of 45 points, the algorithm is able to

−15 −10 −5 0 5 10 15 20

Azimuth [deg]

Fig. 7. Azimuthal (horizontal) and elevation (vertical) angles of the grid of fixational points (blue squares), with respect to the angles estimated by the proposed algorithm (red crosses).

We developed a novel approach for iris segmentation and eye tracking that resorts on the

Once that the pupil center is roughly located and the ellipse that describes the pupil is fitted, the parameters of the pupil ellipse can be exploited to improve the search of the limbus. We developed a transformation of the image to an elliptical domain, that is shaped by the pupil, in order to transform the limbus in a straight line, thus easier to be detected. The points that do not belong to the limbus are removed considering that the border of the superior and inferior eyelids is well described by two parabolas intersecting at the eyelids intersections. The similarity of the projections of iris and pupil allows a proper segmentation even if large parts of the iris are occluded by eyelids. We developed a method that takes into account the orientation and eccentricity of the pupils ellipse in order to fit the limbus ellipse. The iris segmentation algorithm is able to work both on an iris image database and on the images acquired by our system. Since the limbus can be considered a perfect circle oriented in 3D with respect to the image plane, its imaged ellipse is used to compute the gaze direction finding the

Even though the iris segmentation demonstrates a good effectiveness in a large variety of cases and a good robustness to perturbations due to reflections and glasses, the gaze tracking part is in a preliminary implementation, and many improvements can be implemented in the current algorithm. In order to restrain the wrong matching of the pupil center, the pupil search area can be constrained to a circle defined by the pupil points found during the calibration

geometrical characteristics of the projection of the eye on the image plane.

orientation in space of the circle that projects the fitted ellipse.

provide a mean error of ≈ 0.6◦.

Elevation

**7. Discussion and conclusion**

−10 −8 −6 −4 −2 0 2 4 6 8 10

 [deg]

> Portions of the research in this paper use the CASIA-IrisV4 collected by the Chinese Academy of Sciences' Institute of Automation (CASIA).

> This work has been partially supported by the Italian MIUR (PRIN 2008) project "Bio-inspired models for the control of robot ocular movements during active vision and 3D exploration".

#### **8. References**

Leavers, V. (1992). *Shape Detection in Computer Vision Using the Hough Transform*,

The Perspective Geometry of the Eye: Toward Image-Based Eye-Tracking 29

Leedan, Y. & Meer, P. (2000). Heteroscedastic regression in computer vision: Problems with

Li, D., Winfield, D. & Parkhurst, D. (2005). Starburst: A hybrid algorithm for video-based

Lin, C., Chen, H., Lin, C., Yeh, M. & Lin, S. (2005). Polar coordinate mapping method for an

Mackworth, N. & Thomas, E. (1962). Head-mounted eye-marker camera, *J. Opt. Soc. Am.*

Matei, B. & Meer, P. (2006). Estimation of nonlinear errors-in-variables models for computer vision applications, *IEEE Trans. Patt. Anal. Mach. Intell.* 28(10): 1537–1552. Matsumoto, Y. & Zelinsky, A. (2000). An algorithm for real-time stereo vision implementation

Morimoto, C., Koons, D., Amir, A. & Flickner, M. (2000). Pupil detection and tracking using

Ohno, T., Mukawa, N. & Yoshikawa, A. (2002). Freegaze: a gaze tracking system for everydaygaze interaction, *Eye Tracking Research and Applications Symposium*. Parkhurst, D. & Niebur, E. (2004). A feasibility test for perceptually adaptive level of detail

Porrill, J. (1990). Fitting ellipses and predicting confidence envelopes using a bias corrected

Pratt, V. (1987). Direct least-squares fitting of algebraic surfaces, *Computer Graphics* 21: 145–152. Rahib, A. & Koray, A. (2009). Neural network based biometric personal identification with fast iris segmentation, *International Journal of Control, Automation and Systems* 7: 17–23. Reulen, J., Marcus, J., Koops, D., de Vries, F., Tiesinga, G., Boshuizen, K. & Bos, J. (1988).

Robinson, D. A. (1963). A method of measuring eye movemnent using a scleral search coil in a magnetic field, *Bio-medical Electronics, IEEE Transactions on* 10(4): 137–145. Rosin, P. (1993). Ellipse fitting by accumulating five-point fits, *Pattern Recognition Letters*

Rosin, P. L. & West, G. A. W. (1995). Nonparametric segmentation of curves into various

multiple light sources, *Image and Vision Computing* 18(4): 331–335. Nishino, K. & Nayar, S. (2004). Eyes for relighting, *ACM SIGGRAPH* 23(3): 704–711.

*in graphics and visualization*, ACM, New York, NY, USA, pp. 49–56.

kalman filter, *Image Vision and Computing* 8(1): 1140–1153.

*Engineering and Computing* 26: 20–26.

representations, *IEEE Trans. PAMI* 17: 1140–1153.

*2000. Proceedings. Fourth IEEE International Conference on*, pp. 499 –504. Mäenpää, T. (2005). An iterative algorithm for fast iris detection, *in* S. Li, Z. Sun, T. Tan,

eye tracking combining feature-based and model-based approaches, *Computer Vision and Pattern Recognition - Workshops, 2005. CVPR Workshops. IEEE Computer Society*

improved infrared eye-tracking system, *Journal of Biomedical Engineering-Applications,*

of head pose and gaze direction measurement, *Automatic Face and Gesture Recognition,*

S. Pankanti, G. Chollet & D. Zhang (eds), *Advances in Biometric Person Authentication*, Vol. 3781 of *Lecture Notes in Computer Science*, Springer Berlin/Heidelberg,

rendering on desktop systems, *Proceedings of the 1st Symposium on Applied perception*

Precise recording of eye movement: the iris technique part 1, *Medical and Biological*

bilinear constraint, *Int. J. Comput. Vision.* 37(2): 127–150.

Springer-Verlag.

52(6): 713–716.

pp. 127–134.

14: 661–699.

*Conference on*, pp. 79-86.

*Basis and Communicatitons* 17(3): 141–146.

28 Will-be-set-by-IN-TECH

Davé, R. N. & Bhaswan, K. (1992). Nonparametric segmentation of curves into various

Eizenman, M., Frecker, R. & Hallett, P. (1984). Precise non-contacting measurement of eye

Ferreira, A., Lourenço, A., Pinto, B. & Tendeiro, J. (2009). Modifications and improvements on

Fitts, P. M., Jones, R. E. & Milton, J. L. (1950). Eye Movements of Aircraft Pilots during Instrument-Landing Approaches., *Aeronautical Engineering Review* (2(9)): 24–29. Fitzgibbon, A. W., Pilu, M. & Fischer, R. (1996). Direct least squares fitting of ellipses, *Proc. of the 13th International Conference on Pattern Recognition*, Vienna, pp. 253–257. Forsyth, D., Mundy, J., Zisserman, A., Coelho, C., Heller, A. & Rothwell, C. (1991). Invariant

Gath, I. & Hoory, D. (1995). Fuzzy clustering of elliptic ring-shaped clusters, *Pattern*

Halir, R. & Flusser, J. (1998). Numerically stable direct least squares fitting of ellipses,

Harker, M., O'Leary, P. & Zsombor-Murray, P. (2008). Direct type-specific conic fitting and

Hartridge, H. & Thompson, L. (1948). Methods of investigating eye movements, *British Journal*

Haslwanter, T. (1995). Mathematics of three-dimensional eye rotations, *Vision Res.*

Jacob, R. J. K. & Karn, K. S. (2003). Eye Tracking in Human-Computer Interaction and Usability

Kanatani, K. (1996). *Statistical Optimization for Geometric Computation: Theory and Practice*,

Kanatani, K. (1998). Cramer-rao lower bounds for curve fitting, *Graph. Models Image Proc.*

Kanatani, K. & Sugaya, Y. (2007). Performance evaluation of iterative geometric fitting

Kaufman, A., Bandopadhay, A. & Shaviv, B. (1993). An eye tracking computer user interface,

Kyung-Nam, K. & Ramakrishna, R. (1999). Vision-based eye-gaze tracking for human

*Proceedings. 1999 IEEE International Conference on*, Vol. 2, pp. 324–329. Labati, R. D. & Scotti, F. (2010). Noisy iris segmentation with boundary regularization and

reflections removal, *Image and Vision Computing* 28(2): 270 – 277.

*Virtual Reality, 1993. Proceedings., IEEE 1993 Symposium on Research Frontiers in*,

computer interface, *Systems, Man, and Cybernetics, 1999. IEEE SMC '99 Conference*

Research: Ready to Deliver the Promises, *The Mind's eye: Cognitive The Mind's Eye:*

eigenvalue bias correction, *Image Vision and Computing* 26: 372–381.

*Cognitive and Applied Aspects of Eye Movement Research* pp. 573–603.

Javal, E. (1879). Essai sur la Physiologie de la Lecture, *Annales D'Oculistique* .

Elsevier Science Inc., New York, NY, USA.

algorithms, *Comp. Stat. Data Anal.* 52(2): 1208–1222.

descriptors for 3-D object recognition and pose, *IEEE Trans. Patt. Anal. Mach. Intell.*

*Sixth International Conference in Central Europe on Computer Graphics and Visualization*,

Dodge, R. & Cline, T. (1901). The angle velocity of eye movements, *Psychological Review*. Duchowski, A. (2002). A breadth-first survey of eye-tracking applications, *Behav. Res. Methods.*

movements using the corneal reflex, *Vision Research* 24(2): 167–174.

representations, *IEEE Trans. Neural Networks* 3: 643–662.

iris recognition, *BIOSIGNALS09*, Porto, Portugal.

*Instrum. Comput.* 34(4): 455–470.

*Recognition Letters* 16: 727–741.

13(10): 971–991.

pp. 59–108.

*of Ophthalmology* .

35(12): 1727–1739.

60: 93-99.

pp. 120–121.

**2**

*Greece*

**Feature Extraction Based on Wavelet**

**Moments and Moment Invariants in**

G.A. Papakostas, D.E. Koulouriotis and V.D. Tourassis

Recently, there has been an increasing interest on modern machine vision systems for industrial and commercial purposes. More and more products are introduced in the market, which are making use of visual information captured by a camera in order to perform a specific task. Such machine vision systems are used for detecting and/or recognizing a face in an unconstrained environment for security purposes, for analysing the emotional states of a human by processing his facial expressions or for providing a vision based interface in the

In almost all the modern machine vision systems there is a common processing procedure called *feature extraction*, dealing with the appropriate representation of the visual information. This task has two main objectives simultaneously, the compact description of the useful information by a set of numbers (features), by keeping the dimension as low as possible.

Image moments constitute an important *feature extraction method* (FEM) which generates high discriminative features, able to capture the particular characteristics of the described pattern, which distinguish it among similar or totally different objects. Their ability to fully describe an image by encoding its contents in a compact way makes them suitable for many disciplines of the engineering life, such as image analysis (Sim et al., 2004), image watermarking (Papakostas et al., 2010a) and pattern recognition (Papakostas et al., 2007,

Among the several moment families introduced in the past, the orthogonal moments are the most popular moments widely used in many applications, owing to their orthogonality property that comes from the nature of the polynomials used as kernel functions, which they constitute an orthogonal base. As a result, the orthogonal moments have minimum information redundancy meaning that different moment orders describe

In order to use the moments to classify visual objects, they have to ensure high recognition rates for all possible object's orientations. This requirement constitutes a significant operational feature of each modern pattern recognition system and it can be satisfied during

**1. Introduction**

2009a, 2010b).

different parts of the image.

context of the human computer interaction (HCI) etc..

*Department of Production Engineering and Management*

**Machine Vision Systems**

*Democritus University of Thrace,*

## **Feature Extraction Based on Wavelet Moments and Moment Invariants in Machine Vision Systems**

G.A. Papakostas, D.E. Koulouriotis and V.D. Tourassis *Democritus University of Thrace, Department of Production Engineering and Management Greece*

### **1. Introduction**

30 Will-be-set-by-IN-TECH

30 Human-Centric Machine Vision

Ryan, W., Woodard, D., Duchowski, A. & Birchfield, S. (2008). Adapting starburst for elliptical

Shackel, B. (1960). Note on mobile eye viewpoint recording, *J. Opt. Soc. Am.* 50(8): 763–768. Stark, L., Vossius, G. & Young, L. R. (1962). Predictive control of eye tracking movements,

Stiefelhagen, R. & Yang, J. (1997). Gaze tracking for multimodal human-computer interaction,

Taubin, G. (1991). Estimation of planar curves, surfaces and nonplanar space curves defined

Wang, J. & Sung, E. (2002). Study on eye gaze estimation, *IEEE Transactions on Systems,Man*

Wang, J., Sung, E. & Venkateswarlu, R. (2003). Eye gaze estimation from a single image of one

Werman, M. & Geyzel, G. (1995). Fitting a second degree curve in the presence of error, *IEEE*

Wu, W. & Wang, M. (1993). Elliptical object detection by using its geometric properties, *Pattern*

Yin, R., Tam, P. & Leung, N. (1992). Modification of hough transform for circles and ellipses

Yuen, H., Illingworth, J. & Kittler, J. (1989). Detecting partially occluded ellipses using the

Zhu, J. & Yang, J. (2002). Subpixel eye gaze tracking, *IEEE Conference on Automatic Faceand*

*Human Factors in Electronics, IRE Transactions on* 3(2): 52–57.

Tinker, M. (1963). *Legibility of Print*, Iowa State University, Ames, IA, USA.

*IEEE International Conference on*, pp. 1–7.

*Conference on*, Vol. 4, pp. 2617–2620.

*and Cybernetics* 32(3): 332-350.

*Recognition* 26: 1499–1509.

*Gesture Recognition*.

pp. 136 –143.

*Trans. Patt. Anal. Mach. Intell.* 13: 1115–1138.

*Trans. Patt. Anal. Mach. Intell.* 17(2): 207–211.

Yarbus, A. (1959). *Eye movements and vision*, Plenum Press, New York.

detection using 2-d array, *Pattern Recognition* 25: 1007–1022.

hough transform, *Image Vision and Computing* 7(1): 31–37.

iris segmentation, *Biometrics: Theory, Applications and Systems, 2008. BTAS 2008. 2nd*

*Acoustics, Speech, and Signal Processing, 1997. ICASSP-97., 1997 IEEE International*

by implicit equations, with applications to edge and range image segmentation, *IEEE*

eye, *Computer Vision, 2003. Proceedings. Ninth IEEE International Conference on*, Vol. 1,

Recently, there has been an increasing interest on modern machine vision systems for industrial and commercial purposes. More and more products are introduced in the market, which are making use of visual information captured by a camera in order to perform a specific task. Such machine vision systems are used for detecting and/or recognizing a face in an unconstrained environment for security purposes, for analysing the emotional states of a human by processing his facial expressions or for providing a vision based interface in the context of the human computer interaction (HCI) etc..

In almost all the modern machine vision systems there is a common processing procedure called *feature extraction*, dealing with the appropriate representation of the visual information. This task has two main objectives simultaneously, the compact description of the useful information by a set of numbers (features), by keeping the dimension as low as possible.

Image moments constitute an important *feature extraction method* (FEM) which generates high discriminative features, able to capture the particular characteristics of the described pattern, which distinguish it among similar or totally different objects. Their ability to fully describe an image by encoding its contents in a compact way makes them suitable for many disciplines of the engineering life, such as image analysis (Sim et al., 2004), image watermarking (Papakostas et al., 2010a) and pattern recognition (Papakostas et al., 2007, 2009a, 2010b).

Among the several moment families introduced in the past, the orthogonal moments are the most popular moments widely used in many applications, owing to their orthogonality property that comes from the nature of the polynomials used as kernel functions, which they constitute an orthogonal base. As a result, the orthogonal moments have minimum information redundancy meaning that different moment orders describe different parts of the image.

In order to use the moments to classify visual objects, they have to ensure high recognition rates for all possible object's orientations. This requirement constitutes a significant operational feature of each modern pattern recognition system and it can be satisfied during

Feature Extraction Based on Wavelet Moments and Moment Invariants in Machine Vision Systems 33

Table 1. Orthogonal moments' characteristics.

the feature extraction stage, by making the moments invariant under the basic geometric transformations of rotation, scaling and translation.

The most well known orthogonal moment families are: Zernike, Pseudo-Zernike, Legendre, Fourier-Mellin, Tchebichef, Krawtchouk, with the last two ones belonging to the discrete type moments since they are defined directly to the image coordinate space, while the first ones are defined in the continuous space.

Another orthogonal moment family that deserves special attention is the *wavelet moments*  that use an orthogonal wavelet function as kernel. These moments combine the advantages of the wavelet and moment analyses in order to construct moment descriptors with improved pattern representation capabilities (Feng et al., 2009).

This chapter discusses the main theoretical aspects of the wavelet moments and their corresponding invariants, while their performance in describing and distinguishing several patterns in different machine vision applications is studied experimentally.

### **2. Orthogonal image moments**

A general formulation of the *(n+m)th* order orthogonal image moment of a *NxN* image with intensity function *f(x,y)* is given as follows:

$$M\_{nm} = \text{NF} \times \sum\_{i=1}^{N} \sum\_{j=1}^{N} \text{Kernel}\_{nm} \left( \mathbf{x}\_{i\prime} \boldsymbol{y}\_{j} \right) f\left( \mathbf{x}\_{i\prime} \boldsymbol{y}\_{j} \right) \tag{1}$$

where *Kernelnm(.)* corresponds to the moment's kernel consisting of specific polynomials of order *n* and repetition *m*, which constitute the orthogonal basis and *NF* is a normalization factor. The type of Kernel's polynomial gives the name to the moment family by resulting to a wide range of moment types. Based on the above equation (1) the image moments are the projection of the intensity function *f(x,y)* of the image on the coordinate system of the kernel's polynomials.

The first introduction of orthogonal moments in image analysis, due to Teague (Teague, 1980), made use of Legendre and Zernike moments in image processing. Other families of orthogonal moments have been proposed over the years, such as Pseudo-Zernike, Fourier-Mellin etc. moments, which better describe the image in process and ensure robustness under arbitrarily intense noise levels.

However, these moments present some approximation errors due to the fact that the kernel polynomials are defined in a continuous space and an approximated version of them is used in order to compute the moments of an image. This fact is the source of an approximation error (Liao & Pawlak, 1998) which affects the overall properties of the derived moments and mainly their description abilities. Moreover, some of the above moments are defined inside the unit disc, where their polynomials satisfy the orthogonality condition. Therefore, a prior coordinates' transformation is required so that the image coordinates lie inside the unit disc. This transformation is another source of approximation error (Liao & Pawlak, 1998) that further degrades the moments' properties.

The following Table 1, summarizes the main characteristics of the most used moment families.

Table 1. Orthogonal moments' characteristics.

32 Human-Centric Machine Vision

the feature extraction stage, by making the moments invariant under the basic geometric

The most well known orthogonal moment families are: Zernike, Pseudo-Zernike, Legendre, Fourier-Mellin, Tchebichef, Krawtchouk, with the last two ones belonging to the discrete type moments since they are defined directly to the image coordinate space, while the first

Another orthogonal moment family that deserves special attention is the *wavelet moments*  that use an orthogonal wavelet function as kernel. These moments combine the advantages of the wavelet and moment analyses in order to construct moment descriptors with

This chapter discusses the main theoretical aspects of the wavelet moments and their corresponding invariants, while their performance in describing and distinguishing several

A general formulation of the *(n+m)th* order orthogonal image moment of a *NxN* image with

*nm nm i j i j*

where *Kernelnm(.)* corresponds to the moment's kernel consisting of specific polynomials of order *n* and repetition *m*, which constitute the orthogonal basis and *NF* is a normalization factor. The type of Kernel's polynomial gives the name to the moment family by resulting to a wide range of moment types. Based on the above equation (1) the image moments are the projection of the intensity function *f(x,y)* of the image on the coordinate system of the

The first introduction of orthogonal moments in image analysis, due to Teague (Teague, 1980), made use of Legendre and Zernike moments in image processing. Other families of orthogonal moments have been proposed over the years, such as Pseudo-Zernike, Fourier-Mellin etc. moments, which better describe the image in process and ensure robustness

However, these moments present some approximation errors due to the fact that the kernel polynomials are defined in a continuous space and an approximated version of them is used in order to compute the moments of an image. This fact is the source of an approximation error (Liao & Pawlak, 1998) which affects the overall properties of the derived moments and mainly their description abilities. Moreover, some of the above moments are defined inside the unit disc, where their polynomials satisfy the orthogonality condition. Therefore, a prior coordinates' transformation is required so that the image coordinates lie inside the unit disc. This transformation is another source of approximation error (Liao & Pawlak, 1998) that

The following Table 1, summarizes the main characteristics of the most used moment

*M NF Kernel x y f x y*

1 1

*i j*

= =

*N N*

( )( )

= ×∑∑ (1)

, ,

transformations of rotation, scaling and translation.

improved pattern representation capabilities (Feng et al., 2009).

patterns in different machine vision applications is studied experimentally.

ones are defined in the continuous space.

**2. Orthogonal image moments**

under arbitrarily intense noise levels.

further degrades the moments' properties.

kernel's polynomials.

families.

intensity function *f(x,y)* is given as follows:

Feature Extraction Based on Wavelet Moments and Moment Invariants in Machine Vision Systems 35

According to this method and by applying coordinates normalization (Rothe et al., 1996) the

( ) ( )

1, , , 2

where (,) *x y* are the coordinates of the image's centroid, *GMnm* are the geometric moments

*N N n m*

( ) ( ) *N N n m*

additional modifications (Mukundan & Ramakrishnan, 1998) have to be performed in order

By expressing each moment family in terms of geometric moment invariants the corresponding invariants can be derived. For example, Zernike moments are expressed

*nm nmk k i j i j*

−− +

0 0

where n is a non-negative integer and m is a non zero integer subject to the constraints n-

*w with i s k m*

⎧− > = ⎨ =− = − + ≤ ⎩

, 0 1 1, , 0 2

*n k*

1 ! 2 !!! 222

*x x y y f xy*

*x y GM x y f x y* 1 1

0 0

= =

(,)

which are translation invariant. The value of angle θ is limited to *o o* − ≤≤ 45 45

*n sm i*

*n s m ZMI B w GMI*

*km i j*

*n k*

( ) 2

−

⎛ ⎞ + − ⎜ ⎟ ⎝ ⎠ = ⎛ ⎞⎛ ⎞⎛ ⎞ −+ − ⎜ ⎟⎜ ⎟⎜ ⎟ ⎝ ⎠⎝ ⎠⎝ ⎠

= ==

− −

*n m GM GM*

= ⎡ ⎤ − +− ⎣ ⎦

θ

⎡ ⎤ − −− ⎣ ⎦

*x y GM GM*

( ) ( )

θ

*y y x x fxy*

10 01 00 00

( )

= ∑ ∑ (5)

= −− ∑ ∑ (6)

*i j* 2 ,2

( )

+ ⎛ ⎞⎛ ⎞ = ⎜ ⎟⎜ ⎟ ⎝ ⎠⎝ ⎠ ∑ ∑∑ (7)

,

cos sin

cos sin ( , )

 θ *m*

 θ

(3)

(4)

θ

(8)

and

*M N n*

*Geometric Moment Invariants (GMIs)* of the following form, are constructed:

*x y*

− − − = =

00

+

γ

*GMI GM x x y y*

1 11 20 02

μ μ

⎛ ⎞ = ⎜ ⎟ − ⎝ ⎠

*nm*

1 1

− −

0 0

= =

*x y*

μ

= += =

1 2 tan 2

−

1 1

0 0

∑ ∑

*nm*

γ

θ

*nm*

θ.

> 1 π

*i m*

*i m*

( )

*B nk km km*

(Wee & Paramesran, 2007) in terms of GMIs as follows:

*nmk*

μ

and *μnm* are the central moments defined as:

to extent θ into the range *o o* 0 360 ≤ ≤

with

Apart from some remarkable attempts to compute the theoretical image moment values (Wee & Paramesran, 2007), new moment families with discrete polynomial kernels (Tchebichef and Krawtchouk moments) have been proposed, which permit the direct computation of the image moments in the discrete domain.

It is worth pointing out that the main subject of this chapter is the introduction of a particular moment family, the *wavelet moments* and the investigation of their classification capabilities as compared to the traditional moment types.

However, before introducing the wavelet moments, it is useful to discuss two important properties of the moments that determine their utility in recognizing patterns.

### **2.1 Information description**

As it has already been mentioned in the introduction, the moments have the ability to carry information of an image with minimum redundancy, while they are capable to enclose distinctive information that uniquely describes the image's content. Due to these properties, once a finite number of moments up to a specific order *nmax* is computed, the original image can be reconstructed by applying a simple formula, inverse to (1), of the following form:

$$\hat{f}\left(\mathbf{x},\mathbf{y}\right) = \sum\_{n=0}^{n\_{\text{max}}} \sum\_{m=0}^{n} Kernel\_{nm}\left(\mathbf{x},\mathbf{y}\right) \mathcal{M}\_{nm} \tag{2}$$

where *Kernelnm(.)* is the same kernel of (1) used to compute moment *Mnm*.

Theoretically speaking, if one computes all image moments and uses them in (2), the reconstructed image will be identical to the original one with minimum reconstruction error.

### **2.2 Invariant description**

Apart from the ability of the moments to describe the content of an image in a statistical fashion and to reconstruct it perfectly (orthogonal moments) according to (2), they can also be used to distinguish a set of patterns belonging to different categories (classes). This property makes them suitable for many artificial intelligence applications such as biometrics, visual inspection or surveillance, quality control, robotic vision and guidance, biomedical diagnosis, mechanical fault diagnosis etc. However, in order to use the moments to classify visual objects, they have to ensure high recognition rates for all possible object's orientations. This requirement constitutes a significant operational feature of each modern pattern recognition system and it can be satisfied during the feature extraction stage, where discriminative information of the objects is retrieved.

Mainly, two methodologies used to ensure invariance under common geometric transformations such as rotation, scaling and translation, either by image coordinates normalization and description through the geometric moment invariants (Mukundan & Ramakrishnan, 1998; Zhu et al., 2007) or by developing new computation formulas which incorporate these useful properties inherently (Zhu et al., 2007).

However, the former strategy is usually applied for deriving the moment invariants of each moment family, since it can be applied in each moment family in a similar way.

According to this method and by applying coordinates normalization (Rothe et al., 1996) the *Geometric Moment Invariants (GMIs)* of the following form, are constructed:

$$\begin{aligned} \text{GMI}\_{nm} &= \text{GM}\_{00}^{-\prime} \sum\_{\mathbf{x}=0}^{M-1} \sum\_{\mathbf{y}=0}^{N-1} \left[ (\mathbf{x} - \overline{\mathbf{x}}) \cos \theta + (\mathbf{y} - \overline{\mathbf{y}}) \sin \theta \right]^{\prime} \\ &\qquad \left[ (\mathbf{y} - \overline{\mathbf{y}}) \cos \theta - (\mathbf{x} - \overline{\mathbf{x}}) \sin \theta \right]^{m} f(\mathbf{x}, \mathbf{y}) \end{aligned} \tag{3}$$

with

34 Human-Centric Machine Vision

Apart from some remarkable attempts to compute the theoretical image moment values (Wee & Paramesran, 2007), new moment families with discrete polynomial kernels (Tchebichef and Krawtchouk moments) have been proposed, which permit the direct

It is worth pointing out that the main subject of this chapter is the introduction of a particular moment family, the *wavelet moments* and the investigation of their classification

However, before introducing the wavelet moments, it is useful to discuss two important

As it has already been mentioned in the introduction, the moments have the ability to carry information of an image with minimum redundancy, while they are capable to enclose distinctive information that uniquely describes the image's content. Due to these properties, once a finite number of moments up to a specific order *nmax* is computed, the original image can be reconstructed by applying a simple formula, inverse to (1), of the following form:

( ) ( ) *n n*

*f x y Kernel x y M* max 0 0 ˆ , , = =

Theoretically speaking, if one computes all image moments and uses them in (2), the reconstructed image will be identical to the original one with minimum reconstruction error.

Apart from the ability of the moments to describe the content of an image in a statistical fashion and to reconstruct it perfectly (orthogonal moments) according to (2), they can also be used to distinguish a set of patterns belonging to different categories (classes). This property makes them suitable for many artificial intelligence applications such as biometrics, visual inspection or surveillance, quality control, robotic vision and guidance, biomedical diagnosis, mechanical fault diagnosis etc. However, in order to use the moments to classify visual objects, they have to ensure high recognition rates for all possible object's orientations. This requirement constitutes a significant operational feature of each modern pattern recognition system and it can be satisfied during the feature extraction stage, where

Mainly, two methodologies used to ensure invariance under common geometric transformations such as rotation, scaling and translation, either by image coordinates normalization and description through the geometric moment invariants (Mukundan & Ramakrishnan, 1998; Zhu et al., 2007) or by developing new computation formulas which

However, the former strategy is usually applied for deriving the moment invariants of each

moment family, since it can be applied in each moment family in a similar way.

*n m*

where *Kernelnm(.)* is the same kernel of (1) used to compute moment *Mnm*.

discriminative information of the objects is retrieved.

incorporate these useful properties inherently (Zhu et al., 2007).

*nm nm*

= ∑ ∑ (2)

properties of the moments that determine their utility in recognizing patterns.

computation of the image moments in the discrete domain.

capabilities as compared to the traditional moment types.

**2.1 Information description**

**2.2 Invariant description**

$$\begin{aligned} \gamma &= \frac{n+m}{2} + 1, \quad \overline{\chi} = \frac{\text{GM}\_{10}}{\text{GM}\_{00}}, \overline{\chi} = \frac{\text{GM}\_{01}}{\text{GM}\_{00}},\\ \theta &= \frac{1}{2} \tan^{-1} \left( \frac{2\mu\_{11}}{\mu\_{20} - \mu\_{02}} \right) \end{aligned} \tag{4}$$

where (,) *x y* are the coordinates of the image's centroid, *GMnm* are the geometric moments and *μnm* are the central moments defined as:

$$\text{G}M\_{nm} = \sum\_{\mathbf{x}=\mathbf{0}}^{N-1} \sum\_{y=0}^{N-1} \mathbf{x}^n y^m f(\mathbf{x}, y) \tag{5}$$

$$\mu\_{nm} = \sum\_{x=0}^{N-1} \sum\_{y=0}^{N-1} \left(\mathbf{x} - \overline{\mathbf{x}}\right)^{n} \left(y - \overline{y}\right)^{m} f(\mathbf{x}, y) \tag{6}$$

which are translation invariant. The value of angle θ is limited to *o o* − ≤≤ 45 45 θ and additional modifications (Mukundan & Ramakrishnan, 1998) have to be performed in order to extent θ into the range *o o* 0 360 ≤ ≤ θ.

By expressing each moment family in terms of geometric moment invariants the corresponding invariants can be derived. For example, Zernike moments are expressed (Wee & Paramesran, 2007) in terms of GMIs as follows:

$$\text{GZM}\_{nm} = \frac{n+1}{\pi} \sum\_{k=m}^{n} B\_{nmk} \sum\_{i=0}^{s} \sum\_{j=0}^{m} w^{j} \binom{s}{i} \binom{m}{j} \text{GMI}\_{k-2i-j, 2i+j} \tag{7}$$

where n is a non-negative integer and m is a non zero integer subject to the constraints n- |m| even and |m|≤ n and

$$w = \begin{cases} -i & m > 0\\ +i & m \le 0 \end{cases} \quad \text{with} \quad i = \sqrt{-1}, \quad s = \frac{1}{2}(k - m)$$

$$B\_{nmk} = \frac{(-1)^{\frac{(n-k)}{2}} \left(\frac{n+k}{2}\right)!}{\left(\frac{n-k}{2}\right)! \left(\frac{k+m}{2}\right)! \left(\frac{k-m}{2}\right)!} \tag{8}$$

Feature Extraction Based on Wavelet Moments and Moment Invariants in Machine Vision Systems 37

With the above the wavelet basis is constructed by a modified formula of (10) having the form:

 ψ

It has to be noted that the selection of b0 to 0.5 causes oversampling, something which adds significant information redundancy when the wavelet moments are used to reconstruct the initial image, but it doesn't seriously affect their recognition capabilities. Moreover, in order to reduce this affection a feature selection procedure can be applied to keep only the useful

(a)

(b)

Fig. 1. Plots of (a) cubic B-spline and (b) Mexican hat mother wavelet functions.

( ) ( ) *m m mn r rn* /2

0.5 , 0,1,2,3...

*m mm*

0.5 0.5 , 0,1,...,2 +

1

= − 2 2 0.5 (15)

(14)

*m m*

0 0

ψ

0

features, by discarding the redundant ones.

*a a m b b na n n*

= = = = ×× = ×× =

#### **3. Wavelet-based moment descriptors**

In the same way as the continuous radial orthogonal moments such as Zernike, Pseudo-Zernike and Fourier-Mellin ones are defined in a continuous form (9), one can define the *wavelet moments* by replacing the function *gn(r)* with a wavelet basis functions.

$$M\_{nm} = \iint \mathbf{g}\_n(\mathbf{r}) e^{-jm\theta} f\left(\mathbf{r}, \theta\right) r dr d\theta \tag{9}$$

Based on Table 1 it can be deduced that by choosing the appropriate function *gn(r)*, the Zernike, Pseudo-Zernike and Fourier moments are derived. If one chooses wavelet basis functions of the following form

$$
\varphi\_{a,b}(r) = \frac{1}{\sqrt{a}} \varphi\left(\frac{r-b}{a}\right) \tag{10}
$$

where *a b*, ∈ℜ ∈ℜ + are the dilation and translation parameters and ψ ( )⋅ the *mother wavelet* that is used to generate the whole basis.

Two widely used mother wavelet functions are the *cubic B-spline* and the *Mexican hat* functions defined as follows:

$$\text{Cubic B-spline Method } \mathsf{W}\nu(r) = \frac{4d^{n+1}}{\sqrt{2\pi (n+1)}} \sigma\_w \cos \left(2\pi f\_0 (2r - 1)\right) \times e^{\left(\frac{(2r - 1)^2}{2\sigma\_w^2 (n+1)}\right)} \tag{11}$$

where

$$\begin{cases} n = & 3\\ a = 0.697066\\ f\_0 = 0.409177\\ \sigma\_w^2 = 0.561145 \end{cases} \tag{12}$$

$$\text{Maxicam Hat Method} \text{ Moverlet } \psi(r) = \frac{2}{\sqrt{3\sigma}} \pi^{-1/4} \left(1 - \frac{r^2}{\sigma^2}\right) \times e^{\left(-\frac{r^2}{2\sigma^2}\right)} \tag{13}$$

with *σ=1*.

The main characteristic of the above wavelet functions is that by adjusting the *a,b* parameters a basis functions consisting of dilated (scaled) and translated versions of the mother wavelets can be derived.

The graphical presentation of the above two mother wavelet functions for the above set of their parameters, is illustrated in the following Fig.1

Since the a,b parameters are usually discrete (this is mandatory for the case of the resulted moments), a discretization procedure needs to be applied. Such a common method (Shen & Ip, 1999) that also takes into consideration the restriction of *r* ≤ 1 , applies the following relations.

36 Human-Centric Machine Vision

In the same way as the continuous radial orthogonal moments such as Zernike, Pseudo-Zernike and Fourier-Mellin ones are defined in a continuous form (9), one can define the

> () ( ) *jm M g nm n r e f r rdrd* , θ

Based on Table 1 it can be deduced that by choosing the appropriate function *gn(r)*, the Zernike, Pseudo-Zernike and Fourier moments are derived. If one chooses wavelet basis

*a b* ( ) *r b*

 ψ⎛ ⎞ − = ⎜ ⎟

*a a* , 1

Two widely used mother wavelet functions are the *cubic B-spline* and the *Mexican hat*

*a*

π

*w*

ψ

σ

*n a f*0 2 *n*

+

3 0.697066 0.409177

2 1

0.561145

The main characteristic of the above wavelet functions is that by adjusting the *a,b* parameters a basis functions consisting of dilated (scaled) and translated versions of the

The graphical presentation of the above two mother wavelet functions for the above set of

Since the a,b parameters are usually discrete (this is mandatory for the case of the resulted moments), a discretization procedure needs to be applied. Such a common method (Shen & Ip, 1999) that also takes into consideration the restriction of *r* ≤ 1 , applies the following

⎧ = ⎫ ⎪ ⎪ = ⎪ ⎪ ⎨ ⎬ = ⎪ ⎪ ⎪ ⎪ = ⎩ ⎭

*w*

*r fr e*

4 cos 2 2 1

*r e*

= −× ⎜ ⎟

2 1

 π

3

σ

σπ

*r*

ψ

*Cubic B-spline Mother Wavelet* ( ) ( ) ( ) ( )

ψ

*Mexican Hat Mother Wavelet* ( )

their parameters, is illustrated in the following Fig.1

where *a b*, ∈ℜ ∈ℜ + are the dilation and translation parameters and

θ

 θ− = ∫∫ (9)

⎝ ⎠ (10)

ψ

*n n*

*r*

 σ

⎝ ⎠

1/4 2 2

⎜ ⎟ −⎜ ⎟ − ⎝ ⎠ ⎛ ⎞

⎛ ⎞ − ⎜ ⎟ − + ⎜ ⎟ + ⎝ ⎠ = − ×

1 2 1 0

( )⋅ the *mother wavelet*

( ) *w* ( ) *r*

σ

*r*

⎛ ⎞

2 2 2

σ

2 2 2 1

(11)

(12)

(13)

*wavelet moments* by replacing the function *gn(r)* with a wavelet basis functions.

**3. Wavelet-based moment descriptors**

functions of the following form

functions defined as follows:

where

with *σ=1*.

relations.

mother wavelets can be derived.

that is used to generate the whole basis.

$$\begin{aligned} a &= a\_0^{\;\;\;m} = 0.5^m, & m &= 0, 1, 2, 3...\\ b &= b\_0 \times n \times a\_0^{\;\;\;m} = 0.5 \times n \times 0.5^m, & n &= 0, 1, ..., 2^{m+1} \end{aligned} \tag{14}$$

With the above the wavelet basis is constructed by a modified formula of (10) having the form:

$$
\Psi\_{mn}(r) = 2^{m/2} \Psi \left( 2^m r - 0.5n \right) \tag{15}
$$

It has to be noted that the selection of b0 to 0.5 causes oversampling, something which adds significant information redundancy when the wavelet moments are used to reconstruct the initial image, but it doesn't seriously affect their recognition capabilities. Moreover, in order to reduce this affection a feature selection procedure can be applied to keep only the useful features, by discarding the redundant ones.

Fig. 1. Plots of (a) cubic B-spline and (b) Mexican hat mother wavelet functions.

Feature Extraction Based on Wavelet Moments and Moment Invariants in Machine Vision Systems 39

Fig. 2. Six pattern samples of the D1 (1st row), D2 (2nd row), D3 (3rd row) and D4 (4th row)

WMs-2 the Mexican hat mother wavelets respectively).

the samples from the centre classes, having the following form:

The performance of the wavelet moments (WMs) is compared to this of the well-known Zernike (ZMs), Pseudo-Zernike (PZMs), Fourier-Mellin (FMs), Legendre (LMs), Tchebichef (TMs) and Krawtchouk (KMs) ones. In this way, for each dataset, a set of moments up to a specific order per moment family is computed, by resulting to feature vectors of the same length. It is decided to construct feature vectors of 16 moments length which correspond to different order per moment family (ZMs(6th), PZMs(5th), FMs(3rd), LMs(3rd), TMs(3rd), KMs(3rd), WMs(1st) ). Moreover, the wavelet moments are studied under two different configurations in relation to the used mother wavelet (WMs-1 uses the cubic B-spline and

Furthermore, the *Minimum Distance* classifier (Kuncheva, 2004) is used to compute the classification performance of each moment feature vector. This classifier operates by measuring the distance of each sample from the patterns that represent the classes' centre. The sample is decided to belong to the specific class having less distance from its pattern. For the purpose of the experiments the *Euclidean* distance is used to measure the distance of

*Euclidean Distance* ( ) ( )

,

The above formula measures the distance between two vectors the pattern **p**=[p1,p2,p3,…pn] and the sample **s**=[s1,s2,s3,…,sn], which are defined in the Rn space. The following Table 3 and Table 4, summarize the classification rates (18) of the studied moment families for different set of training data used to determine the classes' centres (percent of the entire data

*number of correctclassifiedsamples CRate*

*n*

*i d ps* 2 1

=

*i i*

*totalnumber of samples* = (18)

**p s** = − ∑ (17)

datasets.

– 25%, 50%, 75%, 100%).

Based on the previous analysis and by following the nomenclature of Table 1, the wavelet moments of a *N*x*N* image *f(x,y)* are defined in the unit disc (*r* ≤ 1 ) as follows:

$$\mathcal{W}\_{mnq} = \sum\_{x=1}^{N} \sum\_{y=1}^{N} \Psi\_{mn}(r) e^{-jq\theta} f\left(r, \theta\right) \tag{16}$$

with *r x y y x* 2 2 =+ = , arctan( / ) θ.

The corresponding wavelet moment invariants can be derived by applying the methodologies presented in section 2.2. An important property of the radial moments defined by (9) is that their amplitudes are rotation invariant. The translation invariants are achieved by moving the origin to the mass center (x0,y0) of the image, while scaling invariance is obtained by resizing the image to a fixed size having a predefined area (by a factor *a M area* = 00 / ) by using the zeroth geometric moment order (M00).

The inherent properties of the wavelet moments coming from the wavelet analysis (Strang & Nguyen, 1997), according to which a signal is breaking up into shifted and scaled versions of the base wavelet (mother wavelet), make them appropriate in describing the coarse and fine information of an image. This description has the advantage to study a signal on a *timescale* domain by providing time and frequency (there is a relation between scale and frequency), useful information of the signal simultaneously.

#### **4. Experimental study**

In order to investigate the discrimination power of the wavelet moments, four machine vision experiments have been arranged by using some well-known benchmark datasets. The following Table 2, summarizes the main characteristics of the datasets being used, from different application fields (object, face, facial expression and hand posture recognition). Moreover, Fig. 2, illustrates six pattern samples from each dataset.

Table 2. Characteristics of the benchmark datasets.

Since the wavelet moment invariants are constructed by applying the same methods as in the case of the other moment families and therefore their performance in recognizing geometrical degraded images is highly depended on the representation capabilities of the wavelet moments, it is decided to investigate only the discrimination power of the wavelet moments under invariant conditions.

38 Human-Centric Machine Vision

Based on the previous analysis and by following the nomenclature of Table 1, the wavelet

*W r e f r* 1 1

The corresponding wavelet moment invariants can be derived by applying the methodologies presented in section 2.2. An important property of the radial moments defined by (9) is that their amplitudes are rotation invariant. The translation invariants are achieved by moving the origin to the mass center (x0,y0) of the image, while scaling invariance is obtained by resizing the image to a fixed size having a predefined area (by a

The inherent properties of the wavelet moments coming from the wavelet analysis (Strang & Nguyen, 1997), according to which a signal is breaking up into shifted and scaled versions of the base wavelet (mother wavelet), make them appropriate in describing the coarse and fine information of an image. This description has the advantage to study a signal on a *timescale* domain by providing time and frequency (there is a relation between scale and

In order to investigate the discrimination power of the wavelet moments, four machine vision experiments have been arranged by using some well-known benchmark datasets. The following Table 2, summarizes the main characteristics of the datasets being used, from different application fields (object, face, facial expression and hand posture recognition).

**ID Name Instances D1** COIL (Nene et al., 1996) computer vision 10 12 120

1994) face recognition 40 10 400

Since the wavelet moment invariants are constructed by applying the same methods as in the case of the other moment families and therefore their performance in recognizing geometrical degraded images is highly depended on the representation capabilities of the wavelet moments, it is decided to investigate only the discrimination power of the wavelet

**Classes**

recognition 7 30,29,32,31,30,31,30 213

recognition 10 24 (only the dark

**Instances / Class**

background) 240

**Total**

ψ () ( )

, θ

θ−

= ∑ ∑ (16)

*jq*

moments of a *N*x*N* image *f(x,y)* are defined in the unit disc (*r* ≤ 1 ) as follows:

*N N*

= =

*mnq mn x y*

factor *a M area* = 00 / ) by using the zeroth geometric moment order (M00).

with *r x y y x* 2 2 =+ = , arctan( / ) θ

**4. Experimental study**

**D2** ORL (Samaria & Harter,

**D3** JAFFE (Lyons et al.,

**D4** TRIESCH I (Triesch & von der Malsburg, 1996)

moments under invariant conditions.

Table 2. Characteristics of the benchmark datasets.

1998)

.

frequency), useful information of the signal simultaneously.

Moreover, Fig. 2, illustrates six pattern samples from each dataset.

**Dataset Type Num.**

facial expression

hand posture

Fig. 2. Six pattern samples of the D1 (1st row), D2 (2nd row), D3 (3rd row) and D4 (4th row) datasets.

The performance of the wavelet moments (WMs) is compared to this of the well-known Zernike (ZMs), Pseudo-Zernike (PZMs), Fourier-Mellin (FMs), Legendre (LMs), Tchebichef (TMs) and Krawtchouk (KMs) ones. In this way, for each dataset, a set of moments up to a specific order per moment family is computed, by resulting to feature vectors of the same length. It is decided to construct feature vectors of 16 moments length which correspond to different order per moment family (ZMs(6th), PZMs(5th), FMs(3rd), LMs(3rd), TMs(3rd), KMs(3rd), WMs(1st) ). Moreover, the wavelet moments are studied under two different configurations in relation to the used mother wavelet (WMs-1 uses the cubic B-spline and WMs-2 the Mexican hat mother wavelets respectively).

Furthermore, the *Minimum Distance* classifier (Kuncheva, 2004) is used to compute the classification performance of each moment feature vector. This classifier operates by measuring the distance of each sample from the patterns that represent the classes' centre. The sample is decided to belong to the specific class having less distance from its pattern. For the purpose of the experiments the *Euclidean* distance is used to measure the distance of the samples from the centre classes, having the following form:

$$\text{Euclidean Distance } d \text{(p.s.)} = \sqrt{\sum\_{i=1}^{n} \left(p\_i - s\_i\right)^2} \tag{17}$$

The above formula measures the distance between two vectors the pattern **p**=[p1,p2,p3,…pn] and the sample **s**=[s1,s2,s3,…,sn], which are defined in the Rn space. The following Table 3 and Table 4, summarize the classification rates (18) of the studied moment families for different set of training data used to determine the classes' centres (percent of the entire data – 25%, 50%, 75%, 100%).

$$\text{CRate} = \frac{number \, of \, correct \, classified \, samples}{\text{total number of samples}} \tag{18}$$

Feature Extraction Based on Wavelet Moments and Moment Invariants in Machine Vision Systems 41

moments that are computed up to the pth order is equal to (p+1)2, while in the case of ZMs is (p+1)\*(p+2)/2 due to some constraints. On the other hand the number of wavelet moments that is computed for the pth order is (p+1)2 \* (2p+1+1) due to the third parameter (n) of (16) defined in (14). From this analysis it is obvious that if a common computation algorithm is applied to all the moment families, the time needed to compute the wavelet moments up to

The previous analysis constitutes the first study of the wavelet moments' classification performance in well-known machine vision benchmarks. The experimental results highlight an important weakness of the wavelet moments; the computation of many features for a given order (due to the third parameter), which do not carry enough discriminative information of the patterns. On the other hand, this additional parameter adds an extra degree of freedom to the overall computation which needs to be manipulated appropriately. The usage of a feature selection mechanism can significantly improve the classification capabilities of the wavelet moments by keeping only the useful features from a large pool. In this way, the multiresolution nature of the wavelet analysis can be exploited in order to

Moreover, it has to be noted that a fast and accurate algorithm for the computation of the wavelet moments need to be developed, since their computation overhead is very high, compared to the other moment families, due to the presence of the third configuration

Feng, Z. ; Shang-qian, L.; Da-bao, Wang & Wei, G. (2009). Aircraft recognition in infrared

Kuncheva, L.I. (2004). *Combining pattern classifiers: methods and algorithms*. Wiley-Interscience

Liao, S.X. & Pawlak, M. (1998). On the accuracy of Zernike moments for image analysis. *IEEE Trans. Pattern Analysis and Machine Intelligence*, Vol. 20, No. 12, pp. 1358-1364. Lyons, M.J., Akamatsu, S., Kamachi, M. & Gyoba, J. (1998). Coding facial expressions with

Mukundan, R.; Ong, S.H. & Lee, P.A. (2001). Image analysis by Tchebichef moments. *IEEE*

Nene, S.A.; Nayar, S.K. & Murase, H. (1996). Columbia Object Image Library (COIL-20).

Papakostas, G.A.; Boutalis, Y.S.; Karras, D.A. & Mertzios, B.G. (2007). "A new class of

Zernike moments for computer vision applications. *Information Sciences*, Vol. 177,

Automatic Face and Gesture Recognition (pp. 200-205). Nara, Japan. Mukundan, R. & Ramakrishnan, K.R. (1998). *Moment functions in image analysis*. World

*Transactions on Image Processing*, Vol. 10, No. 9, pp. 1357-1364.

image using wavelet moment invariants. *Image and Vision Computing*, Vol. 27, No. 4,

Gabor wavelets. Proceedings of the 3rd IEEE International Conference on

capture the discriminative information in different discrimination levels.

a specific order (p) is considerable higher.

**5. Discussion – Open issues**

parameter.

**6. References**

pp. 313-318.

Publiching.

Scientific Publisher.

No.13, pp. 2802-2819.

Technical Report No. CUCS-006-96.

Table 3. Classification performance of the moment descriptors.

Table 4. Classification performance of the moment descriptors.

From the above results it is deduced that the percent of the dataset used to determine the classes' centres is crucial to the recognition performance of all the moment families. The performance of the wavelet moments is very low when compared to the other families. This behaviour is justified by the chosen order (1st) that produces less discriminant features. It seems that the existence of the third parameter (n=0,1,…,2m+1) does not add significant discriminative information to the feature vector, compared to that enclosed by the *m* and *q* parameters. As far as the performance of the other moment families is concerned, the experiments show that each moment family behaves differently in each dataset (highest rates: D1(ZMs), D2(LMs), D3(PZMs), D4(PZMs)) with the Pseudo-Zernike moments being the most efficient.

It is worth mentioning that the above rates are not optimized and they can be increased by using a more sophisticated classification scheme (e.g. neural classifier) or by constructing larger or appropriate selected feature vectors.

Besides the classification performance of the compared moment families discussed previously, it is also interesting to analyse their computational load. In almost all the non wavelet moment families (PZMs, FMs, LMs, TMs and KMs) the number of independent moments that are computed up to the pth order is equal to (p+1)2, while in the case of ZMs is (p+1)\*(p+2)/2 due to some constraints. On the other hand the number of wavelet moments that is computed for the pth order is (p+1)2 \* (2p+1+1) due to the third parameter (n) of (16) defined in (14). From this analysis it is obvious that if a common computation algorithm is applied to all the moment families, the time needed to compute the wavelet moments up to a specific order (p) is considerable higher.

### **5. Discussion – Open issues**

40 Human-Centric Machine Vision

 **D1 D2 D3 D4**

**ZMs** 0.6370 0.7235 0.6576 0.6964 0.2004 0.2055 0.0857 0.1014 **PZMs** 0.5912 0.6305 0.6247 0.6657 0.2130 0.2257 0.0906 0.1001 **FMs** 0.5720 0.6000 0.6068 0.6354 0.1837 0.1965 0.0746 0.0872 **LMs** 0.4713 0.5158 0.7770 0.8124 0.2392 0.2547 0.0686 0.0678 **TMs** 0.4688 0.5055 0.7772 0.8073 0.2385 0.2557 0.0689 0.0678 **KMs** 0.5079 0.5915 0.3999 0.4193 0.2090 0.2348 0.0759 0.0823 **WMs – 1** 0.2829 0.2862 0.2252 0.2228 0.1521 0.1616 0.0715 0.0758 **WMs – 2** 0.2723 0.2807 0.2206 0.2219 0.1532 0.1643 0.0682 0.0790

 **D1 D2 D3 D4**

**ZMs** 0.7543 0.8083 0.7289 0.8175 0.2060 0.2723 0.1150 0.2625 **PZMs** 0.6683 0.7417 0.6857 0.7675 0.2385 0.3333 0.1154 0.2958 **FMs** 0.6207 0.7000 0.6396 0.7525 0.2098 0.2723 0.0983 0.2792 **LMs** 0.5457 0.7833 0.8319 0.8975 0.2562 0.3192 0.0681 0.1625 **TMs** 0.5287 0.7833 0.8241 0.8900 0.2719 0.3192 0.0727 0.1583 **KMs** 0.5940 0.7250 0.4206 0.5550 0.2383 0.3146 0.0854 0.2750 **WMs – 1** 0.2887 0.3000 0.2146 0.2425 0.1717 0.1784 0.0844 0.1542 **WMs – 2** 0.2960 0.3083 0.2136 0.2425 0.1702 0.1784 0.0846 0.1500

From the above results it is deduced that the percent of the dataset used to determine the classes' centres is crucial to the recognition performance of all the moment families. The performance of the wavelet moments is very low when compared to the other families. This behaviour is justified by the chosen order (1st) that produces less discriminant features. It seems that the existence of the third parameter (n=0,1,…,2m+1) does not add significant discriminative information to the feature vector, compared to that enclosed by the *m* and *q* parameters. As far as the performance of the other moment families is concerned, the experiments show that each moment family behaves differently in each dataset (highest rates: D1(ZMs), D2(LMs), D3(PZMs), D4(PZMs)) with the Pseudo-Zernike moments being

It is worth mentioning that the above rates are not optimized and they can be increased by using a more sophisticated classification scheme (e.g. neural classifier) or by constructing

Besides the classification performance of the compared moment families discussed previously, it is also interesting to analyse their computational load. In almost all the non wavelet moment families (PZMs, FMs, LMs, TMs and KMs) the number of independent

*25% 50% 25% 50% 25% 50% 25% 50%*

*75% 100% 75% 100% 75% 100% 75% 100%*

**Family Datasets**

Table 3. Classification performance of the moment descriptors.

**Family Datasets**

Table 4. Classification performance of the moment descriptors.

**Moment**

**Moment**

the most efficient.

larger or appropriate selected feature vectors.

The previous analysis constitutes the first study of the wavelet moments' classification performance in well-known machine vision benchmarks. The experimental results highlight an important weakness of the wavelet moments; the computation of many features for a given order (due to the third parameter), which do not carry enough discriminative information of the patterns. On the other hand, this additional parameter adds an extra degree of freedom to the overall computation which needs to be manipulated appropriately. The usage of a feature selection mechanism can significantly improve the classification capabilities of the wavelet moments by keeping only the useful features from a large pool. In this way, the multiresolution nature of the wavelet analysis can be exploited in order to capture the discriminative information in different discrimination levels.

Moreover, it has to be noted that a fast and accurate algorithm for the computation of the wavelet moments need to be developed, since their computation overhead is very high, compared to the other moment families, due to the presence of the third configuration parameter.

### **6. References**

**3**

**A Design for Stochastic Texture Classification**

A texture is a pattern represented on a surface or a structure of an object. Classifying a texture involves pattern recognition. Integration of texture classification is a one step ahead for machine vision application. For the past few decades until now, many researchers from various fields have been trying to develop algorithms to do texture discrimination. Alfréd Haar created the first discrete wavelet transform (DWT) in Haar (1911) which led to the development of the fast Fourier transform (FFT) and the other forms of DWT to detect periodic signals. Haralick had developed grey level co-occurrence probabilities (GLCP) with statistical features to describe textures (Haralick et al., 1973). The hidden Markov model

The classification process starts with image acquisition to retrieve the partial information carried by the test image. The general framework of texture classification is divided into two important stages, which are feature extraction and feature selection. Feature extraction is a process of transforming a texture into a feature domain based on the intensity distribution on a digital image. On the other hand, feature selection is a process of integrating a set of conditional statements and routines, so that the computing system can logically decide

The texture classification methods are the tools which have been used to assist in medical imaging.

In food science, texture classification techniques aid on improving the quality of food and decrease the rate of food poisoning cases. For instance, a non-supervised method of texture classification is proposed to estimate the content of Intramuscular Fat (IMF) in beef on Du et

In remote sensing, texture classification is used to analyze satellite Synthetic Aperture Radar (SAR) imagery for multiple purposes. The analyzed satellite imagery might be used to monitor flood situations (Seiler et al., 2008), large scale construction planning, archeological research, weather forecasting, geological study etc. Texture segmentation is also applied in

(HMM) is another statistical method used for pattern recognition (Rabiner, 1989).

**1. Introduction**

which pixels belong to which texture.

al. (2008) and thus improving the meat quality.

**Methods in Mammography**

**Calcification Detection**

*Universiti Sains Malaysia,*

*Malaysia*

Hong Choon Ong and Hee Kooi Khoo

## **A Design for Stochastic Texture Classification Methods in Mammography Calcification Detection**

Hong Choon Ong and Hee Kooi Khoo *Universiti Sains Malaysia, Malaysia*

### **1. Introduction**

42 Human-Centric Machine Vision

Papakostas, G.A.; Boutalis, Y.S.; Karras, D.A. & Mertzios, B.G. (2009a). Pattern classification

Papakostas, G.A.; Tsougenis, E.D. & Koulouriotis, D.E. (2010a). Near optimum local image

*Imaging Systems and Techniques (IST'10),* pp. 459-462, Thessaloniki – Greece. Papakostas, G.A.; Karakasis, E.G. & Koulouriotis, D.E. (2010b). Novel moment invariants for

Rothe, I., Susse, H. & Voss, K. (1996). The method of normalization to determine invariants.

Samaria, F. & Harter, A.C. (1994). Parameterisation of a stochastic model for human face

Shen, D & Ip Horace H.S. (1999). Discriminative wavelet shape descriptors for recognition of

Sim, D.G.; Kim, H.K. & Park, R.H. (2004). Invariant texture retrieval using modified Zernike

Teague, M. (1980). Image analysis via the general theory of moments. *J. Optical Soc. Am.*, Vol.

Triesch, J. & von der Malsburg, C. (1996). Robust classification of hand postures against

Yap, P.T.; Paramesran, R. & Ong, S.H. (2003). Image analysis by Krawtchouk moments. *IEEE*

Zhu, H., Shu, H., Xia, T., Luo, L., Coatrieux, J.L. (2007). Translation and scale invariants of Tchebichef moments. *Pattern Recognition*, Vol. 40, No. 9, pp. 2530-2542.

complex backgrounds. Proceedings of the 2nd International Conference on Automatic Face and Gesture Recognition (pp. 170-175). Killington, Vermont, USA. Wee, C.Y. & Paramesran, R. (2007). On the computational aspects of Zernike moments.

*Computation*, Vol. 212, No. 1, pp. 162-176.

*Recognition*, Vol. 43, No. 1, pp. 58-68.

Computer Vision (pp. 138 – 142). Sarasota, FL, USA.

2-D patterns. *Pattern Recognition*, Vol. 32, No. 2, pp. 151-165.

*Image and Vision Computing*, Vol. 25, No. 6, pp. 967–980.

*Trans. on Image Processing*, Vol. 12, No. 11, pp. 1367-1377.

moments. *Image and Vision Computing*, Vol. 22, No. 4, pp. 331-342. Strang, G. & Nguyen, T. (1997). *Wavelets and Filter Banks*. Wellesley-Cambridge Press.

376.

70, pp. 920-930.

by using improved wavelet compressed Zernike moments. *Applied Mathematics and*

watermarking using Krawtchouk moments", *IEEE International Workshop on*

improved classification performance in computer vision applications. *Pattern*

*IEEE Transactions on Pattern Analysis and Machine Intelligence*, Vol. 18, No. 4, pp. 366-

identification. Proceedings of the 2nd IEEE Workshop on Applications of

A texture is a pattern represented on a surface or a structure of an object. Classifying a texture involves pattern recognition. Integration of texture classification is a one step ahead for machine vision application. For the past few decades until now, many researchers from various fields have been trying to develop algorithms to do texture discrimination. Alfréd Haar created the first discrete wavelet transform (DWT) in Haar (1911) which led to the development of the fast Fourier transform (FFT) and the other forms of DWT to detect periodic signals. Haralick had developed grey level co-occurrence probabilities (GLCP) with statistical features to describe textures (Haralick et al., 1973). The hidden Markov model (HMM) is another statistical method used for pattern recognition (Rabiner, 1989).

The classification process starts with image acquisition to retrieve the partial information carried by the test image. The general framework of texture classification is divided into two important stages, which are feature extraction and feature selection. Feature extraction is a process of transforming a texture into a feature domain based on the intensity distribution on a digital image. On the other hand, feature selection is a process of integrating a set of conditional statements and routines, so that the computing system can logically decide which pixels belong to which texture.

The texture classification methods are the tools which have been used to assist in medical imaging.

In food science, texture classification techniques aid on improving the quality of food and decrease the rate of food poisoning cases. For instance, a non-supervised method of texture classification is proposed to estimate the content of Intramuscular Fat (IMF) in beef on Du et al. (2008) and thus improving the meat quality.

In remote sensing, texture classification is used to analyze satellite Synthetic Aperture Radar (SAR) imagery for multiple purposes. The analyzed satellite imagery might be used to monitor flood situations (Seiler et al., 2008), large scale construction planning, archeological research, weather forecasting, geological study etc. Texture segmentation is also applied in

A Design for Stochastic Texture Classification Methods in Mammography Calcification Detection 45

not distinctive enough for high frequencies and thus could affect the segmentation result

Texture pattern can also be modeled as a transitional system using HMM. A Markov model assumes a texture pattern have a finite number of states and times. Each probability of a state is determined by the previous probability of the state. Three issues can arise from HMM observations; these are evaluation, decoding, and learning. HMM evaluation is to compare the probabilities of different models which best describe the texture feature. HMM decoding is to decompose and provide an estimated basis of texture patterns based on the HMM observations. The HMM learning searches for which model parameters best describe

The HMM can generate a consistent measurement for texture patterns based on the best probabilities. However, one or more observations which produce undesired probabilities could generate disorderly sequences due to the noisy pattern in the test image (Sonka et al.,

ACF is another method to describe texture patterns. The function helps to search for repeated patterns in a periodic signal. The function also identifies the missing basis of texture patterns hidden under the noisy patterns. In using the ACF, the mean of each image is adjusted before applying the general formula. Thus we are actually computing the normalized auto-covariance function. One can characterize a texture pattern by inferring the

The ACF feature is well demonstrated and distinctive between textures on a three dimensional (3D) graph. In feature selection, inferring the periodicity of a texture feature is done by observing several threshold points of the auto-covariance function and then counting the number of peaks for each threshold in a fixed variation. This may result in a random fluctuation and texture segmentation may fail because of two issues; there is not enough information by taking only one dimensional (1D) threshold to compare and an appropriate set of standard deviation of the distances between peaks are needed to know when the periodicity end. For example, the lagged product estimator and time series estimator are proposed to select ACF feature in Broerson (2005). But to appropriately

Instead of using the Gabor filter for feature extraction, wavelet is well known today as a flexible tool to analyze texture. Wavelet is a function whereby the basic function, namely the mother wavelet is being scaled and translated in order to span the spatial frequency domain. The DWT is done by the product of a corresponding signal generated by a pattern and the complex conjugate of the wavelet, and then integrating over all the distance points conveyed by the signal. In texture analysis, a tree data structure namely a packet wavelet expansion is used to split the signal into smaller packets and expand a chosen band at each

The DWT can become rotation invariant in texture analysis by taking the logarithm of the frequency sub-band. However this will damage the frequency components of the corresponding signal generated by a texture pattern and thus may obtain inaccurate result

GLCP method is a non-parametric solution whereby the textures are described in a discrete domain (Petrou & Sevilla, 2006). GLCP statistics are used to preserve the spatial characteristics of a texture. The selection of certain texture is possible as it is based on the

characterize the texture pattern by its periodicity is still an active area of research.

(Hammouda & Jernigan, 2000).

the texture pattern (Sonka et al., 2007).

periodicity of the pattern (Petrou & Sevilla, 2006).

level of resolution (Petrou & Sevilla, 2006).

(Khouzani & Zadeh, 2005).

2007).

the study of global warming by detecting and estimating the size of the icebergs at the North Pole from time to time (Clausi, 2002a).

Another breakthrough in the advancement of security is in biometric authentication. One of the biometric measures is in using the human eye as an identity of a person. The human iris contains texture patterns which are unique in nature and thus texture classification techniques can be applied on this iris imagery (Bachoo & Tapamo, 2005).

The first section has briefly introduced the texture classification and its applications. We highlighted the importance of texture classification in science. A literature review on the current existing methods is discussed in Section 2 and 3. The main challenge of texture classification is pointed out. Section 4 explains the theories behind our research methodologies. The non-parametric statistics are developed and used to generate statistical features from a given textured image. In the feature extraction section, we will provide a general formulation for our proposed method called cluster coding and explain briefly the parameters being used for texture classification purposes. Furthermore, the solution to solve the misclassification problem is also provided. A practical analysis on the digital mammogram using our proposed method is shown in Section 5. We will conclude with a summary of our contributions in Section 6.

There are four goals in this chapter which are as follows: -

### **2. Feature extraction**

There are several mathematical models which are used to extract texture features from an image. This section briefly discusses the five main feature extraction methods, namely autocorrelation function (ACF), Gabor filter, DWT, HMM, and GLCP. Besides, some combination of methods using different models are also developed to increase the robustness of texture segmentation system. We will briefly discuss two hybrid methods; these are Gabor wavelet and wavelet-domain HMM.

Gabor filter is one of the texture descriptors based on the Fourier transform (FT). A discrete FT is first taken from the test image to generate a two dimensional (2D) sinusoidal signal. For texture recognition, the Gabor filter bank contains a list of sub-bands of different signals generated by different textures. After the sub-bands are determined, the 2D signal of test image will multiply one of the chosen sub-bands and yield only the frequencies that match the sub-band. The product is then transformed back by taking the inverse FT and this leaves only the location of the texture feature which matches the signal. The process continues with each possible sub-band and produces the locations where the same signals occur (Petrou & Sevilla, 2006).

Gabor filter is useful in adapting sinusoidal signals whereby it can be decomposed into a weighted sum of sinusoidal signals. Thus Gabor filter is suitable to decompose textural information. Experiments done by Clausi & Deng (2005) stated that the Gabor filter can well recognize low and medium frequencies, but it produces inconsistent measurements for high frequencies due to the noise in the signal. The feature domain generated by Gabor filters are 44 Human-Centric Machine Vision

the study of global warming by detecting and estimating the size of the icebergs at the

Another breakthrough in the advancement of security is in biometric authentication. One of the biometric measures is in using the human eye as an identity of a person. The human iris contains texture patterns which are unique in nature and thus texture classification

The first section has briefly introduced the texture classification and its applications. We highlighted the importance of texture classification in science. A literature review on the current existing methods is discussed in Section 2 and 3. The main challenge of texture classification is pointed out. Section 4 explains the theories behind our research methodologies. The non-parametric statistics are developed and used to generate statistical features from a given textured image. In the feature extraction section, we will provide a general formulation for our proposed method called cluster coding and explain briefly the parameters being used for texture classification purposes. Furthermore, the solution to solve the misclassification problem is also provided. A practical analysis on the digital mammogram using our proposed method is shown in Section 5. We will conclude with a

a. To show how statistical texture descriptions can be used for image segmentation.

c. To develop algorithms for texture segmentation based on texture features analyzed. d. To illustrate the above objectives using an application on a mammogram image.

There are several mathematical models which are used to extract texture features from an image. This section briefly discusses the five main feature extraction methods, namely autocorrelation function (ACF), Gabor filter, DWT, HMM, and GLCP. Besides, some combination of methods using different models are also developed to increase the robustness of texture segmentation system. We will briefly discuss two hybrid methods;

Gabor filter is one of the texture descriptors based on the Fourier transform (FT). A discrete FT is first taken from the test image to generate a two dimensional (2D) sinusoidal signal. For texture recognition, the Gabor filter bank contains a list of sub-bands of different signals generated by different textures. After the sub-bands are determined, the 2D signal of test image will multiply one of the chosen sub-bands and yield only the frequencies that match the sub-band. The product is then transformed back by taking the inverse FT and this leaves only the location of the texture feature which matches the signal. The process continues with each possible sub-band and produces the locations where the same signals occur (Petrou &

Gabor filter is useful in adapting sinusoidal signals whereby it can be decomposed into a weighted sum of sinusoidal signals. Thus Gabor filter is suitable to decompose textural information. Experiments done by Clausi & Deng (2005) stated that the Gabor filter can well recognize low and medium frequencies, but it produces inconsistent measurements for high frequencies due to the noise in the signal. The feature domain generated by Gabor filters are

b. To develop algorithms for feature extraction based on texture images.

techniques can be applied on this iris imagery (Bachoo & Tapamo, 2005).

North Pole from time to time (Clausi, 2002a).

summary of our contributions in Section 6.

**2. Feature extraction**

Sevilla, 2006).

There are four goals in this chapter which are as follows: -

these are Gabor wavelet and wavelet-domain HMM.

not distinctive enough for high frequencies and thus could affect the segmentation result (Hammouda & Jernigan, 2000).

Texture pattern can also be modeled as a transitional system using HMM. A Markov model assumes a texture pattern have a finite number of states and times. Each probability of a state is determined by the previous probability of the state. Three issues can arise from HMM observations; these are evaluation, decoding, and learning. HMM evaluation is to compare the probabilities of different models which best describe the texture feature. HMM decoding is to decompose and provide an estimated basis of texture patterns based on the HMM observations. The HMM learning searches for which model parameters best describe the texture pattern (Sonka et al., 2007).

The HMM can generate a consistent measurement for texture patterns based on the best probabilities. However, one or more observations which produce undesired probabilities could generate disorderly sequences due to the noisy pattern in the test image (Sonka et al., 2007).

ACF is another method to describe texture patterns. The function helps to search for repeated patterns in a periodic signal. The function also identifies the missing basis of texture patterns hidden under the noisy patterns. In using the ACF, the mean of each image is adjusted before applying the general formula. Thus we are actually computing the normalized auto-covariance function. One can characterize a texture pattern by inferring the periodicity of the pattern (Petrou & Sevilla, 2006).

The ACF feature is well demonstrated and distinctive between textures on a three dimensional (3D) graph. In feature selection, inferring the periodicity of a texture feature is done by observing several threshold points of the auto-covariance function and then counting the number of peaks for each threshold in a fixed variation. This may result in a random fluctuation and texture segmentation may fail because of two issues; there is not enough information by taking only one dimensional (1D) threshold to compare and an appropriate set of standard deviation of the distances between peaks are needed to know when the periodicity end. For example, the lagged product estimator and time series estimator are proposed to select ACF feature in Broerson (2005). But to appropriately characterize the texture pattern by its periodicity is still an active area of research.

Instead of using the Gabor filter for feature extraction, wavelet is well known today as a flexible tool to analyze texture. Wavelet is a function whereby the basic function, namely the mother wavelet is being scaled and translated in order to span the spatial frequency domain. The DWT is done by the product of a corresponding signal generated by a pattern and the complex conjugate of the wavelet, and then integrating over all the distance points conveyed by the signal. In texture analysis, a tree data structure namely a packet wavelet expansion is used to split the signal into smaller packets and expand a chosen band at each level of resolution (Petrou & Sevilla, 2006).

The DWT can become rotation invariant in texture analysis by taking the logarithm of the frequency sub-band. However this will damage the frequency components of the corresponding signal generated by a texture pattern and thus may obtain inaccurate result (Khouzani & Zadeh, 2005).

GLCP method is a non-parametric solution whereby the textures are described in a discrete domain (Petrou & Sevilla, 2006). GLCP statistics are used to preserve the spatial characteristics of a texture. The selection of certain texture is possible as it is based on the

A Design for Stochastic Texture Classification Methods in Mammography Calcification Detection 47

probabilities are only dependent upon the statistical parameters, so we must carefully choose

The self-organizing map (SOM) is an unsupervised single layer artificial neural network (ANN). In the SOM training environment, a digital image is mapped as a grid. A set of neurons will be placed at random grid points where each neuron is stored as a cluster centre (Chen & Wang, 2005). SOM clusters regions which have similar pattern and separates the dissimilar patterns based on a general distance function (Martens et al., 2008). The SOM learning process is similar to the K-means clustering where it iterates until each of the

The advantage of SOM is that the more number of neurons are placed in the grid, the higher classification result will be obtained. However, if the numbers of neurons are too large, SOM may end up with over classification (Martens et al., 2008). On the other hand, the numbers of neurons required is unknown. Furthermore, the classification using SOM may fail at the

In a stochastic texture, it is not possible to determine (a), (b), and (c). To classify this texture, a measure of its density has to be computed. Therefore, the statistical modelling of GLCP is used to extract textural features. A general classification approach, namely the cluster coding is

Given a textured image, the grey level intensities of the pixel set can be retrieved. A grid with equal sized cells is fitted onto the image. For each cell, the grey level intensities of the local neighbourhood pixels are taken into consideration when computing the joint probabilities density function. The Equation 1 shows the statistical feature for a cell, *v*(*i*,*j*) based on the joint

> 1 , 0 1

(,) (, )

∑ ∑

*Fst vi j*

−

*G s t G*

, 0

=

⎝ ⎠

, {1,2,..., 1}

*f k cdf cdf f k G* (2)

*s t*

= −

where *F*(*s*,*t*) represents the frequency of occurrence between grey levels, *s* and *t*. *G* is the

Histogram equalization is applied to achieve higher separability between to different textures in the image. The feature set, *v* is distributed in a frequency histogram in order to

> ⎧ = = ⎨ +∈ − ⎩ 1 ,0

⎛ ⎞ ⎜ ⎟ = ⎜ ⎟

(,)

(1)

*Fst*

proposed. This approach preserves the textural information of a given textured image.

probability density function. *i* and *j* are the positions of a cell on the image.

quantized grey levels which forms a *G×G* co-occurrence matrix.

*k*

calculate the cumulative distribution function, *cdfk* using Equation 2,

−

*k k*

*k*

these parameters, especially the parameters for texture patterns (Diplaros et al., 2007).

cluster centre converges to the centre of the possible texture patterns.

local minimum in training (Abe, 2005).

A texure pattern can be identified by

a. Counting the number of repetition of its primitive,

**4.1 Enhanced two-dimensional statistical features**

b. Average brightness of the primitive, and c. The structural distribution of its primitive.

**4. Methodology**

statistical features. The best statistical features that are used for analysis are entropy, contrast, and correlation (Clausi, 2002a). However, further analysis in Jobanputra & Clausi (2006) shows that correlation is not suitable for texture segmentation. GLCP statistics can also be used to discriminate between two different textures. Boundaries can be created from the shift on statistical feature while moving from one texture to another (Jobanputra & Clausi, 2006). Clausi & Zhao (2002) also proposed grey level co-occurrence linked-list (GLCLL) structure and grey level co-occurrence hybrid histogram (GLCHH) structure in Clausi & Zhao (2003) to this non-parametric solution for storing purpose in order to speed up the computational time for GLCP feature extraction.

Wavelet-domain HMM can be described as a finite state machine in the wavelet domain. The hidden Markov tree (HMT) can identify the characteristics of the joint probabilities of DWT by capturing the scale dependencies of wavelet co-efficient via Markov chains (Fan & Xia, 2003). Since this method is based on the wavelet domain, the disorder sequence of signal generated by the noise in the texture patterns may weaken the cross-correlation between DWT sub-bands (Fan & Xia, 2003). Ming et al. (2008) proposed the wavelet hiddenclass-label Markov random field to suppress the specks of noise, but there are some small blobs of noises still appearing in the results of the segmented imagery.

### **3. Feature selection**

The texture feature selection stage is the most important part of the texture segmentation process because it determines which pixels belong to which texture of an image. The use of parameters of a chosen method is crucial to the output of the result.

K-means algorithm has become a popular clustering method which is used for pattern recognition. Given the number of clusters, K, the algorithm will start at a random K number of centres. Then, each of the centres will group the features using the closest distances or Euclidean distance measures. The locations of the features with the same cluster will determine the new centre for each cluster. The process will then repeat until the centre of each texture class remains the same.

K-means algorithm assumes that all clusters are in spherical shape, but it may return inappropriate result for non-spherical clusters (Clausi, 2002b). In the real life 2D feature set is not always in spherical shape and normally the number of classes is unknown.

Support vector machines (SVM) algorithm is a slow but highly accurate clustering method. The SVM training algorithm was introduced by Boser et al. (1992). The purpose of SVM is to map feature vectors into a higher dimensional feature space, and then creating a separating hyperplane with maximum margin to group the features. Support vectors (SVs) contain highlighted pixels that help to create the margins or boundaries in an image. The higher dimensional space is defined by a kernel function. Some of the popular kernels are shown in Schölkopf & Smola (2002). A combined GLCP and SVM technique has been proposed for two class segmentation with significant result in Khoo et al. (2008).

Expectation-maximization (EM) algorithm is a statistical estimation algorithm used for finding maximum likelihood and estimates of parameters in probabilistic models. The parameters involve means, variances, and weights (Tong et al., 2008). The EM algorithm starts by initializing the parameters to compute the joint probability for each cluster. The algorithm then iterates to re-estimate the parameters and maximize the probability for each cluster until the set of convergence values of probabilities are obtained (Bishop, 2006). The convergence probabilities are only dependent upon the statistical parameters, so we must carefully choose these parameters, especially the parameters for texture patterns (Diplaros et al., 2007).

The self-organizing map (SOM) is an unsupervised single layer artificial neural network (ANN). In the SOM training environment, a digital image is mapped as a grid. A set of neurons will be placed at random grid points where each neuron is stored as a cluster centre (Chen & Wang, 2005). SOM clusters regions which have similar pattern and separates the dissimilar patterns based on a general distance function (Martens et al., 2008). The SOM learning process is similar to the K-means clustering where it iterates until each of the cluster centre converges to the centre of the possible texture patterns.

The advantage of SOM is that the more number of neurons are placed in the grid, the higher classification result will be obtained. However, if the numbers of neurons are too large, SOM may end up with over classification (Martens et al., 2008). On the other hand, the numbers of neurons required is unknown. Furthermore, the classification using SOM may fail at the local minimum in training (Abe, 2005).

### **4. Methodology**

46 Human-Centric Machine Vision

statistical features. The best statistical features that are used for analysis are entropy, contrast, and correlation (Clausi, 2002a). However, further analysis in Jobanputra & Clausi (2006) shows that correlation is not suitable for texture segmentation. GLCP statistics can also be used to discriminate between two different textures. Boundaries can be created from the shift on statistical feature while moving from one texture to another (Jobanputra & Clausi, 2006). Clausi & Zhao (2002) also proposed grey level co-occurrence linked-list (GLCLL) structure and grey level co-occurrence hybrid histogram (GLCHH) structure in Clausi & Zhao (2003) to this non-parametric solution for storing purpose in order to speed

Wavelet-domain HMM can be described as a finite state machine in the wavelet domain. The hidden Markov tree (HMT) can identify the characteristics of the joint probabilities of DWT by capturing the scale dependencies of wavelet co-efficient via Markov chains (Fan & Xia, 2003). Since this method is based on the wavelet domain, the disorder sequence of signal generated by the noise in the texture patterns may weaken the cross-correlation between DWT sub-bands (Fan & Xia, 2003). Ming et al. (2008) proposed the wavelet hiddenclass-label Markov random field to suppress the specks of noise, but there are some small

The texture feature selection stage is the most important part of the texture segmentation process because it determines which pixels belong to which texture of an image. The use of

K-means algorithm has become a popular clustering method which is used for pattern recognition. Given the number of clusters, K, the algorithm will start at a random K number of centres. Then, each of the centres will group the features using the closest distances or Euclidean distance measures. The locations of the features with the same cluster will determine the new centre for each cluster. The process will then repeat until the centre of

K-means algorithm assumes that all clusters are in spherical shape, but it may return inappropriate result for non-spherical clusters (Clausi, 2002b). In the real life 2D feature set

Support vector machines (SVM) algorithm is a slow but highly accurate clustering method. The SVM training algorithm was introduced by Boser et al. (1992). The purpose of SVM is to map feature vectors into a higher dimensional feature space, and then creating a separating hyperplane with maximum margin to group the features. Support vectors (SVs) contain highlighted pixels that help to create the margins or boundaries in an image. The higher dimensional space is defined by a kernel function. Some of the popular kernels are shown in Schölkopf & Smola (2002). A combined GLCP and SVM technique has been proposed for

Expectation-maximization (EM) algorithm is a statistical estimation algorithm used for finding maximum likelihood and estimates of parameters in probabilistic models. The parameters involve means, variances, and weights (Tong et al., 2008). The EM algorithm starts by initializing the parameters to compute the joint probability for each cluster. The algorithm then iterates to re-estimate the parameters and maximize the probability for each cluster until the set of convergence values of probabilities are obtained (Bishop, 2006). The convergence

is not always in spherical shape and normally the number of classes is unknown.

up the computational time for GLCP feature extraction.

**3. Feature selection**

each texture class remains the same.

blobs of noises still appearing in the results of the segmented imagery.

parameters of a chosen method is crucial to the output of the result.

two class segmentation with significant result in Khoo et al. (2008).

A texure pattern can be identified by

In a stochastic texture, it is not possible to determine (a), (b), and (c). To classify this texture, a measure of its density has to be computed. Therefore, the statistical modelling of GLCP is used to extract textural features. A general classification approach, namely the cluster coding is proposed. This approach preserves the textural information of a given textured image.

#### **4.1 Enhanced two-dimensional statistical features**

Given a textured image, the grey level intensities of the pixel set can be retrieved. A grid with equal sized cells is fitted onto the image. For each cell, the grey level intensities of the local neighbourhood pixels are taken into consideration when computing the joint probabilities density function. The Equation 1 shows the statistical feature for a cell, *v*(*i*,*j*) based on the joint probability density function. *i* and *j* are the positions of a cell on the image.

$$w(i,j) = \sum\_{s,t=0}^{G-1} \left( \frac{F(s,t)}{\sum\_{s,t=0}^{G-1} F(s,t)} \right) \tag{1}$$

where *F*(*s*,*t*) represents the frequency of occurrence between grey levels, *s* and *t*. *G* is the quantized grey levels which forms a *G×G* co-occurrence matrix.

Histogram equalization is applied to achieve higher separability between to different textures in the image. The feature set, *v* is distributed in a frequency histogram in order to calculate the cumulative distribution function, *cdfk* using Equation 2,

$$cdf\_k = \begin{cases} f\_k & k=0\\ cdf\_{k-1} + f\_{k'} & k \in \{1, 2, \dots, G-1\} \end{cases} \tag{2}$$

A Design for Stochastic Texture Classification Methods in Mammography Calcification Detection 49

Figure 2 demonstrates the basic concept of the cluster coding algorithm. The classification process is run on a split-and-merge operation. The splitting part is shown in Figure 2(b) and

(a) (b) (c) (d) Fig. 2. Cluster coding conceptual diagram. (a) Desired segments. (b) Feature A coded in 2 clusters. (c) Feature B coded in 2 clusters. (c) Summation of codes from (b) and (c).

Each feature set is split up into coded clusters using the mean shift estimator (MSE). Mean

where *K* is the smoothing kernel function. *δ*1 and *δ*2 are the shifting parameters of *v*(*i*,*j*). The mode seeking of one-dimensional version of feature set is illustrated in Figure 3. In the figure, the dotted line is the splitting point at the position 3.5. The flowchart for cluster

⎡ ⎤ −+ + ⎣ ⎦ = ⎡ ⎤ −+ + ⎣ ⎦

1 2 1 2

δ

 δ

> δ

(4)

δ

shift estimator locate the local maxima of a feature set. The MSE is defined as,

Fig. 3. Line graph of the feature values and the corrresponding MSE.

**4.2 Cluster coding classification**

coding is shown in Figure 4.

2(c), while the merging part is shown in Figure 2(d).

The general definition for histogram equalization (Gonzalez & Woods, 2006), *Hk* is

$$H\_k = \frac{cdf\_k - cdf\_{\min}}{t - cdf\_{\min}} \times (G - 1), \ k \in \{1, 2, \dots, G - 1\} \tag{3}$$

where *cdfk*, a cumulative distribution function is the total number of feature sets and min is the minimum for the set of cumulative distribution function. The design process is shown in Figure 1.

Fig. 1. The flowchart of the GLCP algorithm for feature extraction.

### **4.2 Cluster coding classification**

48 Human-Centric Machine Vision

( 1), 1,2,..., 1 *k*

− = ×− ∈ − − (3)

The general definition for histogram equalization (Gonzalez & Woods, 2006), *Hk* is

*cdf cdf H Gk G*

where *cdfk*, a cumulative distribution function is the total number of feature sets and min is the minimum for the set of cumulative distribution function. The design process is shown in

Start

1.Compute a grey level intensity, *Ii* for each pixel on a test

3.Compute the corresponding grey level quantized value for

2.Set a grey level quantization, *G* value*.*

5.Set the distance between pixel pair (fixed to 1).

If cell reached the pixel block at the bottom right corner of the test image

> Initialize a GLCM(*G×G*) for each angle used.

NO

YES End

1.For each pixel, get the values of quantized grey level at the current

2.Assign the values of the pixel pair in step 1 as entry position to the GLCM, then increase the position in the matrix by one to become the

> Based on the frequencies stored in GLCM, 1.Compute the joint probability density

> > Reset all entries of the GLCMs and the GLCP.

2.Compute all the selected statistical features based on GLCP.

frequency of occurrence between the pixel pair, *F*(*s*,*t*). 3.Repeat step 1 & 2 for other angles used in another GLCM.

function or GLCP.

In the local window/ pixel block:

and the neighbour position.

Fig. 1. The flowchart of the GLCP algorithm for feature extraction.

Move to the next cell (Cell sequence started from top to bottom, then follow by left to right)

First cell start from the pixel block at the upper left corner of the test image

 { } min min

*t cdf*

image.

each *Ii*.

4.Set the cell size, *M×N.*

6.Set the angles to be used.

*k*

Figure 1.

Figure 2 demonstrates the basic concept of the cluster coding algorithm. The classification process is run on a split-and-merge operation. The splitting part is shown in Figure 2(b) and 2(c), while the merging part is shown in Figure 2(d).

Fig. 2. Cluster coding conceptual diagram. (a) Desired segments. (b) Feature A coded in 2 clusters. (c) Feature B coded in 2 clusters. (c) Summation of codes from (b) and (c).

Each feature set is split up into coded clusters using the mean shift estimator (MSE). Mean shift estimator locate the local maxima of a feature set. The MSE is defined as,

$$MSE(v(i,j)) = \frac{K\left[\lfloor v(i,j) - v(i + \delta\_{1}, j + \delta\_{2} \mid )\right] v(i,j)}{K\left[\lfloor v(i,j) - v(i + \delta\_{1}, j + \delta\_{2}) \mid \right]} \tag{4}$$

where *K* is the smoothing kernel function. *δ*1 and *δ*2 are the shifting parameters of *v*(*i*,*j*). The mode seeking of one-dimensional version of feature set is illustrated in Figure 3. In the figure, the dotted line is the splitting point at the position 3.5. The flowchart for cluster coding is shown in Figure 4.

Fig. 3. Line graph of the feature values and the corrresponding MSE.

A Design for Stochastic Texture Classification Methods in Mammography Calcification Detection 51

The region reclassifies if there exist a dominated population surrounding its neighbourhood

NO

YES

1. Search all the same and connected pixels

Check if both the Rule I and Rule II are true.

All the interior pixels are merge to the label of dominant edges.

YES

1. Count the size for each label of the edges

2. Determine the label of dominant edge,

and store in the *interior* vector. 2. Store any counter label of edges in the

*edge* vector.

End

NO

from the *edge* vector.

*edgedmt* among all the edges.

Rule I: The located region has achieved the minimum population, and

region. Figure 5 provides the detail flowchart for this approach.

Start

1.The initial pixel start at the top left of

2. Set the first pixel as centre pixel .

the image.

Check if the centre pixel already visited

Fig. 5. The flowchart of the modified flood fill algorithm for cluster merged.