# Poisson Denoising for Astronomical Images

Authors: Fahad Shamshad; M. Mohsin Riaz; Abdul Ghafoor
Journal: Advances in Astronomy
 (2018)

---

## Abstract

A denoising scheme for astronomical color images/videos corrupted with Poisson noise is proposed. The scheme employs the concept of Exponential Principal Component Analysis and sparsity of image patches. The color space RGB is converted to YCbCr andK-means++ clustering is applied on luminance component only. The cluster centers are used for chromatic components to improve the computational efficiency. For videos, the information of both spatial and temporal correlations improves the denoising. Simulation results verify the significance of proposed scheme in both visual and quantitative manner.


---

## Body

## 1. Introduction

Recent advancements in astronomy and digital systems emphasize the development of more sophisticated image processing algorithm. Acquisition of astronomical images in low photon count region is a dominant source of Poisson noise. The images containing intensity dependent Poisson noise (having same mean and variance value) cannot be accurately modelled with Gaussian distribution.In [1], Poisson noise is converted into Gaussian noise with unit variance using Variance Stabilization Transform (VST). Further improvement is made by applying Bayesian multiscale likelihood models [2] and l1 penalization hypothesis testing [3]. However these schemes [1–3] provide significantly degraded results for high Poisson noise. Close form approximation for unbiased inverse of Anscombe transform yields increased error peak for Poisson noise [4, 5]. Structured dictionary learning approach exploiting the self-similarity of image patches [6] requires significant amount of computational time especially for color images and videos.In [7], multiscale Poisson denoising based on shrinkage operators to attain maximum l2 error is proposed. However l2 distance is not suitable for Poisson corrupted measurements. In [8] spatially adaptive total variation framework is proposed with split Bregman optimization for Poisson denoising to handle high computational load but yields poor results in low photon count regime which is usually the case with astronomical images. Directional lapped orthogonal transform overcomes the limitation of diagonal edges and textures which separable wavelets fail to denoise efficiently [9]. Computationally efficient technique for estimation of unknown noise parameters requires a blank image under same mechanical shutter [10]. In a nutshell, most of above-mentioned techniques either suffer from high computational complexity or give poor results in high Poisson noise.A denoising scheme for astronomical color images/videos corrupted with Poisson noise is proposed. The scheme first converts the color space from RGB to YCbCr and appliesK-means++ clustering on luminance component only. The same cluster centers are used for chromatic components to improve the computational efficiency. Framework of Exponential Principal Component Analysis (EPCA) [11] is employed that better suits high Poisson noise. The proposed algorithm addresses several drawbacks of existing schemes while ingeriting their main strengths. For instance, patch based approach is used to exploit redundancies in image that have been shown to improve performance of image denoising algorithms [6]. Similarly clustering algorithms are used in different denoising algorithms to group similar patches together [5]. Note that careful clustering has significant importance in case of high Poisson noise (where only few photons are acquired by detector). Similarly taking into account temporal correlation of videos also improves the performance of restoration algorithms [12, 13]. Simulation results on different images and videos verify the significance of proposed scheme in both visual and quantitative manner.


## 2. Proposed Technique

LetNβ be the observed image from image acquisition device having dimensions P×Q, where p=1,2,…,P are rows, q=1,2,…,Q are columns, and β∈red, green, blue are color bands. For lϵ{1,…,PQ} let Nβ[l] (column stack representation of Nβ) be the observed pixel at location i whose entries are independent Poisson random variables with true value Iβ[l]>0 which is to be estimated. Given the clean image Iβ[l] likelihood of observing noisy Poisson image Nβ[l] can be written as(1)PNβl∣Iβl=exp⁡-IβlIβlNβlNβl!Noisy image is converted into luminance chrominance (YCbCr) color space by matrix transformation as(2)NYNCbNCr=0.2990.5870.114-0.147-0.2890.4360.615-0.515-0.100NRNGNBLet Nα denote the noisy image in YCbCr color space where αϵ{Y,Cb,Cr}; then matrix MY of vectorized patches is constructed from the NY component by extracting all overlapping patches O of size E×E as(3)MY=N1,1YN1,2Y⋯N1,EYN2,1YN2,2Y⋯N2,EY⋮⋮⋱⋮NO,1YNO,2Y⋯NO,EYwhere Ni,jY represents jth pixel in ith patch for 1≤i≤O and O=(P-E+1)(Q-E+1) (considering image border issues) and let Ii,jY be defined for clean image as Ni,jY for noisy one. Similarly matrices MCb and MCr are built for NCb and NCr, respectively.K-means++ algorithm (having more accuracy and efficiency as compared to conventional k-means) is then applied on matrix MY to obtain clusters c=1,2,…,C and denote Sc for the size of cth cluster. As NY (luminance component) contains most of the valuable information (edges, shades, texture, etc.) of image the same cluster centers are then assigned to Cb and Cr components for computational efficiency of algorithm.Different algorithms have used PCA to project the extracted noisy image patches to low dimensional space. LetIic,jα be the jth pixel of ith patch in cluster c; using PCA it can be written as(4)Iic,jα≈ϕψic,jαwhere ϕ is coefficient matrix of size Sc×L and ψ is dictionary atom matrix of size L×E. Sc is cluster size (patches in cluster c) and L is number of principal axes to be retained. However, in this paper, the framework of [11] is used as it is more suitable for exponentially distributed data. Equation (4) can be written as(5)Iic,jα≈exp⁡ϕψic,jαwhere the exponential operator is element-wise.After neglecting the constant term and applying maximization of negative Poisson likelihood alongwith approximation of (5), the loss function for cth cluster is(6)LossU,V=minϕ,ψ⁡∑i=1Sc∑j=1Eexp⁡ϕψic,jα-Nic,jαϕψic,jαEquation (6) is solved using fast Newton method [14], where ϕ and ψ are initialized randomly. Then ith rows of ϕ (denoted as ϕic,:(t)) at iteration t) for cluster c are updated as(7)ϕic,:t+1=ϕic,:t-exp⁡ϕtψtic,:-Nic,:ψTtϕtXiϕTt-1where Xi=diag(exp(ϕ(t)ψ(t))ic,1,exp(ϕ(t)ψ(t))ic,2,…, exp(ϕ(t)ψ(t))ic,E) is diagonal matrix of size E×E for cluster c and α is omitted for simplicity. Similarly columns of ψ (denoted as ψ:,j(t)) are updated as(8)ψ:,jt+1=ψ:,jt-ϕTt+1Djϕt+1-1ϕTt+1exp⁡ϕt+1ψt:,j-N:,jwhere Dj=diag(exp(ϕ(t+1)ψ(t))1c,j,diag(exp(ϕ(t+1)ψ(t))2c,j+1,…, exp(ϕ(t+1)ψ(t))Sc,E).Row ofϕ and columns of ψ are updated iteratively. Stopping criteria for proposed algorithm are either number of iterations (maximum 20) or a defined error threshold (0.01), whichever is reached first. Equation (6) after solving reduces to(9)Φ^,Ψ^=argminϕ,ψ⁡LossU,Vand denoising estimate for cluster c is (I is clean original image while I^ is denoised image)(10)I^ic,jα=exp⁡Φ^Ψ^ic,jαAfter denoising of all c clusters, denoised patches are projected onto pixels. As each pixel has multiple presentation in different patches, it is done by uniformly averaging all estimates of those patches containing given pixel. Denoising image I^α is then converted to RGB domain I^β as(11)I^RI^GI^B=101.1401-0.395-0.58112.0320I^YI^CbI^CrThe above denoising methodology is extended for Poisson noise corrupted astronomical videos by taking spatial correlation into account. LetNtα be Poison corrupted tth frame of size P×Q in YCbCr color space where t=1,2,…,F. If all overlapping cubes of size E×E×Λ are extracted then total number of cubes will be (P-E-1)(Q-E-1)(F-Λ-1), where Λ represents number of frames. For simplicity denote (P-E+1)(Q-E-1) by O and (F-Λ+1) by Z (total cubes then become O×Z). After concatenating every cube row-wise, patch by patch, in single row matrix form (Ni,j,tα), jth pixel of ith cube in frame t is written as (12)Ni,j,tα=N1,1,1αN1,2,1α⋯N1,E,1αN1,E+1,2α⋯N1,E∗Λ,ΛαN2,1,1αN2,2,1α⋯N2,E,1αN2,E+1,2α⋯N2,E∗Λ,Λα⋮⋮⋯⋮⋮⋯⋮NO,1,1αNO,2,1α⋯NO,E,1αNO,E+1,2α⋯NO,E∗Λ,Λ+1αNO+1,1,2αNO+1,2,2α⋯NO+1,E,2αNO+1,E+1,3α⋯NO+1,E∗Λ,Λ+1α⋮⋮⋯⋮⋮⋯⋮NOZ-1,1,ZαNOZ-1,2,Zα⋯NOZ-1,E,ZαNOZ-1,E+1,Z+1α⋯NOZ-1,E∗Λ,FαNOZ,1,ZαNOZ,2,Zα⋯NOZ,E,ZαNOZ,E+1,Z+1α⋯NOZ,E∗Λ,FαNote that indext is changed after every Eth pixel as these are total entries on front side of cube (size E×E). K-means++ clustering is then applied on these cubes and clusters are denoised using (6) (cubes are vectorized in this case). After denoising these cubes are projected onto pixels by averaging and converted back to RGB color space.


## 3. Results and Discussion

The proposed and state-of-the-art existing Poisson denoising techniques are applied on various astronomical images/videos (taken using Hubble telescope). For images, the proposed technique is compared with K-SVD [15], CBM3D [16], PLOW [17], Anscombe transform [18], and NLPCA [5]. The results on videos are compared with KSVD-3D [12], CBM4D [13], and PURE-LET [19].Figures1–3 show visual comparison of existing and proposed techniques on three astronomical images (Disk, Nebula, and Jupiter) with different noise peaks. It can be observed that the proposed technique is able to reconstruct a better output image as compared to existing techniques (even in low noise peak values). Tables 1 and 2 provide quantitative analysis of existing and proposed techniques using Peak Signal to Noise Ratio (PSNR) and Structural Similarity Index Measure (SSIM) [20]. Note that the proposed technique provides better quantitative values as compared to the state-of-the-art existing techniques.Table 1
PSNR comparison.
Techniques Peak = 0.2 Peak = 0.6 Peak = 1 Disk Infrared Jupiter Disk Nebula Jupiter Disk Nebula Jupiter K-SVD [15] 17.766 19.134 12.750 20.001 22.872 16.997 22.961 23.964 20.278 CBM3D [16] 16.905 20.596 16.59 19.989 21.361 19.243 23.331 22.008 22.205 Plow [17] 20.410 17.135 19.50 23.165 20.674 21.091 24.251 23.767 23.04 Anscombe [18] 18.931 21.072 15.818 23.242 23.521 18.742 25.799 23.889 21.112 NLPCA [5] 21.844 23.216 19.705 24.129 22.172 22.379 26.737 24.801 24.212 Proposed 22.198 24.081 20.392 25.032 24.689 22.614 27.015 25.036 24.627Table 2
SSIM comparison.
Techniques Peak = 0.2 Peak = 0.6 Peak = 1 Disk Infrared Jupiter Disk Nebula Jupiter Disk Nebula Jupiter K-SVD [15] 0.661 0.495 0.362 0.625 0.579 0.428 0.826 0.756 0.569 CBM3D [16] 0.519 0.556 0.451 0.623 0.692 0.525 0.775 0.787 0.645 PLOW [17] 0.732 0.457 0.356 0.779 0.589 0.509 0.839 0.722 0.635 Anscombe [18] 0.731 0.660 0.469 0.774 0.668 0.584 0.844 0.785 0.661 NLPCA [5] 0.752 0.669 0.476 0.788 0.818 0.571 0.853 0.825 0.683 Proposed 0.796 0.681 0.491 0.812 0.699 0.599 0.862 0.839 0.698Figure 1
Disk: noisy input, reference, and denoised images (for peak = 0.2).
(a)
Input (b)
Reference (c)
CBM3D [16] (d)
Plow [17] (e)
Anscombe [18] (f)
K-SVD [15] (g)
NLPCA [5] (h)
ProposedFigure 2
Nebula: noisy input, reference, and denoised images (for peak = 0.2).
(a)
Input (b)
Reference (c)
CBM3D [16] (d)
Plow [17] (e)
Anscombe [18] (f)
K-SVD [15] (g)
NLPCA [5] (h)
ProposedFigure 3
Jupiter: noisy input, reference, and denoised images (for peak = 1).
(a)
Input (b)
Reference (c)
CBM3D [16] (d)
Plow [17] (e)
Anscombe [18] (f)
K-SVD [15] (g)
NLPCA [5] (h)
ProposedFigure4 shows denoising result of proposed and existing techniques with cube size 20×20×3 and noise peak = 1. The reconstructed frames (in Figures 4(i) and 4(j)) using proposed technique are quite similar to the reference frames (in Figures 4(k) and 4(l)). Figure 5 shows the PSNR and SSIM values of individual frames for different noise peaks. It can be observed that the proposed technique provides better quantitative values almost for all frames as compared to the state-of-the-art existing techniques. Table 3 shows the average quantitative results (PSNR and SSIM) of video frames (311–325). For Figure 5, the Saturn video created from Hubble images taken over about 9-hour span is used. Note that the PSNR and SSIM of only 10 frames (from 315 to 325) are shown in the paper, as scenes in video do not change much due to slow motion of Saturn. Consequently, the PSNR and SSIM values of frames are almost constant.Table 3
Average PSNR (dB) and SSIM values of 311–325 frames.
Techniques PSNR SSIM Peak = 0.2 Peak = 0.6 Peak = 1 Peak = 0.2 Peak = 0.6 Peak = 1 KSVD-3D [12] 19.428 22.286 23.718 0.4621 0.6381 0.7221 CBM4D [13] 21.341 23.644 24.692 0.5211 0.6766 0.7620 PURE-LET [19] 20.055 22.145 23.417 0.4892 0.6212 0.6121 Proposed 22.426 24.534 25.528 0.6252 0.7599 0.8312Figure 4
Frames 311 and 318 of a Saturn video.
(a)
Noisy 311 (b)
Noisy 318 (c)
PURE-LET [19] (d)
PURE-LET [19] (e)
KSVD-3D [12] (f)
KSVD-3D [12] (g)
CBM4D [13] (h)
CBM4D [13] (i)
Proposed (j)
Proposed (k)
Reference 311 (l)
Reference 318Figure 5
(a)–(c) PSNR graphs for peaks 0.2, 0.6, and 1; (d)–(f) SSIM graphs for peaks 0.2, 0.6, and 1.
(a)
Peak = 0.2 (b)
Peak = 0.6 (c)
Peak = 1 (d)
Peak = 0.2 (e)
Peak = 0.6 (f)
Peak = 1


## 4. Conclusion

A Poisson denoising scheme in nonlocal framework using EPCA based on the approach of Gaussian mixture models is proposed for astronomical imaging. The images patches sparsity and dictionary learning approach are unified. The method employs EPCA (which is more suitable for Poisson noise). The scheme first converts the color space from RGB to YCbCr and appliesK-means++ clustering on of luminance component only. The same cluster centers are used for chromatic components to improve the computational efficiency. The cubes (in case of color videos) of luminance components (information of both spatial and temporal correlations) are used for improved denoising. Simulation results verify the significance of proposed scheme in both visual and quantitative manner.



---

*Source: 2417939-2018-06-10.xml*