Poisson-Gaussian Noise Reduction Using the Hidden Markov Model in Contourlet Domain for Fluorescence Microscopy Images

In certain image acquisitions processes, like in fluorescence microscopy or astronomy, only a limited number of photons can be collected due to various physical constraints. The resulting images suffer from signal dependent noise, which can be modeled as a Poisson distribution, and a low signal-to-noise ratio. However, the majority of research on noise reduction algorithms focuses on signal independent Gaussian noise. In this paper, we model noise as a combination of Poisson and Gaussian probability distributions to construct a more accurate model and adopt the contourlet transform which provides a sparse representation of the directional components in images. We also apply hidden Markov models with a framework that neatly describes the spatial and interscale dependencies which are the properties of transformation coefficients of natural images. In this paper, an effective denoising algorithm for Poisson-Gaussian noise is proposed using the contourlet transform, hidden Markov models and noise estimation in the transform domain. We supplement the algorithm by cycle spinning and Wiener filtering for further improvements. We finally show experimental results with simulations and fluorescence microscopy images which demonstrate the improved performance of the proposed approach.


Introduction
Digital images are prevalent in our lives as a result of advances in multimedia, internet, computers, and wide spread of portable imaging devices such as consumer digital cameras and camcorders. Image acquisition and processing is also common in medical technology and the service industry. Such demands have consequently led to the advancement of image sensors such as the charge coupled device (CCD) and the complementary metal oxide semiconductor (CMOS). Owing to the development of image sensor hardware, increased spatial resolution has resulted in a decreased sensor pixel size, which causes the expansion of the photon noise effect [1]. In many image applications such as fluorescence microscopy and astronomy, only a limited number of photons can be collected due to various physical constraints such as a light source with low power to avoid phototoxicity, and short exposure time. Therefore, these based on cross validation estimation [29]. Lingenfelter et al. presented the optimal penalty function to produce sparser images for maximum-likelihood solution with Poisson data [30]. Recently, a sparsity-regularized convex optimization algorithm for Poisson noise has been developed [31], and Salmon et al. improved this approach using non-local Principal Component Analysis for Poisson noise [32]. Giryes and Elad proposed a Poisson denoising method using sparse representation modeling by Salmon's method and dictionary learning [33].
Although the Poisson noise removal framework has been exploited in many publications as stated above, there has been very little research on mixed Poisson-Gaussian noise removal. Boulanger et al. proposed an effective denoising method of 3D fluorescence microscopy images for Poisson-Gaussian noise [34]. This method employs the generalized Anscombe transform to stabilize the variance of Poisson-Gaussian noise and achieves noise reduction using a minimizer consisting of an objective non-local energy functional involving spatio-temporal image patches. Poisson-Gaussian unbiased risk estimate-linear expansion of thresholds (PURE-LET) by Luisier et al. is a recent study on this issue [2][19] [35]. This methodology optimizes a thresholding algorithm in the transform domain of the undecimated wavelet transform (UWT) and block discrete cosine transform (BDCT) to denoise images corrupted by mixed Poisson-Gaussian noise. PURE-LET minimizes the noise according to the data-adaptive unbiased estimation of the mean squared error (MSE) by a non-Bayesian framework. This algorithm plays an important role in treating Poisson-Gaussian noise directly. However, PURE-LET approximates the unbiased MSE estimation to minimize the cost function. In this paper, we propose an effective noise removal algorithm for Poisson-Gaussian noise using HMM and contourlet transform. The Poisson-Gaussian distribution model is used to estimate the noise parameters of an image, and these parameters are applied to our proposed denoising algorithm for noise reduction.

Materials and Methods
Fluorescence Microscopy image data set The first set of images was acquired from a Nikon C1 Plus confocal laser microscope at the Medicinal Bioconvergence Research Center at Seoul National University. The data set contained 100 samples with a 512x512 size of fixed HeLa cells, labeled with three fluorescent dyes: Alexafluor555 in the red channel, Alexafluor488 in the green channel, and DAPI in the blue channel. The average of 100 images was used as the baseline for PSNR calculation.
The second data set was obtained from a Nikon A1R confocal laser microscope at the Department of Life Science of Ewha W. University. The data set contained 40 images of fixed HeLa cells of 512x512 size, labeled with two fluorescent dyes: golgin97 in the green channel, and DAPI in the blue channel. The average of 40 images was used as the baseline also.

Noise modeling and estimation
Noise can be classified into signal-dependent noise and signal-independent noise. Signaldependent noise is modeled by a Poisson distribution which is obtained from photon counting and signal-independent noise is normally modeled by a Gaussian distribution [36].
In this section, we follow the Poisson-Gaussian noise modeling of Foi et al. [1]. The observed signal z can be represented as the sum of the original signal y and noise as where x 2 X is the pixel position in the image domain X. The noise model is represented with two mutually independent parts, the signal-dependent Poisson component η p , and the signalindependent Gaussian component η g .
Assume that χ(y(x) + η p (y(x))) has the Poisson distribution, where χ > 0 is a scale factor. The probability distributions of Poisson and Gaussian noise are denoted as follows: wðyðxÞ þ Z p ðyðxÞÞÞ $ PðwyðxÞÞ; and Z g ðxÞ $ Nð0; bÞ; ð2Þ where χ > 0 and b ! 0 are real scalar parameters and P(Á) and N(Á) represent the Poisson and normal distributions. The mean and the variance of the Poisson distribution are the same. Therefore, the variance can be obtained from the intrinsic property of a Poisson distribution as follows: Since E{η p (y(x))} = 0 and χ 2 var{η p (y(x))} = χy(x), the variance can be represented as var {η p (y(x))} = y(x) / χ from Eq (3). Therefore, the Poisson noise component η p has a variance proportional to the signal y(x). That is var{η p (y(x))} = ay(x), where a = χ −1 . On the other hand, the variance of Gaussian noise η g is constant and it will be represented as b.
As a result, the overall noise variance of z in Eq (1) has an affine form, The standard deviation σ becomes sðyðxÞÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ayðxÞ þ b p . In this paper, the minimum and the maximum value of brightness y(x) are normalized to 0 and 1, respectively. Thus the standard deviation ranges from sð0Þ . The Poisson noise parameters are estimated using the noise estimation algorithm in [1] [37]. The noise parameter estimation process is as follows. Firstly, the image in the wavelet domain is segmented into smooth regions using level sets. Secondly, local mean and variance are estimated for each uniform region after the wavelet analysis. Finally, the optimal parametersâ and b are estimated by maximum likelihood (ML) estimation.
The noise variance of z is therefore calculated using Eq (4) and the estimated parametersâ andb.

Contourlet Transform
Signal modeling in the wavelet domain has been employed for an effective image noise reduction. Most natural images, however, are composed of not only discontinuous corners but also smooth curves. While the wavelet can express corners, it cannot capture the smoothness along contours in a compact manner. Hence, new algorithms such as curvelet and contourlet transforms have been developed to improve such shortcomings [38][39]. The curvelet transform, which was developed for better noise removal, consists of a 2-D rotation operation and a frequency domain division based on the polar coordinates [38]. As a result, information of the various directional features in images can be obtained to retain curves and edges compactly. The decomposition into the curvelet is simple in the continuous domain, but becomes problematic in the discrete domain. Geometrical problems such as significant bias in horizontal and vertical directions appear in the generalized rectangular-sampling grid. To overcome such difficulties, Do and Vetterli proposed the contourlet transform, which has similar effectiveness as the curvelet transform, for the direct use in the discrete domain [39].
The contourlet transform consists of two filters: the Laplacian pyramid (LP) filter and the directional filter bank (DFB). The LP filter is used in the existing wavelet filter as a low pass filter to separate high and low frequencies. The directional filter bank provides information on image direction components as described in Fig 1. After passing the LP filter, images are partitioned into high and low frequencies. While the high frequency component is divided according to direction in the DFB, the low frequency component is, after down sampling, divided into low and high frequency bands by the LP filter. DFB filtering is recursively applied to high frequency components. Since the contourlet transform can extract multi-resolution and multidirectional information, the boundaries of natural images can be described effectively. Typically, the LP filter proposed by Burt and Adelson is used [40]. The LP filter produces band-pass images by generating image differences between the down-sampled images and the original images. Although the biggest drawback of the LP filter is over sampling, the filter produces desired band-pass images without mixing frequency components at each pyramid level. While the wavelet filter can mix frequencies in the high-band channel after down sampling, the LP filter does not mix frequencies since the images that only pass the low frequency filter are downsampled [41]. In many cases, the DFB filter is based on the design of Bamberger and Smith [42]. The directional filter produces wedge-shaped frequency division subbands for 2 ℓ at ℓ-level binary tree decomposition.
To evaluate the sparsity between the wavelet transform and the contourlet transform, the histograms of their coefficients are compared with respect to the Lena image after adding Poisson-Gaussian noise (S1 Fig). The wavelet filter and directional filter used for the comparison are 'Daubeichies8' and 'dmaxflat5', respectively. The filter type does not practically affect the performance. Fig 2 shows that the sparsity of the contourlet transform is much higher than that of the wavelet transform. The ratio of zeros in the contourlet coefficients is 4.83%, while that in the wavelet coefficients is 1.88%. Zeros in the contourlet coefficients is almost twice as many as that in the wavelet coefficients.
The contourlet filter bank achieves continuous domain expansion at L 2 ð< 2 Þ using the contourlet transform. The connection between discrete contourlet transform and continuous domain expansion can be established with a new multiresolution analysis structure in a similar method of the link between the wavelet transform and filterbank [43].

The Contourlet HMM for Poisson-Gaussian Noise Reduction
Signal modeling, as well as the noise pdf, is important for effective noise removal. Existing algorithms have applied a generalized Gaussian [25] or Gaussian scale mixtures [6] for signal modeling. In this paper, the HMM method proposed by Crouse et al. is adopted to model the signal pdf for its efficacy [3]. In this section we summarize the HMM signal modeling and present our extension for Poisson-Gaussian noise.
The HMM denoising method by Crouse et al. is based on the wavelet transform and the signal-independent noise model. In this paper, the contourlet transform is employed, whereas the Crouse's method uses the wavelet transform. The contourlet transform results in better sparsity than the wavelet transform as shown in Fig 2. The marginal distribution of natural images in contourlet domain is highly non-Gaussian. Furthermore, the dependency relationship exists between the parent and child coefficients in the contourlet transform, which is an important reason for applying the HMM. As shown in Fig 3(A), the parent coefficient has its four children in two separate directional sub-bands, while parent-children relationship in the wavelet domain is always in the same direction. Fig 3(B) plots the absolute values of the parent and child coefficients in log scale. It demonstrates the dependency between the parent and the child coefficients of the Lena image.
The contourlet HMM utilizes the clustering and persistence properties of contourlet transform coefficients. Clustering and persistence are additional useful properties of contourlet transform besides major properties such as locality, multiresolution and energy compaction. Clustering is the property where adjacent coefficients tend to be large (small), if a contourlet coefficient value is large (small). Persistence shows that the values of contourlet coefficients are very likely to propagate across scales; if a contourlet coefficient is large (small) in one scale, then the parent coefficient is large (small).
Two statistical models can describe contourlet properties in HMM. First, the independent mixture model is used to decompose the marginal probability of each coefficient as a mixture density with a hidden state variable to reflect the non-Gaussian property of the contourlet coefficients. Second, probabilistic graphs are used to represent the dependencies between contourlet coefficients. The hidden Markov chain model captures the horizontal dependencies within each scale and the hidden Markov tree (HMT) represents the vertical dependencies across scale.
The mixture model is composed of two-state zero-mean gaussians: the pdfs of the state variable S, p S (1) and p S (2) = 1 − p S (1), and the variances of Gaussian pdfs in each state. In most of the mixture model applications, the contourlet coefficient W is observable, however, the state variable S cannot be observed, thus appropriately called hidden.
In general, an M-state Gaussian mixture model of random variable W is given as follows.
1. a discrete random state variable S which has the values s 2 1, 2, Á Á Á, M, according to the pdf p S (s).
2. the Gaussian conditional pdfs p W|S (w | S = s), s 2 1, 2, Á Á Á, M. Hence, the probability density function p W is a weighted sum of the Gaussian conditional pdfs: The HMM parameters for the M-state Gaussian mixture model for each contourlet coefficient W i are 1. p S 1 ðmÞ: the pdf for the root node S 1 2. ε mr i;rðiÞ ¼ p S i jS rðiÞ ½mjS rðiÞ ¼ r: the conditional probability that S i is in state m when the given S ρ(i) is in state r 3. μ i,m and s 2 i;m : the mean and variance of contourlet variable W i when the given S i is in state m.
These parameters can be grouped by model parameter vector θ. To estimate the HMM parameters for given data, the expectation maximization (EM) algorithm is employed. The training data w is incomplete; the complete data (w, s) is composed of the training data and the hidden state s. The goal is then to maximize the incomplete log-likelihood function ln p(w | θ), where the EM algorithm performs this challenging maximization with an iterative process separately taking two simple steps, E step and M step. The E step calculates the expected value E S [ln p(w, S | θ) | w, θ l ] in the l th iteration. Then the M step maximizes the log-likelihood function of θ to acquire θ l+1 for the next iterative process. Once the parameters are obtained, then the Bayesian noise removal process can be applied. After the contourlet coefficient w k i of the noise signal and the state s k i are obtained, the conditional estimation for the contourlet coefficient v k i of the original signal is as follows: where s 2 n i is the noise variance in contourlet domain. The conditional mean estimation of v k i can be obtained by using the hidden state probabilities pðS k i jw k ; yÞ as by-products of the EM algorithm through the chain rule of conditional expectation.
Finally the denoised signal is acquired by the inverse transform of the estimated contourlet coefficient.
In what follows, we describe the denoising method for Poisson-Gaussian noise. While the variance of Gaussian noise is constant, the variance of Poisson-Gaussian noise is related to the signal intensity. In other words, each pixel has different noise variance and thus it cannot be controlled using the traditional HMM. Therefore, we extend the contourlet HMM to deal with Poisson-Gaussian noise (PG-HMM). The parameters of Poisson-Gaussian noise a and b defined in Section 2, which are the noise estimates in the image domain, can be estimated by the noise estimation method in [1]. The variance of the noise component in the image domain in Eq (1) is denoted as s 2 e i , where e(x) = η p (y(x)) + η g (x) denotes the sum of Poisson and Gaussian noise which can be calculated by Eq (4). The original signal y i is approximated by the low-pass filtered value of the noisy image. We then need to estimate the noise variance in contourlet domain. Let z i and h i denote the observed pixel values and the contourlet filter coefficients, respectively. The noise component in contourlet domain n i is given by where e i is the noise component in the image domain, h i is the contourlet filter coefficient, and Ã is the convolution operator. Since the noise in contourlet domain is a linear combination of image noise, we need to derive an equation for contourlet noise variance s 2 n i . The noise variance of pixel i, denoted as s 2 e i , is equal to E½e 2 i , since E[e i ] = 0. Now, we derive the noise variance in contourlet domain, s 2 n i , using Eq (10) as follows: Here, the cross product terms of (∑ e i−k h k ) 2 can be eliminated because neighboring noise components are assumed to be statistically independent. Thus, we obtain As shown in Eq (12), the noise variance in contourlet domain is obtained by filtering the estimated noise variances in the image domain using the square of the contourlet filter coefficients. Finally, the denoised image can be obtained through Eq (9) using the HMM parameters and the Poisson-Gaussian noise variance in contourlet domain.
To improve noise removal performance, the Wiener filtering and cycle-spinning methods are cascaded. The Wiener filter is designed to minimize the mean squared error. Cycle spinning for noise removal is a simple and efficient method which can be applied to a shift variant transform [44,45]. The equation for cycle spinning is as follows: S Ài;Àj ðT À1 ðD ½TðS i;j ðzÞÞÞ; ð13Þ where S i,j represents a 2-D circulant shift with i and j shifts in horizontal and vertical directions respectively. K 1 and K 2 are the total number of horizontal and vertical shifts, respectively. T represents a shift variant transform which is the contourlet transform in our case and D represents the presented contourlet transform based HMM denoising algorithm. The variableŷ is the noise removed image obtained by the cycle spinning method. Noise removed images can be obtained repeating 2-D circular shift by K 1 K 2 times. Consequently, the cycle spinning method improves the peak signal to noise ratio (PSNR) by averaging the noise reduced images. A sketch of the GP-contourlet HMM denoising algorithm is illustrated in Fig 4 and presented in Algorithm 1. While the conventional HMM algorithm does not deal with signaldependent noise, our proposed algorithm reduces Poisson-Gaussian noise successfully. Experimental results in the following section show that the proposed method shows the best performance for many images corrupted by strong Poisson noise both visually and in terms of PSNR values.
1. Compute GP noise parameters a and b using Eq (6).
2. Compute the variance in the image domain using Eq (4).

Apply the contourlet transform.
Cycle spinning 4. for i = 1 to K 1 , j = 1 to K 2 do 5. Compute the variance in the transform domain using Eq (11) and the parameters, a and b, obtained at step 1.

Wiener filtering
10. Apply Wiener filtering to the denoised signal for further noise reduction. 11. Obtain the final denoised signal.

Denoising experiments
Implementation of Denoising Algorithm. The Laplacian pyramid (LP) filter and the Directional filter bank (DFB) filter of the contourlet transform used in the experiments are as follows: we adopted 'Daubeiches8' for the LP and 'dmaxflat5' for the DFB [46]. We confirmed experimentally that the numbers of scales and directions in contourlet domain do not influence the performance significantly. We therefore divided the images into five scales and each scale into four directions. Grayscale images Camera man, Lena, Boat, Barbara and Peppers were used as test images with different noise levels. The images are available in S1 File. We first confirmed the validity of the derived Poisson noise variance in contourlet domain, and then evaluated the performance of the PG contourlet HMM noise reduction algorithm using synthetic and real data. Quantitative Evaluation Measure. As a measure of quality, we used the peak signal-tonoise ratio (PSNR), defined as PSNR ¼ 10log 10 ðI 2 max =MSEÞ, where I max is the maximum intensity of the original image and MSE is the mean squared error. The output MSEs of the denoised images were averaged over 10 realizations of pseudo-random noise. The standard deviation of noise was estimated from the noisy image using the noise estimation algorithm proposed by Foi et al. [1].

Results and Discussion
We evaluate the performance of the PG-HMM noise reduction algorithm and apply the proposed method to confocal microscopy images in this section.

The verification of estimated noise variance
To verify the accuracy of estimated noise variance, we performed a computer simulation using the 50th row of the Lena image as an example. We compare the noise variance calculated from multiple realizations of noise, which will be used as a ground truth, the noise variance estimated in the image domain and the theoretical noise variance predicted by Eq (12). The noise variance estimated in the image domain was obtained by Eq (4). In practice, the accuracy of the sample noise variance depends on the number of noise generations. The plots of simulated noise variance illustrated in Fig 5 are results of 100 times and 1000 times of noise generations. The blue solid line, green dashed line, red dotted line and cyan dash-dot line indicate the noise variance by Eq (12), the sample noise variance by Eq (11) from 100 times of noise generation, the sample noise variance by Eq (11) from 1,000 times of noise generation, and the noise variance estimated in the image domain, respectively. As shown in the plot, the result is more accurate as the number of the noise generation increases. In addition, Fig 5 demonstrates that the variance predicted by Eq (12) closely matches the noise variance estimated from multiple noise generation while it is significantly different from the image domain noise variance from Eq (4). The validity of the proposed Poisson noise estimation can be confirmed regardless of test images.

Comparisons with existing denoising algorithms
In this section, we compare our proposed method with various denoising algorithms to show the effective denoising of our method on Poisson-Gaussian noise. In Table 1, we evaluate the performance of each block shown in Fig 4. Here, the total noise variance is fixed and only the ratio of Poisson and Gaussian noise is altered. The PSNRs of HMM based on the Gaussian noise are degraded as the Poisson component is increased. As shown in Table 1, the contourlet transform is more effective than the wavelet transform, and the cycle spinning improves the performance. Next, we compare our approach with the PURE-LET algorithm, which is a state-of-the-art denoising algorithm for Poisson-Gaussian noise, and BM3D, which is one of the most effective methods in Gaussian-based denoising algorithms. We provide comparative results of our method with the BM3D algorithm combined with the generalized Anscombe transform. The experimental results for low photon counts are shown in Table 2. Our method shows the best noise reduction for Cameraman and Peppers in low count images, while BM3D is better for Boat and Lena by a small margin. Fluorescence microscopic images are the results of photonlimited imaging, thus the performance of the algorithm for low-count images is important. In this respect, the experimental results are quite encouraging. For subjective comparison of image quality, we present the denoising results of the test images degraded by simulated Poisson noise with peak intensity 3 and Gaussian noise with standard deviation 0.

Application to fluorescence microscopy images
In this section, we present denoising results of photon-limited fluorescence microscopy images. Our proposed method and the BM3D method are applied to two fluorescence microscopy image sets. The first set of images was acquired from a Nikon C1 Plus confocal laser microscope at the Medicinal Bioconvergence Research Center at Seoul National University. The data set contained 100 images with a 512x512 size of fixed HeLa cells, labeled with three fluorescent dyes: Alexafluor555 in the red channel, Alexafluor488 in the green channel, and DAPI in the blue channel. The average of 100 images is used as the baseline for PSNR calculation. The data set is available in S2 File. The visual quality of the cell image is evaluated from Fig 8. Fig 8A and  8B presents the baseline and the observed data with single acquisition, respectively. Fig 8C and  8D shows the denoised images from the BM3D and the proposed method, respectively. As  Table 3. Our proposed method performs better than the BM3D method, except for the blue channel. The BM3D is preferable in the blue channel as it is particularly effective with periodic textures or flat regions.
The second data set was obtained from a Nikon A1R confocal laser microscope at the Department of Life Science of Ewha W. University. The data set contains 40 images of fixed HeLa cells of 512x512 size, labeled with two fluorescent dyes: golgin97 in the green channel, and DAPI in the blue channel. An average of the 40 images is used as the baseline, which is presented in Fig 9(D). The data set is available in S3 File. To evaluate the denoising performance  for different noise levels, we obtained three image sets with different laser intensities as shown in Fig 9A-9C. The denoising results for the single acquisition image with laser intensity 0.4 are presented in Fig 9D-9F. The PSNR results are provided in Table 4. While the proposed algorithm is more effective than the BM3D when the laser intensity is weak, the BM3D is more effective than our proposed method when the laser intensity increases.

Conclusions
In this paper, an effective denoising algorithm for mixed Poisson-Gaussian noise in low-count images is presented. We applied the contourlet transform for sparse representation of signal, and adopted the HMM for effective signal modeling in the transform domain. The contourlet  transform not only separates of images into high and low frequency components but also provides information about the directional components in the images using the Laplacian pyramid filter and the directional filter bank. The HMM algorithm adopts an independent mixture model to match the non-Gaussian nature of the contourlet coefficients and adopts hidden Markov models to characterize the key dependencies between the contourlet coefficients. Furthermore, this method estimates optimal HMM parameters using the EM algorithm. The Poisson-Gaussian noise variance in contourlet domain is obtained by filtering the noise variance of each pixel with the square of the contourlet filter coefficients. Using the estimated HMM parameters of the signal and noise variances, the signal-dependent noise is reduced through Bayesian estimation. We finally show the experimental results with simulations and fluorescence microscopy images and demonstrate the improved performance of the proposed approach when photon count is limited through extensive comparisons with the traditional Bayesian estimator and state-of-the-art techniques. We demonstrate that the performance of the proposed method, which is based on accurate source pdf modeling with HMM and Gaussian mixture, is comparable to those of high performance denoising methods such as BM3D combined with the modified Anscombe transform or PURE-LET. Our approach has the following advantages. First, noise variance in the transform domain is estimated more accurately using the convolution. Second, the noise reduction performance is superior in very low-count images. Third, the proposed method show good performance for the images with irregular patterns such as cell images. Although the computation time is not low compared with existing methods, it can be shortened by modifying HMM training process. The denoising performance of the proposed method can be improved further by incorporating the nonlocal means algorithm, since our method is based on a point-wise technique. The proposed method can be modified for feature detection or segmentation as well as for denoising 1-D signals.