Frequency-domain adaptive iterative image restoration and evaluation of the regularization parameter Moon Ci Kang Aggelos K. Katsaggelos, MEMBER SPIE Northwestern University McCormick School of Engineering and Applied Science Department of Electrical Engineering and Computer Science Evanston, Illinois 60208-3118 E-mail: [email protected] Abstract. An important consideration in regularized image restoration is the evaluation of the regularization parameter. Various techniques exist in the literature for the evaluation of this parameter, which depend on the assumed prior knowledge about the problem. These techniques eval- uate the regularization parameter either at a separate preprocessing step or by iterating based on the completely restored image, therefore requiring many restorations of the image with different values of the re- gularization parameter. The authors propose a nonlinear frequencydomain adaptive regularized iterative image restoration algorithm. According to this algorithm a regularization circulant matrix is used that corresponds to the assignment of a regularization parameter to each discrete frequency. Therefore, each frequency component is appropriately regularized. The resulting algorithm produces more accurate resuits and converges considerably faster than the algorithm that uses one regularization parameter for all frequencies. The regularization matrix is updated at each iteration step, based on the partially restored image. No prior knowledge about the image or the noise is required. The development of the algorithm is based on a set-theoretic regularization approach, where bounds on the weighted error residual and stabilizing functional are updated in the frequency domain at each iteration step. The proposed algorithm is analyzed theoretically and tested experimen- tally. Sufficient conditions for convergence are obtained in terms of a control parameter, which is specified by the conditions for convergence and optimality of the regularization parameters at each discrete frequency location. Finally the proposed algorithm is compared experimentally with other related algorithms. Subject terms: digital image recovery and synthesis; image restoration; regularization; iterative algorithms. Optical Engineering 33(10), 3222-3232 (October 1994). 1 Introduction A number of techniques providing a solution to the image restoration problem have appeared in the literature.1'2 A clasDueto the imperfection of physical imaging systems and the sification of the major restoration approaches can be found, particular physical limitations imposed on every application for example, in Chapter 1 of Ref. 2. Most existing image where image data are recorded, a recorded image represents a noisy and blurred version of an original scene in many practical situations. The blur may be, for example, due to motion, incorrect focusing, or atmospheric turbulence, while the noise may originate from the image formation process, the transmission medium, the recording process, or any combination of these. The goal of image restoration is to recover the original scene from the degraded observations. Image restoration techniques are oriented toward modeling the degradiations, blur, and noise, and applying an inverse procedure to obtain an approximation of the original image. Paper DIR-Ol received Mar. 2, 1994; accepted for publication June 27, 1994. 1994 Society of Photo-Optical Instrumentation Engineers. 0091-3286/94/$6.OO. 3222/OPTICAL ENGINEERING/October 1994/Vol. 33 No.10 recovery or restoration methods have a common estimation structure, in spite of their apparent variety. This common structure is expressed by regularization theory. Qualitatively speaking, the underlying idea in many regularization approaches is the use of prior information on the solution or the parameters to be estimated. This prior information is combined with the data information and defines a solution by trying to achieve smoothness and yet remain ''faithful" to the data. The parameter that controls the trade-off between fidelity to the data and smoothness of the solution is called the regularization parameter. Its selection is very critical to the overall quality of the restored image. Various approaches exist for evaluating this parameter. Some of these require FREQUENCY-DOMAIN ADAPTIVE ITERATIVE IMAGE RESTORATION knowledge of the noise variance,35 and some do not (crossvalidation,4'6 maximum likelihood4). All of these approaches either require the completely restored image based on a value of the regularization parameter—in order to appropriately adjust the new value of the regularization parameter (resulting in considerable computational overhead)—or treat the evaluation of the regularization parameter as a completely separate initial step. In Refs. 7 and 8 iterative algorithms were introduced and analyzed that evaluate the regularization parameter simultaneously with the restored image, thereby removing the disadvantages of excessive computational load and/or the need for knowledge of the noise variance of existing techniques. Another disadvantage of existing nonstochastic regularization techniques is the use of a single scalar regularization parameter. Its selection, although it is based on an optimality criterion, represents a compromise between regions that need turbation in the solution. A regularization approach replaces an ill-posed problem by a well-posed problem whose solution is an acceptable approximation to the solution of the given ill-posed problem.'2 to be overregularized and regions that need to be under- (3) regularized. Spatially adaptive techniques provide a solution to the restoration of images that are nonstationary in nature, by still using a single regularization parameter but adapting 91 the smoothness In this paper we introduce adaptivity into the restoration process by using a constant smoothness constraint, but assigning a different regularization parameter at each discrete frequency location. We can now ' 'fine-tune' the regularization of each frequency component, thereby achieving more ' desirable restoration and at the same time speeding up the convergence of the iterative algorithm used for obtaining the restored image. We evaluate the regularization parameters simultaneously with the restored image (based on the partially restored image), by extending the results in Refs. 7 and 8. The derivation of the algorithm is based on a set-theoretic regularization approach, as explained in Sec. 2. The convergence analysis of the proposed iterative algorithm is presented in Sec. 3. The analysis of the singular points of the blurring function is shown in Sec. 4. Finally, experimental results and conclusions are given in Secs. 5 and 6, respec- (2) XEQX where Q is a set of signals with certain known properties. The knowledge we want to express by constraining x to Q is that the solution is smooth. This is achieved by defining Q by Q={xI fCxfI2E2} where (I.f( denotes the 12 norm, E is a prescribed constant, and C is a linear operator that has high-pass characteristics. Similarly, the prior information constrains the zero-mean noise n to an ellipsoid Q, given by nnQ{n( 11n112<E2} (4) where E is an estimate on the data accuracy (noise norm). Since n lies in a set, it follows that a given observation y combines with the set Q to define a new set that must contain x. Thus the observation y specifies a set Q that must contain x, i.e., XEQ,dY{X(y—DX)EQfl} (5) Since each of the sets Q, and Q,, contains x, it follows that x must lie in their intersection. The center of one of the ellipsoids that bounds their intersection is chosen as a solution.'°'1' It is given by (DTD+oCTC)x=DTy, tively. (6) (/)2 2 Frequency-Adaptive Regularization 2.1 Formulation The image degradation process can be adequately modeled by a linear blur and an additive white Gaussian noise process.' Then the degradation model is described by y=Dx+n With the set-theoretic approach,'°"3 the prior knowledge constrains the solution to certain sets. The sets to be used are ellipsoids. Therefore, consistency with all the prior knowledge pertaining to the original image serves as an estimation criterion. More specifically, it is assumed that (1) where the vectors y, x, and n represent, respectively, the lexicographically ordered noisy blurred image, the nonstochastic original image, and the additive noise. The matrix D represents the linear spatially invariant distortion. The image restoration problem calls for obtaining an estimate of x given y, D, and possibly some knowledge about the noise process. If the image formation process is modeled in a continuous infinite-dimensional space, Eq. (1) becomes a Fredholm integral equation of the first kind. Then the solution of Eq. (1) is almost always an ill-posed problem.'2 This means that the unique least-squares solution of minimal norm of Eq. (1) does not depend continuously on the data, or that a bounded perturbation (noise) in the data results in an unbounded per- where ci = is the regularization parameter. A posterior test is actually required to make sure that indeed the solution of Eq. (6) satisfies both Eqs. (3) and (4). As already mentioned in the introduction, in this work we propose an adaptive form of Eq. (6) by using a regularization matrix instead of a scalar regularization parameter a. More specifically, we propose to use the following two ellipsoids Q and Q, in a weighted space: Q={xJ JICxIJRER} (7) Q={x Iy—DxE} (8) where P and R are both block-circulant weighting matrices. Then a solution that belongs to the intersection of Q and is given by (DTPTPD+XCTRTRC)iz= DTPTPy, (9) where X = (p/ER)2. Let us define P Tp B, R = P0, and XQ TQ A. Then Eq. (9) can be written as (D TBD + CTBAC)i = DTBy, (10) OPTICAL ENGINEERING I October 1994 I Vol. 33 No. 10 I 3223 KANG and KATSAGGELOS or X0(i)=p(i)D*(i)Y(i) B(DTD+ACTC)=BDy , — (11) since all matrices are block-circulant and they therefore corn- Xk+ 1 = mute. The regularization matrix A is defined based on the set-theoretic regularization as A= My— Dx112(fICxfI2I + A) (12) is a block-circulant matrix used to ensure convergence, and I is the identity matrix; B plays the role of the ' ' 'shaping' matrix'4 for maximizing the speed of convergence at every frequency component, as well as for cornpensating for the near-singular frequency components, as will be described in the following sections. where 2.2 Proposed Iterative Algorithm There are a number of ways for solving Eqs. (6) and (10) D(rn)Xk(rn)2 — 13(/)mI() flIC()Xk(2 + ik(l) C(i){2) Xk(l) (1 +13(i)D*(i)[Y(i)_D(i)Xk(i)] (17) According to this iteration, since 13 and Sk are frequency dependent, the convergence of the iteration can be accelerated, as will become clear from the analysis in the following sections. 3 Optimality and Convergence In this section we first determine the allowable range of each regularization and control parameter in the discrete frequency domain. This is achieved by considering the minimization of when the regularization parameter is known."2 For example, E[IIx — i(A)112] a direct solution in the discrete frequency domain can be where E['] denotes the expectation operator, and (A) the obtained, since the matrices D and C are assumed to be block circulant. However, if the regularization parameter is not known, and also for incorporating prior knowledge about the original image into the solution, iterative algorithms can be employed. Iterative algorithms have been analyzed in detail in Refs. 10 and 1 1 when a, or both bounds E2 and E2 in Eqs. (3) and (4), respectively, are known, when only one bound E2 is known,7 and when no bound is known.8 solution of Eq. (1) as a function of A. Following a procedure similar to the one described in Ref. 8, we obtain the following optimal range of a(l) and i(l) at every discrete frequency location (see Appendix): = iIY() — D(rn)(m)I2 — <a N2 max1(IX(i)121C(i)12) In order to obtain a solution to Eq. (10), successive approximations are used in Refs. 10, 15, and 16, resulting in Xk± I Xk + B[DTy (DTD + AkCTC)xkl ' (18) N2 1Y(m) <— N2 — D(m)(m)J2 where Ak= IIY — DxkII2[fICxkJI2I + kI '. It is mentioned here that the iteration (13) can be also derived from the regularized and equation 0 < kopt(l)N2max (Xk(l)2C(l)2) (DTD + ACTC)i = DTy , (14) by using the generalized Landweber's iteration.'4 Since all matrices in the iteration (13) are block-circulant, the iteration can be written in the discrete frequency domain as X0(i) = p(i)D*(i)Y(i) , Xk+ '(P Xk(l) + p(i)[D*(i)Y(i) (15) where = ('1"2)' 0l1N— 1, 0l2N— 1, Xk+ 1(1) and Y(i) represent the two-dimensional (2-D) DFT of the unstacked image estimate Xk± and the noisy blurred image y, and D(i), C(i), 13(l), and ak(l) represent 2-D DFTs of the 2-D sequences that form the block-circulant matrices D, C, B, and Ak, respectively. Since k is block-circulant, ak(l) is given by I \2 m ') I'?!\ ak(l) ' ak) flIC(n)Xk( + Vf " km) — — — — C(m)Xk(m)2 = pt (20) , where a0(L) and akopt(l) are the values that minimize (18). We turn now to the analysis of the convergence of the frequency-domain adaptive algorithm. The algorithm has more than two fixed points. The first fixed point is the inverse —(ID()I2+ak()IC()I2)Xk(PI 'V (19) min1(JX(/)121C(i)12) (13) ''16) - where k(L) is the 2-D DFT of the sequence that forms k• Therefore, the proposed iteration takes the form 3224 / OPTICAL ENGINEERING / October 1994 / Vol. 33 No. 10 or generalized inverse solution of Eq. (1). Clearly, if D is invertible, then X(i) = Y(i)/D(i) is a fixed point of the iteration (17), as is easily verified. If D is not invertible, then solutions that minimize y — Dx112, that is, satisfy D(i)12 X(i) = D*(i)Y(i), are also fixed points of the iteration (17). For the frequencies for which D(i) = 0, the solution X(l) = 0 is obtained, which is the generalized inverse solution The second type of fixed points are regularized approximations to the original image. Since more than one solution exists to the iteration (17), the determination of the initial condi- tion becomes important. As shown in Sec. 5, it has been verified experimentally that if a ' 'smooth' ' image is . . . used for X0(i), almost identical fixed points result, independently of X0. Sufficient conditions for the convergence of the iteration (17) are derived now, based on the assumption that Xk(l) is in the vicinity of a fixed point and on the validity of the FREQUENCY-DOMAIN ADAPTIVE ITERATIVE IMAGE RESTORATION - 2jY(n) — D(n)Xk(n)I2C(l)I4Xk(l)I2 resulting approximations. The frequency iteration (17) can be rewritten as Xk+ 1(f) = Xk(l) + 3(i)[D*(i)Y(i) (21) —ID()2Xk()—F(Xk())I where the frequency nonlinear factor F (Xk(l)) is equal to F(Xk(l))= m'(m) — D(rn)Xk(rn)2C(l)2Xk(l) [rnC(rn)Xk(rn)2 + k(1)I : Y(n)Hl—Dl2+ k Xk± () — Xk(l) = [1 — 3(l)D(l)2I[Xk(l) — Xk_ ()I - 2Y(n)—D(n)Xk(n)2C(l)I4Xk(1)I2 + k(1)1 . — 1 (01 + 0(h2) , where Mk(l) = 2/Hk(l). A lower bound for Mk(l) independent of the iteration step k is given by 2 ' (30) M(i) — D(i)2 + (24) where min(L) 5 the minimum value over all iterations at a frequency I (also discussed in Sec. 5). The bound M(l) is where F '(Xk(l)) is the derivative of F(Xk(l)) with respect to Xk(l), and 0(h2) is the second-order zero term ofthe iteration step Xk(l) — Xk_ (L) at frequency . In Taylor's expansion (24), we have assumed that the factor Xk(l) — Xk_ 1(L) is very small (first-order zero function); therefore the term Y2[Xk(l) — Xk_ l(l)]2F"(Xk(l)), where F' denotes the second derivative of F, has been omitted. Using Eq. (24), Eq. (23) becomes Xk+ 1(1) — Xk(l) [1 - \\ — 2,jY(n) — strictly positive, since 0D(i)I1 and 0C(i)lE1, and therefore a f3(l) satisfying 0 < 3(l) < M(l) can always be found at every frequency component. The condition Hk(l) > 0 is used in establishing a bound Ofl k(L); it can be rewritten as fD(l)2k(l)2+ = [ r - D(n)Xk(n)2C(i)2]k(i) / \2 + [D(l)2(fC(rn)Xk(rn)2) +fY(n) + k(1)) (25) — If we consider the magnitudes of both sides of Eq. (25), we obtain with the use of the triangular inequality — [2D(l2C(rn)Xk(rni2 + IYQ) rn C(rn)Xk(rn)2 + k(') X[Xk(l)—Xk_J(l)I . (29) , (23) F(Xk(l)) = F(Xk_ 1 (1)) + F '(Xk(l))[Xk(l) - 2 Y(n) - D(n)Xk(n)2IC(i)2C(i)Xk(i)2] . Xk(l) — .,jY(n) — D(n)Xk(n)2C(l)2 1— (28 [C(rn)Xk(rn)I2 O<13(l)<Mk(l) According to the Taylor series expansion, () -- D(n)Xk(n)2C(l)I2 - should be strictly positive, and 3(l) should satisfy —3(l)[F(Xk(l))—F(Xkl(l))J Xk± — jC(rn)Xk(rn)12 + k(l) Rewriting the iteration (21) for two consecutive values of k, we obtain (27 is sufficient for the convergence of the iteration. In order for the inequality (27) to be satisfied, Hk defined as (22) . < The (31) condition for k(l) is m(m)Xk(m)2 + k(') 2Y(n) — — 21 D(l)12 .jC(rn)Xk(rn)l2 — IY()— D(n)Xk(n)121C(l)12 (rnC(m)Xk(m)2 + k(1)) XXk(l)—Xk_l(l) . ' (26) C(l)l [ / 21D(i)12 + 21D(1)12[ C(i)l2 Y(n) D(n)Xk(n)12)2 According to (26) the condition 1 °'P'/D'l2+ Y(n) — D(n)Xk(n)f2IC(l)2 — —— — C(m)Xk(m)I2+k(l) V2 + 8lD(l)l2lC(l)Xk(l)l2lY(n) — D(n)Xk(n)12] COflV(/) . (32) OPTICAL ENGINEERING / October 1994 / VoL 33 No. 10/3225 KANG and KATSAGGELOS Overall, if 0 < 13(1) <M(l) and Eq. (32) is satisfied, the iteration (17) converges. Henceforth, we denote by sed(1) the k(l) used in the iteration (17). 1.6e+13 1.4e+13 1.2e+13 4 Analysis of Singular Points of D(i) le+13 Two conditions have been derived for k(l): Eq. (20), which is based on an optimality criterion, and Eq. (32), resulting from the convergence analysis of the iteration. If both conditions are satisfied, then (Sk(L$e+12 6e+12 4e+12 (33) 2e+12 A used(1) satisfying (33) can be found for almost all discrete frequencies L. However, the lower bound 8"(i) becomes 0 COflV(l)<USed(l)<OPt . 5 greater than the upper bound Pt for frequencies 1 in the neighborhood of the singular or near-singular points of (a) 2.5e+13 D(i), due to the division by D(i) in Eq. (32). The behavior of the algorithm at three types of frequency regions is examined next. 4.1 10 Nuiiiber of iterations 2e+13 Nonsingular-Point Region 1.5e+13 This is the region defined by the discrete frequencies 1 for which '5k (Li le+13 COflV(l)<SOPt . (34) In this case a ''(i) was used, as defined by the relation ased(i) = COflV(i) + .y[Pt — oflv(i)I , (35) 0 10 0< y < 1. The values of tised and mjfl for an arbitrary frequency location 1= (O.71T, O.75'rr) are shown in Figs. 1(a) and 1(b) for SNRs of 20 and 30 dB, respectively. Such a frequency location belongs to the nonsingular-point region. where 5e+12 Number of iterations (b) Fig. 1 Comparison of sed, opt and min with a uniform blurred image (7 x 7) for: (a) SNR = 20 dB at an arbitrary frequency component and (b) SNR = 30 dB at an arbitrary frequency component, at an arbitrary frequency (O.7-rr, O.75'rr). 4.2 Exact-Singular-Point Region Let us now consider the frequencies 1 * for which D(i * ) 0. and therefore ak(1 * ) 0. At these frequencies k(1 'K ) It is this relationship that defines the near-singular-point re- Since according to Eq. (17) we have X0(i) = 3(i)D *(i) gion. In this region, unlike the exact-singular-point region, x is not zero. It is instead very close to the inverse solution, Y(i), it holds that X0(i*) 0. Furthermore, according to the iteration (17), since ak(l) is very small due to the very large value of Xl(l*)=Xo(l*)+13(l*)[D*(l*)Y(l*)_cxk(l*)Xk(l*)] In this region, in order to avoid the inverse solutions, spectral filtering is used. By setting KTK= (DTD + ACTC) and D T = KTg, the Moore-Penrose pseudoinverse solution17 of Eq. (14) is written as = X0(i *) = 0 ,' X(l*)=X0(l*)=0 (36) From the above analysis we see that at the set of frequencies 1*, the generalized inverse solution is obtained, that is, X(i*) = 0, according to Eq. (36). As is clear from Eq. (32), for values of D(i) close to zero the resulting values of oflV(1) become very large. For these frequencies it turns out that Pt (38) where cr , v, , and u denote respectively the singular value of K and the left and right singular vectors of K, and p is the rank of K. The spectral filtering solution is then expressed as14'8 4.3 Near-Singular-Point Region COflV(i) > =:—'—-. i=1°j XSF = s(o) -'- v1 , (39) where s(cr) = 1 if cry> T and 0 otherwise. In our work we (37) 3226 / OPTICAL ENGINEERING / October 1994 / Vol. 33 No. 10 propose to determine 'r by the condition (37). With this simple FREQUENCY-DOMAIN ADAPTIVE ITERATIVE IMAGE RESTORATION spectral-filtering function, the frequency components in the near-singular-point region are treated the same as with the exact singular points. We now include the spectral filtering function B in Eq. (14) following a procedure similar to Strand' 14 Returning to the iteration (13), the restored result at the k'th iteration with zero initial condition is given by Xk= I BKTK =iLrI_(I_BKTK)k]U1Tg (40) Here, B is restricted to have the same set of eigenvectors as KTK, that is, (41) Bv1==p,v, where X1 = a- and we assume that O<p1 X1 <2 for all i. Furthermore, B = F (KTK), where F (X) is a rational fraction and 0< XF(X) <2, for all X. Then Eq. (40) becomes Xk= P T i=1 a-i UT =Rk(Xj)-—-vj (T (42) Equation (43) has the form of Eq. (39); in other words, Rk (X1) should approximate the spectral filtering function s(a-). There are a number of ways to approximate the step truncation function of s(o-1) with the polynomial X[1 — (1 —pX)'i. One way is the minimization of the ripple energy in both the stop band and the passband. In other words, define a polynomial of order m, P(X)=XF(X)=a1X+a2X2+..- +amXm (45) and minimize the ripple energy rT E(T) j P2(X) dX + j [1 — P(X)12 dX 0 -r (47) In order for (l) to satisfy the sufficient condition for convergence of the iteration (29), the values of M(l) in Eq. (30) need to be determined, so that the condition 0 < 3(l) <M(l) is verified at each iteration step. For the determination of M(l), the minimum value of 1sed(1) over all iterations, denoted by min(L)' needs to be known, where sed(1) is defined in Eq. (32) and y = 0.5 is used. Toward that end, the value of s(i) estimated from Eq. (32) is set equal to min(L). If c(L) < ised(1) then l) is set equal to min(L). In all our experiments, it was observed that min(L) was reached after a small number of iterations. In Figs. 2(a) and 2(b), the plots of 13(1) and M(l) — 3(l) are shown, respectively. Since M(l) — 13(1) is positive at all frequencies, 3(l) is less than the convergence bound M(l) at all frequencies, and therefore the sufficient condition for convergence is satisfied. The performance of the restoration algorithm was evaluated by measuring the weighted improvement in signal-tonoise ratio defined in decibels by — SNR10 log10 (44) (46) Therefore, B is a polynomial of D TD with coefficient a, 1443.75D(i)4— 34651D(i)16 + 4504.5D(i)I8 — 30031D(i)f'° + 804.375ID(i)'2 (43) where Rk(X)= 1 —(1 _pX)k = 31.5 — 315ID(i)2+ I Xorig 2 Cong p X0fl11 (48) where Xong denotes the original image, which is available in a simulation, and the restored image at convergence. The weighting matrix P is defined in Eq. (9) and satisfies pTp B [an example of B is shown in Fig. 2(a)]. Since P is a high-pass filter, the high-frequency content of the flumerator and denominator quantities in Eq. (48) is weighted heavier. This is in agreement with the objectives of the adaptive filter we propose. We believe that the metric in Eq. (48) better represents the quality of the image restored by the proposed and nonadaptive algorithms. The values of LNR and the required number of iterations, as well as the estimation of the variance of the noise, a- = y — DxklI2/N2 at convergence, are shown in Table 1 for SNRs of 20 and 30 dB with a uniformly blurred image, for the proposed algorithm as well as for the algorithm in Ref. 8, which uses one regularization parameter for the whole image. From this cornparison the proposed algorithm reaches convergence considerably faster than the algorithm in Ref. 8, and also results in obtained as outlined above. The detailed procedure of deriving the spectral filtering function is shown in Ref. 14. a sharper image in terms of the metric LNR and according to visual inspection. The values of the residual y — DxkII2/N2 and of the error 5 Experimental Results between iteration steps, IIXk — Xk_ 112, for the two SNRs men- The performance of the proposed frequency-domain iterative restoration algorithm is illustrated with artificially blurred images obtained with a 7 X 7 uniform point-spread function (PSF). Some of these results are presented in this section, where a 256- x 256-pixel portrait image was used, the 2-D — Laplacian was used for C, and the criterion Xk112/ 106 was used for terminating the iteration; 3(l), IIXkII2 which is obtained from the analysis of the singular points with T = 0.007, is given by'4 13(i) = F (ID()I2) tioned above are compared in Figs. 3 and 4, respectively. Notice that the plots are scaled appropriately and that the vertical axis in Fig. 4 is logarithmic. It is observed in Fig. 3 that in all cases the residuals converge to a value close to but slightly smaller than the true values. This result is expected, since the solution at convergence [center of an ellipsoid bounding the intersection of Q and Q in Eqs. (3) and (5), respectively] lies inside the ellipsoid Q,,. The plots in Figs. 3 and 4 demonstrate that in all cases the iterative algorithm converges at the same rate. The noisy blurred images are shown in Figs. 5(a) and 5(c) for SNRs of 20 and 30 dB. The corresponding restored images are shown in Figs. 5(b) OPTICAL ENGINEERING/October 1994/Vol. 33 No. 10/3227 KANG and KATSAGGELOS (a) 4 2.5)<10 4 2.0>(10 x1 0 0 4 4 3 5.0)<10 (b) Fig. 2 (a) Values of 13(I) from singular-point analysis; (b) values of M(!) — 13(1). 3228 / OPTICAL ENGINEERING / October 1994 / Vol. 33 No. 10 FREQUENCY-DOMAIN ADAPTIVE ITERATIVE IMAGE RESTORATION Table 1 SNR improvement and estimated variance of noise for various SNRs with a uniform blurred image. - SNR SNR # of iteration Estimated variance of noise frequency adaptive algorithm: with regularization matrix 11 17.90 20 dB 4.7 dB 30 dB 5.6 dB 14 1.871 iteration adaptive algorithm: with one regularization parameter 20 dB 4.2 dB 30 dB 4.9 dB 18.21 89 75 1.902 a,= 22.70 at 20 dB and 2.270 at 30 dB of the noise (which is estimated by the algorithm) or the degree of smoothness of the solution. The partially restored image is used in extracting the required information. Sufficient conditions for the convergence of the algorithm have been derived in the frequency domain. Spectral-filtering resuits have been used in handling the near-singular-point region and in speeding up the convergence of the iteration. Linear constraints can also be incorporated into the iteration. When nonlinear constraints are used, the linearization step in the analysis of the algorithm cannot be applied, although the algorithm has been shown to converge experimentally. 7 Appendix 100 In this appendix, an analysis is presented of the error between the original image and the regularized estimate, following Refs. 4 and 8. It leads to the determination of an optimal range for the parameter k(l) in Eq. (16). Let us denote by 80 60 x the minimum-norm least-squares solution of Eq. (1), that is, x = (D D Ty Then the regularized solution (A) can be written as Residual 40 (A)=(DTD + ACTCy 1DTy 20 (49) =(DTD+ACTC)'DTDx (50) =P(A)x (51) 0 10 5 Number of iterations Fig. 3 Values of IIY — Dxkll2/N2 with a uniform blurred image for SNR=20 dB and 30 dB. where P (A) = (DTD + ACTC)_ 1DTD. The mean squared error (MSE) is given by (52) E[Ifx — (A)Il2] E[fIe(A)112] =E[(x x)TP(A)TP(A)(x —x)] (53) where E[.] denotes the expectation operator. Since E [(A)I = P(A)x, the first term of Eq. (53) is equal to the variance of i(A), while the second term is equal to the bias of the estimate i. Since D and C are block-circulant matrices, the variance and bias term in Eq. (53) can be described in the discrete Fourier transform (DFT) domain as N2 10 Number of iterations Fig. 4 Values of IIXk—Xk_1112 with a uniform blurred image for SNR=20 dB and 30 dB. Var[(A)] = (Iv(rn) — D(rn)X(rn)12) ID(l)12 xIE \1=:l [ID(i)12+a(i)IC(i)12121 , fN2 1 (&NR 4.7 dB) and 5(d) (NR 5.6 dB). In all cases the restored images are very satisfactory, according to the improvement in SNR and visual inspection. 6 Conclusions In this paper we have proposed a regularized iterative image restoration algorithm according to which a restored image and an estimate of the regularization matrix are provided simultaneously at each iteration step. Since the regularization matrix is chosen to be block circulant, the regularization is adapted to the parameters of each discrete frequency component. No prior knowledge is required about the variance (54) N2 IX(l)I2a(l)2IC(l)4 bias[x(A)}= ' (55) [ID()I2+a()C()I212 where L=: (11,12), i= (m1,m2), with OEl1 N— 1, Ol2E N— 1, Om1 N— 1, Om2N— 1, X(i) and Y(rn) represent the 2-D DFTs of the unstacked image estimate i and the noisy blurred image y, and D(i) and C(i) represent 2-D DFTs of the 2-D sequences that form the block-circulant matrices D TD and C Tc respectively. Taking their derivatives with respect to ci() yields OPTICAL ENGINEERING / October 1994 / Vol. 33 No. 10 I 3229 KANG and KATSAGGELOS (a) (b) (c) (d) Fig. 5 Noisy blurred and restored images for 7x7 uniform blur: (a) SNR=20 dB, noisy, blurred; (b) SNR=20 dB, restored image; (c) SNR=30 dB, noisy, blurred; (d) SNR=30 dB, restored image. 3230 / OPTICAL ENGINEERING / October 1994 / Vol. 33 No. 10 FREQUENCY-DOMAIN ADAPTIVE ITERATIVE IMAGE RESTORATION a Var[i(A)J 8a() = 2(—-- : IY(m) N2rn1 x 3a() ID(PI2IC(l)12 ( Y(n) — D(n)Xk(n)12 D(m)X(m)12) <O N2 a bias[i(A)J Since N2 ak(l)(C(m)X(rn)(2+(1) , , (56) we obtain from Eq. (61) that O<kOPt(l)N2 max[Xk(l)I2fC(l)21 N2 a(l)IX(l)J2C(l)4D(l)2>° L=' [D(i)I2+c(i)C(i)I2]3 =2 1 (57) . m From Eqs. (54) and (56) we conclude that the variance is a strictly positive, monotonically decreasing function of a(l) for a(1) > 0 and that Var(oo)] = 0. Similarly, from Eqs. (55) and (57) we conclude that the bias is a positive, monotonically increasing function ofa(l) for a(l) > 0 and that bias[(O)] = 0 and bias[i(oo)I = 11x112. Thus, the total error as a function of ci(L) is the sum of a monotonically decreasing and a monotonically increasing function ofa(l). Therefore, the MSE has either one minimum or one maximum for 0 <a(l) < oo The derivative of the MSE with respect to a(l) is equal to 8E[IIe(A)112] = 8a() x 2 N2 D(l)121C(l)12 — <0 m =I for 0 <a(l) < - N2 Y(m) — D(m)X(m)2 max1[X(/)2C(l)2] ' N2 m= 1 Y(m) — D(m)X(m)2 >0 for — N2 min1[fX(i)(2C(i)2 I <c(l)<oc . (60) Therefore, according to the conditions (56), (57), (59), and (60), the MSE function has a unique minimum for 0 < c(1) < 00; this minimum is attained by the value aOP(l)' which is in the range : =1 M. E. Zervakis and T. M. Kwon, "Robust estimation techniques in regularized image restoration," Opt. Eng. 29(5), 455-470 (May 1990). 6. 5. J. Reeves and R. M. Mersereau, ' 'Optimal estimation of the regularization parameters and stabilizing functional for regularized image restoration," Opt. Eng. 29, 446-454 (May 1990). 7. M. G. Kang and A. K. Katsaggelos, ' 'Simultaneous iterative image restoration and evaluation of the regularization parameter,' ' IEEE Trans. Signal Processing 40, 2329—2334 (Sep. 1992). (59) 8a() i. H. C. Andrews and B. R. Hunt, Digital Image Restoration, PrenticeHall, New York (1977). 2. A. K. Katsaggelos, Ed., Digital Image Restoration, Springer Series in Information Sciences, Vol. 23, Springer-Verlag, Heidelberg (1991). 3. B. R. Hunt, ' 'Application of constrained least squares estimation to image restoration by digital computers," IEEE Trans. Comput. C-22, 8. A. K. Katsaggelos and M. G. Kang, ' 'Iterative evaluation of the regularization parameter in regularized image restoration,' ' J. Visual N2 8E[Ile(A)I2J References 5. Therefore, 8a() Telescope Science Institute, Baltimore, Maryland. 332—336 (July 1992). . IY()— D(rn)X(m)12 — ) (58) 8E[IIe(A)112] Acknowledgment This work is supported in part by a grant from the Space 805—812 (1973). N2 (aII2cI2 (63) 4. N. P. Galatsanos and A. K. Katsaggelos, "Methods for choosing the regularization parameter and estimating the noise variance in image restoration and their relation,' ' IEEE Trans. Image Processing 1, [lD()I2 a(1)C(1)2I3 N2m = 1 (62) Comm. and Image Rep. 3(6), 446—455 (Dec. 1992). 9. A. K. Katsaggelos, "A general formulation of adaptive iterative image restoration algorithms,' ' in Proc. 1986 International Conf on Acoustics, Speech, Signal Processing, pp. 42—47 (Mar. 1986). 10. A. K. Katsaggelos, ' 'Iterative image restoration algorithm,' ' Opt. Eng. 28(7), 735—748 (July 1989). I 1. A. K. Katsaggelos, J. Biemond, R. W. Schafer, and R. M. Mersereau, ' 'A regularized iterative image restoration algorithm,' ' IEEE Trans. Signal Processing 39, 914—929 (Apr. 1991). 12. A. N. Tikhonov and V. Y. Arsenin, Solution of Ill-posed Problems, V. H. Winston and Sons, Washington, DC (1977). 13. F. C. Schweppe, Uncertain Dynamic Systems, Prentice-Hall, Engle- woodCliffs,NJ (1973). 14. 0. N. Strand, ' 'Theory and methods related to the singular function expansion and Landweber's iteration for integral equation of the first kind," SIAM J. Numer. Anal. 11, 798-825 (Sep. 1974). 15. R. W. Schafer, R. M. Mersereau, and M. A. Richards, ' 'Constrained iterative restoration algorithms," Proc. IEEE 69, 432—450 (Apr. 1981). 16. J. M. Ortega and W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in Several Variables, Academic Press, New York (1970). 17. A. Ben-Israel and T. N. E. Greville, Generalized Inverses: Theory and Applications, Wiley, New York (1974). 18. B. J. Sullivan and A. K. Katsaggelos, "New termination rule for linear iterative image restoration algorithm," Opt. Eng. 29, 471—477 (May 1990). N2 Y(rn) - D(m)X(m)2 — N2 max1[IX(l)I2C(l)2] Moon Gi Kang received BS and MS de- grees in electronics engineering from <a() Seoul National University, Korea, in 1986 and 1988, respectively. He received a PhD degree in electrical engineering from N2 fY(rn) — D(m)X(m)2 — <=1 N2 min1[IX(i)121C(i)12] (61) Northwestern University in 1994. He was a research assistant from 1989 to 1993, and since 1994 he has been working as a research fellow in the Department of Electrical Engineering and Computer Science at Northwestern University. His current reOPTICAL ENGINEERING I October 1994 I Vol. 33 No. 10 I 3231 KANG and KATSAGGELOS search interests include regularized iterative image restoration and use of higher order spectra for image restoration. Aggelos K. Katsaggelos received the Diploma degree in electrical and mechanical engineering from the Aristotelian University of Thessaloniki, Thessaloniki, Greece, in 1979, and the MS and PhD degrees, both in electrical engineering, from the Georgia Institute of Technology, Atlanta, Georgia, in 1981 and 1985, respectively. In 1985 he joined the Department of Electrical Engineering and Computer Science at Northwestern University, Evanston, Illi- nois, where he is currently an associate professor. During the 1986—1987 academic year he was an assistant professor at Polytechnic University, Department of Electrical Engineering and Com- 3232 / OPTICAL ENGINEERING / October 1994 / Vol. 33 No. 10 puter Science, Brooklyn, New York. His current research interests include image recovery, processing of moving images (motion estimation, enhancement, very low bit rate compression), and computational vision. Dr. Katsaggelos is an Ameritech Fellow and a member of the Associate Staff, Department of Medicine, at Evanston Hospital. He is a senior member of the IEEE, and also a member of the SPIE, the Steering Committees of the IEEE Transactions on Medical Imaging and the IEEE Transactions on Image Processing, the IEEE Technical Committees on Visual Signal Processing and Communications and on Image and Multi-Dimensional Signal Pro- cessing, the Technical Chamber of Commerce of Greece, and Sigma Xi. He has served as an associate editor for the IEEE Transactions on Signal Processing (1990—1992), and is currently an area editor for the journal Graphical Models and Image Processing. He is the editor of Digital Image Restoration (Springer-Verlag, Heidelberg, 1991) and the general chairman of the 1994 Visual Communications and Image Processing Conference (Chicago).
© Copyright 2024