Citation
Wavelet transform modulus maxima algorithms for medical image processing

Material Information

Title:
Wavelet transform modulus maxima algorithms for medical image processing
Creator:
Jagannath, Saritha
Place of Publication:
Denver, Colo.
Publisher:
University of Colorado Denver
Publication Date:
Language:
English
Physical Description:
x, 67 leaves : illustrations ; 28 cm

Thesis/Dissertation Information

Degree:
Master's ( Master of Science)
Degree Grantor:
University of Colorado Denver
Degree Divisions:
Department of Electrical Engineering, CU Denver
Degree Disciplines:
Electrical Engineering
Committee Chair:
Bialasiewicz, Jan T.
Committee Members:
Radenkovic, Miloje S.
Fardi, Hamid

Subjects

Subjects / Keywords:
Wavelets (Mathematics) ( lcsh )
Diagnostic imaging -- Mathematics ( lcsh )
Imaging systems in medicine -- Mathematics ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaf 67.).
Thesis:
Electrical engineering
General Note:
Department of Electrical Engineering
Statement of Responsibility:
by Saritha Jagannath.

Record Information

Source Institution:
|University of Colorado Denver
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
60423378 ( OCLC )
ocm60423378
Classification:
LD1190.E54 2004m J33 ( lcc )

Downloads

This item has the following downloads:


Full Text
WAVELET TRANSFORM MODULUS MAXIMA ALGORITHMS FOR
MEDICAL IMAGE PROCESSING
Saritha Jagannath
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Electrical Engineering
2004
by


This thesis for the Master of Science
degree by
Saritha Jagannath
has been approved
by
Miloje S. Radenkovic
Date
Hamid Fardi


Jagannath, Saritha (M.S., Electrical Engineering)
Wavelet Transform Modulus Maxima Algorithms for
Medical Image processing
Thesis directed by Professor Jan Bialasiewicz
ABSTRACT
Wavelets are powerful mathematical tools which represent the signal as a
combination of different frequency components, and then study each
component with a resolution matched to its scale. They have an edge over
traditional Fourier methods in analyzing discontinuities and sharp spikes in a
signal. The analysis of signals and images at different scales is an attractive
framework for many problems in signal and image processing. Wavelet
theory provides a mathematically precise understanding of the concept of
multiresolution. The focus of this thesis is on application of wavelets in the
field of medical image processing. Wavelets are used as a tool to analyze
CT scan images. Wavelet analysis is carried out on healthy human tissue
and malignant tissue in an attempt to clearly bring out the differences in the
behaviour between the two. The filter bank algorithm called algorithme a
trous is used for the above analysis. The analysis is carried out by detecting
the edges in the images and looking for a pattern difference between the
malignant tissue and the healthy tissue going across frequency bands.
Horizontal and vertical details are calculated in different frequency bands


and the corresponding modulus maxima images are found. Also Wavelet
packet analysis of healthy image and malignant image is carried out in
different frequency bands and the corresponding edges are found.
This abstract accurately represents the content of the candidates thesis,
recommend its publication.
Signed
Jan T. Bialasiewicz
IV


DEDICATION
I dedicate this thesis to my parents. If not for them I would not have
been where I am today. My deepest gratitude goes to them for
believing that I could actually come this far and get this degree. I truly
thank them for all their support and blessings.


ACKNOWLEDGEMENT
My sincere thanks to my advisor, Jan Bialasiewicz, for introducing me
to the world of wavelets and providing me with an opportunity to
explore the world of wavelets. Also, I wish to thank the University of
Colorado Health Science Center for providing me with kind of material
and support that was required for this thesis.


CONTENTS
Figures..............................................................viii
Tables..................................................................x
Chapter
1. Introduction to Wavelet Transforms...................................1
2. Discrete Dyadic Wavelet Transform....................................8
2.1 Dyadic Wavelet Transform Computed with Algorithme a trous.......11
3. Analysis of Images in Different Frequency Bands....................20
3.1 Medical Image Analysis in Different Frequency Bands................23
4. Introduction to Wavelet Packet Analysis.............................39
4.1 Wavelet Packet Analysis of Medical Images..........................42
5. Conclusion..........................................................48
Appendix
A. MATLAB function to read in the test image..........................50
B. MATLAB function to generate low pass and high pass
decomposition filters..............................................51
C. Algorithme a trous MATLAB computer program.........................52
D. MATLAB function to generate modulus maxima masks
for detail images..................................................59
E. MATLAB function to compute circular convolution....................65
F. MATLAB function to perform left shift of an input vector...........66
References.............................................................67
vii


FIGURES
1.1 Time-frequency boxes (Heisenberg rectangles) representing
the energy spread of Gabor atoms.................................2
1.2 Daubechies, db4, and Symmlet, sym8, Mother Wavelets.............4
1.3 Scaling Functions of the Daubechies, db4, and Symmlet,
sym8 Mother Wavelets..............................................4
1.4 Time-frequency boxes of two wavelets y/u s and ( ...............6
2.1 Filter bank for two-dimensional image decomposition
implementing the Mallat multiresoluton algorithm.................10
2.2 The Spline Wavelet (left) and its Scaling Function (right)......13
2.3 Filter bank for two-dimensional image decomposition
implementing the algorithme a trous............................14
2.4 Filter bank implementing image reconstruction with the
algorithme a trous.............................................16
2.5 Use of phase to determine the Modulus Maxima....................18
3.1 Cross section of neck image with tumor area marked..............25
3.2 Modulus maxima image for scale 2................................26
3.3 Modulus maxima image for scale 4...............................27
3.4 Modulus maxima image for scale 8...............................28
3.5 (a) Cross-section of malignant tissue with tumor area marked,
(b) Modulus maxima image for scale 2, (c) Modulus maxima
image for scale 4, (d) Modulus maxima image for scale 8..........30
3.6 (a) Cross-section of healthy neck, (b) Modulus maxima image
at scale 2, (c) Modulus maxima image at scale 4,
(d) Modulus maxima image at scale 8..............................31
3.7 Modulus maxima image at scale 2'5..............................33
3.8 Modulus maxima image at scale 2'4...............................34
3.9 Modulus maxima image at scale 2'3...............................35
viii


3.10 (a) Cross-section of malignant tissue with tumor area marked,
(b) Modulus maxima image for scale 2'5,(c) Modulus maxima
image for scale 2'4, (d) Modulus maxima image for scale 2'3.....36
3.11 (a) Cross-section of healthy tissue, (b) Modulus maxima image
for scale 2'5, (c) Modulus maxima image for scale 24,
(d) Modulus maxima image for scale 2'3..........................37
4.1 Binary tree of wavelet packet spaces...........................41
4.2 (a) Tree structure depicting three levels of details (b) Modulus
maxima image at nodeVK,' at scale 2, (c) Modulus maxima
image at nodeW23 at scale 4, (d) Modulus maxima image at nodeW37
at scale 8....................................................44
4.3 (a)Tree structure depicting three levels of details of healthy
tissue (b) Modulus maxima image at nodeW,' at scale 2,
(c) Modulus maxima image at node W23 at scale 4,
(d) Modulus maxima image at node W37 at scale 8...................45
4.4 (a) Tree structure depicting three levels of details of healthy
tissue (b) Modulus maxima image at nodeVT,' at scale 25,
(c) Modulus maxima image at node W/,23 at scale 2'4,
(d) Modulus maxima image at node W37 at scale 23...............46
4.5 (a) Tree structure depicting three levels of details of healthy
tissue (b) Modulus maxima image at nodeW,1 at scale 2'5,
(c) Modulus maxima image at node VT23 at scale 2'4,
(d) Modulus maxima image at nodeW37 at scale 2'3................47
IX


TABLES
3.1 Spline filter coefficients at scale 1
22
x


1. Introduction to Wavelet Transforms
Wavelet transforms are an important class of mathematical tools used in
medical signal processing as they help in unveiling hidden information.
Wavelet transform is capable of providing the signal content as a
function of time and frequency. It can be considered as an extension of
the Fourier transform. The major drawback of Fourier transform is that it
has only frequency resolution and no time resolution which means that
although we might be able to determine all the frequencies present in a
signal, we do not know when or where they are present. The Windowed
Fourier transform overcomes this drawback to a certain extent in that it
replaces the sinusoidal wave of Fourier transforms by a sinusoid and a
window which is localized in time however has frequency resolution
problems. Low frequencies can hardly be depicted by short windows and
high frequencies can only be poorly localized by long windows. The
windowed Fourier transform, defined by Gabor, correlates a signal fwith
each atom gui as represented by the following equation:
Sf(u,%)= Jf{t)g[4(t)dt= \f(t)g{t-uY^dt (1.1)
where u represents the time instant and £ represents frequency.


It is a Fourier integral that is localized in the neighborhood of u by the
window g(t-u). The windowed Fourier transform can also be written as a
frequency integral by applying the Fourier Parseval formula

(1.2)
The transform Sf{u,£) thus depends only on the values f(t) and f(a>)in
the time and frequency neighborhoods where the energies of gu,and
gu4 are concentrated (Mallat, 2001). Figure 1.1 shows the time-
frequency boxes of g^andg representing the energy spread of
Gabor atoms.
a
lg (0 I
Figure 1.1 Time-frequency boxes (Heisenberg rectangles)
representing the energy spread of Gabor atoms.
2


The limitations in resolution were the main reasons for invention of
wavelet transforms. The wavelet transform decomposes signals over
dilated and translated windows. The original window centered around
t= 0 and for which the scale s= 1, is called a mother wavelet, which is a
function y/e L2(R) with zero average, i.e.,
OO
jy/(t)dt=0 (1.3)
We obtain a family of time-frequency atoms by scaling y/ by s and
translating it by u(Mallat, 2001), i.e.,
/ \ 1 ft -u^
Vl,At) = -rv\
y/s
(1.4)
V J
The scale s indirectly corresponds to frequency. Larger scale means
dilated wavelet and lower frequencies. Smaller scale means compressed
wavelet and thus higher frequencies. The factor l/Vs is a normalization
factor to ensure that all wavelets have the same energy and this energy
is unity if |^(r| = 1 (Mallat, 2001).
Figure 1.2 shows an example of two different mother wavelets:
Daubechies (db4) and Symmlet (sym8).
3


A I I fv-N
l I II 1
Figure 1.2 Daubechies, db4, and Symmlet, sym8, Mother Wavelets
5 10 15
Symmlel Scaling Function
Figure 1.3 Scaling Functions of the Daubechies, db4, and Symmlet,
sym8, Mother Wavelets
The continuous wavelet transform of a signal f(t) is given by the following
scalar product:

t-u
dt
(1.5)
4


Like windowed Fourier transform, wavelet transform can measure the
time frequency variations of spectral components but it has a different
time-frequency resolution. Wavelets transform correlates f with^i( (. By
applying the Fourier Parseval formula, it can also be returned as
frequency integration:
In time, y/u 5is centered at L/with a spread proportional to s.
Figure 1.4 shows the time-frequency boxes of y/u 5 and y/ When the
scale s decreases, the time support is reduced but the frequency spread
increases and covers an interval that is shifted towards high frequencies
(Mallat, 2001).
(1.6)
5


Figure 1.4 Time-frequency boxes of two wavelets \f/u s and yr
We obtain the continuous wavelet transform of a signal by scaling the
mother wavelet, shifting the wavelet in time, multiplying with the signal
and integrating over all times. The obtained coefficients correspond to
the details of the signal since wavelet is considered as a high pass filter.
To recover f we also need a complement of information corresponding to
Wf(M,s) that is obtained using the scaling function 0. The scaling
function is an aggregation of wavelets at scales larger than 1 (Mallat,
2001). The scaling function 0 can be scaled at any scale s. It is
interpreted as a low pass filter and is given by the equation 1.7.
f-1
Uv
(1.7)
6


The low-frequency approximation of fat scale s is
t u
Lf(u,s )=(/(?),^ ^
(1.8)
If the low frequency approximation was obtained at s=s0, then the
original signal can be recovered using the following equation:
f(t)=^ ]w(,s)* Vs(t)f+^-z/U)* K M
(1.9)
where C = ~^^-dto (1-10)
o >
The first term in equation (1.9) represents the part of f{f) recovered from
its high frequency representation for scales s < s0 and the second term
represents the low frequency part of f(t). The condition:
c = ]MH2rf(y<
o to
is called the wavelet admissibility condition for a real function,
^(/)e L2(r), to be a wavelet (Mallat, 2001).
7


2. Discrete Dyadic Wavelet Transform
In the computer world, an image is a collection of pixels at a certain
resolution. Thus continuous wavelet transform cannot be computed at
arbitrary small scales. To construct a discrete wavelet from a continuous
wavelet, the scale s is discretized but not the translation parameter u.
For fast computations the scale is sampled along a dyadic
sequence The scale parameter!7 is the inverse of resolution! 1.
The dyadic wavelet transform of a signal f e L2(R) is defined by
In computer algorithms, an image is analyzed using discrete filters.
Discrete high pass and low pass filters, called the conjugate mirror
filters, replace continuous wavelets and scaling functions, respectively.
The filters are related to their continuous counterparts as follows:
(2.1)
(2.2)
(2.3)
8


The high pass filter §[/?] is related to the low pass filter h[n] as follows:
g[n]=(-l)' /z[l-n] (2-4)
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 impulse response of a
filter. The dyadic wavelet transform (DWT) analyzes the image at
different frequency bands with different resolutions by decomposing it
into a coarse approximation and detail information with high-pass and
low-pass filters abbreviated as HP and LP, respectively. 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-dimensional
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.
9




v+1
-+d
7+1
'^7+1
Figure 2.1 Filter bank for two-dimensional image decomposition
implementing the Mallat multiresoluton algorithm.
The multiresolution analysis involves a decomposition of the image
space into a sequence of subspaces Vy. More formally, the
approximation of a function at a resolution TJ is defined as an
orthogonal projection on a space V. projection of a function f is defined by a function fs that minimizes the
distancej/-/. ||.
10


As illustrated in Figure 2.1, the following are the steps of the algorithm
(known as Mallats multiresolution algorithm):
First, the rows of tf^ are passed through an HP and LP filter. This
generates two images, which are then down-sampled by two
(indicated by the down-arrow), i.e., 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. This
generates four images: approximation image aj+], the horizontal
detail dhj+], the vertical detail d'j+], and the diagonal detail ddj+v
Typically, quadrature mirror filters such as Daubechies filters are used
for perfect reconstruction. These have a special relationship between the
LP and HP decomposition filters.
2.1 Dyadic Wavelet Transform Computed with
Algorithme a trous
Dyadic wavelet transform can also be computed with a filter bank
algorithm called in French the algorithme a trous. The algorithme a
trous uses quadratic spline filters properly related to the quadratic spline
11


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. A box spline of degree m is a
translation of m+1 convolutions of 1 {0,i] with itself. It is centered at t = 1/2 if
m is even and at t = 0 if m is odd. Its Fourier transform is
<*
e
CO
l£
2
sin
_____2_
co
(2.5)
V 2 ;
with e = 1 if m is even and £ = 0 if m is odd.
From equation 2.2, the relation between the low pass filter and the
scaling function is frequency domain is as follows:
h(co) = V2
(2 co)
(2.6)
where h(co) is the Fourier transform of h[n].
We can construct a wavelet that has one vanishing moment by
choosing g(co) = O{co) in the neighborhood of co = 0. For example
.CO
g(co) = -iyfle sin~ (2-7)
12


The Fourier transform of the resulting wavelet is
W(co) = -j=g
-1(0
-iO){ 1 + f)
. co
sin
____4
co
4* ,
. m-t2
(2.8)
(Mallat, 2001)
The mother wavelet and the scaling function of quadratic spline are
shown in Figure 2.2.
Figure 2.2. The Spline Wavelet (left) and its Scaling Function (right)
In the implementation of this algorithm for image decomposition, all
filtering is done with circular convolutions to prevent the discontinuity
induced by finite-length signals near image boundaries. After each level j
of decomposition, the filters are dilated by inserting 2;-1 zeros between
13


consecutive filter coefficients. This creates holes (trous in French).
There is no down sampling only dilation of the filters.
columns rows
d
j+1
d
h
7+1
Figure 2.3 Filter bank for two-dimensional image decomposition
implementing the algorithme a trous.
One of the decomposition levels of this algorithm is illustrated in Figure
2.3. The approximations are generated by first shifting the columns of
a,j\o the left 2J~l times. Then, the columns are filtered using the LP filter
h[n]. Next, the rows are shifted to the left 2y'1 times and filtered with an
LP filter. The result is the approximation aj+] of the image. The
horizontal details are calculated by first filtering the columns of fusing
the HP filter g[] and then shifting the columns 2 times. The vertical
14


details are calculated by first filtering the rows of fusing the HP filter
and then shifting the rows 2y times.
The following equations used in the MATLAB program give the
convolution formulas that are cascaded to calculate the discrete dyadic
wavelet transform and its inverse:
DWT- aj+][n] = cij hjhj[n] (2-9)
dx;+][n]=aj*gjM (2-10)
dj+in] = aj*^j[n] (2'11)
ajln] = ^(aj+\ *hjhj[n}+d'j+l *gj$n]+dhj+i *<£,[]) (2.12)
Note that Uj = dilated analysis filter and = decimated synthesis
filter yr S[n] is discrete Dirac and xy[n,m]=x[n] y[m] is a separable filter
(Mallat, 2001).
Reconstruction of the image is followed in reverse order from the
decomposition procedure with the exception of shifting that is altered.
After each level of reconstruction, the filters are decimated by inserting
15


2y -1 zeros between the consecutive coefficients. Note that j decreases
at each level of reconstruction.
The synthesis filters HP g[n] and LP h[n]are used and the resultant
approximation and detail images are added together and multiplied by
1/4. The filter bank implementing this reconstruction process is shown in
Figure 2.4.
row
columns
d
h
j+1
d
v
j+
h[n] h[n]
Cl;
Figure 2.4. Filter bank implementing image reconstruction with the
algorithme a trous.
At every decomposition level and at any pixel of the image, we have
available horizontal and vertical detail calculated using the algorithme a
16


trous. These are orthogonal components of the wavelet transform.
Using these components, we calculate the wavelet transform modulus
(WTM) and angle using the following equations:
Mjf(x,y) = J(dhj(x,y))2 + d](x, y))2 (2'13)
Ajf(x,y) = tan~'(dj(x,y)/dj(x,y)) <2-14)
For each pixel location within each scale, a modulus maxima is detected
at location (*,>) if M }f{x,y)> My/(x,,>,) andMjf{x,y)> Mjf{x2,y2).
Locations and (x2,y2) are adjacent locations to (x,y) indicated by
the angle A./(x, y). The angle can be one of eight quantized directions:
0, tt/4, tt/2, 3tt/4, tt, tt/4, -tt/4, -tt/2, or -3tt/4 (Hsung, Lun, & Siu, 1998).
17


Xi/Vi
*2,Y2
Figure 2.5 Use of phase to determine the Modulus Maxima.
Figure 2.5 above is a graphical description of how the angle is used to
locate the two points in which to compare a particular wavelet modulus.
For example, if the angle that corresponds to point (*,>) lies between
the tt/8 and 3tt/8, indicated by dashed lines, the quantized angle
becomes tt/4. The line tangent to this angle, shown by a short line,
correspond to the direction in which the two adjacent positions to be
compared are located. In this example, the wavelet transform modulus
located at (*,>>), indicated by the black square, is compared to the
wavelet transform modulus located at (*,,>,) and(^2,>2), indicated by
18


the two shaded squares. Positive angle values that correspond to the
negative angle values are included in the figure because the MATLAB
programs in Appendix C and Appendix D work with positive angle values
only.
19


3. Analysis of Images in Different Frequency
Bands
Wavelet decomposition provides natural setting for the multi-level edge
detection. Since Wavelet Transform Modulus Maxima (WTMM) provide
t
useful information for pattern recognition, we propose to use it here for
fast feature extraction.
It is possible to analyze an image in different frequency bands. This is
accomplished by assuming that the original image is initially at some
scale V where je Z. With this assumption, the image is decomposed
by going higher in scale from the starting scale. The MATLAB function
scaling (), attached in Appendix B, scales the mother wavelet and
scaling function of quadratic splines. From the scaled functions the
corresponding HP and LP filters are obtained. Using these filters the
image is analyzed in the corresponding frequency band. The scaling of
the wavelet and scaling function is done in the frequency domain for
easy computation purposes. The analysis is carried out using a
threshold, which produces the best visual results.
20


To design dual scaling functions 0 and wavelets ip, which are splines,
we choose h=h. As a consequence, 0 =

condition implies that
g(a>) =
2-
h{(0)
2 \
g (CO)
co
= ~iy[I/'*
(0 \ '
e sin y j
2 n=0
CO
cos
2
2/i
(3-1)
The scaling and wavelet function for the quadratic spline, defined in
frequency domain, are
0{co) =
sin(ry/2)V
(012 J
exp
f -ico^
\ z J
(3-2)
ip(co) =
-1(0
\
sin(ry/4)
V
co! 4
exp
( -ico'
v Tj
(3.3)
The HP and LP decomposition filters are calculated as follows:
( n
h(2jO)) =
0(2 j+'co)
V27 I (3.4)
g(2ja>) =
V2
7+1
A -
\p{l;+l co)
V27 I 0(2 Jco)
(3.5)
where jeZ (Mallat, 2001)
21


Table 3.1 gives the coefficients of quadratic spline filters at scale 1.
n h[n] h[n] g[n] g[n]
-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
Table 3.1 Spline filter coefficients at scale 1.
In this table,
h[n] low pass decomposition filter.
h[n] low pass reconstruction filter.
g[n] high pass decomposition filter.
g[n] high pass reconstruction filter.
The low pass filters h[n] = h[n] are obtained by Inverse Fourier
Transform of equation (3.4) and the high pass filters g[n] and g[n] are
obtained by Inverse Fourier Transform of equations (3.5) and (3.1).
22


3.1 Medical Image Analysis in Different Frequency
Bands
In this section, the above theory of analyzing an image in different
frequency bands is applied to Computed Tomography (CT) scan images.
One of the advantages of CT over standard x-rays is that it clearly shows
soft tissue such as the brain, as well as dense tissue like bone. Wavelet
analysis is carried out at three levels going higher in scale and modulus
maxima image is found at each level. This analysis is carried out in
different frequency bands by assuming that the original image is at some
scale V. First the image is read using the function cal_mm(j, Ivl, image)
attached in Appendix A where j corresponds to scale 2j, Ivl is number of
decomposition levels and image is the image under analysis.
The function 'cal_mm(j, Ivl, image) calls function 'scaling(j) attached in
Appendix B which calculates LP and HP decomposition filters
corresponding to scale 2j. Equations (3.2) and (3.3) are used in the
function to calculate the same. The number of decomposition levels and
image under analysis along with the above obtained filters are passed to
the function 'mm_atrous(lvl,y, h,g,gp) attached in Appendix C. This
function calculates the WTMM, Wavelet Transform Modulus Maxima
23


Masks, for every level. Figure 3.1 gives the CT scan image of cross
section of human neck with tumor area marked.
Modulus maxima are calculated for this image at different scales.
Modulus maxima images at scales 2, 4 and 8 are shown in Figure 3.2,
Figure 3.3, and Figure 3.4, respectively. For this case, the image shown
in Fig.3.1 is assumed to be at scale 1.
24


Figure 3.1. Cross section of neck image with tumor area marked.
25


Figure 3.2 Modulus maxima image for scale 2
26


Figure 3.3 Modulus maxima image for scale 4
27


Figure 3.4 Modulus maxima image for scale 8
28


The tumor area is marked in each case and it is observed that there is
an increase in edges in the tumor area when going from scale 2 to 4.
This is in contrast to tumor free areas where there is a decrease in the
number of edges as we go higher in scales. This is more clearly seen in
Fig.3.5 where the tumor area of the cross section is cut out and analyzed
with an assumption that the original image is at a scale 1. This pattern is
not seen in the healthy tissue. Here the numbers of edges reduce when
the scale is increased from 2 to 4. Figure 3.6 shows the analysis results
of the healthy tissue. In this case, we can see that going higher in scale
the number of edges uniformly decreases.
29


Figure 3.5 (a) Cross-section of malignant tissue with tumor area marked,
(b) Modulus maxima image for scale 2, (c) Modulus maxima image for
scale 4, (d) Modulus maxima image for scale 8.
30


Figure 3.6 (a) Cross-section of healthy neck, (b) Modulus maxima image
at scale 2, (c) Modulus maxima image at scale 4, (d) Modulus maxima
image at scale 8.
31


When analyzed in different frequency band, that is, by assuming that the
original image is at a scale 2j where j ? 0, the edges in the tumor area
reduce with scale unlike the case when j= 0.
Figures 3.7, 3.8 and 3.9 give the modulus maxima images corresponding
to Figure 3.1 at scales 2'5, 2'4 and 2'3 respectively. Smaller scales
correspond to higher frequencies. As we go higher in frequency, the
number of edges detected increase. There is denser pattern observed.
However, there is no difference in pattern between the healthy tissue and
malignant tissue.
Figure 3.10 and Figure 3.11 give the modulus maxima images of
malignant and healthy tissue cross-sections at scales 2'5, 24 and 2'3
32



Figure 3.7 Modulus maxima image at scale 2'5
33


Figure 3.8 Modulus maxima image at scale 2'4
34


Figure 3.9 Modulus maxima image at scale 2'3
35


Figure 3.10 (a) Cross-section of malignant tissue with tumor area marked,
(b) Modulus maxima image for scale 2'5, (c) Modulus maxima image for
scale 24, (d) Modulus maxima image for scale 2'3
36


Figure 3.11 (a) Cross-section of healthy tissue, (b) Modulus maxima
image for scale 2'5, (c) Modulus maxima image for scale 2'4, (d) Modulus
maxima image for scale 2'3.
37


From the above analysis it follows that by wavelet analysis of an image in
different frequency bands, it could be possible
difference in behavior between the healthy tissue
This can be used as a tool for detection of cancer.
to see some kind of
and malignant tissue.
38


4. Introduction to Wavelet Packet Analysis
A wavelet transform is a method of expressing a signal in terms of the
wavelet basis {y/^t)}. This accounts to decomposing the signal space Is
into a direct sum of orthogonal subspaces {V;} and {W;}. If the signal
bandwidth is finite with information upto a resolution level J, a wavelet
transform performs a decomposition of the space V, into a sum of
orthogonal subspaces
Vj =VJ.I=vo_,ew^ew,., = ...=£' v0 <4-1)
Of course this is not the only way to decompose the signal space Is or
V
Wavelet packets, introduced by Coifman, Meyer and Wickerhauser,
generalize the link between multiresolution approximations and wavelets.
(Mallat, 2001).
From the theory of multiresolution analysis, we know that given the basis
functions {0j(t-2jn)} of Vjt { orthonormal basis of Vj+] and WJ+], respectively, and Vj= v,.+lwy.+1.
The relations between these basis functions are represented by the
following equations:
39


(4.2)
/J=oo
Wj+,(t-V+'n) = 4l ^glnty^t-Vn) (4.3)
Where h and g are conjugate mirror filters and
g[/i] = (-1)1 "/j[l-n] (4-4)
So the space Vj can be decomposed into a sum of two orthogonal
subspaces defined by their basis functions given by the equations (4.2)
and (4.3). Using a pair of conjugate mirror filters to decompose a signal
space corresponds to splitting the frequency content of the signal into a
low-frequency and a high-frequency component. In orthogonal wavelet
decomposition procedure, a signal space is decomposed into
approximation space and details space. Next, the approximation space is
further decomposed into next level of approximation and details and so
on. The details space is never decomposed into coarser scales. In wavelet
packet analysis, each details space is also decomposed into two spaces
using the same approach as in decomposition of approximation space. So
in general, wavelet packet decomposition divides the frequency space into
various parts and allows better frequency localization of signals.
40


Instead of decomposing only the approximation spaces VJ into two
orthogonal subspaces, we can decompose the detail spaces W as well.
The recursive splitting of the vector spaces can be represented as a
binary tree shown in Fig.4.1.
iy
w,

w
w
w;
W-
Figure 4.1 Binary tree of wavelet packet spaces
Each node (y,p) corresponds to a space Wf, which admits an orthonormal
basis {y/^(t 2Jn)}. At the root of the tree, we have W0=V0. Analogously,
we can define two orthogonal bases at the children nodes as follows:
y#,(0 = (4.5)
£h[W(,-Vn) (4.6)
Thus far we have been using the combination of {^;(r-2;)}and
{y/j{t-2jn))\o form a basis forV^, and now we have a whole sequence of
41


functions at our disposal. Various combinations of these and their dilations
and translations can give rise to various bases for the function space. So
we have a whole collection of orthonormal bases generated.
4.1 Wavelet Packet Analysis of Medical Images
The transformation of data into wavelet packet basis presents no extra
numerical difficulties. We can simply do a convolution using filters h and g
on the details space as well as on the approximation space. As in the
wavelet transform, we can keep doing the decomposition until we cannot
go any further. On the other hand, we can also choose not to decompose
a particular subspace while decomposing others. So there are many
choices for decomposing a signal. We can keep all the coefficients at all
decomposition levels and generate a table of coefficients of wavelet
packet decomposition.
The MATLAB function wavelet_packet_analysis_of_detail_images(lvl,y,
h,g,gp), given in Appendix D, is used to achieve this. The function
cal_mm(j, Ivl, image) calls function scaling(j) attached in Appendix A
which calculates LP and HP decomposition filters corresponding to scale
2j. Equations (3.2) and (3.3) are used in the function to calculate the
same. The number of decomposition levels and image under analysis
along with the above obtained filters are passed to the function
42


wavelet_packet_analysis_of_detailJmages(lvl,y,h,g,gp)' attached in
Appendix D. This function calculates the WTMM masks for every level of
detail images. The modulus maxima images are obtained by initializing the
detail images at every level as the image to be analyzed for the next level.
Figure 4.1 gives the modulus maxima images of malignant tissue for three
levels of detail images. Figure 4.2 gives the modulus maxima images of
healthy tissue for three levels of detail images. The image is assumed to
be initially at scale 1. Figure 4.3 and Figure 4.4 give the analysis results
for smaller scales corresponding to high frequencies of malignant and
healthy tissue respectively. No difference in pattern between the healthy
tissue and malignant tissue was observed in WTMM masks obtained
using wavelet packet transform of detail images.
43


Figure 4.2 (a) Tree structure depicting three levels of details (b) Modulus
maxima image at nodeW,' at scale 2, (c) Modulus maxima image at
nodeW,3 at scale 4, (d) Modulus maxima image at nodeW,7 at scale 8.
44


Figure 4.3 (a) Tree structure depicting three levels of details of healthy
tissue (b) Modulus maxima image at nodeW,1 at scale 2, (c) Modulus
maxima image at node W? at scale 4, (d) Modulus maxima image at
nodeVK,7 at scale 8.
45


Figure 4.4 (a) Tree structure depicting three levels of details of healthy
tissue (b) Modulus maxima image at nodeW,' at scale 2'5, (c) Modulus
maxima image at node W? at scale 2'4, (d) Modulus maxima image at
node Wj at scale 2'3.
46


Figure 4.5 (a) Tree structure depicting three levels of details of healthy
tissue (b) Modulus maxima image at nodeVT,1 at scale 2'5, (c) Modulus
maxima image at node W* at scale 2'4, (d) Modulus maxima image at
nodeVT37 at scale 2'3.
47


5. Conclusion
An extensive analysis of medical images was carried out using the
Algorithme a trous and quadratic spline wavelets. The images were
analyzed by detecting edges and looking for a pattern difference between
the malignant tissue and the healthy tissue going across frequency bands.
Horizontal and vertical details were calculated in different frequency bands
and the corresponding modulus maxima images were found.
Difference in pattern between the healthy tissue and malignant tissue was
found while going from scale 1 to 2. In case of healthy tissue, the number
of edges decreased from scale 1 to 2 whereas in the malignant tissue
there was an increase in the number of edges. The images were analyzed
at both positive and negative scales by assuming that the original image is
at some scale {2'}y.eZ. Also wavelet packet analysis of healthy image and
malignant image was carried out in different frequency bands and the
corresponding modulus maxima images were found.
The algorithm developed in this thesis was able to detect the differences
in behavior between healthy and malignant tissues only in the modulus
maxima images obtained by wavelet analysis. These differences could not
48


be seen in modulus maxima images obtained by wavelet packet analysis.
Also the differences were seen only when going from scale 1 to 2.
In addition, it should be noted that representing a single CT scan image in
different frequency bands and showing the edges of the obtained images,
we are providing the information not available by looking at the original
image. However, the radiologists have to learn how to use this information
for medical diagnosis.
49


APPENDIX
Appendix A: MATLAB function to read the test image.
function cal_mm(scale,Ivl,image)
A=imread(image);
%Convert RGB image to grayscale
A=double(A);
Xrgb=0.2990*A (:,1)+0.5870*A(:,:,2)+0.1140*A(:,:,3);
NbColors=256;
X=wcodemat(Xrgb, NbColors);
map=gray(NbColors);
%Normalize image
X = X./256;
y = x;
save_orig = y; %Save original image to display later
[nr,nc]=size(y);
%Convert 2-D image into a 1-D vector
orig_ vector = y(1,:);
for i = 2:nr
orig_ vector = [orig_ vector y(i,:)];
end
%Decompose image
%a = last level approximation
%D1_MM = contains all levels of the horizontal details corresponding to
the modulus maxima matrix masks
%D2_MM = contains all levels of the vertical details corresponding to the
modulus maxima matrix masks
%gprime = reconstruction high pass filter
%hprime = reconstruction low pass filter
[h,g,gp]=scaling(scale);
[a, D1_MM, D2_MM, gprime, hprime] = mm_atrous(lvl,y,h,g,gp);
50


Appendix B: MATLAB function to generate low pass and high pass
decomposition filters
%This function is called by cal_mm(scale,lvl,image) to help generate the
low pass and high pass decomposition filters.
%This function generates the high pass and low pass decomposition
filters for particular frequency band.
%lnputs
%j where 2) gives the scale parameter
%Outputs
%h=low pass decomposition filter
%g=high pass decomposition filter
%gp=high pass reconstruction filter
function [g,h,gp]=scaling(j)
d=1;
for n=-2:1:3
syms w a b;
phi=(((sin(w/2))/(w/2))A3) exp(-i*w/2);
psi=(-i*w/4) *(((sin(w/4))/(w/4))A4) *exp(-i*w/2);
H=((sin(2Aj*w)/(2Aj*w))A3)*exp(-
i*2Aj*w)/(((sin(2Aj*w/2)/(2Aj*w/2))A3) *exp(-i*2Aj *w/2));
H=sqrt(2) *H*exp(i*w*n);
G=(-i*2Aj*w/2)*(((sin(2Aj*w/2))/(2Aj*w/2))A4)*exp(-
i*2Aj*w)/(((sin(2Aj*w/2)/(2Aj*w/2))A3) *exp(-i*2Aj *w/2));
G=sqrt(2) *G *exp(i*w*n);
GP=-i*sqrt(2) *exp(-
i*2Aj*w/2)*sin(2Aj*w/2)*(1+(cos(2Aj*w/2))A2+(cos(2Aj*w/2))A4)*exp(i*w*n);
a=-pi;
b=pi;
g(d)=( 1/(2*pi)) *INT(H, w,a,b);
h(d)=( 1/(2*pi)) *INT(G,w, a, b);
gp(d)=( 1/(2*pi)) *INT(GP, w,a, b);
d=d+1;
end
51


Appendix C: Algorithm a trous MATLAB computer program
%mm_atrous(lvl,y, h,g,gp)
%This function is called by cal_mm(scale,lvl,image) to help generate the
%modulus maxima matrix masks.
%This function generates the masked horizontal and vertical details along
%with the last approximation.
%.....................................
%inputs:
%.....................................
%lvl = number of levels of decomposition
%y = original image
%h=low pass decomposition filter
%g=high pass decomposition filter
%gp=high pass reconstruction filter
%..........................................
%Outputs:
%..........................................
%y = last level approximation image
%D1_MM = contains masked horizontal details from every level
%D2_MM = contains masked vertical details from every level
%gprime = reconstruction high pass filter
%hprime = reconsrtuction low pass filter
%Finally the program displays the modulus maxima matrix masks for each
%level.
function [y, D1_MM, D2_MM, gprime, hprime] = mm_atrous(lvl,y,h,g,gp)
%.....................
%Spline wavelet filters
%.....................
hf=double(h)
gf=double(g)
gprime=double(gp)
L=lvl; %L=level of detail
52


[nr,nc]=size(y);
nonzeros = L*2*nr*nc;
O/
/O-------------------------
%Loop for decomposing image
%...........................
for k = 1:L
O/
/o--------------------------------------------
%caicuiates approximations.
%shift either image or latest approximation vertically.
O/
/o--------------------------------------------
for i=1:nc
a{k}(1 :nr, i) = leftshift(y( 1 :nr, i) ',2A(k-1))
O/
/o-----------------------------------------------
%circular convolution of hf with columns from vertically shifted image
O/
/o-----------------------------------------------
a{k}(1:nr,i) = cconv(hf,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));
%..........................................-........
%circular convolution of hf with rows from horizontally shifted image
O/
/o--------------------------------------------------
a{k}(i, 1:nc) = cconv(hf,a{k}(i, 1:nc));
end;
53


0/
/o------------------------------------------
%calculates horizontal and vertical details.
%...........................................
%Circular convolution of high pass filter, gf, with columns from image
or...
%... previous approximation
O/
/o------------* ------------------- -----
for i=1:nc
d 1 {k}( 1 :nr, i)=cconv(gf,y( 1 :nr, i)')
%...............................-.....................
%Shift resultant image vertically. The result is a horizontal detail...
%... image
O/
/o----------------------------------------------------
d1{k}(1 :nr, i)=leftshift(d 1{k}(1 :nr, i) ,2Ak)
end;
%
%Circular convolution of high pass filter, gf, with rows from image or...
%... previous approximation
%
for i=1:nr
d2{k}(i, 1 :nc)=cconv(gf,y(i, 1:nc));
%.................................................
%Shift resultant image horizontally. The result is a vertical detail...
%... image
O/
/o-------------------------*---------------------
d2{k}(i, 1 :nc)=leftshift(d2{k}(i, 1:nc),2Ak);
end;
d{k}=d1{k}+d2{k};
%....-...-
%Dilate filters
%...........
54


hf=[dyadup(hf,2) 0]; %a trous
hprime=hf;
gf=[dyadup(gf,2) 0]; %a trous
gprime = [dyadup(gprime,2) 0];
y=a{k};
end;
%..............................
%Set threshold to zero
Threshold = 0.005;
%......................................
%Calculate the modulus maxima for each level
0/
/o---------------------------------
for k = L:-1:1
wtm{k} = sqrt(d1{k}.A2 + d2{k}.A2); %Wavelet Transform Modulus
[modern,mod_n] = size(wtm{k});
%.......................................
%calculate angle of the wavelet transform vector.
%.......................................
forqt= 1:nr
forr= 1:nc
alpha{k}(qt,r) = atan2(d2{k}(qt,r),d1{k}(qt,r));
if alpha{k}(qt,r) < 0
alpha{k}(qt,r)=2*pi+alpha{k}(qt,r);
end
end
end
%..............................
%calculate modulus maxima image
Re................-......-.....
for r=1 :mod_m
for c=1 :mod_n
mod_max{k}(r,c)=255; Reinitialize modulus maxima array
end
55


end
O/
/o--------------------------
%find local maximum of modulus
%.........-.................
ang=pi/8;
for r=2:mod_m-1
for c=2:mod_n-1
if (alpha{k}(r,c)>=(15*ang)l...
alpha{k}(r,c)<=(1*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c) & wtm{k}(r,c) >= wtm{k}(r+1,c)
mod_max{k}(r, c)=0;
elseif (alpha{k}(r,c)>=( 1 *ang)&...
alpha{k}(r,c)<=(3*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c-1) & wtm{k}(r,c) >= wtm{k}(r+1,c+1)
mod_max{k}(r,c)=0;
elseif (alpha{k}(r, c)>=(3*ang)&...
alpha{k}(r,c)<=(5*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r,c-1) & wtm{k}(r,c) >= wtm{k}(r,c+1)
mod_max{k}(r, c)=0;
elseif (alpha{k}(r,c)>=(5*ang)&...
alpha{k}(r, c)<=(7*ang))&...
wtm{k}(r,c)>=Threshold&...
wtm{k}(r,c)>=wtm{k}(r-1 ,c+1) & wtm{k}(r,c) >= wtm{k}(r-h1,c-1)
mod_max{k}(r, c)=0;
elseif (alpha{k}(r,c)>=(7*ang)&...
alpha{k}(r,c)<=(9*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c) & wtm{k}(r,c) >= wtm{k}(r+1,c)
mod_max{k}(r,c)=0;
elseif (alpha{k}(r,c)>=(9*ang)&...
56


alpha{k}(r,c)<=(11 *ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c-1) & wtm{k}(r,c) >= wtm{k}(r+1,c+1)
mod_max{k}(r, c)=0;
elseif (alpha{k}(r,c)>=( 11 *ang)&...
alpha{k}(r,c)<=(13*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r,c-1) & wtm{k}(r,c) >= wtm{k}(r,c+1)
mod_max{k}(r,c)=0;
elseif (alpha{k}(r,c)>=( 13*ang)&...
alpha{k}(r,c)<=(15*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c+1) & wtm{k}(r,c) >= wtm{k}(r+1,c-1)
mod_max{k}(r,c)=0;
end
end
end
end
%................................................
%Build vertical and horizontal detail matrices that include only those...
%...coefficients that are located at the modulus maxima; otherwise, the...
%...coefficients are set to zero.
O/
/O---------------------------------------*-------
for k = 1:L
d1_mm{k} = zeros(nr,nc);
d2_mm{k} = zeros(nr,nc);
fori = 1:nr
forj = 1:nc
if mod_max{k}(i,j) == 0
d1_mm{k}(i,j) = d1{k}(i,j);
d2_mm{k}(i,j) = d2{k}(i,j);
end
end
end
D1_MM(:,:,k) = d1_mm{k};
D2_MM(:,:,k) = d2_mm{k};
end
57


mm_nonzeros = nnz(D1_MM) + nnz(D2_MM);
percent = mm_nonzeros/nonzeros
O/
/o------------------------------------
%Display modulus maxima matrix mask
%.....................................
figure;
O/
/o------------------------------------
o/
/o------------------------------------
imshow(mod_max{ 1},[]);
xlabel('Modulus Maxima for Level 1');
figure;
%.....................................
%.....................................
imshow(mod_max{2},[]);
xlabel(Modulus Maxima for Level 2');
figure;
%.....................................
0/
/0-------------------------------------
imshow(mod_max{3},[]);
xlabel('Modulus Maxima for Level 3');
58


Appendix D: MATLAB function to generate modulus maxima masks for
details images
% wavelet_packet_analysis_of_detail_images (lvl,y,h,g,gp)
%This function is called by cal_mm(scale,lvl,image) to help generate the
%modulus maxima matrix masks for detail images.
%This function generates the masked horizontal and vertical details along
%with the last details
O/
/o--------------------------------------
%lnputs:
O/
/o--------------------------------------
%lvl = number of levels of decomposition
%y = original image
%h=low pass decomposition filter
%g=high pass decomposition filter
%gp=high pass reconstruction filter
%..........................................
%Outputs:
%..........................................
%y = last level approximation image
%D1_MM = contains masked horizontal details from every level
%D2_MM = contains masked vertical details from every level
%gprime = reconstruction high pass filter
%hprime = reconstruction low pass filter
%Finally the program displays the modulus maxima matrix masks for each
%level.
function [y, D1_MM, D2_MM, gprime, hprime] =
wavelet_packet_analysis_of_detail_images (lvl,y,h,g,gp)
%..........................................
%Spline wavelet filters
%..........................................
hf=double(h)
gf=double(g)
gprime=double(gp)
L=lvl; %L=level of detail
[nr,nc]=size(y);
nonzeros = L*2*nr*nc;
O/
/o----------------------------------------
59


%Loop for decomposing image
%........................
fork = 1:L
%...................................-..........
%calculates approximations.
%shift either image or latest approximation vertically.
O/
/O---------------------------------------------
for i=1:nc
a{k}(1:nr,i) = leftshift(y(1 :nr,i)',2A(k-1))';
O/
/O---------------------------------------------
%circular convolution of hf with columns from vertically shifted image
%..............................................
a{k}(1 :nr, i) = cconv(hf,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));
O/
/o---------------------------------------------
%circular convolution of hf with rows from horizontally shifted image
O/
/o-----------------------------------------
a{k}(i, 1:nc) = cconv(hf,a{k}(i, 1:nc));
end;
%.........................................-
%calculates horizontal and vertical details.
%..............................................
%Circular convolution of high pass filter, gf, with columns from image
or...
%... previous approximation
0/
/o---------------------------------------------
for i=1:nc
d 1 {k}( 1 :nr, i)=cconv(gf,y( 1 :nr, i)')
O/
/O---------------------- -------------------
%Shift resultant image vertically. The result is a horizontal detail...
%... image
%..............................................
d1{k}(1:nr, i)=leftshift(d 1 {k}(1:nr, i)\2Ak)
end;
60


o/
/o-------------------------------------............
%Circular convolution of high pass filter, gf, with rows from image or.
%... previous approximation
%................................................
for i=1:nr
d2{k}(i, 1 :nc)=cconv(gf,y(i, 1:nc));
%............................................
%Shift resultant image horizontally. The result is a vertical detail..
%... image
%............................................
d2{k}(i, 1 :nc)=leftshift(d2{k}(i, 1:nc),2Ak);
end;
d{k}=d1{k}+d2{k};
%............................................
%Diiate filters
%............................................
hf=[dyadup(hf,2) 0]; %a trous
hprime=hf;
gf=[dyadup(gf,2) 0]; %a trous
gprime = [dyadup(gprime,2) 0];
y=a{k};
end;
O/
/o-------------------------------------------
%Set threshold to zero
%............................................
Threshold = 0.005;
%............................................
%Calculate the modulus maxima for each level
for k = L:-1:1
wtm{k} = sqrt(d1{k}.A2 + d2{k}.A2); %Wavelet Transform Modulus
[mod_m,mod_n] = size(wtm{k});
0/
/o------------------------------------------
%calculate angle of the wavelet transform vector.
%............................................
forqt= 1:nr
for r= 1:nc
alpha{k}(qt,r) = atan2(d2{k}(qt,r),d1{k}(qt,r));
if alpha{k}(qt,r) < 0
alpha{k}(qt,r)=2*pi+alpha{k}(qt,r);
61


end
end
end
%..........................-......
%calculate modulus maxima image
O/
/o--------------------------------
for r=1:mod_m
for c=1:mod_n
mod_max{k}(r,c)=255; Reinitialize modulus maxima array
end
end
Re..........................................
Rofind local maximum of modulus
Re..........................................
ang=pi/8;
for r=2:mod_m-1
for c=2:mod_n-1
if (alpha{k}(r,c)>=(15*ang)/...
alpha{k}(r, c)<=( 1 *ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c) & wtm{k}(r,c) >= wtm{k}(r+1,c)
mod_max{k}(r,c)=0;
elseif (alpha{k}(r,c)>=( 1 *ang)&...
alpha{k}(r,c)<=(3*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c-1) & wtm{k}(r,c) >= wtm{k}(r+1,c+1)
mod_max{k}(r,c)=0;
elseif (alpha{k}(r,c)>=(3*ang)&...
alpha{k}(r, c)<=(5*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r,c-1) & wtm{k}(r,c) >= wtm{k}(r,c+1)
mod_max{k}(r, c)=0;
elseif (alpha{k}(r, c)>=(5*ang)&...
alpha{k}(r, c)<=(7*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c+1) & wtm{k}(r,c) >= wtm{k}(r+1,c-1)
mod_max{k}(r,c)=0;
62


elseif (alpha{k}(r,c)>=(7*ang)&...
alpha{k}(r, c)<=(9*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c) & wtm{k}(r,c) >= wtm{k}(r+1,c)
mod_max{k}(r, c)=0;
elseif (alpha{k}(r,c)>=(9*ang)&...
alpha{k}(r,c)<=( 11 *ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r-1 ,c-1) & wtm{k}(r,c) >= wtm{k}(r+1,c+1)
mod_max{k}(r,c)=0;
elseif (alpha{k}(r,c)>=( 11 *ang)&...
alpha{k}(r,c)<=(13*ang))&...
wtm{k}(r,c)>=Threshold &...
wtm{k}(r,c)>=wtm{k}(r,c-1) & wtm{k}(r,c) >= wtm{k}(r,c+1)
mod_max{k}(r,c)=0;
elseif (alpha{k}(r, c)>=( 13 *ang)&...
alpha{k}(r, c)<=( 15*ang))&...
wtm{k}(r,c)>=Threshold&...
wtm{k}(r,c)>=wtm{k}(r-1 ,c+1) & wtm{k}(r,c) >= wtm{k}(r+1,c-1)
mod_max{k}(r,c)=0;
end
end
end
end
%.............................................
%Build vertical and horizontal detail matrices that include only those...
%...coefficients that are located at the modulus maxima; otherwise, the...
%...coefficients are set to zero.
%.............................................
for k = 1:L
d1_mm{k} = zeros(nr,nc);
d2_mm{k} = zeros(nr,nc);
fori = 1:nr
for j = 1:nc
if mod_max{k}(i,j) == 0
d1_mm{k}(i,j) = d1{k}(i,j);
d2_mm{k}(i,j) = d2{k}(i,j);
end
end
end
63


D1_MM(:,:,k) = d1_mm{k};
D2_MM(:,:,k) = d2_mm{k};
end
mm_nonzeros = nnz(D1_MM) + nnz(D2_MM);
percent = mm_nonzeros/nonzeros
%...............................
%Display modulus maxima matrix mask
%.............................
figure;
%................................
imshow(mod_max{ 1},[]);
xlabel('Modulus Maxima for Level 1');
figure;
O/
/O-------------------------------
imshow(mod_max{2},[]);
xlabel('Modulus Maxima for Level 2');
figure;
%................................
imshow(mod_max{3},[]);
xlabel('Modulus Maxima for Level 3);
64


Appendix E: MATLAB function to compute circular convolution
O/ 0/
/O /0"------------------------------------------
%This function does a circular convolution between an input vector, x,
%and a discrete filter, f.
%............................................
function y = cconv(f,x)
m = length(x);
r = length(f);
if r <= m,
%............................................
%lf filter length is smaller than the input vector:
O/
/O'-----------------------------------------
x_extended = [x((m+1-r):m) x];
else
%lf filter length is bigger than the input vector:
O/
/o------------------------------------------
z = zeros(1 ,r);
for i=1:r,
q r*m -r + i-1;
imod = 1 + rem(q,m);
z(i) = x(imod);
end
x_extended = [z x];
end
%................
%Filter the signal
%................
intermediate_y = filter(f, 1 ,x_extended);
O/
/O---------------------------------------------
%Save only the part of the vector that has the circular convolution results:
%.....-.......................................
y = intermediate_y(r+1 :m+r);
65


Appendix F: MATLAB function to perform left shift of an input vector.
%%......................................
%This function does a left shift of input vector x by a length LengthOfShift
%.....................................
function ShiftedVector=ieftshift(x,LengthOfShift)
LengthOfVector=length(Vector);
for i=1:LengthOfShift
Temp=Vector(1);
Temp Vector= Vector(2:LengthOfVector);
Vector=[TempVector Temp];
end
ShiftedVector= Vector
66


REFERENCES
Donoho, D., Duncan, M. R., Huo, X., & Levi, O. (1999). Wavelab 802 for
MATLAB 5.x [ Library of MATLAB software programs]. Retrieved from
http://www-stat.stanford.edu/~wavelab.
Grochenig, K. (1993). Acceleration of the frame algorithm. IEEE
Transactions on Signal Processing, 41(12), 3331-3340.
Hsung, T., Lun, D.P., & Siu, W. (1998). A deblocking technique for block-
transform compressed image using wavelet transform modulus maxima.
IEEE Transactions on Image Processing, 7(10), 1488-1496.
Hwang, W. L. & Mallat, S. (1992). Singularity detection and processing
with wavelets. IEEE Transactions on Information Theory, 38(2), 617-643.
Mallat, S. (2001). A wavelet tour of signal processing (2nd ed.). New York:
Academic Press.
Mallat, S. & Zhong, S. (1992). Characterization of signals from multiscale
edges. IEEE Transactions of Pattern Analysis and Machine Intelligence,
14(7), 710-732.
Polikar, R. (1999). The engineers ultimate guide to wavelet analysis: The
wavelet tutorial. Retrieved October 13, 2002, from the Rowan University
College of Engineering Web site:
http://engineering.rowan.edu/~polikar/WAVELETS/WTtutorial.html.
Shensa, M. (1992). The discrete wavelet transform: Wedding the a trous
and Mallat algorithms. IEEE Transactions on Signal Processing, 40(10),
2464-2482.
67


Full Text

PAGE 1

WAVELET TRANSFORM MODULUS MAXIMA ALGORITHMS FOR MEDICAL IMAGE PROCESSING by Saritha Jagannath A thesis submitted to the University of Colorado at Denver in partial fulfillment of the requirements for the degree of Master of Science Electrical Engineering 2004

PAGE 2

This thesis for the Master of Science degree by Saritha Jagannath has been approved by Miloje S. Radenkovic Hamid Fardi I :L /0 Date

PAGE 3

Jagannath, Saritha (M.S., Electrical Engineering) Wavelet Transform Modulus Maxima Algorithms for Medical Image processing Thesis directed by Professor Jan Bialasiewicz ABSTRACT Wavelets are powerful mathematical tools which represent the signal as a combination of different frequency components, and then study each component with a resolution matched to its scale. They have an edge over traditional Fourier methods in analyzing discontinuities and sharp spikes in a signal. The analysis of signals and images at different scales is an attractive framework for many problems in signal and image processing. Wavelet theory provides a mathematically precise understanding of the concept of multiresolution. The focus of this thesis is on application of wavelets in the field of medical image processing. Wavelets are used as a tool to analyze CT scan images. Wavelet analysis is carried out on healthy human tissue and malignant tissue in an attempt to clearly bring out the differences in the behaviour between the two. The filter bank algorithm called algorithme a trous is used for the above analysis. The analysis is carried out by detecting the edges in the images and looking for a pattern difference between the malignant tissue and the healthy tissue going across frequency bands. Horizontal and vertical details are calculated in different frequency bands lll

PAGE 4

and the corresponding modulus maxima images are found. Also Wavelet packet analysis of healthy image and malignant image is carried out in different frequency bands and the corresponding edges are found. This abstract accurately represents the content of the candidate's thesis. I recommend its publication. Signed Jan T. Bialasiewicz lV

PAGE 5

DEDICATION I dedicate this thesis to my parents. If not for them I would not have been where I am today. My deepest gratitude goes to them for believing that I could actually come this far and get this degree. I truly thank them for all their support and blessings.

PAGE 6

ACKNOWLEDGEMENT My sincere thanks to my advisor, Jan Bialasiewicz, for introducing me to the world of wavelets and providing me with an opportunity to explore the world of wavelets. Also, I wish to thank the University of Colorado Health Science Center for providing me with kind of material and support that was required for this thesis.

PAGE 7

CONTENTS Figures ...................................................................................................... viii Tables ........................................................................................................... x Chapter 1. Introduction to Wavelet Transforms .......................................................... 1 2. Discrete Dyadic Wavelet Transform ......................................................... 8 2.1 Dyadic Wavelet Transform Computed with "Aigorithme a trous" .......... 11 3. Analysis of Images in Different Frequency Bands ................................. 20 3.1 Medical Image Analysis in Different Frequency Bands ......................... 23 4. Introduction to Wavelet Packet Analysis ................................................. 39 4.1 Wavelet Packet Analysis of Medical Images ........................................ 42 5. Conclusion .............................................................................................. 48 Appendix A. MATLAB function to read in the test image ............................................ 50 B. MATLAB function to generate low pass and high pass decomposition filters ............................................................................. 51 C. Algorithme a trous MATLAB computer program .................................... 52 D. MATLAB function to generate modulus maxima masks for detail images .................................................................................... 59 E. MATLAB function to compute circular convolution ................................. 65 F. MATLAB function to perform left shift of an input vector ....................... 66 References ................................................................................................. 67 VII

PAGE 8

FIGURES 1 .1 Time-frequency boxes ("Heisenberg rectangles") representing the energy spread of Gabor atoms .......................................................... 2 1.2 Daubechies, db4, and Symmlet, sym8, Mother Wavelets ...................... 4 1.3 Scaling Functions of the Daubechies, db4, and Symmlet, sym8 Mother Wavelets ............................................................................ 4 1.4 Time-frequency boxes of two wavelets 'If"' and 'lfu,.,, ............................. 6 2.1 Filter bank for two-dimensional image decomposition implementing the Mallat multiresoluton algorithm ................................. 10 2.2 The Spline Wavelet (left) and its Scaling Function (right) ..................... 13 2.3 Filter bank for two-dimensional image decomposition implementing the "algorithms a trous" ................................................... 14 2.4 Filter bank implementing image reconstruction with the "algorithms a trous" ............................................................................... 16 2.5 Use of phase to determine the Modulus Maxima .................................. 18 3.1 Cross section of neck image with tumor area marked .......................... 25 3.2 Modulus maxima image for scale 2 ...................................................... 26 3.3 Modulus maxima image for scale 4 ...................................................... 27 3.4 Modulus maxima image for scale 8 ...................................................... 28 3.5 (a) Cross-section of malignant tissue with tumor area marked, (b) Modulus maxima image for scale 2, (c) Modulus maxima image for scale 4, (d) Modulus maxima image for scale 8 .................... 30 3.6 (a) Cross-section of healthy neck, (b) Modulus maxima image at scale 2, (c) Modulus maxima image at scale 4, (d) Modulus maxima image at scale 8 ................................................... 31 3. 7 Modulus maxima image at scale 2-5 .................................................... 33 3.8 Modulus maxima image at scale 2-4 ..................................................... 34 3.9 Modulus maxima image at scale 2-3 ..................................................... 35 Vlll

PAGE 9

3.1 0 (a) Cross-section of malignant tissue with tumor area marked, (b) Modulus maxima image for scale 25,(c) Modulus maxima image for scale 2-4 (d) Modulus maxima image for scale 2-3 ................ 36 3.11 (a) Cross-section of healthy tissue, (b) Modulus maxima image for scale 2-5 (c) Modulus maxima image for scale 2-4 (d) Modulus maxima image for scale Z3 ............................................... 37 4.1 Binary tree of wavelet packet spaces ................................................... 41 4.2 (a) Tree structure depicting three levels of details (b) Modulus maxima image at node W,' at scale 2, (c) Modulus maxima image at node W2 3 at scale 4, (d) Modulus maxima image at node W3 7 at scale 8 ............................................................................................... 44 4.3 (a)Tree structure depicting three levels of details of healthy tissue (b) Modulus maxima image at node W,' at scale 2, (c) Modulus maxima image at node W2 3 at scale 4, (d)Modulus maxima image at node W3 7 at scale 8 .................................. 45 4.4 (a) Tree structure depicting three levels of details of healthy tissue (b) Modulus maxima image at node W1 1 at scale 2-5 (c) Modulus maxima image at node W2 3 at scale 2-4 (d) Modulus maxima image at node W3 7 at scale 2-3 .............................. 46 4.5 (a) Tree structure depicting three levels of details of healthy tissue (b) Modulus maxima image at node W1 1 at scale 2-5 (c) Modulus maxima image at node W/ at scale 2-4 (d) Modulus maxima image at node W3 7 at scale 2-3 .............................. 47 IX

PAGE 10

TABLES 3.1 Spline filter coefficients at scale 1 ........................................................ 22 X

PAGE 11

1. Introduction to Wavelet Transforms Wavelet transforms are an important class of mathematical tools used in medical signal processing as they help in unveiling hidden information. Wavelet transform is capable of providing the signal content as a function of time and frequency. It can be considered as an extension of the Fourier transform. The major drawback of Fourier transform is that it has only frequency resolution and no time resolution which means that although we might be able to determine all the frequencies present in a signal, we do not know when or where they are present. The Windowed Fourier transform overcomes this drawback to a certain extent in that it replaces the sinusoidal wave of Fourier transforms by a sinusoid and a window which is localized in time however has frequency resolution problems. Low frequencies can hardly be depicted by short windows and high frequencies can only be poorly localized by long windows. The windowed Fourier transform, defined by Gabor, correlates a signal f with each atom g r as represented by the following equation: II,';. Sf(u, q) = J J(t )g :.q (t ):tr = J J(t )g(t-u dt ( 1 .1) where u represents the time instant and ; represents frequency.

PAGE 12

It is a Fourier integral that is localized in the neighborhood of u by the window g(t-u). The windowed Fourier transform can also be written as a frequency integral by applying the Fourier Parseval formula Sf(u,;) = -1 --JJ(m)g,:)! (m}im 27r -00 ., (1.2) The transform SJ(u,;) thus depends only on the values J(t) and ](w)in the time and frequency neighborhoods where the energies of g and u ., are concentrated (Mallat, 2001 ). Figure 1.1 shows the timefrequency boxes of g and g r,y representing the energy spread of Gabor atoms. w lg ( ro) I -==o:=t ==; r-1 .. I I (j(JJ lgu (t) I .., lg Ct) I v,y Figure 1.1 Time-frequency boxes ("Heisenberg rectangles") representing the energy spread of Gabor atoms. 2

PAGE 13

The limitations in resolution were the main reasons for invention of wavelet transforms. The wavelet transform decomposes signals over dilated and translated windows. The original window centered around t=O and for which the scale S=1, is called a mother wavelet, which is a function'I'E L2(R)with zero average, i.e., 00 Jljl(t)dt = 0 (1.3) We obtain a family of time-frequency atoms by scaling 'I' by s and translating it by u (Mallat, 2001), i.e., ) I "j t-u) f/lu.s (t = Fs r l-s(1.4) The scale s indirectly corresponds to frequency. Larger scale means dilated wavelet and lower frequencies. Smaller scale means compressed wavelet and thus higher frequencies. The factor t/ ..[; is a normalization factor to ensure that all wavelets have the same energy and this energy is unity if ll'l'(t (Mallat, 2001 ). Figure 1.2 shows an example of two different mother wavelets: Daubechies (db4) and Symmlet (sym8). 3

PAGE 14

II 05 /' (\ -05 ., Daubeth1esW:noelel IOL ___ ___ __Ji5 SymmleiW
PAGE 15

Like windowed Fourier transform, wavelet transform can measure the time frequency variations of spectral components but it has a different time-frequency resolution. Wavelets transform correlates f with If/, __ ,. By applying the Fourier Parseval formula, it can also be returned as frequency integration: +oo I WJ(u,s) = lf(t )Vt:.q(t = 2Jr l](m)ift:.s (1.6) In time, 'II,_, is centered at u with a spread proportional to s. Figure 1.4 shows the time-frequency boxes of If/,,, and lf/,0_,11 When the scale s decreases, the time support is reduced but the frequency spread increases and covers an interval that is shifted towards high frequencies (Mallat, 2001 ). 5

PAGE 16

!l s !!_ s .. 1\ \ ii' I 1\ II + il s rr -' 1'1\.. s< W)l n u ) (\ vu:, scr II I + q_w_ I Su \jl u,,,S II I \ \ \ 'uu Figure 1.4 Time-frequency boxes of two wavelets lfu.s and If""-'" We obtain the continuous wavelet transform of a signal by scaling the mother wavelet, shifting the wavelet in time, multiplying with the signal and integrating over all times. The obtained coefficients correspond to the details of the signal since wavelet is considered as a high pass filter. To recover fwe also need a complement of information corresponding to Wf (u, s) that is obtained using the scaling function
PAGE 17

The low-frequency approximation of fat scale s is LJ(u,s )= \J(t ), j; u )) (1.8) If the low frequency approximation was obtained at S=s0 then the original signal can be recovered using the following equation: I so ds I J(t) =-JwJ(.,s )* )-2 +-Lf(.,s0)* ,10 (t) C 0 s Cs0 (1.9) 2 where C = J dw 0 (() (1.10) The first term in equation ( 1.9) represents the part of recovered from its high frequency representation for scales s ::; s0 and the second term represents the low frequency part of The condition: C = J 111 w" dw < oo 0 (() (1.11) is called the wavelet admissibility condition for a real function, L2(R), to be a wavelet (Mallat, 2001). 7

PAGE 18

2. Discrete Dyadic Wavelet Transform In the computer world, an image is a collection of pixels at a certain resolution. Thus continuous wavelet transform cannot be computed at arbitrary small scales. To construct a discrete wavelet from a continuous wavelet, the scale s is discretized but not the translation parameter u. For fast computations the scale is sampled along a dyadic sequence {21 Lez. The scale parameter 21 is the inverse of resolution T1 The dyadic wavelet transform of a signal f E L2 (R) is defined by I *(t-uG Wf(u,21 ) = fi! 2J Jt (2.1) In computer algorithms, an image is analyzed using discrete filters. Discrete high pass and low pass filters, called the conjugate mirror filters, replace continuous wavelets and scaling functions, respectively. The filters are related to their continuous counterparts as follows: (2.2) (2.3) 8

PAGE 19

The high pass filter g[n] is related to the low pass filter h[n] as follows: g[n] =(-It" h[l-n] (2.4) 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 impulse response of a filter. The dyadic wavelet transform (DWT) analyzes the image at different frequency bands with different resolutions by decomposing it into a coarse approximation and detail information with high-pass and low-pass filters abbreviated as HP and LP, respectively. 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-dimensional convolutions along the rows and the columns of the image. At any level, the decomposition of the previous level image approximation a 1 is performed with a0 representing the original image. 9

PAGE 20

a. J LP HP LP HP LP dv j+l Figure 2.1 Filter bank for two-dimensional image decomposition implementing the Mallat multiresoluton algorithm. The multiresolution analysis involves a decomposition of the image space into a sequence of subspaces Vi. More formally, the approximation of a function at a resolution r i is defined as an orthogonal projection on a space vi c L2 (R) (Mallat 2001 ). Orthogonal projection of a function f is defined by a function f1 that minimizes the distance III !1 II 10

PAGE 21

As illustrated in Figure 2.1, the following are the steps of the algorithm (known as Mallat's multiresolution algorithm): First, the rows of a j are passed through an HP and LP filter. This generates two images, which are then down-sampled by two (indicated by the down-arrow), i.e., 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. This generates four images: approximation image a J+l, the horizontal detail d J + 1 the vertical detail d)'+ 1 and the diagonal detail d 1 + 1 Typically, quadrature mirror filters such as Daubechies filters are used for perfect reconstruction. These have a special relationship between the LP and HP decomposition filters. 2.1 Dyadic Wavelet Transform Computed with "Aigorithme a trous" Dyadic wavelet transform can also be computed with a filter bank algorithm called in French the "algorithme a trous". The "algorithme a trous" uses quadratic spline filters properly related to the quadratic spline 11

PAGE 22

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. A box spline of degree m is a translation of m+ 1 convolutions of 1 ro. 1J with itself. It is centered at t = Y2 if m is even and at t = 0 if m is odd. Its Fourier transform is 0 OJ smOJ 2 2 m+l withE= 1 if m is even and E = 0 if m is odd. (2.5) From equation 2.2, the relation between the low pass filter and the scaling function is frequency domain is as follows: (2.6) where h(OJ) is the Fourier transform of h[n]. We can construct a wavelet that has one vanishing moment by choosing g(m) = o(m) in the neighborhood of m= 0. For example ,OJ r;:; _,_ (I) g(w) = -iv Le 2 sin-2 12 (2.7)

PAGE 23

The Fourier transform of the resulting wavelet is (2.8) (Mallat, 2001) The mother wavelet and the scaling function of quadratic spline are shown in Figure 2.2. Figure 2.2. The Spline Wavelet (left) and its Scaling Function (right) In the implementation of this algorithm for image decomposition, all filtering is done with circular convolutions to prevent the discontinuity induced by finite-length signals near image boundaries. After each level j of decomposition, the filters are dilated by inserting 2j-1 zeros between 13

PAGE 24

consecutive filter coefficients. This creates holes ("trous" in French). There is no down sampling only dilation of the filters. a. 1 columns rows LP LP HP HP .. aJ+l dl' )+1 1 j+ Figure 2.3 Filter bank for two-dimensional image decomposition implementing the "algorithme a trous". One of the decomposition levels of this algorithm is illustrated in Figure 2.3. The approximations are generated by first shifting the columns of a 1 to the left 2j-l times. Then, the columns are filtered using the LP filter h[n]. Next, the rows are shifted to the left 2i1 times and filtered with an LP filter. The result is the approximation a J+l of the image. The horizontal details are calculated by first filtering the columns of a 1 using the HP filter g[n] and then shifting the columns 2i times. The vertical 14

PAGE 25

details are calculated by first filtering the rows of a J using the HP filter and then shifting the rows 2i times. The following equations used in the MATLAB program give the convolution formulas that are cascaded to calculate the discrete dyadic wavelet transform and its inverse: DWT: (2.9) (2.1 0) (2.11) IDWT: (2.12) Note that u1 = dilated analysis filter u 1 and y1 = decimated synthesis filter y1 o[n] is discrete Dirac and xy[n,m]=x[n] y[m] is a separable filter (Mallat, 2001 ). Reconstruction of the image is followed in reverse order from the decomposition procedure with the exception of shifting that is altered. After each level of reconstruction, the filters are decimated by inserting 15

PAGE 26

2j -1 zeros between the consecutive coefficients. Note that j decreases at each level of reconstruction. The synthesis filters HP g[n] and LP h[n]are used and the resultant approximation and detail images are added together and multiplied by 1/4. The filter bank implementing this reconstruction process is shown in Figure 2.4. row columns h[n] h[n] .. LP LP +. xl 4 g[n] HP g[n] HP Figure 2.4. Filter bank implementing image reconstruction with the "algorithme a trous". At every decomposition level and at any pixel of the image, we have available horizontal and vertical detail calculated using the "algorithme a 16 a. J

PAGE 27

trous". These are orthogonal components of the wavelet transform. Using these components, we calculate the wavelet transform modulus (WTM) and angle using the following equations: M jf(x, y) = (x, y))2 +d). (x, y))2 (2 1 3) A j f ( x, y) = tan -1 ( d; ( x, y) / d J ( x, y)) (2.14) For each pixel location within each scale, a modulus maxima is detected at location (x,y) if MJ(x,y)>MJ(x1,y1 ) andM1J(x,y)>MJ(x2,h) Locations (x1 yJ and (x2 y2 ) are adjacent locations to (x, y) indicated by the angle AJ(x, y). The angle can be one of eight quantized directions: 0, n/4, n/2, 3n/4, n, n/4, -n/4, -n/2, or -3n/4 (Hsung, Lun, & Siu, 1998). 17

PAGE 28

'][ ------;: \----,_ I 0 "'/ / l \ _..--/ /. // \ 9n/8\\ / / 15 n/8 \ / \( / 5n/4 or -3n/4 1 // 7n/4 or -n/4 --""' ,' 11 4 13 n/ 8 or -'J[/2 Figure 2.5 Use of phase to determine the Modulus Maxima. Figure 2.5 above is a graphical description of how the angle is used to locate the two points in which to compare a particular wavelet modulus. For example, if the angle that corresponds to point (x, y) lies between the n/8 and 3n/8, indicated by dashed lines, the quantized angle becomes n/4. The line tangent to this angle, shown by a short line, correspond to the direction in which the two adjacent positions to be compared are located. In this example, the wavelet transform modulus located at (x, y), indicated by the black square, is compared to the wavelet transform modulus located at (x" yJ and (x2 y2), indicated by 18

PAGE 29

the two shaded squares. Positive angle values that correspond to the negative angle values are included in the figure because the MATLAB programs in Appendix C and Appendix D work with positive angle values only. 19

PAGE 30

3. Analysis of Images in Different Frequency Bands Wavelet decomposition provides natural setting for the multi-level edge detection. Since Wavelet Transform Modulus Maxima (WTMM) provide useful information for pattern recognition, we propose to use it here for fast feature extraction. It is possible to analyze an image in different frequency bands. This is accomplished by assuming that the original image is initially at some scale 21 where j E z With this assumption, the image is decomposed by going higher in scale from the starting scale. The MATLAB function 'scaling ()', attached in Appendix B, scales the mother wavelet and scaling function of quadratic splines. From the scaled functions the corresponding HP and LP filters are obtained. Using these filters the image is analyzed in the corresponding frequency band. The scaling of the wavelet and scaling function is done in the frequency domain for easy computation purposes. The analysis is carried out using a threshold, which produces the best visual results. 20

PAGE 31

To design dual scaling functions and wavelets lji, which are splines, we choose h =h. As a consequence, = and the reconstruction condition implies that [ 12 J 2h(m) .w m 211 g(OJ)= =-i..fie-'2 sin OJ:L(cosOJ) g (OJ) 2 n=O 2 (3.1) The scaling and wavelet function for the quadratic spline, defined in frequency domain, are (3.2) (3.3) The HP and LP decomposition filters are calculated as follows: (3.4) (3.5) where j E Z (Mallat, 2001) 21

PAGE 32

Table 3.1 gives the coefficients of quadratic spline filters at scale 1. n h[n] h[n] g[n} g[n] -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 Table 3.1 Spline filter coefficients at scale 1 In this table, h[n]low pass decomposition filter. h[n]-low pass reconstruction filter. g[n]-high pass decomposition filter. g[n]high pass reconstruction filter. The low pass filters h[n] = h[n} are obtained by Inverse Fourier Transform of equation (3.4) and the high pass filters g[n} and g[n] are obtained by Inverse Fourier Transform of equations (3.5) and (3.1 ). 22

PAGE 33

3.1 Medical Image Analysis in Different Frequency Bands In this section, the above theory of analyzing an image in different frequency bands is applied to Computed Tomography (CT) scan images. One of the advantages of CT over standard x-rays is that it clearly shows soft tissue such as the brain, as well as dense tissue like bone. Wavelet analysis is carried out at three levels going higher in scale and modulus maxima image is found at each level. This analysis is carried out in different frequency bands by assuming that the original image is at some scale 21 First the image is read using the function 'ca/_mm(j, /vi, image)' attached in Appendix A where j corresponds to scale 2i, /vi is number of decomposition levels and image is the image under analysis. The function 'ca/_mm(j, /vi, image)' calls function 'scaling(})' attached in Appendix B which calculates LP and HP decomposition filters corresponding to scale 2i. Equations (3.2) and (3.3) are used in the function to calculate the same. The number of decomposition levels and image under analysis along with the above obtained filters are passed to the function 'mm_atrous(lvl,y, h,g,gp)' attached in Appendix C. This function calculates the WTMM, Wavelet Transform Modulus Maxima 23

PAGE 34

Masks, for every level. Figure 3.1 gives the CT scan image of cross section of human neck with tumor area marked. Modulus maxima are calculated for this image at different scales. Modulus maxima images at scales 2, 4 and 8 are shown in Figure 3.2, Figure 3.3, and Figure 3.4, respectively. For this case, the image shown in Fig.3.1 is assumed to be at scale 1. 24

PAGE 35

Figure 3.1. Cross section of neck image with tumor area marked. 25

PAGE 36

Figure 3.2 Modulus m . axlma Image for scale 2 26

PAGE 37

rna image for scale 4 3 3 Modulus maxi Figure 27

PAGE 38

. e for scale 8 I maxima lmag 4 Modu us Figure 3. 28

PAGE 39

The tumor area is marked in each case and it is observed that there is an increase in edges in the tumor area when going from scale 2 to 4. This is in contrast to tumor free areas where there is a decrease in the number of edges as we go higher in scales. This is more clearly seen in Fig.3.5 where the tumor area of the cross section is cut out and analyzed with an assumption that the original image is at a scale 1 This pattern is not seen in the healthy tissue. Here the numbers of edges reduce when the scale is increased from 2 to 4. Figure 3.6 shows the analysis results of the healthy tissue. In this case, we can see that going higher in scale the number of edges uniformly decreases. 29

PAGE 40

Figure 3.5 (a) Cross-section of malignant tissue with tumor area marked, (b) Modulus maxima image for scale 2, (c) Modulus maxima image for scale 4, (d) Modulus maxima image for scale 8. 30

PAGE 41

(a) Figure 3.6 (a) Cross-section of healthy neck, (b) Modulus maxima image at scale 2, (c) Modulus maxima image at scale 4, (d) Modulus maxima image at scale B. 31

PAGE 42

When analyzed in different frequency band, that is, by assuming that the original image is at a scale 2 j where j : 0, the edges in the tumor area reduce with scale unlike the case when j=O. Figures 3.7, 3.8 and 3.9 give the modulus maxima images corresponding to Figure 3.1 at scales 25 24 and 23 respectively. Smaller scales correspond to higher frequencies. As we go higher in frequency, the number of edges detected increase. There is denser pattern observed. However, there is no difference in pattern between the healthy tissue and malignant tissue. Figure 3.1 0 and Figure 3.11 give the modulus maxima images of malignant and healthy tissue cross-sections at scales 25 24 and 23 32

PAGE 43

Figure 3.7 Modulus maxima image at scale 2-5 33

PAGE 44

Figure 3.8 Modulus maxima image at scale 2-4 34

PAGE 45

Figure 3.9 Modulus maxima i mage at scale 2-3 35

PAGE 46

Figure 3.10 (a) Cross-section of malignant tissue with tumor area marked, (b) Modulus maxima image for scale 2 5 (c) Modulus maxima image for scale 24 (d) Modulus maxima image for scale 2 3 36

PAGE 47

. (a) c Figure 3.11 (a) Cross-section of healthy tissue, (b) Modulus maxima image for scale 2-5 (c) Modulus maxima image for scale 2-4 (d) Modulus maxima image for scale 2-3 37

PAGE 48

From the above analysis it follows that by wavelet analysis of an image in different frequency bands, it could be possible to see some kind of difference in behavior between the healthy tissue and malignant tissue. This can be used as a tool for detection of cancer. 38

PAGE 49

4. Introduction to Wavelet Packet Analysis A wavelet transform is a method of expressing a signal in terms of the wavelet basis {'I' /f)}. This accounts to decomposing the signal space L2 into a direct sum of orthogonal subspaces {vi} and { wi }. If the signal bandwidth is finite with information upto a resolution level J, a wavelet transform performs a decomposition of the space V1 into a sum of orthogonalsubspaces (4.1) Of course this is not the only way to decompose the signal space L2 or VJ Wavelet packets, introduced by Coifman, Meyer and Wickerhauser, generalize the link between multiresolution approximations and wavelets. (Mallat, 2001 ). From the theory of multiresolution analysis, we know that given the basis functions {i(t-2in)} of Vi, {i+1(t-2i+1n)}and {'l'i+l(t-2i+ln)J form an orthonormal basis of Vi+l and wi+l' respectively, and Vi= Vi+l Ef> Wi+l. The relations between these basis functions are represented by the following equations: 39

PAGE 50

(4.2) (4.3) n;;:---oo Where h and g are conjugate mirror filters and g[n]= (-It"h[I-n] (4.4) So the space V1 can be decomposed into a sum of two orthogonal subspaces defined by their basis functions given by the equations (4.2) and (4.3). Using a pair of conjugate mirror filters to decompose a signal space corresponds to splitting the frequency content of the signal into a low-frequency and a high-frequency component. In orthogonal wavelet decomposition procedure, a signal space is decomposed into approximation space and details space. Next, the approximation space is further decomposed into next level of approximation and details and so on. The details space is never decomposed into coarser scales. In wavelet packet analysis, each details space is also decomposed into two spaces using the same approach as in decomposition of approximation space. So in general, wavelet packet decomposition divides the frequency space into various parts and allows better frequency localization of signals. 40

PAGE 51

Instead of decomposing only the approximation spaces v, into two orthogonal subspaces, we can decompose the detail spaces w, as well. The recursive splitting of the vector spaces can be represented as a binary tree shown in Fig.4.1. / wo I / / wo 2 WI 2 w22 w23 Figure 4.1 Binary tree of wavelet packet spaces Each node U,p) corresponds to a space Wf, which admits an orthonormal basis {'If! (t2' n)} At the root of the tree, we have W0 = V0 Analogously, we can define two orthogonal bases at the children nodes as follows: (4.5) n=--o::o (!) = L h[n )ltJ' (1-2 J n) (4.6) fi=---OCI Thus far we have been using the combination of { /1-2' n)} and {'If, (12' n)} to form a basis forv,, and now we have a whole sequence of 41

PAGE 52

functions at our disposal. Various combinations of these and their dilations and translations can give rise to various bases for the fun.ction space. So we have a whole collection of orthonormal bases generated. 4.1 Wavelet Packet Analysis of Medical Images The transformation of data into wavelet packet basis presents no extra numerical difficulties. We can simply do a convolution using filters h and g on the details space as well as on the approximation space. As in the wavelet transform, we can keep doing the decomposition until we cannot go any further. On the other hand, we can also choose not to decompose a particular subspace while decomposing others. So there are many choices for decomposing a signal. We can keep all the coefficients at all decomposition levels and generate a table of coefficients of wavelet packet decomposition. The MATLAB function 'wave/et_packet_analysis_of_detail_images(lvl,y, h,g,gp)', given in Appendix 0, is used to achieve this. The function 'cal_mm(j, /vi, image)' calls function 'scaling(})' attached in Appendix A which calculates LP and HP decomposition filters corresponding to scale 2i. Equations (3.2) and (3.3) are used in the function to calculate the same. The number of decomposition levels and image under analysis along with the above obtained filters are passed to the function 42

PAGE 53

'wave/et_packet_analysis_ot_detail_images(lvl,y,h,g,gp)' attached in Appendix D. This function calculates the WTMM masks for every level of detail images. The modulus maxima images are obtained by initializing the detail images at every level as the image to be analyzed for the next level. Figure 4.1 gives the modulus maxima images of malignant tissue for three levels of detail images. Figure 4.2 gives the modulus maxima images of healthy tissue for three levels of detail images. The image is assumed to be initially at scale 1. Figure 4.3 and Figure 4.4 give the analysis results for smaller scales corresponding to high frequencies of malignant and healthy tissue respectively. No difference in pattern between the healthy tissue and malignant tissue was observed in WTMM masks obtained using wavelet packet transform of detail images. 43

PAGE 54

(a) Figure 4.2 (a) Tree structure depicting three levels of details (b) Modulus maxima image at node W,' at scale 2, (c) Modulus maxima image at node W/ at scale 4, (d) Modulus maxima image at node W3 7 at scale 8. 44

PAGE 55

(a) Figure 4.3 (a) Tree structure depicting three levels of details of healthy tissue (b) Modulus maxima image at node W1 1 at scale 2, (c) Modulus maxima image at node W2 3 at scale 4, (d) Modulus maxima image at node W3 7 at scale 8. 45

PAGE 56

/ wo I w.l I / w2 2 w3 2 w6 3 w? 3 (a) Figure 4.4 (a) Tree structure depicting three levels of details of healthy tissue (b) Modulus maxima image at node W1 1 at scale 2-5 (c) Modulus maxima image at node W2 3 at scale 2-4 (d) Modulus maxima image at node W/ at scale 23 46

PAGE 57

I I (a) Figure 4.5 (a) Tree structure depicting three levels of details of healthy tissue (b) Modulus maxima image at node at scale 2 5 (c) Modulus maxima image at node at scale 2-4 (d) Modulus maxima image at at scale 23 47

PAGE 58

5. Conclusion An extensive analysis of medical images was carried out using the Algorithme a trous and quadratic spline wavelets. The images were analyzed by detecting edges and looking for a pattern difference between the malignant tissue and the healthy tissue going across frequency bands. Horizontal and vertical details were calculated in different frequency bands and the corresponding modulus maxima images were found. Difference in pattern between the healthy tissue and malignant tissue was found while going from scale 1 to 2. In case of healthy tissue, the number of edges decreased from scale 1 to 2 whereas in the malignant tissue there was an increase in the number of edges. The images were analyzed at both positive and negative scales by assuming that the original image is at some scale {2 1 tEz Also wavelet packet analysis of healthy image and malignant image was carried out in different frequency bands and the corresponding modulus maxima images were found. The algorithm developed in this thesis was able to detect the differences in behavior between healthy and malignant tissues only in the modulus maxima images obtained by wavelet analysis. These differences could not 48

PAGE 59

be seen in modulus maxima images obtained by wavelet packet analysis. Also the differences were seen only when going from scale 1 to 2. In addition, it should be noted that representing a single CT scan image in different frequency bands and showing the edges of the obtained images, we are providing the information not available by looking at the original image. However, the radiologists have to learn how to use this information for medical diagnosis. 49

PAGE 60

APPENDIX Appendix A: MATLAB function to read the test image. function cal_mm(scale,lvl,image) A=imread(image); %Convert RGB image to grayscale A=double(A); Xrgb=0.2990*A(:,:, 1)+0.5870*A(:,:,2)+0. 1140*A(:,:,3); NbColors=256; X=wcodemat(Xrgb,NbColors); map=gray(NbColors); %Normalize image X= y=X; save_orig = y; %Save original image to display later [nr,nc]=size(y); %Convert 2-D image into a 1-D vector orig_ vector = y( 1, :) ; fori= 2:nr orig_ vector= [orig_ vector y(i,:)]; end %Decompose image %a = last level approximation %01_MM = contains all levels of the horizontal details corresponding to the modulus maxima matrix masks %02_MM = contains all levels of the vertical details corresponding to the modulus maxima matrix masks %gprime = reconstruction high pass filter %hprime = reconstruction low pass filter [h, g, gp ]=scaling( scale); [a, D1_MM, D2_MM, gprime, hprime] = mm_atrous(lvl,y,h,g,gp); 50

PAGE 61

Appendix 8: MATLAB function to generate low pass and high pass decomposition filters %This function is called by cal_mm(scale,lvl,image) to help generate the low pass and high pass decomposition filters. %This function generates the high pass and low pass decomposition filters for particular frequency band. %Inputs %j where i gives the scale parameter %Outputs %h=low pass decomposition filter %g=high pass decomposition filter %gp=high pass reconstruction filter function [g,h,gp]=scaling(j) d=1; for n=-2: 1:3 syms wa b; phi=(((sin(w/2))/(w/2))"'3) exp(-i*w/2); psi=(-i*w/4)*(((sin(w/4))/(w/4))A4)*exp(-i*w/2); H=((sin(2"'j*w)/(2"'j*w))"'3)*exp(i*2"'j*w)/(((sin(2"'j*w/2)1(2"'j*w/2))"'3)*exp(-i*2"j *w/2)); H=sqrt(2) *H*exp(i*w*n); G=(-i*2"'j*w/2)*(((sin(2"'j*w/2))/(2"'j*w/2))A4)*exp( i*2"'j*w)/(((sin(2"'j*w/2)/(2"'j*w/2))"'3)*exp(-i*2"'j *w/2)); G=sqrt(2) *G *exp(i*w*n); GP=-i*sqrt(2) *exp(i*2"'j*w/2) *sin(2"'j*w/2) ( 1 +( cos(2"'j*w/2) )"'2+( cos(2"'j*w/2) )1'14) *exp(i*w*n); a=-pi; b=pi; g(d)=(1/(2*pi))*INT(H,w,a,b); h(d)=(1/(2*pi))*INT(G,w,a,b); gp(d)=(1/(2*pi))*INT(GP, w,a,b); d=d+1; end 51

PAGE 62

Appendix C: Algorithm a trous MATLAB computer program %mm_atrous(lv/,y, h,g,gp) %This function is called by ca!_mm(scale,lvl,image) to help generate the %modulus maxima matrix masks. %This function generates the masked horizontal and vertical details along %with the last approximation. %-------------------------------------------------------%Inputs: %-------------------------------------------------------%/vi = number of levels of decomposition %y =original image %h=low pass decomposition filter %g=high pass decomposition filter %gp=high pass reconstruction filter %----------------------------------------------------------%Outputs: %----------------------------------------------------------%y =last level approximation image %01_MM =contains masked horizontal details from every level %02_MM = contains masked vertical details from every level %gprime = reconstruction high pass filter %hprime = reconsrtuction low pass filter %Finally the program displays the modulus maxima matrix masks for each %/eve/. function [y, D1_MM, D2_MM, gprime, hprime] = mm_atrous(lvl,y,h,g,gp) %----------------------------------------------------------%Spline wavelet filters %----------------------------------------------------------hf=double(h) gf=double(g) gprime=double( gp) L=lvl; %L=Ievel of detail 52

PAGE 63

[nr,nc]=size(y); nonzeros = L *2*nr*nc; 91a---------------------------------------------------------91aLoop for decomposing image 91a----------------------------------------------------------fork= 1 :L 91a---------------------------------------------------------91acalculates approximations. 91ashift either image or latest approximation vertically. 91a----------------------------------------------------------for i=1:nc a{k}(1 :nr,i) = leftshift(y(1 :nr,i)',2/\(k-1 ))'; 91a---------------------------------------------------------91acircular convolution of hf with columns from vertically shifted image 91a----------------------------------------------------------a{k}(1:nr,i) = cconv(hf,a{k}(1:nr,in'; end 91a---------------------------------------------------------91ashift either image or latest approximation horizontally 91a----------------------------------------------------------for i=1:nr a{k}(i, 1 :nc) = leftshift(a{k}(i, 1 :nc),2/\(k-1 )); 91a---------------------------------------------------------91acircular convolution of hf with rows from horizontally shifted image 91a----------------------------------------------------------a{k}(i, 1 :nc) = cconv(hf,a{k}(i, 1 :nc)); end; 53

PAGE 64

horizontal and vertical details. convolution of high pass filter, gf, with columns from image or ... approximation for i=1:nc d1 {k}(1 :nr,i)=cconv(gf,y(1 :nr,i)')'; resultant image vertically. The result is a horizontal detail ... 0/ /o . tmage d1 {k}(1 :nr,i)=leftshift(d1 {k}(1 :nr,i)',2"k)'; end; convolution of high pass filter, gf, with rows from image or ... .. previous approximation for i=1 :nr d2{k}(i, 1 :nc)=cconv(gf,y(i, 1 :nc)); resultant image horizontally. The result is a vertical detail ... d2{k}(i, 1 :nc)=leftshift(d2{k}(i, 1 :nc),2"k); end; d{k}=d1 {k}+d2{k}; filters 54

PAGE 65

hf=[dyadup(hf,2) OJ; %a trous hprime=hf; gf=[dyadup(gf,2) OJ; %a trous gprime = [dyadup(gprime,2) OJ; y=a{k}; end; %---------------------------------------------------------%Set threshold to zero %----------------------------------------------------------Threshold = 0. 005; %----------------------------------------------------------%Calculate the modulus maxima for each level %----------------------------------------------------------fork= L:-1:1 wtm{k} = sqrt(d1{k}/'2 + d2{k}/'2); %Wavelet Transform Modulus [mod_m,mod_nJ = size(wtm{k}); %----------------------------------------------------------%calculate angle of the wavelet transform vector. %----------------------------------------------------------forqt= 1:nr for r = 1 :nc alpha{k}(qt,r) = atan2(d2{k}(qt,r),d1{k}(qt,r)); if alpha{k}(qt,r) < 0 alpha{k}(qt,r)=2*pi+alpha{k}(qt,r); end end end %----------------------------------------------------------%calculate modulus maxima image %----------------------------------------------------------for r=1:mod_m for c=1:mod_n mod_max{k}(r,c)=255; %Initialize modulus maxima array end 55

PAGE 66

end 91a---------------------------------------------------------91afind local maximum of modulus 91a----------------------------------------------------------ang=pi/8; for r=2:mod_m-1 for c=2:mod_n-1 if (alpha{k}(r,c)>=(15*ang)f ... alpha{k}(r,c)<=(1 *ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c) & wtm{k}(r,c) >= wtm{k}(r+ 1,c) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(1 *ang)& .. alpha{k}(r,c)<=(3*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c-1) & wtm{k}(r,c) >= wtm{k}(r+1,c+1) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(3*ang)& .. alpha{k}(r,c)<=(5*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=Wtm{k}(r,c-1) & wtm{k}(r,c) >= wtm{k}(r,c+1) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(5*ang)& .. alpha{k}(r,c)<=(7*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c+1) & wtm{k}(r,c) >= wtm{k}(r+1,c-1) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(7*ang)& .. alpha{k}(r,c)<=(9*ang))& . wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c) & wtm{k}(r,c) >= wtm{k}(r+ 1,c) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(9*ang)& ... 56

PAGE 67

alpha{k}(r,c)<=(11 *ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c-1) & wtm{k}(r,c) >= wtm{k}(r+1,c+1) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(11 *ang)& .. alpha{k}(r,c)<=(13*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r,c-1) & wtm{k}(r,c) >= wtm{k}(r,c+1) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(13*ang)& .. alpha{k}(r,c)<=(15*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c+1) & wtm{k}(r,c) >= wtm{k}(r+1,c-1) mod_max{k}(r,c)=O; end end end end vertical and horizontal detail matrices that include only those ... ... coefficients that are located at the modulus maxima; otherwise, the ... ... coefficients are set to zero. fork= 1:L d1_mm{k} = zeros(nr,nc); d2_mm{k} = zeros(nr,nc); fori= 1:nr for j = 1:nc if mod_max{k}(i,j) == 0 d1_mm{k}(i,j) = d1{k}(i,j); d2_mm{k}(i,j) = d2{k}(i,j); end end end 01_MM(:,:,k) = d1_mm{k}; D2_MM(:,:,k) = d2_mm{k}; end 57

PAGE 68

mm_nonzeros = nnz(D1_MM) + nnz(D2_MM); percent = mm_nonzeros!nonzeros modulus maxima matrix mask figure; imshow(mod_max{1 },[]); xlabei('Modulus Maxima for Level 1 '); figure; imshow(mod_max{2},[]); xlabei('Modulus Maxima for Level 2'); figure; imshow(mod_max{3},[]); xlabei('Modulus Maxima for Level 3'); 58

PAGE 69

Appendix D: MATLAB function to generate modulus maxima masks for details images % wave/et_packet_analysis_of_detail_images (/vl,y,h,g,gp) %This function is called by ca!_mm(scale,/vl,image) to help generate the %modulus maxima matrix masks for detail images. %This function generates the masked horizontal and vertical details along %with the last details %-------------------------------------------------------%Inputs: %-------------------------------------------------------%/vi = number of levels of decomposition %y =original image %h=low pass decomposition filter %g=high pass decomposition filter %gp=high pass reconstruction filter %----------------------------------------------------------%Outputs: %----------------------------------------------------------%y =last level approximation image %01_MM =contains masked horizontal details from every level %02_MM = contains masked vertical details from every level %gprime = reconstruction high pass filter %hprime = reconsrtuction low pass filter %Finally the program displays the modulus maxima matrix masks for each %/eve/. function [y, D1_MM, D2_MM, gprime, wavelet_packet_analysis_ of_ detail_images (lvl,y, h,g,gp) %----------------------------------------------------------%Spline wavelet filters %----------------------------------------------------------hf=double(h) gf=double(g) gprime=double(gp) L=lvl; %L=Ievel of detail [nr,nc]=size(y); nonzeros = L *2*nr*nc; %----------------------------------------------------------59 hprime] =

PAGE 70

%Loop for decomposing image %----------------------------------------------------------fork= 1:L %----------------------------------------------------------%calculates approximations. %shift either image or latest approximation vertically. %----------------------------------------------------------for i=1:nc a{k}(1 :nr,i) = leftshift(y(1 :nr,i)',211(k-1 ))'; %circular convolution of hf with columns from vertically shifted image %----------------------------------------------------------a{k}(1 :nr,i) = cconv(hf,a{k}(1 :nr,in': end %----------------------------------------------------------%shift either image or latest approximation horizontally %----------------------------------------------------------for i=1:nr a{k}(i, 1 :nc) = leftshift(a{k}(i, 1 :nc),211(k-1 )); %----------------------------------------------------------%circular convolution of hf with rows from horizontally shifted image %----------------------------------------------------------a{k}(i, 1 :nc) = cconv(hf,a{k}(i, 1 :nc)); end; %----------------------------------------------------------%calculates horizontal and vertical details. %----------------------------------------------------------%Circular convolution of high pass filter, gf, with columns from image or ... % ... previous approximation %----------------------------------------------------------for i=1:nc d1 {k}(1 :nr,i)=cconv(gf,y(1 :nr,in': %---------------------------------------------------------%Shift resultant image vertically. The result is a horizontal detail ... % ... image %----------------------------------------------------------d1 {k}(1 :nr,i)=leftshift(d1 {k}(1 :nr, i) ',211k) '; end; 60

PAGE 71

91a---------------------------------------------------------91aCircular convolution of high pass filter, gf, with rows from image or ... 91a ... previous approximation 91a----------------------------------------------------------for i=1:nr d2{k}(i, 1 :nc)=cconv(gf,y(i, 1 :nc)); 91a---------------------------------------------------------91aShift resultant image horizontally. The result is a vertical detail ... 0/ /o . 1mage 91a----------------------------------------------------------d2{k}(i, 1 :nc)=leftshift(d2{k}(i, 1 end; d{k}=d1 {k}+d2{k}; 91a---------------------------------------------------------91aDilate filters 91a----------------------------------------------------------hf=[dyadup(hf,2) OJ; 91aa trous hprime=hf; gf=[dyadup(gf,2) 0]; 91aa trous gprime = [dyadup(gprime,2) 0]; y=a{k}; end; 91a---------------------------------------------------------91aSet threshold to zero 91a----------------------------------------------------------Threshold = 0. 005; 91a---------------------------------------------------------91aCalculate the modulus maxima for each level 91a----------------------------------------------------------for k = L:-1:1 wtm{k} = + 91aWavelet Transform Modulus [mod_m,mod_n} = size(wtm{k}); 91a---------------------------------------------------------91acalculate angle of the wavelet transform vector. 91a----------------------------------------------------------for qt = 1:nr for r = 1:nc alpha{k}(qt,r) = atan2(d2{k}(qt,r),d1{k}(qt,r)); if alpha{k}(qt,r) < 0 alpha{k}(qt,r)=2*pi+alpha{k}(qt,r); 61

PAGE 72

end end end modulus maxima image for r=1:mod_m for c=1:mod_n mod_max{k}(r,c)=255; modulus maxima array end end local maximum of modulus ang=pi/8; for r=2:mod_m-1 for c=2:mod_n-1 if (alpha{k}(r,c)>=(15*ang)f ... alpha{k}(r,c)<=(1 *ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c) & wtm{k}(r,c) >= wtm{k}(r+1,c) mod_max{k}(r,c)=D; elseif (alpha{k}(r,c)>=(1 *ang)& .. alpha{k}(r,c)<=(3*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1 ,c-1) & wtm{k}(r,c) >= wtm{k}(r+ 1 ,c+ 1) mod_max{k}(r, c)=D; elseif (alpha{k}(r,c)>=(3*ang)& .. alpha{k}(r,c)<=(S*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r,c-1) & wtm{k}(r,c) >= wtm{k}(r,c+1) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(S*ang)& .. alpha{k}(r,c)<=(l*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c+1) & wtm{k}(r,c) >= wtm{k}(r+1,c-1) mod_max{k}(r,c)=D; 62

PAGE 73

elseif (alpha{k}(r,c)>=(7*ang)& .. alpha{k}(r,c)<=(9*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c) & wtm{k}(r,c) >= wtm{k}(r+ 1,c) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(9*ang)& .. alpha{k}(r,c)<=(11 *ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c-1) & wtm{k}(r,c) >= wtm{k}(r+1,c+1) mod_max{k}(r,c)=O; elseif (alpha{k}(r,c)>=(11 *ang)& .. alpha{k}(r,c)<=(13*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r,c-1) & wtm{k}(r,c) >= wtm{k}(r,c+ 1) mod_max{k}( r, c)=O; elseif (alpha{k}(r,c)>=(13*ang)& .. alpha{k}(r,c)<=(15*ang))& .. wtm{k}(r,c)>= Threshold& .. wtm{k}(r,c)>=wtm{k}(r-1,c+1) & wtm{k}(r,c) >= wtm{k}(r+1,c-1) mod_max{k}(r,c)=O; end end end end 96----------------------------------------------------------96Build vertical and horizontal detail matrices that include only those ... 96 ... coefficients that are located at the modulus maxima; otherwise, the ... 96 ... coefficients are set to zero. 96---------------------------------------------------------for k = 1:L d1_mm{k} = zeros(nr,nc); d2_mm{k} = zeros(nr,nc); fori= 1:nr for j = 1:nc if mod_max{k}(i,j) == 0 d1_mm{k}(i,j) = d1{k}(i,j); d2_mm{k}(i,j) = d2{k}(i,j); end end end 63

PAGE 74

D1_MM(:,:,k) = d1_mm{k}; D2_MM(:,:,k) = d2_mm{k}; end mm_nonzeros = nnz(D1_MM) + nnz(D2_MM); percent = mm_nonzeros/nonzeros 91a---------------------------------------------------------91aDisplay modulus maxima matrix mask 91a----------------------------------------------------------figure; 91a----------------------------------------------------------imshow(mod_max{1 },[]); xlabei('Modulus Maxima for Level 1 J; figure; 91a----------------------------------------------------------imshow(mod_max{2},[]); xlabei('Modulus Maxima for Level 2J; figure; 91a----------------------------------------------------------imshow(mod_max{3},[]); xlabei('Modulus Maxima for Level 3J; 64

PAGE 75

Appendix E: MATLAB function to compute circular convolution function does a circular convolution between an input vector, x, a discrete filter,f. function y = cconv(f,x) m = length(x); r = length(f); if r <= m, length is smaller than the input vector: X_ extended= [x((m+1-r):m) x]; else length is bigger than the input vector: z = zeros(1,r); for i=1 :r, q = r*m -r + i-1; imod = 1 + rem(q,m); z(i) = x(imod); end x_extended = [z x]; end the signal intermediate_y = filter(f, 1 ,x_extended); only the part of the vector that has the circular convolution results: y = intermediate_y(r+ 1 :m+r); 65

PAGE 76

Appendix F: MATLAB function to perform left shift of an input vector, function does a left shift of input vector x by a length LengthOfShift function ShiftedVector=leftshift(x,LengthOfShift) LengthOfVector=length(Vector); fori= 1 :LengthOfShift Temp=Vector(1); Temp Vector= Vector(2:LengthOfVector); Vector=[TempVector Temp]; end ShiftedVector= Vector 66

PAGE 77

REFERENCES Donoho, D., Duncan, M. R., Huo, X., & Levi, 0. (1999). Wavelab 802 for MATLAB 5.x [ Library of MATLAB software programs]. Retrieved from http://www-stat.stanford.edu/-wavelab. Grochenig, K. (1993). Acceleration of the frame algorithm. IEEE Transactions on Signal Processing, 41 (12), 3331-3340. Hsung, T., Lun, D.P., & Siu, W. (1998). A deblocking technique for block transform compressed image using wavelet transform modulus maxima. IEEE Transactions on Image Processing, 7(10), 1488-1496. Hwang, W. L. & Mallat, S. (1992). Singularity detection and processing with wavelets. IEEE Transactions on Information Theory, 38(2), 617-643. Mallat, S. (2001 ). A wavelet tour of signal processing (2nd ed.). New York: Academic Press. Mallat, S. & Zhong, S. (1992). Characterization of signals from multiscale edges. IEEE Transactions of Pattern Analysis and Machine Intelligence, 14(7), 710-732. Polikar, R. (1999). The engineer's ultimate guide to wavelet analysis: The wavelet tutorial. Retrieved October 13, 2002, from the Rowan University College of Engineering Web site: http://engineering.rowan.edu/-polikar/WAVELETSIWTtutorial.html. Shensa, M. (1992). The discrete wavelet transform: Wedding the a trous and Mallat algorithms. IEEE Transactions on Signal Processing, 40(1 0), 2464-2482. 67