KING FAHD UNIVERSITY OF PETROLEUM & MINERALS DHAHRAN 31261, SAUDI ARABIA DEANSHIP OF GRADUATE STUDIES This thesis, written by AZZAM ABDULAZIZ ALFARRAJ under the direction of his thesis adviser and approved by his thesis committee, has been presented to and accepted by the Dean of Graduate Studies, in partial fulfillment of the requirements for the degree of MASTER OF SCIENCE IN MATHEMATICS. Thesis Committee Dr. Mohamed El-Gebeily (Adviser) Dr. Faisal Fairag (Member) Dr. Abdelmalek Zidouri (Member) Dr. Husain Al-Attas Department Chairman Dr. Salam A. Zummo Dean of Graduate Studies Date i ii ➞Azzam Alfarraj 2017 iii Z@YëB @ ..ø Yg. úÍ@ ..h. @Q®Ë@ YÔg áK. YÒm× ..úGAJk ú¯ Èð B@ XAJƒ B@ iv Dedication To my grandfather.. Mohammad Hamad Alfarraj.. My first mentor.. v ACKNOWLEDGMENTS First of all, I thank Allah for all what I have in my life that make it easy for me to finish this degree with easiness. I thank my parents, my wife, my siblings and all members of my family for their unlimited support. Great thanks go to my advisor Dr. Mohamed El-Gebeily for guiding me to be a good Mathematician, for his efforts and for his great patience. I thank Dr. Faisal Fairag for his help throughout my work. I also thank him and Dr. Abdelmalek Zidouri for accepting to be in the committee of my defence. I really thank my friends Jaafar Motasim, Jamilu Hashim, Mohammad Al-Shehri, Mohsen Al-Shahrani for their help throughout my way to the master’s degree. I thank the chairman Dr. Hussain Al-Attas for his efforts and unlimited help. vi TABLE OF CONTENTS ACKNOWLEDGEMENT vi LIST OF FIGURES ix ABSTRACT (ENGLISH) xi ABSTRACT (ARABIC) xii CHAPTER 1 INTRODUCTION 1 CHAPTER 2 LITERATURE REVIEW 3 2.1 Gradient-based Methods . . . . . . . . . . . . . . . . . . . . . . . 3 2.2 Minimization-based Methods . . . . . . . . . . . . . . . . . . . . . 5 2.3 Wavelet Transform-based Methods . . . . . . . . . . . . . . . . . 6 CHAPTER 3 PRELIMINARIES 8 3.1 Wavelet Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 3.2 Vector Valued Inner Products . . . . . . . . . . . . . . . . . . . . 9 3.3 Function Spaces . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.4 Convergence of Sequences of Functions . . . . . . . . . . . . . . . 13 3.5 The Haar System . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.5.1 Dyadic Intervals . . . . . . . . . . . . . . . . . . . . . . . . 14 3.5.2 The Haar Scaling Function and the Haar Wavelet Function 15 3.5.3 Orthogonality of The Haar System . . . . . . . . . . . . . 20 Multiresolution Analysis (MRA) . . . . . . . . . . . . . . . . . . . 22 3.6 vii CHAPTER 4 DISCRETE WAVELET TRANSFORM IN ONE AND TWO DIMENSIONS 25 4.1 Discrete Wavelet Transform in One Dimension . . . . . . . . . . . 26 4.1.1 Fast Wavelet Transform in One Dimension . . . . . . . . . 34 Discrete Wavelet Transform in Two Dimensions . . . . . . . . . . 37 4.2.1 MRA in L2 [0, 1]2 . . . . . . . . . . . . . . . . . . . . . . . 41 4.2.2 Fast Wavelet Transform in Two Dimensions . . . . . . . . 43 4.2 CHAPTER 5 EDGE DETECTION 46 5.1 Detecting Singularities of Functions in One Dimension . . . . . . 47 5.2 Edge Detection Using Discrete Wavelet Transform in Two Dimensions 56 5.3 Edge Detection Using Canny edge detector . . . . . . . . . . . . . 59 5.4 Results and Comparison of Images Without Noise . . . . . . . . . 62 5.5 Results and Comparison of Images with Noise . . . . . . . . . . . 76 5.6 Comparison of Haar with db2 and db3 in edge detection . . . . . 82 5.7 Edge Detection In Colored Images . . . . . . . . . . . . . . . . . . 86 REFERENCES 92 VITAE 96 viii LIST OF FIGURES 3.1 The Haar scaling function . . . . . . . . . . . . . . . . . . . . . . 3.2 (a) The Haar scaling function at level j = 2 and k = 0. (b) j = 2 and k = 1. (c) j = 2 and k = 2. (d) j = 2 and k = 2 15 . . . . . . . 16 3.3 The Haar wavelet function . . . . . . . . . . . . . . . . . . . . . . 19 4.1 The decomposition of an input vector cj . . . . . . . . . . . . . . . 35 4.2 The reconstruction of a vector cj+1 . . . . . . . . . . . . . . . . . . 36 5.1 (a) The function f (x), (b) The detail coefficients for j = 1, (c) The detail coefficients for j = 2 and (d) The detail coefficients for j = 3 53 5.2 Flow chart for wavelet-based singularity detection algorithm . . . 55 5.3 Flow chart for wavelet-based edge detection algorithm . . . . . . . 58 5.4 Flow chart for Canny edge detector algorithm . . . . . . . . . . . 61 5.5 Pepsi can image . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 5.6 Pepsi can edges using dwt2 with Haar . . . . . . . . . . . . . . . . 63 5.7 Pepsi can edges using Canny edge detector . . . . . . . . . . . . . 64 5.8 Times square image . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.9 Times square edge image using dwt2 with Haar . . . . . . . . . . 66 5.10 Times square edge image using Canny edge detector . . . . . . . . 67 5.11 Bird image . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 5.12 Bird edge image using dwt2 with Haar . . . . . . . . . . . . . . . 69 5.13 Bird edge image using Canny edge detector . . . . . . . . . . . . . 70 5.14 Dr Pepper can image . . . . . . . . . . . . . . . . . . . . . . . . . 71 5.15 Dr Pepper can edge image using dwt2 with Haar . . . . . . . . . . 72 ix 5.16 DrPepper can edge image using Canny edge detector . . . . . . . 73 5.17 Toys room image . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 5.18 Toys room edge image using dwt2 with Haar . . . . . . . . . . . . 75 5.19 Toys room edge image using Canny edge detector . . . . . . . . . 76 5.20 Noised Pepsi can image . . . . . . . . . . . . . . . . . . . . . . . . 77 5.21 Noised Pepsi can edge image using dwt2 with Haar . . . . . . . . 78 5.22 Noised Pepsi can edge image using Canny edge detector . . . . . . 79 5.23 Noised Times square image . . . . . . . . . . . . . . . . . . . . . . 80 5.24 Noised Times square edge image using dwt2 with Haar . . . . . . 81 5.25 Noised Times square edge image using Canny edge detector . . . 82 5.26 Pepsi can image . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 5.27 Pepsi can edge image using dwt2 with Haar . . . . . . . . . . . . 84 5.28 Pepsi can edge image using dwt2 with db2 . . . . . . . . . . . . . 85 5.29 Pepsi can edge image using dwt2 with db3 . . . . . . . . . . . . . 86 5.30 Colored Pepsi can image . . . . . . . . . . . . . . . . . . . . . . . 88 5.31 Colored Pepsi can edge image using dwt2 with Haar . . . . . . . . 88 5.32 Colored Pepsi can edge image using Canny edge detector . . . . . 89 5.33 Colored toys room image . . . . . . . . . . . . . . . . . . . . . . . 89 5.34 Colored toys room edge image using dwt2 with Haar . . . . . . . 90 5.35 Colored toys room edge image using Canny edge detector . . . . . 90 x THESIS ABSTRACT NAME: Azzam Abdulaziz Alfarraj TITLE OF STUDY: Edge Detection Using Discrete Wavelet Transform MAJOR FIELD: Mathematics DATE OF DEGREE: December 2017 The edge detection problem plays an important role in many applications. It helps in extracting the main features of an image and specifying its constituent parts. In the literature, this problem has been tackled using different mathematical approaches such as gradient-based detectors and minimization-based approach. Also, the continuous wavelet transform approach is involved in many edge detection algorithms. The cornerstone in all approaches is to look for points where the intensity of the image has a jump discontinuity. Applying the discrete wavelet transform in two dimensions to an image will result in large detail coefficients at the parts that have edges. The objective of this thesis is to propose a rigorously analysed algorithm for edge detection using discrete wavelet transform in two dimensions. To the best of our knowledge, this algorithm has not been treated with rigor in the literature. xi é“CmÌ '@  éË A‚Ó  k , èYK Y« HA  ®J J.¢ ú¯ AÒêÓ @P ðX I.ªÊK èPñ’Ë@ ¬@ ñk ­‚» ú¯ Y«A‚ IJ  KQË@ HA  K. X @ ú¯ .AîD¯ é‚  KñºÖÏ @ YK Ym' ð èPñ’Ë@ ú¯ éJ ƒAƒ B@ @ñmÌ '@ CjJƒ@ è Yë HAJ    èY« Ð@YjJƒAK Im  ®ƒA¾Ë@  ÉJÓ éJ “AK P HAK  . PA®Ó úΫ èYÒJªÖÏ @ HA .  .Ì 'ñ« à @ ‡J.ƒ , éË A‚ÖÏ @  IÓY  jJƒ@ ½Ë Y» .ÉJÊ® JË@ úΫ èYÒJªÖÏ @ ð h. PYJË@ QKñÓ  .' ñÖÏ @ ÉK ñm' éK . PA®Ó ɒJÖÏ @ HAm  . Ë@ ñë HAK  . PA®Ö Ï @ ½ÊK É¿ ú¯ éK ð@QË@ Qm.k . ¬@ñmÌ '@ ­‚»  HAJ  ÓPP@ñk áÓ XY« ú¯ á« IjJ  A® JË@  .' ñÖÏ @ ÉK ñm' AJ®J . £ ñË .XAg ɾ‚  . èPñ’Ë@ èYƒ AëYJ« QªJK úæË@ ©¢® JÖÏ @ HAm  Z@Qk. B@ ú¯ èQJ.» ÉJ“A®K HCÓAªÓ  ø ñm' úæË@ úΫ ɒjJƒ AJKA¯ AÓ èPñ“ úΫ áK YªJ.Ë     Qå” é<Êm    éÓ@ ¬@ñmÌ '@ ­‚ºË . × éJÓPP@ñk Õç' Y®K ñë éËAƒQË@ è Yë ú¯ ¬YêË@ .A¯@ñk  K AÓ ú支 @ I.‚k . áK YªJ.Ë ©¢® JÖÏ @ HAm  .' ñÖÏ @ ÉK ñm' Ð@YjJƒAK. è Yë àA¯ éJË@ Iʓñ  Qå” IÊÔ    . Ë@ @ Yë HAJ  K. X @ ú¯ éÓ@ . IjJ .  « à @ ‡J.‚ ÕË éJÓPP@ñmÌ '@ xii CHAPTER 1 INTRODUCTION The edge detection problem plays an essential role in many applications. In fact, it is a field on its own in image processing [1]. Moreover, it is the first step in many computer vision applications, [2]. It is used mainly for feature extraction [3] or for image segmentation [4], [5]. Edges identify the significant information in an image, which enable the machine to separate its constituent parts. Additionally, it reduces the amount of data stored and preserves the main structural properties of the image for further processing [1], [6]. Other applications include medical imaging [7], [8], face recognition, automatic control of traffic systems, obstacle recognition, etc. [1], [2] Before we start, it is worth noting that, mathematically speaking, an image is a function of two variables u(x, y) defined on a rectangle. If grey-scale images are considered, as in this work, u(x, y) is a grey-scale level that represents the intensity of the image color at the point (x, y), rescaled to be a value in an interval [a, b] for prespecified real numbers a and b. This rescaling process is admissible since 1 images are bounded functions [9]. The cornerstone in edge detection techniques and algorithms is to look for points where the image intensity changes sharply. An edge is a sign for discontinuity [4]. Roughly speaking, it is a segment of connected pixels separating two different regions in the image [1], [4]. Indeed, this definition of an edge is not mathematically rigorous and it may lead to detecting false edges as in the case of textures [9]. This is why in the very early work on this problem, researchers only considered homogeneous (without texture) images [10]. The main challenges in edge detection include detecting false edges, missing edges in low contrast images, high computational time, sensitivity of noise, missing edges of objects with small change of intensity [2]. Amazingly, some detectors face problems in detecting edges of high-definition images. The algorithm proposed in this thesis performs better than most of the edge detectors when applied to high-definition images. The discrete wavelet transform can be applied using plenty of wavelet families such as Haar, Daubechies, Mexican hat, etc. In this thesis we will focus on the Haar wavelet which produced the best results in detecting edges. This is because the abrupt change in its definition matches the discontinuity represented by an edge in an image. Furthermore, it does not mollify the image as other smoother wavelets do. 2 CHAPTER 2 LITERATURE REVIEW 2.1 Gradient-based Methods Gradient-based methods for detecting edges were commonly used in applications since 1960’s. Roberts [11] presented an operator in his Ph.D thesis where he used the following two 2 × 2 kernel matrices    1 0   Dx =    0 −1   0 1  Dy =    −1 0 to compute the spatial gradient components of the image Gx and Gy . Gx and Gy stand for the difference in the x-direction and y-direction respectively. These components are used to get the gradient magnitude and angle as follows: |G| = q G2x + G2y θ = arctan( 3 Gx 3π )− . Gy 4 Then some threshold T is set to mark the edge pixels wherever |G| > T . The angle θ represents the direction perpendicular to the edge line. For further improvements, Sobel [12] and Prewitt [13] presented two kernels as 3 × 3 matrices as follows: Sobel kernel matrices:   −1 0 1      Dx = −2 0 2      −1 0 1     2 1 1      Dy =  0 0 0      −1 −2 −1 Prewitt kernel matrices:   1 1 1     Dy =  0 0 0      −1 −1 −1 −1 0 1      Dx =  −1 0 1     −1 0 1 The gradient magnitude is computed as in Robert’s while the gradient angle is θ = arctan( Gx ). Gy Anyway, all three operators were sensitive to noise. In 1980, Marr and Hildreth [14] introduced an operator that convolves the image u with a Gaussian function G and then applies the Laplacian operator ∇2 to 4 the convolved image. The algorithm then looks for the points where ∇2 (u∗G) = 0 to identify points of local extrema of the gradient. The idea was dominating at that time [1]. In 1986, John Canny presented an approach where edges are marked as maxima in gradient magnitude of a Gaussian-smoothed image [15]. The image u is convolved with a Gaussian function G then the gradient magnitude |∇(u ∗ G)| is computed to detect edges. Bergholm [16] used Canny’s operator and he succeeded to differentiate between shadow contours and perfect ones [17]. Due to its great advantages, Canny edge detector is still used widely after some improvements in nowadays image processing techniques [1]. 2.2 Minimization-based Methods Mumford and Shah introduced a different approach to image segmentation in 1989 [18]. It is based on a minimization problem. They look for a pair (u, K) that minimizes the following functional F (u, K) = Z 2 Ω\K (u − u0 ) dx + α Z 2 Ω\K |∇u| dx + β Z dσ K where Ω = [0, 1]2 , u0 the initial image, K ⊂ Ω is a set of discontinuities and α and β are non-negative parameters. In fact, if we drop the second term and reject 5 the trivial solution, we would get the reduced Mumford and Shah functional E(u, K) = Z 2 Ω\K (u − u0 ) dx + β Z dσ. K This particular case gained some interest and Mumford and Shah [18], Morel and Solimini [19] [20], and Massari and Tamanini [21] devoted some work to it. Indeed, this approach will segment the image into its constituent pieces and K will present the edges in the image. 2.3 Wavelet Transform-based Methods Stephane Mallat et al. [22] implemented Canny’s approach with the conjugate of a convolution kernel θ(x, y), instead of a Gaussian function. Practically, it was computed with two wavelets ψ 1 and ψ 2 that are partial derivatives of the convolution kernel θ [9] ψ1 = − ∂θ ∂x and ψ2 = − ∂θ . ∂y The algorithm computes the wavelet transform components of an image and the result is proportional to the gradient sought after. More explicitly, W u(x, 2j ) = 2j ∇(u ∗ θ̄), 6 where W u(a, 2b ) is the continuous wavelet transform of a function f at scale b and location a. In 2016, Monika published a paper [23] discussing image edge detection using the discrete wavelet transform. The paper was descriptive and it was not clear what algorithm was used. Indeed, she talked about a threshold to get the high detail coefficients but she never specified the threshold or even the way how to arrive at it. In 2005, Chaganti [17] described how to use four kinds of wavelets, Haar, Daubechies, Coifman and biorthogonal, to detect edges. His thesis was also descriptive, as he was an engineering student. Moreover, he did not state a clear algorithm. 7 CHAPTER 3 PRELIMINARIES In this chapter, we will go through many notions and results needed in this thesis. First, we will give the definition of a wavelet. Then, we will give some results of vector valued inner products. Also, we will go through some function spaces needed in our work and the convergence of sequences of functions in these spaces. After that, we will give a brief about the Haar system since our work depends heavily on the Haar system. Finally, we will talk about the notion of multiresolution analysis. 3.1 Wavelet Theory Definition 3.1 (Wavelet) A function ψ(x) is called a wavelet if it satisfies the following: Z |ψ(t)|2 dt < ∞ R 8 (3.1) Z ψ(t)dt = 0 (3.2) R 3.2 Vector Valued Inner Products In this section we let H be a Hilbert space. We will show some definitions and properties of inner products of vectors in H. Now we give the definition of vector valued inner products in H. Definition 3.2 Let F and G be two vectors of elements in H defined as follows:     f1      f   2  F= .  ..        fn ,   g1      g   2  G=  .   ..        gm where fi , gj ∈ H for 1 6 i 6 n and 1 6 j 6 m, then    hf1 , g1 i hf1 , g2 i . . . hf1 , gm i       hf , g i hf , g i . . . hf , g i   2 1 2 2 2 m   hF, Gi =    . .. .. ..   .. . . .       hfn , g1 i hfn , g2 i . . . hfn , gm i (3.3) In fact, (3.3) is obtained by forming (formally) the matrix FG∗ and then applying the inner product to each element of the resulting matrix. Lemma 3.1 The vector valued inner products have the following three properties: 9 (1) hF + E, Gi = hF, Gi + hE, Gi (2) hAF, BGi = AhF, GiB∗ , where A and B are complex matrices and B∗ is the conjugate transpose of B (3) hF, Gi∗ = hG, Fi Definition 3.3 Given a vector G = {g1 , g2 , ..., gn }, we call hG, Gi the Grammian of G. Theorem 3.1 Suppose G = {g1 , g2 , ..., gn }T is a set of vectors in an inner product space H. Then hG, Gi is invertible iff g1 , g2 , ..., gn are linearly independent. Proof. (=⇒:) Assume that hG, Gi is invertible and g1 , g2 , ..., gn are not linearly independent. ∴ ∃ a vector c 6= 0 such that c∗ G = 0 ⇒ hG, c∗ Gi = 0 ⇒ hG, Gic = 0 ⇒ c = 0 which contradicts the assumption. (⇐=:) Assume that g1 , g2 , ..., gn are linearly independent and hG, Gi is not invertible. ∴ ∃ a vector c 6= 0 such that hG, Gic = 0 ⇒ hgi , c∗ Gi = 0, i = 1, 2, ..., n ⇒ c∗ G ⊥ span(G). On the other hand, c∗ G ∈ span(G) ∴ c∗ G = 0 10 Then either c∗ = 0 or G is not linearly independent, which both contradict the assumption. 3.3 Function Spaces In this thesis we will come across many function spaces and I will give a brief definition of each one of them. Definition 3.4 The space of all functions f : Ω −→ R; such that f is continuous on Ω is denoted by C(Ω) and called the space of continuous functions. Definition 3.5 (See [24]) For 1 6 p < ∞, we define Lp (Ω) to be the collection of all functions f for which Z Ω |f (x)|p dx < ∞. Moreover, the Lp -norm of a measurable function f is defined by kf kp = Z p Ω |f (x)| dx 1/p Definition 3.6 (See [24]) A function f is called essentially bounded if there is some M > 0, called the essential upper bound for f , for which |f (x)| 6 M for almost all x ∈ Ω. Definition 3.7 (See [24]) We define L∞ (Ω) to be the collection of all essentially 11 bounded functions. The L∞ -norm of a measurable function f(x) is defined by kf k∞ = ess sup{|f (x)| : x ∈ Ω}. The space L2 (Ω) is a Hilbert space equipped with the inner product hf (x), g(x)i = Z f (x)g(x)dx, f, g ∈ L2 (Ω) Ω where g(x) is the complex conjugate of g(x). Definition 3.8 (Support of a Function) Suppose that f : Ω −→ R, then the support of f(x), denoted by supp(f ), is the closure of the set of all vectors x such that f(x) 6= 0. i.e., supp(f ) = {x ∈ Ω : f (x) 6= 0}. Definition 3.9 (Compact Support) A function f is said to be of compact support if its support can be included in a closed and bounded domain. Definition 3.10 A function f(x) defined on an Ω is C n on Ω with n ∈ N if it is n-times continuously differentiable on Ω. C 0 on I if it is continuous on Ω. f(x) is C ∞ on Ω if it is C n on Ω for all n ∈ N. A function f(x) is Ccn on Ω if it is C n on Ω and compactly supported. Theorem 3.2 (See [25]) If f ∈ L2 (R) and ǫ > 0, then there exists a function g and R > 0, such that R R |f − g|2 dx < ǫ with g(x) = f (x)χ [−R,R] (x). This means that most of the energy of f is in a finite interval. 12 3.4 Convergence of Sequences of Functions Definition 3.11 (Pointwise convergence) Given a function f(x) and a sequence of functions {fn (x)}n∈N defined on an interval I, we say that the sequence {fn (x)}n∈N converges pointwise to f(x) if for every x0 ∈ I the sequence {fn (x0 )}n∈N converges to f (x0 ). We write fn (x) → f (x) pointwise as n → ∞. Definition 3.12 (Uniform (L∞ ) convergence) Given a function f(x) and a sequence of functions {fn (x)}n∈N defined on an interval I, we say that the sequence {fn (x)}n∈N converges uniformly to f(x) if for every ǫ > 0, there exists a number N > 0 such that if n > N, then |fn (x) − f (x)| < ǫ for every x ∈ I. We write fn (x) → f (x) uniformly as n → ∞. Definition 3.13 (Mean-square (L2 ) convergence) Given a function f(x) and a sequence of functions {fn (x)}n∈N defined on an interval I, we say that the sequence {fn (x)}n∈N converges in the mean-square to f(x) if lim n→∞ Z I |fn (x) − f (x)|2 dx = 0. We write fn (x) → f (x) in the mean-square as n → ∞. Theorem 3.3 (See [25]) Given that I is a finite interval, (a) L∞ (I) ⊂ L2 (I). (b) if fn (x) → f (x) uniformly on I, then fn (x) → f (x) in the mean-square on I. 13 Theorem 3.4 (See [25]) Given that I is a finite interval, if fn (x) → f (x) in L∞ (I) or L2 (I), then fn (x) → f (x) pointwise on I. 3.5 The Haar System Alfred Haar, the Hungarian mathematician, itroduced the Haar function in his Ph.D thesis in 1910 [26]. In this section, we will discuss the Haar function and Haar scaling function with some notions, lemmas and theorems that relate to them. First, we will introduce the notion of dyadic intervals. 3.5.1 Dyadic Intervals Definition 3.14 For each pair j, k ∈ Z, the interval Ij,k = [2−j k, 2−j (k + 1)] is called a dyadic interval. The collection of such intervals is known as the collection of dyadic sub-intervals of R. The parameter j will be referred to as the scale level or scale parameter, or simply ”scale”. Note that if j 6= j1 or k 6= k1 then either (i) (ii) Ij,k ∩ Ij1 ,k1 = φ Ij,k ⊆ Ij1 ,k1 or Ij1 ,k1 ⊆ Ij,k (3.4) (3.5) l r Definition 3.15 Given a dyadic interval Ij,k at scale j, we write Ij,k = Ij,k ∪ Ij,k , l r where Ij,k and Ij,k denote the left half and right half of the interval Ij,k and they l r are dyadic intervals at scale j + 1. Indeed, Ij,k = Ij+1,2k and Ij,k = Ij+1,2k+1 . 14 3.5.2 The Haar Scaling Function and the Haar Wavelet Function The Haar scaling function is defined as ϕH (x) := χ[0,1) (x) =     1 if   0, 06x<1 otherwise. Figure 3.1 shows the graph of the Haar scaling function. Figure 3.1: The Haar scaling function j/2 H j Definition 3.16 Let ϕH ϕ (2 x − k) with j, k ∈ Z. j,k (x) = 2 The collec- tion {ϕH j,k (x)}j,k∈Z is called the Haar system of scaling functions. The collection {ϕH j,k (x)}k∈Z is called the system of Haar scaling functions at scale j. 15 Figure 3.2 shows the graphs of the Haar scaling functions at level j = 2. Figure 3.2: (a) The Haar scaling function at level j = 2 and k = 0. (b) j = 2 and k = 1. (c) j = 2 and k = 2. (d) j = 2 and k = 2 Lemma 3.2 For all j, k ∈ Z, H kϕH j,k (x)k2 = kϕ (x)k2 . Proof. 2 kϕH j,k (x)k2 = Z∞ 2 |ϕH j,k (x)| dx = dt = 2j dx , 2 kϕH j,k (x)k2 = Z∞ j [2j/2 ϕH (2j x − k)]2 dx −∞ −∞ Let t = 2j x − k Z∞ H then 2 −j 2 [ϕ (t)] 2 dt = Z∞ −∞ −∞ 16 [ϕH (t)]2 dt = kϕH (x)k22 The system of Haar scaling functions satisfies the following [25] [27]: (i) ϕH j,k (x) is supported on the interval Ij,k . (ii) For every j, k ∈ Z, (iii) For every j, k ∈ Z, R R R 2 |ϕH j,k (x)| dx = 1. −j/2 ϕH . j,k (x)dx = 2 R j/2 H j Lemma 3.3 If ϕH ϕ (2 x − k) with j, k ∈ Z, and supp(ϕH ) = [a, b] j,k (x) = 2 then supp(ϕH j,k ) = Ij,k = [ a+k b+k , j ] 2j 2 Proof. If supp(ϕH ) = [a, b], then we will look for values of x that will make the expression (2j x − k) fall within the interval [a,b]. 2j x − k > a ⇒ x> a+k 2j 2j x − k 6 b ⇒ x6 b+k . 2j Then, supp(ϕH j,k ) = Ij,k = [ Moreover, |Ij,k | = a+k b+k , j ]. 2j 2 b−a 2j We now address the question of which Haar scaling functions have support in [0,1]. Lemma 3.4 The Haar scaling functions at scale j with support in [0, 1] are exj 2 −1 actly {ϕH j,k (x)}k=0 . 17 Proof. First, keep in mind that in the interval [0, 1] we have Ij,k = 2−j [k, k + 1] 2−j k > 0 2−j (k + 1) 6 1 ⇒ ⇒ k>0 (k + 1) 6 2j k 6 2j − 1 ⇒ The Haar wavelet function (or the Haar function) is defined as      1      H ψ (x) = χ[0,1/2) (x) − χ[1/2,1) (x) = −1         0, if if 0 6 x < 1/2 1/2 6 x < 1 otherwise. H Definition 3.17 The collection {ψj,k (x)}j,k∈Z is called the Haar wavelet system. H The collection {ψj,k (x)}k∈Z is called the system of Haar wavelet functions at scale j. Figure 3.3 shows the graph of the Haar wavelet function. 18 Figure 3.3: The Haar wavelet function Note that lemmas 3.2, 3.3 and 3.4 also hold for Haar wavelet functions. The Haar wavelet system satisfies the following [25] [27]: H (i) ψj,k (x) is supported on the interval Ij,k . (ii) For every j, k ∈ Z, (iii) For every j, k ∈ Z, R R R R H |ψj,k (x)|2 dx = H ψj,k (x)dx = R Ij,k R Ij,k H (iv) ψj,k (x) = 2j/2 (χIj,k l (x) − χI r (x)). j,k 19 H |ψj,k (x)|2 dx = 1. H ψj,k (x)dx = 0. 3.5.3 Orthogonality of The Haar System Definition 3.18 (Orthogonal Basis and Orthonormal Basis) A basis {gn (x)}∞ n=1 of a Hilbert space X is orthogonal if hgn , gm i = 0 f or n 6= m. An orthogonal basis is called orthonormal if hgn , gm i = δnm , where δnm is the kronecker δ; δnm =     1 n=m    0, n 6= m. Theorem 3.5 The system of Haar wavelet functions(or the Haar System) H {ψj,k (x)}j,k∈Z is orthonormal on R. Proof. The inner product we have is H hψj,k , ψjH1 ,k1 i = Z∞ H ψj,k (x)ψjH1 ,k1 (x)dx. (3.6) −∞ We show the orthonormality within a scale j first. The properties of dyadic H H intervals tell us that if k 6= k1 then ψj,k ψj,k1 ≡ 0, and thus H H hψj,k , ψj,k i = 0. 1 20 If k = k1 , then H H hψj,k , ψj,k i 1 = Z∞ H H ψj,k (x)ψj,k (x)dx 1 −∞ = Z H |ψj,k (x)|2 dx = 1 Ij,k Now, we show the orthonormality in case of different scales. Without loss of generality, say that j > j1 , then the properties (3.4) and (3.5) give us the following three possibilities: (i) Ij,k ∩ Ij1 ,k1 = φ which results in H hψj,k , ψjH1 ,k1 i = Z∞ H ψj,k (x)ψjH1 ,k1 (x)dx = 0. −∞ (ii) Ij,k ⊆ Ijl1 ,k1 . In this case, ψjH1 ,k1 (x) is identically 1 on Ijl1 ,k1 . This results in H hψj,k , ψjH1 ,k1 i = Z H ψj,k (x)ψjH1 ,k1 (x)dx = Z H ψj,k (x)dx = 0. Ij,k Ij,k (iii) Ij,k ⊆ Ijr1 ,k1 . In this case, ψjH1 ,k1 (x) is identically -1 on Ijr1 ,k1 and on Ij,k . This results in H hψj,k , ψjH1 ,k1 i = Z H ψj,k (x)ψjH1 ,k1 (x)dx Ij,k =− Z H (x)dx = 0. ψj,k Ij,k Theorem 3.6 For any j ∈ Z, the system of Haar scaling functions at scale j is an orthonormal system on R 21 Proof. First, If k = k1 then H hϕH j,k , ϕj,k i = Z∞ H ϕH j,k (x)ϕj,k (x)dx = Z∞ 2 |ϕH j,k (x)| dx = 1. −∞ −∞ H Second, if k 6= k1 then the properties of dyadic intervals tell that ϕH j,k ϕj,k1 ≡ 0, and thus H hϕH j,k , ϕj,k1 i = Z∞ H ϕH j,k (x)ϕj,k1 (x)dx = 0. −∞ H Theorem 3.7 ([25]) For any J ∈ Z, the collection {ϕH J,k (x), ψj,k (x) : j > J, k ∈ Z} is an orthonormal system on R 3.6 Multiresolution Analysis (MRA) H (x)}j,k∈Z = {2j/2 ψ H (2j x − k)}j,k∈Z We saw in Section 3.5 that the collection {ψj,k forms an orthonormal system on R. In this section, we will see how to generalize this construction. In particular, we will see a general framework for constructing wavelet functions ψ(x) such that the collection {ψj,k (x)}j,k∈Z = {2j/2 ψ(2j x − k)}j,k∈Z is an orthonormal system on R. Such collection {ψj,k (x)}j,k∈Z is called a wavelet orthonormal basis on R. Definition 3.19 (Uniformly stable basis) For a given scale j, the basis Φj = 22 {ϕj,k (x) : k ∈ Z} is a uniformly stable basis for the space Sj if there exist two real numbers α and β independent of j and α, β > 0 such that αkck2 6 k X k∈Z cj,k ϕj,k k2 6 βkck2 Definition 3.20 (Multiresolution analysis MRA [28]) A multiresolution analysis is a sequence S = {Sj }j∈Z of subspaces of L2 (R) that satisfies the following properties: (i) The spaces are nested. i.e. Sj ⊂ Sj+1 for every j ∈ Z. (ii) Their union is dense in L2 (R). i.e. (iii) Their intersection is trivial. i.e. T S j∈Z j∈Z Sj = L2 (R). Sj = {0}. (iv) There is a function ϕ such that each Φj = {ϕj,k (x) : k ∈ Z} is a uniformly stable basis for Sj . (v) The spaces arises by scaling. i.e. ϕ ∈ Sj ⇒ ϕ(2·) ∈ Sj+1 . (vi) They are shift-invariant. i.e. ϕ ∈ S0 ⇒ ϕ(· − k) ∈ S0 , k ∈ Z. Definition 3.21 ([28]) For a MRA S, any function ϕ ∈ L2 (R) that satisfies ϕ(x) = X k∈Z ak ϕ(2x − k) (3.7) is called a scaling function of the MRA S. Moreover, the equation (3.7) is called 23 the refinement equation and the sequence {ak }k∈Z is called the refinement coefficients (or refinement mask). 24 CHAPTER 4 DISCRETE WAVELET TRANSFORM IN ONE AND TWO DIMENSIONS In this chapter, we will give the definitions and properties of spaces and operators needed in discrete wavelet transform in one dimension and then we will build on the one-dimension case to go to the two-dimension case. We will discuss the notions of approximation space and detail space and the operators associated with them. Indeed, our interest is in discrete wavelet transform in two dimensions since the images are two dimensional functions. All the work done in this chapter is over the interval [0,1] rather than R since our scope is images which are bounded in domain. Throughout this chapter, we will consider only the Haar system. 25 4.1 Discrete Wavelet Transform in One Dimension Let Sj = span({ϕj,k (x)}k∈Z ) with ϕj,k (x) = 2j/2 ϕ(2j x − k). In discrete wavelet transform in one dimension, we project L2 [0, 1] functions into the approximation space Sj . Pj : L2 [0, 1] −→ Sj hPj f, gi = hf, gi ∀g ∈ Sj (4.1) We can easily show the following properties of the operator Pj (i) Pj is linear. i.e. Pj (αf + h) = αPj f + Pj h. (ii) Pj is an orthogonal projection. i.e. hg, Pj f − f i = 0 ∀g ∈ Sj We will now show how to derive an explicit formula for the operator Pj . Given that Φj = {ϕj,k (x)}k∈Z is a basis for the space Sj , we have the following Pj f = c∗ Φj hf, Φj i = hPj f, Φj i = hc∗ Φj , Φj i = c∗ hΦj , Φj i So, since hf, Φj i = c∗ hΦj , Φj i then, using Theorem 3.1 we get c∗ = hf, Φj ihΦj , Φj i−1 . If Φj is an orthogonal basis, as in the case of Haar system, then 26 hΦj , Φj i−1 = I and c∗ = hf, Φj i. Definition 4.1 (The approximation operator) For every j ∈ Z, the approximation operator Pj of a function f ∈ L2 [0, 1] is Pj f = hf, Φj iΦj Φj is orthogonal If we write the refinement coefficients in Definition 3.21 in a matrix a with elements defined as follows: 1 ak,m = √ am−2k 2 (4.2) we will be able to go from a finer level to a coarser level of approximation subspaces by the relation Φj = aΦj+1 . (4.3) In the following proposition, we will show more properties of the operator Pj . Proposition 4.1 The operator Pj given in Definition 4.1 has the following properties: (i) P2j = Pj . (ii) Pj+1 Pj = Pj Pj+1 = Pj (iii) hPj f, gi = hf, Pj gi ∀f, g ∈ L2 [0, 1] Proof. 27 (i) P2j f = Pj Pj f = Pj hf, Φj iΦj = hf, Φj iPj Φj = hf, Φj iΦj (ii) Pj+1 Pj f = Pj+1 hf, Φj iΦj = hf, Φj iPj+1 Φj = hf, Φj ihΦj , Φj+1 iΦj+1 = hf, Φj ihaΦj+1 , Φj+1 iΦj+1 using the refinement relation (4.7) = hf, Φj iahΦj+1 , Φj+1 iΦj+1 = hf, Φj iaΦj+1 = hf, Φj iΦj = Pj f Pj Pj+1 f = Pj hf, Φj+1 iΦj+1 = hf, Φj+1 iPj Φj+1 = hf, Φj+1 ihΦj+1 , Φj iΦj = hf, Φj+1 ihΦj+1 , aΦj+1 iΦj = hf, Φj+1 ihΦj+1 , Φj+1 ia∗ Φj = hf, Φj+1 ia∗ Φj = hf, aΦj+1 iΦj = hf, Φj iΦj = Pj f 28 using the refinement relation (4.7) (iii) hPj f, gi = hhf, Φj iΦj , gi = hf, Φj ihΦj , gi = hf, hg, Φj iΦj i = hf, Pj gi Definition 4.2 (The detail operator) For every j ∈ Z, the detail operator Qj is defined as Qj = Pj+1 − Pj . (4.4) Moreover, let Wj = Im(Qj ) and it is called the detail subspace. In the following proposition, we will show some properties of the operator Qj . Proposition 4.2 The operator Qj given in Definition 4.2 has the following properties: (i) Q2j = Qj . (ii) Qj Pj = Pj Qj = 0 (iii) hQj f, gi = hf, gi (iv) hQj f, gi = hf, Qj gi (v) Qj g = g ∀g ∈ Wj ∀f, g ∈ L2 [0, 1] ∀g ∈ Wj 29 Proof. (i) Q2j = (Pj+1 − Pj )2 = P2j+1 − Pj+1 Pj − Pj Pj+1 + P2j = Pj+1 − Pj − Pj + Pj using Proposition 4.1(ii) = Pj+1 − Pj = Qj (ii) Qj Pj = (Pj+1 − Pj )Pj = Pj+1 Pj − P2j = Pj − Pj = 0 using Proposition 4.1(ii) (iii) Let g = Qj h for some h ∈ L2 [0, 1], and let f ∈ L2 [0, 1] hQj f, gi = hQj f, Qj hi = h(Pj+1 − Pj )f, Qj hi = hPj+1 f, Qj hi − hPj f, Qj hi = hf, Qj hi − hf, Pj Qj hi = hf, Qj hi = hf, gi 30 using Proposition 4.1(iii) (iv) hQj f, gi = h(Pj+1 − Pj )f, gi = hPj+1 f, gi − hPj f, gi = hf, Pj+1 gi − hf, Pj gi using Proposition 4.1(iii) = hf, (Pj+1 − Pj )gi = hf, Qj gi (v) Let f ∈ L2 [0, 1] hQj g, f i = hg, Qj f i = hg, f i by property (iv) by property (iii) The proposition 4.2 shows that the operator Qj shares the same the properties as Pj . Lemma 4.1 The approximation space Sj+1 can be decomposed as follows: Sj+1 = Sj ⊕ Wj Proof. First, we show that Sj+1 = Sj + Wj . 31 (4.5) (⊆ :) Let f ∈ Sj+1 , then f = Pj+1 f = (Pj+1 − Pj )f + Pj f = Qj + Pj f then Sj+1 ⊆ Sj + Wj . (⊇ :) Since Sj ⊂ Sj+1 and Wj ⊂ Sj+1 , then Sj + Wj ⊆ Sj+1 . So, Sj+1 = Sj + Wj . Now, we show that Sj ⊥ Wj hf, gi = hPj f, gi = hPj f, Qj gi = hQj Pj f, gi = 0. Thus, Sj+1 = Sj ⊕ Wj Now, if we keep applying the decomposition relation 4.5 to each approximation subspace Sl , 0 6 l 6 j we get Sj+1 = S0 ⊕ W0 ⊕ W1 ⊕ · · · ⊕ Wj (4.6) The key idea in discrete wavelet transform is that any function f ∈ Sj+1 can 32 be decomposed into two parts in a coarser level, an approximation part in Sj and a detail part in Wj , as in (4.5). Moreover, it can be decomposed into an approximation part in the coarsest level S0 and details in many levels from 0 to j as shown in (4.6). This decomposition process helps in many applications such as denoising, compression, singularity detection, etc. At this stage, the following question arises: Is there a function ψ such that {ψj,k (x) : k ∈ Z} spans the detail subspace Wj ? The following proposition gives the answer. Proposition 4.3 ([28]) Suppose that the scaling function ϕ has a finite refinement mask {ak }k∈Z and generates an orthonormal MRA S. Let bk = (−1)k a1−k then ψ(x) = P k∈Z bk ϕ(2x − k) is a wavelet such that {ψj,k (x) : k ∈ Z} is an or- thonormal basis for Wj . i.e. for every j ∈ Z, the detail operator Qj on a function f ∈ L2 [0, 1] is given by Qj f = hf, Ψj iΨj = d∗j Ψj where d∗j = hf, Ψj i. If we write the coefficients {bk : k ∈ Z} in Definition 4.3 in a matrix b with elements defined in the same way that we did in (4.2), we get the coarse-fine relation Ψj = bΦj+1 . 33 (4.7) 4.1.1 Fast Wavelet Transform in One Dimension Now, we are going to introduce the fast wavelet transform used in the decomposition and the reconstruction of any function f ∈ L2 [0, 1]. This transform uses the refinement equations introduced previously in (4.2) and (4.7) to get the approximation coefficient cj and the detail coefficient dj at the coarser levels from the approximation coefficient cj+1 at the finer level and vice-versa. Decomposition Given a function f ∈ L2 [0, 1], the projection of f into the space Sj can be written as follows: Pj f = Pj−1 f + Qj−1 f hf, Φj iΦj = hf, Φj−1 iΦj−1 + hf, Ψj−1 iΨj−1 . See Lemma 4.1 Let cTj = hf, Φj i and dTj = hf, Ψj i, we have: cTj−1 = hf, Φj−1 i = hf, aΦj i = hf, Φj iaT = cTj aT ∴ cj−1 = acj (∗) 34 dTj−1 = hf, Ψj−1 i = hf, bΦj i = hf, Φj ibT = cTj bT ∴ dj−1 = bcj (∗∗) The equations (∗) and (∗∗) show how to obtain the coarser approximation and detail coefficients from the finer one. This is a decomposition process. Repeating the decomposition process we get the decomposition diagram shown in Figure 4.1. Figure 4.1: The decomposition of an input vector cj . Reconstruction This process is the reverse of decomposition and hence we keep the same definitions of cj and dj . Recall that Pj+1 f = Pj f + Qj f hf, Φj+1 iΦj+1 = hf, Φj iΦj + hf, Ψj iΨj 35 Then cTj+1 Φj+1 = cTj Φj + dTj Ψj =⇒ cTj+1 Φj+1 = cTj aΦj+1 + dTj bΦj+1 =⇒ cTj+1 Φj+1 = (cTj a + dTj b)Φj+1 =⇒ cTj+1 = cTj a + dTj b =⇒ cj+1 = aT cj + bT dj . The last equation shows how to reconstruct cj+1 from cj and dj . We can start from the coarsest level of approximation and details c0 , d0 , d1 , ..., dj and recursively reconstruct the approximation coefficient cj+1 as depicted in Figure 4.2. Figure 4.2: The reconstruction of a vector cj+1 . The reconstruction process is also known as the inverse discrete wavelet transform. Operation Count In this section, we will calculate the number of operations needed to perform the fast wavelet transform and compare it to the fast Fourier transform (FFT). 36 Assume that the input vector cj is of length N . To get the vector cj−1 using the refinement equation (4.2) we need (l2 − l1 + 1) (l2 − l1 + 1) N multiplication operations and 2 N − 1 addition operations. In total, the calculations needed for cj−1 is 2 (l2 − l1 + 1)N − 1. The same number of operations is needed for dj−1 . So, to get the approximation and detail coefficients at level j − 1 we need 2(l2 − l1 + 1)N − 2 operations. Continuing by induction, we get that the total number of operations needed is 2(l2 − l1 + 1)N (1 + 1 1 1 + + . . . + j−1 ) − 2j < 4(l2 − l1 + 1)N 2 4 2 So, the number of operations needed is O(N ). On the other hand, the (FFT) needs O(N logN ) operations and obviously the fast wavelet transform is optimal [28]. 4.2 Discrete Wavelet Transform in Two Dimensions As we mentioned in the introduction of this thesis, images are functions of two variables. In this section, we will extend the notions and results shown in the previous section to the case of two dimensions. We will build scaling function and associated wavelet functions in two dimensions. Also, we will introduce the approximation subspace and the detail subspaces in three directions, horizontal, vertical and diagonal. The two-dimensional space will be thought of as the tensor product of two one-dimensional spaces. Definition 4.3 (The tensor product) Let V and W be two subspaces of 37 L2 [0, 1], with bases E and F respectively. The tensor product of V and W is defined as V ⊗ W = span{f (x)g(y) : f ∈ E, g ∈ F }. (4.8) A typical element in V ⊗ W is of the form n m X X αij fi (x)gj (y). i=1 j=1 Clearly, V ⊗ W ⊆ L2 [0, 1]2 Theorem 4.1 Let U, V, W be subspaces of L2 [0, 1], then (U + V ) ⊗ W = U ⊗ W + V ⊗ W Proof. Let f ∈ U, g ∈ V and h, r ∈ W , then (⊆ :) (f (x) + g(x))h(y) = f (x)h(y) + g(x)h(y) ∴ (U + V ) ⊗ W ⊆ U ⊗ W + V ⊗ W (⊇ :) Let F ∈ U ⊗ W and G ∈ V ⊗ W such that 38 (4.9) F = Pm Pn F = i=1 j=1 m X n X αij fi (x)gj (y) and G = αij fi (x)gj (y) = i=1 j=1 G= r X s X m X n X Pr i=1 Ps j=1 βij hi (x)gj (y) αij (fi (x) + 0)gj (y) ∈ (U + V ) ⊗ W βij (0 + hi (x))gj (y) ∈ (U + V ) ⊗ W i=1 j=1 βij hi (x)gj (y) = i=1 j=1 r X s X i=1 j=1 ∴ F + G ∈ (U + V ) ⊗ W Corollary 4.2 Let U, V, W, Y be subspaces of L2 [0, 1], then (U + V ) ⊗ (W + Y ) = U ⊗ W + U ⊗ Y + V ⊗ W + V ⊗ Y Definition 4.4 Let Am×n and Br×s be two matrices, then    a11 B a12 B . . . a1n B      a B a B ... a B   21 22 2n   A⊗B =  . .. .. ..   .. .  . .       am1 B am2 B . . . amn B Theorem 4.3 If E is an orthonormal basis for V and F is an orthonormal basis for W , then E ⊗ F is an orthonormal basis for V ⊗ W . Proof. Suppose that E = {ξi }i∈Z and F = {ηj }j∈Z , then E ⊗ F = {ξi ηj }i,j∈Z . First, we show that any function in V ⊗ W can be represented as a linear combination of functions in E ⊗ F . 39 Let G ∈ V ⊗ W , then G= k X l X αij fi gj , i=1 j=1 such that fi = X βr ξr and gj = r∈Ii so fi g j = X s∈Ij X βr βs′ ξr ηs r∈Ii ,s∈Ij = X βrs ξr ηs r∈Ii ,s∈Ij so G= l k X X X αij βrs ξr ηs i=1 j=1 r∈Ii ,s∈Ij = X ′ βrs ξ r ηs r∈Ii ,s∈Ij G = β ′ (E ⊗ F ). hence Now, we show the orthonormality of the elements of E ⊗ F hξi ηj , ξk ηl i = = Z∞ Z∞ ξi (x)ηj (y)ξk (x)ηl (y)dxdy −∞ −∞ Z∞ ξi (x)ξk (x)dx Z∞ ηj (y)ηl (y)dy −∞ −∞ = hξi ξk ihηj , ηl i = δi,k δj,l = 1 ∴ when the elements of E ⊗ F 40 i = k, j = l are orthonormal. βs′ ηs 4.2.1 MRA in L2 [0, 1]2 Let ϕ generate an orthonormal MRA S = {Sj }j∈Z in L2 [0, 1]. Now, we want to build an MRA in L2 [0, 1]2 . Let Sj2 := Sj ⊗ Sj = (Sj−1 ⊕ Wj−1 ) ⊗ (Sj−1 ⊕ Wj−1 ) = (Sj−1 ⊗ Sj−1 ) ⊕ (Sj−1 ⊗ Wj−1 ) ⊕ (Wj−1 ⊗ Sj−1 ) ⊕ (Wj−1 ⊗ Wj−1 ) | {z } 2 Wj−1 2 2 =: Sj−1 ⊕ Wj−1 , 2 where Sj−1 is the approximation space at the coarser level j − 1, Sj−1 ⊗ Wj−1 is the subspace of vertical details, Wj−1 ⊗ Sj−1 is the subspace of horizontal details and Wj−1 ⊗ Wj−1 is the subspace of diagonal details. Orthonormal bases of these spaces are as follows: (See Theorem 4.3) For Sj2 , a basis is {ϕj,k (x)ϕj,k (y) : k, l ∈ Z} = Φj ⊗ Φj = Φ2j , and for Sj−1 ⊗ Wj−1 , a basis is {ϕj,k (x)ψj,k (y) : k, l ∈ Z} = Φj ⊗ Ψj = Ψ21 j , and for Wj−1 ⊗ Sj−1 , a basis is {ψj,k (x)ϕj,k (y) : k, l ∈ Z} = Ψj ⊗ Φj = Ψ22 j , and for Wj−1 ⊗ Wj−1 , a basis is {ψj,k (x)ψj,k (y) : k, l ∈ Z} = Ψj ⊗ Ψj = Ψ23 j .   Ψ21   j    22 . Let Ψ2j =  Ψj      Ψ23 j 41 The projections from L2 [0, 1]2 to the spaces Sj2 and Wj2 are defined as follows: P2j : L2 [0, 1]2 −→ Sj2 P2j = hf, Φ2j iΦ2j Q2j : L2 [0, 1]2 −→ Wj2 Q2j = hf, Ψ2j iΨ2j    Ψ21 j   Ψ21 j             22  = hf, Ψ22 i j  Ψj          23 Ψ Ψ23 j j   Ψ21   j    22 23  22  = (hf, Ψ21 j i hf, Ψj i hf, Ψj i) Ψj      23 Ψj 21 22 22 23 23 = hf, Ψ21 j iΨj + hf, Ψj iΨj + hf, Ψj iΨj . Let cTj = hf, Φ2j i,   21 d1T j = hf, Ψj i, 22 d2T j = hf, Ψj i, d1j       and dj =  d2j .     d3j 42 23 d3T j = hf, Ψj i Refinement Equations Recall that Φj = aΦj+1 and Ψj = bΦj+1 , where a and b are given by ak,m = 1 1 √ am−2k and bk,m = √ bm−2k , hence 2 2 Φ2j = Φj ⊗ Φj = aΦj+1 ⊗ aΦj+1 = (a ⊗ a)(Φj+1 ⊗ Φj+1 )    = (a ⊗ a)Φ2j+1 =: a2 Φ2j+1 | {z } 2  a  (∗) Ψ21   Φj ⊗ Ψj   (a ⊗ b)Φ2j+1   j            2  =  (b ⊗ a)Φ2  22  =  Ψj =  Ψ Ψ ⊗ Φ  j   j  j j+1              23 2 Ψj Ψj ⊗ Ψ j (b ⊗ b)Φj+1   a ⊗ b     2  =  b ⊗ a  Φj+1     b⊗b | {z } b2 =: b2 Φ2j+1 . (∗∗) Equations (∗) and (∗∗) will play an important rule in the decomposition and reconstruction processes. 4.2.2 Fast Wavelet Transform in Two Dimensions After we showed the fast wavelet transform in one dimension, we will now show the fast wavelet transform in the case of two dimensions. The derivation is the 43 same. Decomposition For any function f ∈ L2 [0, 1]2 , we get the approximation coefficient in a coarser level from the finer one as follows: cTj = hf, Φ2j i = hf, a2 Φ2j+1 i = hf, Φ2j+1 ia2T = cTj+1 a2T ∴ cj = a2 cj+1 and we get the detail coefficient as follows: dTj = hf, Ψ2j i = hf, b2 Φ2j+1 i = hf, Φ2j+1 ib2T = cTj+1 b2T ∴ dj = b2 cj+1 This is a decomposition process. Repeating the decomposition process we get the decomposition diagram shown in Figure 4.1. 44 Reconstruction The following equation shows how to reconstruct cj+1 from cj and dj . We can start from the coarsest level of approximation and details c0 , d0 , d1 , ..., dj and recursively reconstruct the approximation coefficient cj+1 as depicted in Figure 4.2. cTj+1 Φ2j+1 = cTj Φ2j + dTj Ψ2j = cTj a2 Φ2j+1 + dTj b2 Φ2j+1 = (cTj a2 + dTj b2 )Φ2j+1 ∴ cj+1 = a2 T cj + b2 T dj 45 CHAPTER 5 EDGE DETECTION Discrete wavelet transform has many applications, including but not limited to compression, denoising and singularity detection. In this chapter, we are interested in image edge detection which can be defined as the process of locating and identifying the sharp discontinuities in an image [17]. First, we are going to consider singularity detection in one dimension. Then, we will extend the work to the case of discontinuities in two dimensions, i.e. discontinuities in images. After that, we will give a brief overview of Canny edge detector and compare the results of our algorithm to Canny’s. In both detectors, we will add noise and observe the effect of noise on edge detection. Furthermore, we will experiment with our algorithm using Haar wavelet, Daubechies wavelets db2 and db3 to illustrate the distinction of the Haar wavelet in edge detection. Then, we will apply our algorithm to colored images and compare the results to Canny’s. 46 5.1 Detecting Singularities of Functions in One Dimension Detecting singularities of functions in one dimension can be done using the discrete wavelet transform by looking for high magnitude detail coefficients. In this work, we are going to use the Haar wavelet because it matches well abrupt changes in signals. Theorem 5.1 will show that the detail coefficients are small when a function f is smooth. As we shall see, it provides the justification for using wavelets in edge detection. Definition 5.1 (Vanishing moment) The function f : Ω −→ R is said to have d vanishing moments if Z xr f (x)dx = 0, r = 0, 1, . . . , d − 1 Ω (5.1) Theorem 5.1 ([25]) For n ∈ N, let the function g be C n (R) and g (n) (x) be L∞ (R). Suppose that the function ψ(x) has compact support and n vanishing moments and that R R |ψj,k (x)|2 dx = 1 for all j, k ∈ Z. Then there exists a constant C > 0 that depends only on n and g such that |hg, ψj,k i| 6 C2−jn 2−j/2 . Proof. First, suppose that ψ(x) is supported in the interval I = [0, a] for some a > 0. Recall that ψj,k (x) is supported in the interval Ij,k = [2−j k, 2−j (k + a)] 47 and |Ij,k | is 2−j a. Denote the center of the interval Ij,k by xj,k and note that xj,k = 2−(j+1) a + 2−j k. Since ψj,k (x) has n vanishing moments, given any polynomial p(x) of degree no greater than n − 1 we have Z p(x)ψj,k (x)dx = 0. (5.2) R Since g(x) is C n (R), for every j, k ∈ Z, g(x) can be expanded about the point xj,k in a Taylor expansion as follows: g(x) = g(xj,k )+(xj,k −xj,k )g ′ (xj,k )+. . .+ 1 (xj,k −xj,k )n−1 g (n−1) (xj,k )+rn (x), (n − 1)! (5.3) where rn (x) = 1 (xj,k − xj,k )n g (n) (ξ) n! for some number ξ between xj,k and xj,k . If x ∈ Ij,k , then xj,k − xj,k 6 2−(j+1) a, which implies that |rn (x)| 6 1 −n(j+1) n 2 a max |g (n) (x)|. x∈Ij,k n! 48 (5.4) Now, we apply (5.2) to (5.3) and compute hg, ψj,k i = Z g(x)ψj,k (x)dx R ! n−1 X 1 (xj,k − xj,k )l g (l) (xj,k ) + rn (x) ψj,k (x)dx = l! R l=0 ! Z n−1 X1 g (l) (xj,k ) (xj,k − xj,k )l ψj,k (x)dx = l! R l=0 Z + rn (x)ψj,k (x)dx Z R = Z rn (x)ψj,k (x)dx, by (5.2). Ij,k Then we apply (5.4) and the Cauchy-Schwarz inequality to get |hg, ψj,k i| = Z rn (x)ψj,k (x)dx Ij,k 1 6 2−n(j+1) an max |g (n) (x)| x∈Ij,k n! Z Ij,k |ψj,k (x)|dx 1 6 2−n(j+1) an max |g (n) (x)||Ij,k |1/2 x∈Ij,k n! Z 2 Ij,k |ψj,k (x)| dx !1/2 1 max |g (n) (x)|2−n(j+1) an 2−j/2 a1/2 n! x∈Ij,k   1 −n n+1/2 −nj −j/2 (n) =2 2 2 a max |g (x)| x∈Ij,k n!   1 −n n+1/2 (n) −nj −j/2 2 a ||g ||∞ . =2 2 n! = with C = 1 −n n+1/2 (n) 2 a ||g ||∞ , Theorem 5.1 is satisfied. n! The above theorem suggests that the magnitude of the detail coefficients that correspond to the smooth parts of a function will be very small compared with 49 the ones corresponding to the non-smooth parts of a function. This result has a crucial implication to edge detection in images, since edges can be interpreted as non-smooth parts of the image and hence can be detected by looking for high detail coefficients. The following example will substantiate this statement. Example 5.2 Consider the function f defined as follows: f (x) =     1    0 if 0 6 x < 1/3 if 1/3 6 x < 1. It is clear that f is continuous and continuously differentiable infinitely many 1 1 1 times on (0, ) and ( , 1) and it is discontinuous at x = . Now, we will calculate 3 3 3 the detail coefficients at level j = 1. Then, we will give a general formula for calculating dj,k for all values of j and k. This will show that, in the subinterval that has the discontinuity, the detail coefficient is nonzero while, in the subintervals that have smooth parts the detail coefficients are zero. H (x) = 2j/2 ψ(2j x − k) Recall first that ψj,k d1,0 d1,1 √ √ √ 2 1 1 H = hf, ψ1,0 i= 2dx + (− 2)dx + 0dx = ( ) 2 − ( ) 2 + 0 = 4 12 6 0 1/4 1/3 3/4 √ √ R R1 H = hf, ψ1,1 i = (0) 2dx + (0)(− 2)dx = 0 + 0 = 0 1/4 R √ 1/2 1/3 R √ 1/2 R 3/4 Now, we are going to deduce a general formula for calculating dj,k . First, let us 50 determine the dyadic subinterval Ij,k′ = [ 1 k′ k′ + 1 ] ∈ [ j, 3 2 2j k′ 1 > j 3 2 ′ 1 k +1 6 3 2j ∴ 1 k′ k′ + 1 , ] that contains the value x = . 2j 2j 3 1 k′ and > j 3 2 2j ′ k 6 3 2j ′ k > −1 3 ⇒ ⇒ ⇒ 1 k′ + 1 6 3 2j 2j k =⌊ ⌋ 3 ′ So, for any value of j we have three cases for k. We will show the value of the detail coefficients in each case. Case I: k < ⌊ 2j ⌋ 3 In this case, x < 1 which implies that f (x) = 1. So 3 dj,k = H hf, ψj,k i = Z H f ψj,k dx = Ij,k Z H ψj,k dx = 0 Ij,k 2j ⌋ 3 1 In this case, x > which implies that f (x) = 0. So 3 Case II: k > ⌊ dj,k = H hf, ψj,k i = Z H f ψj,k dx = Ij,k Case III: k = k ′ = ⌊ In this case, Z Ij,k 2j ⌋ 3 1 ∈ Ij,k′ so we have two subcases: 3 51 H (0)ψj,k dx = 0 Case III(a): 1 l ∈ Ij,k ′ , which gives 3 H dj,k′ = hf, ψj,k ′i = Z Ij,k′ Z1/3 H H ψj,k f ψj,k ′ dx ′ dx = k′ 2j Z1/3 1 k′ 2j/2 dx = ( − j )(2j/2 ) > 0 = 3 2 k′ 2j k′ 1 > j 3 2 since Case III(b): 1 r ∈ Ij,k ′ , which gives 3 H dj,k′ = hf, ψj,k ′i = Z H f ψj,k ′ dx = Ij,k′ =0+ Z l Ij,k ′ Z1/3 H ψj,k ′ dx + Z1/3 H f ψj,k ′ 2k′ +1 2j+1 1 2k ′ + 1 2j/2 = ( − j+1 )(2j/2 ) > 0 3 2 2k′ +1 2j+1 since 1 2k ′ + 1 > j+1 3 2 Figure 5.1 shows the function f (x) and the detail coefficients for the levels j = 1, 2, 3. You can see clearly that the coefficient of the subinterval that contains the discontinuity is high while all other coefficients are zeros. 52 Figure 5.1: (a) The function f (x), (b) The detail coefficients for j = 1, (c) The detail coefficients for j = 2 and (d) The detail coefficients for j = 3 The proposed algorithm in one dimension Now, we present the most important part in this section which is an algorithm to detect singularities using discrete wavelet transform. Suppose that we have a function f ∈ L2 [0, 1], and we want to detect singularities in it. The algorithm proposed is as follows: Algorithm 5.3 (Wavelet-based singularity detection algorithm) 53 1) Apply discrete wavelet transform in one dimension at scale j, chosen a priori, which produces the approximation and detail coefficient vectors cj and dj . 2) Discard the approximation coefficient vector cj . 3) Identify the detail coefficient dj,k that has the maximum modulus in dj and call it m. 4) Choose a threshold T . Here we take T = m . 10 5) Discard all the detail coefficients that are less than T in modulus. 6) The remaining coefficients specify the dyadic subintervals where the singularities are located. The higher the scale j, the more precise the singularity is bracketted since Ij,k = [2−j k, 2−j k + 1]. The reason for step 2 is that the approximation coefficient cj contains the smooth part of the signal, since Sj contains all polynomials Pd of degree d − 1, where d is the vanishing moment of ψ. Singularities are manifested only in the detail coefficients since Wj is orthogonal to these polynomials. Figure 5.2 shows the flow chart of this algorithm. 54 Figure 5.2: Flow chart for wavelet-based singularity detection algorithm 55 5.2 Edge Detection Using Discrete Wavelet Transform in Two Dimensions It is worth noting that images are functions of two variables [29]. Moreover, they are bounded and at the same time have bounded variations [9]. Discrete wavelet transform in two dimensions has many applications in image processing such as compression and denoising. Most importantly for us, it enables us to detect edges by tracking the high detail coefficients in vertical, horizontal and diagonal directions. In this section, we will present an algorithm for image edge detection. In fact, our algorithm is analogous to the one we proposed in the previous section. Edge detection in the presence of noise will also be tried. Images considered in this work, are all high-definition square images of size 2048 × 2048. Indeed, our algorithm outperformed the Canny edge detector, which is the most popular one used in the literature. The highest scale j that can be reached is log2 ( p size of image) − 1 (∗) . The proposed algorithm in two dimensions Suppose that we have an image u ∈ L2 [0, 1]2 and we want to detect the edges in u. The algorithm proposed is as follows: Algorithm 5.4 (Wavelet-based edge detection algorithm) 56 1) Convert the image to a black-and-white image. 2) Apply discrete wavelet transform in two dimensions at scale j according to (∗) to produce the approximation and detail coefficient vectors cj and dj . T  Recall that dj = d1j d2j d3j . (See 4.2.1) 3) Discard the approximation coefficient vector cj . 4) Identify the detail coefficient dlj,k that has the maximum modulus for every l with 1 6 l 6 3 and call them m1 , m2 and m3 , respectively. Then choose the maximum one among them and call it m. 5) Let the threshold be T = m . 10 6) Discard all the detail coefficients in the three directions that are less than T in modulus. 7) Plot the location of detail coefficients after thresholding to get the image edges. It is worth noting that this algorithm is applied using MATLAB software. Figure 5.3 shows the flow chart of this algorithm. 57 Figure 5.3: Flow chart for wavelet-based edge detection algorithm 58 5.3 Edge Detection Using Canny edge detector John Canny proposed his detector in 1986 [15]. In this section, we will give a short description of Canny edge detector. His detector is dominating nowadays due to its great advantages [1]. Given an image u, the algorithm that Canny proposed is as follows: Algorithm 5.5 (Canny edge detector algorithm) 1) Convert the image to a black-and-white image. 2) Convolve the image u with a Gaussian function G. This step is done to smooth the image. i.e. w = u ∗ G. 3) Take the modulus of gradient of the convolved image. i.e. Find |∇w| = |∇(u ∗ G)| 4) Threshold by looking for pixels of high modulus gradient |∇w| > T . 5) Suppress any pixel that is not maximum in its neighborhood. Indeed, here it only compares the pixel with the two neighboring pixels that are in the direction of the tangent θ = arctan wy . This step is called non-maximum wx suppression to thin the edges. 6) Double threshold by setting two thresholds Thigh and Tlow and mark the pixels above Thigh as strong, pixels between Thigh and Tlow as weak and suppress pixels below Tlow . 59 7) Track the true edges by comparing the weak pixels to the 8 neighboring pixels. If a weak pixel is connected to a strong one, then it is preserved. Otherwise, it is considered false edge and, therefore, suppressed. Note that, the Canny edge detector is applied to images using the MATLAB command ”edge(I,’canny’)”. 60 Figure 5.4: Flow chart for Canny edge detector algorithm Figure 5.4 shows the flow chart of this algorithm. 61 5.4 Results and Comparison of Images Without Noise In this section, we compare edge detection by wavelets to that with Canny. We begin by five images without noise. Then we add noise and observe its effect. The first image shows a Pepsi can with a blurred background (Figure 5.5). It is clear that our algorithm outperforms Canny’s. Using our algorithm, we only miss some circles with weak boundaries in the background. For Canny algorithm, the true edges are not clear and it produces many false edge pixels. Figure 5.5: Pepsi can image 62 Figure 5.6: Pepsi can edges using dwt2 with Haar Figure 5.6 shows the Pepsi can edges using Haar wavelet. 63 Figure 5.7: Pepsi can edges using Canny edge detector Figure 5.7 shows the Pepsi can edges using Canny edge detector. 64 The second image that we are going to take is an image of the famous Times Square in New York (Figure 5.8). Looking at the two edge images it is easily seen that our algorithm outperforms Canny’s. Using our algorithm, the edges are much clearer and the constituent parts in the image are easy to distinguish. For example, you can easily see the edges of the cars and also read the word ”Kodak” written vertically which is not the case with Canny’s. Figure 5.8: Times square image 65 Figure 5.9: Times square edge image using dwt2 with Haar Figure 5.9 shows the edge image of the Times Square using Haar wavelet. 66 Figure 5.10: Times square edge image using Canny edge detector Figure 5.10 shows the edge image of the Times Square using Canny edge detector. 67 The third image is of a bird with a blurred background (Figure 5.11). Figure 5.11: Bird image 68 Figure 5.12: Bird edge image using dwt2 with Haar Figure 5.12 shows the edge image of the bird image using Haar wavelet. 69 Figure 5.13: Bird edge image using Canny edge detector Figure 5.13 shows the edge image of the bird image using Canny edge detector. 70 The fourth image is of a Dr Pepper can. (Figure 5.14). Figure 5.14: Dr Pepper can image 71 Figure 5.15: Dr Pepper can edge image using dwt2 with Haar Figure 5.15 shows the edges of the Dr Pepper can image using Haar wavelet. 72 Figure 5.16: DrPepper can edge image using Canny edge detector Figure 5.16 shows the edges of the Dr Pepper can image using Canny edge detector. 73 The fifth image is an image of a toys room (Figure 5.17). Figure 5.17: Toys room image 74 Figure 5.18: Toys room edge image using dwt2 with Haar Figure 5.18 shows the edges of the toys room image using Haar wavelet. 75 Figure 5.19: Toys room edge image using Canny edge detector Figure 5.19 shows the edges of the toys room image using Canny edge detector. 5.5 Results and Comparison of Images with Noise Now we will see the effect of adding the salt and pepper noise to the Pepsi can (Figure 5.20) and Times Square (Figure 5.23) images. Here we use the same 76 algorithm proposed but we lower the threshold to Tnoise = m m , instead of T = . 20 10 Figure 5.20: Noised Pepsi can image 77 Figure 5.21: Noised Pepsi can edge image using dwt2 with Haar Figure 5.21 shows the edges of the Pepsi can noised image using Haar wavelet. 78 Figure 5.22: Noised Pepsi can edge image using Canny edge detector Figure 5.22 shows the edges of the Pepsi can noised image using Canny edge detector. 79 Figure 5.23: Noised Times square image 80 Figure 5.24: Noised Times square edge image using dwt2 with Haar Figure 5.24 shows the edges of the Times Square noised image using Haar wavelet. 81 Figure 5.25: Noised Times square edge image using Canny edge detector Figure 5.25 shows the edges of the Times Square noised image using Canny edge detector. It is clear that our algorithm is not immune to noise but still it is better than Canny edge detector with respect to noise. 5.6 Comparison of Haar with db2 and db3 in edge detection In this section, we will apply our algorithm using Haar, db2 and db3 wavelets to the Pepsi can image and compare the edge images in all cases. This section, will 82 show the distinction of the Haar wavelet in edge detection as we proposed in the beginning of this work. This distinction is due to the better correlation of the Haar wavelet with the characteristic edge jump discontinuity. Figure 5.26 shows the Pepsi can image. Figure 5.26: Pepsi can image 83 Figure 5.27: Pepsi can edge image using dwt2 with Haar Figure 5.27 shows the edges of the Pepsi can image using Haar wavelet. 84 Figure 5.28: Pepsi can edge image using dwt2 with db2 Figure 5.28 shows the edges of the Pepsi can image using db2 wavelet. 85 Figure 5.29: Pepsi can edge image using dwt2 with db3 Figure 5.29 shows the edges of the Pepsi can image using db3 wavelet. It clear that the Haar wavelet is distinct in detecting edges. Moreover, the edge image loses sharpness of edges as we go from db1(Haar) to db2 to db3 an so on. 5.7 Edge Detection In Colored Images In this section, we will directly apply our algorithm to colored images instead of transforming first to a black-white image. The algorithm will be applied to 86 the colored Pepsi can image (Figure 5.30) and to the toys room image (Figure 5.33). Figures 5.31 and 5.32 show the edges of the colored Pepsi can image using Haar wavelet and Canny edge detector, respectively. Figures 5.34 and 5.35 show the edges of the toys room image using Haar wavelet and Canny edge detector, respectively. Recall that a colored image can be thought of as a function u : [0, 1] × [0, 1] −→ R3 where   r       u(x, y) =  g      b and r, g, b are the red, green and blue color intensities at the point (x, y). The idea is to apply the edge detection algorithm to each of the rgb channels representing the color intensities at the image pixels. 87 Figure 5.30: Colored Pepsi can image Figure 5.31: Colored Pepsi can edge image using dwt2 with Haar 88 Figure 5.32: Colored Pepsi can edge image using Canny edge detector Figure 5.33: Colored toys room image 89 Figure 5.34: Colored toys room edge image using dwt2 with Haar Figure 5.35: Colored toys room edge image using Canny edge detector Clearly, from the above images you can see that our algorithm outperforms 90 the Canny edge detector when colored images are processed instead of black-white images. Concluding Remarks 5.6 At the end of this work, we can summarize our findings as follows: ❼ Our algorithm outperforms Canny edge detector when processing high- definition images. ❼ Our algorithm is much more immunized than Canny edge detector when processing noised high-definition images. ❼ Our algorithm is performing the best when using the Haar wavelet. ❼ Our algorithm outperforms Canny edge detector when colored images are processed. Conclusion 5.7 fsd 91 REFERENCES [1] Rashmi, M. Kumar, and S. Rohini, “Algorithm and Technique on Various Edge Detection : A Survey,” Signal & Image Processing : An International Journal, vol. 4, no. 3, pp. 65–75, 2013. [2] S. Bhardwaj and A. Mittal, “A Survey on Various Edge Detector Techniques,” Procedia Technology, vol. 4, pp. 220–226, 2012. [3] M. Ali and D. Clausi, “Using the Canny edge detector for feature extraction and enhancement of remote sensing images,” in IEEE International Geoscience and Remote Sensing Symposium, vol. 5. IEEE, 2001, pp. 2298–2300. [4] Sivakumar P and Meenakshi S, “A REVIEW ON IMAGE SEGMENTATION TECHNIQUES,” International Journal of Advanced Research in Computer Engineering & Technology, vol. 5, no. 3, 2016. [5] P. A. Hajare and P. A. Tijare, “Edge Detection Techniques for Image Segmentation,” International Journal Of Computer Science And Applications, vol. 4, no. 1. 92 [6] R. Maini and H. Aggarwal, “Study and Comparison of Various Image Edge Detection Techniques,” Himanshu Aggarwal International Journal of Image Processing, vol. 3, no. 1, 2009. [7] M. Gudmundsson, E. El-Kwae, and M. Kabuka, “Edge detection in medical images using a genetic algorithm,” IEEE Transactions on Medical Imaging, vol. 17, no. 3, pp. 469–474, jun 1998. [8] Z. Yu-qian, G. Wei-hua, C. Chen Zhen-cheng, T. Jing-tian, and L. Ling-yun, “Medical Images Edge Detection Based on Mathematical Morphology,” in IEEE Engineering in Medicine and Biology 27th Annual Conference, vol. 6. IEEE, 2005, pp. 6492–6495. [9] S. Mallat, A Wavelet Tour of Signal Processing, 2nd ed. Academic Press, 1999. [10] L. S. Davis, “A Survey of Edge Detection Techniques,” Computer Graphics and Image Processing, vol. 4, no. 3, pp. 248–270, 1975. [11] L. Roberts, “Machine Perception of Three-Dimensional Solids,” Ph.D. dissertation, MIT, 1963. [12] I. Sobel and G. Feldman, “A 3x3 isotropic gradient operator for image processing,” a talk at the Stanford Artificial Intelligence Project, 1968. [13] J. Prewitt, “Object Enhancement and Extraction,” in Picture Processing and Psychopictorics, B. Lipkin and A. Rosenfeld, Eds. ch. 1, pp. 75–151. 93 Academic Press, 1970, [14] D. Marr and E. Hildreth, “Theory of edge detection.” Proceedings of the Royal Society of London. Series B, Biological sciences, vol. 207, no. 1167, pp. 187–217, feb 1980. [15] J. Canny, “A Computational Approach to Edge Detection,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-8, no. 6, pp. 679–698, nov 1986. [16] F. Bergholm, “Edge Focusing,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. PAMI-9, no. 6, pp. 726–741, nov 1987. [17] V. Chaganti, “Edge Detection of Noisy Images Using 2-D Discrete Wavelet Transform,” Master’s thesis, Florida State University, 2005. [18] D. Mumford and J. Shah, “Optimal approximations by piecewise smooth functions and associated variational problems,” Communications on Pure and Applied Mathematics, vol. 42, no. 5, pp. 577–685, jul 1989. [19] J.-M. Morel and S. Solimini, “Segmentation of images by variational methods: a constructive approach,” 1988. [20] S. Solimini and J. M. Morel, “Variational methods in image segmentation,” Progress in Nonlinear Differential Equations and Their Applications, 1995. [21] U. Massari and I. Tamanini, “Regularity properties of optimal segmentations,” Journal fur die reine und angewandte Mathematik, vol. 420, pp. 61–84, 1991. 94 [22] S. Mallat and S. Zhong, “Characterization of Signals from Multiscale Edges,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 14, no. 7, pp. 710–732, jul 1992. [23] R. Monika, “Image Edge Detection using Discrete Wavelet Transform,” International Journal of Innovative Research in Computer and Communication Engineering, vol. 4, no. 4, pp. 169–175, 2016. [24] H. L. Royden and P. M. Fitzpatrick, Real Analysis, 4th ed. Pearson, 2010. [25] D. F. Walnut, An introduction to wavelet analysis, 2nd ed. Birkhauser, 2004. [26] A. Haar, “Zur theorie der orthogonalen funktionensysteme,” Mathematische Annalen, vol. 69, no. 3, pp. 331–371, 1910. [27] R. Islam, “Approximation By Wavelets,” Master’s thesis, KFUPM, 2005. [28] K. Urban, Wavelet methods for elliptic partial differential equations. Oxford University Press, 2009. [29] A. Chambolle, R. A. De Vore, N.-Y. Lee, and B. J. Lucier, “Nonlinear wavelet image processing: variational problems, compression, and noise removal through wavelet shrinkage,” IEEE Transactions on Image Processing, vol. 7, no. 3, pp. 319–335, 1998. 95 Vitae ❼ Name: Azzam Abdulaziz Alfarraj ❼ Nationality: Saudi ❼ Date of Birth: 21, Feb 1992 ❼ Email: azzamfr@gmail.com ❼ Permenant Address: P.O. Box 450, KFUPM, Dhahran 31261, Saudi Arabia 96