R E S T O R A T I O N OF I M A G E S D E G R A D E D B Y S Y S T E M S OF R A N D O M I M P U L S E R E S P O N S E by IVAN L. R E V E L A N T A THESIS S U B M I T T E D IN PARTIAL F U L F I L M E N T OF T H E R E Q U I R E M E N T S FOR T H E D E G R E E O F M A S T E R OF APPLIED SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Electrical Engineering We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C OLU MBIA 13 April 1987 © Ivan L. Revelant, 1987 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at The University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Electrical Engineering The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: 13 April 1987 A B S T R A C T The problem of restoring an image distorted by a system consisting of a stochastic impulse response in conjuction with additive noise is investigated. The method of constrained least squares is extended to this problem, and leads to the development of a new technique based on the minimization of a weighted error function. Results obtained using the new method are compared with those obtained by constrained least squares, and by the Wiener filter and approximations thereof. It is found that the new technique, "Weighted Least Squares", gives superior results if the noise in the impulse response is comparable to or greater than the additive noise. 11 C O N T E N T S A B S T R A C T ii LIST OF F I G U R E S v A C K N O W L E D G E M E N T vi 1 I N T R O D U C T I O N 1 2 T H E P R O B L E M 2 2.1 T H E D I S C R E T E C A S E . . 4 2.1.1 The Fourier Space 4 2.2 S C O P E 6 3 V A R I O U S R E S T O R A T I O N M E T H O D S 7 3.1 T H E WIENER F I L T E R 7 3.1.1 Applicability 8 3.2 T H E "GAMMA" FI L T E R 8 3.2.1 Applicability 9 3.3 A N O T H E R APPROXIMATION OF T H E WIENER F I L T E R . . . . 9 3.3.1 Applicability 9 3.4 C O N S T R A I N E D D E C O N V O L U T I O N .. 10 3.4.1 Applicability 12 3.5 W E I G H T E D L E A S T SQUARES ' 13 3.5.1 Generalization of the Method 15 iii 3.5.2 Analysis of the Result 17 3.5.3 Summary 20 3.5.4 Applicability 20 4 R E S U L T S 21 4.1 AN INTERESTING RESULT 22 4.2 COMPARISON WITH O T H E R M E T H O D S 23 4.3 O T H E R CHOICES FOR T H E A L T E R N A T E R O O T . 25 4.4 VISUAL RESULTS 28 5 C O N C L U S I O N 32 B I B L I O G R A P H Y 34 A P P E N D I C E S 35 A M A T H E M A T I C A L C O N V E N T I O N S 36 A . l PARSEVAL'S T H E O R E M 37 A.2 T H E C O N V O L U T I O N T H E O R E M 37 A.3 P O W E R S P E C T R A L DENSITY 38 B D E R I V A T I O N OF T H E C O N S T R A I N E D D E C O N V O L U T I O N F I L T E R 40 C D E R I V A T I O N OF T H E W E I G H T E D L E A S T S Q U A R E S FIL-T E R 44 C l OBTAINING T H E ROOTS 46 iv LIST OF F I G U R E S 3.1 Behaviour of the weighted residual function 18 4.1 Performance of the WLS filter for different values of o\ ....... . 22 4.2 Comparison of all methods using a square lens blur function 25 4.3 Comparison of all methods using a uniform square blur function . . 26 4.4 Performance of WLS using various alternatives to the inverse filter . 27 4.5 Results using oY - 0.001 and o2 = 0.316 29 4.6 Results using o~x — 0 and a2 = 0.316 30 4.7 Results using cr, = 0.001 and cr2 = 0 31 v A C K N O W L E D G E M E N T I wish to thank my supervisor, Professor R. K. Ward, for her guidance and patience during the course of my studies. I would also like to thank my colleague, Mr. Wilfred Bohm, for his encouragement and timely suggestions. Special thanks are due to my parents, Leonardo and Luigia Revelant, for their financial support, love, and encouragement throughout my academic career. vi C H A P T E R 1 - I N T R O D U C T I O N Many situations exist in which an image is obtained which has been degraded so that much of the information contained in the image cannot be perceived. This degradation may take either or both of the following forms: the image may be blurred, or noise may be introduced at various stages in the formation of the image, most commonly at the detection stage. Naturally, it would be desirable to reverse the degradation process as much as possible so as to recover the information which was apparently "lost". This is done by first rendering the degraded image in digital form (if it is not a digital image to begin with) so that it can be processed by a digital computer, then applying an algorithm which is designed to obtain a good estimate of the original image. This paper shall describe a number of commonly-used image restoration methods, present a new technique which the author has developed, and compare the results obtainable with each method. These methods will be applied in particular to restoration problems involving a partially stochastic blur, that is, to images which have been blurred by a system whose fundamental impulse response is known but which is also subject to noise. In this we depart from the usual scope of image restoration work, which assumes that the blur may be modelled in an entirely deterministic fashion. 1 C H A P T E R 2 - T H E P R O B L E M The purpose of the author's research was to find a new, hopefully better means of restoring images which have been degraded by a "noisy" blur in conjuction with additive noise. This degradation process may be modelled as follows [17]: g = (h + n i ) * f + n 2 (2.1) where g is the distorted image, f is the original (undegraded) image, h is a deter-ministic point-spread (blur) function, n i and n 2 are random processes with known statistical properties, and * is the convolution operator.1 In adopting this model, we assume that both the deterministic and stochastic components of the blur are limited in space and shift-invariant, and that the noise sources are uncorrelated with either the blur or the image. We also assume that no statistics are available concerning the original image. The random processes ni and n 2 may represent a variety of real-world phenom-ena. Examples of n i include: • Error in a posteriori estimation of h. Where the blur function is unknown or difficult to model (such as in the case of camera "jitter" or motion blur), one may try to estimate h by studying g. Obviously this can rarely be done perfectly. • Error in modelling h. If, for example, it is known that g was formed by an out-of-focus lens, there will still be some error in determining the degree to 1See Appendix A for an explanation of the notation used in this paper. 2 which the lens is out of focus. Furthermore, the actual transfer function of the lens is unlikely to conform perfectly to its idealistic model. • A natural random process may add itself to the known cause of blur. For example, an unknown amount of dust may appear between the object of interest and the lens, thus altering the blur function in a way that may be modelled statistically. • The blur function may be inherently random, such as in the case of atmo-spheric turbulence. Examples of 112 include: • Noise in the image detector. This is an obvious problem with electronic video equipment, but it also occurs as "graininess" in photosensitive films. • Quantization noise. When an image is rendered in digital form, some roundoff error in the brightness values is inevitable. Hence n 2 must always be assumed to exist when a digitized image is being processed (although it may, in some cases, be considered negligible). We shall be particularly interested in the case in which both n x and n 2 are white Gaussian noise processes. 3 2.1 T H E D I S C R E T E C A S E In order to find a means of restoring a digital image, we first reformulate the degradation process as follows:2 G = B * F + N , * F + N 2 (2.2) Here we have replaced the continuous functions of Equation 2.1 with matrices ob-tained by sampling the corresponding continuous functions on a rectangular tesse-lation. We shall define their dimensions as follows: • F has dimensions M x N; • The blur matrix B and its associated noise N j both have dimensions J x K; • G and the additive noise N 2 both have dimensions (M + J — 1) x (N + K — 1), since the size of G is determined by the convolution of the image and blur matrices. 2.1.1 The Fourier Space In the literature, 3 image restoration problems (not necessarily involving stochastic blur) have often been approached using existing one-dimensional signal processing techniques (which are applied in the time domain), by converting the image to a one-dimensional signal 4 and finding a solution by methods involving the solution of systems of equations. This approach becomes extremely costly when applied 2The convolution operation, as applied to matrices, is defined in Appendix A. 3See, for example, [l], [12], [13), and [17] 4 A "raster scan" is most often used. 4 to images because of the large size of the matrices involved. The cost and com-plexity of restoration may be considerably reduced by considering the problem in the frequency domain, since space-domain convolution operations are replaced by pointwise linear operations in the frequency domain. One way to convert signal processing techniques to the Fourier domain leaves the problem formulated in one dimension, but expresses the restoration filter in terms of circulant matrices. The Fourier version of the filter is then obtained directly using the fact that the eigenvalues of a circulant matrix are equal to the discrete Fourier transform coefficients of the circulant sequence. First used by Hunt [8], this technique is undeniably clever, but we shall find it easier to obtain restoration filters directly in the frequency domain. We may reformulate our problem in the Fourier space as follows: Q = ?QB7 + + M2 (2.3) where §, B, 7, A/i, and M2 are the discrete Fourier transforms of the matrices G, B , F , N i , and N 2 , respectively, of Equation 2.2 (P and Q are discussed below). To make use of Fourier techniques, we must first expand the matrices of Equa-tion 2.2 to a common size, which must be at least (M + J — 1) x (N + K — l). This is done by padding the smaller matrices with zeros. This done, the discrete Fourier transform (DFT) is applied to the degraded image and the blur model to obtain their spectra; a restoration algorithm derived from Equation 2.3 is applied to ob-tain the spectrum of the image estimate T; and the inverse D F T is applied to J to obtain F . Usually, the matrices are expanded to a larger size than is strictly necessary -most commonly a power of 2 on each side - in order to take advantage of fast Fourier transform algorithms. Hence, the size of the operand matrices will, in general, be greater than (M + J — 1) x (N -f K - 1). Henceforward, we shall use the dimensions P x Q for all such matrices. 2.2 S C O P E We shall be expressly concerned here with techniques for the restoration of images, as opposed to their enhancement. This means that we shall only consider techniques which we expect to reverse the degradation process, not those which improve the appearance of images by means of techniques such as edge enhancement, contrast adjustment, smoothing, and so forth. To decide on the success of a restoration method, we need an objective criterion of quality. The most commonly used such criterion is the mean-squared error (MSE) of the estimate, that is, the average of the squares of the differences between each point in F and the corresponding point in F. We shall adopt this criterion due to its popularity and its mathematical simplicity. Another, less objective, criterion is sharpness, which we shall loosely define as the amount of genuine detail visible in the restored image. 6 C H A P T E R 3 - VARIOUS RESTORATION METHODS 3 .1 T H E WIENER FILTER The Wiener filter, or least squares filter, is obtained by minimizing the expected value of the mean-squared error of the estimate: 1 p " l Q _ 1 \ 2 ' minimize E - ^ ^ (Fpq - Fpq) (3.1) _ r V * p=0 q = 0 In general, the solution of (3.1) is a nonlinear function of the degraded image G , and requires knowledge of the joint probability density of F and G , treating both as random fields [13]. To simplify the problem, we restrict F to be a linear function of G , ie. Jpq = Mpq^pq, where M is not a function of either J or Q. Minimization of (3.1) then yields the optimum linear least-squares estimate, which for our problem is: "q P 2 Q 2 | B P J 2 + P2Q^ni(p,c7) + | ^ where Sni and Sn2 are the power spectral densities of N i and N 2 , respectively, and S/f is the power spectral density of the undegraded image F. Note that this means that the original image must be treated as a random field. Also notice that if Sni is zero, that is, if noise is absent from the blur, then the familiar form of the Wiener filter results. The Wiener filter was first applied to images by Helstrom [7] for the case of additive noise only, and by Slepian [14] for the case of an entirely unknown blur (but with known statistics). Equation 3.2 is adapted from Franks [4], where it is applied to one-dimensional signals. Ward and Saleh [17] have used this filter with our imaging model in its space-domain formulation. 3.1.1 Applicability The Wiener filter is the best linear filter in the sense that it minimizes the expected mean-squared error of the estimate. However, it requires some a priori knowledge of F , which is treated as a random field with a known power spectral density, so that it is most often applied to classes of images rather than single pictures. We are interested in restoration of a single image, with no knowledge of its statistics, so the Wiener filter is not strictly applicable to our problem. However, it shall be useful to us as a standard optimum linear filter against which other methods may be compared. It should be stressed that the optimality of this filter is confined to the class of linear filters, and that nonlinear filters have been found which outperform the Wiener filter [17] at the cost of increased complexity of implementation and time required to obtain a solution. 3.2 T H E " G A M M A " F I L T E R In many restoration problems, either the expected power spectrum of the image, or the statistics of the noise, or both, may be unknown. In such cases, although the Wiener filter cannot be considered strictly valid, the following approximation to Equation 3 .2 is often used [13]: n P 2 Q 2 |B P 9 | 2 + r 1 ' Here the noise-to-signal power ratio P 2 Q 2 S n i + SNJSJJ is replaced by the constant T, which is generally chosen by making a rough guess of the overall NSR. (In the 8 literature, where this filter is seen, the noise in the blur is generally absent, so that T just replaces SnJSff.) 3.2.1 Applicability Rather good results (comparable to Wiener) can be obtained with this method by trying different values of T until the best-looking results are obtained. Of course, this adds greatly to the amount of time needed to restore an image. 3.3 A N O T H E R A P P R O X I M A T I O N O F T H E W I E N E R F I L T E R If the noise statistics are known, but not the expected image power spectrum, a "quick-and-dirty" estimate can be obtained by replacing the expected power spectrum of the image in the Wiener filter by the actual power spectrum of the degraded image, as follows: on the assumption that the power spectrum of the degraded image is not far different from that of the original image. We shall refer to Equation 3.4 as the "false Wiener" filter. 3.3.1 Applicability This filter works surprisingly well despite the fact that it is strictly ad hoc. It is (so far) probably the best one-shot linear filter in cases where nothing at all is known about J, and is particularly good in cases where the additive noise is (3.4) P 2Q 2 |3 p , | 2 + P 2Q 2 S n i (?,<?) + Sn2{p,g) 9 small compared to the multiplicative noise, since it then approaches the form of the Wiener filter. 3.4 C O N S T R A I N E D D E C O N V O L U T I O N Constrained Deconvolution was first developed for the one-dimensional case by Phillips [ll] and applied by Twomey and Hunt to the problem of image restoration with additive noise only (nx = 0) [16],[9]. It was derived directly in the frequency domain, again for the case of noiseless blur, by Dines and Kak [3]. I have adapted the method of Dines and Kak to our case of a partially stochastic blur. What follows here is a quick overview of the method; for a detailed derivation, please see Appendix B. The idea of this method is to find an estimate F such that the total computed noise power (known as the "residual") is equal to its expected value: E E k - ( B * F ) P J 2 = £ E D ( N i 4 + N j 2 (3.5) ' p=0 9=0 p=0 ?=0 Some manipulation is required to convert this equation to its frequency domain analogue for the case of white Gaussian noise:1 P - 1 Q - I r " - (M + J - l ) ( N + K - l ) a 2 2 E E p=0 (j=0 I&9 - P Q V J 2 - J K ( 7 1 4 P Q where o\ and o\ represent the variances of the multiplicative noise (noise in the blur) and the additive noise, respectively. Since an infinite number of estimates exist which conform to this constraint, another independent criterion is chosen to find an optimum estimate. For example, 1 Please see Appendix B. 10 one might choose to make the estimate smooth, or to minimize its total energy. The criterion used in this method is to minimize a function of the form: p - i Q - i p=0 17=0 (3.7) since L can very easily be chosen so that minimization of (3.7) is helpful. For example, the Laplacian operator may be chosen: L = 0 1 0 1 -4 1 0 1 0 0 <Z> 0 (3.8) in which case, by minimizing (3.7), we are choosing the version of F with the smallest sum of second differences, and hence a "smooth" F. Note that L is padded with zeros in the same manner as B, so that the convolution in (3.7) is circular. Usually the dimensions of the active part of L are smaller than those of B; if not, it may be necessary to increase the dimensions P and Q to accomodate L * F . Another common choice is to make L a simple impulse: Jpq 1, ifp = q = 0 0, otherwise (3.9) whereby (3.7) minimizes the total energy of F. This is the form of L which we shall be using (see Chapter 4, "Results"). In any case, the method of Lagrange multipliers is used to minimize the frequency-domain form of (3.7), with j constrained by Equation 3.6. This leads to the fol-11 lowing filter: f = P Q & A , 3 1 0 l P q u P 2 Q 2 | £ p g | 2 + P 2 Q 2 | B p , | 2 - mo\ where v is an unknown constant derived from the Lagrange multiplier A. It is impossible to (analytically) solve (3.10) and (3.6) simultaneously to obtain v, so an iterative procedure must be used. A guess is made of the value of v, Equation (3.10) is used to obtain the corresponding J , then this is plugged into (3.6). This is repeated until Equation 3.6 is satisfied. It turns out [13] that the left side of (3.6) is a monotonically increasing function of v in the region of the solution, so we may easily "zero in" on the value of v for which Equation 3.6 is satisfied.2 3.4.1 Applicability Constrained Deconvolution is an example of a nonlinear restoration method. Since its solution cannot be obtained analytically, it is time-consuming compared to a linear method such as the Wiener filter. However, in contrast with the Gamma filter, which also requires that a parameter be adjusted until the "best" solution is found, Constrained Deconvolution obtains its solution automatically, with no need for human judgment to choose an optimum solution. Its advantage over the Wiener filter is that nothing needs to be known a priori about the image to be restored; it is merely necessary to make the assumption that the two noise sources are uncorrelated with each other and with the image, which is usually the case. (Equations 3.6 &; 3.10 are derived for the particular case of stationary white Gaussian noise sources.) 2This is proven in Reference [13] without considering N i . Because the constraint is nonlinear in our case of a noisy blur, we cannot prove that it is still true, but the author has found it to be true in practice. 12 3.5 W E I G H T E D L E A S T S Q U A R E S What follows here is the history of the author's development of the Fourier-domain version of the weighted least-squares (WLS) filter, which began with a modification of Constrained Deconvolution for the particular case of white Gaussian noise. Since the expected value of the residual, E jl2p=o ^q=o Gpq — (B * F) p, j , is dependent on F when o\ ^ 0 (Eq. 3.5), it is possible to minimize the residual itself (and hence, optimistically, the noise in the estimate), rather than a chosen linear function of F , while still keeping the same constraint: that the residual be equal to its expected value. Stating our objective in the frequency domain, we wish to minimize: P - i Q - i p=0 q=0 subject to Equation 3.6. Following a similar derivation to that of the constrained deconvolution filter, we obtain: / p ? _ P 2 Q 2 | B P 9 | 2 + t;JKa 2 l ^ i l J where, again, v is found by iteration with Equation 3.6. It is easily shown that this filter is the same as the constrained deconvolution filter using the minimum energy criterion, unless of course o~i = 0, in which case it degenerates to the simple inverse filter. (The alert reader may also have noticed that this amounts to a determinate "Gamma" filter, ie. it represents a means of choosing T objectively foj: the case of white Gaussian noise.) At first glance, minimization of the residual appears to be a good criterion for finding F. However, it fails to take into account the fact that the noise is signal-13 dependent: it is larger where F is larger. This means that the noise energy is greater in areas3 where F is stronger. Thus, the minimum-residual criterion is weighted in favour of areas where the signal strength is greatest - precisely where there is more noise. Clearly, in order to correct for this behaviour, we should give more weight to the residual at frequencies where J is small. We may do this by minimizing a weighted version of the residual, in such a way that the expected value of this weighted residual is constant for any magnitude of The obvious choice is to weight the residual by its expected value, which leads us to the following criterion: P - l Q - l Q - P O B f 2 minimize > > ~ - (3.12) J ( M + J - l ) ( N + K - l ) c 72 T p 2Q2 p=o g=o J K a f Obviously, it makes no sense to apply the constraint of Equation 3.6, which would simply set this quantity equal to 1. Hence we simply minimize this with respect to f, obtaining:4 ' f = ''T , (3.13) P 2Q 2 ISI JKCT S-PQ07 JKa2|/|- + ( M + J - l ) ( N + K - l ) a 2 / P 2 Q 2 (Here, and henceforward, the subscript pq is omitted from frequency-domain ex-pressions, as it is clear that such expressions are to be evaluated point-by-point.) If a reasonable starting estimate of J is chosen, iteration of Equation 3.13 will converge successfully. However, it is possible to solve this equation analytically using the fact that the phase of f is equal to that of CjB* (from inspection of Equation 3.13). By this means, Equation 3.13 can be rearranged into a quadratic 3Either areas in space, or frequency components. 4Please see Appendix C for a detailed derivation of a generalized form of this filter. 14 equation, the roots of which are: [ 9 - g * ( M + J - l ) ( N + K - l ) g ; J ~ \ P Q S ' P Q J K a 2 £ * j ( 3 ' 1 4 ) Hence, the solution at any point in the spectrum is either the I N V E R S E F I L T E R or a value which is 180° out of phase wi th the inverse filter. In a t r ia l , iteration of Equat ion 3.13 starting wi th the D F T of the undegraded image yielded an estimate which was composed of one or the other of the two roots of (3.14) at each point. Taking a sample of this estimate, it was found that, for the samples taken, the roots obtained were the ones which best fit the constraint of Equat ion 3.6 at each point. We now have a useful, semi-linear filter. The scheme is to obtain both possible roots (3.14) at each point in the spectrum, and choose the one which best fits the point form of the constraint on j : |S - PQB/f = JK.? j/f + (M + J - y - l M ( 3 , 5 ) 3.5.1 Generalization of the Method Minimiza t ion of a weighted residual is the basis of a well-known space-domain technique known as "regression spatial image restoration" [12]. This technique models the degradation of a one-dimensional signal as follows: g = Bf+n (3.16) A minimum-variance unbiased estimate of / can be obtained by minimizing the following function: (g-BfYK-'ig-Bf) (3.17) 15 where K " 1 is the inverse of the autocovariance matrix of the noise. The idea is apparently to minimize the residual, weighted at each point by the variance of the noise. No means has been found to implement this method in the frequency domain, so that its application to pictures is limited due to the enormous size of the matrices involved and the time required to obtain a solution, even when it is formulated in terms of circulant matrices [6]. We now realize that our method is a frequency-domain analogy to the mini-mization of (3.17). In our method, the Fourier version of the residual is minimized, weighted at each point in the spectrum by the power spectral density of the noise, which we define as follows: Sn(p,q) = E[\M\2(p,q)} (3.18) = [DFT{R n }] p g (3.19) For Equation 3.19 to hold, it is necessary for the noise to be stationary, so that R n (the autocorrelation function of the noise) is a function of two dimensions. (Otherwise, it is a function of four dimensions - the location vectors of the two points being correlated.) In any case, observe that the power spectral density of the noise is also equal to the expected value of the residual: Sn{p,g) = E [\gpq - P Q B P 9 J p g | 2 ] . (3.20) This leads us to the following generalized Weighted Least Squares criterion: m i n i m i z e 2 P = 0 g = o sn{p,q) 16 I In the degradation model we have adopted, the additive and multiplicative noise sources are independent of each other and of F, so that (as shown in Appendix A): S n = P 2 Q 2 S n i | J | 2 + S„ 2 (3.22) so that (3.21) leads us to: / = (3 23) p 2 Q , | B | 2 P 2 Q 2 s , , l g - P Q a / 2 and hence: J = { p Q S ' P ^ ^ } ( 3 - 2 4 ) where the correct root is chosen to be that which best fits: | £ - P Q B 7 | 2 = P 2 Q 2 S n i | 7 | 2 - r - S „ 2 . (3.25) Please see Appendix C for the complete derivation of this filter. 3.5.2 Analysis of the Result We have not yet verified that (3.24), the roots of Equation 3.23, actually minimize (3.21). Let us consider this now. Inspection of (3.21) immediately reveals that the quantity to be minimized is a nonnegative real number for all (p,q). Now, if we substitute the first root from (3.24), , (the "inverse filter") into (3.21), we get 0. Substitution of the second root returns a nonzero positive quantity, and we can see A IB I" that, in the limit as | J\ —> oo, (3.21) approaches Since the only finite extrema of (3.21) are given by Equation 3.24, we can only come to one conclusion. Please see Figure 3.1. 17 I ^ - P Q ^ I / X Figure 3.1: Behaviour of the weighted residual function for all estimates 7 at a given point (p, q). (Cross-section through the line 7 = | J|e J a r g^ s*).) In fact, the second root from Equation 3.24 maximizes the weighted error func-tion (3.21). Yet the method described in the previous section works (as will be seen in the next chapter), and works considerably better than the inverse filter, which is the absolute minimum of (3.21). To explain this, we need to consider the inverse filter 7 = pqg, the most primi-tive attempt at image restoration. It is obtained by taking the expectation of both sides of Equation 2.3: E\9\ = E [PQSf + PQ>/iJ + M2] = PQB 7 (3.26) which immediately gives 7 = pqg-This filter performs poorly because in most cases the blur transfer function 8 tends to act as a low-pass filter when applied to 7. Thus, when noise is added to a blurred image, it accounts for most of the higher frequencies in the spectrum of Q. Now, since B is a low-frequency emphasis filter, its inverse is a high-frequency emphasis filter. Due to the difference in the spectra of the blurred image and the noise, such a filter amplifies the noise more than the signal component of Q, 18 so that the "restored" image is very noisy, containing very large high-frequency components. One attempt to circumvent this problem in the application of the inverse filter is described by McGlamery [10]. This method selects a threshold frequency below which the inverse filter may be applied. Above this threshold, the inverse is replaced by the degraded image § . Until now, the selection of this threshold has been a matter of trial and error. Now consider the operation of the "constraint" (Equation 3.25). It selects the root in (3.24) for which the error function \Q — PQS 7\2 is closest to its expected value. We may think of the action of the selection process as follows: "If the estimate 7 which maximizes the weighted residual has a noise power which is closer to the expected quantity than that of the inverse filter (which gives the absolute minimum of the weighted residual), then it appears that the inverse filter fails at that frequency." It is our great fortune that the alternative to the inverse filter given by the other root in (3.24) is small at such points compared to the inverse filter (since their magnitudes are inversely proportional) so that, in effect, our method of selecting between the possible solutions (3.24) is an automatic implementation of McGlamery's thresholding technique. Thus, by a somewhat roundabout method, we have arrived at a determinate means of deciding whether or not to use the inverse filter for each point of the spectrum of the estimate 7: we see whether it fits the constraint (Equation 3.25) — B* S better than the maximizing root pqg.p • If not, we substitute something else. It stands to reason, now, that the optimal "something else" is not the maxi-19 mizing root; after all, it is 180° out of phase with the inverse filter. If our theory is correct, then using, for example, § as the alternative (following McGlamery), or 0, should give better results. This hypothesis has been confirmed; results using various alternatives to the inverse filter are shown in the next chapter. 3.5.3 Summary Our new method, "Weighted Least Squares", boils down to a non-arbitrary means of deciding where to use the inverse filter as an estimate of the undegraded image spectrum. The decision is made by finding which of the two quantities (3.24) best fits the following model of the noise power at each point in the spectrum of the degraded image: \9-?QBT\2 = V2Q2Sni\7\2 + Sn2 If the inverse filter, pqg, fits best, it is used to obtain 7. If not, another quantity is substituted; for example, 0 or p. 3.5.4 Applicabil ity WLS may be used in any image restoration problem which can be described by our model (Equation 2.3) and for which it is possible to model the noise sources by their power spectral densities. Whether the method is effective (gives good results) remains to be seen, and is the topic of the next chapter. 20 C H A P T E R 4 - R E S U L T S The methods presented in the previous chapter were applied to the restoration of an image of size 120 x 120 (the "FACE" of Figure 4.5(a)) which was degraded by a simulated optical system with the following impulse response: h(p) = p^>)' (4.1) where p is the distance from the origin. This represents a diffraction-limited in-coherent optical system of square aperture [2]. This function was sampled and truncated to a size of 9 pixels square, with the scale of p chosen to give 3 complete lobes in all. The resulting matrix (B) was then normalized so that the sum of its elements was equal to 1 (so that convolution of this matrix with the target image would preserve the energy of the image). Hence G had dimensions of 128 x 128, so that no further padding was necessary for the application of the fast Fourier transform. The noise matrices N i and N 2 were generated by stationary white Gaussian pseudo-random number generators with variances o\ and a|, respectively. This system was fine for numerical comparison of the various restoration meth-ods, but for the pictures of Figs. 4.5 - 4.7, a stronger blur was needed to bring out the differences between the various techniques. Therefore, a uniform 9 x 9 step function was used for the photographs. The original image "FACE" is an array of eight-bit positive integers. Its range is (4, 227), its mean value is 106.53, and its sample standard deviation is 43.35. Prior to processing, it is converted to an array of real numbers and padded with its mean (rather than zeros) to size P x Q. 21 4.1 A N I N T E R E S T I N G R E S U L T The author was fortunate to have first tried the WLS filter on an image with a relatively large value of ox. Indeed, the method would have been dismissed entirely on the sole basis of its performance for the case of no multiplicative noise. The performance of the WLS filter is shown in Figure 4.1. Figure 4.1: Performance of the WLS filter for different values of o\, with o\ = 0.1. The filter uses either the actual value of o\ or a fixed minimum value. Note how the mean-squared error of the estimate actually increases as the noise in the blur decreases below a certain point. This is contrary to what we might expect, but can be understood if we refer to Equation 3.13 - as ox approaches 0, Equation 3.13 approaches the inverse filter, which of course fails due to the presence of additive noise. „In order to force this filter to behave well for low values of Oi, we must set a minimum value for ox to be used by the filter. It was found that, for optimum results, this minimum value is approximately proportional to CT2 and, for a given cr2, M S E ( F ) 10-1-I I I ! 1 I I I 0 1E-8 1E-7 1E-6 1E-5 1E-4 0.001 0.01 0.1 2 22 is constant for any quantity of ox below the minimum value. However, the minimum value changes for different types of blur (it seems to drop as the blur becomes broader) and different images, so that, in practice, the WLS filter cannot be used where o\ is small without finding the minimum ax manually for each restoration problem (much as is done with the "Gamma" filter). 4.2 C O M P A R I S O N W I T H O T H E R M E T H O D S In order to assess the value of the WLS filter, the mean-squared error of the estimate using WLS was compared to those obtained using the other techniques presented in the previous chapter. The Wiener filter was selected as the optimum standard, using the actual value of | J | 2 in place of the power spectral density Sff (of course, this would never be used in practice!). It should be stressed that the Wiener filter is not used here as an actual candidate for use in our problem, since we have assumed that nothing is known about the image prior to its restoration. The "Gamma" filter' was chosen as an example of a linear filter which must be "tuned" manually to obtain an optimum result. The "false Wiener" filter is an example of a linear filter which does not require such tuning, and Constrained Deconvolution (CDV) is an example of a nonlinear filter. All other things being equal, we would choose a linear filter over a nonlinear one because it is faster; and any filter which works on its own is preferable to one which must be adjusted manually. From the point of view of implementation, the WLS filter may be considered linear (since it does not require iterations, but simply selects one of two possible 23 linear solutions using a simple criterion). For large values of ox it works on its own, but for smaller values it requires manual adjustment. Thus, to be useful, it must give better results than all other filters (except Wiener) when tuned manually, or it must give better results than the "false Wiener" filter and C D V when the actual value of oi is used by the filter. For the case of o2 — 0, it was found that the "false Wiener" filter always outperforms WLS. This is to be expected, since in this case the false Wiener filter is functionally identical to the Wiener filter and hence gives the optimum linear estimate. When 0\ — 0, a minimum value must always be found for o\ as used by WLS; in such cases, the Gamma filter always gives better results for the same effort. In fact, the author has found that the Gamma filter always gives better results (when properly "tuned") than WLS (and even better than the Wiener filter in some cases). Thus, we shall no longer consider the use of a minimum value of a x with WLS. Figures 4.2 & 4.3 show the performance of all the candidate methods for a fixed value of a 2 (0.1 was chosen as a rough approximation to the amount of additive noise one might expect from quantization error). The MSE of each estimate is shown as a ratio to the MSE of the corresponding Wiener estimate. Constrained Deconvolution is used here with the minimum energy criterion (L is given by Equation 3.9). WLS2 is a variation of WLS which replaces the negative root of Equation 3.13 with 0. As the reader can see from the figures, this variation is superior to the original WLS filter. From these plots, it is apparent that a range of values of ox exists in which 24 M S E (F ) / M S E (Wiener) (dB) 5.0 2.5-0--1.5 False Wietier CDV Gamma I I 1 1 r 1E-8 1E-7 1E-6 1E-5 1E-4 0.001 Figure 4.2: Comparison of all methods using a square lens (sine2) blur function WLS2 gives superior results compared to both C D V and "false Wiener". This is the region in which WLS2 may be considered useful, as it gives the best results of the available "automatic" methods. The Gamma filter still gives better results in this range of o\, but must be adjusted manually. 4.3 O T H E R C H O I C E S F O R T H E A L T E R N A T E R O O T We have already seen that better choices than the "maximizing root", p n ''2 exist as alternatives to the inverse filter for use by WLS. The question is, what would the best alternative root be? Unfortunately, the author has been unable to analytically determine an optimal choice. However, a number of obvious choices present themselves: 25 MSE(F)/MSE (Wiener) (dB) WLS2' 5.0' 2.5-0 --1.5' False Wiener CDV Gamma ' -> Hiug" — — » - - - - • 0 1E-8 1E-7 1E-6 1E-5 1E-4 OM 001 01 " ^ 1 Figure 4.3: Comparison of all methods using a uniform square blur function • Zero, as used by WLS2. This assumes that, where the inverse filter is not selected, the noise-to-signal ratio in the degraded image at those points in the spectrum is greater than unity, so it is best to "leave 'out" those points altogether. • Q, with its phase changed to that of the inverse filter. The phase of the inverse filter is used simply because the Wiener filter always uses it, so that we expect that phase to be the optimum choice. • The following approximation to the Wiener filter: VQ9B* P 2 Q 2 | B | 2 + P 2 Q 2 5 n i (4.2) 26 For each of these alternatives, we may use Equation 3.25 to compare the inverse filter either to the maximizing root or the alternative which we actually intend to use. It has been found that the latter approach gives superior results in all cases tried. Figure 4.4 shows the results obtained using each alternative to the inverse filter. Of these, it is apparent that the Wiener approximation works best for large values of Figure 4.4: Performance of WLS using various alternatives to the inverse filter, using a uniform blur function, with al — 0.1. The M S E of the estimate obtained using each method is displayed as a ratio to the MSE of the "false Wiener" estimate. <7i, and indeed outperforms the false Wiener filter for all o\ > 10~5 (approximately). Q works better for smaller values of a\. All the methods give poor results for very low amounts of noise in the blur. It should be noted that the improvement in M S E obtained by using § instead of 0 (as was used by WLS2) comes at the expense 27 of some detail in the restored image. Hence, if sharpness of the restored image is important, it is preferable to use 0 rather than § for levels of multiplicative noise at which these two alternatives are better than the Wiener approximation. 4.4 V I S U A L R E S U L T S Figures 4.5-4.7 are photographs showing the best results obtainable using each method discussed, applied to a uniform 9 x 9-pixel blur. In Figure 4.5, o\ = 0.001 and a2 = 0.316. Here the Gamma filter gives results as good as the Wiener filter. C D V and the false Wiener filter give somewhat blurred output, while both versions of WLS give a good, sharp image, though noisier than that produced by the Gamma filter. WLS2 may be considered useful here, since ox > 0 i ( m m ) - (This is not true of the original WLS filter.) In Figure 4.6, ox — 0 and cr2 = 0.316. Here the Gamma filter gives results which appear in these photos to be better than those of the Wiener filter. In fact, the Gamma filter has a "grainy" artifact which is not as visible in the photographs as the Wiener filter's artifacts. The false Wiener filter gives terrible results. C D V gives quite good results, and may be considered the best "automatic" filter. Results using WLS are poor, especially in view of the fact that the minimum ox must be used. In Figure 4.7, 0\ = 0.001 and o2 = 0. Again, the Gamma filter gives excellent results. This time, the false Wiener filter is the same as the Wiener filter, so naturally their results match. C D V works rather poorly. WLS is better, but would still not be used since the false Wiener filter works so well. 28 h Figure 4.5: Results using ox = 0.001 and o~2 = 0.316: (a) The original image "FACE", (b) "FACE" blurred and with noise added, (c) Output of Wiener filter, (d) Output of "false Wiener" filter, (e) CDV. (f) WLS. (g) Gamma filter, (h) WLS2. 29 6L c f Figure 4.6: Results using ox = 0 and o2 - 0.316: (a) The original image "FACE", (b) "FACE" blurred and with noise added, (c) Output of Wiener filter, (d) Output of "false Wiener" filter, (e) C D V . (f) WLS. (g) Gamma filter, (h) WLS2. 30 Si my. f K Figure 4.7: Results using <7j - 0.00] and a2 - 0: (a) The original image "'FACE". (b) "FACE" blurred and with noise added, (c) Output of Wiener filter, (d) Output of "false Wiener" filter, (e) CDV. (f) WLS. (g) Gamma filter, (h) WLS2. 31 C H A P T E R 5 - C O N C L U S I O N We have adapted and applied the following image restoration techniques to images degraded by a stochastic blur and additive noise: the inverse filter, the Wiener filter, the "Gamma" filter, the "false Wiener" filter, and Constrained Deconvolu-tion. Investigation of the noise-power constraint, which is the basis of Constrained Deconvolution, has led to the development of a new restoration technique. This technique, Weighted Least Squares, is a modification of the inverse filter whereby the points at which the inverse filter gives good results are determined analyti-cally. Where the inverse filter is predicted to give poor results, one of a number of replacement functions is used in its place. Weighted Least Squares has been found to be useful where the effect of noise added to the blur impulse response is comparable to, or greater than, that due to additive noise. It is applicable to any system which can be described by the degradation model which was presented in Chapter 2, and for which the power spectral densities of the noise sources are known. Four alternatives to the inverse filter were investigated for use by WLS. The first (the "maximizing" root from (3.24)), was shown to be suboptimal. The degraded image spectrum and 0 were found to work well for moderate amounts of blur noise. For larger amounts of this noise, an approximation of the Wiener filter was found to give the best results. Although nonlinear filters may exist which give better results, they have the dis-advantage of requiring considerably more computation time than does the Weighted Least Squares filter, which is "semi-linear" and produces results within a timespan 32 comparable to that required by the fastest linear filters. The author is optimistic that the technique of choosing between alternative solutions of an imaging problem by using the noise-power constraint (Equation 3.6) may be found useful in other methods of image restoration. 33 B I B L I O G R A P H Y [1] H. C. Andrews and B. R. Hunt, Digital Image Restoration. Englewood Cliffs, N.J.: Prentice-Hall Inc., 1979. [2] K. R. Castleman, Digital Image Processing. Englewood Cliffs, N.J.: Prentice-Hall Inc., 1979. [3] K. A. Dines and A. C. Kak, Constrained Least Squares Filtering. IEEE Trans. Acoust. Speech Signal Processing ASSP-25, 1977, pp. 346-350. [4] L. E. Franks, Signal Theory. Englewood Cliffs, N.J.: Prentice-Hall Inc., 1969. [5] R. C. Gonzalez and P. Wintz, Digital Image Processing. Reading, Massachusetts: Addison-Wesley, 1977. [6] L. Guan and R. K. Ward, Restoration of Randomly Blurred Images Using the Wiener and the Minimum Variance Unbiased Techniques. Submitted to: IEEE Trans. Acoust. Speech Signal Processing. (Not yet published) [7j C. W. Helstrom, Image Restoration by the Method of Least Squares. J. Opt. Soc. Amer. 57, March 1967, pp. 297-303. [8] B. R. Hunt, Deconvolution of Linear Systems by Constrained Regression and its Re-lationship to the Wiener Theory. IEEE Trans. Automat. Control. AC-17, 1972, pp. 703-705. [9] B. R. Hunt, The Application of Constrained Least Squares Estimation to Image Restoration by Digital Computer. IEEE Trans. Comput. C-22, 1973, pp. 805-812. [10] B. L. McGlamery, Restoration of Turbulence Degraded Images. J . Opt. Soc. Amer. 57, 1967, p.295. [11] D. L. Phillips, A Technique for the Numerical Solution of Certain Integral Equations of the First Kind. J.ACM 9, 1962, pp. 84-97. [12] W. K. Pratt, Digital Image Processing. New York: John Wiley & Sons, Inc., 1978. [13] A. Rosenfeld and A. C. Kak, Digital Picture Processing. Orlando, Florida: Aca-demic Press, 1982. (Second Edition) [14] D. Slepian, Linear Least-Squares Filtering of Distorted Images. J . Opt. Soc. Amer. 57, July 1967, pp. 918-922. [15] G. B. Thomas and R. L. Finney, Calculus and Analytic Geometry, Part II. Reading, Massachusetts: Addison-Wesley, 1979. (Fifth Edition) 34 [16] S. Twomey, The Application of Numerical Filtering to the Solution of Integral Equa-tions Encountered in Indirect Sensing Measurements. J . Franklin Inst. 297, 1965, pp. 95-109. [17] R. K. Ward and E. A. Saleh, Restoration of Images Distorted by Systems of Random Impulse Response. J . Opt. Soc. Amer. A, Vol. 2, Number 8, August 1985, pp. 1254-1259. 35 A P P E N D I X A - M A T H E M A T I C A L C O N V E N T I O N S If we ignore colour, an image can be simply described as a function of two di-mensions, strictly limited in domain. In this text, a continuous image (such as a photograph or the real image formed by a lens) will be denoted by a normal, lower-case letter, eg. f(x,y). A digital image takes the form of a matrix, which will be denoted by a boldface upper-case letter, and an individual element of such a matrix (a pixel) will be indicated by subscripts, as follows: F zy Where used, vectors will be denoted by lowercase letters, and will be obvious by the presence of an arrow superscript. The dimensions of matrices and vectors are represented by normal uppercase letters corresponding to the matrices' subscripts. For example, a subscript x may take values from 0 to X — 1. The discrete Fourier transform (DFT) of a matrix shall be represented by the script uppercase version of the same letter used to represent the corresponding space domain matrix. It shall be defined as follows: -, M - 1 N - 1 ' - = ^ E E F „ « - * « » * V ) (A.I) I V 1 1 > m = 0 n = 0 The corresponding inverse transform is: M - 1 N - 1 F m n = E E ^"W+f?) (A.2) u = 0 v=0 These definitions are taken from Rosenfeld and Kak [13], as are the following asso-ciated theorems: 36 A . l P A R S E V A L ' S T H E O R E M M - 1 N - 1 M - 1 N - 1 E E f » g : = M N E E (A.3) m=0 n-Q u=0 u=0 A .2 T H E C O N V O L U T I O N T H E O R E M (DFT{F * G } ) M = PQ7 p qQ p q (A.4) Here * denotes circular convolution, which is defined as follows: ! F * G U = ^ E F P , G m - P , „ - g (A.5) where it is understood that G„p + p >jQ + g = Gpq for all integers a, b; in other words, the matrices F and G are extended periodically as required by Equation A.5. This theorem states that circular convolution in the space domain is equivalent to point-by-point multiplication of the corresponding spectra. Throughout this paper, the circular convolution of any two matrices is taken to the same as their linear convolution, since the operands are assumed to have been padded with zeros to reach dimensions equal to those of the output matrix. This is necessary in any case for Fourier methods to be applied to image restoration problems, since all spectra must have the same dimensions if they are to be, for example, added or multiplied together. 37 A . 3 P O W E R S P E C T R A L D E N S I T Y I have chosen to define the power spectral density of a random field as follows: Sn{p,q)=E[\^\ (A.6) This includes the usual definition of 5 n as the discrete Fourier transform of the autocorrelation function of the random field: Sn{p,q) = [DFT{R n } ] p 9 (A.7) which is limited in its application to homogeneous random fields [13]. Equation A.6 is a far more useful definition for use with our discrete degradation model (Equa-tion 2.2), since the noise sources N i and N 2 are limited in space and so cannot be considered homogeneous. Therefore, Equation A.7 cannot be used for our problem. On the other hand, it is easily shown that (A.6) exists for many nonhomogenous random fields (in particular, any field consisting of uncorrelated random variables). Let us consider our special case of zero-mean white Gaussian noise. Take N x ; it consists of JK independent Gaussian noise elements of variance a 2 . The rest of the P x Q matrix is filled with zeros. Since the D F T is linear, we may write: K M = E Z U i M ( A - 8 ) where Mjj is the D F T of a matrix filled with zeros except for the single point Ni ( i , j ) . Now, note that = V(p, 9) (A.9) 38 so that E | M ( p , g ) | 5 E P - l Q - l i=0 j=0 / P - l Q - 1 \ / P - l Q - 1 L E W E E 4 ( M ) \^t=o j=o y \i=o j=o N1(i,i)ra PQ (A.10) (A .11) since the cross terms, being uncorrelated, have an expected value of zero. Now, ^[|Ni(»,j)| 2] o\ if (0 < i < J) AND (0 < j < K) 0 otherwise (A.12) Hence: and, similarly, (M + J - l ) ( N + K - l ) a 2 2 P2Q2 (A.13) (A.14) * * * 39 A P P E N D I X B - D E R I V A T I O N OF T H E C O N S T R A I N E D D E C O N V O L U T I O N F I L T E R We begin with the Fourier version of our degradation model (Equation 2.3). We omit the subscripts pq, since they are clear from the context: g = PQBJ + PQ.A/,J + .A/2 (B.l) From this, we immediately obtain the expression for the total noise power in the degraded image: E E \9-PQB J|2 = EE|PQM> + A/2|2 (B.2) p= 0 g= 0 p= 0 9 = 0 Our constraint will be that the total computed noise power in the degraded image, given an estimate J, shall be equal to its expected value. Although we can hardly expect this to be true at individual points, it is a reasonable expectation for the image spectrum as a whole. In any caser we write: p - i Q - i p= 0 9 = 0 2 A E E E + M2\2 (B.3) p= 0 9 = 0 P - l Q - l = E E £ [ ( P Q M / + M ) ( P Q A / i 7 + A/2)*] (B.4) p= 0 9 = 0 P - l Q - l = E E E [ P 2 Q 2 - V i A / ! * J J * + M2M2*l (B.5) p= 0 9 = 0 since the cross terms, being uncorrelated, have an expected value of 0. This yields our constraint: p - i Q - i P - 1 Q - 1 E E \9 - PQS J|2 = E E [P2QX I J|2 + sn2 p= 0 9 = 0 p= 0 9 = 0 (B.6) 40 Now we wish to minimize expression 3.7. First, we must convert it to the frequency domain: P - l Q - l n 9 P - l Q - l E E M l = E E M i p=0 9=0 p=0 9=0 (since L and F are REAL) p - i Q - l = P Q £ £ l p Q - C ^ l p 9 (by Parseval's theorem) p=0 9=0 P - l Q - l = P 3 Q 3 E E l ^ f (B.7) p=0 9=0 Now let us define: U = P 3 Q 3 ] T £ \£?\2 + ^ E E [ 1 5 - PQB J | 2 - P 2 Q 2 5 n J J | 2 - Sn2] (B.8) p=0 9=0 p=0 9=0 To minimize (3.7) under the constraint (B.6), we must find the 7 which minimizes U} Since U is not an analytic function of 7, we must minimize it separately with respect to the real and imaginary components of 7. Let: A = 5R{/}, C = 5{7} 7 = A + jC (B.9) Then: t/= Ep^EjTo1 1?3QZ\£(A+JC)\2 + \\9-PQB(A + JC)\2 - A P 2 Q 2 5 n i |A + i C | 2 - XSn2] (B.10) Now we differentiate this with respect to each APQ and Cpq, and set the derivatives equal to 0. However, we can immediately see that this results in the same two equations for all pq, so we can drop the sums and just minimize Upq, as follows: Let: u = Upq = P 3 Q 3 \£(A + jC) | 2 + A \Q - PQB {A + JC) f xFor an explanation of the method of Lagrange multipliers, please see eg. [15]. -\P2Q2Sni\A + j C \ 2 - \ S n 2 = P3Q3CiZ*{A2 + C2) + X \gg* - P Q £ 8 * ( A - jC) -PQg*B{A + jC) + P2Q2BB*{A2 + C2) -?2Q2Sni(A2 + C2)-Sn2] (B.ll) C\TJ — = 0 = 2 P 3 Q 3 i i r . 4 + A [ - P Q £ r - P Q £ * B oA +2?2Q2BB*A-2P2Q2SniA] = A [ 2 P 3 Q 3 £ i r + 2AP 2 Q 2 BB* - 2AP 2 Q 2 5 n i ] -\PQ{9B* + Q*B) (B.12) Since gB* + 9*8 = 2%{gB*}, and letting A = PQ/v: PQ3?{£B*} ^ P 2 Q 2 | £ | 2 + P 2 Q 2 | B | 2 - P 2 Q 2 5 n i Similarly, rj r j — = 0 = 2P 3Q 3£rC + X\PQjg8* -PQJ9*8 a A > +2P 2 Q 2 B8*C - 2 P 2 Q 2 5 n i C = C 2 P 3 Q 3 C L * + 2AP 2 Q 2 BB* - 2 A P 2 Q 2 S „ 1 ] (B.13) +\PQj{9B*-Q*B) (B.14) This time, we note that gB* - 9*8 = 2j$s{9B*}, and again let A = PQ/v: = > r = ?Q*{9B*} f B , W P 2 Q 2 | i : | 2 + P 2 Q 2 | B | 2 - P 2 Q 2 5 n i 1 ' ; Since 7 = A + jC, we immediately obtain: 7 = ¥ Q t 9 B * fB 161 ^ P 2 Q 2 | i : | 2 + P 2 Q 2 | S ! 2 - P 2 Q 2 S n i 1 ' ] 42 In Appendix A , we obtained the power spectral densities 5 n i and 5„ 2 for the particular case of white Gaussian noise: = g - ( M + j - y - i ) , } l A l t ) Plugging these quantities into Equations B.6 &; B.16, we obtain Equations 3.6 & 3.10: P - l Q - l P- l Q -l EH\9-PQBI\2 = EE p=0 g=0 p=0 9=0 v P 2 Q 2 | £ | 2 + P 2 Q 2 | S | 2 - J K c r 2 J K c r 2 | J | 2 + ( M + J - 1)(N + K - l ) g | P2Q2 43 A P P E N D I X C - D E R I V A T I O N O F T H E W E I G H T E D L E A S T S Q U A R E S F I L T E R We wish to find an estimate 7 such that the total calculated noise power in the degraded image (the "residual"), weighted by its expected value, is minimized. In Appendix B we obtained an expression for the expected value of the residual using our degradation model: E I P Q - A / X J + ^ I 2 ] = P 2 Q 2 5 n i | 7 | 2 + 5 n 2 (C.l) where, again, the subscripts pq are omitted since it is clear from the context that all operations are carried out point-by-point. Thus, we wish to find the estimate 7 which minimizes the following function: f E ' i W Q f l ' ( c . 2 ) p=0 q=0 P2Q2'S'ni 17\2 + Sn2 This is a sum of positive real numbers, so it is minimized by simply minimizing each term. We note that this is'not an analytic function of J , so we separate 7 into its real and imaginary components (as we did for the derivation of Constrained Deconvolution): 7 = A + jC (C.3) so that we shall now minimize the following function: U [ A C ) _ [ £ - P Q B U + J C)| 2 99* - PQ^B*(A - JC) - PQ9*8{A + JC) p2Q 2 S n i (^ + C2) + Sn 2 P2Q2BB*{A2 + C2) P2Q2SnL(A2 + C2) + Sn, 44 (C.5) which we do by setting its partial derivatives with respect to A and C equal to 0: d£- = Q [P 2Q 2S n i(-4 2 + C 2) + Sn2] [ 2 P 2 Q 2 B B M - P Q £ B * - PQ£*B] dA [ P 2 Q 2 S n i ( A 2 + C 2) + S, "2J [99* - PQ9B*(A - jC) - P Q £ * B ( A +JC)} [2P 2Q 2S n iA] [ P 2 Q 2 5 n i ( ^ 2 + C 2) + 5„ 2] 2 P 2 Q 2 B B * ( ^ 2 + C 2)2P 2Q 2.S n iyl [ P 2 Q 2 5 n i ( ^ + c2) + 5 n 2 ] 2 We can remove the denominator on the (reasonable) assumption that it is nonzero: = > 0 = yl{2P 2Q 2BB* [ p 2 Q 2 5 n i ( A 2 + C 2) + S„ 2] - 2 P 2 Q 2 5 n i | ^ - P Q B ( y l + jC)| 2} ' -PQ{9B* + 9*B) [ p 2 Q 2 5 n i ( ^ 2 + C 2) + Sn2 P Q » { £ B * } P 2 Q 2 S , , , | g - P Q B 7 T (C.7) (C.8) Similarly: dU dC = 0 = c [P2Q2Sni(A2 + C2) + Sn2] [2P 2Q 2BB*C + j P Q £ B * - /PQ£*B] [ P 2 Q 2 5 n i ( A 2 + C 2) + 5 n 2 ] 2 [99* ~ PQ9B*(A - jC) - ?Q9*B(A+jC)} [2P 2Q 25 r a iC] [ P 2 Q 2 5 n i ( A 2 + C 2) + 5„ 2] 2 P 2 Q 2 B B * ( ^ 2 + C 2 ) 2 P 2 Q 2 5 n i C [ P 2 Q 2 5 n i ( y l 2 + C 2) + 5 n 2 ] 2 C {2P 2Q 2BB* [ p 2 Q 2 5 n i ( A 2 + C 2) + Sn - 2 P 2 Q 2 5 n ] \9-?Q8{A + jC)\2} +jPQ{9B* - 9*B) [ P 2 Q 2 S n i ( A 2 + C2) + Sn p2Q21fl|2 P 2 Q 25 , , , l g-PQB7] 2 ^ |0 | P*Q*Sni\t\*+SHN (C.9) (CIO) (C.ll) Since 7 = A + jC, we immediately have: 7 = p2Q2IO|2 _ P 2 Q 2 g » , | g - P Q B ^ I 2 ^ 1 1 P 2 Q 2 S , M | / | 2 + S „ 2 (C.12) 45 C . l O B T A I N I N G T H E R O O T S It was found that iteration of Equation C.12 yields a good estimate of 7 only if 7 is initialized to 7 (the other starting points tried gave poor results). Obviously this is useless in practice. If Weighted Least Squares is to find any practical use, we must find a way to obtain the roots of Equation C.12 analytically. This is done as follows. The trick is to notice that any solution of Equation C.12 must have the same phase as p B * , since the denominator is always real. Thus, we may write: J — | J|gj'arg(5B*) (C.13) (C.14) \g-PQB7\2= ( | £ | - P Q | B | | / | ) 2 (C.15) We already know the phase of 7. To find its magnitude, we rewrite Equation C.12, using Equation C.15, and solve for \7\ as follows: PQISIIBI \7\ = 2 02iDi2 _ r 2Q 2Ml£!- pQ|B!m) ; ^ 1 1 PiQ*Sni\f\i + S,l2 (C.W) \7 P 2 Q 2 | B | 2 ( P 2 Q 2 S N I|J| 2 + S „ 2 ) -P 2QX ( |£ i 2 - 2PQ|B||/||$| + P2Q2|B|2|J|2); = PQ|5 | |B | (P 2 Q 2 5 N I|J| 2 + 5 n 2) P 4 Q 4 | B | 2 5 n i \ 7 f + P 2 Q 2 | B | 2 5 N 2|J| - P 2 Q 2 5„ , \ Q \ 2 \ 7 \ +2P3QX|S||£||JT - P4QX|B|2|J|3 = P 3Q 3| i9 | |B|5 N I|J | 2+PQS„ 2|5 | |B| (C.17) (C.18) 46 => o - P 3 Q X I 5 I | S | | J | 2 + P 2 Q 2 (sn2\B\2-sni\9\2)\f\ ~?QSn2\9\\B\ (C.19) To our delight, the cubic terms cancel! We may thus obtain the solution of Equa-tion C . 1 6 using the quadratic formula: ?2Q2 (Sni\9\2 - Sn2\8\2) \7\ = ± 2 P 3 Q 3 S N I | £ | | B | 0 > 4 Q 4 (Sn2\B\2 - S „ J £ | 2 ) 2 +4P4Q4S N I S N 2 |£| 2 |B 2 P 3 Q 3 S n i | £ | | B sni\9\2 - sn2\B\2 ± \f(sn2\B\2 + sni\9\*y 2YQSni\9\\B\ j 2Sni\9\* - 2 5 n J B l 2 } \ 2 P Q S n i | £ | | B | ' 2 P Q 5 n i | 5 | | B | J ( 151 SnM \ \FQ\B\' P Q S J S I j We can now reconstruct J from Equation C .14: (C.20) (C.21) (C.22) (C.23) \Q\QB* S„2\8\9B* PQ|8||S||B|> P Q 5 „ J 5 | 2 | S | T — / S -S,I2B* \ 7 \ PQS ' PQS„j£* / } (C.24) (C.25) These roots have been verified by substitution into Equation C .16 . The problem remains to select the "best" root from (C .25) . The use of the point form of the constraint: \9 -?QB7\Z = P2Q2Sni\T\2 + Sn2 (C.26) can be justified on the principle that it chooses the root for which the noise power is closest to its expected value. Of course, we have not yet established whether the two roots (C.25) are indeed minima of Equation C .2 ; this is discussed in Section 3.5.2 47 of t h i s paper. * * * 48
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Restoration of images degraded by systems of random...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Restoration of images degraded by systems of random impulse response Revelant, Ivan L. 1987
pdf
Page Metadata
Item Metadata
Title | Restoration of images degraded by systems of random impulse response |
Creator |
Revelant, Ivan L. |
Publisher | University of British Columbia |
Date Issued | 1987 |
Description | The problem of restoring an image distorted by a system consisting of a stochastic impulse response in conjuction with additive noise is investigated. The method of constrained least squares is extended to this problem, and leads to the development of a new technique based on the minimization of a weighted error function. Results obtained using the new method are compared with those obtained by constrained least squares, and by the Wiener filter and approximations thereof. It is found that the new technique, "Weighted Least Squares", gives superior results if the noise in the impulse response is comparable to or greater than the additive noise. |
Subject |
Image transmission Image processing -- Digital techniques |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-07-21 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0065618 |
URI | http://hdl.handle.net/2429/26731 |
Degree |
Master of Applied Science - MASc |
Program |
Electrical and Computer Engineering |
Affiliation |
Applied Science, Faculty of Electrical and Computer Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-UBC_1987_A7 R48_3.pdf [ 3.8MB ]
- Metadata
- JSON: 831-1.0065618.json
- JSON-LD: 831-1.0065618-ld.json
- RDF/XML (Pretty): 831-1.0065618-rdf.xml
- RDF/JSON: 831-1.0065618-rdf.json
- Turtle: 831-1.0065618-turtle.txt
- N-Triples: 831-1.0065618-rdf-ntriples.txt
- Original Record: 831-1.0065618-source.json
- Full Text
- 831-1.0065618-fulltext.txt
- Citation
- 831-1.0065618.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0065618/manifest