Fast Computation of Fourier Integral Operators

We introduce a general purpose algorithm for rapidly computing certain types of oscillatory integrals which frequently arise in problems connected to wave propagation and general hyperbolic equations. The problem is to evaluate numerically a so-called Fourier integral operator (FIO) of the form $\int e^{2\pi i \Phi(x,\xi)} a(x,\xi) \hat{f}(\xi) \mathrm{d}\xi$ at points given on a Cartesian grid. Here, $\xi$ is a frequency variable, $\hat f(\xi)$ is the Fourier transform of the input $f$, $a(x,\xi)$ is an amplitude and $\Phi(x,\xi)$ is a phase function, which is typically as large as $|\xi|$; hence the integral is highly oscillatory at high frequencies. Because an FIO is a dense matrix, a naive matrix vector product with an input given on a Cartesian grid of size $N$ by $N$ would require $O(N^4)$ operations. This paper develops a new numerical algorithm which requires $O(N^{2.5} \log N)$ operations, and as low as $O(\sqrt{N})$ in storage space. It operates by localizing the integral over polar wedges with small angular aperture in the frequency plane. On each wedge, the algorithm factorizes the kernel $e^{2 \pi i \Phi(x,\xi)} a(x,\xi)$ into two components: 1) a diffeomorphism which is handled by means of a nonuniform FFT and 2) a residual factor which is handled by numerical separation of the spatial and frequency variables. The key to the complexity and accuracy estimates is that the separation rank of the residual kernel is \emph{provably independent of the problem size}. Several numerical examples demonstrate the efficiency and accuracy of the proposed methodology. We also discuss the potential of our ideas for various applications such as reflection seismology.Comment: 31 pages, 3 figure

Scale-discretised ridgelet transform on the sphere

We revisit the spherical Radon transform, also called the Funk-Radon transform, viewing it as an axisymmetric convolution on the sphere. Viewing the spherical Radon transform in this manner leads to a straightforward derivation of its spherical harmonic representation, from which we show the spherical Radon transform can be inverted exactly for signals exhibiting antipodal symmetry. We then construct a spherical ridgelet transform by composing the spherical Radon and scale-discretised wavelet transforms on the sphere. The resulting spherical ridgelet transform also admits exact inversion for antipodal signals. The restriction to antipodal signals is expected since the spherical Radon and ridgelet transforms themselves result in signals that exhibit antipodal symmetry. Our ridgelet transform is defined natively on the sphere, probes signal content globally along great circles, does not exhibit blocking artefacts, supports spin signals and exhibits an exact and explicit inverse transform. No alternative ridgelet construction on the sphere satisfies all of these properties. Our implementation of the spherical Radon and ridgelet transforms is made publicly available. Finally, we illustrate the effectiveness of spherical ridgelets for diffusion magnetic resonance imaging of white matter fibers in the brain.Comment: 5 pages, 4 figures, matches version accepted by EUSIPCO, code available at http://www.s2let.or

Parabolic Molecules

Anisotropic decompositions using representation systems based on parabolic scaling such as curvelets or shearlets have recently attracted significantly increased attention due to the fact that they were shown to provide optimally sparse approximations of functions exhibiting singularities on lower dimensional embedded manifolds. The literature now contains various direct proofs of this fact and of related sparse approximation results. However, it seems quite cumbersome to prove such a canon of results for each system separately, while many of the systems exhibit certain similarities. In this paper, with the introduction of the notion of {\em parabolic molecules}, we aim to provide a comprehensive framework which includes customarily employed representation systems based on parabolic scaling such as curvelets and shearlets. It is shown that pairs of parabolic molecules have the fundamental property to be almost orthogonal in a particular sense. This result is then applied to analyze parabolic molecules with respect to their ability to sparsely approximate data governed by anisotropic features. For this, the concept of {\em sparsity equivalence} is introduced which is shown to allow the identification of a large class of parabolic molecules providing the same sparse approximation results as curvelets and shearlets. Finally, as another application, smoothness spaces associated with parabolic molecules are introduced providing a general theoretical approach which even leads to novel results for, for instance, compactly supported shearlets

Curvelets and Ridgelets

International audienceDespite the fact that wavelets have had a wide impact in image processing, they fail to efficiently represent objects with highly anisotropic elements such as lines or curvilinear structures (e.g. edges). The reason is that wavelets are non-geometrical and do not exploit the regularity of the edge curve. The Ridgelet and the Curvelet [3, 4] transforms were developed as an answer to the weakness of the separable wavelet transform in sparsely representing what appears to be simple building atoms in an image, that is lines, curves and edges. Curvelets and ridgelets take the form of basis elements which exhibit high directional sensitivity and are highly anisotropic [5, 6, 7, 8]. These very recent geometric image representations are built upon ideas of multiscale analysis and geometry. They have had an important success in a wide range of image processing applications including denoising [8, 9, 10], deconvolution [11, 12], contrast enhancement [13], texture analysis [14, 15], detection [16], watermarking [17], component separation [18], inpainting [19, 20] or blind source separation[21, 22]. Curvelets have also proven useful in diverse fields beyond the traditional image processing application. Let’s cite for example seismic imaging [10, 23, 24], astronomical imaging [25, 26, 27], scientific computing and analysis of partial differential equations [28, 29]. Another reason for the success of ridgelets and curvelets is the availability of fast transform algorithms which are available in non-commercial software packages following the philosophy of reproducible research, see [30, 31]

Continuous Curvelet Transform: I. Resolution of the Wavefront Set

We discuss a Continuous Curvelet Transform (CCT), a transform f → Γf (a, b, θ) of functions f(x1, x2) on R^2, into a transform domain with continuous scale a > 0, location b ∈ R^2, and orientation θ ∈ [0, 2π). The transform is defined by Γf (a, b, θ) = {f, γabθ} where the inner products project f onto analyzing elements called curvelets γ_(abθ) which are smooth and of rapid decay away from an a by √a rectangle with minor axis pointing in direction θ. We call them curvelets because this anisotropic behavior allows them to ‘track’ the behavior of singularities along curves. They are continuum scale/space/orientation analogs of the discrete frame of curvelets discussed in Candès and Donoho (2002). We use the CCT to analyze several objects having singularities at points, along lines, and along smooth curves. These examples show that for fixed (x0, θ0), Γf (a, x0, θ0) decays rapidly as a → 0 if f is smooth near x0, or if the singularity of f at x0 is oriented in a different direction than θ_0. Generalizing these examples, we state general theorems showing that decay properties of Γf (a, x0, θ0) for fixed (x0, θ0), as a → 0 can precisely identify the wavefront set and the H^m- wavefront set of a distribution. In effect, the wavefront set of a distribution is the closure of the set of (x0, θ0) near which Γf (a, x, θ) is not of rapid decay as a → 0; the H^m-wavefront set is the closure of those points (x0, θ0) where the ‘directional parabolic square function’ S^m(x, θ) = ( ʃ|Γf (a, x, θ)|^2 ^(da) _a^3+^(2m))^(1/2) is not locally integrable. The CCT is closely related to a continuous transform used by Hart Smith in his study of Fourier Integral Operators. Smith’s transform is based on strict affine parabolic scaling of a single mother wavelet, while for the transform we discuss, the generating wavelet changes (slightly) scale by scale. The CCT can also be compared to the FBI (Fourier-Bros-Iagolnitzer) and Wave Packets (Cordoba-Fefferman) transforms. We describe their similarities and differences in resolving the wavefront set

The curvelet transform for image denoising

We describe approximate digital implementations of two new mathematical transforms, namely, the ridgelet transform and the curvelet transform. Our implementations offer exact reconstruction, stability against perturbations, ease of implementation, and low computational complexity. A central tool is Fourier-domain computation of an approximate digital Radon transform. We introduce a very simple interpolation in the Fourier space which takes Cartesian samples and yields samples on a rectopolar grid, which is a pseudo-polar sampling set based on a concentric squares geometry. Despite the crudeness of our interpolation, the visual performance is surprisingly good. Our ridgelet transform applies to the Radon transform a special overcomplete wavelet pyramid whose wavelets have compact support in the frequency domain. Our curvelet transform uses our ridgelet transform as a component step, and implements curvelet subbands using a filter bank of a` trous wavelet filters. Our philosophy throughout is that transforms should be overcomplete, rather than critically sampled. We apply these digital transforms to the denoising of some standard images embedded in white noise. In the tests reported here, simple thresholding of the curvelet coefficients is very competitive with "state of the art" techniques based on wavelets, including thresholding of decimated or undecimated wavelet transforms and also including tree-based Bayesian posterior mean methods. Moreover, the curvelet reconstructions exhibit higher perceptual quality than wavelet-based reconstructions, offering visually sharper images and, in particular, higher quality recovery of edges and of faint linear and curvilinear features. Existing theory for curvelet and ridgelet transforms suggests that these new approaches can outperform wavelet methods in certain image reconstruction problems. The empirical results reported here are in encouraging agreement