Wavelet Domain Image Restoration Using Edge Preserving Prior Models
Wavelet Domain Image Restoration Using Edge Preserving Prior Models Murat Belge and Eric L. Miller Department of Electrical and Computer Engineering Northeastern University Boston, MA 02115
Abstract
In this paper we consider a wavelet-based edge preserving regularization scheme for use in linear image restoration problems. Our efforts build on a collection of mathematical results indicating that wavelets are especially useful for representing functions that contain discontinuities (i.e. edges in two dimensions or jumps in 1D). We interpret the resulting theory in a statistical signal processing framework and obtain a highly exible framework for adapting the degree of regularization to the scale based structure of the underlying image. We demonstrate an ecient algorithm for obtaining the reconstructions from observed data and for choosing the multiple regularization parameters governing the priors.
1 Introduction
In this work, we consider the problem of recovering an unknown image x from measurements y where the two are related via y = Hx + n (1)
with H a known matrix de ning the linear degradation, x and y vectors representing lexiographically ordered pixels in the original and degraded images and n representing additive white Gaussian noise (i.e. n N(0; 2I )). The above problem is typically ill-posed in the sense that the operator H may not be invertible or it may be ill-conditioned. In such a case, the original problem is replaced by a well conditioned one whose solution approximates that of the original. This is usually achieved by incorporating some a priori information about the original image. In the stochastic interpretation of image restoration, the prior information on the image x can be expressed in the form of an a priori probability distribution P (xj) which is assumed to be dependent on a set of parameters . Bayes rule allows us to combine This work was supported by NSF CAREER Grant MIP9623721 and a subcontract GC123920NDG from Boston University under the AFOSR MURI Program on Reduced Signature Target Recognition.
the prior information with the information contained in the data to obtain the a posteriori distribution
P (xjy; ; )/ exp ? 21 2 ky ? Hxk22 ? (x; )
(2)
where the prior distribution on x is assumed to be proportional to expf?(x; )g. Conventional prior distributions are typically expressed as, expf? 212 kLxk22g, which is equivalent to the assumption that x is coming from a Gaussian distribution with covariance matrix (LT L)?1. The restored image x is then found by maximizing the log posterior probability with respect to x. i.e 1 2 x (; ) = arg min x 22 ky ? Hxk2 + (x; ): (3) The quality of the restoration depends on both the form of the prior (i.e. form of (:) in (3)) and the choice of hyperparameters ( ; ). In this paper, we consider a Bayesian image restoration algorithm in which the desired image's subbands are modeled as stationary generalized Gaussian processes. The intuition underlying our modeling is that images typically have spatial structure consisting of smooth areas dispersed with occasional edges. When the image is subjected to a subband decomposition, band pass and high pass pixels of the image are seen to have signi cantly non-Gaussian behavior with a distribution function sharply peaked about zero (contributions coming from smooth areas) and heavy tails which are mainly due to the edges [5]. Gaussian modeling puts too much penalty on the sharp transitions (edges) in the image therefore the restorations based on Gaussian modeling look smooth. Several authors used generalized Gaussian distributions to model subband statistics in the context of image coding and denoising. However, to our knowledge such models have not been used in connection with image restoration yet. In this paper, we develop new multiresolution edge preserving image models for image restoration applications and solve the nonlinear optimization and parameter selection issues associated with the new modeling.
2 Wavelet Transform of Images
The wavelet transform is a time-frequency representation which has good localization in both domains. We restrict our attention to orthonormal separable wavelet transforms and assume that the reader is familiar with the wavelet theory [1]. The following discussion is intended to establish the notation that will be used in the rest of the paper. Following the wavelet literature, elements of the image x(m; n), 1 m; n 2J , are called the nest scale scaling coecients. From these nest scale scaling coecients we obtain the scaling coecients at a lower resolution level which represent a coarser scale version of x and a set of wavelet coecients which represents the horizontal, vertical and diagonal details lost in moving from a ner scale to a coarser scale. By repeating this procedure on the coarse scale coecients k ? 1 times, we can obtain a k-level wavelet decomposition. We denote the scaling coecients at scale j0 and position (m; n) by xj0 ;0(m; n), 1 m; n 2j0 and the wavelet coecients at scale j , j0 j J ? 1, orientation i, 1 i 3 and position (m; n) by xj0;i (m; n). The wavelet transform can be represented as a orthonormal matrix W , which transforms an image x into its wavelet domain representation through matrix multiplication. That is, x^ = W x. By using this, (1) can be transformed into wavelet domain as ?
W HW T W x + W n
Wy = y^ = H^ x^ + n^ ; (4) where W T W = I follows from the orthonormality of the wavelet transformation.
3 A Multiscale Image Model for Use in Image Restoration
Following the work in [5] on image coding and denoising, we model the distribution of wavelet coecients of images by a Generalized-Gaussian (GG) density [9]. The following distribution function serves as the prior model xj;i(m; n) p 1 P (xj;i (m; n)jp; j;i(m; n)) / exp ? p (m; n) ; j;i (5) where 1 p 2 is a parameter which determines the tail behavior of the density function and j;i (m; n) is the scale parameter similar to the variance of a Gaussian density. For simplicity, we will refer to the density in (5) as GG(0; j;i(m; n); p). For p = 1 we have the Laplacian or double exponential density and for p = 2 we have the familiar Gaussian density. Intermediate values of p correspond to increasingly heavier tails. We assume that the scaling coecients xj0;0 , 1 m; n 2J ?L , are i.i.d. with GG(0; j0;0 (m; n); p) density. The speci cation of one parameter for every wavelet coecient results in an image model far too complex to be of use in a restoration procedure. The structure of the model in (5) coupled with the speci cation of the problem in the wavelet domain does
suggest a variety of simpli cations which are of use for the restoration problem. In this work, we consider the following two: 1. Model 1: The variance of the wavelet coef cients decreases exponentially with the scale, i.e, xj;i(m; n) = GG(0; 2? (j ?j0) ; p); i = 1; 2; 3; 1 m; n 2j with j0 the coarsest scale, the scale parameter corresponding to j0 and 0. The rationale behind this model is that it is equivalent to a deterministic modeling of the image as a member of Besov spaces [1]. 2. Model 2: Scaling coecients are i.i.d. with GG(0; j0;0; p) distribution and wavelet coecients at a particular scale are i.i.d. with GG(0; j ; p); j = J ? L; : : :; J ? 1. This model is useful in cases where the variance of the wavelet coecients at di erent scales can not be wellapproximated by a simple exponential law. For the remainder of this paper, we assume that the two parameters and p in the above models are known and xed a priori. Generally, the performance of the regularizer is impacted to a greater extent by the online identi cation of the parameters and it is here where we choose to concentrate our e orts.
4
A Multiscale Bayesian Image Restoration Algorithm
Motivated by the above discussion, we consider a restoration scheme for obtaining an estimate of the image's wavelet coecients x^ which is de ned in terms of the following optimization problem: x^ = arg min fky^ ? H^ x^ k22 + j0 ;0kx^ j0 ;0 kpp ^ +
x
XX
j
i
j;i kx^ j;ikpp g;
(6)
where the notation x^j;i denotes the vector of wavelet(scaling) coecients at scale j and orientation i 2 2 and j;i = p . The cost function in (6) has a unique minimum for a given . This follows from the fact that it is just the summation of an L2 and an L1 norm which are both convex in x^ . The rst order conditions that must be satis ed by x^ is found by di erentiating (6) with respect to x^ . However, the cost function in (6) is not di erentiable in its current form due to the fact that the lp norm term present in the cost function is no di erentiabl at zero. To alleviate this problem, we slightly perturb the lp norm term to ?obtain a difP ferentiable approximation: kxkpp i x2i + p=2 ; where xi is the ith element of the vector x and > 0 is a small constant. Substituting this approximation into (6) and di erentiating with respect to x^ we obtain the following equation # " j;i (7) D = diag (jx^ (m; n)j2 + )1?p=2 j;i H^ T H^ + 2p D x^ = H^ T y^ (8) j;i
optimal(i.e. minimizing the restoration error). The necessary formulas for the computation of Gaussian curvature of the L-hypersurface for the problem in (6) are provided in [8].
4
x 10 8
26.5
7
26
5 Simulation Study
6 25.5 RMSE
Curvature
5 4 3
25 24.5
2 24
1 0 −1
23.5 −1 −2
−1 −3
−1.5 −4 −5
log10 λ
1
−2
−1 −3
−1.5 −4
−2
−5
−2.5 −6
−3
−2
log10 λ2
log10 λ
1
−2.5 −6
−3
log10 λ
2
Figure 1: Curvature of the L-hypersurface and the corresponding RMSE plot. An iterative algorithm based on (8) can be developed to approximate the solution x^ . Starting from an initial extimate x^ (0) , we solve the following equation for x^(k+1) until convergence is achieved (9) H^ T H^ + 2p D(k) x^(k+1) = H^ T y^; where D(k) denotes the diagonal matrix obtained by replacing x^ in (7) by x^(k). It can be shown that the iterative algorithm in (9) converges to the unique minimum of the perturped cost function in (6). The stabilization constant a ects the convergence rate. The bigger the , the faster the convergence. However, using a relatively large has the side e ect of smoothing out the edges in the image. Therefore, should be set to balance the compromise between the convergence rate of the algorithm and sharpness of the resulting solution. Based on our experience, we recommend 1 for problems involving real life images.
4.1
Determination of Regularization Prameters
A noteworthy feature of the mutiscale image restoration algorithm introduced in the previous section is the need for multiple regularization parameters. In this paper, we utilize a multidimensional extension of the popular L-curve method [7] called the \L-hypersurface" to determine mutiple regularization parameters appearing in (6). The L-hypersurface is a plot of the log of residual norm, ky? Hx (p)k22, against the log of side constraint norms kx^ j;i ()kp , j0 j J ? 1, for a range of regularization parameters. It has been argued and shown through numerical experiments that the points on the L-hypersurface where the Gaussian curvature reaches a local maxima are closely tied to the points on the restoration error, kx ? x ()k22 , hypersurface where the restoration error reaches a local minima. The Lhypersurface selection for corresponds to that for which the Gaussian curvature of the L-hypersurface is a maximum. Monte Carlo simulations indicates that [8] the L-hypersurface method produces regularization parameters which are almost as good as the
In this section, we illustrate the performance of our proposed multiscale image restoration algorithm for both real and synthetic images. In the rst example, we used a 2-D Gaussian convolutional kernel, with x = y = 2:0 to blur 256 256 Mandrill image. Zero mean white Gaussian noise was added to set the SNR to 30dB. The top two plots in Figure 2 display the original and the blurred, noisy images. We restored the degraded Mandrill image using three regularization techniques: our proposed multiscale regularization scheme, the Constrained Least Squares (CLS) algorithm with a 2-D Laplacian regularizer [3], and the Total Variation (TV) algorithm [6]. The TV algorithm has been recently developed in the context of edge preserving regularization. In this approach the regularization function is the TV of the image which is simply the l1 norm of the gradient of the image. The TV functional does not penalize the discontinuities as harshly as a quadratic functional and restorations obtained by this method are superior to the conventional techniques especially if the exact image is piecewise continuous. The relevant regularization parameters were determined using either the L-curve or or the L-hypersurface method. For the TV algorithm and our algorithm we used = 1:0 as the stabilization constant. For our multiscale image restoration algorithm, we used the Daubechie's eight tap most symmetrical wavelets [2]. In Figure 2-3 we display the restored Mandrill images. For our multiscale image restoration method we computed two restorations according to the regularization schemes Model 1 and Model 2 described in section 3 with p = 1:0. For Model 1, we used a 5-level wavelet decomposition and set the exponential parameter to = 1:2. For Model 2, we used a 3level wavelet decomposition and set the regularization for the scaling coecients to 10?5 . For the Model 1 case the L-hypersurface was used to determine two parameters corresponding to the scaling coecients and the coarsest scale wavelet coecients. In this case, L-hypersurface is a 2-D function of regularization parameters as seen in Figure 1. Also shown in Fig. q 1 is a plot of the root mean square error (RMSE), 1 ^ k22, as a function of these regularization N 2 kx ? x parameters. Examining these images shows that the L-surface has a distinct extended maxima along which the RMSE is very close to being a minimum. Thus, we see that the restoration algorithm is not overly sensitive to the scaling coecient regularization parameter and locating the correct regularization parameter for the wavelet coecients is more important. In Model 2 restoration, we employed a 3-level wavelet decomposition and assumed that each scale gets a di erent regularization parameter. Based on the insensitivity of the restoration to the scaling coecient regularization parameter we set this value to
10?5. Fig. 3 shows that both the TV algorithm and our algorithm produce restored images visually superior to the CLS algorithm. We also observe that the images restored by our algorithm are a little sharper than that of the TV and the texture-like regions abundant in the Mandrill image (eg. the hairs around the mouth of the Mandrill) are better preserved by our algorithm. The RMSE values are 23.73 for the TV algorithm, 24.11 for the CLS algorithm and 23.75 and 23.84 for the Model 1 and Model 2 restorations respectively. Finally, for this image we see little di erence either in terms of the error norm or in terms of visual quality between the the restorations using our method. The primary advantage of the Model 2 approach is that we no longer need to independently set the parameter. In our nal example, we rst blurred the original Bridge image with a 9 9 uniform motion blur and added white Gaussian noise to the degraded image to set the SNR at 40dB. Figure 4 shows the original and the degraded images. We display the restorations obtained by the TV and the proposed algorithm in Figure 5. For our restoration, we applied the Model 1 regularization scheme with = 1:2, p = 1:2 and a 5-level wavelet decomposition. We used the appropriated 2-D L-hypersurface to determine the relevant regularization parameters. Although, in this case the RMSE values were similar (19.43 for TV and 19.72 for the multiscale algorithm), the two restored images in Figure 5 exhibit vastly di erent visual characteristics. The TV algorithm completely wipes out the small features in the image and produces an overly homogenized restoration resembling an \oil painting" of the original scene. On the other hand, the wavelet-based algorithm is able to reproduce ner detail thereby yielding a more visually appealing restoration.
6 Summary and Conclusions
In this paper, we introduced a wavelet domain multiscale image restoration algorithm for use in linear image restoration problems. Following the recent results in the the area of image denoising and coding, we developed a statistical prior model for the wavelet coecients of images. Our priors are able to capture spatial, scale space characteristics of images accurately. We developed an ecient iterative optimization algorithm to solve the nonlinear optimization problem resulting from using such priors and introduced the L-hypersurface method to choose the multiple hyperparameters governing the structure of our priors. Comparison of our multiscale image restoration algorithm with either the traditional image restoration algorithms or the more recent edge-preserving image restoration algorithms shows that our algorithm can produce restorations which are visually signi cantly better than that of the traditional techniques and at least comparable, if not better, than that of the the edge-preserving algorithms.
References
[1] Y. Meyer, Wavelets and Operators, Cambridge Univ. Press, 1992.
[2] I. Daubechies, Ten Lectures on Wavelets. NewYork: SIAM Press, 1992. [3] H. C. Andrews and B. R. Hunt, Digital Image Restoration, Englewood Cli s, NJ:Prentice-Hall, 1977. [4] P. Charbonnier, L. Blanc-Feraud, G. Aubert and M. Barlaud, \Deterministic edge-preserving regularization in computed imaging," IEEE Trans. Image Processing, vol. 6, no. 2, pp. 298-311, February 1997. [5] E. P. Simoncelli, E. Adelson, \Noise removal via Bayesian wavelet coring," in Proceedings 3rd International Conference on Image Processing, Lausanne, Switzerland, September 1996. [6] L. I. Rudin, S. Osher, E. Fatemi, \Nonlinear total variation based noise removal algorithms," Physica D, vol. 60, pp. 259-268, 1992. [7] P. C. Hansen, \Analysis of discrete ill-posed problems by means of the L-curve," SIAM Review, vol. 34, pp. 561-580, 1992. [8] M. Belge, M. E. Kilmer and E. L. Miller, \Simultaneous multiple regularization parameter selection by means of teh L-hypersurface with applications to linear inverse problems posed in the wavelet domain," to be published in Proceedings of SPIE, San Diego, U.S.A, July 1998. [9] M. Belge, M. E. Kilmer and E. L. Miller, \Wavelet Domain Image Restoration With Adaptive EdgePreserving Regularization," submitted to IEEE Trans. Image Processing.
Figure 2: Top to bottom: Original, blurred, restored by the CLS algorithm.
Figure 3: Top to bottom: restored by TV, restored by Model 1 scheme and restored by Model 2 scheme.
Figure 4: Top: Original Bridge image. Bottom: Blurred, noisy Bridge image.
Figure 5: Top: restored by the TV. Bottom: restored by Model 1 scheme.
Abstract
In this paper we consider a wavelet-based edge preserving regularization scheme for use in linear image restoration problems. Our efforts build on a collection of mathematical results indicating that wavelets are especially useful for representing functions that contain discontinuities (i.e. edges in two dimensions or jumps in 1D). We interpret the resulting theory in a statistical signal processing framework and obtain a highly exible framework for adapting the degree of regularization to the scale based structure of the underlying image. We demonstrate an ecient algorithm for obtaining the reconstructions from observed data and for choosing the multiple regularization parameters governing the priors.
1 Introduction
In this work, we consider the problem of recovering an unknown image x from measurements y where the two are related via y = Hx + n (1)
with H a known matrix de ning the linear degradation, x and y vectors representing lexiographically ordered pixels in the original and degraded images and n representing additive white Gaussian noise (i.e. n N(0; 2I )). The above problem is typically ill-posed in the sense that the operator H may not be invertible or it may be ill-conditioned. In such a case, the original problem is replaced by a well conditioned one whose solution approximates that of the original. This is usually achieved by incorporating some a priori information about the original image. In the stochastic interpretation of image restoration, the prior information on the image x can be expressed in the form of an a priori probability distribution P (xj) which is assumed to be dependent on a set of parameters . Bayes rule allows us to combine This work was supported by NSF CAREER Grant MIP9623721 and a subcontract GC123920NDG from Boston University under the AFOSR MURI Program on Reduced Signature Target Recognition.
the prior information with the information contained in the data to obtain the a posteriori distribution
P (xjy; ; )/ exp ? 21 2 ky ? Hxk22 ? (x; )
(2)
where the prior distribution on x is assumed to be proportional to expf?(x; )g. Conventional prior distributions are typically expressed as, expf? 212 kLxk22g, which is equivalent to the assumption that x is coming from a Gaussian distribution with covariance matrix (LT L)?1. The restored image x is then found by maximizing the log posterior probability with respect to x. i.e 1 2 x (; ) = arg min x 22 ky ? Hxk2 + (x; ): (3) The quality of the restoration depends on both the form of the prior (i.e. form of (:) in (3)) and the choice of hyperparameters ( ; ). In this paper, we consider a Bayesian image restoration algorithm in which the desired image's subbands are modeled as stationary generalized Gaussian processes. The intuition underlying our modeling is that images typically have spatial structure consisting of smooth areas dispersed with occasional edges. When the image is subjected to a subband decomposition, band pass and high pass pixels of the image are seen to have signi cantly non-Gaussian behavior with a distribution function sharply peaked about zero (contributions coming from smooth areas) and heavy tails which are mainly due to the edges [5]. Gaussian modeling puts too much penalty on the sharp transitions (edges) in the image therefore the restorations based on Gaussian modeling look smooth. Several authors used generalized Gaussian distributions to model subband statistics in the context of image coding and denoising. However, to our knowledge such models have not been used in connection with image restoration yet. In this paper, we develop new multiresolution edge preserving image models for image restoration applications and solve the nonlinear optimization and parameter selection issues associated with the new modeling.
2 Wavelet Transform of Images
The wavelet transform is a time-frequency representation which has good localization in both domains. We restrict our attention to orthonormal separable wavelet transforms and assume that the reader is familiar with the wavelet theory [1]. The following discussion is intended to establish the notation that will be used in the rest of the paper. Following the wavelet literature, elements of the image x(m; n), 1 m; n 2J , are called the nest scale scaling coecients. From these nest scale scaling coecients we obtain the scaling coecients at a lower resolution level which represent a coarser scale version of x and a set of wavelet coecients which represents the horizontal, vertical and diagonal details lost in moving from a ner scale to a coarser scale. By repeating this procedure on the coarse scale coecients k ? 1 times, we can obtain a k-level wavelet decomposition. We denote the scaling coecients at scale j0 and position (m; n) by xj0 ;0(m; n), 1 m; n 2j0 and the wavelet coecients at scale j , j0 j J ? 1, orientation i, 1 i 3 and position (m; n) by xj0;i (m; n). The wavelet transform can be represented as a orthonormal matrix W , which transforms an image x into its wavelet domain representation through matrix multiplication. That is, x^ = W x. By using this, (1) can be transformed into wavelet domain as ?
W HW T W x + W n
Wy = y^ = H^ x^ + n^ ; (4) where W T W = I follows from the orthonormality of the wavelet transformation.
3 A Multiscale Image Model for Use in Image Restoration
Following the work in [5] on image coding and denoising, we model the distribution of wavelet coecients of images by a Generalized-Gaussian (GG) density [9]. The following distribution function serves as the prior model xj;i(m; n) p 1 P (xj;i (m; n)jp; j;i(m; n)) / exp ? p (m; n) ; j;i (5) where 1 p 2 is a parameter which determines the tail behavior of the density function and j;i (m; n) is the scale parameter similar to the variance of a Gaussian density. For simplicity, we will refer to the density in (5) as GG(0; j;i(m; n); p). For p = 1 we have the Laplacian or double exponential density and for p = 2 we have the familiar Gaussian density. Intermediate values of p correspond to increasingly heavier tails. We assume that the scaling coecients xj0;0 , 1 m; n 2J ?L , are i.i.d. with GG(0; j0;0 (m; n); p) density. The speci cation of one parameter for every wavelet coecient results in an image model far too complex to be of use in a restoration procedure. The structure of the model in (5) coupled with the speci cation of the problem in the wavelet domain does
suggest a variety of simpli cations which are of use for the restoration problem. In this work, we consider the following two: 1. Model 1: The variance of the wavelet coef cients decreases exponentially with the scale, i.e, xj;i(m; n) = GG(0; 2? (j ?j0) ; p); i = 1; 2; 3; 1 m; n 2j with j0 the coarsest scale, the scale parameter corresponding to j0 and 0. The rationale behind this model is that it is equivalent to a deterministic modeling of the image as a member of Besov spaces [1]. 2. Model 2: Scaling coecients are i.i.d. with GG(0; j0;0; p) distribution and wavelet coecients at a particular scale are i.i.d. with GG(0; j ; p); j = J ? L; : : :; J ? 1. This model is useful in cases where the variance of the wavelet coecients at di erent scales can not be wellapproximated by a simple exponential law. For the remainder of this paper, we assume that the two parameters and p in the above models are known and xed a priori. Generally, the performance of the regularizer is impacted to a greater extent by the online identi cation of the parameters and it is here where we choose to concentrate our e orts.
4
A Multiscale Bayesian Image Restoration Algorithm
Motivated by the above discussion, we consider a restoration scheme for obtaining an estimate of the image's wavelet coecients x^ which is de ned in terms of the following optimization problem: x^ = arg min fky^ ? H^ x^ k22 + j0 ;0kx^ j0 ;0 kpp ^ +
x
XX
j
i
j;i kx^ j;ikpp g;
(6)
where the notation x^j;i denotes the vector of wavelet(scaling) coecients at scale j and orientation i 2 2 and j;i = p . The cost function in (6) has a unique minimum for a given . This follows from the fact that it is just the summation of an L2 and an L1 norm which are both convex in x^ . The rst order conditions that must be satis ed by x^ is found by di erentiating (6) with respect to x^ . However, the cost function in (6) is not di erentiable in its current form due to the fact that the lp norm term present in the cost function is no di erentiabl at zero. To alleviate this problem, we slightly perturb the lp norm term to ?obtain a difP ferentiable approximation: kxkpp i x2i + p=2 ; where xi is the ith element of the vector x and > 0 is a small constant. Substituting this approximation into (6) and di erentiating with respect to x^ we obtain the following equation # " j;i (7) D = diag (jx^ (m; n)j2 + )1?p=2 j;i H^ T H^ + 2p D x^ = H^ T y^ (8) j;i
optimal(i.e. minimizing the restoration error). The necessary formulas for the computation of Gaussian curvature of the L-hypersurface for the problem in (6) are provided in [8].
4
x 10 8
26.5
7
26
5 Simulation Study
6 25.5 RMSE
Curvature
5 4 3
25 24.5
2 24
1 0 −1
23.5 −1 −2
−1 −3
−1.5 −4 −5
log10 λ
1
−2
−1 −3
−1.5 −4
−2
−5
−2.5 −6
−3
−2
log10 λ2
log10 λ
1
−2.5 −6
−3
log10 λ
2
Figure 1: Curvature of the L-hypersurface and the corresponding RMSE plot. An iterative algorithm based on (8) can be developed to approximate the solution x^ . Starting from an initial extimate x^ (0) , we solve the following equation for x^(k+1) until convergence is achieved (9) H^ T H^ + 2p D(k) x^(k+1) = H^ T y^; where D(k) denotes the diagonal matrix obtained by replacing x^ in (7) by x^(k). It can be shown that the iterative algorithm in (9) converges to the unique minimum of the perturped cost function in (6). The stabilization constant a ects the convergence rate. The bigger the , the faster the convergence. However, using a relatively large has the side e ect of smoothing out the edges in the image. Therefore, should be set to balance the compromise between the convergence rate of the algorithm and sharpness of the resulting solution. Based on our experience, we recommend 1 for problems involving real life images.
4.1
Determination of Regularization Prameters
A noteworthy feature of the mutiscale image restoration algorithm introduced in the previous section is the need for multiple regularization parameters. In this paper, we utilize a multidimensional extension of the popular L-curve method [7] called the \L-hypersurface" to determine mutiple regularization parameters appearing in (6). The L-hypersurface is a plot of the log of residual norm, ky? Hx (p)k22, against the log of side constraint norms kx^ j;i ()kp , j0 j J ? 1, for a range of regularization parameters. It has been argued and shown through numerical experiments that the points on the L-hypersurface where the Gaussian curvature reaches a local maxima are closely tied to the points on the restoration error, kx ? x ()k22 , hypersurface where the restoration error reaches a local minima. The Lhypersurface selection for corresponds to that for which the Gaussian curvature of the L-hypersurface is a maximum. Monte Carlo simulations indicates that [8] the L-hypersurface method produces regularization parameters which are almost as good as the
In this section, we illustrate the performance of our proposed multiscale image restoration algorithm for both real and synthetic images. In the rst example, we used a 2-D Gaussian convolutional kernel, with x = y = 2:0 to blur 256 256 Mandrill image. Zero mean white Gaussian noise was added to set the SNR to 30dB. The top two plots in Figure 2 display the original and the blurred, noisy images. We restored the degraded Mandrill image using three regularization techniques: our proposed multiscale regularization scheme, the Constrained Least Squares (CLS) algorithm with a 2-D Laplacian regularizer [3], and the Total Variation (TV) algorithm [6]. The TV algorithm has been recently developed in the context of edge preserving regularization. In this approach the regularization function is the TV of the image which is simply the l1 norm of the gradient of the image. The TV functional does not penalize the discontinuities as harshly as a quadratic functional and restorations obtained by this method are superior to the conventional techniques especially if the exact image is piecewise continuous. The relevant regularization parameters were determined using either the L-curve or or the L-hypersurface method. For the TV algorithm and our algorithm we used = 1:0 as the stabilization constant. For our multiscale image restoration algorithm, we used the Daubechie's eight tap most symmetrical wavelets [2]. In Figure 2-3 we display the restored Mandrill images. For our multiscale image restoration method we computed two restorations according to the regularization schemes Model 1 and Model 2 described in section 3 with p = 1:0. For Model 1, we used a 5-level wavelet decomposition and set the exponential parameter to = 1:2. For Model 2, we used a 3level wavelet decomposition and set the regularization for the scaling coecients to 10?5 . For the Model 1 case the L-hypersurface was used to determine two parameters corresponding to the scaling coecients and the coarsest scale wavelet coecients. In this case, L-hypersurface is a 2-D function of regularization parameters as seen in Figure 1. Also shown in Fig. q 1 is a plot of the root mean square error (RMSE), 1 ^ k22, as a function of these regularization N 2 kx ? x parameters. Examining these images shows that the L-surface has a distinct extended maxima along which the RMSE is very close to being a minimum. Thus, we see that the restoration algorithm is not overly sensitive to the scaling coecient regularization parameter and locating the correct regularization parameter for the wavelet coecients is more important. In Model 2 restoration, we employed a 3-level wavelet decomposition and assumed that each scale gets a di erent regularization parameter. Based on the insensitivity of the restoration to the scaling coecient regularization parameter we set this value to
10?5. Fig. 3 shows that both the TV algorithm and our algorithm produce restored images visually superior to the CLS algorithm. We also observe that the images restored by our algorithm are a little sharper than that of the TV and the texture-like regions abundant in the Mandrill image (eg. the hairs around the mouth of the Mandrill) are better preserved by our algorithm. The RMSE values are 23.73 for the TV algorithm, 24.11 for the CLS algorithm and 23.75 and 23.84 for the Model 1 and Model 2 restorations respectively. Finally, for this image we see little di erence either in terms of the error norm or in terms of visual quality between the the restorations using our method. The primary advantage of the Model 2 approach is that we no longer need to independently set the parameter. In our nal example, we rst blurred the original Bridge image with a 9 9 uniform motion blur and added white Gaussian noise to the degraded image to set the SNR at 40dB. Figure 4 shows the original and the degraded images. We display the restorations obtained by the TV and the proposed algorithm in Figure 5. For our restoration, we applied the Model 1 regularization scheme with = 1:2, p = 1:2 and a 5-level wavelet decomposition. We used the appropriated 2-D L-hypersurface to determine the relevant regularization parameters. Although, in this case the RMSE values were similar (19.43 for TV and 19.72 for the multiscale algorithm), the two restored images in Figure 5 exhibit vastly di erent visual characteristics. The TV algorithm completely wipes out the small features in the image and produces an overly homogenized restoration resembling an \oil painting" of the original scene. On the other hand, the wavelet-based algorithm is able to reproduce ner detail thereby yielding a more visually appealing restoration.
6 Summary and Conclusions
In this paper, we introduced a wavelet domain multiscale image restoration algorithm for use in linear image restoration problems. Following the recent results in the the area of image denoising and coding, we developed a statistical prior model for the wavelet coecients of images. Our priors are able to capture spatial, scale space characteristics of images accurately. We developed an ecient iterative optimization algorithm to solve the nonlinear optimization problem resulting from using such priors and introduced the L-hypersurface method to choose the multiple hyperparameters governing the structure of our priors. Comparison of our multiscale image restoration algorithm with either the traditional image restoration algorithms or the more recent edge-preserving image restoration algorithms shows that our algorithm can produce restorations which are visually signi cantly better than that of the traditional techniques and at least comparable, if not better, than that of the the edge-preserving algorithms.
References
[1] Y. Meyer, Wavelets and Operators, Cambridge Univ. Press, 1992.
[2] I. Daubechies, Ten Lectures on Wavelets. NewYork: SIAM Press, 1992. [3] H. C. Andrews and B. R. Hunt, Digital Image Restoration, Englewood Cli s, NJ:Prentice-Hall, 1977. [4] P. Charbonnier, L. Blanc-Feraud, G. Aubert and M. Barlaud, \Deterministic edge-preserving regularization in computed imaging," IEEE Trans. Image Processing, vol. 6, no. 2, pp. 298-311, February 1997. [5] E. P. Simoncelli, E. Adelson, \Noise removal via Bayesian wavelet coring," in Proceedings 3rd International Conference on Image Processing, Lausanne, Switzerland, September 1996. [6] L. I. Rudin, S. Osher, E. Fatemi, \Nonlinear total variation based noise removal algorithms," Physica D, vol. 60, pp. 259-268, 1992. [7] P. C. Hansen, \Analysis of discrete ill-posed problems by means of the L-curve," SIAM Review, vol. 34, pp. 561-580, 1992. [8] M. Belge, M. E. Kilmer and E. L. Miller, \Simultaneous multiple regularization parameter selection by means of teh L-hypersurface with applications to linear inverse problems posed in the wavelet domain," to be published in Proceedings of SPIE, San Diego, U.S.A, July 1998. [9] M. Belge, M. E. Kilmer and E. L. Miller, \Wavelet Domain Image Restoration With Adaptive EdgePreserving Regularization," submitted to IEEE Trans. Image Processing.
Figure 2: Top to bottom: Original, blurred, restored by the CLS algorithm.
Figure 3: Top to bottom: restored by TV, restored by Model 1 scheme and restored by Model 2 scheme.
Figure 4: Top: Original Bridge image. Bottom: Blurred, noisy Bridge image.
Figure 5: Top: restored by the TV. Bottom: restored by Model 1 scheme.
