Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Restoration of randomly blurred images with measurement error in the point spread function Lam, Edward W. H. 1990

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

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  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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

Comment

Related Items