Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00003971/00001
## Material Information- Title:
- Image denoising with wavelets using Markov random fields
- Creator:
- Ravindra, Vishal C
- Publication Date:
- 2005
- Language:
- English
- Physical Description:
- xi, 84 leaves : ; 28 cm
## Subjects- Subjects / Keywords:
- Wavelets (Mathematics) ( lcsh )
Electronic noise ( lcsh ) Image processing -- Digital techniques ( lcsh ) Electronic noise ( fast ) Image processing -- Digital techniques ( fast ) Wavelets (Mathematics) ( fast ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Bibliography:
- Includes bibliographical references (leaf 84).
- General Note:
- Department of Electrical Engineering
- Statement of Responsibility:
- by Vishal C. Ravindra.
## 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:
- 63755067 ( OCLC )
ocm63755067 - Classification:
- LD1193.E54 2005m R38 ( lcc )
## Auraria Membership |

Full Text |

IMAGE DENOISING WITH WAVELETS USING MARKOV RANDOM
FIELDS by Vishal C. Ravindra B.E., University of Mysore, India, 2001 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 This thesis for the Master of Science degree by Vishal C.Ravindra has been approved Jan Bialaseiwicz Date Ravindra, Vishal C. (M.S., Electrical Engineering) Image Denoising with Wavelets using Markov Random Fields Thesis directed by Professor Jan Bialaseiwicz ABSTRACT This thesis describes a method for the suppression of noise in images via the wavelet transform. The method relies on two measures. The first is a classic measure of the smoothness of the image and is based on an approximation of the local Lipschitz exponent obtained via the wavelet coefficients. The second, relatively new measure takes into account a priori geometrical constraints of the image, which are generally valid for natural images. The smoothness measure and the geometrical constraints are combined in a Bayesian probabilistic framework, and are implemented as a Markov Random Field image model. The manipulation of the wavelet coefficients is consequently based on the obtained probabilities. This abstract accurately represents the content of the candidates thesis. I recommend its publication. Signed m DEDICATION I dedicate this thesis to my parents, my sister, and other family members without whose understanding and support, I would not be in Graduate School. ACKNOWLEDGEMENT My thanks to my advisor, Dr. Jan Bialaseiwicz, for introducing me to the world of Wavelets, and for his guidance and patience with me for this past year. I also wish to thank the staff and faculty of the Graduate School for their support and understanding. CONTENTS Figures............................................................ ix Tables............................................................. xi Chapter 1. Introduction..................................................... 1 2. Brief Review on Wavelets......................................... 3 2.1 Transforms...................................................... 3 2.2 Basis Functions................................................ 3 2.3 Fourier Transforms and the Time-Frequency Wedding.............. 4 2.3.1 Windowed Fourier Transform.................................... 6 2.4 Wavelet Transforms............................................. 8 2.4.1 Wavelet Admissibility Condition............................... 9 2.4.2 Scaling Function.............................................. 9 2.5 Frames and Bases.............................................. 10 2.5.1 Pseudo Inverse............................................. 11 2.5.2 Partial Reconstruction....................................... 12 2.6 Wavelet Frames................................................ 12 2.7 Wavelet Bases................................................. 13 2.8 Continuous Wavelet Transform.................................. 14 vi 2.9 Dyadic Wavelet Transform....................................... 14 2.10 Wavelet Design................................................ 15 2.10.1 Reconstruction Wavelets...................................... 15 2.11 Spline Dyadic Wavelets........................................ 16 3. Algorithme a Trous............................................. 20 4. Singularity Detection and Lipschitz Regularity................. 28 4.1 Wavelet Vanishing Moments.................................... 30 5. Wavelet Decomposition and Notations.............................. 33 5.1 Orthogonal Wavelet Bases and Multiresolution Approximations... 35 6. The Proposed Method.............................................. 38 6.1 Notations....................................................... 39 6.2 Bayesian Probability........................................... 40 7. Markov Random Fields............................................. 43 7.1 Introduction.................................................... 43 7.2 Markov Random Fields and Gibbs Random Fields.................... 44 7.2.1 Gibbs Random Fields........................................... 47 7.2.2 MRF-GRF Equivalence........................................... 48 7.3 Regularity Estimation........................................... 49 7.4 Lipschitz Threshold and Initial Mask Selection.................. 50 8. Conditional Probability Model.................................... 53 Vll 8.1 Stochastic Sampling.............................................. 53 9. Monte Carlo Methods for Stochastic Sampling....................... 55 9.1 Metropolis-Rosenbluth-Rosenbluth-Teller-Teller Algorithm......... 55 9.2 Image Reconstruction............................................ 63 10. Conclusion....................................................... 66 Appendix A. Algorithme a Trous MATLAB Program................................ 67 B. Metropolis Sampler MATLAB Program................................ 77 C. Variance for Lipschitz Estimator Thresholding MATLAB Program........................................................... 79 D. Circular Convolution, Leftshift, and Zero-padding MATLAB Programs.......................................................... 81 E. Notations Used................................................... 83 References........................................................... 84 vm FIGURES Figure 2.1 Fleisenberg box for time-frequency atoms....................... 5 2.2 Biorthogonal spline wavelet................................... 18 2.3 Biorthogonal cubic B-spline scaling function.................. 18 2.4 Decomposition and reconstruction filters for quadratic spline. 19 3.1 The original house 256X256 bitmap image...................... 23 3.2 The noisy house bitmap image (SNR=5dB)...................... 24 3.3 Approximations, horizontal and vertical details for wavelet decomposition for levels 1 and 2............................... 25 3.4 Decomposition filter bank for the Algorithme a Trous.......... 26 3.5 Reconstruction filter bank for the Algorithme a Trous......... 27 7.1 (a) First order neighbourhood (b) second order (c) fifth order, X is the center point (d) single state (e) two state (f) two state, and (g) three states neighbourhoods.............................. 46 7.2 Initial masks for horizontal and vertical components obtained by thresholding Lipschitz estimator with a threshold T= 1.25...... 52 9.1 Metropolis sampler............................................... 58 9.2 Horizontal mask after 1 iteration of Metropolis algorithm....... 59 9.3 Vertical mask after 1 iteration of Metropolis algorithm......... 60 9.4 Modified horizontal details after shrinkage factors have been applied to each coefficient for levels 1 and 2............ 61 IX 9.5 Modified vertical details after shrinkage factors have been applied to each coefficient for levels 1 and 2............... 62 9.6 Approximations, horizontal and vertical details after wavelet reconstruction for levels 1 and 2................................. 64 9.7 Input noisy image, reconstructed image, and the original clean house image...................................... 65 x TABLES Tables E. 1. Notations Used............................................... 83 xi 1. Introduction Although image acquisition techniques yield an ever improving quality, there is a need for post-processing methods to remove noise from images. The problem of suppressing noise in digital images is based on the model y~f+n where, / denotes the true noise-firee pixels, y denotes the observed noisy pixels and n the observed noise. The goal is to find an optimal approximation to/, based on y and possibly knowledge about the noise. The method used in this thesis, is not specific for a particular type of noise, however, it will deal with white noise for simplicity. The ability to suppress noise is valuable in many applications, such as suppressing of motion blur in aerial imaging, and the noise in the visible channel due to fog and clouds. Denoising (noise suppression) is done at the time of acquisition or by post-processing methods. A number of methods for suppression of noise in signals and images based on the wavelet transform have been proposed and worked on in the past few years. The common idea is to compute the wavelet decomposition of the noisy image, and then to modify the wavelet coefficients in a specific way, so that a relatively clean image is obtained after reconstruction. The simplest approach is to replace the coefficients that are supposed to be affected by noise with zero, or by any other suitable value. In the wavelet shrinkage technique (Donoho and 1 Johnstone), the wavelet coefficients that have an absolute value below a threshold are replaced by zero. The other coefficients are reduced in absolute value. In the approach of (Xu, Weaver, et at), the criterion to distinguish noise from meaningful signal is based on the observation that wavelet coefficients of noise have a much weaker correlation between scales than the coefficients of a noiseless image. The method of (Mallat, et at) is based on the assumption that the noiseless image is regular and the noisy image is irregular. It exploits the functions local regularity parameters that can be obtained from the wavelet coefficients. This thesis describes a relatively new method that makes use of the following two aspects: The selection of coefficients is based on an estimation of local regularity. It additionally takes into consideration a priori geometrical information of the image. Whereas most methods divide the coefficients into a noisy and clean category, and apply a different manipulation to each category, this method computes for each coefficient the probability that it is affected (unaffected) by noise. The manipulation of the coefficients is then based on this probability. This method combines the wavelet multiresolution concept with a Markov Random Field (MRF) model. In addition, the method describes a novel Bayesian approach for wavelet coefficients, (Malfait, Roose). 2 2. Brief Review on Wavelets 2.1 Transforms When data of any kind is recorded, it hardly ever reveals its true nature from our standpoint. Transforms are made on this data to perform the task of changing the standpoint. A basic transform in the physical world is, for example, a transform from the time-domain to the frequency-domain and vice-versa. Our ears are built is such a way so as to identify the frequencies and the phases of the sound signal. The frequencies are in turn used to analyze the message, and the phases carry information, which is mostly used to spatial hearing (low frequencies). Transforms thus have a high analytical value. 2.2 Basis Functions The idea of having a basis function is that any signal can be expressed as a weighted sum of a family of functions. These functions are the basis functions of the transform. Depending upon the features of the transform, the basis functions have to fulfill different requirements. Some of the requirements are listed as follows: Orthogonality: From an analysis point of view, it is practical that a set of functions is orthogonal. For a given set of functions, orthogonality on a given interval a functions is zero. Orthogonal basis functions describe the signal explicitly. In other words, there exists only one possible transform for each signal when such functions are used. They are a good choice for compression tasks. Redundancy: Not all transforms aim to analyze the signal afterwards. If we want to protect data against errors, it might be useful to add redundancy. In these cases, the choice of basis functions leads to redundancy, where only certain basis functions exist, and others are forbidden. In the case of an error, the basis function can be corrected to the nearest legal one according to the distance metric of the function space. Compact support: Compact support means that wherever the function is not defined it will have a value of zero. In numerical calculation, this is an advantage as all known zeros speed things up in calculations. Compact support finds applications in wavelet transforms, as it tells about the locality of wavelets in the time domain. 2.3 Fourier Transforms and the Time-Frequency Wedding The Fourier Transform is the most well-known frequency domain analysis tool. Fourier transforms provide a way to convert samples of a standard time-series into the frequency domain. This provides a dual representation of the function, in which certain operations become easier than in the time domain. Applications of Fourier transforms include filtering, image compression, convolution and deconvolution, and computing the correlations of functions. The uncertainty 4 principle states that the energy spread of a function and its Fourier transform cannot be arbitrarily small. In 1946, the physicist Gabor defined elementary time- frequency atoms as waveforms that have a minimal spread in a time-frequency plane. To measure time-frequency information content he proposed decomposing signals over these elementary waveforms. Windowed Fourier transforms and wavelet transforms (two prime examples of time-frequency transforms) are computed by decomposing the signal over different families of time-frequency atoms. A linear time-frequency transform correlates the signal with a family of waveforms that are well concentrated in time and in frequency. These waveforms are called time-frequency atoms, as illustrated in Figure 2.1. Let us consider a general family of time-frequency atoms {fr j^r, where y might be a multi-index parameter. Also, we suppose that The corresponding linear time-frequency operator T associates to any / e L2 (R) 5 (2.1) Tf(r)= //('Â¥; = The Parsevals formula (Mallat, 2001) proves that oo 1 oo T/0) = jfiOfy (0 = \/(<*>)% (co)dco. (2.2) -00 -00 If depends only on the values of / in this neighbourhood. Similarly, if (.f,(f>r ) reveals the properties of / in the neighbourhood ofÂ£ (Mallat, 2001). 2.3.1 Windowed Fourier Transform A Windowed Fourier atom is constructed by translating in time by u and modulating at the frequency Â£, a time widow g: g,^t) = e,9g(f-u) (2.3) It is normalized ||g|| = 1 so that the resulting Windowed Fourier transform of / e L2(/?)is oo S/(u,4) = (f,gu.f)= \f(t)g(t-u)e~,9dt. (2.4) -ao This transform is called the Windowed Fourier transform because the multiplication by g{t u) localizes the Fourier integral in the time window of a neighbourhood around /=. The resolution in time and frequency of the 6 Windowed Fourier transform depends on the spread of the window in time and frequency (Mallat, 2001). For image analysis, the Fourier transform identifies the spectral components present in an image but it does not provide information as to where certain components occur within the image. If we are only interested in stationary signals, the Fourier transform provides complete spectral analysis because all spectral components exist at all times. However, if we are interested in non- stationary signals with transient phenomena, the Fourier transform would only give us the spectral components within the image but not where they are located. The wavelet transform or wavelet analysis is probably the most recent solution to overcome the shortcomings of the Fourier transform. In wavelet analysis a fully scalable modulated window is shifted along the signal and for every position the spectrum is calculated. Then this process is repeated many times with a slightly shorter (or longer) window for every new cycle. The result is a collection of time- frequency representations of the signal, all with different resolutions or frequency bands. This collection of resolutions results in multi-resolution analysis. In case of wavelets we normally do not speak about time-frequency representations but of time-scale representations, where small scales correspond to high frequencies and large scales correspond to low frequencies. In fact, the human visual system performs hierarchical edge detection at multiple levels of resolution and wavelet transforms perform a similar multi-resolution analysis. 7 2.4 Wavelet Transforms Like a windowed Fourier transform, a wavelet transform can measure the time- frequency variations of spectral components, but it has a different time-frequency resolution. To analyze signal structures of very different sizes, it is necessary to use time-frequency atoms with different time supports, it is achieved by varying the scale s. The wavelet transform decomposes signals over dilated and translated wavelets. A wavelet is a function e L (R) with a zero average: f ViOdt = 0 J-oo (2.5) It is normalized 1^1 = 1, and centered in the neighbourhood of t=0. A family of time-frequency atoms is obtained by scaling y/ by s and translating it by u: t-u 1 y/uAt) = -j=Â¥ yJS K s J (2.6) These atoms remain normalized for all scales. The wavelet transform of / at the time u and scale s is dt.={f,Â¥u,s) (2.7) The wavelet transform can be rewritten as a convolution product: wf(u,s)= = (2.8) 8 with (2.10) Larger scales correspond to dilated wavelets and therefore lower frequencies. Smaller scales correspond to compressed wavelets and therefore higher frequencies. Since y/ has a zero average and unity norm, the above wavelet integral measures the variation of / in a neighbourhood of u, whose size is proportional to s. When the scale s goes to zero, it can be proved that the decay of the wavelet coefficients characterizes the regularity of/in the neighbourhood of u (Mallat, 2001). The resulting wavelet coefficients are known as the details of the signal because the wavelet is considered as a high pass filter. In order to perfectly reconstruct the signal from the wavelet coefficients, the low-frequency approximation must also be added to the high-frequency details (Mallat, 2001). 2.4.1 Wavelet Admissibility Condition (Mallat, 2001) 2.4.2 Scaling Function When Wf(u,s) is known only for s < s0, i.e., for high frequencies, to recover / we need a complement of information corresponding to s > s01 i.e., for low frequencies. This is obtained by introducing a scaling function aggregation of wavelets at scales larger than s0, where s0 corresponds to the (2.11) 9 mother wavelet^. The scaling function can be interpreted as the impulse response of a low-pass filter (proof can be found in Mallat, 2001). The low- frequency approximation of/ at the scale s is function ^(/) e LJ{R), to be a wavelet (Mallat, 2001). 2.5 Frames and Bases Frame theory analyzes the completeness, stability and redundancy of linear discrete signal representations. A frame is a family of vectors { characterizes any signal/from its inner products {{/,)} ner Definition (Mallat, 2001): The sequence { (2.12) The original signal can be recovered from the following equation: (2.13) A A lim where CÂ¥ is the admissibility condition, (2.14) When A=B the frame is said to be tight. 10 If the frame condition is satisfied, then U in equation (2.15) is called a frame operator. V6T,t//[] = (/,<*). (2.15) A frame thus defines a complete and stable signal representation, which may also be redundant. When the frame vectors are normalized ||^n|| = 1, this redundancy is measured by the frame bounds A and B. If the frame vectors { independent then it can be proved that A<\< B. The frame is an orthonormal basis if and only if A =5=1. This can be verified by inserting f=$ in equation (2.15) of the definition of a frame (Mallat, 2001). 2.5.1 Pseudo Inverse The reconstruction of / from its frame coefficients Uf[ri\ is calculated with a pseudo inverse. This pseudo inverse is a bounded operator that is expressed with a dual frame. I2(T) is the space of vectors such that ||x||2 = Xnerl^M2 < 00 anc* ^ ImU the image space of all Uf with fe H. Proposition (Mallat, 2001): Ij {(j)n} neT is a frame whose vectors are linearly independent, then ImU is strictly included in I2(r), and U admits an infinite number of left inverses U~': V/ e H,U~'Uf = f. (2-16) 11 The pseudo inverse of a frame operator is related to a dual frame family, which is specified by the following theorem. Theorem 2.1 (Mallat, 2001): Suppose that {tf>n} ner is a frame with frame bounds A and B. Let 2.5.2 Partial Reconstruction Suppose that { to fully recover f The best linear mean square approximation of/computed from these inner products is the orthogonal projection of / on the space V (Mallat, 2001). 2.6 Wavelet Frames Wavelet frames are constructed by sampling the time and scale parameters of a (2.17) P,f = U-'Uf = (Mallat, 2001) (2.18) continuous wavelet transform. A real continuous wavelet transform of / e L2(i?) is defined in previous sections by 12 (2.19) Where, y/ is a real wavelet and (2.20) Intuitively, to construct a frame we need to cover the time-frequency plane with the Heisenberg boxes of the discrete wavelet family. A wavelet y/us has energy in time that is centered at u over a domain proportional to s. Over positive frequencies, its Fourier transform y/us has a support centered at frequency q/s, with a spread proportional to 1/s. To obtain a full cover, we sample s along an exponential sequence {c/}, with a sufficiently small dilation step a> 1. In this thesis, a=2 is used. The time translation u is sampled uniformly at intervals proportional to the scale 2/. Let us denote 2.7 Wavelet Bases One can construct wavelets y/ such that the dilated and translated family of the t -nu0a (2.21) V a / mother wavelet can be expressed as an orthonormal basis of L2 (if). (2.22) 13 Orthogonal wavelets dilated by 27 carry signal variations at resolution 2y (Mallat, 2001). 2.8 Continuous Wavelet Transform (Mallat, 2001) The wavelet transform can be written as a convolution product: with, y/s(t) = r=\f/ Continuous wavelet transforms and windowed Fourier 4 s \ s ) transforms provide translation-invariant representations (Mallat, 2001). 2.9 Dyadic Wavelet Transform In computer programs, an input image is a discrete signal made up of pixels defined within a certain resolution. Therefore, the wavelet transform cannot be computed at an arbitrary fine scale. To construct a translation-invariant wavelet representation, the scale s is discretized but not the translation parameter u. The scale is sampled along a dyadic sequence^7}, to simplify the numerical calculations (Mallat, 2001). The dyadic wavelet transform offe L2(R)is defined by (2.23) 1 f-t (2.24) 14 with, V'2j(t) = V'2j (-0 = 4v Â¥ '-t' \2' j (2.25) 2.10 Wavelet Design A discrete dyadic wavelet transform can be computed with a fast filter bank algorithm if the wavelet is appropriately designed. Let h and g be a pair of finite impulse response (FIR) filters. Suppose that h is a low-pass filter whose transfer function satisfies h(0) = -Jl. We construct a scaling function whose Fourier transform is k h(2-pco) yfl 1 A Q) A 0) V2 2 2 (2.26) The corresponding wavelet y/ has a Fourier transform defined by (Mallat, 2001) (2.27) 1 A 1 'co' = f28 a, h <2; 2.10.1 Reconstruction Wavelets Reconstruction wavelets are calculated with a pair of finite impulse response filters h and g The following Fourier transform then has the following energy: #()n p=\ h(2-pa)) 1 r,o)~ a H-XK-) (2.28) 15 Also, ~ 1 ^ ^(^)=V2 8 co 'a>' {2.29) 2.11 Spline Dyadic Wavelets In the Conjugate Mirror Filter (CMF) algorithm, the continuous wavelet and scaling functions are replaced with discrete high pass and low pass filters respectively. The following equations show the relation between the discrete filters and their continuous counterparts (Mallat, 2001): (2.30) wri k}*-*} (2.31) The high and the low pass filters g[ri\ and h[n], respectively, are related in the following way: g[n] = (-irnh[\-n] (2.32) A box spline of degree m is defined as a translation of m+1 convolutions of 1[(U] with itself. It is centered at t= if m is even and at /= 0 if m is odd. Its Fourier 2 transform is 16 ( co^
ton
Theorem5.1 (proved in Mallat, 2001): Let { Vj}JeZ be a Multiresolution |