RESTORATION OF RANDOMLY BLURRED IMAGES WITH MEASUREMENT ERROR IN T H E POINT SPREAD FUNCTION Edward W . H . Lam B.Sc. Computer Engg., University of Alberta, 1984 A THESIS 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 OF T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F 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 E N G I N E E R I N G We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A August 1990 © Edward Lam, 1990 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 B-itcb^caJL Ksi^\4JLv} y The University of British Columbia Vancouver, Canada Date AlA^t 2S7 jlYlO DE-6 (2/88) Abstract 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 mea-surement 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 in-stalling 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 measure-ment 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. Exper-imental results confirm that the filters we developed perform better than those based on averaging and than those ignoring the PSF 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 Ensemble Averaging 5 1.3.2 - Multi-channel Images Restoration 6 1.3.3 Heuristic Methods 7 1.4 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 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 3.8 Filter Modification for Improved Stability 42 4 Experimental Results 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 80 5.1 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 . . . . 51 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 SNR=40d£ 74 4.10 Restorations of "girl" at SNR=35d£ 77 vi Acknowledgement I would like to thank my supervisor Dr. Rabab K . Ward for her guidance and en-couragement 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, inhomo-geneities 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 orig-inal 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 nonho-mogeneous 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 PSF is constant, a system with a time varying PSF can not be calibrated once and for all before use. Although when a priori knowledge of the PSF is available and the PSF has a fairly simple math-ematical form, it may be possible to estimate the PSF 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 PSF 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 PSF may be derived by utilizing the knowledge about these structures at different orientations [12]. Obtaining the PSF from the degraded image should give more accurate results than that using blind deconvolution since some information about the PSF is available. However when the PSF 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 PSF since the objects from which the PSF are measured often do not have very strong inten-sity levels. This error in the measured PSF 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 PSF. We shall refer to this restoration problem as the "error in PSF 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 PSF are generally ignored or assumed negligible. The problem can be expressed as: g(x,y) = h(x,y)* f(x,y) + v(x,y) (1.1) = JJ f(x-ct,y - P)h(ot,P)dad8+ v(x,y) Chapter 1. Introduction 3 where * 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 differentia-tion 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 corre-spond 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 conven-tional 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 PSF has been studied in the past. The model for such problems can be written as: where h(x,y) and ni(x,y) represents the deterministic and random part of the PSF respectively, and n,2(x,y) is the random detection noise. Denoting the noise terms ni(x,y) * f(x,y) + n2(x,y) as n(x,y), equation (1.2) becomes: 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) esti-mate 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 PSF, they do not account for the measurement error in the PSF. The difference between these models and the PSF error model will be discussed in the next chapter. g(x,y) = (h(x,y) + n1(x,y))* f(x,y) + n2(x,y) = HX,V) * f(x,y) + ni(x,y) * f(x,y) + n2(x,y) (1.2) g(x, y) = h(x, y) * f(x, y) + n(x, y) (1.3) where n(x, y) is a signal dependent noise. 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 PSF models. A l -though the knowledge of the PSF or its statistical properties have been assumed known in these studies, the PSF 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 im-plementation of the "optimal" solution often requires a prohibitively large amount of computational resources. Methods in which the data is compressed before the restora-tion 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 PSF 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 PSF 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 PSF 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 corre-lation 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]. By 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 corre-lated. 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 PSF for each image is not known, but a noisy mea-surement 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 PSF er-ror). In chapter 4, the performance of the filters are tested with simulated degraded images. Comparisons are made between the filters for the "PSF 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 PSF 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 PSF, 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) = Hx.Y2ylh(x-x',y-y',t)f(x',y') + u(x,y,t) = h(x,y,t)*f(x,y) + 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 10 where * 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 PSF is obtained from a point object in the degraded image, and the intensity of the point object(/0) is known, the computation of p(x,y,t) can be described as in the following. Denoting the unnormalized PSF as p'(x,y,t), we can write: where x0, yQ is the coordinate of the point object and S is the width of the support area given to the PSF. Substituting from (2.1), we get, where v'(x,y,t) = u(x0 + x,y0 + y,t) and I0 is the intensity of the point object in the original image. p'(x, y, t) = g(x0 + x,y0 + y) - ! < * < ! , - | < y < ! (2.2) p'(x, y, t) = I0h(x, y, t) + v'(x, y, t) - ! < * < ! > - f < y < f (2.3) The PSF 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 I0 is unknown, its value can be approximated by: & ' ( 3 , y , < ) (2.5) It can be seen from (2.4) that since Ylxy h(x, y, t) is unity, and the noise v'(x,y,t) has zero mean. Averaging (2.3), we obtain (2-6) 2.1 Frequency Domain Representation Transforming equations (2.1) into the frequency domain, g(ux,ujy,t) = h(ux,uy,t)f(ux,uy) + u(ujx,u;y,t) p(u)x,ojy,t) = h(ux,u>y,t) + v(ux,uy,t) (2.7) (2.8) where g(ux,uy,t), h(wx,uy,t), f(ux,uy), u(wx,u>y,t), p(u}x,u>y,t), v(u>x,uy,t) represents the Fourier Transforms for g(x,y,t), h(x,y,t), f(x,y), u(x,y,t), p(x,y,t) and v(x,y,t) respectively. At each frequency pair u)x,u>y the vectors G, P , H, V, U are defined as: G = g(ux,uy,2) g(ux,uy,T) p(ux,uy,l) p(ux,uy,2) p(ux,uy,T) H = h(ux,uy, 1) h(u>x,Uy, 2) h(ux,u;y,T) (2.9) Chapter 2. Point Spread Function Measurement Error Model 12 u(w*,w y ,l) u(u)x,uy,2) v(u>SiUy,l) v(ux,uy,2) u = V = (2.10) u(u>x,uy,T) v(ux,uy,T) so that, G = Hf+U (2-11) P = H + V where / is a scalar quantity. The first equation in (2.11) has a form similar to the classical linear regression prob-lem: where X is the vector of independent variable, 8 is the unknown parameter to be es-timated, 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 PSF is random then the vector of the transfer functions H in general is composed of a deterministic and a random part. Y = X8+v (2.12) H = Hd + Hr (2.13) Chapter 2. Point Spread Function Measurement Error Model 13 where E[H] = HA (2.14) The noise u(x,y,t) is assumed to be zero mean, Gaussian distributed, and uncor-related 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(ux,ujy,i)*u(ux,uy,j)} = o\ for i = j (2.15) = 0 for i 7^ j 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) (2.16) The autocorrelation of v(x,y,t) at time t, Rvv(Ax, Ay, t), is: Rvv(Ax,Ay,t) = E [v(x,y,t) v(x + Ax,y -f Ay,t)} = E [n(x, y, t) n(x + Ax, y + Ay, t) w(x, y) w(x + Ax, y + A?/)] = Rnn{Ax, Ay, t) S2As(Ax, Ay) (2.17) where As denotes a triangular function of width 25 by 2S with a maximum height of 1. Similarly Ruu(Ax, Ay,t) may be found to be I^Rnn(Ax, Ay, t) M2KM(Ax, 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, Rvv(Ax, y + Ay, t) = f ^S[Ax'Ay)Ruu(Ax, Ay, t) (2.19) I'M1 hM{Ax,Ay) In the case when Ruu(Ax, Ay, t) = o\, 6(0,0,/) where 8 is the Dirac delta function, we get: Rvv(Ax, Ay, t) = | ^ ( 0 , 0 , t) (2.20) Thus Rvv(Ax, Ay, t) is proportional to a delta function. The power spectrum of v(x, y, t), Svv(u>x,Uy), being the Fourier transform of Rvv(Ax, Ay, t), is a constant function. It can be concluded then that S2 Svv(ux,uy) = pj^vl = °l (2.21) where S2 is the support area of the PSF and M2 is the size of the image. As with u(x,y,t), v(x,y,t) is uncorrelated from one image to another, E[v(ux,uy,iYv(ux,uy,j)} = c2 for i = j (-,.22) = 0 for i ^ j The cross correlation of u and v can be shown to be E[u(x,y,t)v(x + x',y + y',t + t')] = I0a2v x' = x0,y' = y0,t' = 0 = 0 otherwise where x0,y0 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 15 2.3 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, P{x,y,t) = h(x,y,t) + v(x,y,t) Substituting the second equation into the first, we get (2.25) g{x,y,t) = h(x, y, t) * f(x, y) + na(x, y,t) h(x,y,t) = h(x,y) + n2(x,y,t) where * is the convolution operator h(x, y) is the known average or expected value of the PSF ni(x,y,t) is the additive noise n2(x,y,t) is the random component in the PSF Combining the two equations of (2.25), we obtain g{x,y,t) = h(x,y) * f(x,y) + n2(x,y,t) * f(x,y) + rn(x,y,t) (2.26) For our problem, (2.27) g(x, y, *) = P{x, V, t) * f(x, y) - v(x, y, t) * f(x, y) + u(x,y, t) (2.28) The differences between our PSF 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 "PSF 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 PSF 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 17 3.2 Restoration Based on the Averaged Image Averaging is the most straightforward method for combining the image ensemble. The averaged degraded image and PSF can be defined as: 1 T 1 t 1 T p{x,y) = ™ £ p ( z , < M ) (3-1) 1 t From chapter two equation (2.9), at each spatial frequency, G = Hf + U (3.2) P = H + V 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 (3.3) = (p-v)f + 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 1 2 <4 = \°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 cr2 is the variance of the noise and SJJ is the power spectrum of the input / . In 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/Sffi where 1 T Sffi = 7j;zZS93 iori = 0 1 t = | % | 2 f o r ^ O (3.9) 3.3 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\gi-hj\2 (3.11) / ,=i 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 esti-mator(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 (3.13) Pi = 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 (3.15) Hi = Ui - Vif 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 explana-tory 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 (3.17) °~h + °v Thus plim f 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 P{G\f) (3.19) where G is the observed vector. This estimate is identical to the maximum a posteri-on'(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, 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: max P(f\G) (3.20) which is the same as max P(G\f)P(f) (3.21) P(G1,G2..GT, P1,P2..PT\f, hu h2..hT) = 1 T —=——exp- £ \Pi - hi\2/2o2v (V2irov)T ( \ / 2 ^ u ) 1 -exp-^gi-hif]2^ (3.22) T Chapter 3. Restoration Filters 22 The corresponding log likelihood function is then: C = - |/n(2™*) - ~ £ \Pi - A,f - ± 52\gi-hif\2 (3.23) 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-hi) + ^(gi-fhi) = 0 oh* °l < Rearranging the above equation and letting (i = i..r) (3.26) (3.27) we get Substituting the above into (3.25), hi = Pi + Xf*9i 1 + A | / P (3.28) / = Letting Zplgi + VEgtgi l + A/2 Ep-Pi + A / * + XfEPig* + A 2 | / | 2 E s t a (1 + A | / | 2 ) 2 - i (3.29) 1 T M9'9 = jiY,9i9i Chapter 3. Restoration Filters 23 1 T 1 T M9'V = y & i W MP*P = TfEptlH (3.30) Equation (3.29) can be written as: i = (l + \\f\2)(Mp.g + \fMg.g) Mp.p + A/*MP.S + \fMg.p + \2\f\2Mg.g which simplifies to the quadratic: \pMg.p - f(\Mg.g - Mp.p) - Mp,g = 0 (3.32) with solutions, h = 0 + yJP + \ f2 = 0-y/F+\ (3-33) MQ-„ — \MV*V The value of 6 can be computed from <r2, o~2 and the observed vectors G and P. 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, fx is always positive 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 Mpg (= -f J2Pi9i) i s positive and to choose the negative / if Mpg is negative. The last advice seems intuitively correct. For our case the variables are complex numbers. To determine Chapter 3. Restoration Filters 24 whether fx or / 2 should be chosen, we calculate Mp*g (3.35) and assuming the last three terms are very small, we get: (3.36) And since J2i I hi | 2 is a real number, the phase of Mp*g is approximately equal to the phase of / . Therefore we choose the / whose phase is closer in value to that of Mp*g. However 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 2 is the support area of the PSF, M2 is the size of the image, and I0 is the gray level of the original point object from which the PSF is determined. With A being small, the coefficient of the quadratic term in equation (3.32) (XMp*g) becomes small in 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 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 / = (3.38) 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 ,„„Q* £, (Pi ~ v*i)(pi - 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 becomes Pr+jPi = Hr+jHi-rVr+jVi Gr+jGi = (HT+jHi)(fr+jfi) + UT+jUt Separating the real and imaginary parts, we get (3.41) (3.42) (3.43) Pr = Hr + Vr Pi = Ht + Vi Gr = Hrfr — Hifi + Ur Gi = Hrfi + Hifr + Ui Assuming Hr, Hi, Vr, Vi, Ur, and Ui have independent normal distributions with variances a\r, a\., a^r, o-\., alr, a2., Pr, Pi, GT, Gi, being linear functions of normal variables, will Chapter 3. Restoration Filters 26 have a multivariate normal distribution. Denoting the correlation of two variables x and y (E[xy]) as Sxy, the following can be written for our problem: ^ P r P r = Shrhr ~T ^VrVr SpiPi = Shihi + SViVi SgrPr ~ frShThr ~ fi^hrhi (3-44) < f^frP. = fr$hrhi - fiShihi SgiPr = frShrhi + fiShrhr SgiPi = frShrhr + fiShrhi Combining the above equations, we obtain ^ 9 r P r = fr{SprPr ~ ^ V r V r ) ~ f i ^ P r P i ^flrP, = frSprPi ~ fii^PiPi ~ SviVi) (3.45) ^9iPr ~ frSPrpi 4" fi{SpTPr SvrVr) SgiPi = fr{$piPi ~ $ViVi) + fiSpTPi For 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 Mxy, the following can be written for our variables: MpTpr = f 5Z{ PrtPrt Mpipi = f J2t PitPit Mprpi — f Y2t PrtPit M9rPr = iZl9rtPrt (3-46) Mgrpi = f S ( 9rtPit Mgipr = f X^t 9itPrt MgiVl = hY7t 9hPit Chapter 3. Restoration Filters 27 Since Pr, Pi, G>, and G, have a multivariate normal distribution, the sample moment 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„. MPrVr - Krvr M, VrPi lVrVi f -Af„..„.. + O2 M, PrVi MPiPi - oiiV. MprpT °~vrVr PrPi MgrPr fr MgrVl . fi . MgtPr _MgtPi _ (3.47) 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 fr,fi by solving the above set of linear equations. Adding the first and last row and subtracting the second from the third, (3.48) (3.49) (3.50) 0 ' fr ' MgrPr + MgiPi _ 0 . Mgipr ~ Mgrp, _ where v = MPrPr + MPiPi - o2v 2 - o 2 Vr ViVi ' fr' . fi . Mgrpr+Mgipi MgiPr-MgrPi Therefore, Ht{Pt9t) (3.51) E t M 2 - Tol 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 28 3.5 The Wiener Filter 3.6 Introduction 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 PSF 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), g2(x).. .gr(x) of an original signal f(x), the estimate which minimizes the mean square error between the original and estimated signal When the conditional densities are not known, such an estimate can not be easily ob-tained. 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 (3.52) is known to be [15] [32] the conditional mean of f(x) given <7,(x),i.e. f{x) = E [f(x)\9l(x),g2(x)..gT(x)] V x (3.53) 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 iLi(x) * F(gi,g2.-gT,Pi,P2~PT)] (3.55) t 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) /(<") is equivalent to minimizing mmE[\f(x)-f(x)f] (3.57) /(*) 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. At 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 fol-lowing solution: ' Sgg(l,2) Sgg(2,l) Sgg(2,2) (3.59) . S,f{T) S3a(T,T) where Sgg(i,j) = E[g*gj] and Sgj(i) = E[g*f] and £, i — 1..T are the coefficients of the Wiener filter to be determined. The above equation is the general Wiener solution applied to any model. In our application, the model is, 9i = hif + Ui (3.60) Pi = 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 With u being independent of / and h and self uncorrelated in time, the autocorre-lation of the observed vector G can be written as: SgB{i,j) = E[g^gj] = E[{h*r + um)(hjf + uj)] = Eih'hjSjA+EMxij] (3.61) Chapter 3. Restoration Filters 31 (3.62) Treating the elements in H as constants, Sggihj) = h*hjSff + a l f o r * = 3 and Saj{i) = 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\h2Sjf U KSff h^hiSfj hfaSff + crl L2 h^Sff • h*ThTSff + a2u hySff Each row i of the above matrix can be expressed as: ZfiVhjSff + alS^Lj = h*Sff ZJihlhjLjSfri + alLi = h*Sff Solving for Li in the above equation, Li = ^ ( l - E M o - ) u 3 = h*a Substituting back into (3.65), h*iSff Sff Zjh.Sffih^ + aKa) (3.64) (3.65) (3.66) (3.67) (3.68) Chapter 3. Restoration Filters 32 Solving for a 5 5 (3.69) 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\2 - Tal + £ ff For conventional mo dels (ignoring u,-), the resulting filter would have been / = T E ' P ^ , (3.72) E f N 2 + ^ Since 5/ / is unknown, the solution can be found iteratively as in (3.8), by using the value of YA Sgg as a starting estimate for S/f. 3.6.2 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 = fE[h] E\pi] = E[h] (3.73) Since the noise sources are assumed white and uncorrelated we get Sgg(i,i) = E[g*gi] = E[\hi\2}SIf + ol = E[\hf}Sff + ol (3.74) Chapter 3. Restoration Filters 33 E N 2 ] = £ [ N 2 ] + a 2 = E\\h\2\+ol (3.75) and (3.78) Sgg(hJ) = E[9*9j\ = EMhASjf (3.76) EbM = E[h*h3] (3.77) for i ^ j From (3.74) and (3.75) it is clear that Sgg(i,i) are the same for every i and E [\pi\2] are also the same for every i. Substituting (3.75) into (3.74), and (3.77) into (3.76), we get S„(i,i) = (E[\Pi\2]-o2v)Sff + o2u Sgg(i,j) = E\p*pjSjf] for i^j To obtain Sg*f(i) Sgf(i) = E[gU\ = E[h*]Sff = E[h*]Sff (3.79) Substituting the corresponding values (3.74) (3.76) and (3.79) in (3.59), E[\h\2)S}j + o2u E[h*h2]SfI E{h*MSS} E[\h\2]Sff + o2u E[\h\2]Sff + o2u E [h*] Sff u E [V] Sn _ E [h*] Sff (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\2Sff + olU + E[hihi] S„E^. = £ 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\2] + (T-l)E{h*hj} + jf; (3.82) In (3.82), the value of the filter is in terms of the unknowns E [\h\2] and E [h*hj]. Using (3.78) in we get / = E[\p\2}-al + (T-l)E\p*Pj} + f4 ii To calculate (3.84), we use the sample average 1 for Etf] and the sample variances estimates 1 J^PiPi for E \pi\ and T ( T _ u 2 1Z PmPn for E [p*Pj] Substituting (3.85), (3.86), (3.87) into (3.84) we get, j Hi Pi E i 9i * HplPi - °l + (T- 1 ) ^ Zi Z&k P]Pk + Since 7f zZPiPi + T^T _ ^ zZlZP*jPk = 7f YZp*iPi+ YZYZ p]pk (3.83) (3.84) (3.85) (3.86) (3.87) (3.88) Chapter 3. Restoration Filters 35 (3.89) (3.90) where p = ^ E . P o equation (3.88) can be simplified to, f = \P\2-^l + rf f = • - , f l . .* (3-91) where g = (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 —-fO2, term in the denominator. 3.7 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 36 3.7.1 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: I M W ) I + « M where k(u) — 0 for the regression filter, = -Tal f o r t h e M L nlter> r2 for the conventional Wiener filter. (3.94) = — Tal + Sff\w) * ° r ^ e Wiener filter(H constant). In order to simplify the above notation, we denote f(u) as simply / , so J2l\Pi\2 + k Substituting hi + Vi for p,- and /&,-/ + «,• for and rearranging the above equation, we get where Zjh*vt + Zj\vi\2 + k D f + Zfh*Ui + Zjv*ut D (3.96) D = J2 N 2 + E + E + E H 2 + * ( 3 - 9 7 ) i > < t 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 rep-resents 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: / = / - ^ i - i , 12 " T , . . „ / (3-99) E f K f E n ^ i 2 + E ? ' k - i 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,2 is small, the bias becomes more 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 h2 = 0. 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 h2 = 0, it may be unstable. For the conventional Wiener filter, combining (3.94), (3.96), (3.97), and (3.98) gives, T2 Ef>, l 2 + Sff E f N2 + E f N 2 + ° / = / - ^ T „ •„ . „ r . * f (3-102) Chapter 3. Restoration Filters 38 This estimate is biased for Y h2 ^ 0 since as T —> oo equation (3.102) approaches a value of E f K f E f N 2 + E f k f For E ^ 2 = 0, equation (3.102) becomes / = / T l i z ; ' ' 1 , , / (3.103) / = / - / = 0 (3.104) 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 h2 ^ 0, this estimate is unbiased as T —• oo since the the denominator increases with T while the numerator stays constant. But when E i l ^ t f — 0 the estimate is biased. The estimate is stable for Ytf = 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*vi-> E f h{V*, Yj h*Ui and Yj viui 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 ^ + E f M * + E f ^ (3.107) Chapter 3. Restoration Filters 39 into a single random variable f i , with expected value, = 0 + 0 +2V 2 , = Tal thus n = Tal + n (3.108) 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: = Z\hi\2-Tal regression filter M L filter conventional Wiener filter (3.110) = T,\hi\2 -Tal + jfj Wiener filter From (3.108), (3.109) and (3.110), the expected value of D for each of the filters can be written as: (3.111) D = E i ^ i l 2 + Ta2v + n regression filter = E\hi\2 + n M L filter = Z \hi\2 + Tal + "§fj + n conventional Wiener filter = £ N 2 + # £ + rc 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{\2 = 0 and since n may attain values which 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^«|2 = 0 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 To2. The following table summarizes the behaviour of the filters: biased for stable filter T large regression yes 2nd M L no Ath (least) conv. Wiener yes l s t (most) Wiener(const.H) nof 3rd \ biased on ly if E N 2 = 0. 3.7.2 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 + % ) 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 (3.112) 2 x 2 2 1 2 (3.113) Chapter 3. Restoration Filters 41 Simplifying the notation in (3.112) to / = TT P 9 \p\2 + k Substituting in h + v for p and hf + u for g and rearranging (3.114) gives, (3.114) / = / - h*v + \v\2 + k D f + h*u -f v*u D (3.115) where D = \h\2 + v*h + h*v+\v\2 + k (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 (3.117) k w 0 With 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\2 + k)f + v*u f = f IvP + k V u (3.119) (3.120) |0|2 + k Since both the numerator and denominator approach zero as T increases, stability for all filters is not guaranteed. Chapter 3. Restoration Filters 42 3.8 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 m the denominator for equations (3.51) and (3.71) should be removed as the value of the denominator approaches zero. Therefore, the M L filter which is: / = v P i f H 2 (3-121) E , | p , r -Tal is modified to the following: Zi(p*9i) f = Zi \Pt\2 - Tal(l - e-*>) where D = £ |p,f - Ta2v (3.122) 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*&) /g 1 2 3 \ ZM'-Tal + jfj is modified to: / = ZM'-TaKl-e-^ + jjj a2 where D = ^ \Pi\2 - Tal + —2- (3.124) i bl! 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). / = p's 12 - - |H>) + r ^ i 1 TSjj where D = |pf - 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 PSF 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 PSF is confined to an area of 7 by 7 pixels. Since the exact size of the PSF 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 PSF generated separately from the image, i.e., the correlation between the PSF error and the noise in the image being restored is zero. Simulation for the case where the PSF 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 unrec-ognizable. 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 do-main. 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 SNR 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., _±_YMYM fix v)2 SNR - M 2 y (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 PSF is measured (J0) 4. amount of randomness in the PSF (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 be-tween 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 conven-tional 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 be-come almost the same as that of the regression and M L . This is expected since when T is large, the o-2/Sjf term in the denominator of the Wiener filters become small compared 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 PSF becomes smaller Chapter 4. Experimental Results 4 8 with increasing SNR. Figure 4.2c shows that the difference between the conventional and modified niters decreases as the signal strength of the PSF is increased. This result is also expected since increasing the SNR of the PSF decreases the effect of the PSF measurement error. Figure 4.2d shows that as the PSF 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 PSF. Two explanation for this improvement may be: 1. the least square based estimates become more accurate as the variance of the in-dependent variable is increased (Appendix C), 2. the condition of the problem improves(less ill-conditioned) when the PSF becomes more random since the chances of Ei l^«'|2 w 0 becomes smaller. 4.2.3 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. At 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 PSF 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 PSF becomes more random. Chapter 4. Experimental Results 49 4.3 Modified Filters vs. Averaged Filters 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 accu-rate 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 sta-ble 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 il l-conditioning, the modified filters become much more stable and outperform the conventional filters under most conditions. 2. The improvements of the modified over the conventional filters are most noticeable when the number of observations is large and the error in the measured PSF is large. 3. The difference between the results of M L and Wiener with 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 degraded © « unadjusted o « adjusted 4.1a) ML filter (T = 50,10 = 128, aAh = 10% h(0,0)) Figure 4.1: Comparison of filters adjusted and unadjusted for ill-conditioning Chapter 4. Experimental Results 40 35 30 25 20 15 10 5 Wiener(H fixed) adjusted for ill-conditioning vs. unadjusted 0 25 30 X . 35 40 SNR (dB) 45 degraded x — » unadjusted x x adjusted 4.1b) Wiener filter(constant H) (T = 50,10 = 128, aAh = 10% h(0,0)) Chapter 4. Experimental Results 53 00 160 140 120 100 80 60 40 20 Wiener(H random) adjusted for ill-conditioning vs. unadjusted 0 25 -1 , 1 r—I— I . | T i : :» / : :\ 1 : *.i i • i ' *.» i '• ;* * • ;» ' •' w i : - . i i » - : i ' • :\ i -i i •' ' . i ' •' •.» ' - .» ' •" i i ' •" .» i i / t i : \ • \ ! . * : i \ i ; i \ i i ; I •. 1 • 1 • I : > --1 :1 '- \ • \ • \ . x --* ' • • - . ..-•+ • I ; 1 • I - 1 : i '- \~ • \ : i • i • t 1 - . . 1 ' " • ; + " U 1 1 30 — degraded .+ unadjusted + adjusted 35 40 SNR (dB) 45 50 4.1c) Wiener filter (random H) (r = 50,10 = 128, <rAh = 10% h(0,0)) Chapter 4. Experimental Results 54 Conventional Filters vs. Modified Filters I I I I ! I l l I T 1 1 1—I I I I 0 10° _1 I L _ l ' ' ' ' _1 1 I I I I I I _ l I I I I l _ 10i 102 number of observations 103 degraded o « regression o -o ML(adjusted) * - * conventional Wiener x x const. H Wienerfadjusted) 4.2a) R M S E Response to Variations in the Number of Observations (SNR=40dB, I0 = 128, oAh = 10% h(0,0)) Figure 4.2: Comparison of conventional and modified filters(adjusted) er 4. Experimental Results 55 Conventional Filters vs. Modified Filters 401 , , , fji ' 1 1 1 1 25 30 35 40 45 50 SNR (dB) degraded o regression <, ML(adjusted) * conventional Wiener x const. H Wiener(adjusted) 4.2b) R M S E Response to Variations in the SNR of the Degraded Images (T = 50, I0 = 128, a A h = 10% h(0,0)) Chapter 4. Experimental Results 56 40 35 30 25 20 15 10 5 Conventional Filters vs. Modified Filters 0 50 100 X * x -X degraded 1 1 I • 1 — - — -: \ _ ••» ; — * '*»»» 1 150 200 grey level of point object 250 conventional Wiener const. H Wienerfadjusted) 300 4.2c) RMSE Response to Variations in SNR of the PSF (T = 50, SNR=40dB, crAh = 10% h(0,0)) Chapter 4. Experimental Results Conventional Filters vs. Modified Filters 401 1 1 1 1 1 1 1 — 0 2 4 6 8 10 12 14 16 18 20 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, I0 = 128) .pter 4. Experimental Results 58 Restorations Based on the Averaged Image _i i i i i i 10° 101 10 2 10 3 number of observations degraded r + inverse filter r + Wiener filter y ^ rand. H Wiener filter(adjusted) 4.3a) RMSE Response to Variations in the Number of Observations (SNR=40dB, I0 = 128, oAh = 10% h(0,0)) Figure 4.3: Comparison of averaged filters Chapter 4. Experimental Results 59 40 35 30 25 20 15 10 5 Restorations Based on the Averaged Image 0 25 30 35 1 1 * i - A i / -\ \ \ \ \ \ \ » \ '. \ \ N \ \ • \ * X*. / / / i i < / / - k / \ ,• \ / ;--. V '+ \ \ +. \ \ i~ \ \ \ \ \ \ t \ \ i •-. \ •. \ i '•. V ' •. \ ". \ 1 '- \ * \ \ \ \ 1 1 • • - , \ *.\ 40 SNR (dB) _ degraded .+ inverse filter .+ Wiener filter ..+ rand. H Wiener filter(adjusted) 45 50 4.3b) R M S E Response to Variations in the SNR of the Degraded Images (r = 50, /„ = 128, <7Ah = 10% h(0,0)) Chapter 4. Experimental Results 60 Restorations Based on the Averaged Image 401 ; . ; —i n : 5 -Q l 1 1 1 1 1 50 100 150 200 250 300 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 SNR of the PSF (T = 50, SNR=40dB, oAh = 10% h(0,0)) Chapter 4. Experimental Results 61 5h 01 1 1 i i i i i i i I 0 2 4 6 8 10 12 14 16 18 20 Randomness in H ( % h(0,0)) degraded +-- + inverse filter + + Wiener filter + y rand. H Wiener filter(adjusted) 4.3d) R M S E Response to Variations in the Randomness in H (T = 50, SNR=40dB, I0 = 128) Chapter 4. Experimental Results 62 Adjusted Wiener(constant H) vs. Wiener Based on Averaging GO \ 1 . ' / I 1 "TT ••! i n 1 r i "-r • i — i — i / i n i 1 1 1 1 1 1 1 1 \ \ \ \ / f / i \ i \ i \ i \ 35 \ \ — \ \ t \ \ \ l 1 1 1 1 l 1 l i \ i 1 / * * 1 # i l -30 \ \ - \ \ \ i I t t i -25 i ; / s > \ I y • y • 20 -X. * * V 15 -10 -5 ' X . . X 0 J i i i i i i i 1 1 L_ — L - . 1 -1 .,1 1.1 1 1 1 1, 1 ' ' 1 10° 10 1 10 2 10 3 number of observations _ degraded + Wiener on averaged data x const. H Wienerfadjusted) 4.4a) RMSE Response to Variations in the Number of Observations (SNR=40dB, 10 = 128, aAh = 10% h(0,0)) Figure 4.4: Comparison of averaged and modified filters(adjusted) 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,10 = 128, a A „ = 10% h(0,0)) Chapter 4. Experimental Results 64 Adjusted Wiener(constant H) vs. Wiener Based on Averaging 150 200 grey level of point object 250 300 degraded + + Wiener on averaged data » * const. H Wiener(adjusted) 4.4c) RMSE Response to Variations in SNR of the PSF (T = 50, SNR=40dB, <rAh = 10% h(0,0)) Chapter 4. Experimental Results 65 Adjusted Wiener(constant H) vs. Wiener Based on Averaging w oo s 6 8 10 12 14 16 18 Randomness in H ( % h(0,0)) 20 _ 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, I0 = 128) Chapter 4. Experimental Results 66 b) degraded (SNR=40<f£) c) difference image for b) Figure 4.5: Original and degraded images of "fighter" Chapter 4. Experimental Results 67 e) restored by adjusted M L 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 69 q) restored by adjusted Wiener(rand. H) r) difference image for q) Figure 4.6 Restorations of "fighter" at SNR=40<f£ (cont.) Chapter 4. Experimental Results 70 e) restored by adjusted M L f) difference image for e) Figure 4.7: Restorations of "fighter" at SNR=35d£ 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 SNR=40d£ 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 SNR=40d£ (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 SNR=35d£ Chapter 4. Experimental Results 78 MSm . . mMMM: 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.10 Restorations of "girl" at SNR=35tiB (cont.) Chapter 4. Experimental Results 79 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.10 Restorations of "girl" at SNR=35JB (cont.) Chapter 5 Conclusion In this study, we examine the image restoration problem in which multiple observa-tions of a fixed scene are available. The PSF is assumed to be changing with time so that each observation is distorted differently. A noisy measurement of the PSF 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 indepen-dent variable. For our model, errors in the independent variables are transformed into measurement errors in the PSF. 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 PSF. 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 PSF 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 PSF 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 PSF error have the disadvantage of being unstable when the problem is ill-conditioned. It is found that the difference between the two 8 0 Chapter 5. Conclusion 81 approaches are most obvious with high number of observations, and large PSF errors. When the PSF 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 PSF errors when the problem is ill-conditioned. The results show that under most conditions, accounting for PSF error gives an improved estimate compared to not accounting for PSF error. 5.1 Future Considerations The PSF error model developed in this study is for windowed white noise. One possible extension to the model is for the case in which the PSF is measured separately and not from the degraded image itself. In this case, the PSF 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 PSF 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 PSF other than white noise. For example, in the case where the PSF represents motion blur, the error in the PSF 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 PSF. 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 PSF 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. Am. 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 Deconvo-lution," 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 ultra-sonic images," Ultrasonics, pp.273-278 Nov. 1981. [9] R . N . Colwell Ed., Manual of Remote Sensing 2nd Ed. vol.1, Am. Soc. of Photogram-metry, 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 Proc , 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, NJ : Prentice-Hall, 1979. [14] T . G . Stockham, Jr. T . M . Cannon, and R .B. Ingebretsen, "Deconvolution through digital signal processing" Proc. IEEE 63, pp.678-692, 1975. [15] A . K . Jain, Fundamentals of Digital Image Processing, Englewood Cliffs, NJ : Prentice-Hall, 1989. [16] G. Demoment, "Image Reconstruction and Restoration: Overview of Common Es-timation and Structures and Problems," IEEE 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. Am. 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. Am. 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" Ph.D. Thesis, August 1988. [22] L . Guan and R . K . Ward, "Deblurring random time-varying blur," J . Opt. Soc. Am. vol 6, n o . l l , ppl727-1737, Nov.1989. [23] A . G . Qureshi and M . M . Fahmy, U2-D Kalman filtering for the restoration of stochas-tically blurred images," Proc. ICASSP-88, pp.1024-1027, 1988. [24] P .L . Combettes and H.J . Trussell, "Methods for Digital Restoration of Signals De-graded 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," IEEE 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 Applica-toins, New York: Van Nostrand Reinhold, 1986. [33] Mendel, J . M . , Lessons in Digital Estimation Theory, Pentice Hall, pp. 92, eq.(11-21), 1987. A p p e n d i x A D e r i v a t i o n of 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(e2) / where e2 = \di - hif? i i i i Taking derivation with respect to / and equating to zero, de 2 df* f = hi9i 1 Zih'hi 87 Appendix B Derivation of Complex Wiener Filter The derivation of Wiener Filter for complex signals in (3.59) is as follows: min(e2) where e2 = E\£ \f - f\2} = £ [ | X > < 7 . - / I 2 ] = C£hiSi-fm)(Lhi9i-f) = E E KkgUi - E - E w + n » 3 min\hk\, /_hk(e2) is equivalent to \h\*\g\* + \hk\e->lh> hjgtgj + \hk\e^h^ £ , ¥ f c %*<7fc min\hk\,lhk \ + \h\e-3/-h"9U+\hk\^h"gkr (B.l) 9 & 2 2\hk\9tgk + e " j Z , l f c ( E Wk9i ~ 9lf) + ^ ' ( E Kgkg>[ - gtf*) = 0(B.3) (B.4) Combining B.2 and B.3, 2\hk\glgk = 2 e - j Z ^ ( E Wk9i " 0*/) 88 Appendix B. Derivation of Complex Wiener Filter 89 Solving for hy. \hk\JLh* = hi9k9j ~ 9kf 9k9k The above equation shows that each coefficient hk in the filter depends on all other coefficients. Expressed in a matrix form, or where 9x9\ 9i92 9*92 9292 Li L2 9N9N SggL — S, 9lf 9*2f 9*Nf 9f Sgg(iJ) = E[g*gj] Saf(i) = E[gU] (B.5) (B.6) (B.7) Appendix C Error of Least Square Estimator The error of the ordinary least square(OLS) estimator becomes smaller as the vari-ation of the independent variable increases. The O L S estimate for Yi = BXi + e,- i = 1..N (C.8) 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 x^r (c-9) If the mth element (Xm) is increased by a value of AXm. The numerator of C.9 is increased by AXmem, and the denominator is increased by (AXm)2. Assuming that AXm » em, the denominator increases faster then the numerator and the total error is 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 |
AggregatedSourceRepository | 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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0098397/manifest