Multicomponent Image Restoration, an Experimental Study

Multicomponent Image Restoration, an Experimental Study Arno Duijster, Steve De Backer, and Paul Scheunders IBBT, Vision Lab, University of Antwerp, Universiteitsplein 1, 2610 Wilrijk, Belgium {arno.duijster,steve.debacker,paul.scheunders}@ua.ac.be http://visielab.ua.ac.be

Abstract. In this paper, we study the problem of restoring multicomponent images. In particular, we investigate the effects of accounting for the correlation between the image components on the deconvolution and denoising steps. The proposed restoration is a 2-step procedure, comprising a shrinkage in the Fourier domain, followed by a shrinkage in the wavelet domain. The Fourier shrinkage is performed in a decorrelated space, by performing PCA before the Fourier transform. The wavelet shrinkage is performed in a Bayesian denoising framework by applying multicomponent probability density models for the wavelet coefficients that fully account for the intercomponent correlations. In an experimental section, we compare this procedure with the single-component analogies, i.e. performing the Fourier shrinkage in the correlated space and using single-component probability density models for the wavelet coefficients. In this way, the effect of the multicomponent procedures on the deconvolution and denoising performance is studied experimentally. Key words: multicomponent image restoration, deconvolution, denoising, wavelet transform

1

Introduction

With the evolution of imaging technology, an increasing number of imaging modalities becomes available. In this paper, we focus on what is called multichannel or multicomponent imagery, which are images containing several image components. The most well-known example is a color image, which generally consists of a red, green and blue component. In remote sensing, sensors are available that can generate multispectral or hyperspectral data, involving a few to more than hundred bands. In medical imagery, distinct image modalities reveal different features of the internal body. Examples are MRI images acquired by using different imaging parameters (T1, T2, proton density, diffusion, . . . ), different CT and nuclear medicine imaging modalities. Remote sensed and medical images are often degraded due to limitations such as aperture effects of the camera, motion or atmospheric effects and physical limitations. As a consequence, the images are corrupted by blur and noise.

2

Arno Duijster, Steve De Backer, and Paul Scheunders

The standard imaging model consists of a convolution with some impulse response function and additive noise. The standard way of deconvolution is inverse filtering of the image in the Fourier domain. A regularization is required to handle the noise. This leads to the classical Wiener type filters. In the past, specific multicomponent image restoration was proposed [1–3].

A disadvantage of the Fourier transform is that it does not efficiently represent image singularities, so that only small amounts of shrinkage are allowed to avoid distortion of edges in the image. Recently, 2-step procedures were proposed, involving a Fourier shrinkage and a wavelet based image denoising [4–6]. The first step is a shrinkage in the Fourier domain, in order to exploit the Fourier transform’s economical noise representation, followed by a shrinkage in the wavelet domain, to exploit the wavelet transform’s economical representation of piecewise smooth images.

The wavelet transform offers an efficient representation of spatial discontinuities. It compresses the essential information of an image into a relatively few, large coefficients coinciding with the positions of image discontinuities. Such a representation naturally facilitates the construction of spatially adaptive denoising methods that can smooth noise without excessive blurring of image details. Typically, noise is reduced by shrinking the noisy wavelet coefficient magnitudes. Shrinkage estimators can result from a Bayesian approach, which imposes a prior distribution on noise-free data. Recently, several wavelet based procedures for multicomponent images were proposed that account to some extent for the inter-component correlations, applying wavelet thresholding [7] or Bayesian estimation, using different prior models [8–11].

In this work, we propose to restore multicomponent images by a combination of a Fourier-based and a wavelet-based shrinkage, employing a 2-step procedure. In this study, we want to investigate whether it helps to account for the intercomponent correlations when restoring a multicomponent image. In case of the Fourier shrinkage step, a PCA transformation is applied prior to the Fourier transform, in order to decorrelate the image components, then a componentwise shrinkage is applied, after which the inverse PCA transform is applied. In the wavelet shrinkage step, a multicomponent procedure is applied, where the multicomponent wavelet coefficients are modeled using a multinormal or a multicomponent heavy-tailed distribution (a Gaussian Scale Mixture model). In both cases, the proposed procedure is compared to the single-component analogs, i.e. no PCA and single-component wavelet shrinkage. In the experimental section, conclusions on these comparisons will be drawn. The remaining of this manuscript is organized as follows. In the next section, the 2-step procedure is reviewed. In section 3, we introduce our multicomponent restoration procedure. Finally, section 4 is devoted to experiments and discussion.

Multicomponent Image Restoration

2 2.1

3

The 2-step Restoration Procedure The Image Model

The observed image sample Y of size M × M consists of an unknown signal S, first degraded by a circular convolution (denoted by ~) with a known impulse response H from a linear system H and then corrupted by zero-mean additive white Gaussian noise N , with variance σ 2 Y =H ~S+N = HS + N .

(1)

In the Fourier domain, this becomes Yf = Hf · Sf + Nf .

(2)

Sˆf = Hf−1 (Yf − Nf )

(3)

Inverse filtering

leads to amplified noise components due to singularities at the frequencies where Hf is (almost) zero. Application of a Wiener filter can prevent this. 2.2

Step 1: Fourier Shrinkage

Denote the result of the Wiener filter as X(= SˆWiener ), defined by Xf =

Hf∗ PS Yf |Hf |2 PS + PN

(4)

where PS and PN are the power spectra of the original image and the noise respectively. PS can be estimated by an iterative procedure [12], where an initial estimate of PS is inserted in (4). The initial value is set to PS = max(0, PY −PN ), with PY = |Yf |2 and, in the case of pure white Gaussian noise, PN = M 2 σ 2 . From the resulting Xf a new estimate of PS is set to |Xf |2 . This new estimate can be iteratively inserted in (4). By applying this procedure we obtain an estimate for Xf . Since the Fourier transform does not economically represent singularities, only a small amount of shrinkage is allowed in order to not distort the image edges. Due to this, the estimate X will still be noisy. For this, a wavelet shrinkage step is included. 2.3

Step 2: Wavelet Shrinkage (j,o)

The obtained Wiener filtered estimation X is wavelet transformed. Let xl denote its wavelet coefficient at spatial position l, resolution level j and orientation subband o. The corresponding wavelet coefficients of the noise-free image

4

Arno Duijster, Steve De Backer, and Paul Scheunders (j,o)

(j,o)

and the noise are sl and nl , respectively. Equivalent processing is typically applied to all the wavelet subbands, and hence we shall omit the indexes. Then x = s + n.

(5)

Using the Bayesian approach, a priori knowledge about the distribution of the noise-free data is assumed. The minimum mean squared error (MMSE) estimate is the posterior conditional mean R∞ R∞ Z ∞ s p(x|s)p(s)ds s φ(x − s; σn2 )p(s)ds −∞ sˆ = s p(s|x)dx = R ∞ = R−∞ . (6) ∞ p(x|s)p(s)ds φ(x − s; σn2 )p(s)ds −∞ −∞ −∞ Assuming e.g. a Gaussian prior for the noise-free signal p(s) = φ(s; σs2 ), the above MMSE estimate becomes the Wiener filter: sˆ =

σ ˆs2

σ ˆs2 x. + σn2

(7)

In this work, σ ˆs2 = max(0, σx2 − σn2 ) is calculated over the entire wavelet subband. The noise variance in each wavelet subband σn2 is in general a scaled version of the input image noise variance, where the scaling factors depend on the wavelet filter coefficients (see e.g. [13]). With the orthogonal wavelet families [14], that we use in this paper, the noise variance in all the wavelet subbands is equal to the input image noise variance.

3 3.1

The Multicomponent Restoration Procedure Multicomponent Fourier Shrinkage

To extend the restoration technique from the previous section to multicomponent images, we rewrite (1) in vector-matrix notation. Let S be the vector describing the multicomponent image S = [S 1 , S 2 , . . . , S K ] as a concatenation of multiple single component images S i with K the number of components. We can similarly define Y and N, being respectively the vectors containing the observed multicomponent image and the additive noise. If we assume the blurring operation to be equal for all components and no blurring among different components exist, then the multicomponent blurring operation can be written as H = IK ⊗ H with ⊗ the Kronecker product, IK the identity matrix of size K ×K and H the singlecomponent blurring operation. The multicomponent observation model can be written as Y = HS + N . (8) Following [1], the multicomponent image restoration can, depending on the assumptions made, be performed independently or in decorrelated fashion. In the independent case, all components are deblurred separately. In the decorrelated case, the assumption is made that the spatial correlations of each component, and the spatial cross-correlations between all pairs of components are equal, up

Multicomponent Image Restoration

5

to some scaling factor. We can write the total multicomponent covariance matrix as CS = Ccomp ⊗ CS , where Ccomp is the component covariance matrix of size K × K and CS the spatial correlation of one component, of size M 2 × M 2 . For the noise covariance we assume it to be spatially and component independent CN = σ 2 IM 2 K . Generally, it is assumed that the spatial covariance matrix is diagonalized by the discrete Fourier transform, with the diagonal elements described by the power spectrum PS (CS is a circulant matrix). The multicomponent restoration process can then be performed by first decorrelating the components, i.e. diagonalizing Ccomp , and then performing single-component restoration for each component. Let Upca be the orthogonal transformation diagonalizing the covariance matrix Ccomp by the operation generally referred to as principal component analysis (PCA) or Karhunen-Lo`eve transform T Upca Ccomp Upca = Λ

(9)

with Λ a diagonal matrix containing the eigenvalues of Ccomp on the diagonal. Applying the transformation on the multicomponent image, the covariance matrix CS is then transformed to (Upca ⊗ IM 2 )T CS (Upca ⊗ IM 2 ) = (Upca ⊗ IM 2 )T (Ccomp ⊗ CS )(Upca ⊗ IM 2 ) = Λ ⊗ CS .

(10)

A block diagonal matrix is obtained, removing all cross-component correlations. Applying the discrete Fourier transform on each component fully diagonalizes the multicomponent covariance matrix. 3.2

Multicomponent Wavelet Shrinkage

After Fourier shrinkage, each component of the estimated multicomponent image is wavelet transformed. The wavelet coefficients are concatenated into a vector, obeying x = s + n. (11) ˆ s of the unknown noise-free vector s is given The estimated covariance matrix C by ˆ s = Cx − Cn C (12) where Cx denotes the covariance matrix of the noisy vector x and Cn = σn2 IK is the noise covariance. In this paper, Cx is calculated over the entire image: Cx = hxl xTl il , where it is assumed that the intercomponent correlations are the same for all wavelet ˆ s is a covariance matrix, it needs to be semicoefficients of a subband. Since C positive definite. This is assured by performing an eigenvalue decomposition and clipping possible negative eigenvalues to zero.

6

Arno Duijster, Steve De Backer, and Paul Scheunders

The minimum mean squared error (MMSE) estimate, assuming a Gaussian prior for the noise-free signal p(s) = φ(s; Cs ), becomes the vector-based Wiener filter: ˆ s (C ˆ s + Cn )−1 x . ˆs = C (13) The Gaussian Scale Mixture (GSM) Prior The standard Wiener result is obtained using a multicomponent Gaussian prior model. It accounts for the multicomponent covariance, but it assumes that the marginal densities for the wavelet coefficients are Gaussian. It is well-known that this assumption is not justified: these marginals are symmetric and zero mean, but heavier tailed than Gaussians. Different other priors were proposed to better approximate the marginal densities. The Gaussian scale mixture (GSM) prior [15] models the probability density function p(s) by a mixture of Gaussians Z p(s) = p(s|z)p(z)dz (14) where p(z) is the mixing density, and p(s|z) is a zero mean Gaussian with covariance Cs|z . Under the GSM model, the additive noise model (11) becomes x=s+n=

√

zu + n

(15)

where both u and n are zero-mean Gaussians, with covariances given by Cu and Cn respectively. Then Cs|z = zCu , or, by taking expectations over z with E(z) = 1 , Cs = Cu . GSM densities are symmetric and zero-mean and heavier tailed than Gaussians. They are known to model the shape of the wavelet coefficient marginals more accurate than Gaussians. We apply the GSM to model multicomponent wavelet coefficients. In this way, the prior fully accounts for the inter-component covariances. The Bayes least squares estimate is given by Z ˆs = s p(s|x)ds ZZ ∞ = s p(s, z|x)dzds 0 ZZ ∞ = s p(s|x, z)p(z|x)dzds 0 Z ∞ = p(z|x)E(s|x, z)dz . (16) 0

Since the GSM model s conditioned on z is Gaussian, the expected value within the integral is given by a Wiener estimate E(s|x, z) = zCu (zCu + Cn )−1 x .

(17)

Multicomponent Image Restoration

7

The posterior distribution of z can be obtained, using Bayes’ rule, as p(z|x) = R ∞ 0

p(x|z)p(z) p(x|α)p(α)dα

(18)

with p(x|z) = φ(x; zCu + Cn ) . In [15], the authors motivate the use of the so-called Jeffrey’s prior [16] for the random multiplier z: p(z) ∝ z1 . We also refer to [10] and [15] for further information about the practical implementation of (16).

4

Experiments and Discussion

In all presented techniques, the noise variance is required. Instead of assuming that it is known, the noise variance is estimated as follows: let for each image component x1,HH denote the wavelet coefficients of the first resolution level and orientation subband HH. Then the diagonal elements of the noise covariance are estimated by the classical median estimator [17] σ ˆ2 =

median(|x1,HH |) . 0.6745

(19)

The noise is assumed to be uncorrelated between different image components. After the Fourier shrinkage step, the noise is no longer white, since different frequencies were scaled by different factors. Therefore, before the wavelet shrinkage step, we need to re-estimate the noise variance, at each wavelet scale and orientation separately. However, it can be derived directly from the Fourier shrinkage filter. Denote Wf the Fourier transformed shrinkage filter as estimated in (4): Hf∗ PS Wf = . (20) |Hf |2 PS + PN p Let W = F T −1 ( |Wf |2 ) be the inverse Fourier Transform of the square root (j,o) of the filter’s power spectrum density. Define wl as its wavelet coefficient at location l, scale j and orientation o. Then, at each scale and orientation sX (j,o) (j,o) σ ˆ =σ ˆ |wl |2 . (21) l

In the experimental section, we want to find out whether multicomponent restoration outperforms single-component restoration. For this, the proposed 2-step procedure is applied both in single-component mode as well as in multicomponent mode. The procedures are tested on 2 images of 3 components: the color image ‘Lena’ and the 3-band multispectral remote sensing image ‘San Diego’. The images are blurred using a 3 × 3 boxcar function. On both images, Gaussian white noise is simulated with standard deviations ranging from 4 to 40, with steps of 4.

8

Arno Duijster, Steve De Backer, and Paul Scheunders 33 SF SFSW MF MFMW1 MFMW2

32 31 30

PSNR (dB)

29 28 27 26 25 24 23

5

10

15

20 25 standard deviation of the noise

30

35

40

Fig. 1. The PSNR in function of the noise standard deviation σ for the restored image ‘Lena’.

The following procedures are applied: – SF: the single-component Fourier shrinkage step, without PCA and denoising, using the procedure described in section 2.2. – SFSW: the single-component Fourier shrinkage step, with single-component denoising, using (7). – MF: the multicomponent Fourier shrinkage step, without denoising, using the procedure described in section 3.1. – MFMW1: the multicomponent Fourier shrinkage step, with multicomponent denoising, using the multinormal model (13). – MFMW2: the multicomponent Fourier shrinkage step, with multicomponent denoising, using the multicomponent GSM model (16). In the figures 1 and 2, the results are shown as the PSNR versus the noise standard deviation. For each measurement, 10 experiments with different noise realizations are performed and the average results is shown. Error bars are negligible. When the multicomponent Fourier shrinkage step is applied, a gain of about a half to one dB is obtained, compared to the single-component Fourier shrinkage. This clearly indicates that the multicomponent Fourier shrinkage step does account for the correlations between the image components. Including the wavelet shrinkage step improves the image quality. When applying the single-component

Multicomponent Image Restoration

9

31 SF SFSW MF MFMW1 MFMW2

30 29

PSNR (dB)

28 27 26 25 24 23 22 5

10

15

20 25 standard deviation of the noise

30

35

40

Fig. 2. The PSNR in function of the noise standard deviation σ for the restored image ‘San Diego’.

wavelet shrinkage on top of the single-component Fourier shrinkage, a gain in image quality is observed that increases with higher noise variance. The multicomponent denoising models generate superior results compared to the singlecomponent model. This clearly indicates that, when applying the wavelet shrinkage step, it is favorable to account for the correlations between the image components. When a model is applied that better approximates the probability density of the wavelet coefficients (as in this case the GSM model), better results are obtained. Applying both steps in a multicomponent way is the most favorable for the image quality improvement. In figure 3a the third channel of a part of the multicomponent image ‘San Diego’ is shown. A 3 × 3 boxcar function is used to blur the image and Gaussian white noise with a standard deviation of 30 is added (figure 3b). The figures 3c and d show the results of the single-component Fourier shrinkage without and with denoising, respectively. The results of the multi-component Fourier shrinkage are shown in the last three subfigures: figure 3e shows the image without denoising, while the figures 3f and g show the denoised images, using the multinormal and the GSM model, respectively. In this work, the Fourier shrinkage step was applied with a standard Wiener filter. In [4], a procedure is described where the regularization of the Fourier and wavelet shrinkage are simultaneously optimized. When applying such an optimization, we expect even better results.

10

Arno Duijster, Steve De Backer, and Paul Scheunders

Acknowledgements This work was supported by the Flemish Interdisciplinary institute for Broadband Technology (IBBT), through the project ICA4DT, Image-based Computer Assistance for Diagnosis and Therapy.

References 1. Hunt, B.R., K¨ ubler, O.: Karhunen-Loeve multispectral image restoration, part I: Theory. IEEE Trans. Acoust., Speech, Signal Process. ASSP-32 (1984) 592–600 2. Zervakis, M.E.: Generalized maximum a posteriori processing of multichannel images and applications. Circuits, Systems, and Signal Processing 15 (1996) 233– 260 3. Galatsanos, N.P., Wernick, M.N., Katsaggelos, A.K.: Multi-channel Image Recovery. In: The Handbook of Image and Video Processing. Academic Press (2000) 161–174 4. Neelamani, R., Choi, H., Baraniuk, R.: Forward: Fourier-wavelet regularized deconvolution for ill-conditioned systems. IEEE Trans. Image Process. 52 (2004) 418–433 5. Figueiredo, M.A.T., Nowak, R.D.: An EM algorithm for wavelet-based image restoration. IEEE Trans. Image Process. 12 (2003) 906–916 6. Bioucas-Dias, J.: Bayesian wavelet-based image deconvolution: a GEM algorithm exploiting a class of heavy-tailed priors. IEEE Trans. Image Process. 15 (2006) 937–951 7. Scheunders, P.: Wavelet thresholding of multivalued images. IEEE Trans. Image Process. 13 (2004) 475–483 8. Benazza-Benyahia, A., Pesquet, J.C.: Building robust wavelet estimators for multicomponent images using Steins’ principle. IEEE Trans. Image Process. 14 (2005) 1814–1830 9. Scheunders, P., Driesen, J.: Least-squares interband denoising of color and multispectral images. In: Proc. IEEE Internat. Conf. Image Proc. ICIP. Volume 2., Singapore (2004) 985–988 10. Scheunders, P., De Backer, S.: Wavelet denoising of multicomponent images using a gaussian scale mixture model. In: Proc. Internat. Conf. on Pattern Recognition, ICPR, Hong Kong (2006) 11. Piˇzurica, A., Philips, W.: Estimating the probability of the presence of a signal of interest in multiresolution single- and multiband image denoising. IEEE Trans. Image Process. 15 (2006) 654–665 12. Hillery, A., Chin, R.: Iterative wiener filters for image restoration. IEEE Trans. Signal Process. 39 (1991) 1892–1899 13. Jansen, M.: Noise Reduction by Wavelet Thresholding. Springer-Verlag, New York (2001) 14. Mallat, S.: A Wavelet Tour of Signal Processing. Academic Press, London (1998) 15. Portilla, J., Strela, V., Wainwright, M.J., Simoncelli, E.P.: Image denoising using scale mixtures of gaussians in the wavelet domain. IEEE Trans. Image Process. 12 (2003) 1338–1351 16. Box, G.E.P., Tiao, C.: Bayesian Inference in Statistical Analysis. Addison-Wesley (1992) 17. Donoho, D.L., Johnstone, I.M.: Adapting to unknown smoothness via wavelet shrinking. Journal of the American Statistical Association 90 (1995) 1200–1224

Multicomponent Image Restoration

a.

b.

c.

d.

e.

f.

11

g.

Fig. 3. The third channel of the multicomponent ‘San Diego’ image. From left to right and top to bottom: (a) original scene; (b) blurred and corrupted with noise; (c) restored using SF, (d) SFSW, (e) MF, (f) MFMW1 and (g) MFMW2.