Translator Disclaimer
22 June 2012 Estimation of non-Gaussian noise parameters in the wavelet domain using the moment-generating function
Author Affiliations +
We discuss methods for modeling and removal of noise in astronomical images. For its favorable properties, we exploit the undecimated wavelet representation and apply noise suppression in this domain. Usually, the noise analysis of the studied imaging system is carried out in the spatial domain. However, noise in astronomical data is non-Gaussian, and thus the noise model parameters need to be estimated directly in the wavelet domain.We derive equations for estimating the sample moments for non-Gaussian noise in the wavelet domain. We consider that the sample moments in the spatial domain are known from the noise analysis and that the model parameters are estimated by using the method of moments.



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).1 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).2 However, as demonstrated in our previous work on noise analysis of astronomical video of weak sporadic meteors,3 the contained noise is significantly non-Gaussian. As a consequence, its statistics are not preserved after transformation into the wavelet domain,4 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.5 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.9 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 camera10 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.15K, texp=60s).


Fig. 2

Light image 2g980831.fits (inverted grayscale) without dark frame correction (conditions: T=277.36K, texp=180s).


More sophisticated techniques for dark current suppression are based on statistical modeling of the dark frame. In the spatial domain, Baer11 models the dark current histogram by using the Log-Normal, the Gamma, and the Inverse Gamma distribution. Goesele et al.12 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 Simoncelli13 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, Amirmazlaghani14 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,16 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.


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×710pixels 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 model5 given as

Eq. (1)

where A and B are the sensor material constants, T denotes the temperature in Kelvin, and kβ stands for the Boltzmann constant. This empirical model is valid only over a constrained temperature range. Its parameters can be estimated using linear regression assuming the following model:

Eq. (2)


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·1020).


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

Eq. (3)

where n denotes the pixel values of the acquired dark frames, n¯ is the mean value, and I represents the number of pixels.

We found empirically that a satisfactory fit of the moment temperature dependencies Ml,i(T) is achieved with the exponential model.17 Hence, using Eq. (1) we may write

Eq. (4)

where the coefficients Al and Bl are estimated via minimizing the square error given as

Eq. (5)


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

Eq. (6)


This leads to the following system of linear equations

Eq. (7)

where the matrix C is given by

Eq. (8)

and the column vector b is given by

Eq. (9)


The solution of the equation system z=[z1,z2]T, where Al=ez1 and Bl=z2kβ is given as

Eq. (10)


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 κn(T) evaluated in the spatial domain as

Eq. (11)


Fig. 4

Dependency of the second and the fourth sample moments modeled by an exponential, (a) Dependency of the second moment, A2=3.03·1030, B2=1.17·1019, (b) Dependency of the fourth moment, A4=2.39·1062, B4=1.10·1019.


The computed values of κ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.



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-optical18 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/50mm 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/12mm 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×582pixels 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.


Noise Whitening

For data whitening, we exploit a linear convolution filter that equalizes the nonflat noise spectrum. We derived this filter via inversion of the noise spectrum and thresholding of the spectral samples that are zero or near to zero. Let |S(u,v)| be the modulus of the noise spectrum, where u and v are the spatial frequency in the x and y directions, respectively. The frequency response of the inverse filter Hinv(u,v) is then given by

Eq. (12)

where δ denotes the threshold value.

The noise whitening procedure is applied to the noisy observations y in the frequency domain. For the sake of simplicity, we consider an additive noise model of y given by

Eq. (13)

where x denotes the useful signal and n represents noise introduced into the video data during the acquisition process. The whitening process W{} is implemented via linear filtering, and thus W{y}=W{x}+W{n}. Noise whitening enables us to proceed with the denoising algorithm.


Sample Moments in the Spatial Domain

In this paper, the noise model parameters are estimated via the method of moments. Hence, every noise realization n in the spatial domain is described by the sample moments. The r’th central sample moment is given by

Eq. (14)

where n¯ is the mean value. The sample moments of the noise are computed for every acquired video sequence. The choice of moments depends on the used model.


Undecimated Wavelet Transform

For image denoising, the UWT or the stationary wavelet transform (SWT) is a better choice than the critically sampled discrete wavelet transform (DWT).16 The main reason is that the DWT is shift variant,21 which limits its denoising performance.2 Wavelet shrinkage methods performed on the DWT coefficients usually cause unwanted artifacts around the objects, such as stars.22

For the reasons explained above, we choose the UWT23 for the video frame representation. The UWT is computed using a so-called à trous algorithm24 that produces the same number of wavelet coefficients at each scale (decomposition level). We use the following notation for the respective wavelet subbands: γDξ(v), where ξ in the subscript denotes the decomposition level and the superscript in the parentheses denotes the particular detail subband (v-vertical, h-horizontal, or d-diagonal).


Sample Moments in the Wavelet Domain

As we mentioned above, the noise present in astronomical images is non-Gaussian and thus it is necessary to evaluate the sample moments in the wavelet domain. This may be achieved by using the moment-generating function. This function of the random variable n (representing the analyzed noise) is closely related to the characteristic function4 and is defined by

Eq. (15)


The series expansion of eun suggests that the moment-generating function allows to find all moments of a given distribution.25 Provided that the random variable n has a continuous PDF, the Mn(u) is given by

Eq. (16)

where mk is the k’th moment.

For producing the relations between the moments in the spatial and the wavelet domain, we need to describe the wavelet transform process first. The one-dimensional (1-D) UWT corresponds to convolution filtering of n with the kernel h=[h1,h2hk] while the down-sampling step is omitted. Hence, each wavelet coefficient is computed as the weighted sum of the independent random variables n1,n2nk (noise pixels) given as

Eq. (17)

The moment-generating function MSk(u) of Sk then runs as

Eq. (18)


Eq. (19)


We are going to demonstrate that it is possible to find the sample moments in the wavelet domain by using the values of the sample moments from the spatial domain. As mentioned above, the moment-generating function is closely related with the moments of the distribution. Therefore, the r’th moment may be evaluated using the moment-generating function computed as the r’th derivative with respect to the variable u at u=0 given as

Eq. (20)

Let us consider a zero-mean noise n and its wavelet domain representation N=UWT{n}. In our method, we exploit the second and the fourth moment for noise description. These moments are related via the sample kurtosis as demonstrated in (29). For the sake of simplicity, we assume a short filter such as the Haar filter with the kernel h=[h1,h2]. The previously stated assumptions suggest the following moment relations. The second sample moment M2(N) in the wavelet domain is computed from M2(n) given by

Eq. (21)

Similarly the fourth moment M4(N) computed from M4(n) given by

Eq. (22)

Equations (21) and (22) may be generalized for filters with k coefficients h=[h1,h2hk] given as

Eq. (23)


Eq. (24)


The above equations are demonstrated on the case of the 1-D UWT. Nevertheless, the UWT may be easily extended also to the two-dimensional space. This transform is separable, and thus we carry out the convolution in the row direction and then also in the column direction to obtain the 2-D UWT decomposition. The sample moments of the resulting coefficients are evaluated using the derived equations.


Simplification of the derived equation for the fourth moment

This proposed equation for the fourth moment in the wavelet domain is a little complex, especially for longer wavelet transform filters. We found experimentally that Eq. (24) could be simplified for certain types of impulsive noise. If we consider the salt-and-pepper noise for 8 bpp (bits per pixel) images then the probability of the pixel value flipping to 0 is P(y=0)=ε/2 and the probability of flipping to 255 is P(y=255)=ε/2. If the parameter ε approximately satisfies ε0.05 then

Eq. (25)

As a result, Eq. (24) can be simplified as

Eq. (26)


We tested these findings also on the 16-bpp dark frames26 that were acquired by the SBIG ST-8 astronomical camera as described in Sec. 2. Our experiments indicate that these dark frames may be modeled as white impulsive noise for which the hot pixels do not always reach the maximum value of the dynamic range. Using this model, Eq. (25) is satisfied for all the acquired dark frames within the whole temperature range (from 268.15 to 293.15 K).


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

Eq. (27)

where μ denotes the mean value, the parameter ν presents a generalization in the sense of the shape of the heavy-tailed PDF, and the parameter s controls the width of the PDF. The function Z(s,ν)=2sνΓ(1ν), where Γ(x)=0tx1etdt, normalizes the exponential to a unit area.

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

Eq. (28)


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

Eq. (29)


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)26 must be employed.


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.

LevelWavelet: HaarWavelet: Db4

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.

LevelWavelet: HaarWavelet: Db4

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.

LevelWavelet: HaarWavelet: Db4

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.

LevelWavelet: HaarWavelet: Db4

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.

LevelWavelet: HaarWavelet: Db4


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 entropy28 given by

Eq. (30)

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 Diaconis29 written as

Eq. (31)

where the term in parentheses denotes the interquartile range, N0.75 is the 75th percentile, and N0.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).


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 (JD1), between the model computed by the direct method and the noise histogram (JD2), and between the models produced by the two moment estimation methods (JD12).

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 γD3(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, κ=201, ν=0.29.



Computational Cost

The major advantage of the proposed method is saving in the computational cost. To evaluate the sample moments in the wavelet domain, the direct method requires the UWT to be applied to the extracted noise. The sample moments are then evaluated in the chosen wavelet subband at every decomposition level. By contrast, the proposed method does not require the UWT of the extracted noise that is especially significant for the redundant wavelet representation that we are using. The sample moments of the noise are evaluated in the spatial domain and then converted to the wavelet domain exploiting the derived equations.

To demonstrate this point, we measured the computational time of the two methods under discussion. We used three different image sizes (256×256, 512×512, and 1024×1024). For each image, the moments were evaluated 20 times and the obtained computational times were averaged. The resulting averaged times are displayed in Table 7. The variances were negligible in all cases. The results show that the proposed algorithm saves computational time, especially for images of larger size.

Table 7

The obtained averaged execution time for the proposed method tP¯, and the direct method tD¯ and the ratio of these times for the four-level UWT. Description of the system used: DELL LATITUDE E6500, Intel(R) Core(TM)2 Duo CPU P9500 @ 2.53 GHz processor, 4.00 GB RAM, 32bit Windows Vista™ Business, and Matlab 7.5.0 (R2007b).

Image sizetP¯[s]tD¯[s]tD¯/tP¯[-]


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.27 We applied this noise removal algorithm to multimedia data that were artificially contaminated with a high amount of uniform white noise U(20,20) and also to real astronomical light images contaminated with dark current.


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=UWT{y} denotes the noisy observation (i.e., the acquired data), N=UWT{n} presents noise, and X=UWT{x} is the ideal noise-free image. The conditional mean of the posterior PDF pX|Y(x|y) provides the least square estimation of X. The formula for the MMSE16 estimator runs as

Eq. (32)

where pY|X(y|x) denotes the likelihood function, pX(x) represents the a priori model, and pN(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)

where τ and λ are the GLM parameters of the signal. The fourth moment of Y runs as

Eq. (34)

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)

From Eq. (33), we may derive

Eq. (36)

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 Mk(X)=1Ii=1I[XiE(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)

where x denotes the original image before contamination and x^ is the reconstructed denoised image. However, in most real-world applications, the precontamination image is not available and the RMSEout and RMSEin measures must be estimated in another way. RMSEin is directly related to the noise second moment M2(n)=M2(n)+E2(n), RMSE^in=M2(n). In Sec. 2, we exploit the model of the temperature dependency of dark current for the SBIG camera. The modeled second moment M2(n), the mean value E(n), and the fourth moment M4(n) represent a priori information about the noise component contained in the image.

Fig. 9

Testing images (256×256pixels) denoised by Bayesian shrinkage using the Db4 wavelet and 4 UWT decomposition levels (a) The noise House image [added uniform noise 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 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.


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

Eq. (38)

Using Eq. (38) we may write

Eq. (39)


Eq. (40)

If we consider the reconstructed light image to be given as x^(y)=x+nres, where the noise is suppressed to nres, then RMSEout can be estimated by

Eq. (41)

where E(nres)=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 texp. Then the k’th moment Mktexp(n) normalized to required exposure treq runs as

Eq. (42)


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


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: RMSE^in=1498, proposed: RMSE^out=325, direct: RMSE^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×1024pixels; inverted grayscale) without dark frame correction (T=277.36K, texp=1000s).


Fig. 11

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




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.26

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 ν 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×1024pixels), 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.


This work has been supported by the research project MSM 6046137306 of the Ministry of Education, Youth and Sports of the Czech Republic and by the Czech Grant Agency under Grants No. 205/09/1302 “Study of sporadic meteors and weak meteor showers using automatic video intensifiers cameras,” No. P102/10/1320 “Research and modelling of advanced methods of image quality evaluation,” and by the research program MSM 6840770014 “Research in the Field of Information and Communication Technology.”



A. Postigoet al., “Recent developments in the BOOTES experiment,” in AIP Conf. Proc. 662, 553 –555 (2003). Google Scholar


A. Pizurica, “Image denoising using wavelets and spatial context modeling,” Univ. Gent, Gent, Belgium (2002). Google Scholar


J. Švihlíket al., “Noise analysis of MAIA system and possible noise suppression,” Radioengineering, 20 (1), 110 –117 (2011). RAIOE3 Google Scholar


W. B. Davenport, Probability and Random Processes: An Introduction for Applied Scientists and Engineers, Prentice Hall, Inc., New York (1970). Google Scholar


C. Buil, CCD Astronomy, Willman-Bell, Inc., Richmond (1991). Google Scholar


R. Widenhornet al., “Temperature dependence of dark current in a CCD,” Proc. SPIE, 4669 193 –201 (2002). Google Scholar


GlennS. G., “Apparatus for elimination of dark current utilizing charge coupled devices,” U.S. Patent No. 4142213 (2008).


NobuhiroT., “Imaging device with elimination of dark current,” U.S. Patent No. 5216511 (2008).


A. K. OkyayC. O. ChuiK. C. Saraswat, “Leakage suppression by asymmetric area electrodes in metal-semiconductor-metal photodetectors,” Appl. Phys. Lett., 88 (6), 063506 (2006). APPLAB 0003-6951 Google Scholar


Santa Barbara Instrument Group (SBIG), A Division of Aplegen, Inc., Google Scholar


R. L. Baer, “A model for dark current characterization and simulation,” Proc. SPIE, 6068 606805 (2006). Google Scholar


M. Goeseleet al., “Entropy-based dark frame subtraction,” in Image Processing, Image Quality, Image Capture and Systems Conf. (PICS2001), (2001). Google Scholar


S. W. LyuE. P. Simoncelli, “Modeling multiscale subbands of photographic images with fields of Gaussian scale mixtures,” IEEE Trans. Pattern Anal. Mach. Intell., 31 (4), 693 –706 (2009). ITPIDJ 0162-8828 Google Scholar


M. AmirmazlaghaniH. Amindavar, “Two novel Bayesian multiscale approaches for speckle suppression in SAR images,” IEEE Trans. Geosci. Rem. Sens., 48 (7), 2980 –2993 (2010). IGRSD2 0196-2892 Google Scholar


J. Schmittet al., “Poisson denoising on the sphere: application to the Fermi gamma ray space telescope,” Astron. Astrophys., 517 (A26), 1 –13 (2010). AAEJAF 0004-6361 Google Scholar


M. RaphanE. P. Simoncelli, “Optimal denoising in redundant representations,” IEEE Trans. Image Process., 17 (8), 1342 –1352 (2008). IIPRE4 1057-7149 Google Scholar


A. Ralston, A First Course in Numerical Analysis, 1st Czech ed.Academia, Prague (1978). Google Scholar


M. Rerabeket al., “Space variant point spread function modeling for astronomical image data processing,” Proc. SPIE, 6691 66910T (2007). Google Scholar


K. Fliegelet al., “Meteor automatic imager and analyzer: system design and its parameters,” Proc. SPIE, 8135 81351S (2011). Google Scholar


P. Kotenet al., “Automatic video system for continues monitoring of the meteor activity,” Earth Moon Planets, 108 (1), 69 –76 (2011). 0167-9295 Google Scholar


S. G. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell., 11 (7), 674 –693 (1989). ITPIDJ 0162-8828 Google Scholar


J. L. StarckF. MurtaghA. Bijaoui, Image Processing and Data Analysis: The Multiscale Approach, Cambridge University Press, Cambridge, UK (1998). Google Scholar


Rice University, rice-wlet-tools, Google Scholar


M. Holschneideret al., Wavelets, Time-Frequency Methods and Phase Space: Real Time Algorithm for Signal Analysis with the Help of the Wavelet Transform, Springer-Verlag, Berlin (1989). Google Scholar


M. G. Bulmer, Principles of Statistics, Dover Publications, Inc., Mineola (1979). Google Scholar


J. Švihlík, “Modeling of scientific images using GMM,” Radioengineering, 18 (4), 579 –586 (2009). RAIOE3 Google Scholar


E. P. SimoncelliE. H. Adelson, “Noise removal via Bayesian wavelet coring,” in Proc. 3rd IEEE Int. Conf. Image Process., 379 –382 (1996). Google Scholar


H. Jeffrey, “An invariant form for the prior probability in estimation problems,” Proc. R. Soc. London, 186 (1007), 453 –461 (1946). PRSLAZ 0370-1662 Google Scholar


A. J. Izenman, “Recent developments in nonparametric density estimation,” J. Stat. Assoc., 86 (413), 205 –224 (1991). JSTNAL 0162-1459 Google Scholar


J. ŠvihlíkP. Páta, “Elimination of thermally generated charge in charged coupled devices using Bayesian estimator,” Radioengineering, 17 (2), 119 –124 (2008). RAIOE3 Google Scholar



Jan Švihlík received his MSc (Ing.) from the Czech Technical University in Prague (CTU), in 2005 and PhD from the CTU, in 2008. Now he is an assistant professor at the Institute of Chemical Technology Prague. His research interests include image processing, image denoising, astronomical images, and statistical models of images.


Karel Fliegel received his MSc degree in electrical engineering from the CTU in 2004 and PhD from the CTU, in 2011. Now he is an assistant professor at the CTU. His research interests include image processing, imaging systems, and image and video compression.


Jaromír Kukal obtained the MSc degree from the Institute of Chemical Technology in Prague, in 1978 and PhD from the University of Ostrava, in 2000. Now he is an associate professor at the CTU in Prague. His research interests include PCA, object classification, artificial neural networks, optimization, computer science, C++, and theory of database systems.


Eva Jerhotová obtained the MSc degree from the Institute of Chemical Technology in Prague in 2006, where she is currently finishing her PhD in technical cybernetics. Her research is focused on image denoising and segmentation.


Petr Páta received Mgr (MSc) degree in physics from the Faculty of Mathematics and Physics, Charles University in Prague, in 1996, then he started PhD study in radioelectronics at the Faculty of Electrical Engineering, the CTU in Prague. His PhD thesis was defended in 2002. He is currently working as associate professor with the department of radioelectronics of FEE CTU in Prague.


Stanislav Vítek received his MSc degree in electrical engineering from the CTU in 2002 and PhD from the CTU, in 2009. Now he is an assistant professor at the CTU. His research interests include signal and image processing and analysis, applications of astronomical image data processing, modeling of imaging systems, stereoscopic systems and assistive technology.


Pavel Koten graduated from Faculty of Mathematics and Physics of Charles University in Prague, in 1996. He has been at the astronomical institute since 1996. In 2001, he earned a PhD in astronomy at Charles University. His research fields include photometry, light curves and physical structure of faint television meteors, double station observations of meteors, trajectory computation, and image processing of video records.

© 2012 SPIE and IS&T 0091-3286/2012/$25.00 © 2012 SPIE and IS&T
Jan Svihlik, Eva Jerhotova, Karel Fliegel, Petr Páta, Stanislav Vítek, Jaromir Kukal, and Pavel Koten "Estimation of non-Gaussian noise parameters in the wavelet domain using the moment-generating function," Journal of Electronic Imaging 21(2), 023025 (22 June 2012).
Published: 22 June 2012

Back to Top