Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

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

  • Sejung Yang,

    Affiliation Ewha Institute of Convergence Medicine, Ewha Womans University Medical Center, Seoul, Republic of Korea

  • Byung-Uk Lee

    bulee@ewha.ac.kr

    Affiliation Department of Electronics Engineering, Ewha Womans University, Seoul, Republic of Korea

Abstract

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 applications acquiring images by photon counting have extremely low signal to noise ratio [2]. The aim of this study is to model and remove the noise in a low-count image.

The two predominant sources of noise in digital image acquisition are (a) the stochastic nature of the photon-counting process at detectors, and (b) the intrinsic thermal and electronic fluctuations of the acquisition devices. Traditional denoising algorithms have employed the additive white Gaussian noise modeling to account for the second source of noise, which is signal-independent. Robbins proposed an empirical Bayesian framework to estimate Gaussian noise [3]. Lee, likewise, proposed a two-step empirical Bayesian estimation [4] which estimates the variance of signal from the neighbors of an observed pixel and applies the standard linear least squares (LLS) solution. Malfait and Roose further exploited a methodology to realize the Bayesian approach by applying the Markov random field [5] and Crouse et al. proposed an algorithm which uses hidden Markov models (HMM) to obtain the variance of signal and therefore, to denoise with Bayesian estimation [6]. Simoncelli published a research paper on the Bayesian denoising process [7], and Mihcak et al. applied a maximum a posteriori (MAP) estimation based on exponential marginal prior [8]. The Bayesian least squares-Gaussian scale mixtures (BLS-GSM) algorithm developed by Portilla et al. consists of a multivariate estimator resulting from Bayesian least-squares optimization, assuming Gaussian scale mixtures as a prior for neighborhoods of coefficients at adjacent positions and scales [9]. Block Matching 3D (BM3D), regarded as one of the state-of-the-art methods, improves sparsity by grouping 2D image blocks into 3D data sequences and decreases noise effectively by collaborative Wiener filtering [10]. These algorithms show reasonable performance for signal-independent Gaussian noise. However, the reduction in the size of an image sensor has resulted in a higher impact of signal-dependent noise. Consequently, the mixed Gaussian-Poisson noise model was developed [1]. Since Poisson noise is signal-dependent, it cannot have a constant noise variance, which causes difficulties in designing a denoising algorithm. To overcome the difficulties, variance stabilizing transforms (VSTs) such as the Anscombe transform [11] and the Fisz transform [12] have been introduced. Donoho exploited the Anscombe VST in denoising applications [13]. Later, many noise removal techniques employed the VST with the denoising algorithms based on the Gaussian noise reduction algorithm for Poisson noise removal [1419]. Although the Anscombe VST makes the variance of Poisson noise constant and allows the denoising algorithms based on Gaussian noise to perform effectively, it shows disappointing results in the case of low-count images [1719]. This is because of the inverse transformation rather than the stabilization itself. Bias error is unavoidable after an inverse transformation since the Anscombe transform is a nonlinear transformation. To overcome this bias problem, Mäkitalo and Foi proposed the exact unbiased inverse of the generalized Anscombe transform (GAT) for Poisson-Gaussian noise [2022]. In this paper, we tested the combination of GAT and the existing noise reduction method based on Gaussian noise for comparison with the proposed algorithm.

Another effective method for Poisson noise reduction is an algorithm using the Haar transform. Poisson data after the Haar transform have a binomial distribution for children coefficients given a parent coefficient. This property has been exploited not only in user-calibrated hypothesis testing [23] but in the Bayesian framework as well [18] [2426]. Since the Gaussian noise reduction methods yield higher quality images with the differentiable wavelet transform than with the Haar transform or other shift-invariant transformations, researchers have begun to modify Poisson-intensity estimation in order to incorporate the multi-scale transformation. Kolaczyk developed a shrinkage method using corrected hard/soft threshold based on an arbitrary wavelet transform to handle the nature of burst-like Poisson intensities [27], and Charles and Rasson generalized this approach to operate for various kinds of Poisson data [28]. Nowak and Baraniuk proposed a wavelet shrinkage method where the threshold is locally estimated 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. Signal-dependent 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 (1) where xX 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 signal-independent 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: (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: (3)

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, (4)

The standard deviation σ becomes . 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 to .

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 are estimated by maximum likelihood (ML) estimation. (5) where . Eq (5) can be rewritten as (6)

The noise variance of z is therefore calculated using Eq (4) and the estimated parameters and .

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 [3839].

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 multi-directional 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.

thumbnail
Fig 1. Contourlet Filter bank consisting of Laplacian Pyramid (LP) and Directional Filter Bank (DFB) [39].

https://doi.org/10.1371/journal.pone.0136964.g001

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.

thumbnail
Fig 2. The comparison of sparsity between wavelet and contourlet transforms on the Lena image.

(a) wavelet, (b) contourlet.

https://doi.org/10.1371/journal.pone.0136964.g002

The contourlet filter bank achieves continuous domain expansion at 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.

thumbnail
Fig 3. Parent-child relationship of contourlet coefficients.

(a) a parent coefficient and its pertinent four child coefficients (b) parent-child dependency of contourlet coefficients in log scale.

https://doi.org/10.1371/journal.pone.0136964.g003

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, pS(1) and pS(2) = 1 − pS(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 ∈ 1, 2, ⋯, M, according to the pdf pS(s).
  2. the Gaussian conditional pdfs pW|S (w | S = s), s ∈ 1, 2, ⋯, M.

Hence, the probability density function pW is a weighted sum of the Gaussian conditional pdfs: (7)

The HMM parameters for the M-state Gaussian mixture model for each contourlet coefficient Wi are

  1. : the pdf for the root node S1
  2. : the conditional probability that Si is in state m when the given Sρ(i) is in state r
  3. μi,m and : the mean and variance of contourlet variable Wi when the given Si 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 ES[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 of the noise signal and the state are obtained, the conditional estimation for the contourlet coefficient of the original signal is as follows: (8) where is the noise variance in contourlet domain.

The conditional mean estimation of can be obtained by using the hidden state probabilities as by-products of the EM algorithm through the chain rule of conditional expectation.

(9)

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 , 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 yi is approximated by the low-pass filtered value of the noisy image. We then need to estimate the noise variance in contourlet domain. Let zi and hi denote the observed pixel values and the contourlet filter coefficients, respectively. The noise component in contourlet domain ni is given by (10) where ei is the noise component in the image domain, hi 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 .

The noise variance of pixel i, denoted as , is equal to , since E[ei] = 0. Now, we derive the noise variance in contourlet domain, , using Eq (10) as follows: (11)

Here, the cross product terms of (∑ eikhk)2 can be eliminated because neighboring noise components are assumed to be statistically independent. Thus, we obtain (12)

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: (13) where Si,j represents a 2-D circulant shift with i and j shifts in horizontal and vertical directions respectively. K1 and K2 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 K1K2 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 signal-dependent 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.

thumbnail
Fig 4. Block diagram for the proposed method (CT: contourlet transform).

https://doi.org/10.1371/journal.pone.0136964.g004

Algorihm 1. GP-contourlet HMM.

  1. Compute GP noise parameters a and b using Eq ().
  2. Compute the variance in the image domain using Eq ().
  3. Apply the contourlet transform.

Cycle spinning

  1. 4. for i = 1 to K1, j = 1 to K2 do
  2. 5. Compute the variance in the transform domain using Eq () and the parameters, a and b, obtained at step 1.

Expectation Maximization

  1. 6. Estimate HMM parameter vector using EM algorithm.

Bayesian estimation

  1. 7. Estimate denoised contourlet coefficients using Eq ().
  2. 8. end for
  3. 9. Apply the inverse contourlet transform to recover the denoised image.

Wiener filtering

  1. 10. Apply Wiener filtering to the denoised signal for further noise reduction.
  2. 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-to-noise ratio (PSNR), defined as , where Imax 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.

thumbnail
Fig 5. Plots of noise variances in the wavelet domain for the 50th row of the Lena image; noise variances by Eq (12) (blue solid line), sample noise variances from 100 times of noise generation (green dashed line), 1,000 times of noise generation (red dotted line) and noise variance estimated in the image domain (cyan dash-dot line).

https://doi.org/10.1371/journal.pone.0136964.g005

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.

thumbnail
Table 1. Denoising performance evaluation of each block shown in Fig 4 in terms of PSNR (dB) for various noise ratio of Poisson and Gaussian components with fixed noise variance.

https://doi.org/10.1371/journal.pone.0136964.t001

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 photon-limited 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.3 in Figs 6 and 7. Although all three methods have similar PSNRs, Fig 7 shows that the proposed method provides fewer visual artifacts with a smoother face and preserves edges better than the PURE-LET or BM3D techniques. Our method restores the nose and mouth of the face for a more natural look and does not demolish the boundary as shown in Fig 6. To give the reader an indication of the typical computation times, for Cameraman (512x512) the various denoising methods require time approximately as follows: BM3D with GAT 2.5 s, UWT/BDCT PURE-LET 13.5 s, and the proposed method 10.2 s. The compared algorithms were implemented in MATLAB and BM3D is a pre-compiled execution file. These results are obtained with an Intel Core i7-3770 processor, running at 3.4 GHz.

thumbnail
Fig 6. Denoising results on Lena image (http://en.wikipedia.org/wiki/File:Lenna.png.) with a = 1/3 and b = 0.32.

(a) original image, (b) noisy image (7.48 dB), (c) BM3D (24.80 dB), (d) PURE-LET (24.68 dB), (e) the proposed method (24.72 dB).

https://doi.org/10.1371/journal.pone.0136964.g006

thumbnail
Fig 7. Denoising results on the magnified Lena image with a = 1/3 and b = 0.32.

(a) original image, (b) noisy image, (c) BM3D, (d) PURE-LET, (e) the proposed method.

https://doi.org/10.1371/journal.pone.0136964.g007

thumbnail
Table 2. Comparison of the proposed method with the BM3D method and PURE-LET in terms of PSNR (dB) for various low-count cases.

https://doi.org/10.1371/journal.pone.0136964.t002

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 observed in Fig 8, the proposed algorithm significantly reduces the level of noise and still recovered the details in the green channel. The comparison of PSNR between the proposed method and the BM3D method in each channel is presented in 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.

thumbnail
Fig 8. Denoising results on the first set of HeLa cell images.

(a) an average of 100 images used as the ground truth, (b) single acquisition image, (c) BM3D, (d) the proposed method.

https://doi.org/10.1371/journal.pone.0136964.g008

thumbnail
Table 3. Comparison of the proposed method with the BM3D method in terms of PSNR (dB) for each channel on the first set of HeLa cell images.

https://doi.org/10.1371/journal.pone.0136964.t003

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.

thumbnail
Fig 9. The second set of HeLa cell images and denoising results.

First low: the second set of HeLa cell images with different laser intensities. (a) 0.2, (b) 0.4, (c) 0.8. Second low: denoising results on HeLa cell image with laser intensity 0.4. (d) an average of 40 images used as the ground truth, (e) BM3D, (f) the proposed method.

https://doi.org/10.1371/journal.pone.0136964.g009

thumbnail
Table 4. Comparison of the proposed method with the BM3D method in terms of PSNR(dB) for various laser intensities on the second set of HeLa cell images.

https://doi.org/10.1371/journal.pone.0136964.t004

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.

Supporting Information

S1 File. Test image data set.

Camera man (Figure A), Boat (Figure B), Barbara (Figure C) and Peppers (Figure D)

https://doi.org/10.1371/journal.pone.0136964.s002

(ZIP)

S2 File. The first data set of HeLa cell images.

https://doi.org/10.1371/journal.pone.0136964.s003

(ZIP)

S3 File. The second data set of HeLa cell images.

https://doi.org/10.1371/journal.pone.0136964.s004

(ZIP)

Acknowledgments

The authors would like to thank Professor Dong-min Kang with Ewha Womans University for providing fluorescence microscopy images used in this paper.

Author Contributions

Conceived and designed the experiments: SY, BL. Performed the experiments: SY. Analyzed the data: SY, BL. Wrote the paper: SY, BL. Programmed the software used in analysis: SY.

References

  1. 1. Foi A, Trimeche M, Katkovnik V, Egiazarian K. Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data. IEEE Trans. Image Process. 2008; 17: 1737–1754. pmid:18784024
  2. 2. Luisier F, Blu T, Unser M. Image denoising in mixed Poisson–Gaussian noise. IEEE Trans. Image Process. 2011; 20: 696–708. pmid:20840902
  3. 3. Robbins H. The empirical Bayes approach to statistical decision problems. Ann. Math. Stat. 1964; 1–20.
  4. 4. Lee JS. Digital image enhancement and noise filtering by use of local statistics. IEEE Trans. Pattern Anal. Mach. Intell. 1980; PAMI-2(2): 165–168.
  5. 5. Malfait M, Roose D. Wavelet-based image denoising using a Markov random field a priori model. IEEE Trans. Image Process. 1997; 6: 549–565. pmid:18282948
  6. 6. Crouse MS, Nowak RD, Baraniuk RG. Wavelet-based statistical signal processing using hidden Markov models. IEEE Trans. Signal Process. 1998; 46: 886–902.
  7. 7. Simoncelli EP. Bayesian denoising of visual images in the wavelet domain. Bayesian inference in wavelet-based models: Springer; 1999. 291–308.
  8. 8. Mihcak MK, Kozintsev I, Ramchandran K, Moulin P. Low-complexity image denoising based on statistical modeling of wavelet coefficients. IEEE Signal Process. Lett. 1999; 6: 300–303.
  9. 9. Portilla J, Strela V, Wainwright MJ, Simoncelli EP. Image denoising using scale mixtures of Gaussians in the wavelet domain. IEEE Trans. Image Process. 2003; 12: 1338–1351. pmid:18244692
  10. 10. Dabov K, Foi A, Katkovnik V, Egiazarian K. Image denoising by sparse 3-D transform-domain collaborative filtering. IEEE Trans. Image Process. 2007; 16: 2080–2095. pmid:17688213
  11. 11. Anscombe FJ. The transformation of Poisson, binomial and negative-binomial data. Biometrika 1948; 246–254.
  12. 12. Fisz M. The limiting distribution of a function of two independent random variables and its statistical application. Institute of Mathematics Polish Academy of Sciences. 1955. 138–146.
  13. 13. Donoho DL. Nonlinear wavelet methods for recovery of signals, densities, and spectra from indirect and noisy data. Providence: American Mathematical Society. 1993. 173–205.
  14. 14. Fryzlewicz P, Nason GP. A Haar-Fisz algorithm for Poisson intensity estimation. J Comput. Graph. Stat. 2004; 13: 621–638.
  15. 15. Jansen M. Multiscale Poisson data smoothing. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 2006; 68: 27–48.
  16. 16. Zhang B, Fadili J, Starck J-L, Olivo-Marin J-C. Multiscale Variance-stablizing Transform for Mixed-Poisson-Gaussian Processes and its Applications in Bioimaging; 2007. 233–236.
  17. 17. Zhang B, Fadili JM, Starck J-L Wavelets, ridgelets, and curvelets for Poisson noise removal. IEEE Trans. Image Process. 2008; 17: 1093–1108. pmid:18586618
  18. 18. Lefkimmiatis S, Maragos P, Papandreou G. Bayesian inference on multiscale models for Poisson intensity estimation: Applications to photon-limited image denoising. IEEE Trans. Image Process. 2009; 18: 1724–1741. pmid:19414285
  19. 19. Luisier F, Vonesch C, Blu T, Unser M. Fast interscale wavelet denoising of Poisson-corrupted images. Signal Process. 2010; 90: 415–427.
  20. 20. Mäkitalo M, Foi A. Optimal inversion of the Anscombe transformation in low-count Poisson image denoising. IEEE Trans. Image Process. 2011; 20: 99–109. pmid:20615809
  21. 21. Mäkitalo M, Foi A. A closed-form approximation of the exact unbiased inverse of the Anscombe variance-stabilizing transformation. IEEE Trans. Image Process. 2011; 20: 2697–2698. pmid:21356615
  22. 22. Mäkitalo M, Foi A. Optimal inversion of the generalized Anscombe transformation for Poisson-Gaussian noise. IEEE Trans. Image Process. 2013; 22: 91–103. pmid:22692910
  23. 23. Kolaczyk ED. Nonparametric estimation of intensity maps using Haar wavelets and Poisson noise characteristics, Astrophys. J. 2000; 534: 490–505.
  24. 24. Timmermann KE, Nowak RD. Multiscale modeling and estimation of Poisson processes with application to photon-limited imaging. IEEE Trans. Inf. Theory 1999; 45(2): 846–862.
  25. 25. Kolaczyk ED, Dixon DD. Nonparametric estimation of intensity maps using Haar wavelets and Poisson noise characteristics. The Astrophys. J. 2000; 534: 490.
  26. 26. Lu H, Kim Y, Anderson JM. Improved Poisson intensity estimation: denoising application using Poisson data. IEEE Trans. Image Process. 2004; 13: 1128–1135. pmid:15326854
  27. 27. Kolaczyk ED. Wavelet shrinkage estimation of certain Poisson intensity signals using corrected thresholds. Stat. Sin. 1999; 9: 119–135.
  28. 28. Charles C, Rasson J-P. Wavelet denoising of Poisson-distributed data and applications. Comput. Stat. Data Anal. 2003; 43: 139–148.
  29. 29. Nowak RD, Baraniuk RG. Wavelet-domain filtering for photon imaging systems. IEEE Trans. Image Process. 1999; 8: 666–678. pmid:18267482
  30. 30. Lingenfelter DJ, Fessler JA, He Z. Sparsity regularization for image reconstruction with poisson data. Proc. SPIE, 2009; 7246: 72 460F–72 460F–10.
  31. 31. Harmany Z, Marcia R, Willett R. This is SPIRAL-TAP: Sparse poisson intensity reconstruction algorithms–theory and practice. IEEE Trans. Image Process. 2012; 21(3): 1084–1096. pmid:21926022
  32. 32. Salmon J, Harmany Z, Deledalle CA, Willett R. Poisson noise reduction with non-local PCA. J. Math. Imaging Vis. 2014; 48: 279–94.
  33. 33. Giryes R, Elad M. Sparsity Based Poisson Denoising with Dictionary Learning. IEEE Trans. Image Process. 2014; 23(12): 5057–5069. pmid:25312930
  34. 34. Boulanger J, Kervrann C, Bouthemy P, Elbau P, Sibarita JB, Salamero J. Patch-based nonlocal functional for denoising fluorescence microscopy image sequences. IEEE Trans. Med. Imag. 2010; 29(2): 442–454.
  35. 35. Stein CM. Estimation of the mean of a multivariate normal distribution. Annal. Stat. 1981; 1135–1151.
  36. 36. Starck J-L, Murtagh FD, Bijaoui A. Image processing and data analysis: the multiscale approach. Cambridge University Press. 1998.
  37. 37. Mäkitalo M, Foi A. Noise parameter mismatch in variance stabilization, with an application to Poisson-Gaussian noise estimation. IEEE Trans. Image Process. 2014; 23(12): 5348–5359. pmid:25343760
  38. 38. Starck J-L, Candès EJ, Donoho DL. The curvelet transform for image denoising. IEEE Trans. Image Process. 2002; 11: 670–684. pmid:18244665
  39. 39. Do MN, Vetterli M. The contourlet transform: an efficient directional multiresolution image representation. IEEE Trans. Image Process. 2005; 14: 2091–2106. pmid:16370462
  40. 40. Burt PJ, Adelson EH. The Laplacian pyramid as a compact image code. IEEE Trans Comm. 1983; 31: 532–540.
  41. 41. Po D-Y, Do MN. Directional multiscale modeling of images using the contourlet transform. IEEE Trans. Image Process. 2006; 15: 1610–1620. pmid:16764285
  42. 42. Bamberger RH, Smith MJ. A filter bank for the directional decomposition of images: Theory and design. IEEE Trans. Signal Process. 1992; 40: 882–893.
  43. 43. Mallat S. A wavelet tour of signal processing. Academic press. 1999.
  44. 44. Eslami R, Radha H. The contourlet transform for image denoising using cycle spinning. Asolomar Conf. Signals, Syst. Comput. 2003; 2: 1982–1986.
  45. 45. Coifman RR, Donoho DL. Translation-invariant de-noising: Springer, in Wavelets and Statistics. Springer Lecture Notes in Statistics. 1995; 103: 125–150.
  46. 46. Da Cunha AL, Zhou J, Do MN. The nonsubsampled contourlet transform: theory, design, and applications. IEEE Trans. Image Process. 2006; 15: 3089–3101. pmid:17022272