This image is a denoising experiment using the curvelet transform. The original image is a synthetic noise-free seismogram simulated from a one-dimensional velocity model, courtesy of Felix Herrmann and Eric Verschuur:

---

Image (1) is a zoom-in. In real life, however, seismic data is corrupted by noise. In Image (2), we simulate degradation by adding gaussian white noise with standard deviation equal to one tenth of the total intensity, PSNR = 20.0 dB. Image (3) is the result of standard translation-invariant wavelet thresholding (see WaveLab), PSNR = 30.8 dB. Image (4) is the result of curvelet thresholding, PSNR = 37.6 dB.

The lesson of this experiment is that curvelets are an adequate tool to represent bandlimited wavefronts in an efficient manner, such as those present in reflection seismology data.

For reference, the transform used to generate these images is the FDCT via wrapping, complex-valued and with curvelets at the finest scale. The CurveLab toolbox contains, among others, an 'enhanced denoise' demo which tests the same thresholding algorithm on another image.

Curvelets at work

1. We start with a simple image of size 800 by 600, which we interpret as a matrix with entries one (black) or zero (white). This is a prototypical 'geometrical' image, expected to have an extremely sparse representation in curvelets:
2. We compute the Curvelet transform of the image and put to zero all the coefficients except the 5000 largest ones (in modulus) at the finest scale. This corresponds to a fraction of 0.14 percent of all the coefficients (take this info with a grain of salt since we started with a mostly empty image.) The inverse Curvelet transform is then applied and, after taking the real part, gives the image below. The most negative values are in white and the most positive values in black.

For reference, the transform used to generate these images is the FDCT via wrapping, complex-valued and with curvelets at the finest scale. The CurveLab toolbox contains, among others, a 'partial reconstruction' demo which tests the same thresholding algorithm on a more interesting image.

Other numerical experiments and pictures of curvelets can be found in the CurveLab toolbox.

Rendering of the second dyadic decomposition, by Gabrielle Candes.

The picture on the left shows the second dyadic decomposition (SDD), which is also the frequency partitioning of curvelets. The SDD was invented by Charles Fefferman in 1973 to answer mathematical questions of boundedness of Bochner-Riesz summation multipliers (some convolution operators).

Next, we show examples of curvelets at different scales, both in space and frequency:

Finally, the picture on the right shows the second dyadic decomposition in three dimensions, which is also the frequency partitioning of 3d curvelets.All the pictures are taken from the paper Fast Discrete Curvelet Transforms and refer to the wrapping variant, complex-valued, and with curvelets at the finest scale.

In this thought experiment, we consider the wave equation modeling acoustic, electromagnetic or elastic waves. If the initial condition is chosen localized and directional like a curvelet, then it will remain so for larger times as well---only it will travel along the light (or sound) rays. This is a result we proved in this paper, and remains true even if the index of refraction of the medium is smoothly varying. As a result the curvelet matrix of the Green's function exhibits an optimally sparse structure. This result does not hold for other transforms such as wavelets, Gabor, ridgelets, or contourlets.

This coherence of curvelets under the wave flow has far-reaching consequences. In image (2), we depict a wave front---where the underlying function is nonsmooth---initially in the shape of an ellipse. If we let it evolve under the wave equation, the ellipse will shrink to the fishtail pattern (the outgoing ellipse, if it exists, is not shown). A caustic has formed, and the wavefront seems to have lost some of its smoothness near the tips of the tail, named cusps. This is not a concern in phase-space, or in the curvelet domain, where the cusp pattern is in fact explained as a superposition of curvelets. Curvelets do not see caustics; keep in mind that the wave equation is linear and that the superposition principle is valid.

We have developed a numerical wave solver based on the same sparsity ideas.