3 Filtering

We now apply our digital transforms for removing noise from image data. The methodology is standard and is outlined mainly for the sake of clarity and self-containedness.

Suppose that one is given noisy data of the form

$$\begin{eqnarray*}x_{i,j} = f(i,j) + \sigma z_{i,j}, \end{eqnarray*}$$

where f is the image to be recovered and z is white noise, i.e. $z_{i,j} \stackrel{i.i.d.}{\sim} N(0,1)$. Unlike FFT's or FWT's, our discrete ridgelet (resp. curvelet) transform is not norm-preserving and, therefore, the variance of the noisy ridgelet (resp. curvelet) coefficients will depend on the ridgelet (resp. curvelet) index $\lambda$. For instance, letting F denote the discrete curvelet transform matrix, we have $F z \stackrel{i.i.d.}{\sim} N(0,F F^T)$. Because the computation of F F^T is prohibitively expensive, we calculated an approximate value $\tilde{\sigma}^2_\lambda$ of the individual variances using Monte-Carlo simulations where the diagonal elements of F F^T are simply estimated by evaluating the curvelet transforms of a few standard white noise images.

Let $y_\lambda$ be the noisy curvelet coefficients (y = F x). We use the following hard-thresholding rule for estimating the unknown curvelet coefficients:

$$\hat y_\lambda = y_\lambda \mbox{ if } \vert y_\lambda\vert/\sigma \geq k \tilde{\sigma}_\lambda$$ (6)

$$\hat y_\lambda = 0 \mbox{ if } \vert y_\lambda\vert/\sigma < k \tilde{\sigma}_\lambda.$$ (7)

In our experiments, we actually chose a scale-dependent value for k; we have k = 4 for the first scale (j = 1) while k = 3for the others (j > 1).

Poisson noise

Assume now that we have Poisson data x_i,j with unknown mean f(i,j). The Anscombe transformation (Anscombe 1948)

$$% \tilde{x} = 2\sqrt{x + \frac{3}{8}}$$ (8)

stabilizes the variance and we have $\tilde{x} = 2 \sqrt{f} + \epsilon$ where $\epsilon$ is a vector with independent and approximately standard normal components. In practice, this is a good approximation whenever the number of counts is large enough, greater than 30 per pixel, say.

Figure 6: Top left, part of Saturn image with a Gaussian noise. Top right, filtered image using the undecimated bi-orthogonal wavelet transform. Bottom left and right, filtered image by the à trous wavelet transform algorithm and the curvelet transform.

For small number of counts, a possibility is to compute the Radon transform of the image, and then to apply the Anscombe transformation to the Radon data. The rationale being that, roughly speaking, the Radon transform corresponds to a summation of pixel values over lines and that the sum of independent Poisson random variables is a Poisson random variable with intensity equal to the sum of the individual intensities. Hence, the intensity of the sum may be quite large (hence validating the Gaussian approximation) even though the individual intensities may be small. This might be viewed as an interesting feature as unlike wavelet transforms, the ridgelet and curvelet transforms tend to average data over elongated and rather large neighborhoods.

Gaussian and Poisson noise

The arrival of photons, and their expression by electron counts, on CCD detectors may be modeled by a Poisson distribution. In addition, there is additive Gaussian read-out noise. The Anscombe transformation (Eq. (8)) has been extended to take this combined noise into account. As an approximation, consider the signal's value, s_k, as a sum of a Gaussian variable, $\gamma$, of mean g and standard-deviation $\sigma$; and a Poisson variable, n, of mean m₀: we set $x = \gamma + \alpha n$ where $\alpha$ is the gain. The generalization of the variance stabilizing Anscombe formula is (Murtagh et al. 1995):

 

$$% \tilde x == \frac{2}{\alpha} \sqrt{\alpha x + \frac{3}{8} \alpha^2 + \sigma^2 - \alpha g}.$$ (9)

With appropriate values of $\alpha$, $\sigma$ and g, this reduces to Anscombe's transformation.

Then, for an image containing Poisson and Poisson+Gaussian noise, we apply first respectively the Anscombe and the Generalized Anscombe transform.

These variance stabilization transformations, it has been shown in Murtagh et al. (1995), are only valid for a sufficiently large number of counts (and of course, for a larger still number of counts, the Poisson distribution becomes Gaussian). The necessary average number of counts is about 20 if bias is to be avoided. Note that errors related to small values carry the risk of removing real objects, but not of amplifying noise. For Poisson parameter values under this threshold acceptable number of counts, the Anscombe transformation loses control over the bias. In this case, an alternative approach to variance stabilization is needed. It has been shown that the first step of the ridgelet transform consists in a Radon transform. As a Radon coefficient is an addition of pixel values along a line, the Radon transform of an image containing Poisson noise contains also Poisson noise. Then the Anscombe transform can be applied after the Radon transformation rather than on the original image. The advantage is that the number of counts per pixel will obviously be larger in the Radon domain than in the image domain, and the variance stabilization will be more robust.

Figure 7: Left, XMM/Newton image of the Kepler SN1604 supernova. Right, ridgelet filtered image.

Experiment

A Gaussian white noise with a standard deviation fixed to 20 was added to the Saturn image. We employed several methods to filter the noisy image:

1.

Thresholding of the Curvelet transform.

2.

Bi-orthogonal undecimated wavelet de-noising methods using the Dauchechies-Antonini 7/9 filters (FWT-7/9) and hard thresholding.

3.

À trous wavelet transform algorithm and hard thresholding.

Our experiments are reported in Fig. 6. The curvelet reconstruction does not contain the quantity of disturbing artifacts along edges that one sees in wavelet reconstructions. An examination of the details of the restored images is instructive. One notices that the decimated wavelet transform exhibits distortions of the boundaries and suffers substantial loss of important detail. The "à trous'' wavelet transform gives better boundaries, but completely omits to reconstruct certain ridges. In addition, it exhibits numerous small-scale embedded blemishes; setting higher thresholds to avoid these blemishes would cause even more of the intrinsic structure to be missed.

Further results are visible at the following URL: http://www-stat.stanford.edu/~jstarck.

Figure 7 shows an example of an X-ray image filtering by the ridgelet transform using such an approach. Figure 7 left and right shows respectively the XMM/Newton image of the Kepler SN1604 supernova and the ridgelet filtered image (using a five sigma hard thresholding).

Which transform should be chosen for a given data set?

We have introduced in this paper two new transforms, the ridgelet transform and the curvelet transform. Several other transforms are often used in astronomy, such the Fourier transform, the isotropic à trous wavelet transform and the bi-orthogonal wavelet transform. The choise of the best transform may be delicate. Each transform has its own domain of optimality:

• The Fourier transform for stationary process.
• The "à trous'' wavelet transform for isotropic features.
• The bi-orthogonal wavelet transform for features with a small anisotropy, typically with a width equals to half the length.
• The ridgelet wavelet transform for anisotropic features with a given length (i.e. block size).
• The curvelet transform for anisotropic features with different length and width equals to the square of the length.

Section 5 will show how several transforms can be used simultaneously, in order to benefit of the advantages of each of them.