- Click thumbnails for details -
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
-
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:
-
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.