Image Denoising Using Gaussian Scale Mixtures in the Wavelet domain

1. Image Denoising Using Gaussian Scale Mixtures in the Wavelet domain Image Denoising Using Gaussian Scale Mixtures in the Wavelet domain Alex Cunningham Ben Clarke Dy narath Eang ECE6258 20 November 2008

2. Agenda • Algorithm Overview • Background • Wavelets • Steerable Pyramid Decomposition • Gaussian Scale Mixture • Bayesian Least Square Estimation • Image Denoising Details • Implementation Details • Detailed Algorithm • Detailed Matlab Implementation • Results • Comparison to other denoising techniques • Improvement schemes

3. Algorithm Overview General Denoising Algorithm: • Decompose Image into Steerable Wavelet Pyramid Subbands • For each subband: • Perform local Wiener Estimation • Denoise coefficients using Bayesian Least Squares Estimation • Reconstruct Denoised Subband • Reconstruct image from denoised wavelet subbands Paper Source • J. Portilla, V. Strela, M. Wainwright, and E. P. Simoncelli, “Image denoising using a scale mixture of Gaussians in the wavelet domain,” IEEE Trans. Image Processing, vol. 12, no. 11, pp. 1338–1351, November 2003. Additional Tools • Matlab Pyramid Tools (Computational Vision Lab at NYU)

4. Wavelet Overview Discrete Wavelet Transform (DWT) Continuous Wavelet Transform (CWT) CWT Decomposition: Images Source: Matlab Help

5. Steerable Wavelet Pyramids (Spyr) Overview • Spyrs Combine: • Multiple Levels of DWT • Multiple Orientations of Differential Filters • “Steerability” comes from the set of oriented filters • Can reconstruct image without aliasing (Circle) Image Source: E. Simoncelli, W. Freeman, “The Steerable Pyramid: A Flexible Architecture for Multi-Scale Derivative Computation,” 2nd Ann. IEEE Intl. Conf. on Image Processing. Washington, DC. October, 1995.

6. Gaussian Scale Mixture (GSM) Overview • GSM = A useful statistical model • Problem: Images are NOT spatially homogeneous or scale invariant and their statistics do NOT follow Gaussian distribution • Solution: Use GSM to model wavelet coefficients • Usefulness: Characterize the coefficients’ marginal distribution shape and strong correlation between the neighbors’ amplitudes • Form:

7. Bayesian Least Squares (BLS) Estimation Overview • BLS = A useful estimation method • Problem: Two-step empirical Bayes estimation is suboptimal • Solution: Use BLS to estimate the real coefficients. • Usefulness: Improves the empirical Bayes Estimator and it’s now a single step process

8. There are five main variables: w: noise subband y: observed noisy image subband x: model GSM subband u: estimated real signal subband z: scalar RV used as a multiplier The observation is the signal plus noise – these 3 RV’s are independent Similarly, the covariance matrices follow the same additive rule Taking E{z} = 1, then Eq(6a) can be rearranged Eq (5) Eq (6a) Eq (7) GSM and Covariance Matrices

9. BLS estimates the denoised coefficients by estimating at the center coefficients of all neighborhoods It needs two main things: E[Xc|y,z] p(z|y) Local Wiener Estimation helps compute the 1st main item The 2nd main item can be computed sequentially using Equations 13 and 14. Equations come directly from the author’s paper Eq (8) Eq (9) Eq (10) Eq (11) Eq (12) Eq (14) Eq (13) Bayes Least Squares and Local Wiener Estimation

10. Image Denoising Algorithm Full algorithm:* • 1) Decompose the image into subbands. • 2) For each subband (except the lowpass residual): • a) Compute neighborhood noise covariance, Cw, from the image-domain noise covariance. • b) Estimate noisy neighborhood covariance, Cy. • c) Estimate Cufrom Cw and Cy using (7). • d) Compute A and M (Section III-B: Local Wiener Estimate). • e) For each neighborhood: • i) For each value z in the integration range: • A) Compute E{xc|y,z} using (12). • B) Compute p(y|z) using (14). • ii) Compute p(z|y) using (13) and(4). • iii) Compute E{xc|y} numerically using (8). • 3) Reconstruct the denoised image from the processed subbands and the lowpass residual. (p. 1343) * The algorithm is copied directly from the journal paper, p. 1343.

11. Matlab Algorithm Tree denoise_demo.m Load image of interest Make copy of image and add Gaussian noise of known variance Define number of levels/bands for steerable pyramid (SPyr) decomposition Calculate bound extensions and PSDs Decompose noise-image and noisy image into SPyr structure Perform denoising (denoise_ctrl.m) Reconstruct image from denoised SPyr samples Calculate signal-to-noise ratios, display images denoise_demo.m pyr_decomp.m bound_extension.m denoise_ctrl.m pyr_reconstruct.m pyr_insert.m covariance_calculator.m local_wiener_estimator.m expect_estimator.m pyr_extract.m

12. Denoising Functions denoise_ctrl.m • Extract particular subband (both noise image and noisy image) • Calculate covariance matrices (covariance_calculator.m) • Perform (local_weiner_estimator.m) from covariance matrices to get eigenvalues and eigenvectors • Use eigenvalues/eigenvectors to estimate the real coefficients by find expectation values for each pixel in the neighborhood (expect_estimator.m) • Reinsert denoised subband into pyramid structure covariance_calculator.m • Create neighborhood matrices (for noise image and noisy images) assuming zero-mean Gaussian distributed • Compute signal covariance (Cu) as difference between observed covariance (Cy) and noise covariance (Cw) local_wiener_estimator.m • With assumption of AWGN, perform wiener calculations • Requires computation of eigenvalue/eigenvectors for covariance matrices • Result is that for various covariance matrices, we get wiener matrices expect_estimator.m • Perform equations (12) and (14) • From equation (14), perform equation (13) • Use results of equations (12) and (13) to compute E{xc|y} numerically. Eq (8)* Eq (12) Eq (13) Eq (14)

13. Statistical Results Portilla Basic 512x512 or 256x256, 8-bit grayscale images Scan σ across a range of standard deviations Compute Peak SNR (PSNR) for each of six images Last column σ(PSNR) is the estimated standard deviation of the results. Ours Compute same PSNR Compute Improvement SNR (ISNR) from noisy image to denoised image Portilla PSNR table Our SNR table +σ is the input standard deviation for the AWGN

14. Extreme case: σ (sigma) = 100 Moderate case: σ (sigma) = 10 SNR Improvement of 12.45dB SNR Improvement of 5.87dB Pictorial Results

15. Comparison • Extreme case: • σ (sigma) = 100 • Subjective ranking: • GSM-BLS (Portilla) • Wiener 5x5 • Discrete wavelet transform (DWT) • Wavelet packet (WP)

16. Improvement Scheme #1 • ‘‘Image Denoising using shiftable directional pyramid and CGSM’’ • Authors: An Vo, Truong Nguyen, & Soontron Oraintara • Publication date: May 2007 • Synopsis: Improvement in PSNR values for certain images (Barbara) was achieveable by incorporating CGSM and PDTDFB. • PDTDFB – Pyramidal Dual Tree Directional Filter Bank • CGSM – Complex Gaussian Scale mixture. • *SSIM = Structural Similarity

17. Improvement Scheme #2 • ‘‘Image Denoising using GSM model and BLS estimation over wider range of observation’’ • Author:Hyng Il Koo & Nam Ik Cho • Publication date: April 2008 • Synopsis: Subjective improvement in the quality of the images by considering a larger range of coefficients and using a new statistical model for the p(y|z). • CRF = Conditional random field • Sum-product algorithm = inference algorithm

18. References Original reference • J. Portilla, V. Strela, M. Wainwright, and E. P. Simoncelli, “Image denoising using a scale mixture of Gaussians in the wavelet domain,” IEEE Trans. Image Processing, vol. 12, no. 11, pp. 1338–1351, November 2003. Other denoising papers • A. Vo, T. Nguyen, S. Oraintara, “Image Denoising using Shiftable Directional Pyramid and Scale Mixtures of Complex Gaussians” ISCAS 2007, pp. 4000-4003, May 2007 • H. I. Koo and N. I. Cho, “Image Denoising Based on a Statistical Model for Wavelet Coefficients,” ICASSP 2008, pp. 1269-1272, April 2008. • E. Simoncelli, W. Freeman, “The Steerable Pyramid: A Flexible Architecture for Multi-Scale Derivative Computation,” 2nd Ann. IEEE Intl. Conf. on Image Processing. Washington, DC. October, 1995.