WAVELET BASED IMAGE SEGMENTATION
by
Haleh Zarrini
B.S., University of Urmia, 1991
A thesis submitted to the
University of Colorado Denver
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 noisefree representation of the
edges of the original image, which is provided as an input image to Kmeans
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
Kmeans 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
TABLE OF CONTENTS
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 Kmeans Algorithms..................................................28
3.2 Image Clustering...................................................35
vii
4. Segmentation and its Method
41
5. Conclusion.......................................................53
Appendix
A. Kmeans Clustering Algorithm MATLAB Program......................55
B. Kmeans Clustering Algorithm according to the color property.....57
C. Kmeans 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 Timefrequency boxes representing the energy spread of Gabor toms ... 6
2.2 Morelet wavelet function and Mexicanhat 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 TwoDimensional Wavelet Transform...15
2.6 One (left) and Two (right) Levels TwoDimensional DWT.............16
2.7 The scaling function..............................................18
2.8 The quadratic spline wavelet......................................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 kmeans algorithm................................30
3.3 The Kmeans algorithm on a grayscale image...................... 37
3.4 Color expression by L*a*b concept.................................38
3.5 The Kmeans 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 3D 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 3D 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 Kmeans 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 Kmeans 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 Kmeans algorithm. The
2
Kmeans algorithm is a numerical, unsupervised, nondeterministic,
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 Kmeans algorithm is applied on
grayscale 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, nonstationary 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
timeseries 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 timefrequency plane. To
measure timefrequency 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 timefrequency 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 Timefrequency boxes (Heisenberg rectangles) representing the
energy spread of Gabor atoms.
With wavelets, we usually dont speak about timefrequency but of timescale
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 nonstationary 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) <23>
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 Mexicanhat Wavelet.
7
2.2 Continuous Wavelet Transforms
Continuous Wavelet Transforms (CWT), transforms a continuoustime 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) (24)
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 = iC 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 timescale 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) (28)
with
= s0J/2 xp(s0Jtku0) (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 subsignals 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. Vco = {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,(plin) 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) = 11 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] (220)
Figure 2.4 shows graphical presentation of the Haar scaling function and
wavelet function [2, 19].
1
TI I I1 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 twodimensional 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
downsampled by two (indicated by the downarrow).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 TwoDimensional Wavelet
Transform.
15
Figure 2.6 shows wavelet representation of image. It illustrates decomposition
level of a twodimensional 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 TwoDimensional 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 twodimensional 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 multiscale
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 grayscale. 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) = (<) ard 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 nonorthogonal 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 (Vj
(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 subimages 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 Kmeans Algorithm
There are many clustering algorithms. One of the well known methods that
could be used for clustering is Kmeans 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 kmeans algorithm. Kmeans 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 Kmeans 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 Kmeans 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 kmeans 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 kmeans 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 kmeans 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 grayscale or black and white image
is also called bilevel or monochromic. The grayscale value 0 corresponds to
black and a value of 255 corresponds to white. Such an image is called an 8
bit grayscale image with
where G = {0, 1, 2,....., 255} that can be constructed from the 8bit, 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 Kmeans 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 Kmeans algorithm is used and the clustering is
performed by selecting 3 initial centers for k. These 3 centers were initially
selected randomly and as Kmeans algorithm is performed, the clusters start
to take shape and the centers are replaced after each iteration. Kmeans
clustering is done when there are no more changes in the cluster centers.
Figure 3.3 illustrates the original image and the result of Kmeans 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 Kmeans algorithm on a grayscale image.
37
Now, the kmeans 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 chromaticitylayer (it is
an objective specification of the quality of a color regardless of its luminance);
'a*' indicating where color falls along the redgreen axis and chromaticity
layer; and 'b*' indicating where the color falls along the blueyellow axis. All of
the color information is in the 'a*' and 'b* layers.
3. Apply Kmeans clustering;
Use Kmeans 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 Kmeans 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
Kmeans 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 clusterindex. 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 Kmeans 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 grayscale 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 grayscale 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
600011
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 Kmeans clustering
was selected because the Kmeans method is numerical, unsupervised,
iterative and nondeterministic. This method performs well on image and
texture inputs. The disadvantage with this method is presetting 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 Kmeans 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 xdirection.
(c) component in the ydirection.
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 xdirection, and the
difference between the third and first columns approaches the derivative in
the ydirection. 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 xdirection.
(c) \Gy\ Component in the ydirection. (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 grayscale 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 zcoordinates of
points above a grid in the xy 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 3D 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 Kmeans 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 Kmeans algorithm. The end result of Kmeans algorithm is
the segmented image. This process is represented in the Block Diagram
below:
Image
Modulus
Maxima
Edge
Detection
Select the Image
thont Noise
KMeans 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 Kmeans 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 3D 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 nonstationary 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
noisefree representation of the edges of the original image, which is provided
as an input image to Kmeans 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 Kmeans image segmentation algorithm.
In order to obtain the wavelet transform of the image in a suitable form to get
this transform representation by the modulusangle 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
Kmeans Algorithm MATLAB Program
%The following program clustered a grayscale image which is used house.bmp in this
%process. The Algorithm has been described in chapter 3.
clc
clear all
close all
%..............................
%% Stepl: Read Image
%..............................
he = imread('house.bmp');
%..........................................................
%% 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 kmeans
%..............................
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 Kmeans 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
Kmeans Algorithm on Color Image
%This program has been run on colored image by using kmeans algorithm for segmentation.
clc
clear all
close all
%............................
%% Stepl: Read Image
%............................
he = imread('cell3.jpg');
%.....................................................
%% Step 2: Convert the image to L*a*b* color space
%.....................................................
cform = makecform('srgb2lab');
lab_he = applycform(he, cform);
%...........................................................................
%% Step 3: Use Kmeans 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 Kmeans
%...........................................................................
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
%[X,map] = imread(cameramen.bmp);
[X,map] = imread('house.bmp);
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(k1));
%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(k1));
%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;
from = radian;
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
hd=[dyadup(hd,2) 0]; %a trous
hprime=hd;
gd=[dyadup(gd,2) 0]; %a trous
gprime=[dyadup(gprime,2) 0];
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 Kmeans 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);
%.................................................
% read the image
%.................................................
[X,map] = imread('house.bmp');
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(k1))';
%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(k1));
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
%..........................
hd=[dyadup(hd,2) 0]; %a trous
hprime=hd;
gd=[dyadup(gd,2) 0]; %a trous
gprime=[dyadup(gprime,2) 0];
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 denoising.
alpha = 1.2;
thr = wbmpen(c,1 .sigma,alpha);
% Use wdencmp for denoising 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 DeNoising Algorithm')
%subplot(1,2,2),image(wcodemat(Gd,nbc)); title('DeNoised Image')
%...............................
%% Use kmeans 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
%read image
f = imread('lena.bmp');
%f = imread('man.bmp');
%f = imread('cell.bmp')
%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
%read image
I = imread('house.bmp');
[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
I = imread('house.bmp');
[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);
%..........................................................................................
%Sobel gradient image(horizontal and vertical)
clc
clear all
close all
71
%read image
I = imread('house.bmp');
% 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, 198217, 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 2425, 1999), 357, pp 25432560.
[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,
Linkoping University, PhD Thesis 1982.
73
[9] M.G. Mostafa, T.F. Gharib, Medical Image Segmentation Using a
WaveletBased 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,
Academic Press, OrlandoSanDiego, 2001.
[11] W. Doyle, Operations useful for similarity invariant pattern
recognition,JACM, v. 9, no. 2, April 1962, 259267.
[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. 15321545 ,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 nothreshold histogram
based image segmentation method IEEE Trans. Pattern Recognition,
Volume 35, Issue 10, October 2002, pages 23192322.
[15] Shi Na, Liu Xumin and Guan yong. Research on kmeans Clustering
Algorithm,An Improved kmeans 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. 17711784, 2002.
[18] W.K. Pratt, Digital Image Processing, Wileylnterscience, New York,
1991.
[19] Gilbert G. Walter, Xiaoping Shen, Wavelets and Other Orthogonal
Systems Chapman & Hall 2001.
75
