## 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).^{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 camera^{10} 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.

More sophisticated techniques for dark current suppression are based on statistical modeling of the dark frame. In the spatial domain, Baer^{11} 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 Simoncelli^{13} 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^{14} 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.

## 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{\hspace{0.17em}}\text{\hspace{0.17em}}\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^{5} 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).

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

## (3)

$${M}_{l}=\frac{1}{I}\sum _{i=1}^{I}{({n}_{i}-\overline{n})}^{l},\phantom{\rule[-0.0ex]{2.0em}{0.0ex}}1\le l$$We found empirically that a satisfactory fit of the moment temperature dependencies ${M}_{l,i}(T)$ is achieved with the exponential model.^{17} Hence, using Eq. (1) we may write

## (5)

$$S({A}_{l},{B}_{l})=\sum _{i=1}^{I}{[\mathrm{ln}\text{\hspace{0.17em}}{A}_{l}-\frac{{B}_{l}}{{k}_{\beta}}{T}_{i}^{-1}-\mathrm{ln}({M}_{l,i}({T}_{i}))]}^{2}.$$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

## (6)

$$\frac{\partial S({A}_{l},{B}_{l})}{\partial {A}_{l}}=0,\phantom{\rule[-0.0ex]{2.0em}{0.0ex}}\frac{\partial S({A}_{l},{B}_{l})}{\partial {B}_{l}}=0.$$This leads to the following system of linear equations

where the matrix $\mathit{C}$ is given by## (8)

$$\mathbf{C}=\left[\begin{array}{ll}\sum _{i=1}^{I}({T}_{i}^{-2})& \sum _{i=1}^{I}({T}_{i}^{-1})\\ \sum _{i=1}^{I}({T}_{i}^{-1})& I\end{array}\right],$$## (9)

$$\mathbf{b}=\left[\begin{array}{l}\sum _{i=1}^{I}({T}_{i}^{-1}\xb7{M}_{l,i})\\ \sum _{i=1}^{I}({M}_{l,i})\end{array}\right].$$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

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.

## 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^{18} 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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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{\hspace{0.17em}}\text{\hspace{0.17em}}\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.

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.

## 3.1.

### 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 ${H}_{\mathrm{inv}}(u,v)$ is then given by

## (12)

$$|{H}_{\mathrm{inv}}(u,v)|=\{\begin{array}{rr}|S(u,v){|}^{-1},& |S(u,v)|>\delta \\ \delta ,& |S(u,v)|\le \delta \end{array},$$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

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

### 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

## (14)

$${M}_{r}=\frac{1}{I}\sum _{i=1}^{I}{({n}_{i}-\overline{n})}^{r},\phantom{\rule[-0.0ex]{2.0em}{0.0ex}}1\le r,$$## 3.3.

### 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 UWT^{23} for the video frame representation. The UWT is computed using a so-called à trous algorithm^{24} that produces the same number of wavelet coefficients at each scale (decomposition level). We use the following notation for the respective wavelet subbands: $\gamma {D}_{\xi}^{(v)}$, where $\xi $ 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).

## 3.4.

### 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 function^{4} and is defined by

The series expansion of ${e}^{un}$ 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 ${M}_{n}(u)$ is given by

## (16)

$${M}_{n}(u)={\int}_{-\infty}^{\infty}{e}^{un}p(n)\mathrm{d}n\phantom{\rule{0ex}{0ex}}={\int}_{-\infty}^{\infty}(1+un+\frac{{u}^{2}{n}^{2}}{2!}+\dots )p(n)\mathrm{d}n\phantom{\rule{0ex}{0ex}}=1+u{m}_{1}+\frac{{u}^{2}{m}_{2}}{2!}+\dots .,$$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 $\mathit{h}=[{h}_{1},{h}_{2}\dots {h}_{k}]$ while the down-sampling step is omitted. Hence, each wavelet coefficient is computed as the weighted sum of the independent random variables ${n}_{1},{n}_{2}\dots {n}_{k}$ (noise pixels) given as

The moment-generating function ${M}_{{S}_{k}}(u)$ of ${S}_{k}$ then runs as## (18)

$${M}_{{S}_{k}}(u)={M}_{{n}_{1}}({h}_{1}u){M}_{{n}_{2}}({h}_{2}u)\dots {M}_{{n}_{k}}({h}_{k}u).$$## (19)

$${M}_{{n}_{k}}({h}_{k}u)=[1+{h}_{1}u{m}_{1}+\frac{{({h}_{k}u)}^{2}{m}_{2}}{2!}+\frac{{({h}_{k}u)}^{3}{m}_{3}}{3!}\phantom{\rule{0ex}{0ex}}+\frac{{({h}_{k}u)}^{4}{m}_{4}}{4!}\dots ].$$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

Let us consider a zero-mean noise $n$ and its wavelet domain representation $N=\mathcal{U}\mathcal{W}\mathcal{T}\{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=[{h}_{1},{h}_{2}]$. The previously stated assumptions suggest the following moment relations. The second sample moment ${M}_{2}(N)$ in the wavelet domain is computed from ${M}_{2}(n)$ given by Similarly the fourth moment ${M}_{4}(N)$ computed from ${M}_{4}(n)$ given by Equations (21) and (22) may be generalized for filters with $k$ coefficients $\mathit{h}=[{h}_{1},{h}_{2}\dots {h}_{k}]$ given as## (24)

$${M}_{4}(N)=6{({M}_{2}(n))}^{2}[{h}_{1}^{2}{h}_{2}^{2}+{h}_{1}^{2}{h}_{3}^{2}+\dots +{h}_{1}^{2}{h}_{k}^{2}+{h}_{2}^{2}{h}_{3}^{2}\phantom{\rule{0ex}{0ex}}+{h}_{2}^{2}{h}_{4}^{2}+\dots ]+{M}_{4}(n)\sum _{i=1}^{k}{h}_{i}^{4}.$$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.

## 3.4.1.

#### 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)=\epsilon /2$ and the probability of flipping to 255 is $P(y=255)=\epsilon /2$. If the parameter $\epsilon $ approximately satisfies $\epsilon \le 0.05$ then

## (25)

$$6{({M}_{2}(n))}^{2}[{h}_{1}^{2}{h}_{2}^{2}+{h}_{1}^{2}{h}_{3}^{2}+\dots +{h}_{1}^{2}{h}_{k}^{2}+{h}_{2}^{2}{h}_{3}^{2}+{h}_{2}^{2}{h}_{4}^{2}+\dots ]\phantom{\rule{0ex}{0ex}}\ll {M}_{4}(n)\sum _{i=1}^{k}{h}_{i}^{4}.$$We tested these findings also on the 16-bpp dark frames^{26} 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).

## 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

## (27)

$${p}_{n}(n;\mu ,s,\nu )=\frac{{e}^{-{\left|\frac{n-\mu}{s}\right|}^{\nu}}}{Z(s,\nu )},\phantom{\rule[-0.0ex]{2.0em}{0.0ex}}n\in (-\infty ;\infty ),$$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 $\mu =0$. Then the second and the fourth moment run as

## (28)

$${m}_{2}(s,\nu )=\frac{{s}^{2}\mathrm{\Gamma}\left(\frac{3}{\nu}\right)}{\mathrm{\Gamma}\left(\frac{1}{\nu}\right)},\phantom{\rule[-0.0ex]{2em}{0.0ex}}{m}_{4}(s,\nu )=\frac{{s}^{4}\mathrm{\Gamma}\left(\frac{5}{\nu}\right)}{\mathrm{\Gamma}\left(\frac{1}{\nu}\right)}.$$In accordance with Ref. 27, the parameters estimation task can be simplified by using the kurtosis statistic

## (29)

$${\kappa}_{n}=\frac{{m}_{4}(s,\nu )}{{m}_{2}^{2}(s,\nu )}=\frac{\mathrm{\Gamma}\left(\frac{5}{\nu}\right)\mathrm{\Gamma}\left(\frac{1}{\nu}\right)}{{\mathrm{\Gamma}}^{2}\left(\frac{3}{\nu}\right)}.$$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.

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

Level | Wavelet: Haar | Wavelet: Db4 | ||||||
---|---|---|---|---|---|---|---|---|

Proposed | Direct | Proposed | Direct | |||||

M2 | M4 | M2 | M4 | M2 | M4 | M2 | M4 | |

1st | 8.4 | 189 | 8.4 | 188 | 8.4 | 185 | 8.4 | 184 |

2nd | 8.4 | 205 | 8.3 | 202 | 8.4 | 202 | 8.3 | 200 |

3rd | 8.4 | 209 | 8.0 | 187 | 8.4 | 208 | 7.9 | 182 |

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

Level | Wavelet: Haar | Wavelet: Db4 | ||||||
---|---|---|---|---|---|---|---|---|

Proposed | Direct | Proposed | Direct | |||||

M2 | M4 | M2 | M4 | M2 | M4 | M2 | M4 | |

1st | 1.9 | 198 | 1.9 | 197 | 1.9 | 235 | 1.9 | 233 |

2nd | 1.9 | 58 | 1.9 | 58 | 1.9 | 78 | 1.9 | 67 |

3rd | 1.9 | 23 | 1.9 | 25 | 1.9 | 31 | 2.0 | 27 |

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

Level | Wavelet: Haar | Wavelet: Db4 | ||||||
---|---|---|---|---|---|---|---|---|

Proposed | Direct | Proposed | Direct | |||||

M2 | M4 | M2 | M4 | M2 | M4 | M2 | M4 | |

1st | 9.0 | 243 | 9.0 | 239 | 9.0 | 243 | 9.0 | 241 |

2nd | 9.0 | 243 | 9.0 | 242 | 9.0 | 243 | 9.0 | 240 |

3rd | 9.0 | 243 | 9.0 | 246 | 9.0 | 243 | 9.0 | 246 |

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

Level | Wavelet: Haar | Wavelet: Db4 | ||||||
---|---|---|---|---|---|---|---|---|

Proposed | Direct | Proposed | Direct | |||||

M2 | M4 | M2 | M4 | M2 | M4 | M2 | M4 | |

1st | 3.8 | 73 | 3.8 | 98 | 3.8 | 79 | 3.8 | 100 |

2nd | 3.8 | 50 | 3.8 | 59 | 3.8 | 53 | 3.8 | 60 |

3rd | 3.8 | 44 | 3.8 | 54 | 3.8 | 46 | 3.8 | 54 |

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

Level | Wavelet: Haar | Wavelet: Db4 | ||||||
---|---|---|---|---|---|---|---|---|

Proposed | Direct | Proposed | Direct | |||||

M2 | M4 | M2 | M4 | M2 | M4 | M2 | M4 | |

1st | 1.5 | 7298 | 1.5 | 7318 | 1.5 | 8729 | 1.5 | 8774 |

2nd | 1.5 | 1830 | 1.5 | 1836 | 1.5 | 2615 | 1.5 | 2165 |

3rd | 1.5 | 462 | 1.5 | 458 | 1.5 | 787 | 1.5 | 530 |

## 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^{28} given by

## (30)

$${\mathcal{J}}_{\mathcal{D}}=\sum _{i=1}^{I}\left[{p}_{i}\text{\hspace{0.17em}}\mathrm{ln}\right(\frac{{p}_{i}}{0.5({p}_{i}+{H}_{i})})+{H}_{i}\text{\hspace{0.17em}}\mathrm{ln}(\frac{{H}_{i}}{0.5({p}_{i}+{H}_{i})}\left)\right],$$^{29}written as 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).

Uniform|Proposed | γD3(v) | γD3(h) | γD3(d) | Direct | γD3(v) | γD3(h) | γD3(d) |
---|---|---|---|---|---|---|---|

ν | 2.02 | 2.02 | 2.02 | ν | 2.02 | 1.99 | 2.01 |

s | 4.11 | 4.11 | 4.11 | s | 4.11 | 4.06 | 4.11 |

JD1 | 0.0027 | 0.0036 | 0.0028 | JD2 | 0.0026 | 0.0036 | 0.0027 |

JD12 | 0.0001 | 0.0001 | 0.0000 | JD12 | 0.0001 | 0.0001 | 0.0000 |

Gaussian|Proposed | γD3(v) | γD3(h) | γD3(d) | Direct | γD3(v) | γD3(h) | γD3(d) |

ν | 2.00 | 2.00 | 2.00 | ν | 2.01 | 1.92 | 1.94 |

s | 4.24 | 4.24 | 4.24 | s | 4.25 | 4.10 | 4.19 |

JD1 | 0.0024 | 0.0030 | 0.0026 | JD2 | 0.0024 | 0.0022 | 0.0025 |

JD12 | 0.0000 | 0.0008 | 0.0003 | JD12 | 0.0000 | 0.0008 | 0.0003 |

Real|Proposed | γD3(v) | γD3(h) | γD3(d) | Direct | γD3(v) | γD3(h) | γD3(d) |

ν | 1.87 | 1.87 | 1.87 | ν | 1.60 | 1.71 | 1.62 |

s | 2.64 | 2.64 | 2.64 | s | 2.4 | 2.51 | 2.41 |

JD1 | 0.0105 | 0.0063 | 0.0079 | JD2 | 0.0126 | 0.0074 | 0.0102 |

JD12 | 0.0105 | 0.0034 | 0.0089 | JD12 | 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.

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.

## 5.2.

### 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\times 256$, $512\times 512$, and $1024\times 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 size | tP¯[s] | tD¯[s] | tD¯/tP¯[-] |
---|---|---|---|

256×256 | 0.18 | 0.60 | 3.3 |

512×512 | 0.42 | 2.22 | 5.3 |

1024×1024 | 1.27 | 8.60 | 6.8 |

## 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.^{27} 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=\phantom{\rule{0ex}{0ex}}\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^{16} estimator runs as

## (32)

$$\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},\phantom{\rule{0ex}{0ex}}=\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},$$## (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),$$## (34)

$${m}_{4}(Y)=\frac{{\tau}^{4}\mathrm{\Gamma}\left(\frac{5}{\lambda}\right)}{\mathrm{\Gamma}\left(\frac{1}{\lambda}\right)}+\frac{6{s}^{2}{\tau}^{2}\mathrm{\Gamma}\left(\frac{3}{\nu}\right)\mathrm{\Gamma}\left(\frac{3}{\lambda}\right)}{\mathrm{\Gamma}\left(\frac{1}{\nu}\right)\text{\hspace{0.17em}}\mathrm{\Gamma}\left(\frac{1}{\lambda}\right)}+\frac{{s}^{4}\mathrm{\Gamma}\left(\frac{5}{\nu}\right)}{\mathrm{\Gamma}\left(\frac{1}{\nu}\right)}\phantom{\rule{0ex}{0ex}}={m}_{4}(X)+6{m}_{2}(N){m}_{2}(X)+{m}_{4}(N).$$## (35)

$${\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)}\phantom{\rule{0ex}{0ex}}=\frac{{m}_{4}(Y)-{m}_{4}(N)-6{m}_{2}(N)[{m}_{2}(Y)-{m}_{2}(N)]}{{({m}_{2}(Y)-{m}_{2}(N))}^{2}}.$$## (36)

$$\tau =\sqrt{[{m}_{2}(Y)-{m}_{2}(N)]\frac{\mathrm{\Gamma}\left(\frac{1}{\lambda}\right)}{\mathrm{\Gamma}\left(\frac{3}{\lambda}\right)}}.$$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

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.## 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 | RMSEin | RMSEout|proposed | RMSEout|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

## (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).$$## (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}})},$$## (42)

$${M}_{k}^{{t}_{\mathrm{req}}}(n)=\frac{1}{I}\sum _{i=1}^{I}{\left(\frac{{t}_{\mathrm{req}}{n}_{i}}{{t}_{\mathrm{exp}}}\right)}^{k}={\left(\frac{{t}_{\mathrm{req}}}{{t}_{\mathrm{exp}}}\right)}^{k}{M}_{k}^{{t}_{\mathrm{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.

## 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.^{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 $\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{\hspace{0.17em}}\text{\hspace{0.17em}}\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.

## Acknowledgments

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

## References

*Fermi*gamma ray space telescope,” Astron. Astrophys. 517(A26), 1–13 (2010).AAEJAF0004-6361 Google Scholar

## Biography

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