Citation

Material Information

Title:
Wavelet-based image segmentation
Creator:
Zarrini, Haleh
Place of Publication:
Denver, Colo.
Publisher:
Publication Date:
Language:
English
Physical Description:
xi, 75 leaves : illustrations ; 28 cm

Thesis/Dissertation Information

Degree:
Master's ( Master of Science)
Degree Grantor:
Degree Divisions:
Department of Electrical Engineering, CU Denver
Degree Disciplines:
Electrical Engineering
Committee Chair:
Bialasiewicz, Jan T.
Committee Members:
Papantoni, Titsa
Deng, Yiming

Subjects

Subjects / Keywords:
Wavelets (Mathematics) ( lcsh )
Image processing -- Digital techniques ( lcsh )
Digital images ( lcsh )
Electronic noise ( lcsh )
Digital images ( fast )
Electronic noise ( fast )
Image processing -- Digital techniques ( fast )
Wavelets (Mathematics) ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 73-75).
General Note:
Department of Electrical Engineering
Statement of Responsibility:
by Haleh Zarrini.

Record Information

Source Institution:
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
786172486 ( OCLC )
ocn786172486
Classification:
LD1193.E54 2011M Z37 ( lcc )

Full Text
WAVELET BASED IMAGE SEGMENTATION
by
Haleh Zarrini
B.S., University of Urmia, 1991
A thesis submitted to the
In partial fulfillment
Of the requirements for the degree of
Master of Science
Electrical Engineering
2011

This thesis for the Master of Science
degree by
Haleh Zarrini
has been approved
by
Jan T. Bialasiewicz
TitsaPapantoni
///o Date

Zarrini, Haleh (M.S., Electrical Engineering)
Wavelet Based Image Segmentation
Thesis directed by Professor Jan T. Bialasiewicz
ABSTRACT
Wavelet analysis and its applications have become one of the fastest growing
research areas in recent years. It has been employed in many fields like
Signal and Image Processing, Communications Systems, Biomedical
Imaging, and other areas.
The main contribution of this thesis is the development of a new algorithm,
which performs image segmentation and simultaneously filters out the noise
present in the original image. This algorithm, similarly as the well known
Sobel algorithm, is based on extraction or detection of the edges of an image
to be segmented. However, to detect edges the wavelet transform analysis is
performed using a trous algorithm, which assigns two numbers to every
image pixel at every decomposition level: a horizontal and a vertical detail to
be alternatively represented by a modulus and an angle. Such representation
facilitates obtaining the image edges. It is well known that most of the noise
contained in the original image is extracted in the details of the first level of
the wavelet decomposition. Therefore, the details at the second
decomposition level are almost noise free. The modulus maxima
representation of these details determines the noise-free representation of the
edges of the original image, which is provided as an input image to K-means
algorithm. Summarizing, one can say that the new segmentation algorithm,

developed and investigated in this thesis, consists of the image edge
detection using wavelets modulus maxima of the original image and of the
K-means image segmentation algorithm.
In addition, the thesis gives a comprehensive introduction to wavelets and
their applications to signal/image processing. Moreover, the problem of image
segmentation is introduced and the overview of known segmentation methods
and segmentation applications is given.
This abstract accurately presents the content of the Candidates thesis. I
recommend its publication.
Signed
Jan T. Bialaseiwicz

DEDICATION
I dedicate this thesis to my father Muhammad Arman and to my husband,
Hossein for their continued support and understanding, which gave me an
opportunity to pursue higher education. Special thanks to my two boys, Daniel
and AN, that helped me during the course of my education. My deepest
gratitude goes to them for believing that I could make it to graduate school.

ACKNOWLEDGEMENT
My gratitude goes to my advisor, Dr. Jan Bialasiewicz for introducing me to
the world of wavelets, and for his guidance throughout my research. Through
his guidance, I have taken many interesting signal processing classes which
have greatly helped in understanding wavelet applications. I also wish to
thank Dr. Titsa Papantoni and all the members of my committee for their
participation and support.
i

Figures.................................................................ix
Tables..................................................................xi
Chapter
1. Introduction..........................................................1
1.1 Objectives...........................................................1
1.2 Outlines.............................................................2
2. Brief Review on Wavelet...............................................4
2.1 Wavelet Transform....................................................7
2.2 Continuous Wavelet Transform (CWT)..................................8
2.3 Discrete Wavelet Transform (DWT)....................................9
2.4 Multiresolution Analysis (MRA).....................................10
2.5 Wavelet Based on the Gaussian Function.............................17
2.6 Wavelet Transform Modulus Maxima...................................19
2.7 Algorithme a Trous...............................................21
3. Clustering preview...................................................26
3.1 K-means Algorithms..................................................28
3.2 Image Clustering...................................................35
vii

4. Segmentation and its Method
41
5. Conclusion.......................................................53
Appendix
A. K-means Clustering Algorithm MATLAB Program......................55
B. K-means Clustering Algorithm according to the color property.....57
C. K-means Clustering Algorithm on the color image..................58
D. Algorithme a Trous MATLAB Program................................60
E. MATLAB Program for Segmentation..................................64
F. MATLAB Program for Segmentation by Trasholding...................70
G. MATLAB Program for Segmentation by Sobel gradint.................71
References...........................................................73
viii

LIST OF FIGURES
Figure
2.1 Time-frequency boxes representing the energy spread of Gabor toms ... 6
2.2 Morelet wavelet function and Mexican-hat Wavelet...................7
2.3 The Relationship between scaling and wavelet function spaces......11
2.4 Harr scaling and wavelet Function.................................13
2.5 The Analysis Filter Bank for Two-Dimensional Wavelet Transform...15
2.6 One (left) and Two (right) Levels Two-Dimensional DWT.............16
2.7 The scaling function..............................................18
2.9 Filter bank for two dimensional image decomposition implementing
the algorithme a trous............................................24
2.10 Filter bank implementing image reconstruction with the
algorithme a trous...............................................25
3.1 Illustrate clustering performed based on color property........... 27
3.2 Flowchart of the k-means algorithm................................30
3.3 The K-means algorithm on a gray-scale image...................... 37
3.4 Color expression by L*a*b concept.................................38
3.5 The K-means algorithm on color image..............................40
4.1 Segementation of original imagr beased on Trasholding..............43
4.2 Sobel operators...................................................45
4.3 Segmentation of an image based on sobel operators.................47
4.4 3-D mesh polt for Lena.bmp......................................49
4.5 Modulus Maxima at 3 level.........................................49
4.6 Block Diagram of image segmentation...............................50
IX

4.7 Original image and its 3-D mesh plot...............................51
4.8 Edge detection with Modulus Maxima.........................'....... 51
4.9 Edge Detection with all holes filled...............................52
4.10 Segmentation objects..............................................52
X

LIST OF TABLES
Tables
2.1 Spline filter coefficient used in the algorithm a trous..............23
2.2 Example data of K-means algorithm....................................31
XI

1. Introduction
The subject of image processing has always been expanding and
growing. The term image processing is used for all manipulation on an
image if the output of the manipulation is again an image. The field of
digital image processing refers to processing digital images by a digital
computer. A digital image is composed of a finite number of elements,
each of which has a particular location and intensity value. They are
referred to as picture elements or pixels.
Image segmentation refers to the process of partitioning a digital image
into multiple segments and can be used for object recognition, image
editing, or image database search. The goal of image segmentation is to
cluster pixels into different image regions that correspond to individual
surfaces and objects. The image segmentation algorithm is generally
based on one of two basic properties of intensity values that are
discontinuity and similarity. In the first category, the approach is to
partition an image based on abrupt change in intensity, such as image
edges. In the second category, partitioning an image into regions is based
on similar properties.
1.1 Objective
There are many approaches and techniques for image segmentation. The
object of this thesis is image segmentation using wavelets, which is a new
approach that provides faster and more accurate results with less
calculation.
l

Wavelets are purposefully crafted to have specific properties that make
them useful for signal processing. The wavelet transform has many
properties (multi resolution representation, orthogonality and fast
algorithm). They can be used for a wide variety of fundamental signal
processing tasks, such as compressions, removing noise and enhancing
an image.
In this thesis, I present a different approach to segmenting of an image.
An algorithm was developed to achieve this goal. This algorithm first
applies Wavelet Transform Modus Maxima (WTMM), which used a trous
algorithm to generate the approximation, horizontal and vertical
coefficients of an image. Then the corresponding modulus maxima and
their angular directions are found. The result displayed the edges of the
image. Then the K-means algorithm was utilized to obtain segments and
the object was separated from the background.
1.2 Thesis outline
This thesis is outlined as follows:
Chapter 2 reviews Wavelet Transform; it introduces the Continuous
Wavelet Transform (CWT), District Wavelet Transform (DWT), Multi
resolution Analysis (MRA), Wavelet Transform Modus Maxima (WTMM),
and a trous algorithm.
Chapter 3 introduces clustering, and its applications. One of the
methods that could be used for clustering is K-means algorithm. The
2

K-means algorithm is a numerical, unsupervised, non-deterministic,
iterative method. It is simple and very fast, so in many practical
applications, this method has proven to be very suitable for producing
globular clusters. In this chapter, the K-means algorithm is applied on
gray-scale and on the color images.
Chapter 4 presents the segmentation and its realization by different
methods. Image segmentation is an essential step in many advanced
imaging applications. Accurate segmentation is required for medical
diagnosis, image guided procedures, and 3D rendering. This chapter
describes several algorithms and gives their comparison.
Chapter 5 provides the results and the conclusions based on our
algorithm.
3

2. Brief Review on Wavelets
A wave is defined as an oscillating function of time or space. In terms of
signals or functions are sinusoids or complex exponentials, which has proven
to be extremely valuable in mathematics, science and engineering. A wavelet
is a small wave with energy concentrated in time to give a tool for the analysis
of the transient, non-stationary or time varying phenomena.
Wavelets are functions, which decompose the signal to components of waves
or frequencies and then analyze each component with a resolution to be
matched to its scale. The goal of wavelet research is to create a set of either
general exponential functions or basis functions and transform that will give
useful, efficient and informative description of a signal.
The wavelet transform is capable of providing the signal contained as a
function of time and frequency and can be considered as a continuation of
Fourier transform. Fourier transform provides a way to convert samples of
time-series into the frequency domain. The main disadvantage of the Fourier
transform is that it only has frequency resolution and no time resolution. Even
though we may be able to determine all the frequencies present in a signal,
we dont know where and when they are present. To overcome this, the
windowed Fourier transform replaces the sinusoidal wave of Fourier
Transform by sinusoid and a window which is localized in time with frequency
resolution problems. The windowed Fourier transforms defined by the
physicist Dennis Gabor in 1964; shows elementary time frequency atoms as
waveforms that have a minimal spread in a time-frequency plane. To
measure time-frequency information content, he suggested decomposing
signal over these elementary wavelet.
4

The windowed Fourier transform correlates a signal / with each atom gu,( as
represented by the following equation:
Sf(u,0 = J_+/(t)g^(Odt = /_+/(t)g(t u)e"iftdt (2.1)
where u represents the time instance and Â£ represents frequency. The
multiplication of g(t u) localized the Fourier integral in the time window of a
neighborhood around u The resolution of time and frequency of the
Windowed Fourier transform depends on the spread of the window in time
and frequency.
Figure 2.1 shows the time-frequency boxes of g,$representing the energy spread of Gabor atom by Heisenberg rectangle. This rectangle is centered at (u, 0 and has a time width (St) and frequency width (Sco) as St .8to> K The issue of the windowed Fourier transform was that low frequency can hardly be depicted by short windows and high frequencies can only be poorly contained by long windows. In wavelet analysis a completely scalable modulated window is shifted along the signal. Then this process is repeated many times with a slightly shorter or longer window for every new cycle. This is the important characteristic of wavelets. The result is a collection of a time- frequency representations of the signal, all with different resolutions or frequency bands [10, 13]. 5 Figure 2.1 Time-frequency boxes (Heisenberg rectangles) representing the energy spread of Gabor atoms. With wavelets, we usually dont speak about time-frequency but of time-scale representations, where small scales means correspond wavelet and thus high frequencies. Large scales correspond to dilated wavelets and therefore lower frequencies. However, if we are interested in non-stationary signals, the wavelet transform or wavelet analysis is probably the most recent solution to overcome the shortcoming of the Fourier transforms. 6 2.1 Wavelet Transforms A wavelet transform decompose signals over dilated and translated wavelets of varies frequency components. Suppose that xp(t) is a wavelet function i/> G L2(R), with a zero average: 0(0* = 0 (2.2) It is normalized ||i//(t)|| = L and centered in the neighborhood of t=0. At a scale of s or dilation variable and translation of u the function xpUiS(t) is defined by: tMO = ^ (t1) <2-3> where is the normalizing factor to ensure that all wavelets have the same energy [10]. Figure 2.2 shows an example of two different mother wavelets. Mortet wavelet Mexican hat wavelet Figure 2.2 Morelet wavelet function and Mexican-hat Wavelet. 7 2.2 Continuous Wavelet Transforms Continuous Wavelet Transforms (CWT), transforms a continuous-time signal into its wavelet transform that is a function of continuous shift or translation and a continuous scale. The continuous wavelet transform of a continuous, square integrable function /(t) at a scale s > 0 and the time shift u Â£ R is given by the following formula: Wf(u,s) = Ll/(0^i/>* (^-r)dt = (f>^Pu,s) (2-4) where xp(t) is a continuous function in both time and frequency domain, called the Mother wavelet. For any function / e L2(R), the function/(t) can be obtained by taking the inverse continuous wavelet transform under admissibility condition m = i-C -C duds (2.5) where c, = < + (2 -6) Is called the wavelet admissibility, where r/5(o;) is the Fourier transform of xp(t) [6, 10]. This condition insures that ^(0) = 0, that is the wavelet has no zero frequency components. A wavelet should also have finite energy, that is, /." hKOPdtcoo (2.7) 8 2.3 Discrete Wavelet Transforms A discrete wavelet transform (DWT) is any wavelet transform for which the wavelets are discretely sampled and it is a very important concept for data and image compression. The first DWT was invented by the Hungarian mathematician Alfred Haar. The continuous wavelet transform has disadvantages of redundancy and not being practical with digital computers. Consequently, the scale 5 and shift u parameters are calculated on a discrete grid of time-scale plane which leads to a discrete set of continuous basis function [13]. Our goal is to generate a set of function such that any signal in / e L2(R) can be defined by: /(0 = Sy.fc cij k s~]/2 x/;( s~Jt ku0) = ajk xpjk{t) (2-8) with = s0J/2 xp(s0Jt-ku0) (2.9) where, s = s and u = /cSqU0 for j,k el (s0 >1 is a dilated step and u0 0 is a translation step). The two dimensional set of coefficients ajik is called the discrete wavelet transform of fit). In addition, when an input image is a discrete signal made up of pixels defined within a certain resolution, so the wavelet transform cant be computed at an arbitrary fine scale. To construct a discrete wavelet from a continuous wavelet, the scale 5 is discretized but not the translation parameter u. For fast computations the scale is sampled along a dyadic 9 sequence {2;}, to simplify the numerical calculations [10]. The dyadic wavelet transform of a signal / e L2(R) is defined by: 2.4 Multiresolution Analysis (MRA) Multiresolution or multiscale analysis is the important concept in wavelet analysis. It decomposes the signal into sub-signals of different size and resolution levels. Therefore a signal can be implemented in a multiple resolutions. The different signal scales represent distinct resolution in the same signal. That is why, multi resolution analysis is a better way to perform discrete wavelet transform. The multiresolution analysis is defined as a sequence {Vj: j e Z} of closed subspaces and consists of a sequence of embedded subspaces V_2 c V_i c v0 c c V2 of L2(R). Figure 2.3 illustrated the nesting spaces of L2(R). If the Wj is the difference between any two adjacent scaling subspace Vj and Vj+1, then the Vj+i is obtained by composing Vj and Wj. where denotes the union of spaces and Wj is orthogonal component of Vj and Y/+i- The Vj and Wj are orthogonal to each other. (2.10) VJ+1 = Vj 0 Wj (2.11) 10 Figure 2.3 The relationship between scaling and wavelet function spaces. V0 1 W0 1 Wi The multiresolution follows the following condition: 1. Vj c vj+1, V jEZ 2. V-co = {0} V00 = L2 3. HO e Vj /(2t) G V;+1 4. Va = V0 + ^o V2 = V0 + W0 + W1 5. L2 = V0 + Wi + W2 + - 6. V0 = W.0o + + W_2 + A scaling function { translations and introduced such that = 2 7/2 is an orthonormal basis of the subspace Vj. li If Wj is orthornormal component of Vj (Wj 1 Vj) in subspace then the wavelet function ip(t) introduced such that for each j iPi.M = 2-;/2 ip (7") ,(j,k2) (2.13) is an orthonormal basis of subspace Wj. The scaling function (p(t) = y/2'Znh 2 where, /i^(n) = ((p00,(p-lin) is a scaling wavelet coefficient. Y,n\tl(p(n)\ = 1 for orthogonality. Applying the same analogy for xp(t) we have i/>(t) = yl2Y,nh^,(n) cp(2t n) (2.15) where, h^in) = (xpoiQ, If the wavelet spans the orthogonal complement spaces and the integer translation of the wavelets is orthogonal, then the relation between the wavelet function coefficient and the scaling function coefficient is also as: hxpin) = (l)nAiv(l n) (2.16) where hv(ri), h^(n) can be viewed as low pass (LP) and high pass (HP) filter coefficient respectively. The equation (2.15) is called multiresolution analysis; it means that any expansion function of any subspace can be double- resolution copies of itself. That is the component or function in subspace at lower resolution can be presentable by just a shift at the next higher resolution. 12 Consider the Haar scaling function y(x) , r1 0 < x < 1 *W = lo otAmW*, (2'17> whose filter coefficients h^n), known as low pass filter of this function, are hv(n) = [ 1,1] (2.18) the associated Haar wavelet function can be defined by (1 0 < x < 0.5 i/>(x) = 1-1 0.5 < x < 1 (2.19) v0 elsewhere applying the relation between h^(n) and h^n) in equation (2.16), the high pass filter coefficients are obtained as Mn) = 7=[1.-1] (2-20) Figure 2.4 shows graphical presentation of the Haar scaling function and wavelet function [2, 19]. 1 T-------I----- I I-------1 I * t -if Figure 2.4 Haar Scaling Function (left) and Haar Wavelet Function (Right). 13 In order to analyze an image, a sequence of pixel intensities, corresponding to each row and column of the image, is regarded as a signal that can be filtered or convolved with the impulses response of a filter. The dyadic wavelet transform analyzes the image at different frequency bands with different resolution by decomposing it into a coarse approximation and detail information with high pass and low pass filters abbreviated as HP and LP. Wavelet coefficients at each decomposition level (characterized by scale or frequency band for a particular wavelet used) are calculated with two dimensional separable convolutions. This amounts to one dimension convolutions along the rows and the columns of the image. At any level the decomposition of the previous level image approximation a; is performed with a0 representing the original image. The multiresolution analysis involves a decomposition of the image space to a sequence of subspace vj Officially; the approximation of a function at resolution 2~j is defined as an orthogonal projection on a space vj c L2(R). The computational approach to two-dimensional discrete forward transform is shown in Fig. 2.5, which is also called the analysis filter bank, (known as Mallats multiresolution algorithm). In the figure, the rows of a; are passed through an HP and LP filter. This generates two images, which are then down-sampled by two (indicated by the down-arrow).This is every other column is discarded. The columns of these two output images are passed through an HP and LP filters and every second row is discarded [2, 10]. This generates four images: approximation image a;+1, the horizontal detail df+1, the vertical detail dj+1 and the diagonal detail df+1. 14 Rows Columns > HP 2 * HP d? uj+1 Figure 2.5 The Analysis Filter Bank for Two-Dimensional Wavelet Transform. 15 Figure 2.6 shows wavelet representation of image. It illustrates decomposition level of a two-dimensional image. Here LL, LH, HL and HH are the four sub bands, corresponding to approximation, horizontal, vertical, and diagonal edges respectively. Note that the superscripts of these sub bands indicate the numbers of the level decomposition. When further decomposition is composed, only LL is repeated. LL1 LH1 HL1 HH1 LL2 LH2 LH1 IN -J X HH2 HL1 HH1 Figure 2.6 One (left) and Two (right) Levels Two-Dimensional DWT. 2.5 Wavelet Based on the Gaussian Function The derivative of a Gaussian function has been shown to be useful in image enhancement. The Gaussian function is defined as with expected value =0, variance a and 2.71828182 (Euler's number). Gaussian functions are broadly used in statistics and probability theory, it appears as the density function of the normal distribution, which is a limiting probability distribution of complicated sums, according to the central limit theorem. In signal processing they serve to define Gaussian filters, such as in image processing where two-dimensional Gaussians are used for Gaussian blurs (or Gaussian smoothing) which is typically used to reduce image noise. Mathematically, applying a Gaussian blur to an image is the same as convolving the image with a Gaussian function. The first derivative, which is used for calculating wavelet transforms, is (2.21) (2.22) The Fourier transform of the Gaussian function is (2.23) and the Fourier transform of its first derivative is (2.24) 17 G(n) and G'(n) will be referred to respectively as cubic spline and quadratic spline functions. The mother wavelet and the scaling function of quadratic spline are shown in Figures below. Figure 2.7 The scaling function. Figure 2.8 The quadratic Spline Wavelet. 18 2.6 Wavelet Transform Modulus Maxima The wavelet transforms modulus maxima is one of the multi-scale presentations of a signal or an image. This presentation of a signal is achieved by detecting the local modulus maxima in dyadic wavelet transform across scales. Actually, a signal is equivalent to a function, thus, the edge of the signal can be viewed as such a component, whose value varies suddenly. In a two dimensional signal, for example an image, the edge refers to the sharp variation in color or gray-scale. If we focus the detection on the location of the edge not the width, then only the local maximal point needs to be found, without counting other large derivative points in this area. It corresponds to the skeleton of the edge and is effective to the isolated points. The edge detection conflicts with the noise reduction so the smoothing should not be strong, although it is necessary to smooth the original signal. In a two dimensional signal, the input image f(x,y) and smoothing function 6(x,y) are two variable functions. At the same time, 6{x,y) should have good locality and satisfies II XT.e (*< y) dxdy =1 (2 -25) At this time, function f(x,y) can be smoothed, and we shall describe this by writing (/ 6s)(x,y) if QsXx.y) = Xlll/C* ~u,y- v)0s(u,v) dudv (2.26) where 9s(u,v) = (<) ar|d s > 0 stands for a smoothing scale. The two wavelets planned to be partially derivative of a two dimension smooth function 6(x,y) along x and y directions are 19 (2.27) it necessarily satisfies the conditions Loo Loo ^ (x. y)dxdy = 0 (i = 1,2) (2.28) if xpKx.y) = ^rp1 (L J) and \pj{x,y) i xp2 (|, j), then the two dimension wavelet transforms of an image f(x,y) with respect to ^(x.y) and xp2(x,y) become the gradients of the image at multiples of dyadic scales and are obtained respectively by In other words, the calculation of the local maximal modulus of the derivative of a smoothing function along the gradient direction is equivalent to the computation of the local maxima modulus of a wavelet transform. The next step will be to determine the gradient direction of a digital signal as well as its local maximum along this direction. If the W1f{2j,x,y) and W2f(V,x,y) are defined to be the horizontal and vertical detail coefficients of the function f{x,y) smoothed at scale s = 2j respectively, then the wavelet transform modulus is bounded by The latest equations (2.30) and (2.31) will require less computation time for obtaining diagonal detailed coefficients. The local modulus maxima of W1f{s,x,y) = / ipl(x,y) and W2f(s,x,y) = / ip2(.x,y) (2.29) Mf(7), x, y) = y/\W1f(V,x,y')\2 + \W2f{2>,x,y)\2 (2.30) and the phase angle between them is (2.31) 20 Wxf{lj,x,y) along x, for constant y, are detected at points in which the edges of the image are located, of the horizontal variation of a function f{x,y) smoothed at scale of 2>. Equally, the local modulus maxima of W2f{li,x,y) along y, for constant*, are detected at points of vertical variation at a scale of 2j. They are where the modulus image Mf(2j,x,y) is locally maximum along the phase of degree Af(2j,x,y) in eight quantized directions: 0, |, n, p -1, ~ degrees for which 0 degree is the low value from the horizontal to the higher one at 90 degrees [7, 17]. 2.7 Algorithme a Trous A fast dyadic wavelet transform can be computed with a filter bank algorithm called in French algorithme a trous. It was originally developed for music synthesis using non-orthogonal wavelets but it is now of particular interest in image processing applications [3]. The algorithme a trous uses quadratic spline filters properly related to the quadratic spline wavelet (or the derivative of a Gaussian function) and cubic spline scaling function (or the Gaussian function). They have been shown to be optimum filters for image processing. During decomposition, at every scale, the filters are dilated by inserting 2; 1 zeros (Where j is the current decomposition level) between each of the coefficients of the analysis filters used in the a trous algorithm. This creates holes (a trous in French). Again the filter of the a trous algorithm is decimated at each level of reconstruction by putting 2; 1 zeros between each of the coefficients of the synthesis filters. Only horizontal and vertical details, dj1 21 and dj, respectively as well as an approximation, a7, of the image at each level are created. When compressing an image using Wavelet Transform Modulus Maxima, there are always occurring some changes in the edges. To prevent these changes, the mother wavelet should be symmetrical at a given scale for the edge location matching the Wavelet Transform Modulus Maxima. For this reason, the quadratic spline filters are chosen. Mallat suggested the quadratic spline filers are designed as follows: If 0 and 0 are the Fourier transforms of the scaling function (p and wavelet function i|/ respectively. The" translation of m + 1 convolutions of l[0,i] with itself, is given by: The low pass and high pass filters in the neighborhood of o> = 0 are designed as (2.32) with e = 1 if m is even and e = 0 if m is odd. This box spline is centered at t = ^ if m is even and at t = 0 if m is odd. V2e tez (cos 771+1 (2.33) (2.34) where h((o), g((o) are the Fourier transform of h[n], g[ri). 22 The Fourier transform of the resulting wavelet is /. \ # it i v i/i T 2 TV . 1 (0)\^(0)\ -ia> -ta)(1+Â£) (sinj\ 'K) = ^8(1)p{1) = e (-V-j (2.35) with a box spline of degree m + 2 at center of t = ^ . The quadratic spline filter coefficients are shown in Table 2.1. Here /i(n) and g(n) are the low pass and high pass analysis filters while h(ri) and g{n) are the low pass and high pass synthesis filter respectively [2,8, and 10]. Table 2.1 Spline filter coefficients used in the a trous algorithm. n h(ri) h(n) gin) gin) -2 -0.0442 -1 0.1768 0.1768 -0.3094 0 0.5303 0.5303 -0.7071 -0.9723 1 0.5303 0.5303 0.7071 0.9723 2 0.1768 0.1768 0.3094 3 0.0442 23 Figure 2.9 shows the decomposition filter banks for the two dimensional a trous algorithm. The approximation is generated by first shifting the columns of a.j to the left, 2;-1 times. Then, the columns are filtered using the low pass filter h[n] and the rows are shifted to the left 2; 1 times and filtered with a low pass filter. The result is the approximation a;+1 of the image. The vertical details are calculated by first filtering the rows of a; using the high pass filter and then shifting the rows, 2j times. The horizontal details are calculated by first filtering the columns of a; using the high pass filter, g[n] and then shifting the columns, 2j times. Columns Rows aj+1 dj+1 Figure 2.9 Filter bank for two dimensional image decomposition implementing the algorithme a trous. 24 To reconstruct the original image, reverse the order from the decomposition procedures with the exception that the shifting is altered. At each level of reconstruction, the filters are decimated by inserting 2; 1 zeros between each coefficient. The output of the horizontal detail coefficients and high pass filter, the vertical detail coefficients and high pass filter and the sum of the approximation and low pass filters, are added up. The results are then multiplied by ^ . Rows Columns aj+i Figure 2.10 Filter bank implementing image reconstruction with the "algorithme a trous". 25 3. Clustering preview Clustering is grouping of similar objects or is the allocation of a set of observations into subsets, so that observations in the same cluster are similar in some area or category. Clustering is a usual technique for statistical data analysis and a method of unsupervised learning used in many fields, including image analysis, data mining, pattern recognition, and other applications. The term partitioning is sometimes used as synonyms for clustering, this term is frequently used for approaches outside the traditional bounds of cluster analysis. For example, the term partitioning is often used in connection with techniques that divide images into sub-images and that are not strongly connected to clustering. The goal of clustering is to determine the essential grouping in a set of unlabeled data [16]. Generally, clustering has many applications which could be classified as follows: Marketing: classifying the customers to groups based on their action, buying patterns and their needs according to some specification, and the latest buying habits; Biology: grouping of animals and plants based on their specification; Data Mining: data finding based on previous data structures; City Mapping and Planning: identifying groups of houses based on their house type and geographic location; Earthquake Studies: finding the potential hazardous area based on previous data; 26 Web: document grouping or viewer alliance to the relevant sites; Speech Recognition: building code book from specifics and categorizing the speech according to the dialogue or speech compression; Image Segmentation: classifying the medical or Satellite images. Figure 3.1 illustrates clustering based on color property with a simple example. The MATLAB program kmn.m is explained in Appendix B. In this case, I can easily identify the two clusters into which the data can be divided. Figure 3.1 Clustering based on color property. 27 3.1 K-means Algorithm There are many clustering algorithms. One of the well known methods that could be used for clustering is K-means algorithm. It was introduced by James MacQueen in 1967. In this algorithm we select /reenters for Â£ clusters and every cluster center is determined by calculating the mean of all data in that cluster; this is why it is called the k-means algorithm. K-means algorithm is widely used for the clustering large set of data. This method is a partitioning clustering algorithm that classifies the given data objects into k different clusters through the iterative algorithm converging to a local minimum. The algorithm consists of two separate stages. The first stage selects k centers randomly, where the value of /ris initially determined. The second one is to take each data object to the nearest center. Euclidean distance determines the distance between each data object and the cluster centers. The Euclidean distance d(*i,yj) between one vector* = (x1,x2,x3,...,xn) and another vector y = (y\,y2>yz -,yn) can be obtained as follows: When all the data objects are included in some clusters, the first step is completed and an early grouping is finished. Recalculating the average of the early form clusters is done by iterative processes that continue repeatedly until the criterion function becomes the minimum. The criterion function is dixi.yi) = [Â£?=10; y*)2]1^ (3.1) 2 (3.2) where Cj is cluster center and is data objects. E is the sum of the squared error of all objects in the data set. 28 The Euclidean distance x,a) Ci is the criterion function, which is used for determining the nearest distance between each data objects x[j) and cluster center Cj. In simple words, the value of E has to be minimum to make sure the clustering has been done properly [15]. The process of K-means algorithm proceeds as follows: a: we select k initial centers for k clusters from the data, which we may identify as Cj C2........ Ck This selection could be done randomly or based on original data distribution. b: as we sweep the data (pixels for image), we associate each pixel with one of these k clusters. This association is based on the distance between the pixel and the cluster center, which should be the shortest compared to other cluster centers. c: when the association is established and the pixel becomes a member of the cluster, the calculation begins to calculate the new center of the cluster based on a new added member. The calculation is done simply by obtaining the mean of the members coordinates in every direction. In this case, we obtain the mean of x and y values for all new members in each cluster. d: the above steps will be repeated for all new members and new centers will be calculated as the new members join the clusters, the operation stops as the centers of all clusters do not change position anymore. 29 The K-means algorithm can be illustrated in the following simple flowchart. f......*.........i j Number of / / clusterK_______j I Centroid Distance objects to centroids I Grouping based on minimum distance Flrure 3.2 Flowchart of the k-means algorithm. 30 The numerical example below in table 3.1 is provided to clarify this simple iteration. Suppose we have four students and each one has two attributes or features as shown in the table below. Our goal is to group these objects (students) into K=2 groups based on the specified two features (grade and age). Table 3.1 Example data of k-means algorithm. Object attribute 1 (X): Age attribute 2 (Y): Grade Student A 6 1 Student B 7 1 Student C 8 3 Student D 10 4 Each student represents one point with two attributes (X,Y); we can represent it as a coordinate in an attribute space as shown in figure below. Iteration 0 o < 0 2 4 6 8 10 12 X: Age 31 1. Suppose we select student A and student B as the first centroids. Let cx = (6,1) and c2 = (7,1) notate the coordinate of these centers. Iteration 1 X: Age 2. We calculate the distance (Euclidean distance) between cluster centroid and each object, then we have a distance matrix at iteration 1 as: D! [0 1 2.83 5 ] c1 = (6,1) Group 1 ll 0 2.23 4.24J C2 = (7,1) Group 2 Each column in the distance matrix symbolizes the object. The first row of the distance matrix corresponds to the distance of each object to the first centroid and the second row is the distance of each object to the second centroid. For example, distance from student C = (8,3) to the first centroid 6^(6,1) is yj(8 6)2 + (3 l)2 = 2.83, and its distance to the second centroid c2 = (7,1) is V(8 7)2 + (3 l)2 = 2.23 etc. 3. We assign each object based on the minimum distance. Thus, student A is assigned to group 1, student B to group 2, student C to group 2 and 32 student D to group 2. The element of the group matrix below is 1 if and only if the object (student) is assigned to that group. ri ri 0 0 0] Group l lo 1 1 lJ Group 2 4. Knowing the members of each group, we now compute the new centroid of each group based on these new memberships. Group 1 only has one member thus the centroid remains in ^(6,1). Group 2 now has three members, thus the centroid is the average coordinate among the three members: C2 = ^7++1 f --+i.+4^ = (8.3,2.6). Iteration 2 X: Age 5. The next step is to compute the distance of all objects to the new centroids. Similar to step 2, we have a distance matrix at iteration 2: D2 = \ 0 1 2.83 5 Ci = (6,1) Group 1 b.80 2.06 0.5 2.20J C2 = (8.3,2.6) Group 2 6. Similar to step 3, we assign each object based on the minimum distance. According to the new distance matrix, we move the student B to Group 1 33 while all the other objects remain the same. The Group matrix is shown below G2 110 0] Group 1 0 0 1 lJ Group 2 7. Now we repeat step 4 to calculate the new centroids coordinates based on the clustering of the previous iteration. Groupl and group 2 both have two members, thus the new centroids are Cx = - (6.5,1) and Ci = (8^0^) = (9)35) Iteration 3 X: Age 8. Repeat step 2 again and we have new a distance matrix at iteration 3 as D3 r 0.5 0.5 2.5 4.61] C1 = (6.5,1) Group 1 U.91 3.20 1.11 2.06J C2 = (9,3.5) Group 2 9. Again, we assign each object based on the minimum distance. r3 rl 1 0 0] Group 1 lo 0 1 lJ Group 2 The obtained result is G2 = G3 Comparing the grouping of the last iteration and this iteration reveals that the objects do not move groups anymore. Thus, 34 the computation of the k-means clustering has reached its stability, no more iteration is needed, and we achieved the final grouping. 3.2 Image Clustering The image /, can be defined as a rectangular matrix (called image matrix) with image rows x and image columns y. A row number together with a column number defines a small image area called a pixel (from picture element), which is assigned a value that represents the brightness of the pixel. A gray-scale or black and white image is also called bi-level or monochromic. The gray-scale value 0 corresponds to black and a value of 255 corresponds to white. Such an image is called an 8- bit gray-scale image with where G = {0, 1, 2,....., 255} that can be constructed from the 8-bit, which represents 28 or 256 possible values. On the other hand, sometimes 256 values may not be enough for displaying all of the information an image contains. Typically, a color image cannot be represented by the method shown in Equation (3.3). In that case, we can extend the image / by using more planes such as: 1 = fix,y) (3.3) f{x,y) G G (3.4) / = f(x,y,n) (3.5) with n representing the plane index [1]. 35 The following process will be performed for image segmentation in MATLAB using K-means clustering on the sample image (house.bmp ). The program is attached in Appendix A. First the image is read by using the Matlab command imread, then the image input data is brought into a 256x256 matrix. Each element of this matrix represents an identical pixel in the image and will carry a number representing the gray level or intensity level of the same pixel. By using the reshape command, the image matrix (256x256) is converted to a column matrix (1x65536). The K-means algorithm is used and the clustering is performed by selecting 3 initial centers for k. These 3 centers were initially selected randomly and as K-means algorithm is performed, the clusters start to take shape and the centers are replaced after each iteration. K-means clustering is done when there are no more changes in the cluster centers. Figure 3.3 illustrates the original image and the result of K-means clustering by selecting 3 clusters. In the first cluster, the final position of center (ct) is detremined to be 71.43 and all the pixels relatively close to this center are collected in this cluster. It can be seen that darker pixels are selected for the first cluster. In the second cluster c2 = 120.68 and the gray pixels are selected. In the third cluster C3 = 188.87 and all light gray or white pixels are sleeted. 36 VdteMM(1 Oipct* civrt* 2 OljUtt rtuitw S Figure 3.3 The K-means algorithm on a gray-scale image. 37 Now, the k-means algorithm on the RGB or color image is to be performed. The colored images are maintained by a combination of three main colors which are Red, Green and Blue. The spectrum of intensity for these colors will create a vast range of natural colors. The following steps are taken to process the color image: 1. Read image; 2. Convert image from RGB color space to L*a*b color space; 255 Figure 3.4 Color expression by L*a*b concept. The L*a*b* space consists of: a luminosity layer 'L*'; a chromaticity-layer (it is an objective specification of the quality of a color regardless of its luminance); 'a*' indicating where color falls along the red-green axis and chromaticity- layer; and 'b*' indicating where the color falls along the blue-yellow axis. All of the color information is in the 'a*' and 'b* layers. 3. Apply K-means clustering; Use K-means to cluster the objects into three clusters using the Euclidean distance metric. Label every pixel in the image and create images that segment the original image by color. 38 Figure 3.5 shows K-means clustering on the color image. The file cell3.jpg (which can be found on the web) was read by using the Matlab command imread' that means it is changed into a 256x256 matrix. Then converting the color image from the RGB form to L*a*b color space by using the makecforni Matlab command. The segmentation was performed using the K-means algorithm for clustering in three levels. The program is attached in Appendix C. The first output figure (left) shows the labeling of every pixel in the image with its cluster-index. The second figure represents cluster 1 with a center coordination of Cj (150.89, 106.92) that displays the background. The second cluster C2 (173.52, 66.71) shows the purple cells and the last image (right) with C3 (129.96, 123.15) represents the pink cells. 39 Input Image Image labeled by cluster index Objects in cluster 1 Objects in duster 2 Objects m cluster 3 Figure 3.5 The K-means algorithm on color image 40 4. Segmentation Image segmentation is the process of assigning a label to every pixel in an image such that pixels with the same label share certain visual characteristics. Segmentation refers to the process of partitioning a digital image into multiple segments (sets of pixels). Segmentation precision determines the final success or failure of computerized analysis procedure. Image Segmentation is one of the most difficult tasks in image processing; for this case considerable care should be taken to improve the probability of accurate segmentation. The goal of segmentation is to simplify or change the representation of an image into something that is more meaningful and easier to analyze. The result of image segmentation is a set of segments that collectively covers the entire image. Each of the pixels in a region is similar in respect to some characteristics or computer properties, such as color, intensity or texture [6]. 4.1 Segmentation Method There are many approaches and techniques for image segmentation. In this section several techniques are presented. The simplest method of implementation is amplitude segmentation. This is a variant of a tresholding technique. Tresholding techniques separate the pixels into baskets based on a set of threshold values for some identifying characteristics. In this method the identifying characteristic is the pixel based value. This can be either gray-scale or color intensity amplitudes. Here the 41 pixels amplitude is compared to a set of predetermined thresholds and placed in the relevant region. Suppose that the gray-scale histogram shown in Fig. 4.1(b) corresponds to an image, f{x,y) in such a way that object and background pixels have gray levels grouped into two dominant modes. One obvious way to extract the objects from the background is to select a threshold T that separates these modes. Then any point (x,y) for which f(x,y) > T is called an object point; otherwise the point is called a background point. When/(x,y) is the gray level of point (x,y), a threshold image g(x,y) is defined as if f(x,y)>T if f{x,y) (4.1) Thus pixels labeled one correspond to objects, whereas pixels labeled zero correspond to the background [6, 11, and 14]. Segmentation is then accomplished by scanning the image pixel by pixel and labeling each pixel as object or background, (using im2bW command) depending on the gray level of that pixel, it is greater or less than the value of T. Figure 4.1(c), shows the result of global threshold with T between the maximum and minimum gray level. The product obtained using T = 111 (with mean (ff command) segmented the original image in Fig.4.1 (a). The program is located in Appendix F. This method would be useful in solid color images, line art or text. Amplitude tresholding performs poorly with images or when the input is of a texture type even though it is the least expensive. 42 (a)Orginal Image (b):Histogram of Original Image 60001---------1------------------------------- Thereshol Level = 111 Figure 4.1 (a) Original image, (b) Image histogram, (c) Segmentation of original image. 43 The second method of segmentation is clustering and K-means clustering was selected because the K-means method is numerical, unsupervised, iterative and non-deterministic. This method performs well on image and texture inputs. The disadvantage with this method is pre-setting the k value (number of clusters), which makes it difficult to find the optimal k value. However, the clear advantage is that with a large amount of data, it is faster and will produce tighter clusters. The clusters are flat and they do not overlap. It is well suited to generating globular clusters. The K-means algorithm has been completely described in chapter 3. The third method of segmentation involves edge detection. Edges show immense local changes of intensity in an image and typically occur on the boundary between two different regions within an image. It is usually used to produce a line drawing of a scene from an image showing important features (e.g., corners, lines, curves). Several approaches are suitably utilized for edge detection. One of the common methodologies used in edge detection is Sobel operators. In this method, the gradient of the image intensity is calculated at each point, giving the direction of the largest possible increase from light to dark and the rate of change in that direction. The result therefore shows how "sharply" or "smoothly" the image changes at that point and therefore how likely it is that part of the image represents an edge, as well as how that edge is likely to be oriented. Figure 4.2 illustrates a 3x3 region of the image (the Z's are gray- level value) and Sobel operators. 44 Zi z2 Z3 z4 Zs Ze z7 N 00 Zg (a) -1 -2 -1 0 0 0 1 2 1 (b) 1 0 1 -2 0 2 -1 0 1 (c) Figure 4.2 (a) 3x3 region of an image (b) component in the x-direction. (c) component in the y-direction. The mathematical approach is to use masks of size 3x3 and the following formulas: Gx = (z7 + 2z8 + z9) (Zi + 2z2 + z3) (4.2) Gy = (z3 + 2z6 + z9) (zx + 2z4 + z7) (4.3) In this formulation, the difference between the first and the third rows of 3x3 image regions approximates the derivative in the x-direction, and the difference between the third and first columns approaches the derivative in the y-direction. A weight value of 2 is used to achieve some smoothing by giving more importance to the center point. At each point in the image, the resulting gradient approximations can be combined to give the gradient magnitude, using: 45 (4.4) V/ = magiVf) = [G$ + Gy2p
the gradient's direction can also be calculated:
a(x,y) = tan_1(|i) (4.5)
where, a is 0 for a vertical edge which is darker on the left side.
Figure 4.3 shows the Sobel gradient of the image (house.bmp). The operator
uses two 3x3 masks which are convolved with the original image to calculate
approximations of the derivatives, one for horizontal changes, and one for
vertical. If we define 4.3(a) as the original image, then Fig. 4.3(b) and Fig.
4.3(c) are two images which at each point contain the horizontal and vertical
derivative approximations. Figure 4.3(d) illustrates the resulting gradient
approximations (horizontal and vertical) that are combined and gives the
gradient total magnitude image. The Sobel method has a major drawback of
being very sensitive to noise. Another weakness is the size of the masks that
are fixed and cannot be adapted to a given image [4, 6].
46

Original image
Finding horizontal edge with SOBEL Method
Figure 4.3 (a) Orignal image, (b) \GX(Component in the x-direction.
(c) \Gy\ Component in the y-direction. (d) Gradient image,V/ = \GX\ + |Cy|.
47

Another method used for edge detection is Wavelet Transform Modulus
Maxima (WTMM). This method is for fast feature extraction. By using the
Modulus Maxima in a two dimensional dyadic wavelet transform across
scales, edges are detected. Modulus Maxima and a trous algorithm have
been explained thoroughly in chapter 2.
Figure 4.4 shows a 256x256 pixel gray-scale image of Lena and the mesh
diagram of this image is also shown. Mesh plot is the basic command for
three dimensional graphing and it defines a surface by the z-coordinates of
points above a grid in the x-y plane using straight lines to connect adjacent
points.
Figure 4.5 illustrates the Modulus Maxima analysis for this image (lena.bmp)
at different scales. After the image is read, the high pass (HP) and low pass
(LP) filters of the corresponding scale 2; where j El are calculated. The
coefficients in Table 2.1 are the spline filter coefficients used in the a trous
algorithm. Then, the wavelet transform modulus and angle of the wavelet
transform modulus for each corresponding pixel location is generated. It is
also the direction where the function has the sharpest variation. This is done
for each level of decomposition. The wavelet transform modulus Mf(2j,x,y)
and angle Af(2j,x,y) are calculated by equations of (2.30) and (2.31). The
outcome of this analysis of detected edges in three levels is shown. WTMM
has great advantages of being resistant to noise and faster than Sobel
method. The Matlab program is attached in Appendix D.
48

Figure 4.4 3-D mesh plot for Lena.bmp.
Modulus Maxima
Amplitude coeficent at level 1
Modulus Maxima
Amplitude coeficent at level 2
Modulus Maxima
Amplitude coeficent at level 3
/
V
Figure 4.5 Modulus Maxima at 3 levels.
49

The algorithm represented in this thesis, which perform image segmentation
is a combination of K-means and a trous algorithms. First, I applied table (2.1)
values to set LP and HP filters in a trous algorithm, then horizontal and
vertical details were calculated. The wavelet transform modulus and angle of
the wavelet transform modulus for each corresponding pixel location was
generated. The image edges were found. Next, the output will be provided as
an input image to K-means algorithm. The end result of K-means algorithm is
the segmented image. This process is represented in the Block Diagram
below:
Image

Modulus
Maxima

Edge
Detection
Select the Image
thont Noise
K-Means Image
Algorithm Segmentation
Figure 4.6 Block Diagram of image segmentation.
Figure 4.7 shows the selected image. The first step is to read the image.
The second step is to initialize a trous filter, which means calculating LP
and HP decomposition filters corresponding to a scale of2j. In the third
step, Modulus Maxima is calculated. The outcome of this process is
shown in Fig. 4.8. As you can see the second image has the lowest noise
that is why it is selected.
In Fig. 4.9, the edges holes will be filled by using the imfilt command. It
is necessary to have clear and solid edges.
In the last step, I apply K-means algorithm (by assuming k=2), which
provides a clear segmented image.
50

Figure 4.10, shows the image divided into two segments. This can
separate the image and can effectively classify the foreground and
background. The Matlab program is attached in Appendix E.
Orifpnal Input knag*
Figure 4.7 Original image and its 3-D mesh plot.
Modulus Maxims
AmpMuda coaUcent st Laval 1
Modulus Mawna
AmpMuds coaftcant at Laval 2
Modulus Maxima
Amphtuda coatScant at Laval 3
Figure 4.8 Edge detection with Modulus Maxima
51

Detecting Edge
Figure 4.9 Edge Detection with all holes filled.
SEGMENT 2
Figure 4.10 Segmented object.
52

5. Conclusion
It is perceived that the wavelets are an important tool for analyzing and
processing non-stationary signals. In this thesis, a possibility of wavelet
analysis application for obtaining image segmentation algorithm less sensitive
to noise included in original image has been considered. The image
segmentation was introduced as one of the most fundamental and important
procedures in image processing, which is typically used to locate objects and
boundaries in images.
In this thesis, a few segmentation techniques were examined. Initially,
tresholding technique was investigated but it was found that even though it
works well for solid color images, it is unsuitable for textures. Next, Sobel
operators application to edge detection was demonstrated but this approach
turned out to be very sensitive to noise.
It is well known that most of the noise contained in the original image is
extracted in the details of the first level of its wavelet decomposition.
Therefore, the details at the second decomposition level are almost noise
free. The modulus maxima representation of these details determines the
noise-free representation of the edges of the original image, which is provided
as an input image to K-means algorithm. Therefore, the new segmentation
algorithm, developed and investigated in this thesis, consists of the image
edge detection using wavelets modulus maxima of the original image and of
the K-means image segmentation algorithm.
In order to obtain the wavelet transform of the image in a suitable form to get
this transform representation by the modulus-angle pair, the a trous wavelet
53

decomposition algorithm was employed. This algorithm assigns two numbers
to every image pixel at every decomposition level: a horizontal and a vertical
detail, which can be alternatively represented by a modulus and an angle.
Such representation enables easy detection of the image edges.
The main contribution of this thesis is the development of this new algorithm,
which performs image segmentation and simultaneously filters out the noise
present in the original image.
54

APPENDIX A
K-means Algorithm MATLAB Program
%The following program clustered a gray-scale image which is used house.bmp in this
%process. The Algorithm has been described in chapter 3.
clc
clear all
close all
%..............................
%..............................
%..........................................................
%% Step 2: Convert the image input data to column matrix
%..........................................................
column_matrix = double(he);
[mrows, ncols] = size(column_matrix);
column_matrix = reshape(column_matrix,mrows*ncols,1);
%..............................
%% Step 3: Use k-means
%..............................
Ncluster = 3;
% repeat the clustering 3 times
[clusterjdx cluster_center] = kmeans(column_matrix,Ncluster,...
'distance','sqEuclidean',...
'Replicates',3);
pixeljabels = reshape(cluster_idx,256,256);
segjmages = cell(1,Ncluster);
for k = 1 :Ncluster
55

segmentedjmages = he;
for i = 1:mrows;
for j = 1:ncols;
if pixel_labels(i,j) ~= k;
segmentedjmages(i,j) = k;
else segmentedjmages (i,j) = 256;
end
end
end
segjmages {k}= segmentedjmages;
end
%..................
%% Ploting Section
%..................
figure
subplot(2,3,1:3),imshow(he), title('lnput Image');
subplot(2,3,4),imshow(segJmages{1}), title('Objects in cluster 1');
subplot(2,3,5),imshow(segJmages{2}), title('Objects in cluster 2');
subplot(2,3,6),imshow(segJmages{3}), title('Objects in cluster 3');
56

APPENDIX B
Kmn.m Matlab Program
%ln this program I have clustered the image according to the color property.
clc
clear all
close all
%Display image
X = [randn(100,2)+ones(100,2);...
randn(100,2)-ones(100,2)];
opts = statset('Display','final');
%use K-means algorithm
[idx.ctrs] = kmeans(X,2);
% the red pixcels are shown as o
plot(X(idx==1,1 ),X(idx==1,2),'ro','MarkerSize',10)
Centers = ctrs
hold on
% the blue pixcels are shown as
plot(X(idx==2,1),X(idx==2,2),'b*','MarkerSize',10)
plot(ctrs(:,1),ctrs(:,2),'kx',...
'MarkerSize',12,'LineWidth',2)
plot(ctrs(:,1),ctrs(:,2),'ko',...
'MarkerSize',12,'LineWidth',2)
legend('Cluster 1'.'Cluster 2','Centroids',...
'Location','NW')
57

APPENDIX C
K-means Algorithm on Color Image
%This program has been run on colored image by using k-means algorithm for segmentation.
clc
clear all
close all
%............................
%............................
%.....................................................
%% Step 2: Convert the image to L*a*b* color space
%.....................................................
cform = makecform('srgb2lab');
lab_he = applycform(he, cform);
%...........................................................................
%% Step 3: Use K-means to cluster the objects in the 'a*b' into three clusters
%...........................................................................
ab = double(lab_he(:,:,2:3));
nrows = size(ab,1);
ncols = size(ab,2);
ab = reshape(ab,nrows*ncols,2);
nColors = 3;
% repeat the clustering 3 times to avoid local minima
[clusterjdx cluster_center] = kmeans(ab,nColors,...
'distance','sqEuclidean',... 'Replicates',3);
%...........................................................................
%% Step 4: Label every pixel in the image using the results from K-means
%...........................................................................
58

pixeljabels = reshape(cluster_idx,nrows,ncols);
%% Create Images that Segment the input Image by Color,
segmentedjmages = cell(1,3);
rgbjabel = repmat(pixeljabels,[1 1 3]);
for k = 1:nColors
color = he;
color(rgbjabel ~= k) = 0;
segmented_images{k} = color;
end
%...................................
%% Step 4: Ploting Section
%...................................
Figure
subplot(2,4,1:4),imshow(he), title('lnput Image');
text(size(he,2),size(he,1 )+15,...
Sample Microscopic Cellular Image file ...
'FontSize',7,'HorizontalAlignment7right');
subplot(2I4,5),imshow(pixelJabels,[]), title('lmage labeled by cluster index');
subplot(2,4,6),imshow(segmentedJmages{1}), title('Objects in cluster 1');
subplot(2,4,7),imshow(segmentedJmages{2}), title('Objects in cluster 2');
subplot(2,4,8),imshow(segmentedJmages{3}), title('Objects in cluster 3');
59

APPENDIX D
Edge Detection with Modulus Maxima
% This program has been run on an image by using Modulus Maxima for segmentation,
clc
clear all
close all
%..............................................
%input:
%..............................................
% L=number of level of decomposition
% y=original image
% h=low pass decomposition filter
%g= high pass decomposition filter
%gp=high pass decomposition filter
%..............................................
%output:
%..............................................
%y=last level approximation image
%gprime=reconstruction high pass filter
%hprime= reconstruction low pass filter
%finaly the program displays the modulus maxima matrix for each level
%.........................................................
% Read the Image and Initialize Stage
imshow(X,map)
%change data type of read image for applying computation
house =double(X); y=house;
60

% Initialize the composition level
decomp_times = 3;
%.........................................
%Spline wavelet filters
%.........................................
hf=[.125 .375 .375 ,125].*sqrt(2);
hprime = hf;
gf=[.5 -.5].*sqrt(2);
gprime = [-0.03125 -0.21875 -0.6875 0.6875 0.21875 0.03125].*sqrt(2);
save_orig = y;
L=decomp_times;
%L=level of detail
[nr,nc]=size(y);
hd=hf;
gd=gf;
%........................
% Loop for decomposing Input Image
%........................
for k = 1:L
%.............................
%calculate approximations
%.............................
%shift either image or latest approximation vertically
for i=1:nc
a{k}(1:nr,i) = leftshift(y(1:nr,i),2A(k-1));
%circular convolution of hd with colums from vertically shifted image
a{k}(1:nr,i) = cconv(hd,a{k}(1:nr,i));
end
%shift either image or latest approximation horizontally
for i=1 :nr
a{k}(l,1 :nc) = leftshift(a{k}(l,1:nc),2A(k-1));
%circular convolution of hd with rows from horizontally shifted image
61

a{k}(l,1:nc) = cconv(hd,a{k}(l,1:nc));
end;
%.............-.............
%calculate details d1 and d2
%...........................
%shift image or previous approximation vertically
for i=1 :nc
d1 {k}(1 :nr,i)=cconv(gd,y(1 :nr,i));
d1{k}(1 :nr,i)=leftshift(d1 {k}(1 :nr,i),2Ak);
%d1 comes from convolving columns of shifted image
end;
%shift image or previous approximation horizontally
for i=1 :nr
d2{k}(l,1 :nc)=cconv(gd,y(l,1 :nc));
d2{k}(l,1 :nc)=leftshift(d2{k}(l,1 :nc),2Ak);
%d2 comes from convolving rows of shifted image
end;
end;
%.............................................................
%% Calculate the Modulus Maxima for each level
%.............................................................
to = degrees;
for k = L:-1:1
wtm{k} = sqrt(d1{k}.A2 + d2{k} A2); %The Wavelet Transform Modulus
TMP = atan (d1{k}./(0.0000001+d2{k}));
Ammfk} = unitsratio(to, from) TMP; % The phase
[mod_m,mod_n] = size(wtm{k});
%.............................................................
%%calculate a Trous filters
%.............................................................
62

hprime=hd;
y = a{k};
end
%..............................................
%Display image
%..............................................
figure;
for k = 1 :decomp_times
subplot(1 ,decomp_times,k)
imshow(wtm{k},[]);
title({Modulus Maxima;[Amplitude coeficent at level ,int2str(k),]})
end
63

APPENDIX E
%This program has been run on an image by using Modulus Maxima and K-means clustering
%for segmentation.
clc
clear all
close all
%..................................
%input:
%..................................
% L=number of level of decomposition
% y=original image
% h=low pass decomposition filter
%g= high pass decomposition filter
%gp=high pass decomposition filter
%..................................
%output:
%y=last level approximation image
%gprime=reconstruction high pass filter
%hprime= reconstruction low pass filter
%finaly the program displays the modulus maxima matrix for each level
decomp_times = 3;
%............................
%Spline wavelet filters
%............................
hf=[.125 .375 .375 .125].*sqrt(2);
hprime = hf;
gf=[.5 -,5].*sqrt(2);
64

%gprime = [0 0 -.5 .5 0 0].*sqrt(2);%[-0.03125 -0.21875 -0.6875 0.6875 0.21875
0.03125].*sqrt(2);
gprime = [-0.03125 -0.21875 -0.6875 0.6875 0.21875 0.03125].*sqrt(2);
%.................................................
%.................................................
imshow(X,map); title('Original Input Image')
house = double(ind2gray(X,map));
y=house;
save_orig = y;
L=decomp_times;
%L=level of detail
[nr,nc]=size(y);
hd=hf;
gd=gf;
%.................................................
%Loop for decomposing House
%.................................................
for k = 1:L
%.............................
%calculate approximations
%.............................
%shift either image or latest approximation vertically
for i=1 :nc
a{k}(1:nr,i) = leftshift(y(1:nr,i)',2A(k-1))';
%circular convolution of hd with colums from vertically shifted image
a{k}(1:nr,i) = cconv(hd,a{k}(1:nr,i)')';
end
%shift either image or latest approximation horizontally
for i=1 :nr
a{k}(i,1:nc) = leftshift(a{k}(i,1 :nc),2A(k-1));
65

%circular convolution of hd with rows from horizontally shifted image
a{k}(i,1 :nc) = cconv(hd,a{k}(i,1 :nc));
end;
%..........................
%calculate details d1 and d2
%..........................
%shift image or previous approximation vertically
for i=1 :nc
d1 {k}(1 :nr,i)=cconv(gd,y(1 :nr,i)')';
d1{k}(1:nrJ)=leftshift(d1{k}(1:nrJ)\2Ak)';
%d1 comes from convolving columns of shifted image
end;
%shift image or previous approximation horizontally
for i=1 :nr
d2{k}(i,1 :nc)=cconv(gd,y(i,1 :nc));
d2{k}(i,1 :nc)=leftshift(d2{k}(i,1 :nc),2Ak);
%d2 comes from convolving rows of shifted image
end;
%..........................
%calculate a Trous filters
%..........................
hprime=hd;
y = a{k};
end;
for k = L:-1:1
wtm{k} = sqrt(d1 {k}.A2 + d2{k} A2); % Wavelet Transform Modulus
tmp = atan (di{k}./(o.000000l+d2{k})); % The phase
[mod_m,mod_n] = size(wtm{k});
End
66

%...............................................
figure;
for k = 1 :decomp_times
subplot(1 ,decomp_times,k)
imshow(wtm{k},[]);
title({'Modulus Maxima';['Amplitude coefficent at Level ',int2str(k),"]})
end
%..........................
%% Edge Detection
%..........................
AA = im2uint8(mat2gray(log(1+double(wtm{1}))));
DD = im2uint8(mat2gray(log(1+double(wtm{2}))));
HH = im2uint8(mat2gray(log(1+double(wtm{3}))));
HISt = histeq(DD,256);
%...............................................
%...............................................
figure;
subplot(1,2,1);
% Histogram of Maxima Amplitude Level 2
imhist(DD);
title('Madulus Maxima Amplitude Level 2 Matrix Histogram');
subplot(1,2,2);
imhist(HISt);
title('Madulus Maxima Amplitude Level 2 Matrix After Enhanced Contrast Histogram')
%%
SA =DD;
G = imfill(SA);
%..................................................
%% Convert the image input data to column matrix
%..........................................................
he = G;
column_matrix = double(he);
[mrows, ncols] = size(column_matrix);
67

column_matrix = reshape(column_matrix,mrows*ncols,1);
%%.....................................................
wname = 'coif2'; lev = 3;
nbc = size(map,1);
[c,s] = wavedec2(G,lev,wname);
% Estimate the noise standard deviation from the
det1 = detcoef2('compact',c,s,1);
sigma = median(abs(det1 ))/0.6745;
% Use wbmpen for selecting global threshold
% for image de-noising.
alpha = 1.2;
thr = wbmpen(c,1 .sigma,alpha);
% Use wdencmp for de-noising the image using the above
% thresholds with soft thresholding and approximation kept,
keepapp = 1;
Gd = wdencmpCgbl'.c.s.wname,lev,thr,'s',keepapp);
%figure;
%subplot(1,2,1),imshow(G);title('Sourse of De-Noising Algorithm')
%subplot(1,2,2),image(wcodemat(Gd,nbc)); title('De-Noised Image')
%...............................
%% Use k-means algorithm
%...............................
Ncluster = 2;
% repeats the clustering 2 times
[clusterjdx cluster_center] = kmeans(column_matrix,Ncluster,...
'distance','sqEuclidean',...
'Replicates',2);
[m,n]=size(he);
pixeljabels = reshape(cluster_idx,m,n);
segjmages = cell(1,Ncluster);
TR=mean(mean(SA));
for k = 1:Ncluster
segmentedjmages = he;
for i = 1 :mrows;
68

for j = 1 :ncols;
if pixel_labels(i,j) ~= k;
segmentedjmages(i,j) = 0;
else segmentedjmages (i,j) = X(i,j);
end
if G(i,j)>=TR+TR/3.5;
EDg(iJ) = 0;
else
EDg(i,j) = 256;
end
end
end
segjmages {k}= segmentedjmages;
end
%............................................
%Display image
figure;
imshow(EDg);title('Detecting Edge')
%............................................
figure;
for k = 1 :Ncluster
subplot(1,Ncluster,k);
imshow(segJmages{k},[]);
title(['SEGMENT ',int2str(k)],'Color','b');
end
69

APPENDIX F
Tresholding Matlab Program
%This program has been run on an image by using Global Thresholding,
clc
clear all
close all
%calculate threshold T
T = mean(mean(f));
g = im2bw(f,T/255);
Hst = imhist(f);
%output
Y=1:50:max(Hst);[m ,n]=size(Y); X = T*ones(n);
%display image
subplot(2,2,1);imshow(f);title('(a):Orginal Image','FontSize',13);
subplot(2,2,[2 4]);area(Hst);title(' (b):Histogram of Original lmage','FontSize',13 );
hold on; line(X,Y,'LineWidth',2,'Color',[.9 .2 ,5]);figure(gcf);
xlabel(['Thereshol Level = \int2str(T)]);
subplot(2,2,3);imshow(g);title(' (c):Segmented Image','FontSize',13);
70

APPENDIX G
Sobel Operators Matlab Program
%This program has been run on an image by using Sobel gradient,
clc
%Horizontal edge with Sobel methode
clear all
close all
[gv,t] = edge(l,'sobel','horizontal');
%display image
subplot(1,2,1);imshow(l);title('Original Image','fontsize', 12);
subplot(1,2,2);imshow(gv);title('Finding horizontal edge with SOBEL Method','fontsize', 12);
%..........................................................................................
%Vertical edge with Sobel methode
clear all
close all
[gv,t] = edge(l,'sobel','vertical');
%display image
subplot(1,2,1 );imshow(l);title('Original Image','fontsize', 12);
subplot(1,2,2);imshow(gv);title('Finding vertical edge with SOBEL Method','fontsize', 12);
%..........................................................................................
clc
clear all
close all
71

% use Sobel
BW1 = edge(l,'sobel');
%display image
subplot(1,2,1 );imshow(l);title('Original Image','fontsize', 12);
subplot(1,2,2);imshow(BW1);title('SOBEL Method','fontsize', 12);
%..............................................................
72

REFRENCES
[1] Thomas klinger Image processing with lab view and IMAQ vision
Prentice Hall 2003
[2] S. Mallat, Multiresolution Approximation and Wavelet Orthornormal
Bases of L2(R), Trans. American Math. Soc., 1989
[3] M. Shensa, The Discrete Wavelet Transform: Wedding the a trous and
Mallat algorithms, IEEE Transaction on a signal processing, Vol. 40, No
10, October 1992.
[4] Marr, D. Hildreth E. Theory of Edge Detection Proc. R.
Soc.land.B207, 198-217, 1980
[5] N G Kingsbury, Image processing with complex wavelets, Phil. Trans.
Royal Society London A. September 1999, Special issue for the
discussion meeting on Wavelets: the key to intermittent information?
(held Feb 24-25, 1999), 357, pp 2543-2560.
[6] Rafael C. Gonzalez, Richard E. Woods, Digital Image Processing
prentice Hall 2002.
[7] T. Hsung, D. Lun, W. Siu, A Deblocking Technique for Block-
Transform compressed image using Wavelet Transform Modulus
Maxima, IEEE Trans. On Image Processing, Vol 7, No. 10, October 1998
[8] Knutson, H. Filtering and Reconstruction in Image Processing,
73

[9] M.G. Mostafa, T.F. Gharib, Medical Image Segmentation Using a
Wavelet-Based Multiresolution EM Algorithm IEEE International
Conference on Industrial Electronics Technology & Automation, Dec.
2001.
[10] S. G. Mallat. A Wavelet Tour of Signal Processing, Third Edition,
[11] W. Doyle, Operations useful for similarity invariant pattern
recognition,JACM, v. 9, no. 2, April 1962, 259-267.
[12] S.G.Change,B.Yu and M.vetterli, Adaptive Wavelet Thresholding For
Image Denoising And Compression IEEE Trans. Image Processing, vol.9
, no pp. 1532-1545 ,SEP.2000
[13] C.S.Burrus R.A.Gopinath H.Guo Introduction To Wavelet And
Wavelet Transforms". NewJeresy: Pretice Hall,1998.
[14] N. Bonnet, J. Cutrona and M.Herbin. A no-threshold histogram-
based image segmentation method IEEE Trans. Pattern Recognition,
Volume 35, Issue 10, October 2002, pages 2319-2322.
[15] Shi Na, Liu Xumin and Guan yong. Research on k-means Clustering
Algorithm,An Improved k-means Clustering Algorithm. IEEE Third
International Symposium on Intelligent Information Technology and
Security Informatics, Capital Normal University
[16] A. Proch azka A. Gavlasov a. Wavelet Transform in Classification of
Biomedical Image. In IFMBE European Conference on Biomedical
Engineering, 2005.
74

[17] Lei Zhang, Paul Bao, Edge detection by scale multiplication in
wavelet domain, Pattern Recognition Letters 23, pp. 1771-1784, 2002.
[18] W.K. Pratt, Digital Image Processing, Wiley-lnterscience, New York,
1991.
[19] Gilbert G. Walter, Xiaoping Shen, Wavelets and Other Orthogonal
Systems Chapman & Hall 2001.
75