1.

Introduction

Noise analysis and noise suppression play an important role in imaging processing applications. This paper is motivated by the need of noise reduction in astronomical images. For this type of images, noise contamination represents a significant problem, since these data are usually acquired at very low levels of illumination. The image data used in this paper originate from an international (Czech-Spanish) experimental system called burst observer optical transient exploring system (BOOTES).¹ This system has been in service since 1998 as the first Spanish robotic telescope for sky observation. This system is one of six fully operational similar systems in the world. The main goal of this project is to observe extragalactic objects and detect a new optical transient (OT) of gamma-ray burst (GRB) sources.

Within the denoising algorithm, we decided to use the wavelet transform for video frame representation. A lot of papers that deal with noise suppression and modeling in the wavelet domain consider noise to be Gaussian, uncorrelated, and independent. For Gaussian noise, it is relatively simple to estimate its variance. One of the most widely used estimators is the median absolute deviation (MAD).² However, as demonstrated in our previous work on noise analysis of astronomical video of weak sporadic meteors,³ the contained noise is significantly non-Gaussian. As a consequence, its statistics are not preserved after transformation into the wavelet domain,⁴ and the noise model parameters have to be estimated for every level of the wavelet decomposition, which is computationally expensive. To overcome this drawback, we derive equations that allow for finding the $k$’th sample moment for every decomposition level by evaluating the sample moments only in the spatial domain.

The possible usage of these proposed equations is in the process of astronomical light images acquisition, where the dark current contaminates the acquired light image data. The dark current originates from the thermally generated charge in the crystalline lattice of the charge-coupled device (CCD) sensor (or another type of a photosensitive device) and is present in the camera output even when an astronomical camera is darkened. For very low temperatures of the CCD sensor (approximately 173 K), the generated charge is negligible.⁵ Hence, astronomical cameras are usually cooled utilizing the Peltier effect in order to partially suppress the thermal charge. The temperature dependency of the dark current is discussed in Refs. 5 and 6.

Dark current may also be eliminated by a specific configuration of the CCD chip. A number of configuration methods have been developed recently.^7,8 Okyay et al.⁹ study the effect of the varying electrode area asymmetry on the leakage behavior of metal-semiconductor-metal photodetectors (MSM-PDs). The authors demonstrate an effective method of dark current suppression by the means of the asymmetric electrode area design and the appropriate biasing scheme in MSM-PDs.

Another group of methods for dark current suppression is based on postprocessing of the acquired images. Dark current may be removed from the image by using dark frame subtraction. The dark frame represents a mapped thermally generated charge that is acquired by a darkened camera at the same conditions (i.e., exposure time and temperature) as the corresponding light image. An example of the dark frame acquired by the SBIG ST-8 camera¹⁰ is shown in Fig. 1. However, light images (see Fig. 2) are frequently acquired without the corresponding dark frames. In this case median filtering is usually used for dark current elimination. Unfortunately, this simple method often causes significant corruption of the imaged astronomical objects.

Fig. 1

Dark frame (inverted grayscale) acquired by SBIG ST-8 (conditions: $T=288.15\mathrm{K}$, $t_{\exp }=60\mathrm{s}$).

Fig. 2

Light image 2g980831.fits (inverted grayscale) without dark frame correction (conditions: $T=277.36\mathrm{K}$, $t_{\exp }=180\mathrm{s}$).

More sophisticated techniques for dark current suppression are based on statistical modeling of the dark frame. In the spatial domain, Baer¹¹ models the dark current histogram by using the Log-Normal, the Gamma, and the Inverse Gamma distribution. Goesele et al.¹² propose a dark frame subtraction method based on the entropy. This method is based on the assumption that the dark current signal increases the entropy of the resulting image, and thus the entropy is minimal when the dark current signal is not present at all.

For its sparsity and multi-scale nature, the wavelet transform is used by many researchers in various signal processing areas such as image denoising or compression. Lyu and Simoncelli¹³ model natural photographic images represented in a multi-scale basis using the Gaussian mixture model. The authors proposed algorithm for additive Gaussian noise suppression on this framework. In the field of geoscience, Amirmazlaghani¹⁴ introduced a Bayesian-based speckle-suppression method that employs the two-dimensional (2-D) generalized autoregressive conditional heteroscedasticity (2-D-GARCH) model of the wavelet coefficients obtained for log-transformed SAR images. Nevertheless, wavelet-based denoising plays an important role also in astronomy. In Ref. 15, Schmitt, Starck et al. study the large area telescope (LAT) representing the main instrument of the Fermi gamma-ray space telescope. Their two main scientific objectives—studying of the Milky Way diffuse background and detection of point sources—were complicated by the lack of photons. Hence, they proposed a powerful method of Poisson noise removal on the sphere that is efficient for low count Poisson data. This method uses the isotropic undecimated wavelet transform (IUWT) and the curvelet transform as spherical multiscale transforms. In general, the undecimated wavelet transform (UWT) gives outstanding results in denoising,¹⁶ and thus we decided to use this representation in this paper.

The paper is organized as follows. Section 2 discusses image acquisition using the SBIG ST-8 camera and modeling of the real dark frames. Section 3 is dedicated to acquisition and analysis of real data using MAIA system. This section also describes noise extraction and the whitening process applied to the extracted noise due to its undesirable behavior. At the end of this section, the UWT is described and the relation between the sample moments in the spatial and the wavelet domain is derived using the moment-generating function. Section 4 describes the used noise PDF model and the related moment equations. Section 5 contains the comparison of the two discussed methods for estimating the sample moments in the wavelet domain: the proposed method and the direct method that evaluates the moments directly from the noise wavelet decomposition. The main difference of these methods lies in computational cost that is experimentally evaluated. Section 6 is dedicated to performance of the proposed method in the denoising algorithm based on Bayesian shrinkage.

2.

Dark Frame Acquisition for Modeling the Temperature Dependency

Dark current (represented by the dark frame) is a highly undesirable component of astronomical images and video frames. It may be modeled as white impulsive noise. This section describes the statistical characteristics of dark frames in the spatial domain and also the temperature dependency of the sample moments. The obtained dependency allows us to evaluate the moments for an arbitrary temperature.

For analyzing the temperature dependency of dark current, we acquired a set of dark frames using the SBIG ST-8 camera with a sensor size of $510\times 710\text{pixels}$ and an exposure time of 60 s. By exploiting the temperature control functionality of the Peltier device, a set of 100 images is collected for each temperature within the range from 268.15 to 293.15 K with the step of 5 K. The 100 images acquired at each temperature were averaged in order to suppress dark current fluctuations.

The temperature dependency of the dark current may be modeled using the exponential model⁵ given as

The estimated model for the observed data is presented in Fig. 3. The exponential model is obtained from Eq. (1) and captures the temperature dependency of the first sample moment (i.e., the mean value).

Fig. 3

Model of the dark current dependency on temperature (the first sample moment) for the observed data ($A=6070$, $B=2.05·{10}^{-20}$).

For PDF model parameters estimation, it is necessary to evaluate the second and fourth sample moments of dark current at each measured temperature in the spatial domain. The temperature dependencies of these moments are modeled using the exponential models to capture the prior information about the CCD sensor in the spatial domain. The $l$’th sample central moment is given by

We found empirically that a satisfactory fit of the moment temperature dependencies $M_{l,i}(T)$ is achieved with the exponential model.¹⁷ Hence, using Eq. (1) we may write

Minimization of $S(A_l,B_l)$ is carried out using the partial derivatives with respect to the variables $A_l$ and $B_l$ given as

This leads to the following system of linear equations

The solution of the equation system $\mathit{z}={[z_1,z_2]}^T$, where $A_l=e^{z_1}$ and $B_l=z_2{k}_{\beta }$ is given as

The second and fourth sample moment evaluated for the mean dark frames at different temperatures are displayed in Fig. 4. The quality of the fit may be assessed by exploiting the sample kurtosis ${\kappa }_n(T)$ evaluated in the spatial domain as

Fig. 4

Dependency of the second and the fourth sample moments modeled by an exponential, (a) Dependency of the second moment, $A_2=3.03·{10}^{30}$, $B_2=1.17·{10}^{-19}$, (b) Dependency of the fourth moment, $A_4=2.39·{10}^{62}$, $B_4=1.10·{10}^{-19}$.

The computed values of ${\kappa }_n(T)$ are displayed in Fig. 5. Due to great oscillations, the polynomial approximation is not suitable. These oscillations also suggest that the polynomial models generate negative values of the sample kurtosis.

Fig. 5

Sample kurtosis in the spatial domain computed using the modeled moment dependencies.

3.

Acquisition and Analysis of Real Video Data

The video data used in this experiment come from the meteor automatic imager and analyzer (MAIA) project. A number of papers are dedicated to this system. The electro-optical¹⁸ characteristics of this system are covered by Ref. 19, its noise characteristics are discussed in Ref. 3, and the overall description of this system can be found in Ref. 20. The system for video acquisition consists of the Pentax SMC FA $1.4/50\mathrm{mm}$ fast autofocus lens, the Philips (Photonis) XX1332 2nd Generation image intensifier, and the JAI CM-040GE monochrome progressive scan camera with GigE Vision interface and the Pentax H1214-M $1.4/12\mathrm{mm}$ fast lens. The camera in the test setup was delivering an uncompressed video stream of 61.15 frames per second with the resolution of $776\times 582\text{pixels}$ and 256 grayscale levels.

We carry out the noise analysis similarly as for the testing data in Ref. 3. A real video sequence depicts the selected area of interest from the night sky. Approximately 100 frames do not contain any moving object and depict the same objects. And thus these frames represent multiple realizations of the same random process. It is possible to average these frames to obtain the averaged frame with suppressed noise [see Fig. 6(a)]. The average frame is then subtracted from each noisy frame to extract several noise realizations for our experiments. We assume noise to be stationary (of constant statistical properties) from frame to frame, at least for a certain number of consecutive video frames. In Ref. 3, we found that mainly for low illumination levels, noise from our system is signal-dependent.

Fig. 6

(a) Averaged video frame with inverted grayscale, (b) Amplitude spectrum of the extracted noise from the MAIA system.

However, the main problem revealed by the noise analysis is that the extracted noise component is correlated. This means that the noise samples are not independent and the noise spectrum is not flat [as displayed in Fig. 6(b)]. Hence, we need to implement the whitening process prior to the application of the proposed algorithm.

4.

Marginal Noise Model

For modeling of non-Gaussian noise in the wavelet domain, we exploit an exponential marginal model—the generalized Laplacian model (GLM). This model is commonly used for modeling of filtered images such as the result of the wavelet transform. The GLM is defined as

The parameters of this PDF may be estimated by using the system of moment equations.²⁷ For the sake of simplicity, let us consider the noise $n$ to have $\mu =0$. Then the second and the fourth moment run as

In accordance with Ref. 27, the parameters estimation task can be simplified by using the kurtosis statistic

The GLM model was chosen for our experiments because of its simplicity and flexibility to different types of noise. In the cases that the GLM cannot be used for the particular type of noise, another, more complex model such as Gaussian mixture model (GMM)²⁶ must be employed.

5.

Comparison of the Direct and the Proposed Method for Moments Estimation

This section is focused on the comparison of the direct method and the proposed method for sample moments evaluation in the wavelet domain. The direct method estimates the sample moments directly from the wavelet decomposition of noise, whereas the proposed method exploits the derived equations. The exact choice of the equations is driven by the length $k$ of the convolution filters. In our experiments, we exploit two wavelets: the Haar wavelet and the Daubechies 4 (db4) wavelet. The moments for the Haar 2-tap filters are computed by using Eqs. (21) and (22), whereas for the Daubechies 8-tap filters, we need to use the generalized formulations in Eqs. (23) and (24).

We consider that the analyzed noise has the same moment values for every wavelet subband at a given decomposition level. Therefore a single moment value may be evaluated at every decomposition level (from the approximation coefficients). This statement holds for white noise. If that noise is not white, the whitening process described in Sec. 3.1 may be applied.

The comparison of the two methods for different types of noise (both real and simulated) is presented as follows. The computed moments for uniformly distributed white noise are presented in Table 1 for salt-and-pepper white noise in Table 2 for Gaussian noise in Table 3 for real whitened noise in Table 4 and for real dark current in Table 5. The parameters of the simulated noise components were chosen with the aim to obtain similar values of the computed moments for both methods.

Table 1

The summary of the evaluated sample moments for uniform white noise U(−5,5) for the input moments M2(n)=8.4, M4(n)=126.

Table 2

The summary of the evaluated sample moments for salt-and-pepper noise distributed from 0 to 20 and with ε=0.01 for the input moments M2(n)=1.9, M4(n)=760.

Table 3

The summary of the evaluated sample moments for Gaussian white noise N(0,3) for the input moments M2(n)=9.0, M4(n)=244.

Table 4

The summary of the evaluated sample moments for whitened real noise extracted from our system for the input moments M2(n)=3.8, M4(n)=166.

Table 5

The summary of the evaluated sample moments for real dark current represented by dark frame captured at T=268.15  K and texp=2.5  s for the input moments M2(n)=1.5, M4(n)=2917.

5.1.

Estimated Marginal Models

Table 6 contains the estimated parameters of the GLM in the wavelet domain for different types of noise, with the exception of the generated salt-and-pepper noise that cannot be fitted by this model, This type of noise requires a more general model, such as the Gaussian mixture model. The goodness-of-fit is measured by the Jeffrey divergence (i.e., the symmetrical version of the Kullback-Leibler divergence), which is an empirical measure of similarity between distributions (PDFs) based on the relative entropy²⁸ given by

Eq. (30)

$${\mathcal{J}}_{\mathcal{D}}=\sum _{i=1}^I\left[p_i\ln \right(\frac{p_i}{0.5(p_i+H_i)})+H_i\ln (\frac{H_i}{0.5(p_i+H_i)}\left)\right],$$

where $I$ stands for the number of PDF samples, $p$ denotes the model PDF, and $H$ represents the histogram of the noise in the wavelet domain. This histogram is normalized and its bin width is optimized according to Freedman and Diaconis²⁹ written as

Eq. (31)

$$B_W=2(N_{0.75}-N_{0.25})I^{-\frac{1}{3}},$$

where the term in parentheses denotes the interquartile range, $N_{0.75}$ is the 75th percentile, and $N_{0.25}$ is the 25th percentile.

Table 6

The estimated parameters of the noise GLM in the wavelet domain (for the Haar wavelet and three decomposition levels).

$\text{Uniform}|\text{Proposed}$	$\gamma D_{\mathbf{3}}^{(v)}$	$\gamma D_{\mathbf{3}}^{(h)}$	$\gamma D_{\mathbf{3}}^{(d)}$	Direct	$\gamma D_{\mathbf{3}}^{(v)}$	$\gamma D_{\mathbf{3}}^{(h)}$	$\gamma D_{\mathbf{3}}^{(d)}$
$\nu$	2.02	2.02	2.02	$\nu$	2.02	1.99	2.01
$s$	4.11	4.11	4.11	$s$	4.11	4.06	4.11
${{\mathcal{J}}_{\mathcal{D}}}_1$	0.0027	0.0036	0.0028	${{\mathcal{J}}_{\mathcal{D}}}_2$	0.0026	0.0036	0.0027
${{\mathcal{J}}_{\mathcal{D}}}_{12}$	0.0001	0.0001	0.0000	${{\mathcal{J}}_{\mathcal{D}}}_{12}$	0.0001	0.0001	0.0000
$\text{Gaussian}|\text{Proposed}$	$\gamma D_{\mathbf{3}}^{(v)}$	$\gamma D_{\mathbf{3}}^{(h)}$	$\gamma D_{\mathbf{3}}^{(d)}$	Direct	$\gamma D_{\mathbf{3}}^{(v)}$	$\gamma D_{\mathbf{3}}^{(h)}$	$\gamma D_{\mathbf{3}}^{(d)}$
$\nu$	2.00	2.00	2.00	$\nu$	2.01	1.92	1.94
$s$	4.24	4.24	4.24	$s$	4.25	4.10	4.19
${{\mathcal{J}}_{\mathcal{D}}}_1$	0.0024	0.0030	0.0026	${{\mathcal{J}}_{\mathcal{D}}}_2$	0.0024	0.0022	0.0025
${{\mathcal{J}}_{\mathcal{D}}}_{12}$	0.0000	0.0008	0.0003	${{\mathcal{J}}_{\mathcal{D}}}_{12}$	0.0000	0.0008	0.0003
$\text{Real}|\text{Proposed}$	$\gamma D_{\mathbf{3}}^{(v)}$	$\gamma D_{\mathbf{3}}^{(h)}$	$\gamma D_{\mathbf{3}}^{(d)}$	Direct	$\gamma D_{\mathbf{3}}^{(v)}$	$\gamma D_{\mathbf{3}}^{(h)}$	$\gamma D_{\mathbf{3}}^{(d)}$
$\nu$	1.87	1.87	1.87	$\nu$	1.60	1.71	1.62
$s$	2.64	2.64	2.64	$s$	2.4	2.51	2.41
${{\mathcal{J}}_{\mathcal{D}}}_1$	0.0105	0.0063	0.0079	${{\mathcal{J}}_{\mathcal{D}}}_2$	0.0126	0.0074	0.0102
${{\mathcal{J}}_{\mathcal{D}}}_{12}$	0.0105	0.0034	0.0089	${{\mathcal{J}}_{\mathcal{D}}}_{12}$	0.0105	0.0034	0.0089

Table 6 contains the values of Jeffrey divergence computed as the similarity measure between the model calculated by the proposed method and the optimized histogram of noise wavelet coefficients (${{\mathcal{J}}_{\mathcal{D}1}}_{}$), between the model computed by the direct method and the noise histogram (${{\mathcal{J}}_{\mathcal{D}2}}_{}$), and between the models produced by the two moment estimation methods (${{\mathcal{J}}_{\mathcal{D}12}}_{}$).

The models from Table 6 are depicted in Fig. 7. Both the models and the optimized histograms are plotted in the logarithmic scale to better visualize the quality of fit.

Fig. 7

The estimated GLMs of diagonal details at decomposition level 3 $\gamma D_3^{(d)}$ for the direct and the proposed method modeling; (a), (b) uniform white noise, (c), (d) Gaussian white noise, and (e), (f) real whitened noise.

According to the central limit theorem (CLT), a weighted sum of Gaussian random variables again produces a Gaussian variable. If we apply the same process to a uniform variable, the conditions of the CLT may be satisfied. A lot of pseudorandom number generators are based on this principle. The satisfaction of the CLT for the filtered uniformly distributed variable is well visible in Table 1, where the kurtosis approaches 3 practically for all decomposition levels. However, the decomposed astronomical data (dark frames) are strongly non-Gaussian. The kurtosis is usually considerably greater than 3 for all decomposition level. Table 5 contains the second and fourth sample moments of the decomposed dark frame and illustrates its strong non-Gaussianity. Figure 8 contains the typical histogram with high value of kurtosis of the decomposed dark frame.

Fig. 8

GLM of wavelet band along with histogram (decomposed dark frame), fourth decomposition level, vertical details, $\kappa =201$, $\nu =0.29$.

6.

Testing of the Derived Equations Using Bayesian Shrinkage

We have demonstrated that the proposed method surpasses the direct method in terms of computation cost. In this section, we experimentally verify that the proposed method does not negatively impact the denoising results. For its suitability, we selected the Bayesian shrinkage algorithm based on the minimum mean square error (MMSE) estimator.²⁷ We applied this noise removal algorithm to multimedia data that were artificially contaminated with a high amount of uniform white noise $\mathcal{U}(-20,20)$ and also to real astronomical light images contaminated with dark current.

6.1.

Multimedia Data

As discussed above, the parameters of non-Gaussian noise are not preserved by convolution filtering. From the various types of non-Gaussian noise that we generated, uniform noise was selected because of the perfect fit of the wavelet-based GLM. We consider image $x$ contaminated by additive noise as defined in Eq. (13). Hence, the additive model in the wavelet domain is given as $Y=X+N$, where $Y=\mathcal{U}\mathcal{W}\mathcal{T}\{y\}$ denotes the noisy observation (i.e., the acquired data), $N=\mathcal{U}\mathcal{W}\mathcal{T}\{n\}$ presents noise, and $X=\mathcal{U}\mathcal{W}\mathcal{T}\{x\}$ is the ideal noise-free image. The conditional mean of the posterior PDF $p_{X|Y}(x|y)$ provides the least square estimation of $X$. The formula for the MMSE¹⁶ estimator runs as

Eq. (32)

$$\begin{gathered} \widehat{X}(Y)={\int }_{-\infty }^{+\infty }{p}_{X|Y}(x|y)x\mathrm{d}x=\frac{{\int }_{-\infty }^{+\infty }{p}_{Y|X}(y|x)p_X(x)x\mathrm{d}x}{{\int }_{-\infty }^{+\infty }{p}_{Y|X}(y|x)p_X(x)\mathrm{d}x}, \\ =\frac{{\int }_{-\infty }^{+\infty }{p}_N(y-x)p_X(x)x\mathrm{d}x}{{\int }_{-\infty }^{+\infty }{p}_N(y-x)p_X(x)\mathrm{d}x}, \end{gathered}$$

where $p_{Y|X}(y|x)$ denotes the likelihood function, $p_X(x)$ represents the a priori model, and $p_N(x)$ stands for the noise model. Both these random variables are modeled by the wavelet-based GLM. Similarly to Ref. 30, we define the theoretical central moments of $Y$. The second moment is given by

Eq. (33)

$$m_2(Y)=\frac{{\tau }^2\mathrm{\Gamma }\left(\frac{3}{\lambda }\right)}{\mathrm{\Gamma }\left(\frac{1}{\lambda }\right)}+\frac{s^2\mathrm{\Gamma }\left(\frac{3}{\nu }\right)}{\mathrm{\Gamma }\left(\frac{1}{\nu }\right)}=m_2(X)+m_2(N),$$

where $\tau$ and $\lambda$ are the GLM parameters of the signal. The fourth moment of $Y$ runs as

Eq. (34)

$$\begin{gathered} m_4(Y)=\frac{{\tau }^4\mathrm{\Gamma }\left(\frac{5}{\lambda }\right)}{\mathrm{\Gamma }\left(\frac{1}{\lambda }\right)}+\frac{6s^2{\tau }^2\mathrm{\Gamma }\left(\frac{3}{\nu }\right)\mathrm{\Gamma }\left(\frac{3}{\lambda }\right)}{\mathrm{\Gamma }\left(\frac{1}{\nu }\right)\mathrm{\Gamma }\left(\frac{1}{\lambda }\right)}+\frac{s^4\mathrm{\Gamma }\left(\frac{5}{\nu }\right)}{\mathrm{\Gamma }\left(\frac{1}{\nu }\right)} \\ =m_4(X)+6m_2(N)m_2(X)+m_4(N). \end{gathered}$$

The GLM parameters of the noise-free signal are estimated from second and fourth moment of the observed signal and noise using the kurtosis formula

Eq. (35)

$$\begin{gathered} {\kappa }_X=\frac{\mathrm{\Gamma }\left(\frac{5}{\lambda }\right)\mathrm{\Gamma }\left(\frac{1}{\lambda }\right)}{{\mathrm{\Gamma }}^2\left(\frac{3}{\lambda }\right)} \\ =\frac{m_4(Y)-m_4(N)-6m_2(N)[m_2(Y)-m_2(N)]}{{(m_2(Y)-m_2(N))}^2}. \end{gathered}$$

From Eq. (33), we may derive

Eq. (36)

$$\tau =\sqrt{[m_2(Y)-m_2(N)]\frac{\mathrm{\Gamma }\left(\frac{1}{\lambda }\right)}{\mathrm{\Gamma }\left(\frac{3}{\lambda }\right)}}.$$

The values of the moments used in the above equations are estimated from the data using the sample moments. The $k$’th central sample moment of $X$ is given by $M_k(X)=\frac{1}{I}\sum _{i=1}^I{[X_i-E(X)]}^k$.

Figure 9 depicts two testing images (Brada and House) from our image database. These images (8 bpp and uncompressed tiff format) were contaminated with additive uniform white noise and subsequently denoised by the Bayesian estimator using both the direct method and the proposed method for estimation of the GLM parameters. Table 8 contains the root mean square error (RMSE) measure computed as

Eq. (37)

$${\mathrm{RMSE}}_{\text{out}}=\sqrt{\frac{1}{I}\sum _{i=1}^I{(x-\widehat{x})}^2},$$

where $x$ denotes the original image before contamination and $\widehat{x}$ is the reconstructed denoised image. However, in most real-world applications, the precontamination image is not available and the ${\mathrm{RMSE}}_{\text{out}}$ and ${\mathrm{RMSE}}_{\text{in}}$ measures must be estimated in another way. ${\mathrm{RMSE}}_{\text{in}}$ is directly related to the noise second moment $M_2^{\prime }(n)=M_2(n)+E^2(n)$, ${\widehat{\mathrm{RMSE}}}_{\text{in}}=\sqrt{M_2^{\prime }(n)}$. In Sec. 2, we exploit the model of the temperature dependency of dark current for the SBIG camera. The modeled second moment $M_2(n)$, the mean value $E(n)$, and the fourth moment $M_4(n)$ represent a priori information about the noise component contained in the image.

Fig. 9

Testing images ($256\times 256\text{pixels}$) denoised by Bayesian shrinkage using the Db4 wavelet and 4 UWT decomposition levels (a) The noise House image [added uniform noise $\mathcal{U}(-20,20)$], (b) the denoised House image by using the direct method (c) the denoised House image using the proposed method, (d) the noisy Brada image [added noise $\mathcal{U}(-20,20)$], (e) the denoised Brada image using the direct method, and (f) the denoised Brada image using the proposed method.

Table 8

The computed RMSE measure for denoising via Bayesian shrinkage using the Db4 wavelet and four UWT decomposition levels for the multimedia images.

Image	${\mathrm{RMSE}}_{\text{in}}$	${\mathrm{RMSE}}_{\text{out}}|\text{proposed}$	${\mathrm{RMSE}}_{\text{out}}|\text{direct}$
House	11.55	6.13	6.19
Brada	11.55	4.06	4.10

Since we consider $y=x+n$, the second sample moment of $y$ (a noncentral moment in this case) is given as

Eq. (38)

$$m_2(y)={\int }_{-\infty }^{\infty }{y}^2[{\int }_{-\infty }^{\infty }{p}_n(n)p_x(y-n)\mathrm{d}n]\mathrm{d}y=m_2(x)+2E(x)E(n)+m_2(n).$$

Using Eq. (38) we may write

Eq. (39)

$$M_2^{\prime }(y)=M_2^{\prime }(x)+2E(x)E(n)+M_2^{\prime }(n),$$

Eq. (40)

$$E(y)=E(x)+E(n).$$

If we consider the reconstructed light image to be given as $\widehat{x}(y)=x+n_{\mathrm{res}}$, where the noise is suppressed to $n_{\mathrm{res}}$, then ${\mathrm{RMSE}}_{\text{out}}$ can be estimated by

Eq. (41)

$${\widehat{\mathrm{RMSE}}}_{\text{out}}=\sqrt{M_2^{\prime }(n_{\mathrm{res}})}=\sqrt{M_2^{\prime }[\widehat{x}(y)]-M_2^{\prime }(x)-2E(x)E(n_{\mathrm{res}})},$$

where $E(n_{\mathrm{res}})=E(n)$ denotes the mean value of the residual noise. Furthermore, if the modeled temperature dependency of the sample moments belongs to dark frames acquired at an arbitrary exposure time, we can normalize it to any selected exposure time. We consider a dark frame with the exposure time $t_{\exp }$. Then the $k$’th moment $M_k^{t_{\exp }}(n)$ normalized to required exposure $t_{\mathrm{req}}$ runs as

Eq. (42)

$$M_k^{t_{\mathrm{req}}}(n)=\frac{1}{I}\sum _{i=1}^I{\left(\frac{t_{\mathrm{req}}{n}_i}{t_{\exp }}\right)}^k={\left(\frac{t_{\mathrm{req}}}{t_{\exp }}\right)}^k{M}_k^{t_{\exp }}(n).$$

The RMSE measure in Table 8 illustrates that the proposed algorithm gives practically the same results as the direct method for multimedia images.

6.2.

Astronomical Data

We applied the Bayesian shrinkage also to an astronomical light image (containing the object such as stars, nebulae, etc.) contaminated by dark current. We exploited the following process:

• 1. Read the temperature and the exposure time of the acquired light image in the image header.

• 2. Find the appropriate values of sample moments from the approximated temperature dependencies (see Sec. 2). If necessary, normalize the value of the sample moments in accordance with the exposure time value.

• 3. Convert the sample moments into the wavelet domain using the derived equations.

• 4. Estimate the parameters of the wavelet-based PDF model.

• 5. Transform the light image contaminated by dark current into the wavelet domain using the UWT.

• 6. Apply the Bayesian estimator to the detail wavelet subbands at all decomposition levels.

• 7. Apply the inverse UWT to the altered coefficients.

The astronomical light image of the Messier 42 nebula and the same image corrected by the MMSE estimator is depicted in Figs. 10 and 11, respectively. The estimations of the RMSE measure for light astronomical data 3m42-d03.sbg.dat (five UWT decomposition levels, Db4, wavelet) are: ${\widehat{\mathrm{RMSE}}}_{\text{in}}=1498$, proposed: ${\widehat{\mathrm{RMSE}}}_{\text{out}}=325$, direct: ${\widehat{\mathrm{RMSE}}}_{\text{out}}=333$. Hence, the direct and the proposed algorithm produce practically the same results in the sense of RMSE. In both cases, dark current is considerably suppressed and the objects become visible and detectable.

Fig. 10

The 3m42-d03.sbg.dat light image ($512\times 1024\text{pixels}$; inverted grayscale) without dark frame correction ($T=277.36\mathrm{K}$, $t_{\exp }=1000\mathrm{s}$).

Fig. 11

The 3m42-d03.sbg.dat light image (inverted grayscale) denoised using the MMSE estimator.

7.

Conclusion

The proposed algorithm based on the moment-generating function can be used for evaluation of the noise model parameters for noise removal in the wavelet domain. Generally, the proposed method is capable of finding all sample moments in the wavelet domain by converting the sample moments originally evaluated in the spatial domain. This is possible owing to the fact that the second sample moment is preserved in accordance with Parseval’s theorem.²⁶

Table 6 illustrates behavior of the proposed noise modeling algorithm. For Gaussian white noise, the model parameters are preserved in the wavelet domain. On the case of uniform white noise, we may observe the effect of the CLT on a weighted summation of independent random variables. (The value of shape parameter $\nu$ is approaching 2.)

The GLM was chosen for modeling of several types of noise in the wavelet domain. This model is widely used for filtered images modeling, and moreover, the parameters estimation is simple when we use the method of moments. As a quality criterion for model parameters estimation, the Jeffrey divergence was chosen. The Jeffrey divergence was evaluated for the optimized histogram of the wavelet coefficients and the estimated models computed using both the proposed and the direct method. For further methods comparison, this measure was also computed for the models produced by these two methods. The results indicate that if noise is not white (i.e., is correlated), the whitening process can be used.

The proposed method and the direct method were compared with respect to computational efficiency (see Table 7). The direct method exploits two separate undecimated wavelet decompositions: the UWT of the image data and the UWT of the extracted noise, whereas the proposed method uses the UWT only for the image data and estimates the noise model parameters in the spatial domain. The difference in computation cost between these two methods increases with the size of the analyzed image. For an image of the size ($1024\times 1024\text{pixels}$), the proposed method runs approximately seven times faster than the direct method.

Performance of the proposed algorithm was also tested using the Bayesian shrinkage. This denoising method exploits the noise and the image PDF model for image estimation. The proposed algorithm was tested for multimedia images contaminated with uniform white noise and a real astronomical image. The results show that the direct and the proposed method cause practically the same denoising results.