1. Introduction
In digital imaging systems, images are deteriorated by random noise coming from the complementary metal-oxide semiconductor/charge-coupled device (CMOS/CCD) image sensor [
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12,
13,
14,
15,
16]. Compared with the additive signal-independent noise model, the signal-dependent noise model is more accurate at characterizing random noise of the CMOS/CCD image sensor [
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12]. Most researchers assumed that the signal-dependent noise model of the digital imaging sensor is cond as a Poisson-Gaussian noise model and the validity of the Poisson–Gaussian noise model was certified by CMOS sensors from Nokia camera phones, CCD sensors from Fujifilm cameras and CMOS sensors from Canon cameras [
1,
2,
3,
4,
5,
6,
7,
8,
9,
10,
11,
12]. The Poisson-Gaussian noise model is composed of a signal-dependent term accounting for photon noise (Poisson) and a signal-independent term accounting for the remaining noise in the readout data (Gaussian) [
1], as shown in (1).
In (1),
m and
n are the horizontal and vertical coordinates of the image pixel, respectively,
y(
m,
n) is the noisy pixel value and
x(
m,
n) is the noiseless value. Further,
is the zero-mean signal-independent Gaussian noise,
is the signal-dependent Poisson noise and the mean and variance of
are equal to
. Moreover,
a is the variance of noise that is associated with analog gain of the CMOS/CCD image sensor and
b is the variance of noise that is related to readout noise of the CMOS/CCD image sensor [
1]. In this study, noise parameters,
a and
b, are estimated.
Because the Gaussian noise and Poisson noise are uncorrelated, (2) can be deduced from (1).
In (2), is the variance of noise .
Parameter estimation is a vital step in achieving a noiseless image. Many algorithms for estimating parameters of the Poisson–Gaussian noise model were presented [
1,
2,
5,
6,
7]. Owing to the linear relationship between the pixel intensity and the variance of noise, as verified in (2) and the intensity similarity of pixels in the homogeneous region, the main focus of most estimation algorithms is to acquire the parameters of the signal-dependent noise by finding the homogeneous regions [
1,
5,
6]. In [
1], analysis and smoothing in wavelets were used to ascertain the local estimation pairs of the pixel intensity and noise variance. Then, the maximum-likelihood approach was employed to fit the global noise model function. In [
5], image patches were classified based on their intensity and variance for finding the homogeneous regions that best represent the noise. Then, the weights of the cluster of connected patches were calculated based on the degree of similarity to the noise model. In [
6], the true parameter values of Poisson-Gaussian noise were estimated by searching for intersections of the unitary variance contours. Furthermore, methods in other related works did not strive to find the homogeneous regions to address this problem. For example, expectation-maximization was employed to estimate the parameters of Poisson-Gaussian noise in [
2]. In [
7], the generic Poisson-Gaussian noise model was simplified to a Gaussian-Gaussian noise model and the least squares method was used to estimate the noise model parameters.
The general idea of the existing estimation algorithms based on finding homogeneous regions is to first calculate the noise variance and noiseless intensity of every homogeneous region; then, the pairs of noise variances and noiseless intensities obtained from all homogeneous regions are used to fit the linear or piecewise linear curve; finally, the noise parameters from the fitting curve are acquired.
With consideration of effectiveness to detect and denoising of images, the proposed estimation method is also based on finding homogeneous regions to determine the noise parameters. However, different from the existing estimation algorithms, the proposed algorithm estimates the noise parameters by deriving the mean and variance of the mixed Poisson noise samples that are composed of the Poisson noise of all homogeneous regions in every block image and by building the mapping function among noise parameters to be estimated—variance of Poisson-Gaussian noise and that of Gaussian noise corresponding to the stitched region in every block image.
In this study, the input image is divided into 16 blocks. Then, all homogeneous regions in every block image are detected and denoised in the wavelet domain. Next, all the denoised homogeneous regions are pieced together to form a new stitched image in every block image and histogram analysis is used to obtain the intensities and the corresponding number of every intensity in the stitched noiseless image of every block image. Thus, the mixed Poisson noise samples corresponding to the stitched image in every block image can be obtained according to the definition of the signal-dependent Poisson noise in (1). Next, the mean and variance of the mixed Poisson noise samples in every block image are deduced. The mapping function among noise parameters to be estimated—variance of Poisson-Gaussian noise and that of Gaussian noise corresponding to the stitched region of every block image—is constructed accordingly. Finally, the unbiased estimations of noise parameters are obtained from the mapping functions of 16 block images.
The remainder of the paper is organized as follows.
Section 2 presents the proposed algorithm. The experimental results and performance comparison with other state-of-the-art parameter estimation approaches are reported in
Section 3.
Section 4 concludes our study.
2. Proposed Method
The flowchart of the proposed method is shown in
Figure 1. It is comprised of the following eight main steps.
Step ①: Block the input noisy image. With consideration of computational efficiency, the input noisy image is divided into 16 large blocks.
Step ②: Detect the homogeneous regions of every block in the wavelet domain. Owing to the good performance of detecting the homogeneous region in the wavelet domain in [
1], each block image is decomposited into four sub-band images
LL,
HL,
LH and
HH, as shown in
Figure 2.
The LL sub-band image consists of wavelet coefficients that are obtained by performing low-pass wavelet filtering in the row and column of the block image; the HL sub-band image consists of wavelet coefficients that are obtained by performing high-pass wavelet filtering in the row and low-pass wavelet filtering in the column of the block image; the LH sub-band image consists of wavelet coefficients that are obtained by performing low-pass wavelet filtering in the row and high-pass wavelet filtering in the column of the block image; and the HH sub-band image consists of wavelet coefficients that are obtained by performing high-pass wavelet filtering in the row and column of the block image. With respect to the orthogonality and regularity, the Daubechies wavelet basis (db6) is employed in this study. With consideration of complexity and accuracy of computation, the standard deviation of the 5 × 5 slipping window in the LL sub-band is calculated and compared with the threshold shown in (3) to ascertain the homogeneous region.
In (3), W{LL} is the wavelet coefficient of the 5 × 5 slipping window in the LL sub-band, is the mean of the 5 × 5 slipping window in the LL sub-band, denotes the standard deviation of the 5 × 5 slipping window in the LL sub-band and δ represents the threshold of . If the standard deviation of the 5 × 5 slipping window is less than δ, this region is homogeneous and the coordinates of the central wavelet coefficient in the 5 × 5 slipping window are recorded.
Furthermore, the pixel intensities in the homogeneous region are close and the maximal intensity difference of two pixels in the homogeneous region can be set to 15, according to the literature [
17]. The setting of threshold
δ of
in (3) is explained in
Figure 3. The base value of the homogeneous region is the minimal value of this homogeneous region. Because the
LL sub-band image is obtained by performing low-pass filtering between adjacent pixels in the row and column, the range of
LL wavelet coefficient corresponding to the homogeneous region is from base value corresponding to this homogenous region to base value +7.5. That is, the value range of all the
LL wavelet coefficients in the 5 × 5 window is from the base value of this homogeneous region to base value +7.5. Through calculation, the value of
in the 5 × 5 window will reach its maximum when there are 13
LL wavelet coefficients equal to base value +7.5 and the remaining 12
LL wavelet coefficients equal to base value. As a result,
δ in (3) is set to the maximum of
, that is,
δ = 3.75.
Step ③ Combine all homogeneous 5 × 5 windows to form the stitched sub-band image. According to the coordinates of the central wavelet coefficient in the 5 × 5 slipping window, all homogeneous windows in the
LL sub-band of each block are extracted and stitched to form a new stitched image. The same operation will be performed for
LH,
HL and
HH sub-bands and the stitched images in
LH,
HL and HH sub-bands can be obtained accordingly. The extraction and combination process is depicted in
Figure 4. The stitched wavelet coefficients in
LL,
LH,
HL and
HH wavelet sub-band can be reconstructed as a stitched image in every block, denoted as
ySi,
i = 1~16.
Step ④ Parameter estimation of Gaussian noise of the stitched image in every block. The median absolute deviation (MAD) is used to estimate the standard deviation of Gaussian noise of the stitched image
ySi [
7] as shown in (4).
In (4),
(
i = 1~16) is the estimated standard deviation of the Gaussian noise of the stitched image
ySi in every block,
Median(●) is the MAD value and
is the
HH sub-band wavelet coefficient of the stitched image in every block shown in
Figure 4. It can be determined from (1) that
,
i = 1~16.
Step ⑤ Denoise the stitched sub-band image of every block. In order to obtain the denoised stitched sub-band image, the mean values of the 5 × 5 slipping window in LH, HL and HH sub-bands—which are located in the corresponding position as the 5 × 5 slipping window in the LL sub-band—are calculated as (5). Then, all the wavelet coefficients of the 5 × 5 slipping window in LH, HL and HH sub-bands are replaced by these three mean values.
In (5), are the mean values of the 5 × 5 slipping window in the LH, HL and HH sub-bands and are the wavelet coefficients of the 5 × 5 slipping window in LH, HL and HH sub-bands of the noisy image.
After all the homogeneous regions are denoised, the block image in the wavelet domain is translated into a temporal image, denoted as Si, i = 1~16.
Step ⑥ Calculate the variance and mean of the mixed Poisson noise samples. According to (1), the Poisson-Gaussian noise corresponding to the noiseless stitched region
Si is denoted as
(
i = 1~16); the mixed Poisson noise corresponding to the noiseless stitched region
Si is denoted as
(
i = 1~16); and the Gaussian noise corresponding to the noiseless stitched region
Si is denoted as
(
i = 1~16). The forming process of the mixed Poisson noise
(
i = 1~16) corresponding to the noiseless stitched image
Si is shown in
Figure 5.
It can be seen from
Figure 5 that the histogram analysis offers the pixel value
x (
x = 0~255) and the corresponding pixel number (denoted as
ni(x) (
i = 1~16,
x = 0~255)) of every noiseless stitched image (denoted as
Si (
i = 1~16)). Because the parasitic Poisson noise is signal-dependent and
, as shown in (1), the Poisson noise samples
(
i = 1~16,
x = 0~255) and the corresponding sample size
ni(x) (
i = 1~16,
x = 0~255) can be obtained according to histogram analysis. All the noise samples
constitute the mixed Poisson noise sample
(
i = 1~16) corresponding to
Si (
i = 1~16) and the sample size of
can be represented as
Ni (
i = 1~16). As a result,
(
i = 1~16) and the unbiased estimations of the mean of
(
i = 1~16) can be obtained as (6).
In (6), is the unbiased estimation of the mean of the mixed Poisson noise sample and is the unbiased estimation of mean of Poisson noise samples .
The unbiased estimations of the variance of (i = 1~16) can be obtained as (7).
Substituting with (6), the unbiased estimations of the variance of can be obtained as (8).
In (8), is the unbiased estimation of variance of the mixed Poisson noise sample and is the unbiased estimation of the variance of Poisson noise samples .
Step ⑦ Build the mapping function among the noise parameters to be estimated, variance of and variance of . can be obtained by calculating the difference between the noisy stitched image ySi and the noiseless stitched image Si, i = 1~16. The unbiased estimation of variance of (i = 1~16) can be calculated by (9).
In (9), is the unbiased estimation of the mean of and denotes the unbiased estimation of the variance of .
As the Poisson noise and Gaussian noise are irrelevant, the variance of is the sum of the variance of and the variance of , as shown in (10).
From (7) to (10), let and . Thus, (11) can be obtained.
As a result, the estimation of parameter ai in every block can be acquired from (12).
In (12), the values of ai (i = 1~16) are positive.
Step ⑧ Perform unbiased parameter estimation of 16 block images. From (7) and (12), the parameters a and b in each block image are acquired. To improve the estimation accuracy, the unbiased estimations of a and b in 16 blocks are calculated as (13).