Mixed Poisson Gaussian noise reduction in ﬂuorescence microscopy images using modiﬁed structure of wavelet transform

Fluorescence microscopy is an important investigation tool of discoveries in the ﬁeld of biological sciences where the imaging phenomena are limited by the noise. This paper introduces the integration of biorthogonal wavelet ﬁlters along with mixed Poisson-Gaussian unbiased risk estimate (MPGURE) based subband adaptive thresholding function for the restoration of low photon count microscopy images. The proposed algorithm consists of four steps. In the ﬁrst step, variance stabilization transform along with a multi-scale Wiener ﬁltering approach is used to ﬁlter out the noise and blurring effect. In the second step, deconvolved images are further decomposed by the biorthogonal wavelet ﬁlters. The modiﬁed wavelet subband structure is used for the identiﬁcation of noisy and noise-free subbands. In the next stage, different noisy coefﬁcients are thresholded using the MPGURE-based thresholding operation. The different thresholded images are combined along with different optimum coefﬁcients. Finally, inverse variance stabilization transformation is applied to obtain the ﬁnal restored output. Performance of the proposed algorithm is tested on 14 different benchmark image data sets with performance evaluation measures like signal-to-noise ratio, peak signal-to-noise ratio, mean structural similarity index measur, and correlation coefﬁcient. Simulation results of the proposed algorithm claim better results than other state-of-the-art techniques.


INTRODUCTION
Fluorescence microscopy imaging is an enormously powerful research and diagnostics tool nowadays for microscopic discoverers which precisely empowers visualization of biological structures and sub-microscopic details at a super-resolved scale. The innovation and research of labelling subtle details with Fluorescence protein allow a deeper view on the cellular world of biology for the study of DNAs, cells, and organisms [1]. Now, the access to the dynamics of living cells, organisms, and tissues at a super-resolution scale is possible with Fluorescence microscopy imaging. Fluorescence microscopy imaging has several applications in the field of neurology, biomedical, and material science. In the field of modern cell biology, it is widely used to analyse biolog-This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. IET Image Processing published by John Wiley & Sons Ltd on behalf of The Institution of Engineering and Technology ical structures such as the protein-protein interactions and its sub-cellular localization. Now, the live-cell imaging is also possible with Fluorescent microscopy [2]. The biological researcher needs a clear noise-free image for interpretation and analysis purpose. The observed images under Fluorescent microscopy are degraded by mixed Gaussian-Poisson noise (MPG) due to several factors such as low light conditions, limited exposure time, and low signal-to-noise ratio (SNR) [3]. These problems can be handled by devoted noise removing algorithms which rely on the MPG model and concreted to the high noise level. The proposed method provides better denoised images used by biological researchers to make meaningful quantitative measurements about a living specimen.
wileyonlinelibrary.com/iet-ipr of a weighted fidelity and a sparse regularization term. It uses a wavelet frame transform on graphs to keep feature details. In [5], the authors have introduced a non-blind multi-frame superresolution reconstruction based on mixed Poisson-Gaussian noise (MPGSR). The variational Bayesian approach and application to image deblurring is introduced for MPG noise using posterior mean estimation [6]. Mostly image restoration algorithms designed for an additive white Gaussian noise (AWGN), which gives an insignificant performance under Poisson statistics [7,8]. Total variation methods [9] are efficient for piece-wise constant specimens. The TV-MAP (total variation maximum a posterior) approach [10] works in the Bayesian framework but results in the artefacts like an amplification of noise. The nonlocal filtering methods uses a blending of similar patches from the full image but suffers from major Halo artefact. The NLM (non-local mean) Poisson scheme [11] uses Newton's gradient descent method for optimization but computationally inefficient due to its parametric adjustments and non-tunable regularization. The non-local principle component analysis (NLPCA) scheme [12] uses patch clusters and principle component analysis (PCA)-based estimator along with optimization but produces over smoothing effect. The dictionary-based methods, such as image denoising via sparse and redundant representations [13] uses the training of a dictionary on its degraded content using the K-SVD (Kmeans singular value decomposition) algorithm. K-SVD [14] shows a great performance for a big database but a training set that limits the performance. The problem of image denoising addressed based on the non-negative matrix factorization and sparse representation over redundant dictionary [15]. In other applications, the multi-temporal recovery of quantitative data is the extension of dictionary learning approach [16]. Nonnegative matrix factorization and error correction methods are combined to recover the original information covered by the clouds and the accompanying shadows [17]. The BM3D (Blockmatching and 3D filtering) and block matching [18] uses a nonlocal wavelet-based filtering approach and suffers from major artefacts for low photon count statistics. The stochastically connected random field (SRF) [19] is an extension of block matching but unable to retain the fine structural details under Poisson statistics. In photon-limited microscopy applications, the restoration methods discussed above show sub-optimal performance under low input SNR. So interpretation from such inseparable and heavily degraded observations with low SNR is the most challenging problem. This demands and motivates us to design an effective restoration method to efficiently resolve the said problem.
This paper proposes a very pre-eminent image restoration technique with a subband adaptive thresholding scheme which is non-iterative. The variance stabilization transform (VST) and multi-scale Wiener filtering approach filters-out the effect of blurring caused by microscopic components. A biorthogonal wavelet transform (BWT)-based image decomposition and Gaussian model-based adaptive selection of decomposition levels provides different frequency components related to MPG noise. The proposed biorthogonal wavelet struc-ture is used to split the noisy and noise-free wavelet coefficients. The mixed Poisson-Gaussian unbiased risk estimate (MPGURE) which is a more accurate approximation to MSE (mean square error). Here, the minimization of MPGURE is done in order to evaluate the signal adaptive threshold. Noisy coefficients are thresholded with adaptive thresholding and further optimization framework provides a better estimate of noise-free image. We extensively validated our proposed method on 14 different image data sets which illustrates that, the efficiency of proposed method in terms of SNR, correlation coefficient (CR), and mean structural similarity index measure (MSSIM) is very good. Hence, the proposed algorithm is very efficient and provides a better reconstruction of an image.
This paper is organized as follows: Section 2 explores an overview of the proposed image restoration mechanism with a block diagram. The detailed qualitative and quantitative analysis along with the microscopic data set used are described in Section 3. Finally, in Section 4 overall outline has been summarized as conclusions.

PROPOSED METHODOLOGY
The detailed block diagram of the proposed image restoration mechanism is shown in Figure 1.

Image formation model
In Fluorescent microscopy imaging, the observed images are degraded by two major sources of noise, namely, the Poisson and Gaussian. However, due to electronic and intrinsic thermal vibrations of the microscopic components and lack of infinite precision for photo-detector produces distortion in the observed specimen. It can be modelled as an AWGN [9]. Besides, the highly random nature of photon arrivals causes image degradation which is modelled as Poisson noise [20]. The combined effect is termed as MPG. The observed specimen further degraded by point spread function (PSF) of the microscope [21]. We incorporated an MPG noise model for microscopy images under low photon count statistics [9]. The proposed algorithm considers the accurate modelling of image formation it is expected that it will also perform well in the real experiments. Observed image under the microscope can be represented as where Y is the observed image, X is the original noise-free image, H is the PSF of the microscope, * indicates the convolution operation, ℙ(.) is the realization of Poisson noise, and is the scaling factor which controls the strength of Poisson noise. ℕ(0, 2 ) represents the AWGN with 0 mean and variance 2 .

VST and multi-scale Wiener filtering approach
The VST is used as a pre-possessing step to remove the data-dependency of the noise and to get the stable noise variance from the whole image. It is a non-linear operation that transforms a Poisson distribution into an approximately standard Gaussian distribution. VST fixes the variance of an image affected by Poisson noise and makes it to constant [22].
The observed image Y has variance which is directly proportional to true image intensity. To stabilize the variance of Y , we used variance stabilization transformation as, Z = T (Y ) which makes Var (Z ) asymptotically constant irrespective of image intensity as, where b is the normalizing factor, sgn() is the sign function, '⋅' indicates multiplication, and c is the convergence rate. We have considered optimal value of the convergence rate c as 0.0625 [23]. The deconvolution with a multi-scale Wiener filtering approach filters-out the effect of blurring and to provide better superlative resolution [21,24].

Adaptive selection of decomposition level using Gaussian model and BWT
The approach used for selecting the decomposition level is given in Figure 2. The input image after first-level decomposition contains four different subbands, namely, LL, LH , HL, and HH. The LL 1 indicating the LL subband obtained after first-level decomposition. The Gaussian model-based threshold d t calculation is done for each LL 1 subbands. Here, the higher value of d t indicating the more the level of noise. If d t > T , where T is pre-defined threshold, then second-level decomposition of LL 1 subband is obtained. Again out four different subbands LL 2 can be considered for calculating the threshold d t . If the condition d t < T is satisfied then s is the accurate level of decomposition. Otherwise, we will go for the thirdlevel decomposition, and the process repeats until the condition satisfied.
The BWT which has the property of orthonormality and linear phase. The perfect reconstruction filters (h, g) and (h,ḡ) with high pass filter h and low pass filter g arranged in infinite cascading yields two scaling and wavelets functions. The scaling function Φ and inverse scaling functionΦ at time scale t can be expressed as

FIGURE 2 Selection of decomposition levels
where n is the different scale levels.
The wavelet function Ψ and inverse wavelet functionΨ can be expressed as Here, the Gaussian model is used to find out an exact number of decomposition levels required for wavelet decomposition. The LL subband having Gaussian distribution with a bell-shaped blob and its width is controlled by the standard deviation. Higher the value of standard deviation indicating more the level of noise. The further decomposition of an image reduces the value of standard deviation indicating lower the level of noise. Decomposition is done until the value of standard deviation is less than the pre-defined level of threshold. The different value of pre-defined level are considered to find an optimum value of the pre-defined threshold. The selection of optimum threshold is varying the value of thresholds from The value of threshold d t is calculated for different subbands with s = 1, 2, 3, … , N as where is the standard deviation indicating a spread of blob, 0 is the average value of LL subband, and LL(i, j ) indicates the value of frequency coefficients. The scaling function of the wavelet representation is given as where f (x, y) is the input image and Φ j 0 ,m,n (x, y) is the scaling function. The wavelet function of horizontal HL representation of images are given as where j ,m,n (x, y) is the wavelet function for HL subband.
The only odd number rows are selected for a noise-free estimate of HL subband is shown in Figure 4. The size of HL subband is M × N , after selecting alternating rows the size ofĤL becomes M ∕2 j × N where j is the level of decomposition. The wavelet function of vertical LH representation of images are After selecting alternative odd number of rows and columns, we are getting noise-free estimate as The even number of rows and columns are selected for noisy subbands indicated byL H andHL.
The HH subbands contain detail diagonal information along with noise so thresholding operation is directly performed on it. The wavelet function of diagonal HH representation is given as The three-level decomposition using traditional way and modified subband wavelet method The total image after decomposition can be represented as where i = H , V , D is the multi-resolution representation of scaling and wavelet coefficients. The wavelet and scaling function can be represented as The HL, LH , and HH are the subbands with scaling and wavelet function can be represented as, The noisy and noise-free estimates denoted by S 1 and S 2 are evaluated using Equations (11) and (12). The three-level wavelet decomposition in a traditional way and modified method is shown in Figure 5.

Subband adaptive threshold calculation
The approximation to MSE has been made by estimating MPGURE as it is impossible to evaluate MSE directly without knowledge of true image [25]. The MPGURE can be calculated as where N is the total number of pixels within subband,S is the estimated subband, S n is the noisy subband, ‖.‖ 2 2 indicates the square of L 2 norm,S − is the function of real valued subband, ' ⟨⋅⟩ ' indicates the inner product between vectored subbands, is the gain for Poisson process, is the standard deviation, and ∇ is the divergence operator.
The first-order Taylor's series expansion of estimated realvalued function subband can be calculated as,S − ≈S -S . Hence, the modified Equation (17) can be written as The evaluation of some terms in Equation (18) such asS − , S , and 2S become practically unrealistic. In photon-starved applications, the Poisson noise is the dominant source of noise so the effect of Gaussian noise can be guardedly neglected [9] with approximation of 2 = 1 [25,26], which results as where n (S n ) is the difference between noisy and noise-free subbands for nth wavelet coefficient can be expressed as where c is the error between nth estimated and noisy coefficients, s is the total number of coefficients. The best estimation for standard deviation̄is the median value of absolute deviation with noisy wavelet coefficients can be calculated as [27] where S 2 is the noisy estimate evaluated using Equation (12). Similarly, 2 can be calculated with HH subband using Equation (21). The S 2 and HH are identically independent subband. Hence, total standard deviation (̄) can be calculated as,̄= 1 + 2 . For microscopic images, the proposed method uses a universal threshold as Th = 3̄which is calculated empirically.
For every noisy coefficient S i j , we fixed a square kernel A i j at the centre of noisy coefficient, and the size of corresponding kernel A i j is kept constant as S × S where S is an odd number and positive. We have considered d i j as The estimation of noise-free wavelet coefficientsS i j can be represented asS , T s is the value of threshold which is initially set as Th, and S i j is the noisy wavelet coefficients. The error between noisy and noisefree subbands indicated by n (S n ) for nth wavelet coefficient is given as The difference between noisy and noise-free subbands for nth wavelet subband can be obtained by applying a threshold to individual wavelet coefficients, The first-order partial derivative of wavelet coefficient is calculated with Equations (24) and (26) as The ‖ n (S n )‖ 2 2 term used in Equation (31) can be calculated as The expectation of MSE in terms of n (S n ) can be obtained as where s is constant which depends upon the variance of subband, ‖.‖ 2 2 indicates the L 2 norm square for the inside term, E[.] is the expectation operator, and n (S n ) is the difference between noisy and noise-free subband expressed in Equation (20). The MPGURE estimate for subband can be approximated using Equations (27), (28), and (29) as This is the unbiased risk estimate for MSE between noisy and noise-free subbands. The expectation of MPGURE is same as expectation of MSE which completely relies on noisy observations [25]. The MSE can be calculated as The theoretical proofs of Equation (32) can be found easily in [25,26]. We calculated the subband adaptive threshold which minimizes the MSE. The optimal threshold can be obtained as where T s and d s indicating optimal values of threshold and neighbouring kernel, respectively. The estimated subbandsS 2 and diagonal subbandH H are considered for thresholding operation. First, we calculate the optimal threshold of T s and neighbouring kernel d s using Equation (33). All the frequency coefficients of every subband are shirked with optimal threshold and kernel using Equation (23).
The low-frequency coefficients are less susceptible to MPG hence excluded from thresholding [28]. The noise-free subbands obtained after performing threshold operation are further used for reconstruction.
The rows of S 1 andS 2 are alternatively combined in order to obtain S 5 which is an approximation of LH subband. The columns of S 1 andS 2 are alternatively combined in order to obtain the approximation for HL subband as where ⨀ indicates alternative combining of rows as the first row of S 1 is considered and then after second row ofS 2 is considered. So alternative selection of rows is considered. The approximation to LH subband is done as

Wavelet reconstruction and MPGURE minimization
The individual reconstructed images are obtained with a threelevel reverse BWT. The optimal coefficients obtained with MPGURE minimization. The final reconstructed output can be expressed as an affine combination of individual reconstructed outputs R 1 , R 2 , … , R l . The special case of linear expansion which is an affine expansion of reconstructions (AER) can be expressed aŝ whereR is final reconstructed output and l represents the total number of optimized coefficients 1 , 2 , … , l . The optimum coefficients are calculated by minimizing the MSE given as subjected to constraints as the sum of all optimized coefficients are equal to 1 where Y s is the noisy image. The final denoised image with inverse VST asŶ whereR is the reconstructed image andŶ is the final denoised image.

RESULTS AND DISCUSSION
The proposed algorithm is implemented in MATLAB R2018a. It is run on Intel i3, 2 GHz PC with 8 GB RAM and Windows 10, 64-bit operating system. We have compared the results of the proposed method with competing algorithms such as MPGSR [5], NLPCA [12], PNLW (Poisson Non-local Wiener estimator) [29], TV-MAP [30], PURELET [26], and Un-modified wavelet method (UW). For better comparison, we have provided results of proposed method with traditional way of decomposition is termed as UW. The results of proposed method with modified wavelet structure (proposed) are discussed in this section. We have validated our proposed method on 14 different image data sets which include images from different imaging modalities. The qualitative and quantitative analysis of the results illustrates, the proposed algorithm gives better denoising performance under MPG noise.

Description of the database and performance metrics
The performance of the proposed algorithm has been tested on different image data sets taken from the yeast resource center public image repository (YRCPIR) [31]. The data set contains 532,182 TIFF images for sub-cellular localization of proteins and classes of proteins. For simulation, we selected 14 different sets of microscopic images tagged with various types of Fluorescent proteins. Each data set contains more number of different microscopic specimen images limited to a particular type of Fluorescent labels. These images taken under different light conditions.
To provide a quantitative evaluation of the proposed scheme, we used the average value of peak signal-to-noise ratio (PSNR) [32], MSSIM [32], and CR [32] as an indirect performance metric to evaluate the performance of the proposed scheme.

Parameters tuning
For parameter setting, we have smoothened the images with a Gaussian blurring with variance 2 ∈ [0 5] which is the approximation PSF. We have a convoluted input image with different sizes of Gaussian kernels such as 3 × 3, 7 × 7, and 9 × 9 for the deconvolution process. For noise degradation, we have quantized incoming photons count into different levels which represent an observed image. The MPG noise level can be controlled with parameters of and 2 . The strength of Poisson noise can be controlled with ∈ [0 12]. The different values of has chosen to control the strength of the noise. The AWGN noise with mean equal to '0' and variance 2 represented as ℕ(0, 2 ) is added to noise-free image. We have chosen a small value of variance 2 as its impact is secondary in low photon count statistics. For performance comparison, we have applied a deconvolved image as an input to all the state-of-the-art methods except the PURE-LET method.

Quantitative analysis
The SNR results are shown in Figure 6(a) for different methods. It clearly illustrates when photon count is low, the proposed algorithm provides better SNR results than competing algorithms. As the photon count increases towards the higher value, the more improvement in the output SNR is observed and finally becomes constant. The comparison of proposed method along with different competing methods in terms of the correlation coefficient is given in Figure 6(b). The proposed algorithm provides a better value of correlation coefficient which is close to unity which is a better result than other competing algorithms. The MSSIM results from Figure 6(c) clearly illustrates, a proposed algorithm achieves better MSSIM improvement than other competing methods. The MSSIM improvement provided by all methods varies with the structural dynamics of microscopic images. The average input PSNR value for 14 different data sets labelled with different types of Fluorescent proteins are reported in Table 1. All other competing algorithms provide better results with improved PSNR but this improvement is only marginal. Among all comparative methods, the PURE-LET method provides better PSNR improvement than all other methods. However, the proposed method gives better PSNR improvement than the existing five state-of-the-art algorithms. The proposed algorithm with modified wavelet structures gives better results than traditional way of wavelet decomposition. The results of proposed method in all tables are highlighted in bold face. The average value of MSSIM and correlation coefficient for proposed method along with competing methods are reported in Tables 2 and 3. The proposed method having a better MSSIM and correlation coefficient results than five existing state-of-the-art methods. In Tables 1, 2

Qualitative analysis
The visual quality analysis results obtained for Fluocell, Subcellular microtube, Fluorescent protein, and Saccharomyces cerevisiae species are shown in Figure 7. The original and its corresponding noisy images are also shown for better comparison purpose. For Fluocell, results are presented in the first column of Figure 7. The NLPCA algorithm has maintains the structural details present in an image to some extent but more affected by  noise. The images obtained with NLPCA method are given in the third row of Figure Figure 7. The sub-cellular microtube contains a small circular nucleus along with tiny thread-like structural elements. The NLPCA and TV-MAP techniques failed to preserve significant tiny thread-like structural elements and more effect of noise is observed. The restored results of PNLW are better than other existing comparative methods. The PNLW method reduces the noise to some extent but still, noise is present in the images. The MPGSR techniques have failed to retain fine structural elements and more suspected of noise as shown in the sixth row of Figure 7. The PURE-LET algorithm results in the loss of the tiny structural elements of high biological significance. The restoration results of the proposed algorithm are shown in last row of Figure 7 which demonstrates, the proposed method efficiently pre-serves the tiny thread-like structural elements at the same time it removes the effect of mixed noise to a great extent.
Fluorescent protein images contain two circular-shaped small tags that are of biological significance. The third column of Figure 7 represents the restoration results of different methods for Fluorescent protein images. The results obtained with the NLPCA method are more susceptible to noise hence show the limitation. The results of PNLW technique are better than other competing techniques. However, a small amount of noise is present. The TV-MAP method lay down the noise for maintaining cell structure as shown in fifth row of Figure 7. The results obtained with MPGSR technique completely brighten the complete portion of the cell and more suspected of noise. It is shown in the sixth row of Figure 7. The PURE-LET method can remove noise effectively but its overblurring effect causes smoothing hence fails to preserve the finest constituents of high significance. The two circular-shaped Fluorescent tags are diminished by PURE-LET method. The proposed approach has shown the potential to preserve fine details and removing noise very efficiently. This can be easily visualized from last row of Figure 7 with two small dots on Fluorescent protein.
The subsequent results obtained for Saccharomyces cerevisiae species are shown in the fourth column of Figure 7. The NLPCA and TV-MAP techniques can maintain structural details but more susceptible to noise as shown in third and fifth row of Figure 7, respectively. The PNLW method preserves the smallest cellular constituents like three circular shape dots present in an image while losing the ability to remove noise, effectively as shown in the fourth row of Figure 7. The more improvement in denoising is given by the MPGSR method as shown in Figure 7. The PURE-LET approach has totally disappeared the fine structural details like three small dots as shown in sixth row of Figure 7. The results obtained in last row of Figure 7 shows, the proposed algorithm has preserved small three    dots along with microscopic sub-details. The results obtained for PHC-C2DH-U373 cell along with five different methods are shown in last column of Figure 7. All denoising results validate, the proposed algorithm has good edge preservation capability of cells while it retaining the low-intensity details.

CONCLUSION AND FUTURE WORKS
This paper puts forward a new image restoration scheme for highly degraded microscopic images. The proposed technique solves two important problems. First, a VST provides a constant variance of the noisy image. Multi-scale Wiener filtering technique reduces the effect of blurring caused by the PSF of microscope. Second, the BWT-based image decomposition and Gaussian model-based adaptive selection of decomposition levels are used which provides the detailed coefficients related to noise. The modified wavelet subband structure with MPGUREbased thresholding operation eliminates the effect of noise. Furthermore, MPGURE-based optimization framework also provides a better estimate of noise-free image. The performance of proposed technique is verified with different evaluation measures like SNR, PSNR, MSSIM, and CR. The simulation result shows, the proposed algorithm gives better results than other existing state-of-the-art algorithms. Here, the proposed algorithm can be optimized in future by more accurate modelling of the PSF of microscope for deconvolution process.