RESTORATION OF RANDOMLY BLURRED IMAGES Ling Guan B. Sc. (Electronic Engineering) Tianjin University, China M . A . Sc. (Systems Design Engineering) University of Waterloo 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 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 D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y O F G R A D U A T E S T U D I E S D E P A R T M E N T O F 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 O F B R I T I S H C O L U M B I A August 1988 © Ling Guan r 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 The University of British Columbia Vancouver, Canada Date T % K , • / r „ DE-6 (2/88) Abstract This thesis considers the restoration of images blurred by random point spread functions. The stochastic point spread functions are represented by the sum of a deter-ministic blur function plus an image-dependent noise term. The stochastic properties of such systems are studied. Several linear filters are derived. These filters are based on the following criteria: the Wiener, the minimum variance unbiased, and the constrained least squares. A filter based on the geometrical mean methodology is also presented. The popular Wiener filter, and the minimum variance unbiased filter have already been applied to this problem for the one-dimensional case. We generalize these filters to the two-dimensional case. These filters are found to give good restoration results when the blur and/or the additive noise effects are not severe. When either or both of the above conditions are not satisfied, numerical instability may occur. The constrained least squares filter does not suffer from such ill conditioning. However, the performance of this method needs to be manually controlled in order to yield good restoration results. The geometrical mean filter is formed by combining the modified Wiener filter and the constrained least squares filter. This filter gives the best restoration amongst the four linear-based filters because it does not have the shortcomings the other filters have. A nonlinear filter, the maximum a posteriori filter is also derived and implemented. Since the random blur effect is a nonlinear stochastic process, we expect the nonlinear filtering scheme to simulate the reverse of the blurring process more closely. The superiority of this filter over the modified Wiener filter is demonstrated by examples. In addition, the restoration of random time-varying blurred images is also discussed. ii The Wiener criterion is adopted. The solution of this problem involves the processing of a huge amount of multi-image data. To circumvent this problem we present two approaches to realize this Wiener filter. One is based on the Karhunen-Loeve transfor-mation, and the other is a direct Wiener filter implementation. These two approaches give better results than that given by processing any single frame individually. Since image restoration involves computation of very large systems, The computa-tional effort is an important factor to consider in the implementation of all restoration algorithms. By means of the circulant matrix approximation and the fast Fourier transform, our methods are implemented in the frequency domain. The computational efforts are reduced significantly. iii Table of Contents Abstract ii List of Figures viii List of Tables yii Acknowledgement yiii 1 Introduction 1 1.1 The Need for Digital Image Restoration 1 1.2 Conventional Image Restoration 3 1.3 Distortion by Random Point Spread Function 5 1.4 Objective and Scope of the Thesis 7 2 M o d e l of R a n d o m Distort ion - 10 2.1 Introduction 10 2.2 The Model of Random Distortion 10 2.3 Image Statistics 16 2.4 Noise Statistics 17 2.4.1 Deterministic image field 18 2.4.2 Random image field 22 2.5 Summary 24 3 Restorations by the Wiener A n d the M V U Techniques 25 iv 3.1 Introduction 25 3.2 Restoration by the Wiener Filter 27 3.2.1 Filter Discussion 27 3.2.2 Examples 31 3.3 Restoration by the M V U Filter 38 3.3.1 The Filter Derivation 38 3.3.2 Examples 41 3.4 Some Remarks 45 3.5 Summary 46 4 The Constrained Least Squares and the Geometrical Mean Filters 48 4.1 Introduction 48 4.2 Restoration by the Constrained Least Squares Estimation 49 4.2.1 The Filter Derivation 49 4.2.2 The Algorithm and Experiment Results 53 4.3 The Geometrical Mean Filter 57 4.4 Summary 59 5 Restoration by the Nonlinear M A P Filter 60 5.1 Introduction 60 5.2 The Estimation Procedures 62 5.3 The Computational Algorithm 67 5.4 Experiment Results . 70 5.5 Generalization to Two-dimensions 79 5.6 Generalization to Nonlinear Recording System 80 5.7 Summary 81 v 6 D e b l u r r i n g Random Time-varying B l u r 83 6.1 Introduction 83 6.2 Problem Formulation 84 6.3 Restoration Based on the Average Image 87 6.4 Restoration Based on the Time-averaged Spatial Autocorrelation of the Image 88 6.5 Restoration Based on Image Samples 90 6.5.1 Discretization of the Problem 90 6.5.2 The Wiener Filter Criterion 93 6.5.3 The Time-Space-Separable Case 95 6.5.4 Direct Wiener Filter Implementation 100 6.6 Experiment Results 104 6.7 Summary 121 7 Conclusions 123 7.1 General Discussion 123 7.2 Summary of Results 124 7.3 Suggestion for Future Work 127 References 128 Appendices 134 A Frequency Domain Computation of the M V U F i l t e r i n g 134 B Computing VW^ 139 C Computing VW2 141 vi D The Alternative Expression of the Wiener Filter 143 vii List of Figures 3.1 The image Face 33 3.2 The sine2 random point spread function: a) the mean of the blur, b) the random function 33 3.3 Distorted version of Figure 3.1 34 3.4 Restoration of Figure 3.3 by standard Wiener filter 34 3.5 Restoration of Figure 3.3 by: a) the modified Wiener filter, b) the mod-ified parametric Wiener filter 35 3.6 Random blur: noise confined to the width of the blur function 36 3.7 Comparison between the modified Wiener filter and the Backus-Gilbert technique: a) The original image; b) The blurred image; c) Restora-tion by the Wiener-based filter; d) Restoration by the Backus-Gilbert technique 37 3.8 Distortion and restorations of the image Number: a) The original im-age; b) The blurred image; c) Restoration by the MVU-based filter; d) Restoration by the Wiener-based filter 42 3.9 Filter performance as ajj changes 43 3.10 Filter performance as ax changes 43 3.11 Filter performance as a2 changes 44 3.12 Filter performance as the image size changes . . 44 viii 4.13 The performance of the constrained least squares filter: a) distortion; b) restoration by the Wiener-based filter; c) restoration by the constrained least squares filter: constraint satisfied; d) restoration by the constrained least squares filter: manual adjustment 56 4.14 The restoration by the geometrical mean filter 59 5.15 An one dimensional example, a) The original signal; b) the blurred signal; c) Restoration by the M A P filter; d) Restoration by the Wiener-based filter 74 5.16 Image F16: distortion and restorations, a) The original image; b) the blurred image; c) Restoration by the M A P filter; d) Restoration by the Wiener-based filter 75 5.17 Image Face: distortion and restorations, a) The original image; b) the blurred image; c) Restoration by the M A P filter; d) Restoration by the Wiener-based filter 76 5.18 Subimages of Face: the right eye. a) The original image; b) the blurred image; c) Restoration by the M A P filter; d) Restoration by the Wiener-based filter 77 5.19 The difference images of the right eye. a) The difference between Figure 5.18 a) and 5.18 b); b) The difference between Figure 5.18 a) and 5.18 c); c) The difference between Figure 5.18 a) and 5.18 d); d) The difference between Figure 5.18 a) and the noisy version of 5.18 c); 78 6.20 Image Eye 107 6.21 The four distorted versions of Eye 108 ix 6.22 Restoration of Figure 6.21 by jointly processing the data, a) Restoration . by Karhunen-Loeve transformation; b) Restoration by direct Wiener filter 109 6.23 a) Individual restoration of Figure 6.21 a); b) Restoration based on time-averaged image 109 6.24 The difference pictures of example 1: a) The difference between Figure 6.20 and 6.21 a); b) The difference between Figure 6.20 and 6.22 a); c) The difference between Figure 6.20 and 6.22 b); d) The difference between Figure 6.20 and 6.23 a) 110 6.25 Image Face 112 6.26 The four distorted versions of Face 113 6.27 Restoration of Figure 6.26 by jointly processing the data, a) Restoration by Karhunen-Loeve transformation; b) Restoration by direct Wiener filter 114 6.28 a) Individual restoration of Figure 6.26 a); b) Restoration based on time-averaged image 114 6.29 The difference pictures of example 2: a) The difference between Figure 6.25 and 6.26 a); b) The difference between Figure 6.25 and 6.27 a); c) The difference between Figure 6.25 and 6.27 b); d) The difference between Figure 6.25 and 6.28 a) 115 6.30 Image F16 117 6.31 The four blurred versions of F16 118 6.32 Restoration of Figure 6.31 by jointly processing the data, a) Restoration by Karhunen-Loeve transformation; b) Restoration by direct Wiener filter 119 x 6.33 a) Individual restoration of Figure 6.31 a); b) Restoration based on time-averaged image 119 6.34 The difference pictures of example 3: a) The difference between Figure 6.30 and 6.31 a); b) The difference between Figure 6.30 and 6.32 a); c) The difference between Figure 6.30 and 6.32 b); d) The difference between Figure 6.30 and 6.33 a) 120 xi List of Tables 3.1 RMS errors of the Wiener technique 32 3.2 RMS errors of the B - G technique 36 3.3 RMS errors of the M V U filter . . . 41 4.4 RMS errors of the constrained least squares filter 55 5.5 RMS errors of the M A P technique 72 6.6 RMS errors of the 1st example in time-varying case. 106 6.7 RMS errors of the 2nd example in time-varying case I l l 6.8 RMS errors of the 3rd example in time-varying case 116 xii Acknowledgement I would like to thank my supervisor Dr. Rabab K. Ward for her encouragement, guidance, and all the help she gave me throughout this research. Her dedication is sincerely appreciated. I would also like to thank the other members of my supervisory committee: Dr. Cyril Leung and Dr. Harvey Richer, for their valuable time and discussions. Special thanks to Mr. Robert Ross, the computer system manager, for the many helps he provided during the research. The research assistantship from the Natural Science and Engineering Research Council of Canada, and the teaching assistantship from the Department of Electri-cal Engineering, University of British Columbia are gratefully acknowledged. On the personal side, I would like to thank my wife, Weiwei, for her patience and support during my doctoral program. xiii TO Ml/ PARSMTS x i v Chapter 1 Introduction 1.1 The Need for Digital Image Restoration Considerable effort has been dedicated in recent years to the problem of digital image restoration. The image restoration problem arises due to the necessity to improve images that are not of optimum utility for the purposes at hand. Millions of images are created every day. Most of them are of very fine quality. Some of them are of lesser quality; and of these that are of lesser quality a certain subset are of such importance, or are so unique, that it is feasible to consider the techniques by which the image quality may be enhanced or by which the degrading phenomena may be removed and the undegraded image restored. The areas of space imagery, biomedical imagery, industrial radiographs, photorecon-naissance images, television, remote sensing, and several multispectral or other esoteric forms of mapping scenes or objects onto a two-dimensional format, are all likely can-didates for digital image restoration. Yet many nonnatural images are also subject to digital image processing. By nonnatural images one might refer to two-dimensional formats for general data presentation for more efficient human consumption. Thus range, range-rate planes, range-time planes, voiceprints, sonargrams, etc., may also find themselves subject to general two-dimensional enhancement and restoration. Restoration is an image processing technique which attempts to simulate the reverse of the degrading process, remove the blur effect caused by the point spread function of 1 Chapter 1. Introduction 2 the optical system and clear the noise involved in the image formation. The mathematics of the image restoration problem can be posed in either discrete or continuous variable form. Thus , many image restoration schemes can be implemented either by digital computers or by analog (optical) computers. This duality gives rise to an obvious question: Since an image by nature is usually an analog (optical) quantity, why should there be such a growth in interest in digital implementation? Several facts lie in the heart of the answer to the question. In particular we note the following: 1. Increasingly, the natural form in which an image is formed and acquired is not analog but of discrete or digital form. A number of electronic image sensors, such as channel-plate photomultipliers, are wholely discrete. Further, digital data transmission is increasingly popular so that even if the original image is not discrete, there are genuine bonuses to be gained in sampling and conversion to discrete form, particularly for bandwidth reduction, noise immunity, etc. 2. The past ten years have seen a dramatic decrease in the price/performance ca-pabilities of digital hardware. Computers that are considered medium scale by today's measures would have been considered large- or super-scale computers ten years ago, and today's super-scale computers are capable of performing in seconds image processing operations that would have taken hours ten years ago. 3. Equally dramatic have been the innovations in processing algorithms. The dis-covery of fast transform and fast convolution algorithms have resulted in several orders of magnitude improvement in the computation times required in many image restoration and enhancement processes. 4. Increasing sophistication in digital image input/output systems has resulted in Chapter 1. Introduction 3 eliminating what used to be a major source of headaches in digital image pro-cessing. Modern scanners, digitizers, and displays are accurate, reliable, and increasingly less expensive. The usage of interactive, on-line video displays is in-creasing; there are several installations combining interactive video displays with time-shared or multiprogrammed computers, giving a flexibility in processing and problem turnaround that cannot be equaled even in the most well-organized op-tical systems. 5. The digital computer is inherently more flexible and is more readily used in operations that involve nonlinearities and/or decision-making process. Some of the most exciting work in image restoration involves the use of basis functions, functional decompositions, and nonlinear and/or iterative techniques; only with great difficulty can these operations be done optically, if at all. 1.2 Conventional Image Restoration A general image formation equation is usually written in the following mathematical form g{x,y) = j j h(x,y;x1,y1\f(xi,y1))dxidy1 (1.1) where / represents an unknown object image, g represents the corresponding measured image, h represents the point spread function of the image formation system. If the imaging system is linear, then g(x,y) = f f /i(x,y;x 1 ,y 1 )/(x 1 ,y 1 )(ix 1 dy 1 (1.2) J J — OO Furthermore, it is possible to envision an image formation system that acts uniformly across the image and the object planes; hence the point spread function is independent Chapter 1. Introduction 4 of position. Such an image formation system is said to possess a shift-invariant point spread function. The output is then a function only of the difference in variables in the coordinate system, and the corresponding equation for image formation is, 9{x,y) = f I h{x - xi,y - yt; f{x1,yl))dxldyl (1.3) J J —OO If we assume that the imaging system is both linear and shift-invariant, and that the point spread function h is deterministic, we obtain the conventional image formation equation. 9{x,y) = J J^Hx - xi,y-yi)f(x1,y1)dxidyi + n(x,y) (1.4) = h{x,y) * f{x,y) + n(x,y) where * represents convolution. In the above equation, the term n(x, y) is an additive noise which arises due to errors in transmission, or the uncertainty inherent in the electronics of the imaging sensors and detectors. It is equation (1.4) that the majority of the restoration methods have been devoted to. The earliest attempts toward image restoration of the conventional class (1.4) were based on the concept of inverse filtering in which the degrading transfer function is inverted to yield a restored image [1] [2] and [3]. It should not* be surprising that the inverse filtering performs poorly in the presence of noise since the filter design ignores the noise process. Various other techniques have been applied to image restoration ever since. These techniques give better restoration quality than the inverse filter does. Among them, the important methods includes the Wiener filter [4] [5] [6], the Homomorphic filter [7], the geometrical mean filter [7], the constrained least squares filter [8], the maximum entropy filter [9], and the pseudo-inverse filter [10] [11]. The powerful Kalman filtering technique has also been introduced to restoration of the conventional class [12] [13] [14]. A comprehensive discussion of restoration methods under the above "usual" conditions can be found in[l5], [16], and [17]. Chapter 1. Introduction 5 1.3 Distortion by Random Point Spread Function Although the conventional distortion model sets up the basis for most of the popu-lar restoration methods, there are situations when the degradation caused by the image formation system is much more complicated than that represented by the model (1.4). In this thesis we address an image restoration problem in which the unpredicted effects occur in the imaging system itself, i.e., the point spread function of the imaging system is itself random. Randomness of the pupil function of an optical system is not uncommon in the real world. Such phenomena may result from: o Random amplitude and phase fluctuations of the pupil function of the optical sys-tem. These may result from external effects such as optical propagation through atmospheric turbulence, or dust particles on the optical surface. © Underwater imaging or communication, and more generally any scattering prob-lems involving random media [18]. © Scanning image-formation system (scanning microscope, microdensitometer, or laser printer) for which the width of the beam is random^ © Random vibration of the optical instrument (camera) relative to the object. © X-ray imaging, where the image formed by a phosphor screen-film system results from the stochastic amplification and the random scattering of quanta. All such effects make the impulse response function of the imaging system a random function. Early attempts addressing the issue of restoration under random blurring effects can be traced back to Slepian's work in 1967 [19]. Major contributions since include: Chapter 1. Introduction 6 • Frank (1969) applied the Wiener filter to communications through random dis-persive channels [20]. • Von Heide (1979) developed restoration schemes allowing for random parameters in the system's impulse response function but assuming that all randomness obey Gaussian statistics [21]. • Kassam et al (1980) introduced a minmax approach to solve this problem [22]. • Other authors investigated the problem of restoration when the point spread function is unknown by the so-called blind deconvolution [23] [24]. Recently, Ward and Saleh in a series of publications [25-28] presented several restora-tion techniques which restore images suffering from random degradation. The schemes they introduced include the digital Wiener filter, the minimum variance unbiased (MVU) filter, (both are of one-dimensional), and the Backus-Gilbert filter. The Wiener filter has been generalized to the two-dimensional image restoration by Guan and Ward [29]. A nonlinear filter based on the maximum a posteriori methodology by Guan and Ward [30] [31] and a Kalman filter approach by Qureshi and Fahmy [32] have also been reported in the literature. Due to the introduction of the uncertainty in the image formation system, the conventional distortion model is no longer suitable for the purpose of restoration of randomly distorted images. In [25] Ward and Saleh suggested a modeling scheme based on the separation of the stochastic part of the point spread function from the deterministic part. Their proposed Wiener and M V U based filters are both derived using this model. The concept of the model is simple, but its applicability is not. Another important factor in image restoration is the cost of computation. The restoration of a digital image of 100 x 100 pixels (a very small one!), using linear based Chapter 1. Introduction 7 filters in the space domain, would require solving a system of linear equations with size 10000 x 10000. This involves computations of O(1006), obviously a computationally prohibitive approach. Thus we must seek efficient methods. 1.4 Objective and Scope of the Thesis In this thesis we will investigate the restoration of images distorted by a random blur function using the digital computer. We will show that the significance of the model proposed by Ward and Saleh goes far beyond the application of the filters they developed. In fact, application of this modeling scheme enables most of the restora-tion algorithms established under the "usual" conditions to be effectively used in the restoration of random blurred images. We will propose several iterative algorithms to solve the problem based on this model. As in the conventional restoration case, linear filters are utilized whenever possible. The reason for the choice of linear filters is due to their simple mathematical representation, economical computational efforts, and the restoration quality they can provide when the nonlinear phenomena are not severe. Hence four linear filters are modified and presented in this work: the Wiener filter, the M V U filter, the constrained least squares filter, and the geometrical mean filter. The geometrical mean filter we consider here is different from the filter introduced in [ 7 ] in that, instead of using the inverse filter, we use the constrained least squares fil-ter as part of the filter. Since the random distortion constitutes nonlinear features to the problem, the expensive but more powerful nonlinear filters are expected to yield higher quality restoration. Therefore, we construct a nonlinear filter, the maximum a posteriori (MAP) filter. In addition to the above work we will also address an image restoration problem in which, the point spread function is not only random, but also time-varying. This Chapter 1. Introduction 8 problem is closely related to, but more difficult to solve than, the problem mentioned before because here we need to restore the object image from a series of recorded images. The Wiener filter criterion is adopted to perform the restoration. Two different approaches are used to realize the filter. One is based on a Karhunen-Loeve-type transformation which is applicable for the case where the spatiotemporal correlation of the random processes is separable. The other is applicable under the more general conditions. We also consider the complexity of the computational algorithms. We will show that under certain statistical conditions, by using the circulant matrix approximation, and the fast Fourier transform, (FFT) , the computational complexity of restoring an „ M x M image will be reduced significantly. The reduction per iteration is from 0 ( M 6 ) to 0(M2logM) for the linear based filters except the M V U filter. The reductions for the nonlinear M A P filter and in the time-varying case will be indicated in the respective chapters. The thesis is arranged in seven chapters. In Chapter 2, the model proposed by Ward and Saleh will be presented, and some useful statistics will be given. Chapter 3 deals with the modified Wiener-based filter, and the M V U filter. Chapter 4 discusses the the constrained least-squares filter and the geometrical mean filter. The presentation of the important nonlinear filter, M A P , will be covered in Chapter 5. In Chapter 6, we discuss the restoration of images degraded by a random time-varying impulse response. We will set up a model to represent the degradation, and use the Wiener filter criterion to perform the restoration. In Chapter 7, we conclude the thesis and point out some possible future work. In this chapter we have used the terms "degradation", "distortion", and "blur" to describe the effects of the imaging formation procedure without distinguishing the subtle differences amongst them. We will continue to use these terms interchangingly Chapter 1. Introduction throughout the thesis. Chapter 2 M o d e l of R a n d o m Distortion 2.1 Introduction As indicated in Chapter 1, Ward and Saleh proposed a simple and efficient scheme to model the randomness of the blur function which enables the generalization of many linear filters to the random blur problem. The details of their original work was im-plemented for modeling the one-dimensional image (or signal) formation problem. Al -though they indicated that their modeling scheme could be generalized to the two dimensional case, the details of the generalization has not been done. In this chapter, the model proposed by Ward and Saleh is extended to the two-dimensional image formation case. The important statistics associated with the model for the two-dimensional case are derived. The statistics include the correlation and the covariance matrices of / , Rf and Kf, and the correlation and'covariance matrices of the noise, Rn and Kn. 2.2 The M o d e l of R a n d o m Distortion The general model of an image formation system in which the point spread function is random can be written as sKx>y) = 1 1 h ( x - x u v - y i ) f { x u y i ) d x i d y i + n2{*,y) (2.5) J J —OO -- h * / + n 2 10 Chapter 2. Model of Random Distortion 11 where h is the random point spread function and n 2 is the additive noise. In (2.5) we have assumed that the optical system is linear shift-invariant. Ward and Saleh suggested that the random point spread function can be treated as the sum of the expected value of the blur function and a zero mean random fluctuation around the expected value. Thus the function h is represented by the formula h = h + ni (2.6) where , h = E[h] is the expected value of h, and is the random fluctuation. This separation leads to the imaging equation g = h*f + n (2.7) where n = ni * / + ri2 Then the blurred image g consists of two parts, the deterministic term h * f, and the stochastic term n = ni * f + n2. Thus, we can regard the distortion in expression (2.7) as the effect of the degradation of a deterministic blur h plus an additive image dependent noise n. The deterministic blurs we shall use are the sine2, Gaussian and Motion blurs. Please note that these blur functions do not contain negative values. Thus the random blurs will not have significant zero-crossings. Although estimation methods based on the solution of the integral equation (2.7) are possible, our efforts are directed towards digital processing. It is therefore expedient to put the problem in discrete form right away. In the discrete domain, we replace the function f(x,y) by a matrix F of dimension N x N; the mean of blur h by a matrix H of dimension L x L\ and the functions Chapter 2. Model of Random Distortion 12 g(x,y), n(x,y), and n2(x,y) by matrices G, N, and N2 of size M x M, respectively, where M = N + L - 1. The structure of our problem is such that the matrices are of block Toeplitz form. This enables us to use the circulant matrix approximation and F F T in the filter im-plementation in the following chapters. We augment columns and rows of zeros to the image F to extend it to a matrix of M x M. Thus the original value, i(i,j<N fii = (2.8) 0. otherwise For the sake of processing, the two-dimensional data are converted into one-dimensional data by the usual stacking technique [16]. As a result of the modification / ( i , y), g(x,y), n(x,y), and U2(x,y) may be represented by vectors / , g, n, and n2, of length M 2 . 1 In addition the functions h(x,y) and n%(x,y) are replaced by block Toeplitz matrices H and Ni of dimensions M 2 x M2. The discretization follows the following definition: fujw = / ( i x A x , j v A y ) , 9im.i, = g{ix&x,ivAy), nu,i, = n{ixAx,ivAy), Nk,., = n2(icAx,iyAy), {H]ij = h[{ix - jx)Ax,{iv - jv)Ay], [Ni]ij = ni[[iK-jt)&x,(iv-jy)Ay]. (2.9) in which j = M{jv - 1) + jx, i = M{iv - 1) + i, 1 In this thesis we use the same symbols for the continuum and the discrete vector forms of / , g, n, and ii?. To distinguish between the two forms we shall use the two independent variables (x, y) when we refer to the continuum form, e.g., f(x,y) is the continuum form, and / is its discrete (lexicographically ordered) vector form. Chapter 2. Model of Random Distortion 13 where t. = 1,2,... M , *„ = 1,2,... M , _i« = l,2,-... M , = 1,2,... M . In equations (2.9) A x and Ay are the pixel lengths in the x and y coordinates, respec-tively. Equations (2.7) may now be written in the matrix vector form g = Hf + n (2.10) where n = Nif + n2 Since the columnwise ordering scheme is used, the structure of the vectors / , g, n, and n 2 are, respectively, fi h M 9 = 9i 92 9M f\f,M 9l,l 9 M, I 91,2 (2.11) (2.12) 9M,M Chapter 2. Model of Random Distortion 14 i i . i n = ni n2 "•1,2 NM,M (2.13) and U2 = H i [n2]2 \n2}M ( n 2 ] l , l [«2]M ,1 [ « 2 ] l , 2 (2.14) [^2]M,M In (2.11) to (2.14) fk, gk, nk, and [n2]k represent the fcth column in matrix F, G, N, N2, respectively. The structure of the matrix H is H = H2 Hi . H2 H2 Hr Hr 0 . . . Hx Hr • H2 H\ (2.15) Chapter 2. Model of Random Distortion 15 where Hi hj2 hji hj2 hir 0 hjr . hjt hj2 0 hji hjr . . hj2 hji and r = M-N + 1 In equation (2.15), matrix H is not only block Toeplitz, but also block circulant. This is possible due to the result of the augmentation of the object image F. We also write matrix Ni as a block matrix W W {N1)2M (2.16) where Kivo.-.li.M [ ( ^ i ) . , y ] 2 , M For our problem, we assume that Ni is stationary in space. Then matrix Ni here is Chapter 2. Model of Random Distortion 16 also block Toeplitz, i.e. {Nfrj = (iVi),-i,,-i and for all relevant i,j,k,l. Equation (2.10) is the distortion model we are seeking. This model sets up the foundation for the restoration algorithms discussed in the following chapters. 2.3 Image Statistics Most of the restoration algorithms involve statistical consideration of the noise and signal processes. The important statistics include the mean vector, the covariance matrix, and the correlation matrix. The correlation matrix of the object image is denoted by Rf. Since the object / is a lexicographically ordered image vector, then Rf = E[ffT) (2.17) [Rf]i,i [RAi,2 ••• [Rf]lM [RAM,I • • • \RI\MM _ where each matrix is of size M x M, [Rfh = Elfifi], (2-18) fi, fi are the tth and j'th row partition of the ordered vector / , and the superscript T denotes transposition. Since [R/]ij — [Rf]ji, then (2.17) is a symmetric matrix. The Chapter 2. Model of Random Distortion 17 diagonal blocks [R/]u describe the variation of each row partition (column of the original matrix F) with itself. Blocks off the diagonal describe the "correlation" of columns in the matrix separated by |t — j\. Likewise, the covariance matrix is defined by Kf = E[(f - J)(f - ff) (2.19) \Kf]lA \Kf]i,2 • • • \Kf}lM _ lKfh,l \Kf)2,2 • ' 1 [Kf]M,M where J is the mean vector of the object image / . The mean vector J can be either a constant vector, or a variable vector. • If / 7^ A*) a constant vector, the model is nonstationary. However, the quantity / — / may still be stationary. Therefore the correlation matrix Rf is not block Toeplitz, but the covariance matrix Kf may be block Toeplitz. • If 7 = A4, a constant vector, the underlying random process for / is stationary. In this case both Kf and Rf are block Toeplitz. It will be mentioned, whenever necessary, what model assumption is used in filter design. 2.4 Noise Statistics As far as our problem is concerned, the detection noise n2 is assumed independent of / and the randomness in the point spread function. The mean vector of the noise is: E[n] = E[Nif + n2] (2.20) = ElNifl + Eim] Chapter 2. Model of Random Distortion 18 It is also assumed that the object / is uncorrelated with N\, and Nt is a zero mean random fluctuation. Then we have: E[n] = £ [ / V i ] £ [ / ] + E[m\ = E[n2] (2.21) To derive the correlation matrix Rn and the covariance matrix Kn of the noise n, two cases need to be considered. The first case is when the image / is deterministic, (as assumed by the M V U filter, and the constrained least squares filter); the second is when / is a sample of a stochastic process (as assumed by Wiener filter, for example). These two considerations are discussed separately. 2.4.1 Deterministic image field If the image / is regarded as deterministic, the correlation matrix of n is Rn{f) = E(nnT) = £ [ ( 7 ^ / + n 2 )( iV 1 / + n 2) T] (2.22) where Rn{f) indicates that the correlation matrix of n is dependent on / . Since Ni and n 2 are independent of each other, it is obvious that equation (2.22) can be written as Rn(f)^E[(Nlf)(Nlf)T] + E(ninl) It is easy to show that Rn{f) is a matrix with the following structure (2.23) Rn(f) = [i*n(/)]l.i [ JM / ) l u [Rn(f)h,l [Rn(f)]*.2 [Mf)} IM (2,24) [Rn{f)]M,l ••• \Rn{f)]l,M where [Rn{f)]i,j is a submatrix. To find the elements of {Rn(f)]ij, we consider the Chapter 2. Model of Random Distortion 19 (i — 1)M + /th row vector of matrix Nx: (2.25) in which [(JVi)i,,]i is the /th row in the submatrix {Ni)iiq (referring to (2.16)). Then the entries in Rn{f) can be expressed by [{Rn{f))i,m}i,n = ^ { [ ( A r i ) ( , - i ) M + J ] T [ ( i V 1 ) ( m _ 1 ) M + „ / ] } + ^ { [ n 2 ] J , , [ n 2 ] n i m } (2.26) = £ { / T [ ( i V 1 ) ( , - _ 1 ) M + l ] r [ ( J V l ) ( » - l ) W + » / l } + ^{[n>]l1i[»2]n,m} = / T ^ { [ ( / V 0 ( i - 1 ) M + I ] T [ ( i V 1 ) ( m - i ) M + n ] } / - r- Je{ [ r i 2 ] J i 4 n 2 ] „ ) m } Substituting (2.25) into 7^{[(iV 1) ( i_ 1 ) M + ,] T[(/V 1) ( m_ 1 ) M + n]} yields ^{[(JVl)(,--l)M +l]T[(^l)(m-l)M+n]} E E m ^ M r m u M ... E i ^ ^ f m u M Hence (2.27) / r f ? { [ ( J V l ) ( . - - l , « + l ] r ( ( J V l ) ( m - 1 ) M + n ] } / - [A7 • • • fL\ X ' ^{[(^i)i.i]r[(^i)».iU} ••• EiKN^JKNj^Ml} . « { [ (^U in ( iV lU} ^{[(^l).^]n(^)m^]n> E E /,T^{[(^l).,jri(^l)m,,]n}/ a * 1 M (2.28) Chapter 2. Model of Random Distortion 20 where / , and / , are the qth and the sth columns of F, respectively. Furthermore, the term f t m N i k M N i U M f , has a structure similar to that of / T £ { [ ( J V 1 ) ( i . l M # + i ] T [ ( J V 1 ) ( m _ 1 ^ + n l } / . if we regard subvectors fq and [(iVx)»,,Ji as a single unit. Hence, following a similar procedure as from (2.26) to (2.28), / / , ^{[(JVi ) < I , ] n( J v i )m, . ]n}/ . can be written as ffE{[W^]f[WmAn}f. = [fl,q • • • fM>q) X fl,3 ^{[(JVl)«,,U[(JVl)m,.]n,l} • fM,9 = E E / p , ^ { [ ( ^ l ) , , , ] « , P [ ( i V l ) m , , ] n , r } / r , . (2.29) r p Substituting (2.28) and (2.29) into (2.26) gives the explicit expression of [(Rn{f))i,m]i,n, {(Rn(f))i,m)l,n = E E E E ^ ^ O M U ^ W V } / , , / , r a p q +E{[n2]i,i{n2}nim} (2.30) where the notation [A.jjfc.j indicates the (A;,/)th entry in the (t,j)th submatrix of a general matrix A which can represent any of the matrices mentioned in (2 .30) . As for Rn2 the correlation matrix of the detection noise n2, we have i,m]',n = ^ { [ « 2]z , i ( « 2 ] n , m } ( 2 - 3 1 ) The matrices Rn2 and Rn{f) may therefore be computed from knowledge of the corre-lation functions of the noises n i ( x , y ) and n2(x,y). Chapter 2. Model of Random Distortion 21 Since Ni arises from a stationary process, we have the following relations E{[{N1)iigW1A(NlU,\n+l,r} = ^{[(^l)^ll,[(^l)m,.l»,r} Thus the matrix Rn{f) has a block Toeplitz structure. A simple but useful case is that when n2 is zero mean white noise. Hence (2.20) and (2.23) become E(n) = 0 (2.32) and Rnif) = E[(N1f)(N1f)T] + o\l (2.33) respectively, where o2 is the standard deviation of n2, and I is the identity matrix. For simplicity, define C ( , , ' m n ) = ^ { [ ( ^ O J - D M ^ I K ^ I ) ^ - ! ) ^ ] } (2.34) Ni may also be a white noise process. The latter arises wh^n estimating the blur function. The estimate suffers from measurement errors which may be modeled by white noise. For discussion on identification of blur function please refer to [37]. When Ni is white, the matrices C^' ' , m n ^ have elements where C\ is the standard deviation of Ni and [6ij]u is the Dirac Delta function 1 if i = j and k = I (2.35) 0 otherwise Chapter 2. Model of Random Distortion 22 Hence the tjth entry of the correlation matrix Rn{f) is [{Rn{f))i,m]l,n = ^iEE fp,qfm-i+p,n-l+q + °\bx-m,l-n (2.36) P « Another useful quantity in image restoration is the covariance matrix of n, Kn. When the point spread function is random, Kn is also dependent on / . Therefore we would rather write it as Kn(f). Under the assumption that n is of zero mean, the covariance matrix Kn(f) is equal to the correlation matrix Rn{/)- Hence when both N\ and ri2 have white statistics, the tjth entry of Kn(f) is also given by [{Kn[f))i,m]l,n - E E fp,qfm-i+p,n-l+q + ^l^i-m.l-n (2.37) P 1 which is the same as equation (2.26). 2.4.2 Random image field It is often convenient to regard an image as a sample of a stochastic process. If the image / is assumed to be a sample of a stochastic process, the correlation matrix of the noise is somehow different from that of the assumption that the object image is a deterministic signal. We shall start from the second equation in (2.26) [{Rn{f))i,m]i„ = ^{ / r [ (^i ) ( , - i )A / + / ] T [ ( iVi) ( m - i )M + n/ ]} + £ { [ n 2 ] M [ n 2 ] „ , m } (2.38) The first term in (2.38) can be expanded to ^ { / T [ ( ^ 1 ) ( i - 1 ) M + < ] T [ ( ^ l ) ( m - l ) M + n ] / } -I T E M m,l In [ ( J V l W l M m . l l n = ^ E E / J [ ( ^ ) « l [ ( J V i U / . } s q = EEE{fqT{(N1)i,q)r[(Nx)mAM [{Nl)iM]f[{Nl)mM]n fl fhl (2.39) » 9 Chapter 2. Model of Random Distortion 23 Furthermore, the term E{fqT\{Ni)iiq]f[{Ni)miS]nfs} is expanded to fl,q T • [(^'l).-, < r]l,l[(JVl)m,.]n 1M / l , a fM,q [(^l),- , ,]j^[(JVl)m,.]n,l • • • [(JVl).- , ,] l^[(JVl)m,.]n,M / M , » f P = E E ^ { [ ( ^ ) ^ ] « , [ ( ^ l ) - n , . ] - i , r / r , . / r , . } .(2.40) f P Since Ni and / are uncorrelated, we have ^ { [ ( ^ l ) ^ ] l , [ ( J V l ) « . . ] » . r / p r f / r . . } = ^ { [ ( ^ 1 ) ^ ] l , [ ( ^ l ) « , . ] » . r } J E ? t / p ^ / r > . ] (2-41) Substituting (2.39) and (2.40) into (2.38) gives the explicit expression of \(Rn(f))iim\iin, [(Rn(f))iim]iin = EEEE-5{|(-Viy.*[(-ViWvWM/rj r s p q } (2-42) By a similar argument as in the case of deterministic image field, we can show that the matrix Rn{f) is block Toeplitz. For the case when n x and n2 are uncorrelated white noise, an entry in Rn{f) is simplified to [{Rn{f))i,m]i,n = o\{M - \m - t | ) ( M - |n - l\)[{Rf)i,m]i,n + ^ ( 2 . 4 3 ) The covariance matrix here is also the same as the correlation matrix when the mean of the noise E[n] = 0. Chapter 2. Model of Random Distortion 24 2.5 Summary In this chapter the modeling scheme proposed by Ward and Saleh was generalized to the two-dimensional case. By means of this model we have converted the problem of restoration of images distorted by a system of random point spread functions to a problem similar to the conventional restoration. The random blur function is separated in two parts; the mean of the blur which acts as the deterministic blur function, and a zero mean random fluctuation around the mean which, together with the additive noise, and the object image, forms an image dependent noise. The expected value, and the correlation and the covariance matrices of the image dependent noise are obtained for various feasible situations. The material considered in this chapter establishes the foundation for the following chapters. Chapter 3 Restorations by the Wiener And the M V U Techniques 3.1 Introduction Linear filters have made great contributions in image restoration of the conventional class because of their simple formulation and economical computational efforts. Linear filters have also been introduced to the random blur problem [18-28,31]. Among them, Ward and Saleh [25] introduced restoration methods to the one dimensional problem, based on modifications of the Wiener and the minimum variance unbiased (MVU) esti-mation methodologies. Iterative schemes were adopted. The statistics of the unknown signals were assumed not to be completely known. Ward and Saleh [28] also discussed the same problem via the Backus-Gilbert technique. But there, the randomness in the blur was confined in its width to that of the deterministic part of the blur. In this chapter we extend the modified Wiener and M V U "filters [25] to the two-dimensional case. The noise in the blur here does not have to be confined in its width to that of the deterministic part of the blur. Since image restoration necessitates solv-ing huge systems of linear equations, whose computation in the space domain is pro-hibitively expensive, we shall have to resort to different methods of solution. We shall use circulant matrix approximation and the fast Fourier transform (FFT) , whenever possible. The computation of the Wiener filter in the Fourier domain using circulant approximation is a straightforward matter. That is one of the main reasons for the popularity of that filter in image restoration. The main drawback of the Wiener filter 25 Chapter 3. Restorations by the Wiener And the MVU Techniques 26 is the assumption that the correlation matrix of the unknown object image (determined by averaging over the ensemble of the objects) is known. In most practical cases, this is unrealistic, and one resorts to an approximation which hopefully might lead to the correct restoration. The M V U estimator does not require such knowledge, but its computation in the space domain is also infeasible cost-wise. Unfortunately using the circulant matrix approximation and the Fourier domain analysis in the same manner as in the Wiener filter does not help since this leads to the simple but not so desired inverse filter. Here we shall design more elaborate schemes so that the M V U estimator is realized. The solution will be partly in the Fourier domain and partly in the space domain. The cost is still much higher than that of the Wiener filter, but this at least makes the estimation feasible. Experimental results indicate that our proposed iterative Wiener and M V U filters are very worthwhile in that they outperform the standard Wiener and M V U filter (based on neglecting the randomness in the point spread function). We compare the performance of the iterative Wiener filter with that of the Backus-Gilbert technique [27]. Our findings are that the visual qualities of the restorations by the two filters are the same but the Wiener filter gives smaller root mean squared irrors than those based on the Backus-Gilbert technique. We also find experimentally that the MVU-based filter provides better results than the Wiener-based filter when the image is small and the distortion is not severe. Otherwise the M V U based filter may run into convergence problem. Therefore the practical choice between the two is the Wiener-based filter since it is computationally robust and inexpensive. Chapter 3. Restorations by the Wiener And the MVU Techniques 27 3.2 Restoration by the Wiener Filter The Wiener filter is one of the earliest methods that was applied to image restora-tion. The Wiener filter also finds extensive applications in control, signal processing, communication, and many other science and engineering fields. The reason for the popularity of the Wiener filter is due to its simple but powerful criterion and also because its solution has closed form. Using this method to restore randomly blurred images in the continuous domain was introduced by Slepian [19]. Ward and Saleh proposed discrete implementation of the Wiener filter [25] to the restoration of the one-dimensional randomly distorted signals. In this section, we generalize the work in [25] to the two-dimensional image restoration problem. 3.2.1 Filter Discussion The Wiener filter method is based on finding the linear filter L such that: f = Lg (3.44) The filter L is selected such that the mean squared error e2 between the unknown object / and the restored image / is minimized, i.e., mine 2 = wnE{\\f - f\\2} (3.45) JU Lt = mmE{Tr[(f-Lg)(f-Lg)T)} where Tr[A] indicates the trace of matrix A. Substituting g — Hf + n into (3.45) and manipulating the results yields e2 = E{Tr[(f-L(Hf + n))(f-L(Hf + n)f}} (3.46) = [{Tr[ffT - L(HffT + nfT) - (ffTHT + fnT)LT + L(H ffTHT + nfTHT + HfnT + nnT)LT}} Chapter 3. Restorations by the Wiener And the MVU Techniques 28 Notice that Tr[A] — Tr[AT] and since the trace is linear, it interchanges with the expectation operator E[ ]. Then e2 = Tr{Rf - 2LHRf + LHRfHTLT + LRnLT} (3.47) In the conventional case the noise is signal independent, and has zero mean. Differenti-ating (3.47) with respect to L and setting the result equal to zero, the mean-square-error solution is L — R/HT[HRfHT + Rn}'1 (3.48) In the random blur case the following equation E[n] = EliNif + ni)} = 0 (3.49) is true since in Chapter 2 we assumed that Ni is a zero mean process, and the measure-ment noise n2 has zero mean in general. Hence, there is no difference in the derivation of the filter, but the dependence of the noise on the image / should be specified by denoting Rn by Rn(f)- Thus the modified Wiener filter is L = RfHT[HR/HT + Rn(f)}-1 (3.50) where Rn(f) = E[nnT) = E\{NJ + n2)(NJ + n2)T] (3.51) Here equation (2.42) should be used to represent Rn{f) since / is assumed to belong to a stochastic process. If the noise sources Ni and n2 are white and uncorrelated with each other, the entries in Rn{f) is given by (2.43) 4 {Rn{f)i,m]i,m = {M - \m - i\){M - \n - l\)0l[{Rf)i,m\i,n + o\6i-.l<n-n (3.52) The major difficulty with the Wiener filter is that Rf is usually not known. At-tempts to circumvent this difficulty using a minmax approach have been advanced [22]. Chapter 3. Restorations by the Wiener And the MVU Techniques 29 Another approach is to find an estimate R{f) of Rf based on the estimated / , e.g., by using the scheme described in [33], {{R{f))i,m]l,n = T-pY1^2fp,gfp+rn-i,q+n-l (3.53) P I Because Rf is needed to compute / , an iterative algorithm is necessary. Using (3.44), (3.50) and (3.53), we write / r + 1 = R(fr)HT[HR(fr)HT + Rn(fr))-lg (3.54) where fr and R(fr) are the estimates of / and Rf in the rth iteration. From now on the subscripts indicating iterations shall be dropped. This is for the purpose of clarity of presentation. As mentioned in Chapter 1, space domain computation of (3.54) is very expensive, but since every term on the right-hand-side is block Toeplitz, we can approximate those matrices by circulant matrices, and the filter L (see (3.50)) is approximated by L « R f c H j [ H c R f c H j + J U / ) ] " 1 (3-55) Equation (3.55) is the circulant approximation to the Wiener filter. In (3.55) Rfc, Rnc{f), and H c are the circulant approximations of Rf, Rn{f), and H, respectively. The resulting system can be efficiently computed by F F T and iteratively. Since all circulants have the same set of eigenvectors, then it is obvious that sums, products, and inverses of circulants are also circulants. Equation (3.55) may be rewritten as: L « TM2S.fk*H\kHkfk*H + An]"1 J~\ (3.56) where A/ , Ah, and A n are the eigenvalue matrices of Rfc, Hc, and Rnc{f), respec-tively, IM2 is the two-dimensional D F T matrix, and * denotes complex conjugate. In Chapter 3. Restorations by the Wiener And the MVU Techniques 30 equation(3.56), the fact was used that if C is a block circulant with eigenvalue matrix A then A = IM* C 7Ml which implies A* = TM^C7!^ Equations (3.56) and (3.54) give: = A/A# [AjjA/A# + b.n\~l?M\g (3.57) Since equation (3.57) is a relation in diagonal matrices and JM*1 a n < i ^Af*9 a r e recognizable as the two-dimensional DFT's of / and g, respectively, then the matrix equation (3.56) can be written in terms of Fourier quotients: L(WJ,W*) = H ^ ' ' ^ , r (3.58) where Wn(u>i,u>k) is the power spectrum of Rn(f), and Wy(u>i,cn/fc) is the power spec-trum of Rj, H(wi,u>k) is the D F T of H, and * indicates complex conjugate. We also calculate the Fourier transform of the distorted image g, G(u/|,w f c), and multiply (3.58) by G(u>j,u;fc). Finally, we do the inverse Fourier transform to obtain the restored image fik To be more general, the Wiener criterion used for the two-dimensional restoration is extended so as to become parametric. Thus we need to add a parameter in the standard Wiener filter (3.50) in order to optimize the restored images. To be explicit, we rewrite (3.50) in its parametric form. L = RfHT[HRfHT + ^Rn(f)]~1 (3.60) Chapter 3. Restorations by the Wiener And the MVU Techniques 31 where 7 is the optimum parameter determined experimentally. Taking the discrete Fourier transform on both sides of (3.60), we obtain L ( W l , W f c ) = , (3.61) The inverse Fourier transform yields the restored image / M J |/r(W||W4)|» + 7 { ^ 3 S (3.62) 3.2.2 E x a m p l e s In order to test the two-dimensional Wiener filter, we use an image Face, which consists of 120 x 120 pixels, shown in Figure 3.1. The numerical values for the pixels are from 0 to 255. A random uncorrected blur function was generated by superimposing white Gaussian noise, Nlt on the deterministic sine2 shaped blur H. The sine2 blur function is known as the impulse response of an ideal diffraction-limited incoherent imaging system of a square aperture [37], h(x,y) = S i n 2 [ i r x ] 3 i n 2 { w y ] [irx]2 [ny]2 The standard deviation of the noise Ni is o\. The shape of the random blur is shown in Figure 3.2, where 3.2 a) is the deterministic blur and its noisy counterpart is 3.2 b). The additive noise n 2 is assumed white Gaussian, with standard deviation <T2, and is uncorrected with Ni, the randomness in the blur. The distorted image is given in Figure 3.3 for which C\ = 0.0005, and a 2 = 1.25, Although o\, the strength of the noise Ni, is small, its contribution to the noise in g can be very large, due to the relationship shown in equations (3.52) and (3.53). Chapter 3. Restorations by the Wiener And the MVU Techniques 32 A straightforward linear Wiener estimate, based on neglecting the effect of the impulse response noise TVi, is shown in Figure 3.4. The magnification of noise by the linear Wiener filter is obvious. This estimate serves as an initial guess for the iterative approach of Eq.(3.54), in which both noise sources are included. We first used the modified Wiener filter (3.59), that is 7 = 1., then the parametric Wiener filter (3.62) to do the restoration. Experiments show that the restoration is optimized when 7 is equal to 2. See Figure 3.5 a) and 3.5 b). The parametric Wiener filter is seen to result in a smoother picture. We can see that the linear Wiener elimi-nates the blur effect but, the restoration quality is obscured by oscillating noise. On the other hand, the iterative filters give sharp and smooth restorations. For quantitative measurement, the root mean squares (RMS) errors of each image in the example are listed in Table 3.1. Table 3.1: RMS errors of the Wiener technique RMS Blurred image 14.49 Linear Wiener 33.81 Modified Wiener 11.29 Parametric Wiener 10.34 Chapter 3. Restorations by the Wiener And the MVU Techniques 33 Chapter 3. Restorations by the Wiener And the MVU Techniques 34 Figure 3.3: Distorted version of Figure 3.1 • • Figure 3.4: Restoration of Figure 3.3 by standard Wiener filter Chapter 3. Restorations by the Wiener And the MVU Techniques 35 a) b) Figure 3.5: Restoration of Figure 3.3 by: a) the modified Wiener filter, b) the modified parametric Wiener filter We also compared the performance of the modified Wiener filter and the Backus-Gilbert technique [28]. The results are shown in Figures 3.6 and 3.7. Figure 3.6 is the random blur used for this comparison. According to the assumptions of the Backus-Gilbert technique, the noise is confined to the width of the blur function. Figure 3.7 a) is the same as Figure 3.1, but with a smaller scale. The blurred image is shown in Figure 3.7 b) with the strengths of the noise components being 0\ = .006 and er2 = 3., respectively. (Since the effect of the noise in the blur function is magnified by the convolution with / and since in this case ni is confined to the width of the blur, bigger 0\ and a 2 can be selected). Figure 3.7 c) and 3.7 d) are the results of the modified iterative Wiener and the technique in [28], respectively. It is found that the modified Wiener filter gives smaller RMS errors than the Backus-Gilbert technique (see Table 3.2), but the visual qualities of the restorations by the two methods look the same. Chapter 3. Restorations by the Wiener And the MVU Techniques Figure 3.6: Random blur: noise confined to the width of the blur function Table 3.2: RMS errors of the B - G technique RMS Blurred image 11.67 Modified Wiener 6.54 B-G technique 7.63 Chapter 3. Restorations by the Wiener And the MVU Techniques 37 Figure 3.7: Comparison between the modified Wiener filter and the Backus-Gilbert technique: a) The original image; b) The blurred image; c) Restoration by the Wiener-based filter; d) Restoration by the Backus-Gilbert technique. Chapter 3. Restorations by the Wiener And the MVU Techniques 38 3.3 Restoration by the M V U Filter When the unknown / is fixed, the solution of the minimum variance unbiased method is equivalent to the solution of the weighted least-squared estimate with the weight matrix being the covariance matrix of the noise. The rationale of the weighted least-squared criterion is to select a solution that is consistent with the assumption that a function of the noise term n is as small as possible given the data g. The rationale behind the M V U criterion is that the estimated signal has the same mean as the original signal and has the smallest variance amongst all such estimates. 3.3.1 The Filter Derivation The M V U filter seeks the solution min nTn = {g - Hf)TKn{g - Hf) (3.63) where Kn is the covariance matrix of n. The solution of (3.63) is derived by differen-A A tiation with respect to / and solving for / . For the usual image restoration problem wherein the covariance matrix Kn is known a priori the M V U estimator is J = [HTK-lH\-lHTK~lg (3.64) When Kn is diagonal with equal elements, the estimator becomes the pseudo-inverse filter / = H~1g (3.65) For our problem Kn(f) depends on the unknown / . If Kn(f) is unknown a priori, then equation (3.64) is no longer the M V U estimate, and determination of the M V U filter is very complicated. Actually the determination of the M V U filter in this case is quite similar to the M A P filter which we shall discuss in the next chapter. Chapter 3. Restorations by the Wiener And the MVU Techniques 39 For our problem, when the noise sources nx and n2 are white, the entries in /^(/Jbecomes (2.36) [Kr,(/),,m]j,n = <J\ E E fp.qfp-i+m.q-l+n + o\f>i-m,l-n P « Assuming Kn(f) is known a priori, then the M V U estimate is f^\HTK-\f)H]-'HTK-\f)g (3.66) to solve (3.66) we employ iterative technique. We select initial guess / (°) for the un-known image / , from which we compute a estimate for the covariance matrix J F S T „ ( / ( 0 ^ ) . With the covariance matrix known, we subsequently use equation (3.66) to determine the M V U estimate plK We use that estimate to compute a new covariance matrix Kn(fW), and so on. The scheme follows the following algorithm: f{m+i) = \HTK-\f^)H\-lHTK-l0^)g m = 0,1,... (3.67) The initial guess we used is the M V U estimate based on neglecting the /-dependent part of the covariance matrix Kn{f), i.e., using instead Kn2, of which f^ = [HTK-21H]-1HTK-21g ' (3.68) If the algorithm converges, then the result converges to a value / that satisfies the nonlinear equation f = [HTK-\J)H\^HTK-\J)g (3.69) Although the mathematical derivation of the M V U filter under the above conditions is not hard, its application to the two-dimensional signal, or image restoration requires more elaborate schemes than the Wiener filter. The problem with the M V U is that, if we apply F F T with circulant matrix approximation, the resulting filter is the inverse filter [15], see Appendix A for proof. This is not what we want. Thus the proposed Chapter 3. Restorations by the Wiener And the MVU Techniques 40 solution will be basically in the space domain. However, since the statistical assumption of the problem enables us to use circulant approximation and F F T , these shall be exploited to our advantage only when the basic information is not canceled out by their use. Savings on computations, though not as significant as with the Wiener filter, are still achieved. The procedure to compute the MVU filter involves a combination of space domain and Fourier domain manipulation. To begin with, we rewrite equation (3.69) in its equivalent form HT K~1(f)Hf = HT K~1(f)g (3.70) Next we define T as T = HTK~\f)H (3.71) and g as g = HTK-l{f)g (3.72) Then (3.70) becomes Tf = g. (3.73) Because both H and Kn(f) can be approximated as circulant,-The F F T could be ap-plied to solve (3.71), (3.72) and (3.73). But this results in the inverse filter and the useful information carried by Kn(f) is canceled. Instead, F F T and circulant approx-imations shall only be used to invert Kn(f) and compute approximations for (3.71) and (3.72). (See Appendix A for details). Then (3.73) is solved in the space domain. Equations (3.71) and (3.72) could also be formed in the space domain, but it is more efficient and does not constitute much loss in information to form them in the frequency domain. The complexity of inverting Kn(f) and forming (3.71) and (3.72) is reduced significantly, from 0 ( M 6 ) to 0 (M 2 l ogM) . The Appendix A shows that the resulting system of linear equations is the approximation to (3.73), which is block Toeplitz and Chapter 3. Restorations by the Wiener And the MVU Techniques 41 contains all the data useful to compute / . To solve this system, the Toeplitz system solver designed by W . F . Trench [34] and extended to the block Toeplitz system by H. Akaike [35] can be applied. The computation of this method is 0 ( M 6 ) instead of 0 ( M 6 ) . 3.3.2 Examples Due to the limitation of facilities and the complexity of computations, 0 ( M 6 ) for an M x M image, the M V U estimator was applied to the small 32 x 32 image, Number, shown in Figure 3.8 a). A Gaussian shaped blur function was used. -fx 2 + v2) H(x, y) = Cexp \ , V I (3.74) where C is a constant for scaling. The image Number was distorted by uncorrelated white Gaussian noise sources as discussed in the example of Section 3.2.2, but with Oi = 0.0005, and o2 = 2.55. Figure 3.8 b) gives the distorted image. Restorations by the modified iterative M V U and the Wiener filters are presented in Figure 3.8 c) and 3.8 d), respectively. The RMS errors for this example are fisted in Table 3.3. The Wiener approach gives slightly better result than the M V U in this example. Table 3.3: RMS errors of the M V U filter RMS Blurred image 19.49 M V U 16.49 Modified Wiener 15.78 ' - " - J - -For comparisons on the RMS for am 10 x 10 image, among linear Wiener, and the modified iterative Wiener and M V U filters, when one of the parameters, an, a l 5 and d, is changing, see Figures 3.9 - 3.11. Figure 3.12 shows the RMS vs. M , (the size of image processed is M x M ) , for Wiener and M V U filter. The results demonstrate that the iterative M V U filter outperforms the iterative Wiener filter on small images. For larger images though, there are indications that the MVU-based filter may run into convergence problems. This is not true of the Wiener-based filter. 01568 a) b) 0IS68 InBB^fl^B^fl^B^BHfl 1 • 1 1 H J - 1 B otsee c) d) Figure 3.8: Distortion and restorations of the image Number: a) The original image; b) The blurred image; c) Restoration by the MVU-based filter; d) Restoration by the Wiener-based filter. Chapter 3. Restorations by the Wiener And the MVU Techniques 43 Figure 3.10: Filter performance as ot changes Chapter 3. Restorations by the Wiener And the MVU Techniques 44 Figure 3.12: Filter performance as the image size changes Chapter 3. Restorations by the Wiener And the MVU Techniques 45 3.4 Some Remarks In the process of restoring randomly distorted images, some technical details should be considered carefully. Firstly, in order to better approximate the noise correlation matrix Rn{f) and the covariance matrix K„(f) by circulant matrices, we need to take into account the statistical nature of the noise. In the case that the noise is image independent, it possesses no predictable correlation beyond a distance d. That is , \i-j\>d (3.75) Typically d is about 20-30 pixels, see [15]. But when the noise is image dependent, the distance d is related to the size of the image under consideration. If the distance d is not properly chosen, restoration may give incorrect results. If the distance chosen is too small, the approximated Rn{f) tends to become an identity matrix, and useful information is lost. On the other hand, if the distance chosen is too large, some artifacts are introduced, which may result in an unstable solution. We determined the distance d by experiment. We found that d is dependent on the size of the image, and satisfies d« (0.05 ~ 0.2)M ' (3.76) for a matrix of M x M. The above discussion is also valid for the covariance matrix. The second technical consideration involved is in computing the M V U estimator. From the discussion in Section 3.3.1, we know that the computation of Tf = g is still expensive, 0 ( M 5 ) . Although there exist some better tools to solve the system, e.g. the Toeplitz solver introduced by R.R. Bitmead [36], with complexity of 0 (M 4 l ogM) , more powerful algorithms need to be explored. Otherwise, the M V U estimator remains a noneconomical solution. Finally, we discuss the criteria to judge the restoration quality. Visual judgement Chapter 3. Restorations by the Wiener And the MVU Techniques 46 is important in image processing. But the criteria for visual judgement are very sub-jective. Viewers vary in their preference between smooth images at the expense of resolution and high resolution images regardless of noise. The criterion we choose is to seek restorations having reasonably good resolution without much noise, i.e. a trade off between high resolution and noise reduction. The judges of visual quality are mainly myself and the research supervisor. When in doubt, fellow students were asked to speak out their judgement. The students were told of criterion for our judgement. We also use the RMS errors as quantitative measurements. The visual judgements and the RMS measurements have been used in this chapter, and will be used throughout this thesis. Sometimes, we need more detailed information to judge the quality, e.g., the last example in Chapter 5, and the examples in Chapter 6. Then we will show the difference pictures. 3.5 S u m m a r y For restoration of images distorted by systems of random impulse response and additive noise, the Wiener filter and the M V U estimator are modified, extended and applied to the two-dimensional case. The randomness in the blur is not assumed to be necessarily confined to the width of the impulse response. It is shown that the algorithms for computing the modified Wiener filter can be implemented in a straight-forward manner using circulant approximation and the F F T . On the other hand, in order to make the computation of the expensive M V U estimator feasible, a combination of space domain and frequency domain manipulations has to be appropriately applied. Although the scheme for the M V U estimator is developed for the more general prob-lem (random blur with additive noise), it is also applicable to the classical restoration Chapter 3. Restorations by the Wiener And the MVU Techniques 47 problem (deterministic blur with additive noise). The simulation examples presented confirm that these filters work well for the restoration of the random blur problems. The computationally expensive MVU-based estimator was found to do better than the Wiener-based filter for smallimages, when the width of the point spread function is not too high. Thus for small images the M V U -based estimator is recommended. Otherwise the MVU-based filter may run into con-vergence problems. In which case, the very robust Wiener-based filter is recommended. The newly derived filters always give better results than their linear counterparts (based on ignoring the noise in the impulse noise). The modified parametric Wiener filter gives pleasing (smoother) images than the standard modified Wiener filter. The latter also gives better results than those based on the Backus-Gilbert technique [28]. Chapter 4 The Constrained Least Squares and the Geometrical Mean Filters 4.1 Introduction In this chapter we discuss the application of the constrained least squares technique and the geometrical mean methodology to the random distortion problem. The reason for choosing to apply the constrained least squares filter to our case is that we found that the restoration results of the modified Wiener and the M V U filters tend to be noisy when the blur effect and/or the noise effect become severe. The constrained least squares technique can be formulated so as to remedy that. Another advantage of using this method is that it allows the inclusion of a variety of constraints. The constrained least squares filter can be implemented in the Fourier domain. Hence the computation of this technique is inexpensive. We applied this filter to our restoration problem. It is found that this filter does give smooth solutions. As we shall see in the following discussion, although the constrained least squares fil-ter solves the problem of noisy restoration, satisfactory results have only been achieved by manual adjustment of a parameter. The geometrical mean methodology is also based on the adjustment of parameters. Since the Wiener-based filter recovers the blur effect well but yields a noisy solution while the constrained least squares filter gives a smooth result, a combination of these two filters should result in a restoration better than either. Therefore it is worthwhile trying the geometrical mean approach based on the Wiener and the constrained least squares techniques. 48 Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 49 4.2 Restoration by the Constrained Least Squares Estimation The constrained least squares estimation was first introduced to digital image restoration by B.R. Hunt [8]. The reason behind the introduction of this technique is to overcome the inherent ill conditioning in the image restoration problem. Numeri-cal solutions to the integral equations of the first kind, equation (1.4), are obscured by wildly oscillating noise. The gross oscillations are usually recognized as being not con-sistent with a priori knowledge that the solution is smooth (or at least smoother than the noise oscillations actually obtained). Therefore, the solution should encompass a smoothness measure of some kind. The smoothness measure Hunt chose was the inner product of the second difference of the solution vector, which is quantitatively sensitive to the relative smoothness of the solution vector. 4.2.1 The Filter Derivation The derivation of the constrained least squares filter involves minimization of some linear operator applied on the object / subject to the constraint that the sum of the noise residuals is known. Mathematically the problem is written as mm(Cf)T(Cf) (4.77) subject to (g - Hf)T(g - Hf) = nTn (4.78) In (4.78) the value nTn is assumed known or measurable a posteriori. Depending on the requirements of the quality of the restoration the linear operator C in (4.77) can take a variety of forms. For instance, Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 50 • C — I. i.e., a minimum norm solution on / . The answer to the restoration problem here subject to the noise norm constraint is derived to be: / = [HTH + !l]HTg (4.79) where 7 will be explained later. C = [finite difference matrix]. It may be chosen to minimize the second (or higher-order) difference energy of the estimated object. In this case, for the second difference, C = Cl®C1 (4.80) with -2 1 1 - 2 1 0 1 -2 0 0 1 - 2 1 1 -2" which is a tridiagonal matrix. This operator results in a smooth picture. Such operator constraint guarantees that the object estimate / does not oscillate wildly in the constrained solution. • C = [eye model]. One may desire that the restoration be appealing to a human from a perceptual viewpoint. In this case C is a block circulant whose p rope r t i e s in the Fourier domain match the spatial frequency response of the p sychophys i c s of the human visual system [38]. • C = R-}ll2R-1'2. This case results in a parametric Wiener filter. Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 51 Since our motivation of using the constrained least squares is to eliminate the noisy phenomena in the solution, we choose the second difference matrix as our linear oper-ator. Using the method of Lagrange multipliers, the unconstrained objective function is W{f) = (CffiCf) + A[(f7 - Hf)T(g - Hf) - nTn) (4.81) Taking derivatives of the objective function with respect to / and setting the result equal to zero, = 2[C Cf — XH (g — Hf)] = 0 (4.82) df Solving for the / that provides the minimum of the objective function yields / = {HTH + ^CTC)-lHTg (4.83) where 7 = 1/A. Termination of the algorithm is governed by satisfying condition (4.78). Because the right hand side of (4.78) is not available, an approximation must be selected. A reasonable approximation is derived from an unbiased estimate of the variance of n [8] 2 nTn - p,2 & = M 2 - (4.84) where M = ^ £ I > ; (4-85) is the unbiased estimate of the mean of the noise n. Therefore the noise term can be estimated as nTn = M2d2 + fj,2 (4.86) Hunt proved that the satisfaction of equation (4.78) can be achieved automatically using an optimization routine [8]. This is because the function 4>h) = {g-H})T{g-Hf) (4.87) Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 52 is monotonically increasing as a function of 7, when 7 > 0. When this method is applied to our random blur restoration problem, the solution (4.83) will differ. General stationarity assumptions on noise statistics are not relevant to us here. We restrict the noise Nt and n 2 to the class of the white random process. Differentiating (4.81) with respect to / and considering that n = Nif + n 2 is function of / result in = 2[CTCf - XHT(g - Hf) - XN^NJ + n2)] = 0 (4.88) of Hence the estimated image is / = [HTH - NfNt + 7 C r C ] - 1 [ t f T f f + N?n3] (4.89) Since Ni and n 2 are noise, it is impossible to obtain the exact solution to the constrained least squares method. However, the statistics of the noise which are assumed known can be utilized to compute an approximate solution. We first consider the term iNrfn2. We know that Ni is an matrix of M2 x M2 and n2 is an vector of length M 2 . The tth element in N?n2 is \N?n3]i = Y l W M i « 0 - (4-90) since Ni and n2 are uncorrelated processes, and both have zero mean. Next we consider the term N[N\. The tjth element in the matrix product is Since Ni is block Toeplitz and the elements in N\ are assumed to be generated from a uncorrelated white process, equation (4.91) is approximately equal to zero unless i — j. When i = j, the element can be estimated following a procedure similar to equations (4.84) to (4.86). Thus [NfN^u » o\M2 (4.92) Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 53 In equation (4.92) we discarded /x since Ni is known to have zero mean. The above discussion indicates that the matrix NfNi can be approximated by a diagonal matrix in which the diagonal entries are identical, and equal to o\M2. Thus NfNi w o2M2I (4.93) where J is the identity matrix. Substituting (4.90) and (4.93) into (4.89) yields f^[HTH~o\M2I + 1CTC}-1HTg (4.94) Since 7 is unknown in (4.94) we choose 7 so that the resultant / satisfies the constraint (4.78). To compute an estimate to the noise norm nTn for (4.78), we have mentioned that n has a zero mean (as a result of the statistics of Ni and n 2). Hence equation (4.86) becomes nTn = M2d-2 (4.95) where » , = [JRn(/)]n = a ;E .E / ' i + ^ (4-96) is the variance of noise n using the assumption of the deterministic object image indi-cated in Chapter 2. Since / is unknown in (4.96), the estimated / will be used instead. 4.2.2 The A l g o r i t h m and Experiment Results As with the Wiener filter, numerical solution of the constrained least squares filter in the space domain is extremely lengthy due to the size of the problem. Since C\ is a Toeplitz matrix, the structure of matrix C — Ci<g>Ci is block Toeplitz. We already know that H is a block Toeplitz matrix. Therefore C and H are approximated by circulant Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 54 matrices, and the computation is carried out in the Fourier domain. Restoration by the constrained least squares filter using the F F T is summarized in the following algorithm. Algorithm 1. Choose and fix an initial value for 7; 2. Compute the F F T of H, C, and g, then compute p(l, , , 1 _ H'{us,uv)G{ut,u)v) 1 " " j " | / f ( a ; . , W r ) | » - M » d ? + ' y | C ( a ; . , W i r ) | » ^ where * indicates complex conjugate. A A 3 . Compute the inverse F F T of F(u)x,u>v) to obtain / ; 4. Compute ^ ( 7 ) ; 4>{l) = {9- Hf)T(g - Hf). 5. Estimate nTn by (4.95); 6. If ^ ( 7 ) — nTn\ < a, where a is a predefined positive constant which determines accuracy, stop the iteration. Otherwise, use a zero finding algorithm, such as bisection procedure to improve 7. In his original work [8], Hunt suggested a Newton-Raphson algorithm to find the correct value of 7 such that equation (4.78) is satisfied since he proved that, in the con-ventional image restoration, the function ^(7) is monotonically increasing as a function of 7. Therefore there is only one unique value of 7 such that equation (4.78) is satis-fied. In our problem the condition of monotonicity may be violated due to the negative term —o\M2 in equation (4.94). There might exist more than one solution. Hence the solution derived by the bisection procedure might be a suboptimal solution. Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 55 The modified constrained least squares filter was tested on the image shown in Figure 3.1. The image is blurred by a random impulse response with a Gaussian shaped mean plus white additive noise. The distorted version is given in Figure 4.13 a). Here the standard deviation for noise sources n\ and n 2 are 0\ = 0.00075 and cr2 = 8.0, respectively. The restoration given by the modified Wiener filter is shown in Figure 4.13 b). The noisy effect in this restored image is very severe. The constrained least squares filter is also applied to this example. The result is given in Figure 4.13 c). Because of the second difference operator, C tends to smooth the picture when condition (4.78) is satisfied, the restored image looks very smooth, but resolution is lost. If we control 7 manually, resolution of the image is improved and the restoration does not look noisy, as shown in Figure 4.13 d). The RMS errors obtained in this experiment are summarized in Table 4.4. For completeness, the RMS errors by applying the conventional Wiener filter and the parametric Wiener filter are also included in the table. To further study the results, the difference images of the distortion and three restorations are given in Figure 4.13A. The observation of the difference images is agreeable with the discussion above. The advantage of using constrained least squares filter is that this approach is very simple. The only required statistic & is easy to obtain. Table 4.4: RMS errors of the constrained least squares filter RMS Blurred image 15.57 Conventional Wiener 116.47 Modified Wiener 17.82 Modified parametric Wiener 16.34 CLS: constraint satisfied 14.11 CLS: manual adjustment 11.46 Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 56 a) b) c) d) Figure 4.13: The performance of the constrained least squares filter: a) distortion; b) restoration by the Wiener-based filter; c) restoration by the constrained least squares filter: constraint satisfied; d) restoration by the constrained least squares filter: manual adjustment. Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters c) d) Figure 4.13A. The difference pictures: a) The difference between Figure 3.1 and 4.13 a); b) The difference between Figure 3.1 and 4.13 b); c) The difference between Figure 3.1 and 4.13 c); d) The difference between Figure 3.1 and 4.13.d -Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 57 4 . 3 The Geometrical Mean Filter The first geometrical mean approach to the restoration problem was put forward by E.R. Cole in his Ph.D. dissertation [7]. His aim was to parameterize the ratio of the effect on the restoration of the inverse filter to the Wiener filter. The motivation for such a parameterization is the desire to de-emphasize the low frequency dominance of the Wiener filter while avoiding the early singularity of the inverse filter. One such parameterization might yield an estimate of the object image as where 0 < a < 1 and the filter is defined in the scalar or diagonal spaces. For the random blur problem we have shown, in the previous discussions, that the modified Wiener filter gives a restoration which recovers the blur effect but somehow increases the noise effect, whereas the generalized constrained least squares filter has the opposite function. An intuitive solution to this problem is to design a filter similar to the geometrical mean filter shown in (4.98). Thus designing a filter B such that For convenience, we call this filter the geometrical mean filter. While the filter above may not be obvious in its current form, we can develop its Fourier diagonalized equivalent: f = [Inverse filter)01 [Parametric Wiener filter)1 ag (4.98) B = [Constrained least squares filter)01 [Wiener filter) l-a (4.99) (4.100) The algorithm to implement the geometrical mean filter is the following: 1. Use the algorithm in Section 4.2 to find a constrained least squares solution /(!), and the optimum value of 7 in the sense that condition (4.78) is satisfied. Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 58 2. Use the Wiener-based filter presented in Chapter 3 to yield another restoration /(2)-3. Use /(x) to compute b\ and /(2) to compute Rj and Rn{f)-4. Apply the discrete Fourier transform to H, C, Rf, and Rn{f) to obtain the required quantities in equation (4.100). Choose a a value. 5. Form B(ux,uy) by equation (4.100). Then compute the restoration in Fourier domain by F(ut,wv) = B{ux,uv)G{ux,wv) (4.101) the 7 value used here is the optimum value for the constrained least squares filter. 6. Use inverse Fourier transform to obtain a geometrical mean restoration. 7. If the restoration quality is unsatisfactory, change the a value and go through steps 5 and 6 again. In Cole's work there are two parameters, a, and 7 in the parametric Wiener filter, are both adjusted manually. Our approach here is different. W? first find the "best 7 " for the constrained least squares filter automatically, then find the parameter a in the geometrical mean filter manually. We test the obtained geometrical mean filter on the example given in the previous section. The restoration is given in Figure 4.14. The a value for this result is 0.35. The RMS error is 11.562, which is the same as the result by the manual adjustment of the constrained least squares filter. But this result does not have the dark patchs which can be seen in the restoration by the constrained least squares filter. Chapter 4. The Constrained Least Squares and the Geometrical Mean Filters 59 Figure 4.14: The restoration by the geometrical mean filter 4.4 Summary In this chapter we generalize the constrained least squares technique to the random restoration problem. We also considered the application of a geometrical mean method-ology to the addressed problem. The algorithms for computing the constrained least squares filter and the geometrical mean filter can be implemented in straightforward manners using circulant approximation and the F F T . The simulation examples presented confirm that the constrained least squares filter performs in a consistent manner with the o priori knowledge that the solution to restoration problem is smooth. B y manipulating a parameter, the geometrical mean filter provides satisfactory results. Chapter 5 Restoration by the Nonlinear MAP Filter 5.1 Introduction The modified iterative methods derived in Chapters 3 and 4 become linear filters when the noise in the blur is abscent, i.e. when applied to the "conventional" model. Linear methods have accounted for most of the image restoration solutions to the practical applications of real world problems. The reason is simple: Linear methods can be computed in a straightforward and economical fashion. Conversely, nonlinear meth-ods mostly require much more elaborate and costly computational procedures. In spite of this difficulty, nonlinear methods should not be overlooked, for two reasons. Firstly, the capabilities of computers and numerical algorithms are continuously increasing. Secondly, nonlinear methods must be considered for comprehensiveness and correct-ness; many nonlinear solutions result from a more concise and "accurate description of the image restoration problem. The necessity for nonlinear restoration methods arise in several ways. The first situation in which nonlinear image restoration may arise is due to inherent existing constraints in the problem. Secondly, the nonlinearity may directly come about from the nature of the image formation and recording system. The nonlinearities in the recording media are ignored or approximated in the discussion of linear filters. A third situation in which a nonlinear image restoration problem can arise is due 60 Chapter 5. Restoration by the Nonlinear MAP Filter 61 to the way the problem is formulated. Formulation of the problem using descriptions-mathematical or probabilistic-that do not simplify to linear forms creates a nonlinear problem in image restoration. Typically such formulation arise from the use of op-timization criteria that are more complex than the quadratic or the mean-squares criteria. In this chapter we discuss how to apply the traditional nonlinear restoration meth-ods to the random distortion problem. In our case the reason to obtain nonlinear filters is because the noise component n is a function of the object / . Therefore, using a non-linear restoration procedure which does not necessarily become linear when applied to the conventional problem is expected to reverse the distortion more correctly, and thus to result in a better restoration quality. Our aim here is to adopt a nonlinear statistical method: the maximum a posteriori (MAP) method to our problem. The M A P filter was originally used in the conventional image restoration problem so as to deal with sensor nonlinearities [15]. The conventional distortion model with sensor nonlinearities is given by g{x, y) = s{h{x, y) * /(x, y)} + n(x, y) ^ (5.102) where s{.} indicates the sensor nonlinearities. The scheme we attempt is based on the model given in Chapter 2. Since the noise in this model is related to the unknown / , the derived filter is more complicated than that of the conventional model with sensor nonlinearities. Chapter 5. Restoration by the Nonlinear MAP Filter 62 5.2 T h e Est imat ion Procedures In this section, the performance of the MAP estimator is investigated. Based on the Bayes' formula „(/,„) = 'MM!! (5. 1 03) where p(.) represents the probability density function, pdf, the MAP estimator finds A the estimated value / which satisfies p{f\g) = maxp(/|<7) (5.104) Since an arbitrary assumption about the unknown / and the noise Ni and n 2 results in an extremely complicated problem, we will assume that both the unknown / and the noise belong to Gaussian stochastic processes. We shall show later that the estimation algorithm derived under the Gaussian assumption works well on practical examples. The pdf of a Gaussian process / is given by P(f) = 1TT^T^^P{-\{f-l)TK]\f-1)} (5.105) where / and Kf, the mean vector and the covariance matrix c-f / , are assumed to be ensemble statistics. Finally, the conditional a priori probability p{g\f) can be derived. Since both Nx and n 2 are Gaussian with zero means, the conditional pdf of g given / is also Gaussian, and given by = , „ , M , t ( t ^ * P { - \ i 9 - Hf)TK-l{f){g - Hf)} (5.10G) (2-T) a 2 where Kn is the covariance matrix of n introduced in Chapter 2. For simplicity, the discussion about the MAP filter will be limited to the one di-mensional case; generalization to the two-dimensional case is straightforward. In the Chapter 5. Restoration by the Nonlinear MAP Filter 63 one-dimensional case the matrices H and N\ are both Toeplitz. It has been shown that the covariance matrix Kn{f) has Toeplitz structure and the elements in the correlation matrix is given by [25] [Kniflk^f A*'f+ (5.107) where the entries in matrices A*' are [A%m = E W i A N ^ } (5.108) Furthermore, if N\ and n 2 are white noise, and uncorrelated with each other, the elements of Kn(f) in (5.107) can be written explicitly, [*»(/)].•> = <>\ E M i - i + i + (5-109) i Substituting (5.105) and (5.106) into (5.103) results in l{9)\Kn{f)\> 2 + (g-Hf)TK-1(f)(g-Hf)}} where g(g) = (2n)M\Kf\^p(g) is a function of g. Because In x is a monotonic function in terms of x, and because In p(f\g) is much easier to deal with, the logarithm of (5.110) will be considered. For the sake of clarity, we will denote \np(f\g) by W(f). Taking the logarithm of (5.110) yields W(f) = \np(f\g) (5.111) = -Q.h{\n\Kn{f)\ + {f-J)TKj\f--f) + {g - Hf)TK-\f){g - Hf)} - \nq(g) Chapter 5. Restoration by the Nonlinear MAP Filter 64 Equation (5.111) is the objective function (of the M A P estimator), which is to be maximized. Unlike the modified Wiener and the M V U filters, which can be computed in a closed form, the nonlinear M A P filter can only be realized by an optimization algorithm. The algorithm we choose is based on the steepest ascent method. We first differentiate W{f) with respect to / to obtain the gradient vector V W ( / ) VW(f) = -Kj1[f-J)-0.5{ where Sin |gw(/)l , d\(g-Hf)TK^(f)(g-Hf)}] { df df ( JL a ah V 3 /M J (5.113) is the gradient with respect to / . To maximize W(f) in equation (5.111), gradient (5.112) can be computed, equated to zero, and solved to obtain a necessary condition for a maximum. If there exists a value of / such that equation (5.104) is satisfied, then the value is called the M A P estimate (for maximum a posteriori) and designated by the notation /MAP- Assuming that / M A P exists, then setting V W ( / ) equal to zero and rearranging (5.112) yields f - 7 os FT <dln\K»V)\ , d[{g - H f)T K-\f){g - H f)} /MAP = f - 0.5if/{ — + — }\f=fMAp (5-114) This is a nonlinear matrix equation for / M A P since / M A P appears on both sides. If we ignore the randomness in the point spread function, then the problem is simplified to the conventional form. Specifically, the statistical quantities are n = n 2 Kn{f) = Kn2 Chapter 5. Restoration by the Nonlinear MAP Filter 65 which are not object dependent any more. This fact leads to d\n\Kn(f)\ df = 0 (5.115) and (5.116) Therefore, equation (5.112) equal to zero has a explicit linear solution: W = [HTK^H + I f / 1 ] - 1 Jf r / f - a 1 S + [HTK~2lH + Kj^Kj1/ (5.117) Equation (5.117) is the same as the second form of the Wiener filter expression given in Appendix D provided This result is not surprising since it is known that in a conventional linear system, the M M S E and the M A P estimates assuming Gaussian statistics are the same [39]. In the case where the point spread function is random, solving the matrix equations (5.112) cannot be done in a closed form. As mentioned before we choose the iterative procedure, steepest ascent, to find the solution. The steepest ascent algorithm requires an explicit expression for the gradient VW{f). Therefore the next step is to find the two derivatives in equation (5.112). To make the notation shorter, we use V W i and VW2 to represent the two derivative items in equa-tion (5.112), i.e. .7 = 0. d\n\Kn{f)\ df (5.118) and VW2 = d[(g-Hf)TK-l(f)(g-Hf)} df (5.119) Chapter 5. Restoration by the Nonlinear MAP Filter 66 respectively. It is shown in Appendix B that / a[y„(/)],y \ Ofi (5.120) df, where Kn1(f) is the inverse of Kn(f). If NY is white and uncorrelated with n 2 , the derivatives in (5.120) are obtained by differentiating \Kn(f)]ij with respect to ft,I = 1,2 , . . . ,M, ' fi+i-i + fi-i+j, if l + i-j <M, I- i + j > 0; /,+<_;,' if l + i - j < M , l - i + j <0; (5.121) /,_,+,-, if / + i - j > M,I - i + j > 0; 0, otherwise. Since the entries of Kn(f) are functions of / , the entries in the inverse matrix of Kn(f) are also functions of / , and the derivation of VW2 involves differentiation of the inverse of Kn{f). Details of the computation is covered in Appendix C. We write the result in the following equation: \VWt\i = ^[^K-'img-Hf)], -{g - Hf)TK-\f)^P-K-\f){g - Hf) (5.122) ofi where the subscript / indicates the /th entry. Substituting (5.120) and (5.122) into (5.112), we obtain the /th entry of the gradient vector VW(f) |v^( / ) l , = - [ t f / M / - 7)]< - 0.5[VWV+W 2 ], = - [ ^ ( / - J i i i - o ^ D C f / l k (5 .123) d\Kn(f)h dfi +{HTK-\f){g - Hf)\t + 0.5(g ~ Hf)T^(f)^HK-l(f){g - Hf) ofx Chapter 5. Restoration by the Nonlinear MAP Filter 67 In the steepest ascent approach, the gradient vector VW(f) gives the current searching direction. If the fcth estimation has been accomplished, then the (k + l)th estimated value of / is computed by /fc+i = fk + <*kVW{fk) (5.124) where a* is a control factor to accelerate the convergent rate. A golden section search is selected to find the best a at the current step. Due to the nonlinearity of the problem, the computation and storage requirements of the M A P estimator are extensive. The following section discusses this difficulty, and suggests a solution. 5.3 The Computational Algorithm Since (5.124) is the formula to compute the M A P solution, the only question left concerns the computation of the gradient function VW(f). The following characteris-tics of V W ( / ) are of interest. • It can be found, by tracing the previous section, that although we have assumed Toeplitz structure for the matrix H and stationary statistics for the random processes / , Ni and n 2 , the derivation of the M A P filter does not need these assumptions. However for general space-variant image formation, with random process that are not stationary, the evaluation of (5.123) requires matrix multi-plications and consumes an impossible amount of computer time for any image with a typical number of samples, say a 256 x 256 image. • If the image formation process is linear, shift-invariant, as we have formulated in the previous section, then the general matrix H becomes block Toeplitz. Likewise, Chapter 5. Restoration by the Nonlinear MAP Filter 68 if the random process which generates the ensemble, from which the object image / is drawn, is a process which is stationary, and if the noise is also stationary, then the matrices Kf and Kn(f) can also be replaced by block Toeplitz matrices. All the block Toeplitz matrices can then be approximated by circulant matrices. • Assuming the matrices if, Kn(f), and Kf are approximated by circulants then all the matrix operation in (5.123) can be computed as convolutions, hence im-plemented by F F T . The computation of the M A P filter by the steepest ascent method with a Fourier domain approach is summarized in the following algorithm. Algorithm 1. Set k = 0, and Construct an initial guess f0. 2. Compute Kj1{fk - f) by the F F T . This is the first term in (5.112) 3. Compute the Fourier transform of g, fk, H, and K„(f) by the F F T , respectively, where Wn{ui) is the power spectrum of n, (n = Nif + n2). Invert (5.125) back into the space domain to obtain a vector yk, and form no*) (5.125) yk = K-1(fk)(g-Hfk) (5.126) 4. Compute v*(w,-) = /T(<*)n(w.-) (5.127) Chapter 5. Restoration by the Nonlinear MAP Filter 69 where * denotes complex conjugate. Compute the inverse F F T of (5.127), it yields vk = HtK;1(fk)(g-Hfk) (5.128) 5. Compute d[Kn(f)]/d[fk]i for / = 1,2, . . . , M , and transform the result into the Fourier domain d[fi * = 1 ,2 , . . . ,M where 7M{-} indicates the discrete Fourier transform. (5.129) 6. Multiply Yjt(u;,) by \d[Kn{uJi)}i to form / = 1 ,2 , . . . ,M (5.130) Apply the inverse F F T to (5.130) to obtain zi = ^^K-\fk){g-Hfk), d[f, /.= 1 ,2 , . . . ,M k i (5.131) 7. Compute the components of VW2 (the last two terms of (5.123)) by [VW a ] , = -2M - ylkzlk (5.132) = -2\HTK-\f){g - Hf)\i -(g- Hf)TK~1(f)^^j^-K~1(f)(g - Hf) 8. Compute — 0.5VWj, the second term in (5.123) using (5.133) 9. Form the derivative vector VW(fk) by VW(fk) = -Kj\fk - /) - -{VWX + VW2) (5.134) Chapter 5. Restoration by the Nonlinear MAP Filter 70 10. Use the golden section search to find a proper value of a f c , and generate a new estimate fk+i by fk+i = fk + cckVW{fk) (5.135) where fk+i > 0. 11. If the convergence criterion l / f c + i - / * ! < « (5.136) is satisfied, where e is a predefined positive constant, stop restoration process. Otherwise, repeat step 2 through step 10. By applying this algorithm, the reduction on computations is from 0 ( M 4 ) to 0(M2logM) in the one-dimensional case. 5.4 Exper iment Results The mathematical derivation of the M A P filter presented in the previous sections applies to one-dimensional signals, and with a few modifications, as will be shown in the next section, to the two-dimensional signals. However, the computational effort for the two-dimensional signals is very large and beyond our existing facilities. For instance, programmed in F O R T R A N and run on an VAX11/750, the restoration of a two-dimensional image of size 64 x 64 would spend more than 5 hours of C P U time. On the other hand, to process the same image distorted by motion blur, we only need approximately 14 minutes. Thus the examples attempted are for the one dimensional signals and for pictures distorted by motion blur. The M A P filter was applied to the restoration of images distorted by random motion blur with the presence of additive detection noise. The results were compared with Chapter 5. Restoration by the Nonlinear MAP Filter 71 those given by the modified Wiener estimator [25]. The first example is one dimensional. The original signal is plotted in Figure 5.15 a). The numerical values in this signal are from 0.0 to 1.0. This signal is blurred by a random impulse response generated through superposition of white noise Ni on a Gaussian shaped point spread function and additive noise n2. The blurred version is shown in Figure 5.15 b). The standard deviation of the Gaussian shaped PSF is 2.0. The random statistics used here are ot — 0.02 and o2 = 0.03, respectively. Restoration by the M A P and, the modified Wiener estimators are given in Figures 5.15.c) and 5.15 d), respectively. It is seen that the M A P estimator recovers the lost information, such as the first narrow spike, which is not recovered by the modified Wiener estimator. In the following examples, the M A P estimator is applied to restore random motion blurred pictures. Each pixel in the pictures is represented by the usual integer code value in the range (0, 255) which corresponds to an 8-bit uniform quantization of the intensity values. The random motion blur effect was generated by superimposing white Gaussian noise on the impulse response of motion blur, where the scalar K is used for normalization in order to model a lossless linear system, i.e. Y^iH(i) — 1. The standard deviation of iVi is o\. The additive noise n2 is also white Gaussian noise with standard deviation o2, and n2 is uncorrelated with N\. The first of these examples is an 120 x 120 image, which is part of the picture of an F16 jet fighter, as shown in Figure 5.16 a). The distortion parameters used in this example are K = 9, o\ = 0.006 and o2 = 4.0. The distorted image is given in Figure 5.16 b). Restorations by the M A P estimator and by the modified Wiener estimator [25] are shown in Figure 5.16 c) and 5.16 d), respectively. Since the blur effect and the noise effect are both severe in the distortion, the modified Wiener filter cannot clear (5.137) Chapter 5. Restoration by the Nonlinear MAP Filter 72 the noise completely. Restoration by the ordinary Wiener, (ignoring the presence of the noise in the PSF) gave RMS error = 44.7. Thus the ordinary Wiener filter just makes things worse. The last example is a 116 x 116 image "Face", as in Figure 5.17.a). The distortion parameters used in this example are K = 13, C\ — 0.008 and o2 — 5.1. The distorted image is given in Figure 5.17 b). Restorations by the M A P estimator and by the modified Wiener estimator are shown in Figure 5.17 c) and 5.17 d), respectively. The ordinary Wiener estimator gave RMS = 62.0. To further differentiate the results given by the M A P , and the modified Wiener filters, Figure 5.18 show the subimages of Figure 5.17, the right eye of Face, and Figure 5.19 a) to c) show the difference images of Figure 5.18 b), c) and d), with the original Figure 5.18 a), respectively. From the difference images, it might seem that the Wiener filter more correctly reverses the blur effect than the M A P filter. However, this observation is not true. If we add white noise to the image restored by the M A P filter, Figure 5.18 c), so that the resulting RMS error is equal to that of the Wiener restored image Figure 5.18 d), then the difference image is shown in Figure 5.19 d). This resulting difference picture looks almost the same and h&s the same variance as Figure 5.19 c), which is the difference picture of the Wiener restored image. Thus we can conclude that the deblur functions of the two filters are comparable, but the M A P filter clears much more noise than the modified filter does. The RMS errors for the three examples are tabulated in Table 5.5. Since the phenomenon of motion blur is that of one dimensional distortion, the restoration was actually carried out row by row. From the examples given, the superi-ority of the M A P estimator over the Wiener estimator is obvious. The disadvantage is that the M A P estimator requires larger amounts of computation and working memory space. Chapter 5. Restoration by the Nonlinear MAP Filter 73 Table 5.5: RMS errors of the M A P technique RMS Test 1 Test 2 Test 3 Blurred image 0.183 27.72 20.75 Modified Wiener 0.160 24.65 24.36 M A P technique 0.141 19.12 14.41 Although the convergence of the algorithm has not been proved, our experience is that convergence always occurred. For our examples, we used signals of length 128 and the width of the blur function K varied from 3 to 15. The choice of e in (5.136) is important. In [40] the authors discussed one way to find e, but their scheme is not suitable for gradient descent methods. Cooper and Steinberg [41] suggested determining e experimentally. We followed their suggestion and found that values for e is in the range 0.1M < e < 0.2M where M is the length of the signal. Further iterations do not provide better visual results. The algorithm can converge to the global maxima is the solution set is convex, otherwise it may stop at a local maxima. It is noticed that the relevance of the root mean squared error criterion used by the M A P estimator is questionable, given the enigmatic nature of the human visual system. The motivations for using the RMS are, of course, historic precedent and mathematical tractability. Better criteria need to be developed to facilitate the measurement of restoration quality. Chapter 5. Restoration by the Nonlinear MAP Filter 74 Figure 5.15: An one dimensional example, a) The original signal; b) the blurred signal; c) Restoration by the M A P filter; d) Restoration by the Wiener-based filter. Chapter 5. Restoration by the Nonlinear MAP Filter 75 Figure 5.16: Image F16: distortion and restorations, a) The original image; b) the blurred image; c) Restoration by the M A P filter; d) Restoration by the Wiener-based filter. Chapter 5. Restoration by the Nonlinear MAP Filter 76 Figure 5.17: Image Face: distortion and restorations, a) The original image; b) the blurred image; c) Restoration by the M A P filter; d) Restoration by the Wiener-based filter. Chapter 5. Restoration by the Nonlinear MAP Filter 77 Figure 5.18: Subimages of Face: the right eye. a) The original image; b) the blurred image; c) Restoration by the M A P filter; d) Restoration by the Wiener-based filter. Chapter 5. Restoration by the Nonlinear MAP Filter 78 c) d) Figure 5.19: The difference images of the right eye. a) The difference between Figure 5.18 a) and 5.18 b); b) The difference between Figure 5.18 a) and 5.18 c); c) The difference between Figure 5.18 a) and 5.18 d); d) The difference between Figure 5.18 a) and the noisy version of 5.18 c); Chapter 5. Restoration by the Nonlinear MAP Filter 79 5.5 Generalization to Two-dimensions Although the MAP filter was formulated for the one-dimensional case, it can be generalized to the two-dimensional case. The mathematical derivation of the MAP filter in the two-dimensional case is principally the same, but there are few terms which should be changed to forms suitable for the two-dimensional representation. The first change is the expression for the elements in the derivative of the covariance matrix Kn(f) with respect to / shown in equation (5.121). This should be modified to [**n(/)ltf fr-i+l,l-} + l if r + i- 1 < A f , r - i + 1 > 0,s + j - 1 < M,s -j + 1 > 0 + / r + » - l , « + j ' - l / r - i + M - i + l ^ ' fr+i-l,»+j-l if r - i + 1 > 0,s - j + 1 > 0,r + i - 1 > M o r r - t + l > 0,3 - i + l > 0,3 + j ' - l > A f r - f t - 1 < M,s + j - 1< M,r - » + 1 <0 or r + i -l<M,s + j - 1 < M,s- j + 1<0 0 otherwise (5.138) The second change is that equation (5.120) should be modified to its corresponding two dimensional form (5.139) Chapter 5. Restoration by the Nonlinear MAP Filter 80 The last change to be made is that the Fourier computation formulae obtained in the algorithm should be changed to two-dimensions. Thus the Fourier argument u> is replaced by the pair (u>x,wv). 5.6 Generalization to Nonlinear Recording System As mentioned earlier the M A P filter was first developed for the conventional image restoration system to circumvent the sensor nonlinearities problems. Our previous discussion on the random distortion problem can also be generalized to incorporate the sensor nonlinearities. The distortion model for such a recording system has been given in (5.102), we present the discrete version in the following g — S{Hf + Nif} + n2 (5.140) where S{.} indicates the nonlinearities. When <T1( the standard deviation of Nx, is small, Taylor's expansion can be used to separate the stochastic term N\f from the nonlinearities. In the one dimensional case, taking the first two terms of the Taylor's expansion yields g = S{Hf) + SiNif + n 2 (5.141) where Si, is a diagonal matrix Sb = with as du u=6i 0 dS du Chapter 5. Restoration by the Nonlinear MAP Filter 81 Under the statistical assumptions of the linear recording system discussed in Section 4.2, the vector W(f) = lnp{f\g) is W(f) = lnp(f\g) (5.142) = -0 .5{ In | / f B ( / ) | + ( / - 7 ) r V ( / - 7 ) + (g - S{Hf})TK?{f){g - S{Hf})} - \nq(g) where Kn(f) = EliS^r + n2){ShNJ + n 2) T] Following a procedure similar to that in Section 4.2 and 4.3, it can be shown that the computation of V W ( / ) involves only F F T and convolutions. Thus we can use an optimization algorithm, such as the steepest ascent, to solve the restoration problem. The computational complexity for this case can be shown to be of 0(M2logM) per iteration. 5.7 S u m m a r y The application of the nonlinear maximum a posteriori filter to the restoration of images distorted by random point spread functions, and additive Gaussian detection noise is discussed in this chapter. Due to the nonlinear nature of the problem, a nonlinear restoration approach such as this M A P filter is expected to more closely simulate the reverse process of the distortion. The solution is approached iteratively. Unlike the linear based filters, a closed form for the M A P filter can not be derived. We propose an optimization procedure to compute the restoration algorithm. One advantage of this optimization procedure is that the constraint of positiveness of the radiant energy can be incorporated. Other constraints may be added if desired. Chapter 5. Restoration by the Nonlinear MAP Filter 82 We test the M A P filter and compare it to the modified Wiener filter for randomly distorted signal and for several randomly motion blurred images. Values of the RMS errors are listed. In all cases, the nonlinear M A P filter gives better results than the modified Wiener filter. The MAP.technique formulated in this chapter applies to general random distortion problem, but the computational algorithm of the filter is restricted to the linear, shift-invariant point spread function, and stationary image and noise covariance statistics. For the general case, it is very difficult to make the M A P estimate a practical and economical alternative to the linear filters. Chapter 6 Deblurring Random Time-varying Blur 6.1 Introduction In this chapter we wish to address an image-restoration problem in which the point spread function of the imaging system not only is random, but also varies with time (because of, for example, random fluctuation of the system's pupil function) while the object image remains fixed. It is not impractical that there exists random fluctuation in the pupil function of an optical system. This phenomenon occurs in cases such as imaging through turbulent atmosphere, imaging under water or in a dusty atmosphere, or with a randomly vibrat-ing scanning camera. The effect of atmospheric turbulence on imaging has long been paid attention to due to its important applications in astronomy and in remote-sensing [42], [43]. The attempts to restore images distorted by the atmospheric turbulence have been reported in the literature [44-52]. In [26] the problem of restoring signals distorted by random time-varying impulse response has been studied via the minimum variance unbiased methodology. This method is computationally burdensome, consumes a lot of time, and may suffer from instability due to ill conditioning. In [27] the Backus-Gilbert technique has been applied to the problem of restoring images distorted by random and time-varying blur. The results of the latter study are good, but they unfortunately de-pend on the adjustment of a parameter and are thus subjective. Our purpose here is to explore some other restoration techniques for this random time-varying blur problem. 83 Chapter 6. Deblurring Random Time-varying Blur 84 6.2 Problem Formulation The random time-varying distortion problem is modeled by g{x,y;t) - / / h(x - xuy - y1;t)f(x,y)dxldyi + n2[x,y;t) —oo = h ( x , y ; « ) * /(x,y) +n2(x,y;<) (6.143) where h(x,y;i) is the random time-varying point spread function. We have also al-lowed for the measurement noise n2(x, y; t) being time dependent. The recorded image g(x,y; t) is now time-varying and random. It is recorded in the time interval [0, T]; the object / is to be estimated by assuming that some statistical properties of the stochastic processes h(x,y;t) and n2(x,y; t) are known. In particular, we assume that h(x ,y;£) is stationary in the weak sense in time and space and has a known mean Furthermore, h(x,y;<) is assumed to be ergodic. Similarly, n2(x,y;t) is assumed to have zero mean and an autocorrelation function h[x,y) = E[h(x,y;t)} (6.144) and a known spatioternporal autocorrelation function Rh{x, y; T) = E[h{0,0; t)h(x, y; t + r)] (6.145) Rn2{x,y,T) = E[n2(0,0;t)n2(x,y;t + r)] (6.146) We may write equation (6.143) in the alternative form g(x, y; t) = h(x, y) * f{x, y) + n(x, y; t) (6.147) where n(x,y;t) = ni(x,y;f) * /(x,y) + n2{x,y;t) (6.148) Chapter 6. Deblurring Random Time-varying Blur 85 The function n 1(x,y;t) is derived from m{x, y; t) = h{x,y; t) - h{x, y) (6.149) which may be regarded as the time-varying noise in h(x, y; t). The observed available data to be processed for such a problem is extremely large, therefore it is necessary to compress the data without losing useful information. Three possible approaches are given below: 1. A time-averaged image may be constructed on which the estimations are based. 2. A time-averaged spatial correlation is computed from which the estimations are calculated; all other data are ignored. 3. Samples of the image at a finite number of time instants (frames) are measured; the estimations are based on the obtained sequence of images. In the following sections we shall discuss estimations that are based on the com-pressed data in approaches 1 and 2 briefly; our concentration will be on the third approach. By using approach 3, we restore the unknown object image from the observed multiframe images. This problem, as well as the problems of restoration of multi-spectral satellite images, color images, and multiple sensor images, all fall under the subject of the restoration of multi-channel images. Other techniques on restoration of multi-channel images can be found in [53-56]. Because of the multiframe nature of the data, the size of the problem is expected to be extremely large for all practical cases. In this chapter, we approach this difficulty by using two different schemes based on the Wiener criterion. The first employs a Chapter 6. Deblurring Random Time-varying Blur 86 Karhunen-Loeve-like block-diagonalizing transformation that makes the transformed frames uncorrelated. This technique, which is valid when the spatiotemporal correlation of the noise is separable in space and time, reduces the size of computation to that of a single frame. This technique is similar to an approach recently developed for the one-dimensional problem via the minimum variance unbiased scheme [26], and also to a method developed for the restoration of multispectral images [53] in that all the three methods use the Karhunen-Lo^ve transformation. The second scheme requires that the structure of the correlation matrix of the noise be block Toeplitz. This case is valid under the usually assumed condition that the noise is stationary in time and in space, and is a less strict condition than the above time-space-separable. The computation algorithm here is also different from the above case, and requires more time. Here it will be shown that the solution may be reduced to that of solving a series of Toeplitz system of linear equations in the Fourier domain. By this the size of the problem is reduced considerably. However the computational effort of this scheme is still greater than that in the time-space-separable case, using the Karhunen-Loeve transformation. We shall use numerical examples to show that these schemes' are capable of produc-ing high quality results. The computations will be done in the frequency domain using F F T and circulant matrix approximation since this requires far fewer computations than those required for space domain manipulation. Chapter 6. Deblurring Random Time-varying Blur 87 6.3 Restoration Based on the Average Image The mechanism of restoration of the random time-varying blur problem based on the average image 1 fT 9(x>y) = TJ0 9{x,y;t)dt (6.150) is straightforward. Although such a estimation does not utilize all the information contained in the image g(x, y;t), it has the advantage of simplicity. By substituting (6.143) into (6.150), g{x,y) can be related to /(x, y) by g(x, y) = h{x, y) * f(x, y) + n(x, y) (6.151) n(x,y) = ni(x,y) * /(x,y) + n 2(x,y) (6.152) where n 2(x,y) is the time-averaged measurement noise 1 fT n2(x>y) = ji JQ n2{x,y;t)dt (6.153) and 1 fT rti(x, y) = h(x, y) - h(x, y) = — / nt(x,y; t)dt (6.154) 1 Jo in which h(x,y) is the time averaged impulse response = ^ I*h{x,y;t)dt (6.155) l Jo The noise ni(x,y) represents the effect of the finiteness of the averaging time T. For if T is sufficiently long, then, because the process h(x,y;() is ergodic in the mean, h(x,y) approaches h(x, y), and n 1(x,y) approaches 0. Equation (6.151) is then in the conventional image restoration form. Yet now the noise n is linearly dependent on the unknown object / . This deconvolution problem is similar to deconvolution in signal-dependent noise, except that the noise dependent on / rather than on / x h. Chapter 6. Deblurring Random Time-varying Blur 88 It is recognized that the restoration based on average image can be addressed by the extended linear techniques: the Wiener filter, the M V U filter, or the constrained least squares filter discussed in Chapter 3 and in Chapter 4. 6.4 Restoration Based on the Time-averaged Spatial Autocorrelation of the Image Another function on which the estimation may be based is the time-averaged spatial autocorrelation function 1 fT r9{x,y) = - Jo Rg{x,y;t)dt (6.156) where Rg{x,y,t) = I I g{xi,yi;t)g(x1 + x,y1 + y;t)dxidy1 (6.157) J J-OO is the spatial autocorrelation of the image g(x,y;t) at time t. Alternatively, the time average of the squared absolute value of the image may be recorded as 1 fT Sg{ux,uv) = - ^ \G{ug,uv;t)\2dt - (6.158) where G{ux,wy,t) = J J^og{x,y;t)ei^M'"'Uxdy (6.159) is the Fourier transform of g(x,y;t). The function rg(x, y) may be computed from the function Sg(u>x, w„) by an operation of inverse Fourier transform. In the limit T —• oo the function Sg(u>x, u>y) represents the spatial power spectral density of the image. In the absence of the measurement noise, i.e. n 2(x,y;t) = 0 for all t 6 [0, T], equation (6.158) implies that = \F{ujx,uv)\2Sh{ujx,uy), (6.160) Chapter 6. Deblurring Random Time-varying Blur 89 where Sh(u>x,wv) = ^ / r |H( W l ,uy,*) | 2 <fc, (6.161) 1 Jo with H(ux,<jjy;t) being the Fourier transform of h(x,y;t). Equation (6.160) implies that rg[x, y) = r h (x, y) * rf(x, y) (6.162) where rj,(x,y) <-> is the time-averaged autocorrelation function of the stochastic process h(x,y;t) and r/(x,y) is the autocorrelation function of the object image / : /oo f{xu yi)f{xi + x,yi + y)dxldyl (6.163) -oo Note that, in the limit T —• oo, r\(x, y) approaches the known autocorrelation function Rh(x,y;0). The difference between r^(x,y) and i2j,(x,y,0) may be regarded as noise. For a finite time T and in the presence of measurement noise, equation (6.156) is modified rg(x, y) = Rh(x, y; 0) * r,(x, y) + n 0(x, y) (6.164) where no(x,y) is noise due to the effects of measurement noise and the finiteness of the observation time T. The function r/(x,y) may be computed from rg(x,y) by deconvolution. Note that the noise n 0 (x ,y) is again a function of the unknown image / , this time having a quadratic instead of a linear dependence on / . Recovering the image / from its auto-correlation function r/(x,y) is not a easy task [48], yet the reward of this method, as compared with the previous one (based on measurement of the averaged image), may be great. As an explanation, we write the autocorrelation function i?h(x,y;0) as Rh{x, y; 0) = h{x, y)h{0) + c h(x, y) (6.165) where Ch(x,y) is the autocovariance function. Typically, cj,(x,y) is a much narrower function than the function h(x,y) itself. Therefore, we can think of equation (6.165) Chapter 6. Deblurring Random Time-varying Blur 90 as a system whose transfer function h(0)H(ujx,uv) + C\\(UJX,UV), (where H(wx,wv) and Ch(wx,u>„) are the Fourier transform of h(x,y) and Ch(x,y), respectively), contains a wide-band component that allows high spatial frequencies of the image to be recovered. This is the basis of the well-known technique of speckle interferometry [47-51]. Although many algorithms have been developed for solving the signal-recovery prob-lem in speckle interferometry, the problem of the dependence of the noise due to finite measurement time on the unknown function / has, to our knowledge, not been ad-dressed. 6.5 Restoration Based on Image Samples In this section we shall first discretize the random time-varying problem in the first subsection, and throughout the rest of this chapter, we shall use the matrix vector description of the linear, shift-invariant system for image formation. 6.5.1 Discretization of the Problem For time-varying imagery, the discretization of such a model begins with the as-sumption that S frames of the image g(x, y;t); s = 1,2,... S are assumed to be observed during the time interval [0, T], The original object f(x,y) is to be estimated from the observed frames g(x,y;t,)'s. The sampling of this model for digital computation [15] leads to S sets of linear equations M M g{j,k,ta) = £ h ( j - m , f c - / ; i s ) / ( m , / ) 5 = 1,2,... S (6.166) m=l( = l and these sets of equations are replaced by their matrix vector description g{tt) = H ( t , ) / + n2{ta) 5 = 1,2,... ,S (6.167) Chapter 6. Deblurring Random Time-varying Blur 91 where the vectors / , g(tt), and n2(ts) are the lexicographic orders of the two-dimensional sample matrices of size M x M into the one-dimensional vectors of length M2. All of the matrices H ( £ , ) , are block Toeplitz matrices of size M2 x M2, and thus may be approximated for computational ease, by circulant matrices. For our time-varying imaging system, we assume that all the observed image frames are perfectly registered with one another. By "perfectly registered" we mean that all coincident pixels in the sampled image arrays correspond to the same point on the object being imaged. Thus if we pick the coordinate (j, k) and project the pixels at this coordinate in different frames back through the optical system to the object, all these pixels will project back through the optics to the same point on the object. Our approach is based on stacking the vectors that represent the sampled frames so that one long vector containing all the measured data is formed. We adopt the following notation 9 = 9(h) 9{tt) 9(ts) n2 = n2{ti) n2(t2) n2(ts) (6.168) to symbolize the stacking (or the lexicographic order) of all the S vectors of length M 2 into the single vector of length SM2. The matrices of degradation for the multiframe system is also stacked to form matrix H such that: Chapter 6. Deblurring Random Time-varying Blur 92 H(*i) H(t a ) H - (6.169) H ( « 5 ) Obviously, H is of dimensions SM2 x M2. Thus the multiframe time-varying imaging equation model takes the matrix vector form 9 — H / + n 2 (6.170) where the components are as represented in equations (6.168) and (6.169). Without loss of generality, we write the time-varying random blur functions H(£„) as H(t.) = H + Niit.) s = l,2,... ,S (6.171) Here H is the mean part of the random blur which is time-invariant and Ni{tt) is the zero mean random fluctuation around the mean which varies with time. Using (6.171) equation (6.170) can be rewritten as g = Hf + n n = Nif + n 2 (6.172) where H H N!(t2) H (6.173) H Nx(ts) Chapter 6. Deblurring Random Time-varying Blur 93 The submatrices H and Ni(t,) are black Toeplitz with similar structures as H and Ni shown in Chapter 2. The discretized elements in the above equations are from their continuous counterparts by in which where fjxjy = / O x A z , J v Ay) , [»('•)].•-,<, = ^ t ' . A i , ,ivAy;tt), (»(*.)].•...•, = n(jxAx ,iyAy;t,), = n2(ixAx,iyAy;t„), = h[(is- Jx)Ax, (ty -- ni[[ix-•jx)Ax, (iy --jy)Ay;t, 1)+;;, i = M(iy - 1) + »', ** = 1,2,. . . M , iv = 1,2,... . M , ix = i ,2 , . .. M, iv = 1,2,.. . M. (6.174) It is equation (6.172) that we shall use in developing the restoration algorithm. 6.5.2 The Wiener Filter Criterion Our restoration procedures are presented in this section. Both approaches are based on the Wiener filter criterion. The first procedure is based on the assumption that the spatiotemporal correlation of each of the noise sources N\ and n2 is separable in time and space. The second procedure requires that the structure of the correlation matrix of the noise n be block Toeplitz. In both cases, we assume that if , the mean of the point Chapter 6. Deblurring Random Time-varying Blur 94 spread function, in (6.167) is known and that the covariances of the noise elements Ni and n2 are known. The restoration of the object / in the model (6.172) by means of the Wiener esti-mation has been given in Chapter 3 / = RfHT[Rn(f) + HR/HT]~1g (6.175) A where / is the estimat of the object image / , Rf is the correlation matrix of / , i.e., Rf = E[ffT) (6.176) In the random time-varying distortion case, the correlation matrix Rn of the object dependent noise n is Rn(f) = RklA){f) *i u )(/) (6.177) Rks's)(f) Here, the submatrices #(/'*)(/) are determined from the cross-correlation /#•*>(/) = E[n{ti)n{h)T) (6.178) By using a similar procedure as in Section 2.4.2, it can be shown that the elements in Rn'kHf) i S g i v e n by l(R{J'k)(f)km}l,n = E E E E ^ ^ l f e ) ) M k p I ( ^ ( ^ ) ) m , ] n , r } £ [ / p , , / r , ] r a p q +E{[n2{tj)}l,i[n2{tk)}n,m} (6.179) It can also be shown that the matrices Rn{f) and each of {R^'k^(f)} are block Toeplitz as long as the noise rii(x,y;t) and also n2(x,y;t) are each stationary in time and in space. Chapter 6. Deblurring Random Time-varying Blur 95 Although the solution of the multiframe estimation problem based on the Wiener filter (3.50) is straightforward, the dimensionality of this problem for even moderate S and M will be large enough to pose a computational problem. Indeed for practical image sizes and when the number of frames is not small the computational and the active memory requirements are prohibitively large. However, we can alleviate this problem in certain special but practical cases. 6.5.3 The Time-Space-Separable Case Under time-space-separable conditions, we show that the formulation and compu-tational effort of the multiframe problem can be reduced to that of the order of a single frame. By the term " time-space-separable ", we mean that the correlation functions of each separate as a function of time and another function of space. Conditions under which impulse response functions are separable are not uncommon; the assumption of separability is widely used. This occurs in imaging systems for which the fluctuations of the pupil function are spatially homogeneous, stationary, and cross-spectrally pure, i.e., the temporal spectra at all points and the cross spectra at all pair of points have the same shape. Such a model has been used to model atmospheric turbulence [52]. Spatiotemporal separability implies that Ni(t,) and n 2 ( £ 8 ) have correlation matrices that are also separable: £ { [ ( t f i ( * ; ) ) < J i * [ ( ^ i M ) m , . ] n ^ (6-180) and £{M*/)]l,.-M*fc)]n,m} = {R[%A(Ri S)h*nkn (6.181) The superscripts (t) and (s) indicate time (frames) and space, respectively. In equation Chapter 6. Deblurring Random Time-varying Blur 96 (6.180) the condition of stationarity has been used. Since n 2(x, y; i) represents detection noise which is normally known to be uncorrelated in time and in space, each of R^^ and JRJ^ may be assumed to be proportional to an identity matrix. In order to make use of the property of time-space-separability, the Wiener filter estimation in Eq. (3.50) must be expressed in a different form. It has been reported in the literature [15] [59] that if the inverse of Rn, R/, and H exist, then an equivalent form of the Wiener filter is / = [HTR?{f)H + R]l}-lRTK\f)9 (6-182) In Appendix D we show that (6.182) remains true if the condition on the matrix H is relaxed to that of a full column rank, and thus still holds for our problem. Substituting (6.180) and (6.181) into (6.177) and (6.179) yields Rn(f) = ® R['\f) + J#> ® R['\ (6.183) where [ ( i? i ^ ( / ) ) , . ^ (6-184) p q r a and ® denotes the Kronecker-matrix product. f To formulate and reduce the computational effort of the problem to that of a single frame we consider a Karhunen-Loeve-type transformation [57]. Consider a transforma-tion on the image g defined by the following formula: g' = {A ® IM)g (6.185) where A is an S x S matrix and IM is an M x M identity matrix. We similarly transform the noise n : ri = (A ® IM)n (6.186) Chapter 6. Deblurring Random Time-varying Blur 97 Therefore the blurred image is represented in an equivalent transformed form: g' = H'f + n' (6.187) where H1 = [A ® IM)H (6.188) In the transformed domain, it is easy to show that the correlation matrix of the noise n' is: R'n(f) = E[n'n'T] (6.189) = E{{(A®IM)n)[(A®IM)n}T} = (A®IM)E[nnT}(AT®IM) = {A®IM)Rn{f){AT®IM) Substituting (6.183) into (6.189) yields R'n(f) = (A®IM)[R[t)®R{a)(f)+Rit)®Ri,))(AT®IM) (6.190) = (A® IM)[RM ® R[a)(f)](AT ® ITM) + (A® IM)[Rit] ® RP}(At ® ITM) = [AR?AT]®R[l){f) + {ARit)AT)®R2s) In order to use the transform we must assume either ni(x,y;t) or n2(x,y;t) to be temporally uncorrelated. Thus either R[^ — I or R^ = / , where I is the 5 x 5 identity matrix. Since the measurement noise is normally assumed uncorrelated between frames, we select R? = I Hence (6.190) is simplified to R'n(f) = [AR^AT\ ® R['\f) + (AAT) ® R(2s) (6.191) Chapter 6. Deblurring Random Time-varying Blur 98 The matrix R^, by virtue of being a correlation matrix, is positive definite and sym-metric. Thus, the transformation matrix A can be selected such that ARX'A (6.192) where A is a diagonal matrix. A = in L22 A S S j It is easy to identify that the diagonal entries in A are the eigenvalues of R*i \ and the row vectors in the matrix A are the corresponding eigenvectors [58]. We normalize the eigen vectors so that AAT = I Then equation (6.191) is reduced to K(/) = A ® i 2 S 8 ) ( / ) + ^ ^ 8 ) (6.193) (6.194) The above expression represents a block-diagonal matrix. The ith nonzero diagonal submatrix is [R[t{f)]u = A4tR['){f) + Ri') (6.195) The estimation problem (6.182) in the Karhunen-Loeve transform domain becomes , f = [R^R'-HnS' + R]l\-lR'TR'-l{f)g' (6.196) By substituting (6.188) and (6.195) into (6.196), and after some manipulations we obtain f = [H TJ2 c}TiH + RjrW E e , - iy («,-)] (6.197) Chapter 6. Deblurring Random Time-varying Blur 99 where T, = [AURPW + RI'Y1 (6.198) s c« = T,AU (6.199) i = i and g'(tj) is the j th transformed image frame. Equation (6.197) is now recognized as similar to that of the single frame restoration problem. The computations now require inverting matrices of size M 2 x M 2 instead of SM2 x SM2. Even though an eigenvalue problem has to be solved, substantial reductions in needed memory as well as computation are achieved. More important is that under the condition of spatial stationarity for the spatial correlations in the time-space-separable case, i?[*' is block Toeplitz, and thus can be approximated by circulant matrix Ru . Similarly if, T,-, and Rf can be approximated by circulant matrices i f c , r , c , and Rfc. Hence all calculations in equation (6.197) can be performed by Fourier computation. Indeed, the two dimensional Fourier domain form of the Wiener restoration in (6.197) has the form *,/ x _ ^ * ( ^ x , ^ ) E j =iC J T i ( u ; I , u ; t > ) G , ( u ; g , u ; t / ; t < , ) - | i j ( w . > W r ) | » I ^ l c ? r < K | W f ) + ^ > . , W f ) 1 6 - 2 0 0 ) where Wf(ux,uv) is the power spectrum of the image / . Although T,- is not block Toeplitz, its inverse r,- 1 = AuRl']{f) + Rii] is block Toeplitz. Therefore the approximation is relevant, and the ease of the feasibility of the solution is now achieved. By doing so the computational effort is reduced from 0 ( 5 3 M 6 ) to 0{SM2logM). If we take iZj*' = 7, then the resultant estimate is similar to the one shown in (6.197), except for that Ti = [Rif)(f) + AiiRt')\-1 (6.201) Chapter 6. Deblurring Random Time-varying Blur 100 where Aj, is the tth eigenvalue of R2 , and the matrix A is composed of the eigenvectors of Ri l). 6.5.4 Direct Wiener Filter Implementation In this subsection, an approach based on the direct form of the Wiener filter (3.50) is implemented to restore time-varying blur. Conditions for the application of this approach, as mentioned before, is broader than those posed for the Karhunen-Loeve transformation used in the time-space-separable case. The only assumptions required here are that the submatrices R%'^ in the noise correlation matrix Rn{f) are all block Toeplitz and Rn{f) is also block Toeplitz. Such conditions are true under the usual assumptions of stationarity in time and in space. With reference to equation(3.50), if we denote the Wiener filter by L, the direct form of L is L = RfH T[Rn(f) + HRfR T}- 1 (6.202) where (6.203) In (6.203) Li is the tth submatrix of L. Taking the transpose of (6.202) yields L T = [Rn{f) + HRfH T]~ lHRf (6.204) Premultiplying both sides of (6.204) by [ # „ ( / ) + HRfH T] results in [/^(/) + HRfH T]L T = HRf (6.205) Chapter 6. Deblurring Random Time-varying Blur 101 Equation (6.177) shows the structure of the matrix Rn{f). Inserting the definition of H in (6.173) into HR/H gives RRfH T = HRfH^" HRfH^ HRfH T HRfH T HRrR T (6.206) which contains S x S identical submatrices HR/H T. For the sake of simplicity, we introduce the following notations: A = HR/H T E = HRf Then (6.205) is explicitly written as: RgV{f) + A R^(f) + A Rg>lHf) + A Rg#(f) + A E = E E (6.207) Expression (6.207) represents a system of linear equations in SM 2 x M 2 unknowns. The entries in Lk can be calculated row by row from (6.207). But the size of the system is very huge, SM 2 x SM 2 with SM 2 unknowns, and must be repeatedly solved M 2 times. Space domain calculation is extremely expensive. The alternative solution is, again, to make use of the Toeplitz property of the system. The system matrix in (6.207) is a three level block Toeplitz matrix and its normal solution requires three level (or three dimensional) F F T . However, the approximation of the system matrix of (6.207) by circulant matrix is not a good approximation here. But the circulant matrix approximation does apply to its submatrices. Chapter 6. Deblurring Random Time-varying Blur 102 We approximate the matrices H, Rf, R^,k^(f), j,k = 1,2,... S, by the circulant matrices Hc, RfC, and \RH'k^{f)]c- The approximated version of (6.207) is where and <f>21 <f>22 L\ = Ee (6.208) • • • <f>SS . LTS Ec <t>fk = Hc(Rf)eHj + {R^(f))c (6.209) Ee = Hc{Rf)c (6.210) <pki and Ec are both block circulant. Since this system has multiple unknowns, if we find a simple solution to any of them, the rest can be solved following the same procedure. Let Lu be defined as the tth column of Lf and let e,- be the t'th column of the right hand side Ec. Then the solution of Lu involves solving the system of equations <t>n <f>u • • • <f>2l <t>22 L2i = e, (6.211) • • • • 4>ss . Lsi . e, Since all the submatrices in the syst em matrix of equation (6.211) are now block circulant, we use Fourier transform to find the unknowns. For obtaining the discrete Fourier transforms we rewrite (6.211) as: 1=1 5 = 1,2,... 5 (6.212) Chapter 6. Deblurring Random Time-varying Blur 103 Since every term in (6.212) is a function of two variables in space and since <i>uLu represents convolution, the discrete Fourier transform of (6:212) is: s £ * w ( w « , w y ) . C K ( w 1 B > w , ) = £,-(w„w„); k = l,...S; x,y = l , . . . M (6.213) i=i Thus for frequency pair (wx,uv), we can form a system of linear equations. * u ( w „ w y ) £ 1 , - ( w s , i i ; I , ) + . . . + *is(w*>w»)"Cst(w*.Wy) = <fi(w«,w„) (6.214) Thus for each frequency pair (wx, wtf) we must solve a system of linear equatoins of size 5 x 5 , i.e., in 5 unknowns. To solve the system (6.211) we need to solve M 2 such linear systems. The computational burden, although it appreas large, is still greatly reduced form that of solving (2.211) in SM2 unknowns in the space domain. Since solving an N x TV system of linnear equations requires computations of 0(NS), we have sofar reduced the computations from 0 ( 5 3 M 6 ) to that of 0 ( 5 3 M 2 + SM2logM). Further reduction can be achieved by the observation that the form of the matrix in (6.214) is Toeplitz. Thus the computational effort is further reduced to 0 ( 5 2 M 2 + SM2logM). The inverse Fourier transform on £H(UJU, w„), / = 1,... 5; u, v = 1,... M provides the solution L«, / = 1,... 5. The other columns of Lf in (6.208) can be similarly calculated. Since the system matrix in (6.208) is invariant with respect to each column of L T , the matrix in (6.214) is also invariant. Therefore, matrix factorization in (6.214) is only needed once when computing the first column of LT. This factorization result is kept to compute all the subsequent columns [ £ l t £ 2« • • • £si], 1• = 1> 2, • • • M. Chapter 6. Deblurring Random Time-varying Blur 104 6.6 Experiment Results We have developed the theory of restoration of images distorted by a random time-varying point spread function in the previous section based on a complete multiframe model. In the following we will look at some examples. The pixels in the digitized images are represented by integer values in the range (0, 255). The average impulse response function h(x,y) used in all the examples is of Gaussian shape since this is known to model the effect of atmospheric turbulence. (x2 + v2) h(x,y) = Cexp{K J y '} (6.215) In (6.215) C is a scalar to ensure that E/i(x , y ) = l The noise we use in the examples is time-space-separable. In that case both of the Wiener filter implementations discussed in the previous section can be applied to the same problem. We assume that the impulse response noise is uncorrelated in space, i.e., R[f) = I (6.216) but correlated in time, in accordance with a first order Markovian process, so that = o-lX^ A < 1 (6.217) where A is the correlation coefficient between any two successive frames. Using (6.216) the spatial part of the correlation matrix of the noise given by equation (6.179), there-fore, becomes [(R{i\f)km)i,n = (M - |m - i\)(M - \ n - l\)[(Rf)i,m}i,n (6-218) Chapter 6. Deblurring Random Time-varying Blur 105 Under the conditions discussed in Section 6.5.2 this matrix is block Toeplitz and sym-metric. The additive noise represents detection noise, and thus is uncorrelated in time and space. Thus where <72 is the standard deviation of n 2(x,y;i). Now we have completely defined the noise properties for the examples. The distorted images g(tk) are generated from / and H by using equation (6.172) and adding the simulated noise that satisfied the properties discussed above. The Markovian noise ni(x,y,tK) is generated by means of the first order autoregressive model: where e(x, y; k) is an M x M matrix of uncorrelated random numbers of variance o\ = (1 - A 2 )a 2 . Equation (6.220) is applied successively starting from a sufficient large value A:0 = 45 so that the transient effects die away. The additive noise n2(J*) is accomplished by generating uncorrelated random numbers. 3 In our experiment, three examples are tested. Four frames (S — 4) are generated for each of the examples. The first example is an 56 x 56 image, "Eye", shown in Figure 6.20. The impulse response function extends over 9 x 9 pixels. The resulting blurred images are shown in Figure 6.21 a) - d) , respectively. In this example the noise statistics used are: 0\ = 0.0006, o2 = 5.1, and A = 0.8, respectively. Restoration based on the time7space-separable model by applying the Karhunen-Loeve transformation is given in Figure 6.22 a). The estimation of the direct Wiener filter application discussed in Section 6.5.4 is shown in Figure 6.22 b). For comparison, every frame is restored individually by the modified Wiener filter [25] (which is the /4S) = / (6.219) «i(«•»!/;'*) = Afi!(a:,y;tfc_i) + e(x,y; fc), * = -ko,... 1,2,... S (6.220) Chapter 6. Deblurring Random Time-varying Blur 106 Wiener filter applied to the restoration of a single frame distorted by random point spread function and ignoring the knowledge about its time variation). The restored image of one of these frames, the upper left image of Figure 6.21, is shown in Figure 6.23 a). We also show the restoration based on time-averaged image in Figure 6.23 b). The RMS errors are listed in Table 6.6. From the comparison, it is obvious that Table 6.6: RMS errors of the 1st example in time-varying case. RMS Frame 1 Frame 2 Frame 3 Frame 4 Blurred image 16.35 16.68 17.11 16.46 Individual 18.60 18.90 18.32 18.73 Time-averaged image 19.22 K - L transform , 13.91 Direct Wiener 11.34 any of the restorations based on the joint processing of the multiframes is superior to any of those based on individual processing of any single frame, and is also superior to restoration based on time-averaged image. We also found that the direct Wiener filter gave the best restoration in this example since this restoration is of high resolution and looks smooth. Another way to judge the quality of restoration is to compare the difference picture among the distortion and the restorations. Such difference pictures are given in Figure 6.24. Figures 6.24 a), b), c), d) correspond to Figures 6.21 a), 6.22 a), 6.22 b), and Figure 6.23 a), respectively. The observations certainly agree with the direct comparisons. Chapter 6. Deblurring Random Time-varying Blur 107 Figure 6.20: Image Eye Chapter 6. Deblurring Random Time-varying Blur 108 a) b) c) d) Figure 6.21: The four distorted versions of Eye Chapter 6. Deblurring Random Time-varying Blur 109 Figure 6.23: a) Individual restoration of Figure 6.21 a); b) Restoration based on time-averaged image. Chapter 6. Deblurring Random Time-varying Blur 110 a) b) -c) d) Figure 6.24: The difference pictures of example 1: a) The difference between Figure 6.20 and 6.21 a); b) The difference between Figure 6.20 and 6.22 a); c) The difference between Figure 6.20 and 6.22 b); d) The difference between Figure 6.20 and 6.23a). Chapter 6. Deblurring Random Time-varying Blur 111 Figure 6.25 shows the second example which is an 116 x 116 image "Face". The i randomly blurred images are given in Figure 6.26 where the point spread function consists of 13 x 13 pixels, and the noise statistics used are: 0\ = 0.00045, Oj '= 10.0, and A = 0.8, respectively. Restorations based on the Karhunen-Loeve transformation and from the direct re-alization of Wiener filter are shown in Figures 6.27 a), and 6.27 b), respectively. The modified Wiener filter restoration of the 1st frame is shown is Figure 6.28 a) and the restoration based on time-averaged image is given in Figure 6.28 b). Individual restora-I tion results in a very noisy image in this example. The difference pictures of frame 1 in distortion of Figure 6.26 and the three restorations of Figures 6.27 a), b) and 6.28 a) are given in Figures 6.29 a), b), c), and d), respectively. It is shown that the di-i rect Wiener filter gives a slightly better (smoother) result than the Karhunen-Loeve transformation. We list the RMS errors of this example in Table 6.7. Table 6.7: RMS errors of the 2nd example in time-varying case. RMS Frame 1 Frame 2 Frame 3 Frame 4 Blurred image 15.88 15.74 16.03* 16.12 Individual 20.74 20.26 20.94 21.02 Time-averaged image 20.91 K - L transform 13.76 Direct Wiener 12.44 \ I Figure 6.25: Image Face Chapter 6. Deblurring Random Time-varying Blur 113 c) d) Figure 6.26: The four distorted versions of Face Chapter 6. Deblurring Random Time-varying Blur 114 Figure 6.28: a) Individual restoration of Figure 6.26 a); b) Restoration based on time-averaged image. Chapter 6. Deblurring Random Time-varying Blur 115 c) d) Figure 6.29: The difference pictures of example 2: a) The difference between Figure 6.20 and 6.26 a); b) The difference between Figure 6.20 and 6.27 a); c) The difference between Figure 6.20 and 6.27 b); d) The difference between Figure 6.20 and 6.28.a). Chapter 6. Deblurring Random Time-varying Blur 116 The previous two examples suggest that using the direct Wiener filter realization may be preferable to that of the Karhunen-Loeve-transformation approach. But this suggestion is not always true. The next example shows that the direct Wiener filter may have some difficulties when dealing with images which have very sharp, regular, and clustered edges. The original image is presented in Figure 6.30, which has 120 x 120 pixels, and is a subimage of the F16 jet fighter. The image is stochastically distorted by the Gaussian shaped blur function of width 9 x 9 , with the statistics being ox — 0.0003, o2 = 5.1, and A = 0.8, respectively. The distorted multiframes are given in Figure 6.31. Again, this blurred image is restored by all the four methods mentioned above, and the results are presented in Figures 6.32 a), 6.32 b) and Figure 6.33 a), 6.33 b), respectively. Although the direct Wiener filter still gives the smallest RMS error, see Table 6.8, the visual quality (Figure 6.32 b)) of the characters U S A I R F O R C E in the image is not as sharp as that restored by the Karhunen-Loeve transformation (Figure 6.32 a)). The RMS errors for this test are summarized in Table 6.8. Table 6.8: RMS errors of the 3rd example in time-varying case. RMS Frame 1 Frame 2 Frame 3 Frame 4 Blurred image 19.86 19.54 20.32 19.75 Individual 20.74 20.31 20.16 19.89 Time-averaged image 22.32 K - L transform 17.10 Direct Wiener 14.80 Chapter 6. Deblurring Random Time-varying Blur 117 This observation is conformed by comparison of the difference pictures. The differ-ence pictures of frame 1 in distortion and the three restorations are shown in Figures 6.34 a), b), c), and d), respectively. The characters U S AIR FORCE in the difference picture Figure 6.34 c) ( restored by the direct Wiener filter) can still be seen. Figure 6.30: Image F16 Chapter 6. Deblurring Random Time-varying Blur 118 c) d) Figure 6.31: The four blurred versions of F16 Chapter 6. Deblurring Random Time-varying Blur 119 Figure 6.33: a) Individual restoration of Figure 6.31 a); b) Restoration based on time-averaged image. Chapter 6. Deblurring Random Time-varying Blur 120 Figure 6.34: The difference pictures of example 3: a) The difference between Figure 6.30 and 6.31 a); b) The difference between Figure 6.30 and 6.32 a); c) The difference between Figure 6.30 and 6.32 b); d) The difference between Figure 6.30 and 6.33 a). Chapter 6. Deblurring Random Time-varying Blur 121 6.7 Summary This chapter considered the restoration of images blurred by a random time-varying point spread function and additive noise. One feature distinguish this restora-tion problem - its multiframes nature. The restoration was based on the Wiener criterion. The model was formulated in such a way that a closed form solution to this problem was obtained in a straightforward fashion. The problem then became that of a large size problem. The problem size was circumvented. Two versions were introduced. Under the assumption of time-space-separability, it was possible to use a Karhunen-Loeve transformation to decorrelate the frames and reduce the size of the problem to that of the order of a single frame. When the noise is not necessarily separable a solution based on the direct imple-mentation of the Wiener filter was presented. This solution is valid as long as the correlation matrix of the noise has block Toeplitz structure. This is not a strict condi-tion and is valid for such a case as of space and time stationarity of the noise. Using circulant matrix approximation and FFT, whenever possible, it was shown that the order of computations is reduced to 0(S2M2 + SM2logM) in the frequency domain instead of 0(S5M6) in the space domain. The proposed filters were applied to some examples. The results were compared with those based on restoration of a single frame, and those based on time-averaged images. The results demonstrated that both of the new filters give superior results to those processed by the other two techniques. This is to be expected since the latter does not take, or does not fully take, into account the correlation of the noise amongst the frames. When the filters were applied to problems for which the conditions of the two versions were satisfied, the direct implementation version of the Wiener filter was Chapter 6. Deblurring Random Time-varying Blur 122 found to give smaller R M S errors than those using the Karhunen-Loeve transformation. The discrepancies between the results of the two filters when applied to the same problem are to be expected. These are due to the effects of approximations of the different circulant matrices involved. They are also due to the different numerical computational procedures of the two methods, e.g. the direct Wiener filter solves a system of linear equations after finding the Fourier transform of the system, while the Karhunen-Loeve transformation method involves the Fourier transform after the transformation. Chapter 7 Conclusions 7.1 General Discussion This thesis discusses the solution of the image restoration problem in which the images are distorted by a system of random linear shift-invariant point spread functions. This problem is different from the conventional image restoration problem in that the point spread function of the imaging system is random. To tackle this problem we adopted the modeling scheme which makes this problem equivalent to a deterministic blur in the presence of additive image-dependent noise. This problem differs from the usual image-dependent noise restoration problem in that here the noise is convolved with the unknown input image and not multiplied by it. The implementation of the filters required iterative solutions: we began with a guess for the unknown image and computed the noise properties; we then treated the problem as that of estimation with additive noise of known properties and found an optimum estimate for the unknown image; that estimate was then used as a better guess and so on. We generalized four linear filters, the Wiener filter, the M V U filter, the constrained least squares filter, and the geometrical mean filter, and proposed one nonlinear fil-ter, the M A P filter. We also proposed the application of the Wiener criterion to the restoration of images blurred by an optical system with random time-varying point spread function. The large dimensionality of the matrices and the vectors involved 123 Chapter 7. Conclusions 1 2 4 has been shown to cause a serious problem in the implementation of the filters. How-ever, approximation by circulant matrices, and the use of the F F T circumvented this difficulty. 7.2 Summary of Results In Chapter 2, we presented the modeling scheme for the image formation problem whose pupil function contains stochastical uncertainties. Our modeling scheme is based on the assumption that the point spread function of the optical system can be separated into two parts, the deterministic part which is the mean of the point spread function, and the stochastic part which is a zero mean random fluctuation around the mean. The random fluctuation is assumed to be a stationary process. Following these assumptions, we have converted the addressed problem into a format similar to the conventional image restoration problem. We studied the statistical properties of the model, and derived the general expressions for the correlation and the covariance matrices of the object dependent noise. We have shown that the expressions are quite complicated, but the matrices have block Toeplitz structure under practical conditions. As a special case, the correlation and covariance matrices were derived under the assumptions that the uncertainty and the additive noise encountered at the recording stage are white processes, and uncorrelated with each other. In Chapter 3, we generalized two linear filters, the Wiener filter, and the minimum variance unbiased filter to the addressed problem. The algorithm to implement these filters were iterative. The realization of the Wiener filter and the M V U filter for the one-dimensional case has been reported in the literature, our contribution here was to apply them to the two-dimensional case. It is found that both filters are better than Chapter 7 . Conclusions 125 their linear counterpart (by ignoring the uncertainty in the point spread function). The M V U filter is found to function better than the modified Wiener filter when the noise is not severe and the width of the point spread function is within a moderate range, otherwise it may run into convergence problems. Besides, the computation of the M V U filter is too expensive. On the other hand, the modified Wiener filter is more robust and computationally economical. The modified Wiener filter was also compared with the Backus-Gilbert technique. The former gave a slightly better restoration. The experimental results showed that while the information was well restored, both filters did not completely eliminate the noise effect in the restorations when the blur and/or noise effect are severe. In Chapter 4, the constrained least squares filter and a filter based on the geometri-cal mean methodology were presented. The motivation behind applying the constrained least squares technique to our problem is that of circumventing the noise effect. In order to do so, the second order difference matrix was chosen as the linear smoothing opera-tor. The application of the constrained least squares filter did yield smooth estimation but a manual adjustment of the Lagrange multiplier resulted in more pleasing visual restoration. By applying the geometrical mean technique, we Combined the modified Wiener filter and the constrained least squares filter into one filter. This filter possesses the desired effects of these two filters, and yielded good results. Chapter 5 described a nonlinear filter based on the maximum a posteriori method-ology. The reason this filter was investigated is because the random distortion problem involves nonlinear features, hence a nonlinear restoration process should provide su-perior results. Unlike the above modified linear filters the mathematical realization of this filter does not have a closed form; the M A P filter had to be realized by an opti-mization procedure. An advantage for using an optimization algorithm is that other nonlinear constraints, e.g., the positiveness of the radiant energy, may be incorporated. Chapter 7. Conclusions 126 We chose the steepest ascent approach for the optimization technique. Like all the nonlinear solvers, the implementation of the M A P algorithm involved a huge amount of computations, if done in the space domain. To avoid this, we designed an algorithm which computes the optimization procedure in the Fourier domain. The M A P filter was examined against the modified Wiener filter for the one-dimensional signal and for motion blurred images. Our conclusion is that the nonlinear M A P filter yields better restoration than the modified Wiener filter. In Chapter 6, we proposed solutions for the problem of image restoration under random time-varying blur. This problem is closely related to, but more difficult to solve than the time-invariant problem. This is because the problem of extracting a fixed object image from a continuum of differently blurred images is constrained by the availability of a large amount of data that contains only finite amount of information. Because it is impossible to process the entire continuum of the available data, some form of compression must be employed. We discussed three schemes of compression: replacing the time-varying images (l) with a single image-the average image, (2) with a single image-the time-averaged spatial correlation function, or (3) with a finite number of images (frames) representing images at pre-selected time instants. Our computer implementation was based on the data compression scheme number 3. Restoration based on this scheme presents a problem of data size. Two approaches were described. Both are based on the Wiener criterion. The first one utilized a Karhunen-Loeve-type transformation to handle the case in which the noise processes have time-space-separable statistics. Using this transformation the frames became decorrelated, and the amount of computation was eventually reduced to that of a single frame. The second one used the direct Wiener filter methodology, but the computational algorithm had to be designed to avoid the three dimensional Fourier transform. The proposed approach utilized the two-dimensional F F T and the solution of linear equations. Although the Chapter 7. Conclusions 127 computational effort of the second approach is greater than the first one using the Karhunen-Loeve deconvolution, we still managed to reduce the size of the problem considerably. The applicability of the second approach is broader than the first one. 7.3 Suggestion for Future Work In this thesis, we have developed several image restoration methods to solve the ran-dom distortion problem, and used numerical examples to demonstrate their practicality and usefulness. However more work could be done in this area. In Chapter 2, we generalized the M V U filter to the two-dimensional restoration problem using a space-frequency approach. But the solution to the Toeplitz system Tf = g has not been completely derived. Since the application of our work on the M V U filter is not limited to the random restoration problem, it is worthwhile finding a fast solution to the second half of the algorithm. We discussed the modification of the M A P filter to the problem of sensor nonlin-earity, and gave the formulae in Chapter 4. However the implementation has not been done yet. We suggest that the implementation of such a M A P filter be completed. One powerful linear filter is the Kalman filter. The Kalman filter has been intro-duced to conventional image restoration successfully, e.g., [12], [13], and [14]. Since our modeling scheme is iteratively oriented, the inherent iterative procedure of the Kalman filter may add some useful and interesting flavor to the category of random restoration. Although some preliminary work on the random restoration problem by the Kalman filter [32] has been put forward, more about this issue still needs to be explored. References [1] A. Marechal, P. Croce, and K. Dietzel, "Amelioration du contrast des details des images photographiques par filtrage des frequencies spatiales," Opt. Acta, 5, pp.256-262, 1958. [2] J.L. Harris, Sr., "Image evaluation and restoration" J. Opt. Soc. Am. 56, pp.566-574, 1966. [3] P.F. Mueller and G.O. Reynolds, "Image restoration by removal of random media degradations," J. Opt. Soc. Am 57, pp.1338-1344, 1967. [4] C.W. Helstrom, "Image restoration by the method of least squares," J. Opt. Soc. Am 57, pp.297-303, 1967. [5] J.L. Harris, Sr., "Potential and limitations of techniques for processing linear motion-degraded imagery," in Evaluation of Motion Degraded Images, U.S. gov-ernment Printing Office, Washington D.C., pp.131-138, 1968. [6] J.L. Horner, "Optical spatial filtering with the least-mean-square-error filter," J. Opt. Soc. Am 59, pp.553-558, 1969. [7] E.R. Cole, "The removal of of unknown image blurs by homomorphic filtering," Ph.D. dissertation, Department of Electrical Engineering, University of Utah, Salt Lake City, June, 1973. [8] B.R. Hunt, "The application of constrained least squares estimation to image 128 References 129 restoration by digital computer,"" I E E E Trans, on Computers, C-22, No.9, pp.805-812, 1973. [9] B.R. Frieden, "Restoring with maximum likelyhood and maximum entropy," J . Opt. Soc. A m 62, pp.511-518, 1972. [10] N.D.A. Mascarenhas and W.K. Pratt, "Digital image restoration under a regres-sion model," I E E E Trans. Circuits and Systems, CAS-22, pp.252-266, 1977. [11] W . K . Pratt, "Pseudoinverse image restoration computational algorithms" in Op-tical Information Processing Vol. II, G.W. Stroke, Y . Nesterikhin, and E.S. Bar-rekette, Eds., Plenum Press, New York, 1977. [12] J .W. Woods and C H . Radewan, "Kalman filtering in two dimensions," I E E E Trans. Inform. Theory, IT-23, pp.473-482, 1977. [13] J .W. Woods and V . K . Ingle, "Kalman filtering in two dimensions: further re-sults," I E E E Trans. Acoust., Speech, and Sig. Processing, ASSP-29, pp.188-197, 1981. [14] B.R. Suresh and B.A. Shenoi, "New results in two dimensional Kalman filtering with application to image restoration," I E E E Trans. Circuits and Systems, CAS-28, pp.307-318, 1981. [15] H . Andrews and B.R. Hunt, Digital Image Restoration, Englewood Cliffs, NJ: Prentice-Hall, 1977. [16] W . K . Pratt, Digital Image Processing, New York: Wiley, 1978. [17] A. Rosenfeld and A. Kak, "Picture Processing by Digital Computers" 2nd ed. New York: Academic, 1976. References 130 [18] V . Twerskey, "On propagation in random media of discrete scatters," R. Bellman ed., Proc. Symp. Appl. Math., Vol. 16 AMS Providence, RI 1964. [19] D. Slepian, "Least-squares filtering of distorted images," J . Opt. Soc. Am. A 57, pp.918-922, 1967. [20] L . Frank, Signal Theory, Englewood Cliffs, NJ: Prentice-Hall, 1969. [21] K. Von Heide, "Least squares image restoration" Opt. Commun. 31, pp.279-284,1979. [22] S.A. Kassam, T . L . Lim, and L . J . Cimini, "Two dimensional filters for signal processing under modeling uncertainties," IEEE Trans. Geosci. Rem. Sensing, GE-18, pp.331-336, 1980. [23] T . G . Stockham, Jr. T . M . Cannon, and R.B. Ingebretsen, "Deconvolution through digital signal processing," Proc. IEEE 63, pp.678-692, 1975. [24] J .B. Morton and H . C . Andrews, "A posteriori method of image restoration," J . Opt. Soc. Am. 69, pp.280-290, 1979. [25] 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. [26] 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. [27] R . K . Ward and B . E . A . Saleh, "Image restoration under random time-varying blur," Appl. Opt. 26, pp.4407-4412, 1987. [28] 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 131 [29] L . Guan and R.K. Ward, "Restoration of randomly blurred images by the Wiener filter" I E E E Trans. Acoust. Speech and Signal Process, in press. [30] L . Guan and R.K. Ward, "A maximum a posteriori approach to the restoration of randomly distorted signals," Proc. ICASSP-88, pp.1770-1773, 1988. [31] L . Guan and R.K. Ward, "Restoration of random blurred images by the maximum a posteriori technique," Submitted to IEEE Trans. Acoust. Speech and Signal Process. [32] A . G . Qureshi and M . M . Fahmy, "2-D Kalman filtering for the restoration of stochastically blurred images," Proc. ICASSP-88, pp.1024-1027, 1988. [33] A . V . Oppenheim and R.W. Schafer, "Digital Signal Processing," Englewood Cliffs, N J : Prentice-Hall, 1975. [34] W . F . Trench, "An algorithm for the inversion of finite Toeplitz matrix," J . SIAM Appl. Math. 12, pp.515-522, 1964. [35] H . Akaike, "Block Toeplitz inversion," J . SIAM Appl. Math. 24, pp.234-241, 1973. [36] R.R. Bitmead, "Asymptotically fast solution of Toeplitz and related systems of linear equations" Linear Algebra and Its Applications, Vol. 34, pp.103-116, 1980. [37] K . R . Castleman, Digital Image Processing, Englewood Cliffs, N J : Prentice-Hall, 1979. [38] B.R. Hunt, "Digital image processing," I E E E Proc, Vol. 63, pp.693-708, 1975. [39] H.L. Van Trees, Detection, Estimation, and Modulation Theory Vol. I, New York: Wiley, 1967. References 132 [40] H.J . Trussell and M.R. Civanlar, "The feasible solution in signal restoration" I E E E Trans. Acoust. Speech, and Signal Process. ASSP-32, pp.201-212, 1984. [41] L . Cooper and D. Steinberg, Introduction to Methods of Optimization, Philadel-phia: W.B. Saunders Company, 1970. [42] 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. [43] E . Blackman and R. Barakat, "Statistics of the optical transfer function: corre-lated random amplitude and random phase Effects," J. Opt. Soc. Am. 69, pp.544-548, 1979. [44] B . L . McGlamery, "Restoration of turbulence-degraded images," J . Opt. Soc. Am. 57, pp.293-297, 1967. [45] J .L . Horner, "Optical restoration of images blurred by atmospheric turbulence using optimum filter theory" Appl. Opt. 9, pp.167-171, 1970. [46] K . T . Knox and B.J. Thompson, "Recovery of images from atmospherically de-graded short-exposure photographs," Astrophys. J . Lett. 193, L45-L48, 1974. [47] J .H . Shapiro, "Optimum adaptive imaging through atmospheric Turbulence," Appl. Opt. 13, pp.2609-2613, 1974. [48] J . Fienup, "Space object imaging through the turbulence atmosphere" Opt. Eng. 18, pp.529-534, 1979. [49] G . J . M . Aitkin and D.L. DeSaulniers, "Restoration of atmospherically degraded images using complex spectral ratios" Opt. Commun. 28, pp.26-29, 1979. References 133 [50] P. Nisenson, A. Stachnik, C. Papliolios, and P. Horowitz, "Data recording and processing for speckle image reconstruction" Proc. Soc. Photo-Opt. Instrum. Eng. 243, pp.88-94, 1980. [51] J . G . Walker, "Object reconstruction from turbulence-degraded images," Opt. Acta 28, pp.1017-1019, 1981. [52] R. Barakat and P. Nisenson, "Finite exposure time astronomical speckle transfer function," Opt. Acta 30, pp.1405-1416, 1983. [53] B.R. Hunt and O. Ktibler, "Karhunen-Loeve multispectral image restoration, Part I: Theory" I E E E Trans. Acoust. Speech, and Signal Process. ASSP-32, pp.592-600, 1984. [54] N.P. Galatsanos and R.T . Chin, "Digital restoration of multi-channel images," Proc. ICASSP-87, pp.1244-1247, 1987. [55] N.P. Galatsanos and R.T . Chin, "Restoration of multi-channel images by Kalman filtering," Proc. ICASSP-88, Section 41.M10.18, 1988. [56] D .C . Ghiglia, "Space-invariant deblurring given N independently blurred images of a common object" J . Opt. Soc. Am. A 1, pp.398-402, 1984. [57] K. Fukunaga, Statistical Pattern Recognition, New York: Academic, 1972. [58] G . Stewart, Introduction to Matrix Computing, New York: Academic, 1973. [59] E . L . Hall, Computer Image Processing and Recognition, New York: A c a d e m i c , 1979. Appendix A Frequency Domain Computation of the M V U Filtering The formulation of equation (3.73) requires the inversion of the covariance matrix Kn{f), and the computation of the triple matrix multiplication T = H T K~1(f)H (A .221) and the right hand side g = H T K - \ f ) g (A .222) Space domain computation is too expensive, hence we approximate the matrices by circulants, and compute them in frequency domain. The circulant approximation of T is Te = H j K ^ { f ) H c (A .223) where the subscript c denotes the circulant approximation of a matrix; AH and AKn are diagonal matrices, in which the diagonal elements are the eigenvalues of H c and Kn(f), respectively, and 7M3 is the two-dimensional discrete Fourier matrix. Similarly for g, gc = HjK-el[f)Hc ' (A.224) •= A # 7M i 7M3 A*-* 7M7 g = J A ^ A ^ A ^ 134 Appendix A. Frequency Domain Computation of the MVU Filtering 1 3 5 Further Fourier domain manipulation yields f = T~lg (A.225) Therefore the filter become an inverse filter, and useful information in Kn(f) is canceled. In order to avoid the cancellation of Kn(f), we shall obtain the approximated version of T from Tc, and g from gc, and then solve Tf = g in the space domain. By doing so, the computation of forming T and g is reduced from 0 ( M 6 ) to 0 ( M 2 log M ) . For a linear shift-invariant imaging system, the matrix H is block Toeplitz with structure ' Hx H2 Ht 0 • H2 H = H2 Hr (A.226) Appendix A. Frequency Domain Computation of the MVU Filtering 136 where • hi,i hj,r 0 h),r for j = 1,2, • • • M, and r = M — N + 1. For the sake of Fourier domain computation, matrix H is approximated by a block circulant Hc, Ha Hci Ha (A.227) HCr H. c2 Hc H, 0 HC2 Hcl HCr Hcr_i . Hci where Hcj = [Hj | Hja Appendix A. Frequency Domain Computation of the MVU Filtering 137 hj,i hj,2 hjti 0 k3,r h3,r h3,r is the circulant approximation to submatrix Hj. Inserting (A227) into the expression HjK^(f)Hc yields H j K - c \ f ) H e = B T K ^ ( f ) B | B * K £ { f ) B . R T a K ^ { f ) H | R T a K ^ { f ) R a (A.228) In equation (A.228), the upper left submatrix H T K n e 1 ( f ) R contains the useful infor-mation about T. Rewriting H in (.4.227) as Hi Hia H 2 Hia Hi H X a 0 H 2 H 2 i x • H = H r H r Hr H r " ra Hi Hu H2 H2a Hi H^ (A.229) Appendix A. Frequency Domain Computation of the MVU Filtering 138 which can be equivalently written as H = [Hi Hu • • • H M HJ Ma (A.230) Sustituting (A230) into # T #"-(/)#, the upper left submatrix of HjK-cl(f)He results in HlK^(f)Hi HfK;}(f)Hia HlK-HDHi KK-c\f)Hia HMK~cl(f)Hi HMK^(f)Hla HMaK-c\f)Hi HMaK^{f)Hu From (A.231) we form a matrix , HTiK-e\f)Hi HTK-}{f)H = ' (A.231) BfK~l (f)HM HiK~el (f)HMa HlK^ (f)HM H\aK^ (f)HMa HMK-el(f)HM HMK~cl(f)HMa HMaTK^l(f)HM HMaK^(f)HMa (A.232) HLK~cl(f)Hi ••• HMK-cHf)HM where the (t,j')th submatrix in (A.232) is the (2t — 1,2/ — l)th" submatrix in (A.231). Equation (.A.232) is the approximation to the original T matrix. Following a similar procedure, the approximation to the right hand side g is formed by ' HjK-}{f)g HTK^(f)g = HlK-}[f)g (A.233) HTMK-}{f)g The approximation to T in (A.232) is an block Toeplitz matrix which can be solved by the Toeplitz solver designed by Akaike [35]. Appendix B Computing VWi In order to find the explicit expression for VWj , we first investigate the differentia-tion of the logarithm of the determinant of a matrix function In |A(x)| with respect to its sole variable x. Here all the entries in the matrix A(x) are assumed to be functions of x. It is easy to see that d\n\A{x)\ _ 1 d\A(x)\ dx ~ \A(x)\ dx * • 4> Since |A(x)| is a function of a,y ( the ij th. entry of the matrix A{x)), the chain rule gives ~dx~ = ^ ~daT~dx~ ( B 2 3 5 ) Expanding |A(x)| according to Laplace's expansion theorem, it yields |A(x)| = £ > * A ; * (B.236) = aijAij + J^^Aik (B.237) where Aik is the cofactor of a,-*. Since Aik is not a function of a t J when j ^ k, it yields, = -£-1^ + J2aikAik] = Ay (B.238) Therefore equation ( £ . 2 3 5 ) can be simplified to 139 Appendix B. Computing VWi 140 Substituting (5.239) into (5.234), and considering that r£fy\ is the jtth entry of A(x)~x, the inverse matrix of A(x), it results in din |^4(z)| d a i > If A(x) is symmetric, so is A(x)~1, and (B.240) can be written as (B.240) dx dx (B.241) Since the /th entry in vector VW\ is dIn\Kn(f)\/dfi, straightforward application of equation (5.241) yields dln\Kn(f)\ = £ £ [ t f ; i ( / ) ] / t * " W ] ' > dh Therefore we obtain the derivative VWX ( 3[K„(/)|,y ^ 3/l v 8/M y (B.242) (B.243) which is equation (5.120). Appendix C Computing VW2 Since Kn(f) is a function of / , the computation of VW2 involves the differentiation of K~1(f) with respect to / . Therefore W i = m-BWm-m (c.244) df = -2HTK-\f){g - Hf) -(g- i f / ) T ^ ( / ) ( g - Hf) (C.245) The algorithm to compute a K ^ a } ^ can be demonstrated by the differentiation of the inverse of a matrix function A(x) with respect to its sole variable x. The procedure starts with differentiating the identity equation A[x)A[x)-l = I (C.246) Taking derivatives on both sides of (C.246) with respect to x yields d[A(x)A{x)~1 dx = 0 (C.247) Expansion of (C.247) results in **®-A[zr* + A { x ) d - ^ = 0 (C.248) dx dx Therefore * ^ = - A { x ) - i * m A [ x ) - i (c.249) dx dx Normally, A[x)~x is very difficult to obtain analytically, hence (C.249) provides a so-lution to compute the derivative of A(x)~1. 141 Appendix C. Computing WV2 142 It is trivial, from (C.249), to identify that where / , is the /th member of the image / . Inserting (C.250) into (C.245) yields the expression for the /th entry in VW2 [ W 2 ] , = -2\HTK-\f){g - if/)], - (g- Hff K ~ \ f ) d - ^ 1 K-\f){g - Hf) (C.251) which is equation (5.122). Appendix D The Alternative Expression of the Wiener Filter In this appendix, we show that the Wiener filter in equation (3.50) can be equiva-lently written as equation (6.182) so that the Karhunen-Loeve-type transformation is relevant to apply to the the restoration of images distorted by a time-varying point spread function. We first rewrite the Wiener filter in equation (3.50): L = RfHT[Rn + HRfHT}-1 (D.252) By inversion lemma, we know that + HRfH7)-1 = R-1 - R~lH[HTR~lH + R]l\~lHTR-1 (D.253) as long as Rn and Rj are of full rank, and H is of full column rank. The block Toeplitz structure of H guarantees that H is of full column rank. Since i?» and Rj are correlation matrices of the noise n and the object / , respectively, the full rank condition is normally satisfied. Substituting (Z?.253) into (D.252) yields L = RfHT{R-1-R-1H[HTR-1H + Rj1}-1HTR-i} = RjHTR~l — RfHT R~lH[HT R~lH + Rj1]"1 HT R~l = RjHTR? - Rj\HTRZlH + Rj1 - R-fl)\HTR-lH + Rjl}~1HTR-1 = RjHTR~l - Rf\HTR~lH + R]l\\HTR~lH + R]l]~lHTR-1 +RfRj1[HTR~1H + R]l\-lHTR-x = \HTR-lH + R]l]'lHTR-1 (D.254) 143 Appendix D. The Alternative Expression of the Wiener Filter 144 The last equation is easily identified as the filter in equation (6.182) if we replace Rn by Rn(f).
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Restoration of randomly blurred images
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Restoration of randomly blurred images Guan, Ling 1988
pdf
Page Metadata
Item Metadata
Title | Restoration of randomly blurred images |
Creator |
Guan, Ling |
Publisher | University of British Columbia |
Date Issued | 1988 |
Description | This thesis considers the restoration of images blurred by random point spread functions. The stochastic point spread functions are represented by the sum of a deterministic blur function plus an image-dependent noise term. The stochastic properties of such systems are studied. Several linear filters are derived. These filters are based on the following criteria: the Wiener, the minimum variance unbiased, and the constrained least squares. A filter based on the geometrical mean methodology is also presented. The popular Wiener filter, and the minimum variance unbiased filter have already been applied to this problem for the one-dimensional case. We generalize these filters to the two-dimensional case. These filters are found to give good restoration results when the blur and/or the additive noise effects are not severe. When either or both of the above conditions are not satisfied, numerical instability may occur. The constrained least squares filter does not suffer from such ill conditioning. However, the performance of this method needs to be manually controlled in order to yield good restoration results. The geometrical mean filter is formed by combining the modified Wiener filter and the constrained least squares filter. This filter gives the best restoration amongst the four linear-based filters because it does not have the shortcomings the other filters have. A nonlinear filter, the maximum a posteriori filter is also derived and implemented. Since the random blur effect is a nonlinear stochastic process, we expect the nonlinear filtering scheme to simulate the reverse of the blurring process more closely. The superiority of this filter over the modified Wiener filter is demonstrated by examples. In addition, the restoration of random time-varying blurred images is also discussed. The Wiener criterion is adopted. The solution of this problem involves the processing of a huge amount of multi-image data. To circumvent this problem we present two approaches to realize this Wiener filter. One is based on the Karhunen-Loeve transformation, and the other is a direct Wiener filter implementation. These two approaches give better results than that given by processing any single frame individually. Since image restoration involves computation of very large systems, The computational effort is an important factor to consider in the implementation of all restoration algorithms. By means of the circulant matrix approximation and the fast Fourier transform, our methods are implemented in the frequency domain. The computational efforts are reduced significantly. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-10-11 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
IsShownAt | 10.14288/1.0098294 |
URI | http://hdl.handle.net/2429/29108 |
Degree |
Doctor of Philosophy - PhD |
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_1989_A1 G82.pdf [ 11.38MB ]
- Metadata
- JSON: 831-1.0098294.json
- JSON-LD: 831-1.0098294-ld.json
- RDF/XML (Pretty): 831-1.0098294-rdf.xml
- RDF/JSON: 831-1.0098294-rdf.json
- Turtle: 831-1.0098294-turtle.txt
- N-Triples: 831-1.0098294-rdf-ntriples.txt
- Original Record: 831-1.0098294-source.json
- Full Text
- 831-1.0098294-fulltext.txt
- Citation
- 831-1.0098294.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0098294/manifest