Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00002080/00001
## 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:
- 2004
- 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 )
## Auraria Membership |

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. 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.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
(3.4) |

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 |