A Woiselle
email: arnaud.woiselle@cea.fr

J-L Starck
email: jsarck@cea.fr

J Fadili

3D Curvelet transforms and Astronomical Data Restoration

Keywords: 3D transform, curvelet, wavelet, multi-resolution representations, astrophysics, galaxy distribution

. The first one, called BeamCurvelet transform, is well designed for representing 1D filaments in a 3D space, while the second one, the RidCurvelet transform, is designed to analyze 2D surfaces. We show that these constructions can be useful for different applications such as filament detection, denoising or inpainting. Hence, they could lead to alternative approaches for analyzing 3D cosmological data sets, such as catalogs of galaxies, λCDM simulation or weak lensing data.

Introduction

Sparse representations such as wavelets or curvelets have been very successful for 2D image processing. Impressive results were obtained for many applications such as compression (see [START_REF] Chandrasekaran | Representation and compression of multidimensional piecewise functions using surflets[END_REF] for an example of Surflet compression; the new image standard JPEG2000 is based on wavelets rather than DCT like JPEG), denoising [START_REF] Starck | Very high quality image restoration by combining wavelets and curvelets[END_REF][START_REF] Starck | The curvelet transform for image denoising[END_REF][START_REF] Hennenfent | Seismic denoising with nonuniformly sampled curvelets[END_REF], contrast enhancement [START_REF] Starck | Gray and color image contrast enhancement by the curvelet transform[END_REF], inpainting [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis[END_REF][START_REF] Fadili | Inpainting and zooming using sparse representations[END_REF] or deconvolution [START_REF] Starck | Deconvolution based on the curvelet transform[END_REF][START_REF] Starck | Wavelets and curvelets for image deconvolution: a combined approach[END_REF]. Curvelets [START_REF] Starck | The curvelet transform for image denoising[END_REF][START_REF] Candès | New Tight Frames of Curvelets and Optimal Representations of Objects with C2 Singularities[END_REF], Bandelets [START_REF] Peyre | Discrete Bandelets with Geometric Orthogonal Filters[END_REF] and Contourlets [START_REF] Do | The contourlet transform: an efficient directional multiresolution image representation[END_REF] were designed to well represent edges in an image while wavelets are especially efficient for isotropic feature analysis.

In 3D, the separable Wavelet transform (decimated, undecimated, or any other kind) and the Discrete Cosine transform are certainly the most known decompositions [START_REF] Selesnick | Video denoising using 2D and 3D dual-tree complex wavelet transforms[END_REF][START_REF] Dima | Semiautomatic quality determination of 3D confocal microscope scans of neuronal cells denoised by 3D wavelet shrinkage[END_REF][START_REF] Chen | Breast volume denoising and noise characterization by 3D wavelet transform[END_REF]. The DCT is mainly used in video compression, but has also been used in denoising [START_REF] Rusanovskky | Video Denoising Algorithm in Sliding 3D DCT Domain[END_REF]. A lot of effort has been made in the last five years to build sparse 3D data representations, which represent better geometrical features contained in the data. The 3D beamlet transform [START_REF] Donoho | Fast X-Ray and Beamlet Transforms for Three-Dimensional Data[END_REF] and the 3D ridgelet transform [START_REF] Starck | Analysis of the spatial distribution of galaxies by multiscale methods[END_REF] were respectively designed for 1D and 2D features detection. Video denoising using the ridgelet transform was proposed in [START_REF] Carre | 3D fast ridgelet transform[END_REF]. The 3D fast curvelet transform [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF] consists in paving the Fourier domain with angular wedges in dyadic concentric squares using the parabolic scaling law to fix the number of angles depending on the scale, and has atoms designed for representing surfaces in 3D. The Surflet transform [START_REF] Chandrasekaran | Surflets : A sparse representation for multidimensional functions containing smooth discontinuities[END_REF] -a d-dimensional extension of the 2D wedgelets [START_REF] Donoho | Wedgelets: nearly minimax estimation of edges[END_REF][START_REF] Romberg | Multiscale wedgelet image analysis: fast decompositions and modeling[END_REF] -has recently been studied for compression purposes [START_REF] Chandrasekaran | Representation and compression of multidimensional piecewise functions using surflets[END_REF]. Surflets are an adaptive transform estimating each cube of a quad-tree decomposition of the data by two regions of constant value separated by a polynomial surface. Another possible representation uses the Surfacelets developed by Do and Lu [START_REF] Lu | 3-d directional filter banks and surfacelets[END_REF]. It relies on the combination of a Laplacian pyramid and a d-dimensional directional filter bank. Surfacelets produce a tiling of the Fourier space in angular wedges in a way close to the curvelet transform, and can be interpreted as a 3D adaptation of the 2D contourlet transform. This transformation has recently also been applied to video denoising [START_REF] Lu | Multidimensional directional filter banks and surfacelets[END_REF].

3D multiscale transforms and cosmological data set

Different statistical measures have been used in the cosmological literature to quantitatively describe the cosmic texture [START_REF] Martinez | Wavelet analysis of the multifractal character of the galaxy distribution[END_REF], i.e. the large-scale structure of the universe showing intricate patterns with filaments, clusters, and sheet-like arrangements of galaxies encompassing large nearly empty regions, the so-called voids. Wavelets have been used for many years [START_REF] Escalera | New evidence for subclustering in the Coma cluster using the wavelet analysis[END_REF][START_REF] Slezak | Objective detection of voids and high density structures in the first CfA redshift survey slice[END_REF][START_REF] Martinez | Wavelet analysis of the multifractal character of the galaxy distribution[END_REF], and it has been shown that denoising the galaxy distribution using wavelet instead of the standard Gaussian filtering, allows us to better preserve structure at different scales, and therefore better constrain our cosmological models [START_REF] Martínez | Morphology of the Galaxy Distribution from Wavelet Denoising[END_REF].

Noise is also a problem of major concern for N-body simulations of structure formation in the early Universe and it has been shown that using wavelets for removing noise from N-body simulations is equivalent to simulations with two orders of magnitude more particles [START_REF] Romeo | N-body simulations with two-orders-of-magnitude higher performance using wavelets[END_REF][START_REF] Romeo | A wavelet add-on code for new-generation n-body simulations and data de-noising (jofiluren)[END_REF].

Finally, 3D walelets, ridgelet and beamlet were also used in order to extract statistical information from galaxy catalogs [START_REF] Starck | Analysis of the spatial distribution of galaxies by multiscale methods[END_REF] and compare our data set to simulations obtained from different cosmological models.

Why new transforms ?

These 3D transforms all aim at representing the data using a minimal number of active coefficients, and by construction are better adapted to capture a specific kind of pattern. For the wavelet transform the pattern is smooth and isotropic, while for the DCT it is oscillating in all directions. All previously mentioned 3D transforms, except the beamlet transform, use plate-like functions, useful to represent surfaces in a 3D volume. The beamlet is therefore the only existing geometric decomposition allowing an efficient detection of filaments in 3D. Two relatively different implementations have been proposed. One [START_REF] Donoho | Fast X-Ray and Beamlet Transforms for Three-Dimensional Data[END_REF] suffers the lack of any reconstruction algorithm, and the other one [START_REF] Starck | Analysis of the spatial distribution of galaxies by multiscale methods[END_REF] is only optimal for the detection of filaments with a specific filament size. This has motivated the design of the FABT transform in biology [START_REF] Berlemont | Detection of curvilinear objects in biological noisy image using feature-adapted fast slant stack[END_REF], but this decomposition is also limited, since its optimality is only for filament of a given width. None of these transforms presents the nice scaling properties similar to those of 2D curvelet transform.

This paper

We propose in this paper a new transform, the BeamCurvelet transform, which is well adapted for the detection of filaments of different sizes and widths. A minor variation in its construction leads to another transform with plate-like elements, the RidCurvelet transform. Following the notations in [START_REF] Donoho | Tight frames of k-plane ridgelets and the problem of representing objects that are smooth away from d-dimensional singularities in R n . Applied Mathematics[END_REF], a more standard name for these representations would be "Local-k plane Ridgelet Bases in n-D : LRB(k,n)", thus leading respectively to LRB [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF][START_REF] Candès | Ridgelets: the key to high dimensional intermittency[END_REF] and LRB [START_REF] Berlemont | Detection of curvilinear objects in biological noisy image using feature-adapted fast slant stack[END_REF][START_REF] Candès | Ridgelets: the key to high dimensional intermittency[END_REF]. Both constructions have interesting scaling properties, and offer exact reconstruction. They can therefore be used for different applications such as denoising or inpainting. This paper is organized as follows : In the second and third section, we show how we can extend the 2D curvelet transform to 3D, leading to two new decompositions, the BeamCurvelet and the RidCurvelet. In section 4, we investigate different approaches for BeamCurvelet /Rid-Curvelet denoising. The last section presents applications how these two transforms can be used for inpainting.

The BeamCurvelet transform

The 2D Curvelet transform

In order to understand our construction for 3D Curvelets, we first recall the simpler but very similar 2D case.

The first generation curvelets [START_REF] Starck | The curvelet transform for image denoising[END_REF][START_REF] Candès | Curvelets, multiresolution representation, and scaling laws[END_REF] were built using the isotropic undecimated wavelet transform [START_REF] Starck | Astronomical Image and Data Analysis. Astronomical image and data analysis[END_REF] and the ridgelet transform [START_REF] Candès | Ridgelets: the key to high dimensional intermittency[END_REF].

The ridgelet transform is useful for representing global lines crossing an image on its full length the same way the wavelets represent isolated isotropic singularities and shapes. This property is obtained by the Radon transform, which transforms lines into points. A ridgelet function is indexed by a scale parameter a, a position b and an orientation θ.

Let γ = (a, b, θ) ∈ R * + × R × [0, 2π[. Given a wavelet function ψ, we define a ridgelet ψ γ by ψ γ = a -1/2 ψ ((x 1 cos θ + x 2 sin θ -b)/a) .
The ridgelet transform is implemented as a line extraction in Fourier domain, using the projectionslice theorem. The process is shown in Figure 1. The main drawback of this transform is that

Original image FFT2D F F T 1 D -1 WT1D

Angle Angle

Wavelet scales lines span the whole image and thus aren't well localized in space. The idea was then to make it local and multiscale by applying it blockwise on a multiscale isotropic wavelet transform. The key property is a law forcing the size of the blocks from one scale to the next, following a parabolic scaling which states that the number of blocks to get in a scale is downsized by a factor four each other scale (the size of the blocks is doubled each other scale on an isotropic undecimated wavelet transform, and divided by two on a decimated one). Figure 2 shows the global process described above. The ridgelet transform is implemented through a Radon transform in Fourier domain, followed by a one-dimensional wavelet transform. As in two dimensions, the 3D first generation Curvelet transforms we develop in this paper are based on the Radon transform applied to localized blocks of a given size of a 3D spectral decomposition of the data. There are two ways of extending the Radon transform in three dimensions, which lead to the two transforms described below. The first one is obtained by projecting along 2D planes (3D partial Radon transform), which leads to the BeamCurvelets, and the second one by projecting only along 1D lines (3D Radon transform), which leads to the RidCurvelets.

Radon Transform

Original image

The 3D Continuous BeamCurvelet Transform

In order to separate the signal into spectral bands, we use a filter-bank. Let N * be the set of strictly positive integers. Given a smooth function ψ ∈ L 2 (R 3 ) : ∀s ∈ N * , ψ 2s = 2 6s ψ(2 2s •) extracting the frequencies around |ξ| ∈ [2 2s , 2 2s+2 ], and a low-pass filter ψ 0 for |ξ| ≤ 1. We get a partition of unity :

| ψ0 (ξ)| 2 + s>0 | ψ2s (ξ)| 2 = 1 (1)
Let P 0 f = ψ 0 * f and ∆ s f = ψ 2s * f , where * is the convolution product. We can represent any signal f as (P 0 f, ∆ 1 f, ∆ 2 f, ...).

We tile the spatial domain of each scale s with a set Q s of regions Q of size 2 s :

Q = Q(s, k 1 , k 2 , k 3 ) = k 1 2 s , k 1 + 1 2 s × k 2 2 s , k 2 + 1 2 s × k 3 2 s , k 3 + 1 2 s ⊂ [0, 1] 3 (2)
with smooth windows w Q localized near Q, and verifying Q∈Q s w 2 Q = 1, with

Q s = Q(s, k 1 , k 2 , k 3 )|(k 1 , k 2 , k 3 ) ∈ [0, 2 s ) 3 . (3)
Each element of the frequency-space w Q ∆ s is transported to [0, 1] 3 by the transport operator

T Q : L 2 (Q) → L 2 ([0, 1] 3 ) applied to f ′ = w Q ∆ s f T Q : L 2 (Q) → L 2 ([0, 1] 3 ) (T Q f ′ )(x 1 , x 2 , x 3 ) = 2 -s f ′ k 1 + x 1 2 s , k 2 + x 2 2 s , k 3 + x 3 2 s . ( 4
)
For each scale s, we have a space-frequency tiling operator g Q , the output of which lives on [0, 1] 3

g Q = T Q w Q ∆ s . (5)
We can apply a 3D Beamlet transform [START_REF] Donoho | Fast X-Ray and Beamlet Transforms for Three-Dimensional Data[END_REF][START_REF] Donoho | Beamlets and multiscale image analysis[END_REF] on each block of each scale, by projecting on the beamlet functions :

β σ,κ 1 ,κ 2 ,θ 1 ,θ 2 (x 1 , x 2 , x 3 ) = σ -1/2 φ((-x 1 sin θ 1 + x 2 cos θ 1 + κ 1 )/σ, (6)
(x 1 cos θ 1 cos θ 2 + x 2 sin θ 1 cos θ 2 -x 3 sin θ 2 + κ 2 )/σ).
where σ is the Beamlet scale parameter, (θ 1 , θ 2 ) the orientation parameter and (κ 1 , κ 2 ) the location parameter, which is two dimensional because the beamlet transform integrates the data over one dimension through the partial Radon transform (see section 2.3). φ ∈ L 2 (R 3 ) is a smooth function satisfying the following admissibility condition

s∈Z φ 2 (2 s u) = 1, ∀u ∈ R 2 . (7)
Finally, the BeamCurvelet transform of a 3D function f is

BC f = T Q w Q ∆ s f, β σ,κ 1 ,κ 2 ,θ 1 ,θ 2 : s ∈ N * , Q ∈ Q s . (8)
As we can see, a BeamCurvelet function is parametrized in scale (s, σ), position (Q, κ 1 , κ 2 ), and orientation (θ 1 , θ 2 ). The following sections describe the discretization and the effective implementation of such a transform.

Discretization

For convenience, and as opposed to the continuous notations, the scales are now numbered from 0 to J, from the finest to the coarsest. As seen in the continuous formulation, the transform operates in four main steps.

First the frequency decomposition is obtained by applying a 3D wavelet transform on

the data with a wavelet compactly supported in Fourier space like the pyramidal Meyer wavelets with low redundancy [START_REF] Starck | Image reconstruction by the wavelet transform applied to aperture synthesis[END_REF], or using the 3D isotropic à trou wavelets. 2. Each wavelet scale is then decomposed in small cubes of a size following the parabolic scaling law, forcing the block size B s with the scale size N s according to the formula

B s N s = 2 s/2 B 0 N 0 , (9)
where N 0 and B 0 are the finest scale's dimension and block size.

3. Then, we apply a partial 3D Radon transform on each block of each scale. This is accomplished by integrating the blocks along lines at every direction and position. For a fixed direction (θ 1 , θ 2 ), the summation gives us a plane. Each point on this plane represents a line in the original cube. We obtain projections of the blocks on planes passing through the origin at every possible angle. 4. At last, we apply a two-dimensional wavelet transform on each Partial Radon plane.

Steps 3 and 4 represent the Beamlet transform of the blocks. The 3D Beamlet atoms aim at representing filaments crossing the whole 3D space. They are constant along a line and oscillate like φ in the radial direction. Arranged blockwise on a 3D isotropic wavelet transform, and following the parabolic scaling, we obtain the BeamCurvelet transform. Repeat over all its parallels to get the dark plane and apply a 2D wavelet transform within that plane. Repeat for all the directions to get the 3D Beamlet transform. See the text (section 3) for a detailed explanation and implementation clues.

Wavelet transform

Original datacube

(θ1, θ2) (θ 1 , θ 2 ) d ir e c ti o n s (θ 1 , θ 2 ) d ir e c ti o n s 3D Beamlet transform

Algorithm summary

As for the 2D Curvelets, the 3D BeamCurvelet transform is implemented effectively in the Fourier domain. Indeed, the integration along the lines (3D partial Radon transform) becomes a simple plane extraction in Fourier space, using the d-dimensional projection-slice theorem, which states that the Fourier transform of the projection of a d-dimensional function onto an mdimensional linear submanifold is equal to an m-dimensional slice of the d-dimensional Fourier transform of that function through the origin in the Fourier space which is parallel to the projection submanifold. In our case, d = 3 and m = 2. Algorithm 1 summarizes the whole process. Given a Cube of size N × N × N, a cubic block of length B s at scale s, and J + 1 scales, the redundancy can be calculated as follows : According to the parabolic scaling, ∀s > 0 : B s /N s = 2 s/2 B 0 /N 0 . The redundancy induced by the 3D wavelet tansform is

R w = 1 N 3 J s=0 N 3 s , (10)
with N s = 2 -s N for pyramidal Meyer wavelets (the case used on experiments described in sections 4 and 5), and thus B s = 2 -s/2 B 0 according to the parabolic scaling (see equation 9). The partial Radon transform of a cube of size B 3 s has a size 3B 2 s × B 2 s to which we apply 2D decimated orthogonal wavelets with no redundancy. There are (ρN s /B s ) 3 blocks in each scale because of the overlap factor (ρ ∈ [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF][START_REF] Berlemont | Detection of curvilinear objects in biological noisy image using feature-adapted fast slant stack[END_REF]) in each direction. So the complete redundancy of the transform using the Meyer wavelets is

R = 1 N 3 J-1 s=0 ρ N s B s 3 3B 4 s + N 3 J N 3 = 3ρ 3 J-1 i=0 B s 2 -3s + 2 -3J (11) = 3ρ 3 B 0 J-1 s=0 2 -7s/2 + 2 -3J (12)
= O 3ρ 3 B 0 when J → ∞ (13) R(J = 1) = 3ρ 3 B 0 + 1 8 (14) R(J = ∞) ≈ 3.4ρ 3 B 0 ( 15
)

Inverse BeamCurvelet Transform

Because all its components are invertible, the BeamCurvelet transform is invertible and the reconstruction error is comparable to machine precision. Algorithm 2 details the reconstruction steps.

An example of a 3D BeamCurvelet atom is represented in Figure 5. The BeamCurvelet atom is a collection of straight smooth segments well localized in space. Across the transverse plane, the BeamCurvelets exhibit a wavelet-like oscillating behavior.

The 3D RidCurvelet transform

The continuous transform

As referred to in 2.2, the second extension of the curvelet transform in 3D is obtained by using the 3D Ridgelet transform [START_REF] Candès | Ridgelets: the key to high dimensional intermittency[END_REF] Apply a 3D inverse isotropic wavelet transform end is given by :

ρ σ,κ,θ 1 ,θ 2 (x 1 , x 2 , x 3 ) = σ -1/2 φ 1 σ (x 1 cos θ 1 cos θ 2 + x 2 sin θ 1 cos θ 2 + x 3 sin θ 2 -κ) , (18)
where σ and κ are respectively the scale and position parameters, and φ ∈ L 2 (R) satisfies equation 7. The global RidCurvelet transform of a 3D function f is then

RC f = T Q w Q ∆ s f, ρ σ,κ,θ 1 ,θ 2 : s ∈ N * , Q ∈ Q s . (19)

Discretization

The discretization is made the same way, the sums over lines becoming sums over the planes of normal direction (θ 1 , θ 2 ), which gives us a line for each direction. The 3D Ridge function is useful for representing planes in a 3D space. It is constant along a plane and oscillates like φ in the normal direction. The main steps of the Ridgelet transform are depicted in figure 6.

Implementation

The RidCurvelet transform is also implemented in Fourier domain, the integration along the planes becoming a line extraction in the Fourier domain. The overall process is shown in Figure 7, and Algorithm 3 summarizes the implementation.

Radon transform 1D Wavelet transform

Sum over the planes at a given direction Repeat over all its parallels to get the (θ 1 , θ 2 ) line and apply a 1D wavelet transform on it. Repeat for all the directions to get the 3D Ridgelet transform. See the text (section 3) for a detailed explanation and implementation clues.

(θ 1 , θ 2 ) position (θ 1 , θ 2 ) direction wavelet scales (θ 1 , θ 2 ) direction

Properties

The RidCurvelet transform forms a tight frame. Additionally, given a 3D cube of size N × N × N, a block of size-length B s at scale s, and J + 1 scales, the redundancy is calculated as follows :

The Radon transform of a cube of size B 3 s has a size 3B 2 s × B s , to which we apply a pyramidal 1D wavelet of redundancy 2, for a total size of 3B each scale because of the overlap factor (ρ ∈ [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF][START_REF] Berlemont | Detection of curvilinear objects in biological noisy image using feature-adapted fast slant stack[END_REF]) in each direction. Therefore, the complete redundancy of the transform using many scales of 3D Meyer wavelets is

R = J-1 s=0 6B 3 s ρ N s B s 3 + 2 -3J = 6ρ 3 J-1 s=0 2 -3s + 2 -3J (20)
R = O(6ρ 3 ) when J → ∞ (21) R(J = 1) = 6ρ 3 + 1/8 (22) R(J = ∞) ≈ 6.86ρ 3 . ( 23
)

Inverse RidCurvelet Transform

The RidCurvelet transform is invertible and the reconstruction error is comparable to machine precision. Algorithm 4 details the reconstruction steps.

An example of a 3D RidCurvelet atom is represented in Figure 8. The RidCurvelet atom is composed of planes with values oscillating like a wavelet in the normal direction, and well localized due to the smooth function used to extract blocks on each wavelet scale.

Denoising

Introduction

In sparse representations, the simplest denoising methods are performed by a simple thresholding of the discrete curvelet coefficients. The threshold level is usually taken as three times the noise standard deviation, such that for an additive gaussian noise, the thresholding operator kills all noise coefficients except a small percentage, keeping the big coefficients containing information. The threshold we use is often a simple κσ, with κ ∈ [START_REF] Candès | Ridgelets: the key to high dimensional intermittency[END_REF][START_REF] Candès | Curvelets, multiresolution representation, and scaling laws[END_REF] to 0.27% and 6.3 • 10 -5 false detections. Sometimes we use a higher κ for the finest scale [START_REF] Starck | The curvelet transform for image denoising[END_REF]. Other methods exist, that estimate automatically the threshold to use in each band like the False Discovery Rate (see [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF][START_REF] Miller | Controlling the false-discovery rate in astrophysical data analysis[END_REF]). The correlation between neighbor coefficients intra-band and/or inter-band may also be taken into account (see [START_REF] Sendur | Bivariate shrinkage with local variance estimation[END_REF][START_REF] Sendur | Bivariate shrinkage functions for wavelet-based denoising exploiting interscale dependency[END_REF]). In this paper, in order to evaluate the different transforms, we use a κσ Hard Thresholding in our experiments.

Algorithm

Due to the high redundancy of the proposed transforms, and because of the huge size of 3D data, the transforms were implemented in a filtering-oriented way, in order to spare the machine resources, and enable easy multi-threading. The two new Curvelet transforms operate blockwize, and when there is no intra-block correlation taken into account in the denoising process, each block can be treated independently and on a different processor. Therefore, we never have the full transform in memory, only the wavelet transform of the data, and the Curvelet transform of one block (times the number of CPUs if working on a cluster). Algorithm 5 summarizes the denoising process. Using this algorithm, the memory used to filter a cube of any size using the Curvelet transforms is about twice that of the isotropic Meyer wavelet transform, as the size of one block is usually negligible compared to the size of the data. As an example, a data cube of size 128 × 128 × 128 (8MB with 32bit floating points) transformed with the BeamCurvelet and a block size of 17 with full 50% overlapping requires 3.3GB of memory, while it runs under 70MB with the low memory filtering algorithm. The transform and reconstruction with the RidCurvelet (BeamCurvelet) transform take respectively about 25(260) and 35(460) seconds without overlapping on a laptop (using a single core at 2.16GHz). With full overlapping, it takes about six times more time, but as the code is parallelized with OpenMP, it can use all the CPUs of the computer to lower the execution time dramatically.

Experiments

We have performed a few denoising experiments to emphasize the applicability domain and specificity of each transform.

Toy model : Structure detection

First of all, in order to have a quantitative indicator of the efficiency of the two curvelet transforms to approximate and detect specific structures, we transform two cubes, one containing a plane and the second a filament. Table 1 contains the maximum coefficients of the transformed cubes, showing the level of detection of each structure in a given space. As we can see, the RidCurvelet has the greatest coefficient for the plane and the BeamCurvelet for the filament. In many applications (denoising, inpainting, morphological component analysis, ...), the process is a thresholding of the coefficients at a level fixed by the parameters of the algorithm and the noise level. Having many coefficients above a given level (or much energy), is an indication of effectiveness.

In order to see the gain when using one transform rather than another, we threshold the RidCurvelet transform of the cube containing the plane at the maximum BeamCurvelet detection level (which is also far above the wavelet one). This means that, had we used another transform besides the RidCurvelet at the same detection level, we would not have detected anything. The reconstruction is shown in figure 9 (top row). We have also performed the same experiment with the BeamCurvelet transform of the cube containing the filament with a threshold equal to the level where the three other tested transforms don't detect anything. As expected, there is much information kept by the most adapted transform, even at a very high level of thresholding.

Structure denoising

Another way to see the power of each transform when associated to the right structures is to denoise a synthetic cube containing plane-and filament-like structures. Figure 10 shows a cut and a projection of the test cube containing parts of spherical shells and a spring-shaped filament. We denoise this cube using wavelets, RidCurvelets and BeamCurvelets. As shown in figure 11, the RidCurvelets denoise correctly the shells but poorly the filament, the BeamCurvelets restore Table 2: PSNR of the denoised synthetic cube using wavelets, RidCurvelets or BeamCurvelets the helix more properly while slightly underperforming for the shells, and wavelets are poor on the shell and give a dotted result and misses the faint parts of both structures. The PSNRs obtained with each transform are reported in Table 2. Here, the Curvelet transforms did very well for a given kind of features, and the wavelets were better on the signal power. In Starck et al. [START_REF] Starck | Very high quality image restoration by combining wavelets and curvelets[END_REF], in the framework of 3D image denoising, it was advocated to combine several transforms in order to benefit from the advantages of each of them.

Combined Filtering

Let x be the unknown data, and z a white gaussian noise map of unitary variance. We observe y = x + σz, where σ 2 is the noise variance. A combined filtering algorithm has been developed by Starck et al. [START_REF] Starck | Very high quality image restoration by combining wavelets and curvelets[END_REF], estimating x as x as described below, using K dictionaries Φ k , each having a fast transform and reconstruction. We denote by Φ T k the forward transform. Let α k = α k,i be the coefficient list indexed by i obtained by applying the k-th transform α k = Φ T k y. Let Ω y k be the list of coefficients of y in Φ k which have an absolute value greater than κσ : Ω y k = i : |α k,i | > κσ , κσ being the final threshold level (usually 3σ). The Combined Filtering method consists in minimizing

min x K k=1 Φ T k x ℓ 1 s.t. ∀k ∈ [1, K] : Φ T k x -Φ T k y Ω k ℓ ∞ ≤ σ 2 , (24)
where (α k ) Ω k denotes the restriction to atoms in Ω k . The algorithm proposed in [START_REF] Starck | Very high quality image restoration by combining wavelets and curvelets[END_REF] solves this minimization problem by iteratively transforming the solution in each dictionary, comparing the coefficients with the noisy data's coefficients, and soft-thresholding with a decreasing level.

The main drawback of this algorithm is that it requires to have in memory all coefficients of all transforms during the iterations in order to evaluate the right-hand part of the functional. With the size of 3D data and the redundancy of the proposed transforms, this requires a huge memory size. Hence, using this method, we lose the advantages of the low-memory implementation that make our new transforms useable on a laptop computer with very low memory consumption. Therefore, we propose an alternative approach for the Combined Filtering.

The idea is that the residual between the data and the filtered cube must not be sparse in any of the chosen transforms, and thus leads to small coefficients. It can be interpreted as "There is no structure (in the sense of detection with Φ k ) in the residual". We aim to solve the optimization problem

min x y -x 2 ℓ 2 s.t. ∀k ∈ [1, K] : Φ T k (y -x) ℓ ∞ < κσ (25)
where σ is the noise estimated on the data (or the energy of the faintest detail we wish to recover). This problem can be solved using Algorithm 6. With this new algorithm, the memory required Algorithm 6: The Combined Filtering Data: The observed data y, the noise level σ, an initial threshold level k (0) σ, and a number of iterations N. Result: The estimate x (N) of the unknown x. Let ST λ be the soft-thresholding operator with threshold λ. Let Φ k , k ∈ [0, K) be the K dictionaries to combine. Let κ = 3 be the final threshold level (when multiplied by σ). begin

x (0) = 0 // Initial estimate for n = 0 to N -1 do u (0) = x (n) λ (n) = (k (0) σ -κσ)(N -1 -n)/(N -1) + κσ // Current threshold level for k = 0 to K -1 do u (k+1) = u (k) + Φ k ST λ (n) Φ T k (y -u (k) ) // New estimate x (n+1) = u (K) end
to treat data of size 128 3 with wavelets, BeamCurvelets and RidCurvelets with overlapping is about 110MB, while it would need about 4GB with the old algorithm if it were optimally coded (not using transforms as black boxes, but interacting with them). In our current implementation, it uses twice that.

We used this Combined Filtering method to denoise the synthetic data in Figure 10, and compared to the single-transform denoising PSNRs (see Table 2), we get a PSNR of 45.4dB, which is almost 2dB higher than the best transform alone. Figure 12 shows the result of the new CFA.

λCDM denoising

In astrophysics, we use N-body simulations when making an initial homogeneous universe evolve with gravity to form structures, and during this process, there is a noise present in the low density areas. Simulations of Cold Dark Matter (λCDM) exhibit formation of clusters and filamentary, the density of the filaments being a thousand times lower than the clusters. As the noise is important at low density, because of material particles of a given mass, we have to apply a denoising between the time iterations of the simulation in order to improve it. In the following example, we show a denoising of a simulated λCDM cube obtained by a simulation using the RAMSES code [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. a new high resolution code called ramses[END_REF] with cosmological parameters as follows : Ω m = 0.3 (i.e. 30% of standard matter), Ω λ = 0.7 (70% of Dark Matter, which is interacting gravitationally only), the Hubble parameter h = 0.7 (the rate of expansion of the universe) and σ 8 = 0.9 (the dark matter density fluctuations at present time at scale 8MPc). The noise added to the data has a standard deviation (σ = 0.025) comparable to the filaments amplitude (0.02), which is very small compared to the total amplitude of the data (≈ 100), as the most dense clusters are five orders of magnitude higher than the filaments. We first tried denoising the data using the undecimated wavelet transform, which correctly recovers the high-level clusters, but it fails to separate the filamentary structure from the noise. Facing this problem, we used the BeamCurvelet transform instead. The result, shown in figure 13, although recovering the filaments, reveals strong artifacts surrounding the high-density clusters, and therefore is useless. This phenomenon is not due to the noise, but to the form of the BeamCurvelet atoms. Fine-scale BeamCurvelet atoms have an oscillating behavior, and although these oscillations are very low (order 1e -4 compared to the atom's amplitude), they become a serious problem when structures contained in the data have a huge dynamic range, as in our λCDM simulations where clusters are 10 5 larger than filaments. When the atom is used to approximate a high density object, and when we look round it at scales 10 5 smaller, the oscillations become dominant. In order to solve this problem, we must use conjointly the wavelets and the BeamCurvelets combined in the algorithm we presented in section 4.3.3.

The results of denoising the λCDM data using the new combined filtering algotithm (Alg. 6) are shown in figure 14. As expected, the result is far better than using one transform alone, and we restore the filaments very well. The oscillating structure has also naturally completely disappeared from the residual.

3D Inpainting

Introduction

When measuring quantities on sky observations, we often have to apply a mask to the data to cover inappropriate data, due to the stacking on the line of sight. For example, on the full-sky microwaves images, we have to mask out the whole galactic plane (roughly a cylinder on the sphere) because it is very polluted by the foreground, the Milky Way. In the case of 3D data, there can be a whole cone missing, a sphere of a given radius around the center, or any local area hidden by something. This problem has already been addressed with sparse representations in 2D, for example in weak lensing data [START_REF] Pires | FASTLens (FAst STatistics for weak Lensing) : Fast method for Weak Lensing Statistics and map making[END_REF]. We investigate, in this section, how the 2D sparse inpainting methods can be applied in 3D, using the new 3D Curvelet transforms.

The Inpainting problem

Inpainting is the process of recovering missing parts in altered data from the still available clues. Let x be our three-dimensional data cube with missing data indicated by the mask M. The available data is y = Mx. Given a dictionary Φ, we are trying to recover x from the observed y and the mask M. This is an ill posed inverse problem. To get a consistent solution, one must seek regularized solutions. One such regularization is to suppose that x is sparse in one dictionary of atoms Φ, which means that x can be represented by a few atoms from Φ. Therefore, we want to solve the following inpainting problem :

arg min x Φ T x 0 s.t. y -Mx 2 < ǫ(σ), (26)
where Φ + stands for the analysis operator (forward transform). To solve this problem, Algorithm 7 is applied (see [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis[END_REF]).

Example 1 : Random missing voxels

We tried to push the degradation process to see how the BeamCurvelet performed. We applied different masks to the λCDM data, with an increasing percentage of randomly missing voxels. The data are strictly positive (matter density); to avoid high dynamic range issues, in the whole inpainting section we work on the logarithm of the data, but for more clarity, we display all figures with a linear scale, truncated at about 10 times the level of the filaments; the strong clusters are thus saturated. Figure 15 Algorithm 7: The Inpainting Process Data: A mask M, the observed data x = My, y being the unknown cube , the number of iterations N, an initial threshold level k (0) . Result: The estimate x (N) of y. Let HT λ (n) be the hard-thresholding operator with threshold λ (n) . MAD stands for Median absolute deviation. begin x (0) = y ∀i ∈ [0, N) : inpainting process and the transform used (the BeamCurvelet transform) are very well adapted to recover this kind of structure. Animations of the inpainting process and 3D representations of the results of this experiment and the following ones can be seen online at http://arnaud.woiselle.fr in the research pannel.

k (i) = k (0) (N -1 -i)/(N -1) for n = 0 to N -1 do λ (n) = k (n) • MAD Φ T x (n) x (n+1) = Φ HT λ (n) Φ T x (n) + M y -Mx (n)

Example 2 : Missing blocks

In order to see the gain when using one transform rather than another in inpainting when the missing data is spatially coherent, we applied a mask to our λCDM cube, with 50% randomly missing voxels, and an additional 3D checkerboard array of missing boxes of 343 contiguous voxels each (7 × 7 × 7). We applied the inpainting algorithm to restore the data using either the RidCurvelet transform, or the BeamCurvelet transform. The results are shown in Figure 16. The reconstruction is very good in both cases, but as the BeamCurvelets are more adapted to the structure of the data than the RidCurvelets, their reconstruction is better : the L 2 and L 1 norms of the errors are lower with the BeamCurvelts, and there is visually less residual noise in the reconstruction, which can be understood as the atoms (planes) of the RidCurvelet transform must spread around the data's filaments more than the atoms (filaments) of the BeamCurvelet transform to approximate the data.

Conclusion

We have presented two novel multiscale geometric decompositions of a three dimensional volume. The RidCurvelet represents well surfaces and the BeamCurvelet is well adapted for 22 analyzing filamentary features. Preliminary results have shown that they produce interesting results for denoising and inpainting. An important aspect we did not treat is the noise property for denoising applications. Here we have considered only Gaussian noise. However, most 3D astronomical data sets require that we consider Poisson noise. A solution would be to consider non linear 3D decompositions, using both BeamCurvelet and the Multiscale Variance Stabilization, in a similar way to what has been done for the 2D curvelet transform Poisson denoising [START_REF] Zhang | Wavelets, ridgelets, and curvelets for Poisson noise removal[END_REF].

Aknowledgements

We wish to thank Romain Teyssier for providing us the λCDM simulated data [START_REF] Teyssier | Cosmological hydrodynamics with adaptive mesh refinement. a new high resolution code called ramses[END_REF] used in sections 5 and 4.

This work is supported by Sagem DS.

Figure 1 :

 1 Figure 1: Scheme showing the main steps of the ridgelet transform in 2D : lines passing through the origin are extracted from the Fourier transform of the image; a wavelet transform is applied on each of these lines.

FFT2DF

Figure 2 :

 2 Figure 2: Scheme showing the main steps of the curvelet-99 transform : the ridgelet transform is applied to each block of each scale of an isotropic 2D wavelet transform.

Figure 3 Figure 3 :

 33 Figure 3: Schematic view of a 3D Beamlet transform. At a given direction, sum over the (θ 1 , θ 2 ) line to get a • point.Repeat over all its parallels to get the dark plane and apply a 2D wavelet transform within that plane. Repeat for all the directions to get the 3D Beamlet transform. See the text (section 3) for a detailed explanation and implementation clues.

Figure 4 :

 4 Figure 4: Global flow graph of a 3D BeamCurvelet transform.

Figure 5 :

 5 Figure 5: Examples of a BeamCurvelet atoms at different scales and orientations. These are 3D density plots : the values near zero are transparent, and the opacity grows with the absolute value of the voxels. Positive values are red/yellow, and negative values are blue/purple. The right map is a slice of a cube containing these three atoms in the same position as on the left. The top left atom has an arbitrary direction, the bottom left is in the slice, and the right one is normal to the slice.

Figure 6 :

 6 Figure 6: Overview of the 3D Ridgelet transform. At a given direction, sum over the normal plane to get a • point.Repeat over all its parallels to get the (θ 1 , θ 2 ) line and apply a 1D wavelet transform on it. Repeat for all the directions to get the 3D Ridgelet transform. See the text (section 3) for a detailed explanation and implementation clues.

Figure 7 :

 7 Figure 7: Global flow graph of a 3D RidCurvelet transform.

Figure 8 :

 8 Figure 8: Examples of RidCurvelet atoms at different scales and orientation. The rendering is similar to that of figure 5. The right plot is a slice from a cube containing the three atoms shown here.

Figure 9 :

 9 Figure 9: Top row : a plane, and the reconstruction after thresholding the RidCurvelet coefficients at 8.2 (see table1). Bottom row : a filament, and its reconstruction after the thresholding of the BeamCurvelet coefficients at 1.31.

  Figure 9: Top row : a plane, and the reconstruction after thresholding the RidCurvelet coefficients at 8.2 (see table1). Bottom row : a filament, and its reconstruction after the thresholding of the BeamCurvelet coefficients at 1.31.

Figure 10 :

 10 Figure 10: From left to right : a 3D view of the cube containing pieces of shells and a spring-shaped filament, a slice of the previous cube, and finally a slice from the noisy cube. Wavelets RidCurvelets BeamCurvelets Shells & spring 40.4dB 40.3dB 43.7dB

Figure 11 :

 11 Figure 11: From left to right : a slice from the filtered test-cube described in section 4.3.2 (orignial in figure 10) by the wavelet transform (isotropic undecimated), the RidCurvelets and the BeamCurvelets.

Figure 12 :

 12 Figure 12: A slice and a 3D view of the noisy cube on Figure 10 after denoising with the presented Combined Filtering.

Figure 13 :

 13 Figure 13: The central slice of (a) the original λCDM data cube, (b) the noisy data, (c) the data recovered by Wavelet Hard Thresholding, and (d) recovered by BeamCurvelet thresholding.

Figure 14 :

 14 Figure 14: Denoised λCDM data using the Combined Filtering (see Algorithm 6).

=Φ

  HT λ (n) Φ T (I -M)x(n)

Figure 15 :

 15 Figure 15: Central slice of the masked λCDM data with 20, 50, and 80% missing voxels, and the inpainted maps. The missing voxels are dark red.

Figure 16 :

 16 Figure 16: First row : original central frame of the λCDM data cube, and degraded version with missing voxels in red. Bottom row : the filtered results using the RidCurvelets (left) and the BeamCurvelets (right). Missing voxels are dark red.

  instead of the Beamlets. A three-dimensional ridge function Algorithm 2: The Inverse BeamCurvelet Transform Put the obtained Fourier plane to the block, such that the plane passes through the origin of the block with normal angle (θ 1 , θ 2 ) -Apply a 3D IFFT -Add the block to the wavelet scale at the position (B x , B y , B z ), using a weighted function if overlapping is involved if the scale number is even then according to the parabolic scaling :

Data: An initial block size B, and the BeamCurvelet coefficients : series of wavelet-space planes indexed by a scale, angles (θ 1 , θ 2 ), and a 3D position (B x , B y , B z ) Result: The reconstructed data cube begin for all scales from the finest to the second coarsest do Create a 3D cube the size of the current scale (according to the 3D wavelets used in the forward transform) for each block position (B x , B y , B z ) do Create a block B of size B × B × B for each plane (θ 1 , θ 2 ) indexed with this position do -Apply an inverse 2D wavelet transform -Apply a 2D FFT -B = 2B (in the undecimated wavelet case) B = B/2 (in the pyramidal wavelet case)

  2 s × 2B s = 6B 3 s . There are (ρN s /B s ) 3 blocks in

	Original
	datacube
	3D	(θ1, θ2)
	Ridgelet
	transform
	position	wavelet scales
	directions	directions
	(θ1, θ2)	(θ1, θ2)
		Wavelet transform

  Put the obtained Fourier line to the block, such that the line passes through the origin of the block with the angle (θ 1 , θ 2 ) -Apply a 3D IFFT -Add the block to the wavelet scale at the position (B x , B y , B z ), using a weighted function if overlapping is involved if the scale number is even then

	Data: A data cube x and an initial block size B Result: RidCurvelet transform of x begin Apply a 3D isotropic wavelet transform for all scales from the finest to the second coarsest do Cut the scale into small cubes of size B for each block do Apply a 3D FFT Extract lines passing through the origin at every angle (θ 1 , θ 2 ) for each line (θ 1 , θ 2 ) do apply an inverse 1D FFT apply a 1D wavelet transform to get the RidCurvelet coefficients if the scale number is even then according to the parabolic scaling : B = 2B (in the undecimated wavelet case) B = B/2 (in the pyramidal wavelet case) according to the parabolic scaling : B = 2B (in the undecimated wavelet case) B = B/2 (in the pyramidal wavelet case) Apply a 3D inverse isotropic wavelet transform , which corresponds respectively Algorithm 3: The RidCurvelet Transform end

end Algorithm 4: The Inverse RidCurvelet Transform Data: An initial block size B, and the RidCurvelet coefficients : series of wavelet-space lines indexed by a scale, angles (θ 1 , θ 2 ), and a 3D position (B x , B y , B z ) Result: The reconstructed data cube begin for all scales from the finest to the second coarsest do Create a 3D cube the size of the current scale (according to the 3D wavelets used in the forward transform) for each block position (B x , B y , B z ) do Create a block B of size B × B × B for each line (θ 1 , θ 2 ) indexed with this position do -Apply an inverse 1D wavelet transform -Apply a 1D FFT -

Algorithm 5 :

 5 Curvelet denoising Algorithm Data: The observed data, the noise level σ and a threshold level κσ. Separate the following loop on multiple CPUs for each block (s, B x , B y , B z ) of all scales except the coarsest do Extract the block +1 block Apply a 3D Beamlet or Ridgelet transform to the block +1 transformed block Hard Threshold the coefficients at κσ Apply a 3D Beamlet or Ridgelet inverse transform to the block Put the block into the new wavelet scale, or add it with a weighing window if overlapping is involved Apply an inverse 3D wavelet transform end

	Result: The denoised data.			// Memory usage
	begin
	Apply a 3D isotropic Meyer wavelet transform		1 Wavelet transform
	Duplicate it and set it to zero		+1 Wavelet transform
		RidCurvelets BeamCurvelets Undecimated Wavelets Decimated Wavelets
	Plane	11.78	8.11	5.06	3.13
	Filament	1.02	2.22	1.30	1.16

Table 1 :

 1 Maximum value of the coefficients of a transformed cube containing either a plane or a filament, using a given representation.


Keywords: Data Structure Generalized N3-tree (octree), Space partitioning, Empty space compaction

Ideas [x] Rendering only dependent on what is visible ray-tracing approach [x] Load only needed data, at the needed resolution occlusion + LOD (level of details) ray-guided streaming [x] Reuse loaded data as much as possible cache mechanism on GPU [x] Minimize computation, minimize memory transfers Tree Data Management (space partitioning) to store and organize data (octree or generalized N3-tree, + SDK example kd-tree) Cache System on GPU : LRU mechanismn (least recently used) (to get temporal coherency) Data Production Management : on host, GPU, or hybrid mode Goal : produced data are kept in cache on GPU Visit algorithm : traverse your data (loaded in cache) as could be done for rendering Renderer (hierarchical volume ray-casting, cone tracing, emission of requests, brick marching)

File format

By level of resolution and by data types 1 -nodes (spatial structure) 2 -bricks of voxels no header, no meta-file…

  mechanism that continuously add details :  Re-use noise computation at previous resolution level noise : sum of octaves (frequencies) compute only one octave by level of resolution store them in GPU cache reproduce the sum : F(level N) = currentNoise(level N) + previousNoise(level N-1) Finally : replace "noise computation" by "data fetch in cache" Hypertexture (shape + solid texturing) • F(x)==0 : outside the object • F(x)==1 : strictly inside the object • 0<F(x)<1 : inside a fuzzy region Idea : add details on armor soldier • Distance function : sphere • Thickness : fuzzy region (Perlin noise) • Current work : optimizations -Stop adding noise components if density is outside traditional CG rasterized scenes (or compositing) -renderer takes color and depth buffer as input and updates them with voxels OpenGL interoperability Voxel-based Instancing Goal (not yet, in study/progress) Have several GigaVoxels entities (forest) -Need proxy geometry -Render projected 2D BBox -Use GLSL fragment shader -(rays start + direction) Instancing Goal Environment to study, test, profile, debug… Kind of "generic" API : browsers, editors, etc… Tools : transfer function, performance monitor, etc… 3D models importer (+ GLSL shaders editor) Pipeline » : dynamic library (.dll, .so) => kind of « effect » Idea : plug this "node" in a scene graph à la OpenSceneGraph (see the "ppu" project) or Ogre3DNot yet finished : it has been "polluted" with Qt to have USER custom editors. inline PTX assembly in CUDA. // Target ISA Notes : %clock64 requires sm_20 or higher. // // The constraint on output is "=l" : // -'=' modifier specified that the register is written to // -'l' modifier refers to the register type ".u64 reg" asm volatile ("mov.u64 %0,%clock64;" : "=l"(result) : : ); return result; } Record device function for Rendering kernel (1 thread per pixel), then compute average time on all pixels French ANR Research Project RSA Cosmos : planetariums [2010-2014] "Real-Time and Interactive Galaxy for Edutainment" I. Visualize static galaxies II. Data amplification (procedural generation) III. Animate galaxies Optimisations use OpenGL billboards to stop rays (VBO) multi-GPUs ? time budget ?, priority on bricks ? RTIGE PART 3 -The library -[ 1 ] -Core library -[ 2 ] -User custom API The library Goal -Capitalize on the knowledge of the team -Continue research on the subject -Hide the complexity of the underlying « core » library (cache mechanism, data structure management, etc…) -Give USER access to customize the « data production » and the Cuda 4.x and Compute Capability SM at least 2.0 [ 1 ] -Common API : « hidden » generic mechanisms Tree Data Management (space partitioning) to store and organize data (octree or generalized N3-tree, + SDK example kd-tree) Cache System on GPU : LRU mechanismn (least recently used) (to get temporal coherency) Data Production Management : on host, GPU, or hybrid mode Goal : produced data are kept in cache on GPU Visit algorithm : traverse your data (loaded in cache) as could be done for rendering Renderer (hierarchical volume ray-casting, cone tracing, emission of requests,