R E S T O R A T I O N OF R A N D O M L Y B L U R R E D IMAGES W I T H M E A S U R E M E N T E R R O R IN T H E POINT SPREAD F U N C T I O N Edward W . H . Lam B.Sc. Computer Engg., University of Alberta, 1984 A T H E S I S S U B M I T T E D IN P A R T I A L F U L F I L L M E N T O F T H E REQUIREMENTS FOR T H E DEGREE OF Master of Applied Science IN T H E F A C U L T Y OF G R A D U A T E STUDIES D E P A R T M E N T OF E L E C T R I C A L ENGINEERING We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y O F BRITISH C O L U M B I A August 1990 © Edward Lam, 1990 In presenting' this degree at the thesis in partial fulfilment of University of British Columbia, I agree freely available for reference copying of department publication this or of and study. thesis for scholarly by this his or her purposes representatives. Ksi^\4JLv} y The University of British Columbia Vancouver, Canada Date DE-6 (2/88) AlA^t 2S7 jlYlO that the may be It thesis for financial gain shall not B-itcb^caJL requirements I further agree permission. Department of the that Library by the understood be an advanced shall make it permission for extensive granted is for head that allowed without of my copying or my written Abstract The restoration of images degraded by a stochastic, time varying point spread function(PSF) is addressed. The object to be restored is assumed to remain fixed during the observation time. A sequence of observations of the unknown object is assumed available. The true value of the random P S F is not known. However, for each observation a "noisy" measurement of the random P S F at the time of observation is assumed available. Practical applications in which the P S F is time varying include situations in which the images are obtained through a nonhomogeneous medium such as water or the earth's atmosphere. Under such conditions, it is not possible to determine the P S F in advance, so attempts must be made to extract it from the degraded images themselves. A measurement of the P S F may be obtained by either isolating a naturally occurring point object in the scene, such as a reference star in optical astronomy, or by artificially installing an impulse light source in the scene. The noise in the measurements of point spread functions obtained in such a manner are particularly troublesome in cases when the light signals emitted by the point object are not very strong. In this thesis, we formulate a model for this restoration problem with P S F measurement error. A maximum likelihood filter and a Wiener filter are then developed for this model. Restorations are performed on simulated degraded images. Comparisons are made with standard filters of the classical restoration model(ignoring the P S F error), and also with results based on the averaged degraded image and averaged PSF's. Experimental results confirm that the filters we developed perform better than those based on averaging and than those ignoring the P S F measurement error. n Table of Contents Abstract ii Acknowledgement vii 1 Introduction 1 1.1 The Nature of the Restoration Problem 2 1.2 Restoration of Stochastic Distortion 4 1.3 Restoration Based on Image Ensembles 5 1.3.1 5 1.4 Ensemble Averaging 1.3.2 - Multi-channel Images Restoration 6 1.3.3 7 Heuristic Methods Objective and Scope 7 2 Point Spread Function Measurement Error Model 9 2.1 Frequency Domain Representation 11 2.2 Signals Statistics 12 2.3 Comparison with the Stochastic Distortion Model 15 3 Restoration Filters 16 3.1 Introduction 16 3.2 Restoration Based on the Averaged Image 17 3.3 The Regression Filter 18 3.4 The Maximum Likelihood Filter 20 iii 3.4.1 M L Estimate with H as a Fixed Vector 21 3.4.2 M L Estimate with H as a Random Process 25 3.5 The Wiener Filter 28 3.6 Introduction 28 3.6.1 Wiener Estimate with H as a Fixed Vector 30 3.6.2 Wiener Estimate with H as a Random Process 32 3.7 3.8 Discussion of Filters 35 3.7.1 Regression, M L and Wiener Filter(H constant) 36 3.7.2 Averaged Filters and Wiener(H random) 40 Filter Modification for Improved Stability 4 Experimental Results 42 44 4.1 Experimental Set Up 44 4.2 Discussion of Results 45 4.2.1 Adjusting for 111-Conditioning 46 4.2.2 Conventional vs. Modified filters 46 4.2.3 Comparison of Filters Based on Averaging 48 4.3 Modified Filters vs. Averaged Filters 49 4.4 Image Results 49 4.5 General Conclusions 49 5 Conclusion 5.1 80 Future Considerations 81 References 83 Appendices 87 iv A Derivation of Complex Regression Formula B Derivation of Complex Wiener Filter C Error of Least Square Estimator List of Figures 4.1 Comparison of niters adjusted and unadjusted for ill-conditioning 4.2 Comparison of conventional and modified filters(adjusted) 54 4.3 Comparison of averaged filters 58 4.4 Comparison of averaged and modified filters (adjusted) 62 4.5 Original and degraded images of "fighter" 66 4.6 Restorations of "fighter" at SNR=40dB 67 4.7 Restorations of "fighter" at SNR=35<£B 70 4.8 Original and degraded images of "girl" 73 4.9 Restorations of "girl" at S N R = 4 0 d £ 74 4.10 Restorations of "girl" at S N R = 3 5 d £ 77 vi . . . . 51 Acknowledgement I would like to thank my supervisor Dr. Rabab K . Ward for her guidance and encouragement throughout this research. I am also deeply indebted to her for her efforts in writing and encouraging me to publish the material of this thesis. I would also like to thank my colleague Frank Wan for his help and advice throughout this two year period, and Rob Ross, for his assistance in the usage of the computer facilities. vn Chapter 1 Introduction In a practical imaging system the achievable resolution of images are usually far poorer than that predicted by the theoretical diffraction limit. Factors responsible for the reduction in resolution include, inherent aberrations in the imaging system, inhomogeneities in the medium through which the image is received, and motion or misfocusing of the camera with respect to the scene objects. Resolution enhancement through means of deconvolution has found applications in a wide range of areas such as astronomy [1] [2] [3] [4], biomedical sonagraphy [6] [7] [8] and remote sensing [9]. The general area of image restoration deals with the problem of recovering an original image from its degraded and noisy observation(s). In this study, we address the restoration problem in which multiple observations of a fixed scene are available, and the degradation process is stochastic. Such a problem arises when imaging through nonhomogeneous media such as water or the earth's atmosphere. The success of the restoration depends on the accuracy of the knowledge of the point spread function(PSF) as well as the restoration algorithm itself. Unlike the case in which the P S F is constant, a system with a time varying P S F can not be calibrated once and for all before use. Although when a priori knowledge of the P S F is available and the P S F has a fairly simple mathematical form, it may be possible to estimate the P S F analytically, in many practical situations, the distortion is unknown or is too complex to be determined mathematically. In such cases, one may have to resort to blind deconvolution techniques [8] [14]. One potentially more accurate method of determining the P S F however, is to extract it from 1 Chapter 1. Introduction 2 the distorted images themselves. This is possible when the original image contains point objects(such as isolated stars in astronomical images). If the images are known to contain edges or lines, the P S F may be derived by utilizing the knowledge about these structures at different orientations [12]. Obtaining the P S F from the degraded image should give more accurate results than that using blind deconvolution since some information about the P S F is available. However when the P S F is determined from the degraded image, the accuracy is limited because the measurement is degraded by the same random noise as that affecting the rest of the image. This noise may cause significant error in the P S F since the objects from which the P S F are measured often do not have very strong intensity levels. This error in the measured P S F has either been ignored by researchers in the past or has been dealt with via averaging. In this study, we explore the potential benefits of other filters which take into account the error in the observed P S F . We shall refer to this restoration problem as the "error in P S F measurements" problem. In chapter 2, we develop the mathematical model of this problem. Solutions based on the maximum likelihood and the Wiener criteria are developed in chapter 3. In the remainder of this chapter, we give a brief review of related models and restoration methods. 1.1 The Nature of the Restoration Problem Although this problem is different from the classical restoration problem, it shares the general characteristics and limitations encountered in the classical model. The classical image restoration problem deals with the restoration of a single image [10] [11]. Errors in the P S F are generally ignored or assumed negligible. The problem can be expressed as: g(x,y) = h(x,y)* = JJ f(x,y) f(x-ct,y + v(x,y) - P)h(ot,P)dad8+ (1.1) v(x,y) Chapter 1. Introduction where 3 * denotes the convolution operation g(x,y) is the observed image f(x, y) is the original image h(x,y) is the point spread function v(x,y) is random noise In general the above problem is known to be an ill-posed problem [16]. A "perfect" solution to such problems is generally not possible due to the following limitations [11]: 1. The problem may be singular. If the image contains components that are orthogonal to the point spread function, the convolution of the two gives zero result. 2. The problem may be ill-conditioned. The restoration problem is one of differentiation since it involves finding the solution from its integral. Differentiation problems are known to be noise sensitive. A small error in the degraded image may correspond to a large change in the input. 3. The solution is not unique. In the presence of random noise, mathematically there is no unique association between the original and the degraded image. When the noise and the distortion are small, statistical estimation techniques may be applied to obtain the "best" result. When the noise or the distortion is large, the result often becomes unacceptable due to the ill-conditioned nature of the problem. In such cases "regularized" solutions must be sought. Such solutions usually represent a trade off between the fidelity of the observed data and a second criterion such as requiring the solution to be a smooth image. Many linear and nonlinear methods have been developed in the past for the conventional restoration problem [11] [15] including the Wiener filter and its derivatives, the Chapter 1. Introduction 4 constrained least square filter, the maximum entropy filter, the pseudo inverse filter and the maximum likelihood filter. 1.2 Restoration of Stochastic Distortion Restoration with stochastic P S F has been studied in the past. The model for such problems can be written as: g(x,y) = = (h(x,y) + n (x,y))* f(x,y) + n (x,y) 1 2 (1.2) H ,V) * f( ,y) + ni(x,y) * f(x,y) + n (x,y) X x 2 where h(x,y) and ni(x,y) represents the deterministic and random part of the P S F respectively, and n,2(x,y) is the random detection noise. Denoting the noise terms ni(x,y) * f(x,y) + n (x,y) as n(x,y), equation (1.2) becomes: 2 g(x, y) = h(x, y) * f(x, y) + n(x, y) (1.3) where n(x, y) is a signal dependent noise. Methods of dealing with such models have been developed in the past [17] [20] [22]. Solutions based on the LMMSE(Wiener) and minimum variance unbiased(MVU) estimate have been derived in [17]. A restoration technique based on the Backus-Gilbert technique has also been developed in [20]. Several techniques including the Wiener filter, the constrained least square filter, the geometric mean filter, the nonlinear maximum a posterior(MAF) have been developed in [21]. Restoration based on the Projection into Convex Sets method has been developed in [24]. Although these models take into account the stochastic nature of the P S F , they do not account for the measurement error in the PSF. The difference between these models and the P S F error model will be discussed in the next chapter. Chapter 1. Introduction 5 1.3 Restoration Based on Image Ensembles The problem of restoring an image from multiple observations has been addressed in numerous previous studies for both deterministic and stochastic P S F models. A l though the knowledge of the P S F or its statistical properties have been assumed known in these studies, the P S F measurement errors have not been taken into account. In the following section, we outline some of the previous work on restoration based on multiple observations. 1.3.1 Ensemble Averaging With multiple observations available, researchers are faced with the task of finding a method which can take advantage of the information in the large amount data while keeping the storage and computational requirements at a manageable level. The implementation of the "optimal" solution often requires a prohibitively large amount of computational resources. Methods in which the data is compressed before the restoration is performed are often more practical even though an inferior estimate is expected. The most straight forward compression method is the averaging of the images. In application to images distorted by atmospheric turbulence, multiple short exposure images of a single original image are often obtained. The observations are often combined by direct averaging, and restoration is then based on the averaged P S F and the averaged degraded image. Labeyrie [5] has suggested the method of speckle interferometry. With this method, instead of averaging the observations directly, averaging is performed on the autocorrelation of the images. This way, the square magnitude spectrum of the P S F and the degraded image are obtained. The magnitude spectrum of the original image can then be recovered from the results. The advantage of this method is that while the averaged spectrum of the P S F goes to zero at a frequency far below the diffraction Chapter 1. Introduction 6 limit, the averaged square magnitude spectrum retains its accuracy up to much higher frequencies. The drawback of the method is that only the magnitude of the original image is recovered, and the phase information must be retrieved separately. 1.3.2 Multi-channel Images Restoration Restoration methods which utilize all available data have been developed in the past [18] [22] [26] [27] [28]. Multichannel restoration techniques take advantage of the correlation of the signal (or noise) amongst channels, giving an improved estimate compared to restoring each channel individually. Due to the large size of the problem, solving the problem in conventional manners is infeasible in most practical applications. Galatsanos and Chin [27] addressed the multichannel restoration problem using the Wiener(LMMSE) estimate. Although their model is developed for restoring multiple original images, the algorithm can be adopted to the single image case. Guan and Ward [22] applied the L M M S E to restore multiple observation of randomly blurred images. In their model a single original image is being restored, the noise is correlated with the input as a result of the random nature of the point spread function. In both of the above studies, the problem with large size of the data due to large matrices is managed by transforming the problem into the frequency domain. Another approach for reducing the size of the multichannel problem makes use of Karhunen-Loeve type transforms [18] [22] [26]. B y applying such a transform, each of the channels is de-correlated from the rest. Thus the problem is reduced to that of the size of a single channel case. This approach assumes that the problem is time-space separable, i.e. the correlation matrix of the noise (or signal) can be expressed as a product of a time function and a space function. In the above studies the noise amongst the different observations is assumed correlated. If this correlation amongst the different frames is absent, these methods reduce to Chapter 1. Introduction 7 the Wiener filter applied on the averaged observations. 1.3.3 Heuristic Methods The amount of data is expected to remain very large even after its reduction by the techniques described above. For large image ensembles, the computations still remain excessive. Thus heuristic techniques are sought in many practical applications. Ward [19] restored stochastically blurred ensembles by first deblurring individual images through the Backus-Gilbert technique, then combining the result with the maximum likelihood estimator at each pixel in the time domain. If the noise correlation amongst the different frames is not present, this method reduces to averaging the individually deblurred images. Ghiglia [25] deblurred multiple, independently blurred images using the constrained least square approach. His model assumes uncorrelated noise. The filter is formulated in the frequency domain and has great similarity with the ones used in this study. His study shows that good results are obtainable by using a few images only. 1.4 Objective and Scope In this study, we address the problem of restoring an original image from multiple degraded observations. The exact P S F for each image is not known, but a noisy measurement of it is available from each image. A model for this problem is developed in chapter 2. In chapter 3, restoration filters are derived based on the maximum likelihood technique and the Wiener criterion. The results are compared to the regression filter and the conventional Wiener filter derived for the conventional model(ignoring P S F error). In chapter 4, the performance of the filters are tested with simulated degraded images. Comparisons are made between the filters for the " P S F measurement model" and the conventional model. Comparison are also performed with results obtained based Chapter 1. Introduction on averaging of the multiple observations. Chapter 2 Point Spread Function Measurement Error Model In this study, we deal with the restoration problem in which multiple observations of a fixed scene are available. The distortion is assumed to be time variant. The images measured are different because each observation is distorted by a different distortion pi-ocess. For each image, a measurement of the P S F is assumed available because of the assumed presence of a point object in the scene. In this chapter we formulate a degradation model for this problem. The model will take into account the measurement error in the P S F , which by large has been neglected in the conventional dealings of this problem. The statistical properties of the noise will be presented. A brief comparison to some of the existing restoration models will then be given. We assume that the digitized degraded images are sampled on a regularly spaced grid, and the pixels are perfectly registered. In other words, the same pixel on each of the images corresponds to the same point in the scene. We also restrict our attention to space invariant distortions only. Under the above assumptions, we formulate the problem directly in the discrete space time domain. For an input of size M by M and an ensemble of size T, we have: g(x,y,t) = H .Y2 h(x-x',y-y',t)f(x',y') x = yl h(x,y,t)*f(x,y) + u(x,y,t) + u(x,y,t) (2.1) p(x,y,t) = h(x,y,t) + v(x,y,t) x,y = L . M , t = 1..T 9 « Chapter 2. Point Spread Function Measurement Error Model where 10 * denotes the convolution operation g(x, y, t) is the measured degraded images at time t p(x,y,t) is the measured point spread function at time t h(x,y,t) is the true unknown point spread function at time t f(x,y) is the original image to be estimated u(x,y,t) is the detection noise in g(x,y) at time t v(x,y,t) is the detection noise in p(x,y) at time t When the P S F is obtained from a point object in the degraded image, and the intensity of the point object(/ ) is known, the computation of p(x,y,t) can be described as in the 0 following. Denoting the unnormalized P S F as p'(x,y,t), we can write: p'(x, y, t) = g(x + x,y + y) 0 0 - ! < * < ! , (2.2) - | < y < ! where x , y is the coordinate of the point object and S is the width of the support area 0 Q given to the P S F . Substituting from (2.1), we get, p'(x, y, t) = I h(x, y, t) + v'(x, y, t) 0 - ! < * < ! > (2.3) - f < y < f where v'(x,y,t) = u(x + x,y + y,t) and I is the intensity of the point object in the 0 0 0 original image. The P S F p(x,y,t) is then, p(x,y,t) yp'(x,y,t) h(x,y,t) + v(x,y,t) where v(x,y,t) yv'(x,y,t) (2.4) Chapter 2. Point Spread Function Measurement Error Model 11 When I is unknown, its value can be approximated by: 0 (2.5) &'(3,y,<) It can be seen from (2.4) that since Yl h(x, y, t) is unity, and the noise v'(x,y,t) has xy zero mean. Averaging (2.3), we obtain (2-6) 2.1 Frequency Domain Representation Transforming equations (2.1) into the frequency domain, g(u ,uj ,t) = h(u ,u ,t)f(u ,u ) p(u) ,oj ,t) = h(u ,u> ,t) + h(w ,u ,t), f(u ,u ), x y x where g(u ,u ,t), x y y x y x y x x x the Fourier Transforms for g(x,y,t), + u(uj ,u; ,t) y x (2.7) y v(u ,u ,t) y x (2.8) y u(w ,u> ,t), p(u} ,u> ,t), v(u> ,u ,t) represents y x y h(x,y,t), f(x,y), x y u(x,y,t), x y p(x,y,t) and v(x,y,t) respectively. A t each frequency pair u) ,u> the vectors G, P , H, V, U are defined as: x y p(u ,u ,l) h(u ,u , p(u ,u ,2) h(u>x,Uy, 2) x g(u ,u ,2) x y x x y y G = y 1) (2.9) H = g(u ,u ,T) x y p(u ,u ,T) x y h(u ,u; ,T) x y Chapter 2. Point Spread Function Measurement Error Model u(w*,w ,l) v(u> Uy,l) u(u) ,u ,2) v(u ,u ,2) Si y x y x y (2.10) V = u = v(u ,u ,T) u(u> ,u ,T) x 12 x y y so that, G = P = Hf+U H +V (2-11) where / is a scalar quantity. The first equation in (2.11) has a form similar to the classical linear regression problem: Y = X8+v (2.12) where X is the vector of independent variable, 8 is the unknown parameter to be estimated, and v is the noise vector. The major differences between our model and the linear regression model are that in our model, 1. the independent variables H are corrupted by the noise V, 2. all variables are complex variables 2.2 Signals Statistics We shall consider the general case where the unknown vector of PSF's may be fixed or random. In our model, if the P S F is random then the vector of the transfer functions H in general is composed of a deterministic and a random part. H = Hd + H r (2.13) Chapter 2. Point Spread Function Measurement Error Model 13 where E[H] = H (2.14) A The noise u(x,y,t) is assumed to be zero mean, Gaussian distributed, and uncorrelated in space and time. The variance of u(x,y,t) is denoted by <j\ and is assumed known. This is a popular assumption since in practice, this statistic may be measured with a constant intensity reference image. In the frequency domain, the variance of the noise is found to be constant over all frequencies with a value of o\. And since u{x, y, t) is independent for each of the images, E[u(u ,uj ,i)*u(u ,u ,j)} x y x y = o\ for i = j = for i 7^ j (2.15) 0 The noise v{x,y,t) is proportional to a windowed version of u(x,y,t). As a matter of fact, both noise sources v(x,y,t) and u(x,y,t) are windowed versions of an infinite noise source n(x,y,t). The size of the window for v(x,y,t) is S by S and that for u(x,y,t) is M by M. Denoting the window function by w(x,y), v(x,y,t) = n(x,y,t)w(x,y) The autocorrelation of v(x,y,t) at time t, R (Ax, Ay, t), is: vv R (Ax,Ay,t) (2.16) = E [v(x,y,t) v(x + Ax,y -f Ay,t)} vv = E [n(x, y, t) n(x + Ax, y + Ay, t) w(x, y) w(x + Ax, y + A?/)] = Rn {Ax, n Ay, t) S A (Ax, Ay) 2 s (2.17) where As denotes a triangular function of width 25 by 2S with a maximum height of 1. Similarly R (Ax, uu Ay,t) may be found to be I^R (Ax, nn Ay, t) M K (Ax, 2 M Ay) (2.18) Chapter 2. Point Spread Function Measurement Error Model 14 where AM has a width of 2M by 2M with a maximum height of 1. Thus, R (Ax, y + Ay, t) = vv In the case when R (Ax, f I'M ^ [ ' R (Ax, h {Ax,Ay) S Ay, t) Ax Ay) uu 1 (2.19) M Ay, t) = o\, 6(0,0,/) where 8 is the Dirac delta function, we uu get: R (Ax, Ay, t) = | ^ ( 0 , 0 , t) vv Thus R (Ax, Ay, t) is proportional to a delta function. The power spectrum of v(x, y, t), vv being the Fourier transform of R (Ax, S (u> ,Uy), vv Ay, t), is a constant function. It can vv x (2.20) be concluded then that S 2 S (u ,u ) vv where S x y = pj^vl is the support area of the P S F and M 2 (2.21) = °l is the size of the image. As with 2 u(x,y,t), v(x,y,t) is uncorrelated from one image to another, E[v(u ,u ,iYv(u ,u ,j)} x y x = y c 2 = 0 for i = j for i^ j (-,.22) The cross correlation of u and v can be shown to be E[u(x,y,t)v(x + x',y + y',t + t')] = I a 2 0 = where x ,y 0 0 0 v x' = x ,y' = y ,t' = 0 0 0 otherwise is the coordinate of the point object. If u(x,y,t) lies outside the region of interest (the region being restored) in the image, then the correlation is zero, i.e. E [u{x,y,t)v(x\y',t')] =0 V x,y,t (2.24) If the point object lies inside the region of interest, this correlation must be taken into consideration. In the derivation of the filters in the next sections, for simplicity, we assume the point object to be outside the region being restored. Chapter 2. Point Spread Function Measurement Error Model 2.3 15 Comparison with the Stochastic Distortion Model Models used in previous studies for stochastic distortions [19] [22] have similar form to the one presented above. In these previous models, where g{x,y,t) = h(x, y, t) * f(x, y) + n (x, y,t) h(x,y,t) = h(x,y) + n (x,y,t) a (2.25) 2 * is the convolution operator h(x, y) is the known average or expected value of the P S F ni(x,y,t) is the additive noise n (x,y,t) is the random component in the P S F 2 Combining the two equations of (2.25), we obtain g{x,y,t) = h(x,y) * f(x,y) + n (x,y,t) 2 * f(x,y) + rn(x,y,t) (2.26) For our problem, (2.27) P{x,y,t) = h(x,y,t) + v(x,y,t) Substituting the second equation into the first, we get g( , y, *) = P{x, V, t) * f( , y) - v(x, y, t) * f(x, y) + u(x,y, t) x x (2.28) The differences between our P S F error model(2.28) and previous model(2.26) are that 1. in our model the independent variable p(x,y,t) is correlated with the noise term since it is correlated with v(x,y,t) whereas in the stochastic distortion model, h is not correlated with the noise term 2. p(x,y,t) is a function of time whereas h(x,y) is constant. Chapter 3 Restoration Filters 3.1 Introduction In order to compare the performance of the " P S F measurement error model" to that of the regular restoration model, in this chapter, we derive two niters. The first one will be the maximum likelihood(ML) filter. Its results will be compared to the conventional regression(minimum least square) filter. The second will be a LMMSE(Wiener) filter and its results will be compared to that of the conventional Wiener filter. Two different versions of each filter will be examined. The first one treats the values of the P S F as fixed numbers; the second as random drawings from random variable. For comparison with restoration based on the ensemble averages, we also define the inverse and Wiener filter for the averaged images. The filters will be developed and applied in the frequency domain on a pixel by pixel basis. One major advantage of all of these filters aside from their simplicity, is that they utilize only the first and second order statistics of. the images which can be updated as each image is received. Thus the actual image ensemble does not have to be retained, allowing arbitrary large ensembles to be processed. After the filter derivations are presented, a brief discussion will be given on their expected performance. 16 Chapter 3. Restoration Filters 3.2 17 Restoration Based on the Averaged Image Averaging is the most straightforward method for combining the image ensemble. The averaged degraded image and P S F can be defined as: 1 T t 1 1 ™£p(z,<M) t T p{x,y) = 1 (3-1) From chapter two equation (2.9), at each spatial frequency, G = Hf + U P = H +V (3.2) Denoting the mean of the vectors G, P, H, U, V as g, p, h, u, and v, we can write p = h+v g = hf + u = (p-v)f (3.3) +u Since the variances of the noise u and v are reduced by a factor of T by the effect of averaging, we have: 2 <4 = 1 2 \°l (3.4) Any conventional filter developed for single image restoration may be applied to the averaged image. For comparison with the regression and Wiener filter to be developed in the next section, we define the inverse and the Wiener filter. The inverse filter of the averaged image can be defined at each frequency as: Chapter 3. Restoration Filters 18 The conventional Wiener filter for single image restoration [11] can be defined at each spatial frequency as: ~ W ! + kls„ ( 3 ' 6 ) where cr is the variance of the noise and SJJ is the power spectrum of the input / . In 2 our case, h is replaced by p, so the Wiener filter for the averaged image is: 1 = WToWsj; ( 3 - 7 ) Since SJJ is usually unknown in real life, an iterative approach is adopted as in numerous previous studies [17] [21]. Utilizing the averaged power spectrum of the degraded images as an initial estimate for Sfj, an estimate / can be calculated. A new S/j can then be found from the power spectrum of / as follows. U (3-8) = \P? + °l/Sf fi where 1 7j;zZ 93 t |%| T S ffi = = 3.3 iori = 0 S 1 for^O 2 (3.9) The Regression Filter As shown in chapter two equation (2.11), the errors in measurement model in the frequency has the form of a typical linear regression problem with error in the independent variable. The difference is that in our case all variables are complex variables. The errors in independent variable model has long been a topic of interest to statisticians and econometricians [29] [30]. The usual regression model can be written as: 9 i = hif + u, i = l..T (3.10) Chapter 3. Restoration Filters 19 where T is the number of observations, / is the parameter to be estimated, h and g are the independent and observed variables, and u is measurement noise. The ordinary regression model finds the estimate which minimizes the mean square of the residual errors. The residual error is the difference between the estimated (jr, and the observed g;. Thus the problem is: T mm^2\ -hj\ / ,=i 2 gi (3.11) The solution is shown to have the form(Appendix A ) : / LiM^. = (3.12) where * denotes the complex conjugate. This estimate is the best linear unbiased estimator(BLUE) [33] under the conditions that: 1. the noise is statistically independent of the variable h 2. the noise has an expected value of zero, and is serially independent with a constant variance For the errors in variables model J Pi = (3.13) hi + Vi where p,- is the observed independent variable and u,- is the corresponding measurement noise, condition 1 is not met and the estimate f = Ikl& (3.14) EiP*Pi is not the B L U E estimator. This can be shown by rewriting (3.13) as: 9i = Pif + "t Hi = Ui - Vif (3.15) Chapter 3. Restoration Filters 20 Although u and v are mutually and serially uncorrelated, the covariance of p and n is E\p*n] - E\p*]E[n] = E[(h*+ v')(-vf+ u)] - E[h']E[-vf+ u] (3.16) Since this covariance is nonzero, a dependency exists between the noise n and the explanatory variable p in (3.15). The consequence of this dependency is that the straightforward application of linear regression would yield a "biased" estimate of / . Furthermore this bias will not disappear as the sample size becomes infinitely large. That is, the estimate (3.14) is "inconsistent". Since the noise sources are distributed independently of one another and of the true values of the independent variable, it can be shown that the limiting value to which / in (3.14) tends in probability is 2 plimf = -£—f °~h Thus plim f (3.17) + °v f but is in fact an underestimate of / . If the variance of v is say, 10 percent of the variance of the true transfer function h, then straightforward linear regression (3.14) would underestimate / by 10 percent even for large sample sizes. The ratio (3.18) is known in social science literature as the "reliability ratio" [30]. In the following sections we develop the maximum likelihood estimator and the Wiener filter for this problem. 3.4 The Maximum Likelihood Filter The maximum likelihood(ML) estimate is a useful and powerful nonlinear estimate. It gives an estimate of the unknown input / which maximizes the conditional density Chapter 3. Restoration Filters 21 function of G given / , i.e. max (3.19) P{G\f) where G is the observed vector. This estimate is identical to the maximum a posterion'(MAP) estimate when the the input distribution of / is constant. The M A P estimate chooses a value of / which maximizes the a posteriori density function of / conditional on G, max (3.20) P(f\G) which is the same as max P(G\f)P(f) (3.21) When P(f) is unknown(as in our case), we have to be content with maximizing P{G\f) which is the maximum likelihood method. When G is a joint Gaussian function, the M L is known to be a linear estimate. There are two cases to be considered when finding the M L estimate to our error in independent variables model. The first one regards the true independent variables /i;, i = 1..T as constants; the second treats the h\s as samples of a random variable. In the following sections the two approaches will be discussed. 3.4.1 M L Estimate with H as a Fixed Vector With Ui and v, being independent zero mean Gaussian random variables, the joint density function for the observed vectors G and P will have the following multivariate normal distribution: P(G ,G ..G , 1 2 P ,P ..P \f, T 1 1 —=——exp(V2iro ) T v 2 h h ..h ) = T u 2 T £ \ Pi - hi\ /2o 2 T 1 T 2 v -exp-^gi-hif] ^ 2 (\/2^ ) u (3.22) Chapter 3. Restoration Filters 22 The corresponding log likelihood function is then: C = - |/n(2™*) - ~ £ \ - A,f Pi ± 52\ -hif\ (3.23) 2 gi Taking partial derivatives with respect to / * and equating to zero, dC T (3.24) Denoting the M L of a variable / by / , (3.25) The above result indicates that the M L estimate for the error model is the same as the linear regression with hi, the true unknown variable, being the independent variable. The M L estimate hi of the variable hi can be found by taking the derivative of (3.23) and equating the result to zero. ^ = ^(Pi-h ) oh* °l + ^(g -fh ) i i < (i = i..r) =0 i (3.26) Rearranging the above equation and letting (3.27) we get hi = Pi + Xf*9i 1 + A|/P (3.28) Substituting the above into (3.25), / = Zplgi + VEgtgi Ep-Pi + A / * + XfE g* (1 + A | / | ) Pi l + A/2 2 Letting 1 9'9 M = T jiY,9i9i 2 + A |/| E sta 2 2 -i (3.29) Chapter 3. Restoration Filters 23 1 T 9'V = 1 y & i W M* = TfEptlH T M P P (3.30) Equation (3.29) can be written as: i (l + \\f\ )(M . + 2 = p g M. p p + A/*M. + \fM . P S g p \fM . ) g g + \ \f\ M . 2 2 g g which simplifies to the quadratic: \pM . g p - f(\M . g g - M.) - M, p p p g = 0 (3.32) with solutions, h = 0 + yJP + \ f = 0-y/F+\ (3-33) 2 M -„ — \M * Q V V The value of 6 can be computed from <r , o~ and the observed vectors G and P. 2 2 However / has two possible values and the problem as to which / is the correct answer arises. For the case when / is a real(not complex) number, it is clear from inspection of equation (3.33) that irrespective of whether 8 is positive or negative, f is always positive x and fi is always negative. Kendall [31] has advised taking the positive root to maximize the likelihood function. Johnston [29] suggested choosing the / depending on the sign of the relation between hif (the true value of gi), and /i;(the true value of pi). If pi and gi are both of the same sign, we expect / to be positive, and if they are of opposite signs, we expect / to be negative. Thus the advice is to choose the positive / if M pg (= -f J2Pi9i) i positive and to choose the negative / if M s pg is negative. The last advice seems intuitively correct. For our case the variables are complex numbers. To determine Chapter 3. Restoration Filters whether fx or / 2 24 should be chosen, we calculate M * p g (3.35) and assuming the last three terms are very small, we get: (3.36) And since J2i I hi | is a real number, the phase of M * is approximately equal to the phase 2 p g of / . Therefore we choose the / whose phase is closer in value to that of M * . However p g in the experimental results we did not have any luck with this method. The results were incomprehensible. In fact the restored image consists of noise only. To understand why, we had another look at the equations and a deeper look at the characteristics of our variables. It is found that the value of A in our application is a small value. (3.37) where S is the support area of the P S F , M 2 2 is the size of the image, and I is the 0 gray level of the original point object from which the P S F is determined. With A being small, the coefficient of the quadratic term in equation (3.32) (XM * ) becomes small in p g comparison to the other coefficients in the equation, and the estimator has ill-conditioned characteristics. This can be demonstrated by considering the extreme case when A = 0. When the solution is calculated from equation (3.32) we get /= (3.38) which is the correct answer. However calculations done using the latter(3.33) give 0 = —co and fi,h = —00,4-00. When the value of A is small(but nonzero), the result is finite, but a small error in the value of 6 causes a large change in the solutions. Thus the solution of (3.33) becomes very unstable. Because of this unstable characteristic of Chapter 3. Restoration Filters 25 this estimator, we abandoned this approach and looked for an approximation of / from equation (3.25) that is less noise sensitive. Rewriting (3.25) as: i _ TJ(P* - Vj)9i £, (Pi ~ v*i)(pi - ,„„ * Q Vi) in the limit as the number of observations increases, the M L estimate can be approximated by: 3.4.2 M L Estimate with H as a Random Process In this section, we derive a M L estimate for / by treating the true independent variables hi as random samples drawn from a fixed Gaussian distribution. Rewriting every complex variable as the sum of its real and imaginary part, the system P = H +V G = Hf + U (3.41) becomes Pr+jPi G +jGi = = r H +jHi-rV +jVi r (H +jHi)(f +jfi) T (3.42) r r + U +jU T t Separating the real and imaginary parts, we get Pr = H + V Pi = H + Vi G = Hf Gi = H fi + Hif + Ui r r r t r r r (3.43) — Hifi + U r r Assuming H , Hi, V , Vi, U , and Ui have independent normal distributions with variances r r r a\ , a\., a^ , o-\., al , a ., P , Pi, G , Gi, being linear functions of normal variables, will 2 r r r r T Chapter 3. Restoration Filters 26 have a multivariate normal distribution. Denoting the correlation of two variables x and y (E[xy]) as S , the following can be written for our problem: xy Sh h = ^PrPr r ~T r ^V V r SpiPi = Shihi + SgrPr ~ <^ffrP. = T - r = = frShrhi frSh h r r S Vi Vi ( -44) 3 ~ fi^hrhi r fr$h hi SgiPr SgiPi frSh h r fiShihi + fiSh hr r + fiShrhi Combining the above equations, we obtain ^9rPr ^flrP, fr{SprPr ~ ^ V = r ~ V ) r fi^PrPi frSprPi ~ fii^PiPi ~ SviVi) = (3.45) ^9iPr SgiPi For ~ frS pi 4" fi{Sp Pr Pr T Sv V ) r r + fiSp fr{$piPi ~ $ViVi) = TP i a multivariate normal distribution, it is known that the sample variances and co- variances are the maximum likelihood estimates of the variances and covariances, and the sample means are the maximum likelihood estimates of the means. Denoting the sample moment of any two variables x and y as M , the following can be written for xy our variables: Mp p = f 5Z{ Pr Pr T r M t t = f J2t PitPit pipi Mp pi — f Y2t PrtPit r = iZl9rtPrt M 9rPr Mg pi = f S ( 9r Pit r t Mp r = f X^t 9itPrt Mg iVl = hY7t 9hPit gi (3-46) Chapter 3. Restoration Filters 27 Since P , Pi, G>, and G , have a multivariate normal distribution, the sample moment r of any two variables is the maximum likelihood estimate of the correlation of these two variables. Substituting the sample moments for the corresponding correlation in (3.45), we obtain M PrVr - -M„. VrVi K vr M,PrVi PiPi rPr f -Af„..„.. +O M,VrPi M Mg l r Mp p r - oi . T 2 °~v V r r fr M . fi . M grVl gtPr _M PrPi iV (3.47) gtPi _ Since the estimate for a variable derived from a function of maximum likelihood estimates is itself a maximum likelihood estimate [33], we can solve for f ,fi r the above set of linear equations. by solving Adding the first and last row and subtracting the second from the third, 0 M ' fr ' _ 0 where + grPr . Mg ipr v = M PrPr + M M giPi ~ M , grp - o v Vr ' fr' . fi . Mg +Mg rpr Mg -Mg iPr _ -o 22 PiPi (3.48) 2 ViVi (3.49) ipi (3.50) rPi Therefore, Ht{Pt9t) Et M 2 - Tol (3.51) Equation (3.51) is the maximum likelihood estimate of / when h is considered random. It can be easily shown that this estimate is consistent, i.e., plim f = f. Notice that this is the same as the approximation(3.40) obtained in the previous section when the h'.s are considered fixed and the number of observations approaches infinity. Chapter 3. Restoration Filters 3.5 The Wiener Filter 3.6 Introduction 28 The Wiener filter has been a popular tool in image restoration due to its simplicity and effectiveness. In this chapter, the Wiener filter for the P S F error model will be presented. Unlike the regression or M L methods, the Wiener filter requires the knowledge of the autocorrelation function of the true unknown / . In general, given T observations g\(x), g (x).. .gr(x) of an original signal f(x), 2 the estimate which minimizes the mean square error between the original and estimated signal (3.52) is known to be [15] [32] the conditional mean of f(x) given <7,(x),i.e. f{x) = E [f(x)\ (x),g (x)..g (x)] 9l 2 T Vx (3.53) When the conditional densities are not known, such an estimate can not be easily obtained. In such cases, the estimate is restricted to be a linear function of the observed signal. The Wiener filter yields the estimator which is optimal amongst the class of all linear estimators which minimizes the least square errors as described in (3.52). It is also known that when the joint distributions of / and gi i = 1..T is Gaussian, the result of this linear estimate is the same as that obtainable from maximizing the conditional densities. In the conventional problem, the observations available are only those oi gi, i = 1..T. The estimate is taken as T (3.54) Chapter 3. Restoration Filters 29 and the values of the filter coefficients L\s are sought. In our problem, the observations are not only the g^s, but also the p^s. Thus in the most general case, our estimate / should be a linear function F of a combination of g[s and p\s, i.e., T f = zZ i i( ) L t x * F(gi,g .-gT,Pi,P2~PT)] (3.55) 2 This is no longer a straight forward problem and the uniqueness of the solution is no longer guaranteed. For our purpose then, we shall follow the procedure as in the conventional case. We shall find the linear estimate based on the g' s as in (3.54) which minimizes the { square error (3.52). The answer then will be as in the conventional function of Sjg^s and Sgg^s. But unlike the conventional case, for our problem these values are not known since the system(the /it's) are not known. Thus for the values of the unknowns we shall substitute estimates based on the observations and statistics of the p;'s. From Parseval's energy theorem it can be shown that minimizing the mean square error in the frequency domain: minE[\f(u)-f(u)\ ) /(<") 2 (3.56) mmE[\f(x)-f(x)f] /(*) (3.57) is equivalent to minimizing in the space domain. As in the M L and regression filters, the Wiener will be applied in the frequency domain independently at each spatial frequency. A t each frequency u, we seek f(co) such that (3.56) is satisfied. Assuming /(w) = Ylf Li(u)gi(uj) at each frequency, the filter minimizes the mean square error min E L T EWftH-/(«)l a • (3.58) Chapter 3. Restoration Filters 30 It can be shown (Appendix B) that the minimization of the above leads to the following solution: ' S (l,2) gg S (2,l) S (2,2) gg gg (3.59) S (T,T) . S, {T) 3a where S (i,j) gg f = E[g*gj] and S j(i) = E[g*f] and £, i — 1..T are the coefficients of the g Wiener filter to be determined. The above equation is the general Wiener solution applied to any model. In our application, the model is, 9i = Pi = hif + Ui (3.60) hi + Vi As with the M L estimate, we consider two approaches, the first treats hi, i — 1..T as fixed values, the second treats the h'-s as samples from a random variable. The derivation for both approaches is presented in the following. 3.6.1 Wiener Estimate with H as a Fixed Vector W i t h u being independent of / and h and self uncorrelated in time, the autocorrelation of the observed vector G can be written as: S {i,j) gB = E[g^ ] = E[{h*r + u )(hjf + uj)] = Eih'hjSjA+EMxij] gj m (3.61) Chapter 3. Restoration Filters 31 Treating the elements in H as constants, Sggihj) = * jSff h + l h f a o r *= 3 (3.62) and S j{i) = a E[9if] = E[(h*f + u*)f} = (3.63) KSss For the special case in which the noise u, is uncorrelated in time, substituting (3.62) and (3.63) into (3.59) we get: h\h Sjf 2 h^hiSfj hfaSff + crl U KSff L h^Sff 2 (3.64) • h* h S T T ff + a hySff 2 u Each row i of the above matrix can be expressed as: ZfiVhjSff + alS^Lj ZJihlhjLjSfri + alLi = h*S = h*S ff (3.65) ff Solving for Li in the above equation, Li = ^ ( l - E M o - ) u = (3.66) 3 (3.67) h*a Substituting back into (3.65), h*iSff Sff (3.68) Zjh.Sffih^ + aKa) Chapter 3. Restoration Filters 32 Solving for a (3.69) 5 5 From (3.67) and (3.69) the filter coefficients are then, U = —~ (3.70) Since in reality the value of h* and E j hj are unknown, we use the sample estimates for these variables. The final Wiener estimate is then, E / \Pi\ - Tal + £ 2 ff For conventional mo dels (ignoring u,-), the resulting filter would have been / = ' P E f N 2 E T ^ , (3.72) + ^ Since 5 / / is unknown, the solution can be found iteratively as in (3.8), by using the value of 3.6.2 YA S gg as a starting estimate for S/f. Wiener Estimate with H as a Random Process In this case, in the model (3.60), hi is the sample at time i of the random variable h. Since E [f,] and E [u-\ are equal to zero, taking expected values of (3.60), we get E M = E\pi] = fE[h] E[h] (3.73) Since the noise sources are assumed white and uncorrelated we get S (i,i) gg = E[g*gi] = E[\hi\ }S = E[\hf}S 2 ff If + ol + ol (3.74) Chapter 3. Restoration Filters 33 E N ] 2 = £[N ]+a = E\\h\ \+ol 2 2 (3.75) 2 and S gg(hJ) EbM = E [9*9j\ = EMhASjf (3.76) = E[h*h ] (3.77) 3 for i ^ j From (3.74) and (3.75) it is clear that S (i,i) are the same for every i and E [\pi\ ] are 2 gg also the same for every i. Substituting (3.75) into (3.74), and (3.77) into (3.76), we get S„(i,i) = S (i,j) = E\p*pjSj ] gg (E[\ \ ]-o )S 2 +o 2 Pi v ff 2 u (3.78) for f i^j To obtain S *f(i) g S (i) = E[gU\ gf = E[h*]S = E[h*]S ff (3.79) ff Substituting the corresponding values (3.74) (3.76) and (3.79) in (3.59), E[\h\ )S j + o 2 2 } E{h*MS u E[h*h ]S 2 E[\h\ ]S 2 S} ff E [h*] S fI +o ff 2 u u E[\h\ ]S 2 ff +o 2 u E [V] S n _ E [h*] S ff (3.80) Chapter 3. Restoration Filters 34 For the special case when E [h*hj] — E [h*hk] V i ^ j\ I ^ k, the system in (3.80) can be written as: \h\ S + olU + E[hihi] S„E^. 2 ff = £ W] Sj, (3.81) Therefore all the L^s are the same. Letting Z , = L and from (3.80) we get, E [h*\ E[\h\ ] + (T-l)E{h*h } + 2 j jf (3.82) ; In (3.82), the value of the filter is in terms of the unknowns E [\h\ ] and E [h*hj]. Using 2 (3.78) in (3.83) we get /= E[\p\ }-al + (T-l)E\p* } 2 + f ii Pj (3.84) 4 To calculate (3.84), we use the sample average 1 for Etf] (3.85) for E \pi\ (3.86) and the sample variances estimates 1 J^PiPi and T ( T _u 2 1Z PmPn for E [p* ] Pj (3.87) Substituting (3.85), (3.86), (3.87) into (3.84) we get, j Hi Pi E i 9i * HplPi - °l + (T- (3.88) 1 ) ^ Zi Z&k P]Pk + Since 7f zZPiPi + ^ T T _ ^ zZlZP*jPk = 7f YZp*iPi+ YZYZ p]pk Chapter 3. Restoration Filters 35 (3.89) (3.90) where p = ^ E . P o equation (3.88) can be simplified to, f f = = •-, \P\ -^l 2 f l . + .* rf where g = (3-91) (3-92) Comparing (3.91) with the Wiener filter of the averaged image (3.7), we see that the two are identical except for the 3.7 —-fO , 2 term in the denominator. Discussion of Filters In the previous sections, we have derived the maximum likelihood filter and the Wiener filter for the error in variables model. In the following, we predict the behaviour of the filters by examining their accuracy and stability. A filter is unstable when there are one or more spatial frequencies at which the noise becomes overly amplified by the restoration process. This occurs when the denominator of the filter equation is small in comparison to the noise terms numerator. The value of the noise in the numerator depends on the variance of the noise sources as well as the number of observations. In the following discussion we assume the variance of the noise to be fixed while examining the filters for large and small number of observations. Chapter 3. Restoration Filters 3.7.1 36 Regression, M L and Wiener Filter(H constant) From the last sections, we see that the regression, M L and Wiener filter(with constant H) can be combined into the following general form: IM )I + « M W where k(u) — 0 for the regression filter, = f o r -Tal t h e M L > (3.94) nlter for the conventional Wiener filter. = 2 — Tal + r Sf \w) f *° r ^ e Wiener filter(H constant). In order to simplify the above notation, we denote f(u) J2l\ i\ 2 P Substituting hi + Vi for p,- and /&,-/ + «,• for Zjh*v + Zj\vi\ 2 t as simply / , so + k and rearranging the above equation, we get + k D f + Zfh*Ui + Zjv*u t D (3.96) where D = J2 N i 2 + E > + E < + E H t 2 + * ( - ) 3 9 7 Equation (3.96) indicates that the second and third terms on the R.H.S are two error sources, the first one represents an attenuation on the estimated / and the second represents additive noise. In the absence of noise, the second error term is zero, so the first error term should be minimized to obtain the least biased result. In the presence of noise, the first error term should be chosen so as to compensate for the value of the second term. The compensation is done by setting the k term to a nonzero value. The Wiener filter can be interpreted as the one with "optimal" k value in terms of expected square error. Chapter 3. Restoration Filters 37 Number of Observations Large We first examine the case when the number of observations T is large. With the assumption that h is uncorrelated with the noise sources v and u and that v and u are uncorrelated with each other, the following approximations can be made: Tjh'vi « 0 Tjhivt « 0 (3.98) Tjv*Ui « 0 Substituting (3.94), (3.98) into (3.96) and (3.97), the following can be written for the regression filter: Ef Kf "T,.. „ / / = / - ^ i - i , 12 (3-99) En^i +E?'k-i 2 a Equation (3.99) indicates that the regression estimate is biased, and the amount of bias increases linearly with / . When the value of E f h, is small, the bias becomes more 2 severe. When E f h] = 0, (3.99) becomes / = / - / = 0 (3.100) thus the estimate is biased, but the filter remains stable for E f h = 0. 2 For the M L estimate, combining (3.94), (3.96), (3.97) and (3.98) gives: f = f - ^ f (3-101) which indicates the estimates is unbiased. However for E h = 0, it may be unstable. 2 For the conventional Wiener filter, combining (3.94), (3.96), (3.97), and (3.98) gives, Ef>,l + 2 2 / = / - ^ T „ EfN •„ . „ r . 2 T * f + E f N + Sff ° 2 (3-102) Chapter 3. Restoration Filters 38 This estimate is biased for Y h ^ 0 since as T —> oo equation (3.102) approaches a value 2 of / = / For E ^ 2 z T l i ;E' ' f K f , , / EfN (3.103) 1 2 + Ef k f = 0, equation (3.102) becomes (3.104) / = /-/= 0 which indicates the estimate is stable. For the Wiener filter with constant H , (3.94), (3.96), (3.97) and (3.98) gives el f = f-—Ju _/ (3.105) For Y h ^ 0, this estimate is unbiased as T —• oo since the the denominator increases 2 with T while the numerator stays constant. biased. The estimate is stable for Ytf But when E i l ^ t f — 0 the estimate is = 0, since under this condition, equation (3.105) is becomes / = / - / = 0 (3.106) Number of Observations Not Large In the case when the number of observations T is not large, the terms Yj h* i-> v E f h{V*, Yj h*Ui and Yj i i v u m (3.96) can no longer be ignored. We analyze the stability of each of the estimates by examining (3.97). We combine the random terms in the equation, E f ^ +EfM* +E f ^ (3.107) Chapter 3. Restoration Filters 39 into a single random variable f i , with expected value, thus = 0 + 0 +2V , = Tal 2 (3.108) n = Tal + n where n is real, random with zero mean. It can be seen from (3.94) and (3.97) that the denominator for each of the filters are identical except for a constant bias. In another words, the value of the denominator for each of the filters can be expressed as: D = n + <t> (3.109) where <j> is a real constant defined as follows: regression filter = M L filter Z\hi\ -Tal 2 (3.110) conventional Wiener filter = T,\hi\ -Tal + jfj Wiener filter 2 From (3.108), (3.109) and (3.110), the expected value of D for each of the filters can be written as: D = E i ^ i l + Ta = E\hi\ = Z \hi\ + Tal + "§fj + = £ N 2 2 v + regression filter n M L filter +n 2 2 2 + # £+ rc n (3.111) conventional Wiener filter Wiener filter In (3.111) D is real and all terms on the R.H.S. are non-negative except for n which may be positive, negative or zero. When \h{\ = 0 and since n may attain values which 2 are very near zero, then the M L estimate becomes infinite and the filter is unstable. Chapter 3. Restoration Filters 40 The Wiener filter may also be unstable when E t l^«| = 0 2 a n d when Sff is very large. The regression and the conventional Wiener filters are expected to be stable due to the presence of To . 2 The following table summarizes the behaviour of the filters: biased for filter T large regression yes 2nd ML no A (least) conv. Wiener yes l Wiener(const.H) nof \ biased on ly if E N 3.7.2 stable th s t (most) 3rd 2 = 0. Averaged Filters and Wiener(H random) The estimates based on averaging(3.5), (3.8) and the Wiener with H as a random variable(3.91) has the following general form: /(«) = |p»P + % ) (3.112) where k{u>) =0 for the inverse filter on the averaged image, = <?l/Sff(u) for the Wiener filter on the averaged image, = —<7? + o\j'Sjf(uj) for the Wiener filter with H random with 2 2 x 2 1 2 (3.113) 41 Chapter 3. Restoration Filters Simplifying the notation in (3.112) to P9 \p\ + k (3.114) / = TT 2 Substituting in h + v for p and hf + u for g and rearranging (3.114) gives, h*v + \v\ + k 2 /=/ - D f + h*u -f v*u D (3.115) where D = \h\ + v*h + h*v+\v\ + k 2 2 (3.116) As T increases, since u and v have zero mean, the following approximations can be made, h*v » 0 hv* & 0 h*u w 0 k (3.117) w 0 W i t h the above approximations, when / i is not zero, as T increases, equation (3.115) becomes / = / (3-118) thus, all estimates are unbiased. For h = 0 equation (3.115) becomes: (\v\ + k)f + v*u IvP + k 2 f = f V u |0| + 2 k (3.119) (3.120) Since both the numerator and denominator approach zero as T increases, stability for all filters is not guaranteed. Chapter 3. Restoration Filters 3.8 42 Filter Modification for Improved Stability While the M L and Wiener filters are expected to be more accurate due to their unbiased nature, they are expect to be less stable compared to the conventional Wiener and regression filter. This is confirmed in preliminary simulations results. Thus in the following we modify the M L and Wiener filters to improved their stabilities. It is desirable to keep the characteristics of the conventional filters when the denominator of the filter is near zero while keeping the accuracy of the modified filters when the denominator is not near zero. In another words, the term —Tal the denominator for equations (3.51) and m (3.71) should be removed as the value of the denominator approaches zero. Therefore, the M L filter which is: / = PifH v E,|p,r (3-121) 2 -Tal is modified to the following: f Zi(p*9i) = Zi \ \ - Tal(l - e-*>) 2 Pt where D = £ |p,f - Ta (3.122) 2 v Thus when D —• oo, the M L filter results as desired. But when D —> 0, the regression filter results which is known to be stable. In the M L filter(3.121), when D —* 0 the filter becomes unstable. The Wiener filter(constant H) which is f = Ei(p*&) ZM'-Tal /g 1 2 3 \ + jfj is modified to: / = ZM'-TaKl-e-^ + jjj a 2 where D = ^ \Pi\ - Tal + — 2 2 i b l! (3.124) Chapter 3. Restoration Filters 43 Thus when D —» oo the Wiener filter(constant H) results. When D —» 0, (3.124) becomes the conventional Wiener filter which is stable. For the Wiener filter with random H , we alter equation (3.91) so that it becomes the same as the conventional Wiener filter based on the averaged image(3.7). / = where D = 12 - ' - |H>) + p s r ^ i TSjj 1 | p f - l<r„ + 2 (3.125) Again when D —> oo (3.125) becomes the Wiener filter(random H) and when D —> 0, (3.125) becomes the conventional Wiener filter on the averaged image. Chapter 4 Experimental Results To test the effectiveness of our model, we apply the filters described in the previous chapter to ensembles of simulated images. The results of the experiment are discussed in this chapter. 4.1 Experimental Set Up The original images used for the simulation are shown in Figure 4.5a(64 by 64 pixels) and Figure 4.8a(128 by 128 pixels). The P S F is composed of a deterministic and a random part. The deterministic part is a fixed Gaussian shape, and the random part is generated by white Gaussian noise. Such transfer function resembles those observed in short exposure photographs [4], or that found in images distorted by random vibration. The size of the P S F is confined to an area of 7 by 7 pixels. Since the exact size of the P S F is assumed unknown, a region of size 10 by 10 is used for the measured PSF, i.e., the region containing measurement error is 10 by 10 pixels. The results presented here are obtained with the P S F generated separately from the image, i.e., the correlation between the P S F error and the noise in the image being restored is zero. Simulation for the case where the P S F is contained within the degraded image itself (by installing a point object in Figure 4.5a has also been tested, the results were found to be similar to the ones presented here. As each degraded image is generated, the statistics are updated, thus the actual images are not stored, so storage requirement is very much reduced. 44 Chapter 4. Experimental Results 45 The results for the Wiener filter for the averaged image are obtained by iterating five times equations (3.8). Five iteration are used for all the other Wiener based filters as well. In order to avoid division by zero, at any spatial frequency where the absolute value of the denominator of any filter is less than a threshold value of l O - 1 0 , the result at that spatial frequency is set to zero. The small threshold value is chosen so that it will not impact the performance of the filters. 4.2 Discussion of Results For convenience, in the following discussion, the filters that are based on the averaged data, i.e., inverse(3.5), Wiener(3.7), and Wiener with H random(3.91) will be referred to as the "averaged" filters. The filters derived for conventional models i.e., regression(3.14) and conventional Wiener filter(3.72) will be referred to as the "conventional" filters and the ML(3.51) and Wiener with fixed H(3.71) will be referred to as the "modified" filters. In the simulation, it is found that the restoration results are not alway stable due to the ill-conditioned nature of the problem. The difference between stable and unstable results are usually very obvious. Unstable results are usually very bad and often unrecognizable. When comparing two sets of stable results, the difference in visual quality may not always be obvious, and the "difference picture" (obtained by subtracting the original from the restored image) are sometime required to see the difference. In the following, the results of the M L , Wiener (const ant H) and the Wiener(random H) filter before and after adjusting for ill-conditioning are first compared. Comparison are than made between the modified and the conventional filters, amongst the different averaged filters, and between modified and averaged filters. The root mean square error(RMSE) used for measuring quality of the results are Chapter 4. Experimental Results 46 calculated as the difference between the restored and the original image in the time domain. The R M S E for the degraded image represents the average R M S E for the degraded ensemble. 4.2.1 Adjusting for Ill-Conditioning We first compare the M L , Wiener(constant H), and Wiener(random H) filter(3.51, 3.71, 3.91) to their corresponding counterparts adjusted for ill-conditioning(3.122, 3.124, 3.125). Figures 4.1a to 4.1c show the R M S E of the results as the S N R is varied. The SNR is defined as the ratio of the power of the original signal / over the power of the noise u, i.e., _±_Y Y M SNR - M M 2 y fix v) 2 (4 1) °~l The figures show that adjusting for ill-conditioning improves the stability of the filters. In 4.1a and 4.1b, the M L and Wiener(constant H) filters are stable for SNR of 40 dB or more. After adjusting for ill-conditioning, they remain stable for SNR=30 dB. When the unadjusted filters become stable, they are marginally better than the adjusted filters, although the difference is quite small. Figure 4.1c shows that the restoration with the Wiener(random H) filter is not very successful. Even though the adjusted version is much improved over the unadjusted version, the results are still very much unstable. Since the adjusted version of the filters show improved stability without losing too much accuracy, they will be used in the comparison with the results from the conventional filters and averaged filters in the following sections. 4.2.2 Conventional vs. Modified filters To see the difference between the modified and the conventional filters, the results of M L and the Wiener(constant H) filters are compared to that of the regression and Chapter 4. Experimental Results 47 conventional Wiener filter respectively. The results are shown in Figures 4.2a to 4.2d. The figures show the R M S E of the restorations as the following parameters are varied: 1. number of observations (T) 2. signal to noise ratio (SNR) 3. intensity levels of the original point object from which the P S F is measured (J ) 0 4. amount of randomness in the P S F (cr\ ) h From the figures, it can be concluded that the results of the modified filters are superior to those obtained from the conventional filters. Figure 4.2a shows that the difference between the conventional and modified filters become larger as the number of observations increases. The results of the regression filter are very similar to those of the conventional Wiener filter; and the results of the M L filters are very similar to those of the Wiener(constant H) filter. Although the Wiener based filters are expected to give the lowest R M S E values, when the number of observations is not large, the results show the that the Wiener filters give slightly higher R M S E than the regression and M L filters. This may be caused by the fact that since the true power spectrum of the original image is unknown the results of the Wiener filters are not as accurate as they would have been. As the number of observations becomes large, the results of the Wiener based filters become almost the same as that of the regression and M L . This is expected since when T is large, the o- /Sjf term in the denominator of the Wiener filters become small compared 2 to the other terms. In Figures 4.2b, 4.2c and 4.2d, the results are obtained for T = 50. The results of the regression and conventional Wiener filter are almost identical to that of the M L and Wiener (const ant H) filter. In Figure 4.2b, as the SNR is increased, the difference between the conventional and modified filters decreases. This is expected since the error in the P S F becomes smaller Chapter 4. Experimental Results 48 with increasing S N R . Figure 4.2c shows that the difference between the conventional and modified niters decreases as the signal strength of the P S F is increased. This result is also expected since increasing the S N R of the P S F decreases the effect of the PSF measurement error. Figure 4.2d shows that as the P S F becomes more random, the degraded image becomes increasingly distorted. The restoration results do not become worst, but actually becomes improved with increasing randomness in the P S F . Two explanation for this improvement may be: 1. the least square based estimates become more accurate as the variance of the independent variable is increased (Appendix C), 2. the condition of the problem improves(less ill-conditioned) when the P S F becomes more random since the chances of Ei l^«'| 2 w 4.2.3 0 becomes smaller. Comparison of Filters Based on Averaging Figures 4.3a to 4.3d show the R M S E of the images restored by the filters based on averaging (the inverse(3.5), Wiener(3.7) and the random H Wiener(3.125) filter). From the figures it can be concluded that the Wiener filter outperforms the inverse and random H Wiener in all cases. The random H Wiener filter is found to be very unstable. In Figure (4.3a), all three niters appear unstable at the given values of the parameters. Figure (4.3b) shows the performance of the inverse and Wiener filters improves as the SNR is increased. A t high SNR(50 dB) the difference between the results are not distinguishable. In Figure (4.3c), although the filters are expected to improve as the P S F signal strength is increased, the results are still unstable, and such conclusion can not be made. Figure (4.3d) shows that the results of all the filters are improved with increasing randomness in H (especially for the random H Wiener filter). This may be caused to the fact that the problem becomes less ill-conditioned as the P S F becomes more random. Chapter 4. Experimental Results 49 Modified Filters vs. Averaged Filters 4.3 In Figures 4.4a to 4.4d the results of the Wiener filter(constant H) and the Wiener filter based on the averaged image are compared. In Figure 4.4a, 4.4c and 4.4d the modified filter consistently outperforms the average based filter. Figure 4.4b shows that the modified filter gives lower R M S E than the averaged filter except for the case when SNR=25 dB, where the averaged filter did not fail as badly as the modified filter. 4.4 Image Results The restored images for SNR=40dB and SNR=35dB are shown in Figures 4.5 to 4.10. The quality of the restored images reflects that observed from the R M S E mea- surements. Comparison between the adjusted and unadjusted versions of the M L and constant H Wiener filters show that at 40 dB, the adjusted version (4.6f, 4.61) is slightly less accurate compared to the unadjusted ones (4.6d, 4.6j). At 35 dB, the unadjusted versions(4.7c, 4.7i) are unstable. Comparison between the conventional(4.6a, 4.6g) and the modified(adjusted version) filters(4.6e, 4.6k) shows the modified filters give improved results. No visible difference is observed between the results from the regression based filter (4.6a, 4.6c, 4.6e, 4.7a, 4.7c, 4.7e) and Wiener based filters (4.6g, 4.6i, 4.6k, 4.7g, 4.7i, 4.7k). The results from the average based filters(4.6m, 4.6o, 4.6q) are not as accurate as the non-average based filters(4.6a to 4.6k). The same conclusions can be made from the "girl" image shown in Figures 4.8 to 4.10. 4.5 General Conclusions The experimental results generally agree with that predicted in the last chapter. Restorations with the M L and modified Wiener filter are more accurate but less stable compared to the regression and conventional Wiener filter. The following general Chapter 4. Experimental Results 50 conclusions can be made from the results: 1. After adjusting for ill-conditioning, the modified filters become much more stable and outperform the conventional filters under most conditions. 2. T h e improvements of the modified over the conventional filters are most noticeable when the number of observations is large and the error i n the measured PSF is large. 3. T h e difference between the results of M L and Wiener w i t h H fixed is found to be small(especially when the number of observations is large). 4. Restorations based on the averaged image do not give as accurate results as that obtainable by using the conventional or modified filters. 5. Comparison of the methods based on averaging shows the Wiener filter gives better results than the inverse filter and the random H Wiener filter (which is not very unstable). Chapter 4. Experimental Results 51 ML adjusted for ill-conditioning vs. unadjusted © o degraded « unadjusted « adjusted 4.1a) ML filter (T = 50,1 = 128, a 0 Ah = 10% h(0,0)) Figure 4.1: Comparison of filters adjusted and unadjusted for ill-conditioning Chapter 4. Experimental Results 40 Wiener(H fixed) adjusted for ill-conditioning vs. unadjusted 35 30 25 20 15 10 5 X. 0 25 30 35 40 SNR (dB) — x x degraded » unadjusted x adjusted 4.1b) Wiener filter(constant H) (T = 50,1 = 128, a 0 Ah = 10% h(0,0)) 45 Chapter 4. Experimental Results 53 Wiener(H random) adjusted for ill-conditioning vs. unadjusted 160 1 , 1 r—I— iI :. i• i• i' •' ii / 1 140 i - ' * | : : » 120 ' ' i 100 ' i i ' 00 -.i : i •' •' •" '.i •.» .» i .» i t •" \ • \ :\ ii ;\ ii 80 ; I •. 1 • 1 • I * : > - : -11 i '- \ • \ • \ 60 40 ;» :\ -i • ! . *.» w i : / T *.i '• : ' :» :\ ;* .x • I ; 1 • I - 1 : i - *'••-. - . 20 0 25 '- \~ • \ : i • i • t ; ..-•+ . 1 ' " • " + U 1 1 30 1 35 40 SNR (dB) — degraded .+ unadjusted + adjusted 4.1c) Wiener filter (random H) (r = 50,1 = 128, <r = 10% h(0,0)) 0 Ah 45 50 Chapter 4. Experimental Results 54 Conventional Filters vs. Modified Filters 0 10° I I _1 I I I L _ l ! ' I l ' ' l I ' T 1 1 1—I I I I _1 1 I I I I I I 10i 10 _l I I 2 number of observations degraded o o * x « regression -o ML(adjusted) - * conventional Wiener x const. H Wienerfadjusted) 4.2a) R M S E Response to Variations in the Number of Observations (SNR=40dB, I = 128, o 0 Ah = 10% h(0,0)) Figure 4.2: Comparison of conventional and modified filters(adjusted) I I l _ 10 3 er 4. Experimental Results 55 Conventional Filters vs. Modified Filters 401 , fji 25 , , ' 1 1 1 1 30 35 40 45 50 SNR (dB) degraded o regression <, ML(adjusted) * conventional Wiener const. H Wiener(adjusted) x 4.2b) R M S E Response to Variations in the SNR of the Degraded Images (T = 50, I = 128, a 0 A h = 10% h(0,0)) Chapter 4. Experimental Results 56 Conventional Filters vs. Modified Filters 40 1 I 1 • 1 — - — 35 30 25 20 - 15 10 _ : \ ••» 5 0 50 ; — * '*»»» 1 100 150 200 250 grey level of point object degraded X * conventional Wiener x -X const. H Wienerfadjusted) 4.2c) RMSE Response to Variations in SNR of the PSF (T = 50, SNR=40dB, cr Ah = 10% h(0,0)) 300 Chapter 4. Experimental Results Conventional Filters vs. Modified Filters 401 0 1 1 1 1 1 1 2 4 6 8 10 12 1 — 14 16 Randomness in H ( % h(0,0)) - degraded o regression o ML(adjusted) * conventional Wiener x const. H Wiener(adjusted) 4.2d) RMSE Response to Variations in the Randomness in H (T = 50, SNR=40dB, I = 128) 0 18 20 .pter 4. Experimental Results 58 Restorations Based on the Averaged Image _i 10° 10 101 2 number of observations r r y degraded + inverse filter + Wiener filter ^ rand. H Wiener filter(adjusted) 4.3a) RMSE Response to Variations in the Number of Observations (SNR=40dB, I = 128, o 0 Ah = 10% h(0,0)) Figure 4.3: Comparison of averaged filters i i ii i 10 3 Chapter 4. Experimental Results 59 Restorations Based on the Averaged Image 40 1 1 i i - A 35 / / \ \ \ \ 30 - i 25 * \ i \ \ \ \ » '. \ N \ / \ • \ * X*. < k - / \ / / / ,• / \ ;- . V '+ 20 - i~ \ \ +. \ 15 \ \ \ \ t \ \ \ \ •-.•.\ \ \ 10 i i '•. V ' •. \ ". \ '- \ *\\ \\ -,\ *.\ 1 5 0 25 1 1 • 30 35 40 • 45 SNR (dB) _ degraded .+ inverse filter .+ Wiener filter ..+ rand. H Wiener filter(adjusted) 4.3b) R M S E Response to Variations in the S N R of the Degraded Images ( r = 50, /„ = 128, <7 = 10% h(0,0)) Ah 50 Chapter 4. Experimental Results 60 Restorations Based on the Averaged Image 401 ; . ; n —i : 5Ql 50 1 1 1 1 100 150 200 250 grey level of point object ++ + degraded —+ inverse filter + Wiener filter ,. rand. H Wiener filter(adjusted) 4.3c) R M S E Response to Variations in S N R of the P S F (T = 50, SNR=40dB, o Ah = 10% h(0,0)) 1 300 Chapter 4. Experimental Results 61 5h 01 1 1 0 2 4 i 6 i i i i 8 10 12 14 i 16 Randomness in H ( % h(0,0)) +-+ + degraded + inverse filter + Wiener filter rand. H Wiener filter(adjusted) y 4.3d) R M S E Response to Variations in the Randomness in H (T = 50, SNR=40dB, I = 128) 0 i 18 I 20 62 Chapter 4. Experimental Results Adjusted Wiener(constant H) vs. Wiener Based on Averaging 1 .'/ \\ \\ f \\ l \ 1 1 35 — \\ 1 t 1 \\ l 1 \ l \\ \\ Ii 30 t t \ ii I 1 "TT ••! i n 1 r / "-r• i / ; / s 1 1 1 1 1 1 1 1 - * 1 i l - > \I • 20 - 1 ni \ \ \ \ \ / * # 25 GO i—i—i/i i i i i i i y • y * * V X. 15 - 10 - 5 ' X.. X 0 10° J i i i i i i 1 i 10 1 L_ — L - . 1 -1 .,1 1.1 10 1 1 1 1 2 number of observations _ degraded + Wiener on averaged data const. H Wienerfadjusted) x 4.4a) RMSE Response to Variations in the Number of Observations (SNR=40dB, 1 = 128, a 0 Ah = 10% h(0,0)) Figure 4.4: Comparison of averaged and modified filters(adjusted) 1, 1 ' ' 1 10 3 Chapter 4. Experimental Results 63 Adjusted Wiener(constant H) vs. Wiener Based on Averaging SNR (dB) degraded .+ Wiener on averaged data const. H Wiener(adjusted) 4.4b) R M S E Response to Variations in the SNR of the Degraded Images (T = 50,1 = 128, a „ = 10% h(0,0)) 0 A 64 Chapter 4. Experimental Results Adjusted Wiener(constant H) vs. Wiener Based on Averaging 150 200 250 grey level of point object + » degraded + Wiener on averaged data * const. H Wiener(adjusted) 4.4c) RMSE Response to Variations in SNR of the PSF (T = 50, SNR=40dB, <r Ah = 10% h(0,0)) 300 Chapter 4. Experimental Results 65 Adjusted Wiener(constant H) vs. Wiener Based on Averaging w oo s 6 8 10 12 14 16 Randomness in H ( % h(0,0)) _ degraded .+ Wiener on averaged data •x const. H Wiener(adjusted) 4.4d) R M S E Response to Variations in the Randomness in H (T = 50, SNR=40dB, I = 128) 0 18 20 Chapter 4. Experimental Results b) degraded (SNR=40<f£) 66 c) difference image for b) Figure 4.5: Original and degraded images of "fighter" Chapter 4. Experimental Results e) restored by adjusted M L 67 f) difference image for e) Figure 4.6: Restorations of "fighter" at SNR=40c?B k) restored by adjusted Wiener(const.H) 1) difference image for k) Figure 4.6 Restorations of "fighter" at SNR=40<z5 (cont.) Chapter 4. Experimental Results q) restored by adjusted Wiener(rand. H) r) difference image for q) Figure 4.6 Restorations of "fighter" at SNR=40<f£ (cont.) 69 Chapter 4. Experimental Results e) restored by adjusted M L 70 f) difference image for e) Figure 4.7: Restorations of "fighter" at S N R = 3 5 d £ Chapter 4. Experimental Results 71 g) restored by conventional Wiener h) difference image for g) i) restored by Wiener(const.H) j) difference image for i) k) restored by adjusted Wiener(const.H) 1) difference image for k) Figure 4.7 Restorations of "fighter'' at S N R = 3 5 d £ (cont.) er 4. Experimental Results o) restored by Wiener of average p) difference image for o) q) restored by adjusted Wiener(rand. H) r) difference image for q) Figure 4.7 Restorations of "fighter" at SNR=35<f£ (cont.) Chapter 4. Experimental Results 73 a) original "girl" image b) degraded (SNR=40di?) c) difference image for b) d) degraded (SNR=35(iB) e) difference image for d) Figure 4.8: Original and degraded images of "girl" Chapter 4. Experimental Results a) restored by regression b) difference image for a) c) restored by M L d) difference image for c) e) restored by adjusted M L f) difference image for e) Figure 4.9: Restorations of "girl" at S N R = 4 0 d £ Chapter 4. Experimental Results g) restored by conventional Wiener h) difference image for g) i) restored by Wiener(const.H) j) difference image for i) k) restored by adjusted Wiener(const.H) 1) difference image for k) Figure 4.9 Restorations of "girl" at S N R = 4 0 d £ (cont.) Chapter 4. Experimental Results 76 m) restored by inverse of average n) difference image for n) o) restored by Wiener of average p) difference image for o) q) restored by adjusted Wiener(rand. H) r) difference image for q) Figure 4.9 Restorations of "girl" at SNR=40<1D (cont.) > c) restored by M L d) difference image for c) e) restored by adjusted M L f) difference image for e) Figure 4.10: Restorations of "girl" at S N R = 3 5 d £ Chapter 4. Experimental Results 78 MSm g) restored by conventional Wiener i) restored by Wiener(const.H) . . mMMM: h) difference image for g) j) difference image for i) k) restored by adjusted Wiener(const.H) 1) difference image for k) Figure 4.10 Restorations of "girl" at SNR=35tiB (cont.) Chapter 4. Experimental Results o) restored by Wiener of average 79 p) difference image for o) q) restored by adjusted Wiener(rand. H) r) difference image for q) Figure 4.10 Restorations of "girl" at SNR=35JB (cont.) Chapter 5 Conclusion In this study, we examine the image restoration problem in which multiple observations of a fixed scene are available. The P S F is assumed to be changing with time so that each observation is distorted differently. A noisy measurement of the P S F is assumed available for each of the images. A model has been developed for this restoration problem. In the frequency domain, the model has the form of a linear regression problem but with errors in the independent variable. For our model, errors in the independent variables are transformed into measurement errors in the P S F . In the past, averaging has been a popular method for combining image ensembles. In this study, we explore two approaches which give improved solutions to those obtainable from restorations based on the averaged image. The first approach ignores the errors in the P S F . In this case, restoration can be performed by straight forward linear regression methods. However the results from this approach is theoretically biased. The second approach takes into account the P S F error and overcomes this bias. The maximum likelihood filter and the Wiener filter are developed. Restorations were performed with the filters on simulated images. The results show that restorations which take into account the P S F error are superior to the ones which ignore this error, these in turn are superior to the ones based on averaging. However, the methods which take into account the P S F error have the disadvantage of being unstable when the problem is ill-conditioned. It is found that the difference between the two 80 Chapter 5. Conclusion 81 approaches are most obvious with high number of observations, and large P S F errors. When the P S F error becomes too excessive, the filters that account for the error become unstable due to ill-conditioning. To improve the stability of these niters, these filters are modified so that they become equal to the ones ignoring P S F errors when the problem is ill-conditioned. The results show that under most conditions, accounting for P S F error gives an improved estimate compared to not accounting for P S F error. 5.1 Future Considerations The P S F error model developed in this study is for windowed white noise. One possible extension to the model is for the case in which the P S F is measured separately and not from the degraded image itself. In this case, the P S F noise will not be confined to a small area but will extend over the full size of the image. The model will remains generally the same but the P S F noise will become more significant so that the potential gain from applying this model may become larger. It may be possible to build similar models to account for inaccuracies in the P S F other than white noise. For example, in the case where the P S F represents motion blur, the error in the P S F will likely involve uncertainties in the length of the blur instead of random additive noise. In this case, the model will probably be more complex since the power spectrum of the error is not constant over all frequencies. Such a model, once developed, will be very useful as motion blur is one of the most common type of distortions. Many improvements and generalizations of the model used in this study are possible. For instance, the filters developed for this study have not accounted for the correlation between the additive noise in the degraded images and the noise in the P S F . When the PSFs lies inside the region of the image being restored, this correlation is equal to zero. Filters that take into account this correlation have still to be developed. Another Chapter 5. Conclusion 82 improvement would be to generalize the model to include PSFs that are correlated in time or space. Although the covariance of the random part of the P S F is zero for this study, in many practical applications this may not be true. Filters that will take this correlation into account have to be developed. And finally besides the M L and Wiener filter developed in this study, many other linear and nonlinear estimation techniques may be applied to this error model. These all remain as problems for candidate future work. References [1] D . L . Fried, "Optical resolution through a randomly inhomogeneous medium for very long and very short Exposures," J . Opt. Soc. A m . 56, pp.1372-1378, 1966. [2] R. Barakat and P. Nisenson, "'Finite exposure time astronomical speckle transfer function," Opt. Acta 30, pp.1405-1416, 1983. [3] J . G . Walker, "Object reconstruction from turbulence-degraded images," Opt. Acta 28, pp.1017-1019, 1981. [4] J . W . Goodman, "Statistical Optics, New York: Wiley, 1985. [5] A . Labeyrie, "Attainment of Diffraction Limited Resolution in Large Telescopes by Fourier Analysis Speckle Pattern in Star Images," Astron. h Astrophys., 6, No. 1, 85-87, 1970. [6] M . Fatemi and A . C . Kak, "Ultrasonic B-scan Imaging: Theory of Image Formation and a Technique for Restoration," Ultrasonics Imaging 2 pp.1-47, 1980. [7] W . Vollmann, "Resolution Enhancement of Ultrasonic B-Scan Images by Deconvolution," I E E E Trans. Sonics and Ultrasonics. SU-29, no. 2, pp.78-83, March 1983. [8] W . Kohl and R. Schwarz, "A simple method of realizing the deconvolution of ultrasonic images," Ultrasonics, pp.273-278 Nov. 1981. [9] R . N . Colwell Ed., Manual of Remote Sensing 2nd E d . vol.1, A m . Soc. of Photogrammetry, 1983. 83 References 84 [10] W . K . Pratt, Digital Image Processing, New York: Wiley, 1978. [11] B . R . Hunt, "Digital image processing," I E E E P r o c , Vol. 63, pp.693-708, 1975. [12] A . Rosenfeld and A . Kak, "Picture Processing by Digital Computers,'''' 2nd ed. New York: Academic, 1976. [13] K . R . Castleman, Digital Image Processing, Englewood Cliffs, N J : Prentice-Hall, 1979. [14] T . G . Stockham, Jr. T . M . Cannon, and R . B . Ingebretsen, "Deconvolution through digital signal processing" Proc. I E E E 63, pp.678-692, 1975. [15] A . K . Jain, Fundamentals of Digital Image Processing, Englewood Cliffs, N J : Prentice-Hall, 1989. [16] G . Demoment, "Image Reconstruction and Restoration: Overview of Common Estimation and Structures and Problems," I E E E Trans. Acoust. Speech, and Signal Process. ASSP-37, pp.2024-2036, 1989. [17] R . K . Ward and B . E . A . Saleh, "Restoration of images distorted by systems of random impulse response," J . Opt. Soc. A m . A 2, pp.1254-1259, 1985. [18] R . K . Ward and B . E . A . Saleh, "Restoration of images distorted by systems of random time-varying impulse response," J . Opt. Soc. A m . A 3, pp.800-807, 1986. [19] R . K . Ward and B . E . A . Saleh, "Image restoration under random time-varying blur," Appl. Opt. 26, pp.4407-4412, 1987. [20] R . K . Ward and B . E . A . Saleh, "Deblurring random blur," I E E E Trans. Acoust. Speech and Signal Process. ASSP-34, pp.1494-1498, 1987. References 85 [21] L . Guan, "Restoration of Randomly Blurred Images" P h . D . Thesis, August 1988. [22] L . Guan and R . K . Ward, "Deblurring random time-varying blur," J . Opt. Soc. A m . vol 6, n o . l l , ppl727-1737, Nov.1989. [23] A . G . Qureshi and M . M . Fahmy, 2-D Kalman filtering for the restoration of stochasU tically blurred images," Proc. ICASSP-88, pp.1024-1027, 1988. [24] P . L . Combettes and H . J . Trussell, "Methods for Digital Restoration of Signals Degraded by a Stochastic Impulse Response" I E E E Trans. Acoust. Speech, and Signal Process. ASSP-37, pp.393-401, 1989. [25] D . C . Ghiglia, "Space-invariant deblurring given N independently blurred images of a common object," J . Opt. Soc. A m . A 1, pp.398-402, 1984. [26] B . R . Hunt and 0 . Kubler, "Karhunen-Loeve multispectral image restoration, Part I: Theory," I E E E Trans. Acoust. Speech, and Signal Process. ASSP-32, pp.592-600, 1984. [27] N.P. Galatsanos and R . T . Chin, "Digital restoration of multi-channel images," I E E E Trans. Acoust. Speech, and Signal Process. ASSP-37, pp.415-421, 1989. [28] M . I . Sezan and H . J . Trussell, "Use of a priori Knowledge in Multispectral Image Restoration," I E E E IC-ASSP 1989. [29] Johnston, J , Econometric Methods, McGraw-Hill 1963. [30] Fuller, W . A . , Measurement Error Models, John Wiley & Sons, Inc. 1987. [31] M . G . Kendall and A . Stuart, The Advanced Theory of Statistics, Griffin, London, vol2 1961. References 86 [32] N . Mohanty, Random Signal Estimation and Identification: Analysis and Applicatoins, New York: Van Nostrand Reinhold, 1986. [33] Mendel, J . M . , Lessons in Digital Estimation Theory, Pentice Hall, pp. 92, eq.(11-21), 1987. Appendix A D e r i v a t i o n o f C o m p l e x Regression F o r m u l a The derivation of general regression formula for complex numbers in (3.12) is as follows: min(e ) / 2 where e 2 = \di - h if? i i i Taking derivation with respect to / and equating to zero, de 2 df* f 1 h = i9i Zih'hi 87 i Appendix B Derivation of Complex Wiener Filter The derivation of Wiener Filter for complex signals in (3.59) is as follows: min(e ) 2 where e 2 = = E\£ \f - f\ } 2 £[|X><7.-/I ] 2 C£ iSi-f )(L i9i-f) = = h E E KkgUi » m h - E - E w + n 3 min\h \, /_h (e ) is equivalent to 2 k k \h\*\g\* + \h \e-> > hjgtgj + \h \e^ ^ £ , lh min\h \,lh k k \ k ¥ f c %*<7fc (B.l) \h\e- - "9U+\h \^ "g r + 9 & 2 h k 3/ h h k 2\hk\9tgk + e " j Z , l f c k ( E Wk9i ~ 9lf) + ^ ' ( E Kg g>[ - gtf*) = 0(B.3) k (B.4) Combining B.2 and B.3, 2\h \glg k k =2e- j Z ^ ( E Wk9i " 0*/) 88 Appendix B. Derivation of Complex Wiener Filter 89 Solving for hy. \h \J * i9k9j = Lh k ~ 9kf h 9k9k The above equation shows that each coefficient hk in the filter depends on all other coefficients. Expressed in a matrix form, 9x9\ 9i92 9*92 9292 Li 9lf L 9*2f 2 (B.5) 9N9N 9* f N or SggL — S,9f (B.6) where S (iJ) = E[g* ] S (i) = E[gU] gg af gj (B.7) Appendix C Error of Least Square Estimator The error of the ordinary least square(OLS) estimator becomes smaller as the variation of the independent variable increases. The O L S estimate for Yi = BXi + e,- (C.8) i = 1..N can be written as follows: B E f XjYj E t XiXi Z?Xj(0Xi + ei) E i XiXi Thus the error in the estimate is: V N Ye (c x^r If the m th element (X ) m increased by AX e , m AX m » m is increased by a value of AX . m - 9) The numerator of C.9 is and the denominator is increased by (AX ) . 2 m Assuming that e , the denominator increases faster then the numerator and the total error is m decreased. 90
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Restoration of randomly blurred images with measurement...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Restoration of randomly blurred images with measurement error in the point spread function Lam, Edward W. H. 1990
pdf
Page Metadata
Item Metadata
Title | Restoration of randomly blurred images with measurement error in the point spread function |
Creator |
Lam, Edward W. H. |
Publisher | University of British Columbia |
Date Issued | 1990 |
Description | The restoration of images degraded by a stochastic, time varying point spread func-tion(PSF) is addressed. The object to be restored is assumed to remain fixed during the observation time. A sequence of observations of the unknown object is assumed available. The true value of the random PSF is not known. However, for each observation a "noisy" measurement of the random PSF at the time of observation is assumed available. Practical applications in which the PSF is time varying include situations in which the images are obtained through a nonhomogeneous medium such as water or the earth's atmosphere. Under such conditions, it is not possible to determine the PSF in advance, so attempts must be made to extract it from the degraded images themselves. A measurement of the PSF may be obtained by either isolating a naturally occurring point object in the scene, such as a reference star in optical astronomy, or by artificially installing an impulse light source in the scene. The noise in the measurements of point spread functions obtained in such a manner are particularly troublesome in cases when the light signals emitted by the point object are not very strong. In this thesis, we formulate a model for this restoration problem with PSF measurement error. A maximum likelihood filter and a Wiener filter are then developed for this model. Restorations are performed on simulated degraded images. Comparisons are made with standard filters of the classical restoration model(ignoring the PSF error), and also with results based on the averaged degraded image and averaged PSF's. Experimental results confirm that the filters we developed perform better than those based on averaging and than those ignoring the PSF measurement error. |
Subject |
Imaging systems -- Image quality |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-28 |
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. |
DOI | 10.14288/1.0098397 |
URI | http://hdl.handle.net/2429/29629 |
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 |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1990_A7 L35.pdf [ 16.32MB ]
- Metadata
- JSON: 831-1.0098397.json
- JSON-LD: 831-1.0098397-ld.json
- RDF/XML (Pretty): 831-1.0098397-rdf.xml
- RDF/JSON: 831-1.0098397-rdf.json
- Turtle: 831-1.0098397-turtle.txt
- N-Triples: 831-1.0098397-rdf-ntriples.txt
- Original Record: 831-1.0098397-source.json
- Full Text
- 831-1.0098397-fulltext.txt
- Citation
- 831-1.0098397.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-0098397/manifest