Antoine Guerin

Jean-Pierre Rossetti

Didier Hantz
email: dider.hantz@ujf-grenoble.fr

Michel Jaboyedoff
email: michel.jaboyedoff@unil.ch

Estimating rock fall frequency in a limestone cliff using LIDAR measurements

Keywords: rock fall, frequency, hazard, risk, LiDAR, terrestrial laser scanner

Terrestrial Laser Scanner has been used to detect rock falls which have occurred in a limestone cliff during some years, in the difficult configuration of the Subalpine Chains. In a rock wall of width 750 m and height 200 m, 130 rock falls larger than 0.1 m 3 have been detected for a period of 1180 days. The distribution of the rock fall volumes is well fitted by a power law, with an exponent which is compatible with the exponent found for the 120 km long cliffs of the Grenoble area. But the spatialtemporal frequencies given by the two analyses are very different. The number of rock falls larger than 1 m 3 , which occur per century and per hm 2 , is about 150 times larger for the bedded limestone of Sequanian stage than for the massive limestones of Tithonian and Barremian stages.

Introduction

Estimating rock fall frequency is needed to characterize a diffuse rock fall hazard [START_REF] Picarelli | Hazard characterization and quantification[END_REF][START_REF] Fell | A framework for landslide risk assessment and management[END_REF][START_REF] Hantz | Quantitative assessment of diffuse rock fall hazard along a cliff foot[END_REF]. Up to now this frequency is determined from historical inventories. The minimal volume detected in these inventories can be relatively small when the rock blocks fall on a road or railway from a cut slope, but it is larger when they fall from a high rock cliff [START_REF] Dussauge-Peisser | Probabilistic approach to rock fall hazard assessment: potential of historical data analysis[END_REF]. In the last years, terrestrial laser scanner (TLS) has been used to detect rock falls in coastal cliffs [START_REF] Rosser | Terrestrial laser scanning for monitoring the process of hard rock coastal cliff erosion[END_REF][START_REF] Dewez | Five seasons of chalk cliff face erosion monitored by terrestrial laser scanner: from quantitative description to rock fall probabilistic hazard assessment[END_REF] and in high mountain walls [START_REF] Ravanel | Quatre années de suivi de la morphodynamique des parois rocheuses du massif du Mont Blanc par laserscanning terrestre. Images et modèles 3D en milieux naturels[END_REF], where the rock fall frequency is relatively high. We have used this method to detect rock falls in a limestone cliff of the Subalpine Ranges, which threatens the Grenoble town.

Description of the cliff and measurements

The Mont Saint-Eynard (1308 m) is located 4 km to the North of the Grenoble centre and towers above a residential area of the town (Figures 1 and2). Its South-East face consists of, from top to bottom: a 120 m high limestone cliff (Tithonian and upper Kimmeridgian stages); a 100 m high forested slope of marl and marly limestone (Kimmeridgian stage); a 240 m high limestone cliff (Sequanian stage); a 300 m high forested talus slope, covering marl and marly limestone of the Orfordian stage. This paper describes the results obtained for the Sequanian cliff.

The measurement place was located at the foot of the talus slope, on a protection embankment at an elevation of 580 m. The inclined distance to the cliffs ranges between 625 m and 900 m. Note that no place was found closer to the cliffs. Photographs and laser measurements were carried out in August 2009 and November 2012.

The laser scanner technology also called LiDAR (Light Detection And Ranging), is based on the acquisition of a point cloud using a time-of-flight distance measurement of an infrared laser pulse which reflects on the topography. The raw data consist of the x, y, z coordinates of each reflection point and the intensity of the reflected pulse. The y axis corresponds to the outward axis of the LiDAR camera and x and z axes are parallel to the sides of the LiDAR scene (x is roughly horizontal). We have used two Optech systems: ILRIS-3D in 2009 and ILRIS-LR in 2012. The main characteristics of these systems are given in Table 1. It can be seen that ILRIS-LR has a higher repetition rate, allowing a greater number of points to be measured for a given period of time. 7 mm @ 100 m 7 mm @ 100 m Raw angular accuracy 8 mm @ 100 m 8 mm @ 100 m Field of view 40° x 40° 40° x 40° Minimum step size 0,001146° (20 µrad) 0,001146° (20 µrad) Maximum density 2 cm @ 1000 m 2 cm @ 1000 m Rotational speed 0,001 to 20°/sec 0,001 to 20°/sec Beam diameter 22 mm @ 100 m 27 mm @ 100 m Beam divergence 0,009740° (170 µrad) 0,014324° (250 µrad) Laser wavelength 1535 nm 1064 nm Integrated camera 3,1 MP 3,1 MP Table 1. Main characteristics of the LiDAR used.

Parameter

ILRIS-3D ILRIS-LR

The point spacing was about 20-30 cm in 2009 and about 10-13 cm in 2012. According to the distances given above and the accuracy given in Table 1, the expected accuracy of our measurements ranges from about 5 cm for the closest points to 7,5 cm for the farthest ones. Two scans were taken to cover a cliff width of about 750 m.

Data analysis

The softwares 3DReshaper Application and Cloud Compare has been used to process the point clouds.

Cleaning the raw point cloud

Vegetation has a lower reflectance than the rock making up the cliff. Thus, a reflectance threshold has been chosen to remove most of the points corresponding to vegetation.

Georeferencing

Georeferencing the LiDAR point clouds was made by registering these with a DEM (1 m x 1 m) using Lambert 2 extended (x,y) coordinates and NGF IGN 69 leveling. The DEM is shown in Figure 3. Then the coordinate system has been changed in order to easily determine the width and the thickness of the fallen compartments. The width is defined horizontally, parallel to the cliff (x direction), the thickness is defined horizontally, perpendicular to the cliff direction (y direction) and the height is parallel to the z axis. The positive direction is inside the cliff for the y axis and towards the East side for the x axis.

Meshing the 2012 point clouds

The more recent point clouds ( 2012) have been transformed in two meshes (polyhedrons), made up of 1.6 million (left scene) and 1.1 million (right scene) of triangles (Figures 4).

The average distances between the vertex of the polyhedrons are respectively 26 cm and 36 cm, and the vertex numbers 800 000 and 570 000. Note that these numbers are about ten times less than the initial numbers of points of the clouds. This reduction is associated to the noise reduction process we used in 3DReshaper and is also necessary for numerical reasons. Registration of the 2012 point cloud with the corresponding mesh allows estimating the roughness of the rock surface at the scale of the triangles making up the mesh. It appears that about 50 % of the points are closer than 1 cm from the mesh, 90 % are closer than 3 cm and 99 % are closer than 7 cm. Note that, in the later case, most of the deviations larger than 7 cm are located on vegetation areas.

Registration of the point clouds

As georeferencing with the DEM was not precise enough, the 2012 meshes and the point clouds acquired in 2009 have been registered (fitted) together in order to put them in the same coordinate system. Ideally, the deviations between these objects should be due only

Detection of the rock falls

A rock fall is considered certain when a positive deviation is observed and the comparison of the 2009 and 2012 photos shows that a rock fall has occurred between these dates. Such a comparison is shown in Figure 5. When a positive deviation is observed and the comparison of the 2009 and 2012 photos shows that no rock fall has occurred, the deviation is considered to be a false rock fall. False rock falls are due to the conditions (c) or (d) expressed in the later paragraph. When a positive deviation is observed and the photo comparison can't conclude if a rock fall has occurred or not, the deviation is considered to be an improbable rock fall if the conditions (c) or (d) occur or a probable rock fall if these conditions don't occur. This situation occurs more and more when the extent of the deviation zone decreases. The reason is that small rock falls are difficult to observe on photographs. Figure 6 shows the proportions of certain, probable, improbable and false rock falls as a function of the rock fall volume. [START_REF] Hantz | An historical, geomechanical and probabilistic approach to rock-fall hazard assessment[END_REF] analyzed the cumulative distribution of rock fall volumes between 10 2 and 10 7 m 3 , occurred in the 120 km long limestone cliffs of the Grenoble area, which include the Mont Saint-Eynard cliff. They found that a power law well describes the distribution, with an exponent absolute value of 0.55 ± 0.11. The values obtained for the Grenoble area and for the Mont Saint-Eynard appear to be compatible with a value of 0.6 ± 0.11. Concerning the spatial-temporal rock fall frequency, the number of rock falls larger than 1 m 3 , which occur per century and per hm 2 , is about 75 for the Mont Saint-Eynard and 0,5 for the Grenoble area. These frequencies correspond to different rock masses. The cliffs of the Grenoble area consist mainly of massive limestone of Tithonian and Barremian stages, whereas the studied cliff consists of bedded limestone of Sequanian stage. Although the latter is included in the Grenoble area cliffs, it represents a tiny part of it.

Size-frequency relation

Discussion

Conclusion

Terrestrial Laser Scanner can be used to detect rock falls which have occurred in high rock walls during some years, in the difficult configuration of Subalpine Chains. In a rock wall of width 750 m and height 200 m, 130 rock falls larger than 0.1 m 3 have been detected for a period of 1180 days. The distribution of the rock fall volumes is well fitted by a power law, with an exponent (b = 0.65) which is compatible with the exponent found for the 120 km long cliffs of the Grenoble area (b = 0.55 ± 0.11). But the spatial-temporal frequencies given by the two analyses are very different. The number of rock falls larger than 1 m 3 , which occur per century and per hm 2 , is about 150 times larger for the bedded limestone of Sequanian stage than for the massive limestone of Tithonian and Barremian stages.

Figure 1 .

 1 Figure 1. Localisation of the Mont Saint-Eynard.

  After cleaning, the two point clouds consisted of 1.5 Mpt (left scene) and 1.2 Mpt (right scene) in 2009 and 7.1 (left) and 5.8 Mpt (right) in 2012. They are shown in Figure 2. The average distance between the points was 21-29 cm in 2009 and 10-13 cm in 2012.

Figure 2 .

 2 Figure 2. Point clouds measured in 2009 (high) and 2012 (low).

Figure 3 .

 3 Figure 3. Digital Elevation Model used for georeferencing.

Figure 4a .

 4a Figure 4a. Mesh for the left scene 2012 and rock fall detected.

Figure 4b .

 4b Figure 4b. Mesh for the right scene 2012 and rock fall detected.

  to rock falls occurred between 2009 and 2012. But in reality, there are other causes of deviations: (a) measurement accuracy; (b) the 2009 measurement points do not correspond to the 2012 ones and are not exactly on the triangles defined by the 2012 vertex (due to the curvature and the roughness of the rock surface); (c) the later cause is accentuated in areas where the triangles are large (this situation occurs particularly near the limits of the mesh); (d) vegetation element which has not been removed; (e) earth slide due to the impact of an overlying rock fall.

Figure 5 .

 5 Figure 5. Rock fall proven by comparison of 2009 and 2012 photographs (81 m 3 ).

Figure 6 .

 6 Figure 6. Proportions of certain, probable, improbable and false rock falls as a function of the rock fall volume.

Figure 7

 7 Figure 7 shows the cumulative distribution of the rock fall volumes larger than 0.1 m 3 . It appears that the distribution is well fitted by a power law: b aV N

Figure 7 .

 7 Figure 7. Cumulative distribution of the rock fall volumes.


A Woiselle

J-L Starck

J Fadili

3-D Data Denoising and Inpainting with the low-redundancy Fast Curvelet Transform

Keywords: 3-D Curvelets, sparsity, denoising, inpainting, Morphological Component Analysis, video deinterlacing

, which makes it more suitable to applications including massive data sets. We report an extensive comparison in denoising with the Curvelab implementation as well as other 3-D multi-scale transforms with and without directional selectivity. The proposed implementation proves to be a very good compromise between redundancy, rapidity and performance. Secondly, we exemplify its usefulness on a variety of applications including denoising, inpainting, video de-interlacing and sparse component separation. The obtained results are good with very simple algorithms and virtually no parameter to tune.

Introduction

Since the introduction of the multiscale point of view in image processing, an important gap has been crossed with the development of directional transforms such as the curvelet transform [START_REF] Candès | Curvelets -a surprisingly effective nonadaptive representation for objects with edges[END_REF].

For instance, the latter provide a non-adaptive (near-)optimally sparse representation of smooth curvilinear structures in images, while their wavelet predecessors obtained as tensor-product of 1-D wavelets fail to efficiently represent these anisotropic structures. In fact, the curvelet transform provides atoms that are well localized in space and frequency and exhibit a strong directional selectivity. This transform has found a wide spectrum of applications including 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]. With the increasing computing power and memory storage capabilities of computers, it becomes feasible to analyze 3-D data as a volume and not only slice-by-slice, which would mistakingly miss the 3-D geometrical nature of the data. In order to process such massive data, there is an urgent need to design new 3-D multiscale directional transforms. In the literature, some 3-D representations have already been developed, in addition to the natural extension of separable wavelets to 3-D. For example, 2-D 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] have been extended to d-D Surflets for compression purposes [START_REF] Chandrasekaran | Surflets : A sparse representation for multidimensional functions containing smooth discontinuities[END_REF][START_REF] Chandrasekaran | Representation and compression of multidimensional piecewise functions using surflets[END_REF]. d-dimensional wedgelets correspond to an adaptive transform based on a 2 d -tree spatial partitioning of the data, where each leaf of the 2 d -tree has two constant regions separated by a smooth polynomial manifold of codimension 1. Surfacelets have been proposed in [START_REF] Lu | 3-D directional filter banks and surfacelets[END_REF] as a non-adaptive transform to efficiently represent singularities along 2-D smooth surfaces. They correspond to a 3-D extension of contourlets, and similarly to contourlets, their implementation relies on a filterbank decomposition of the 3-D data in both scale and angle. They produce a tiling of the 3-D frequency domain quite similar to the one of curvelets although not strictly faithful to the sharp frequency localization of curvelets. The latter were first designed and characterized in a 2-D continuous setting [START_REF] Candès | Curvelets -a surprisingly effective nonadaptive representation for objects with edges[END_REF][START_REF] Candès | New tight frames of curvelets and optimal representations of objects with C 2 singularities[END_REF], then discretized and extended to higher dimensions [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF][START_REF] Candès | Fast discrete curvelet transforms[END_REF]. The 2-D curvelet representation was formally proved to yield a (near) optimal sparse representation of piecewise smooth (C 2 ) images away from singularities along C 2 edges. More generally, the curvelet transform would provide efficient representation of d-D smooth objects with singularities along smooth manifolds of codimension 1. The 3-D implementation of the fast curvelets was proposed in [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF][START_REF] Candès | Fast discrete curvelet transforms[END_REF] with a public code distributed (including the 2-D version) in Curvelab, a C++/Matlab toolbox available at www.curvelet.org. This 3-D fast curvelet transform has found applications mainly in seismic imaging, for instance for denoising [START_REF] Ma | Three-dimensional curvelets for coherent vortex analysis of turbulence[END_REF] and inpainting [START_REF] Herrmann | Non-parametric seismic data recovery with curvelet frames[END_REF]. However, a major drawback of this transform is its high redundancy factor, of approximately 25. As a straightforward and somewhat naive remedy to this problem, the authors in [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF][START_REF] Candès | Fast discrete curvelet transforms[END_REF] suggest to use wavelets at the finest scale instead of curvelets, which indeed reduces the redundancy dramatically to about 5.4 (see Section 2.3 for details). However, this comes at the price of the loss of directional selectivity of fine details. On the practical side, this entails poorer performance in restoration problems compared to the full curvelet version. Note that directional selectivity was one of the main reasons curvelets were built at the first place.

In this paper, we present a new implementation of the fast curvelet transform (FCT) with a special emphasis on the 3-D case that has a reduced redundancy, of only 10, without sacrificing the directional selectivity property at the finest scale. Our implementation provides a very good tradeoff between the full FCT and its less-redundant counterpart using only wavelets at the finest scale to reduce the redundancy, as suggested in [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF][START_REF] Candès | Fast discrete curvelet transforms[END_REF]. In Section 2, we present our new implementation of the FCT, and compare it to the one implemented in Curvelab. In Section 3, we show how this new curvelet construction can be wisely used for many different applications including denoising, inpainting, video deinterlacing and sparse component separation.

Low redundancy 3-D fast curvelets 2.1 Background

The FCT of a 3-D object consists of a low-pass approximation subband, and a family of curvelet subbands carrying the curvelet coefficients indexed by their scale, position and orientation in 3-D. These 3-D FCT coefficients are formed by a proper tiling of the frequency domain following two steps (see Figure 1):

• Cartesian coronization or multiscale separation: first decompose the object into (Cartesian) dyadic coronae based on concentric cubes;

• Angular separation: each corona is separated into anisotropic wedges of trapezoidal shape obeying the so-called parabolic scaling law (to be defined shortly). The 3-D FCT coefficients are obtained by an inverse Fourier transform of applied to each wedge appropriately wrapped to fit into a 3-D rectangular parallelepipeds.

In the rest of this section, we will rehearse the main details of these two steps of the 3-D FCT, and describe the modifications we propose in the implementation to reduce its redundancy.

Symmetry

The FCT isolates, in the frequency domain, oriented and localized 3-D wedges. There is a symmetry on a 3-D Cartesian grid : the cube has six faces which can be processed in a similar way. Let ω = (ω 1 , ω 2 , ω 3 ) ∈ [-1/2, 1/2] 3 be a frequency in the 3-D Shannon band. Exploiting the above symmetry, we will only focus on the subspace

{ω 1 > 0, |ω 2 /ω 1 | < 1, |ω 3 /ω 1 | < 1} which is a squared-based pyramid.
The five other parts can be dealt with exactly in the same way by symmetry around the origin and exchange of axes.

Cartesian coronization

The multiscale separation is achieved using a 3-D Meyer wavelet transform [START_REF] Starck | Image reconstruction by the wavelet transform applied to aperture synthesis[END_REF][START_REF] Mallat | A wavelet tour of signal processing[END_REF], where the Meyer wavelet and scaling functions are defined in Fourier domain with compactly supported Fourier transforms.

Let's denote ψ j as the Meyer wavelet at scale j ∈ {0, • • • , J -1}, and φ J-1 the scaling function at the coarsest scale. Denote M j = ψj = 2 -3j/2 ψ(2 -j •) and M J = φJ-1 = 2 -3(J-1)/2 φ(2 -(J-1) •) their Fourier transforms (see Section 2.2.1 for a formal definition). M J is a lowpass and the wavelet functions {M j } 0≤j<J is a family of bandpass frequency localized windows that form a uniform partition of the unity. Applied to a 3-D object, the family {M j } 0≤j<J separates it into Cartesian coronae (annuli), and M J selects its low frequency content. This coarsest subband is kept unaltered, and after an inverse Fourier transform, provides us with the first curvelet coefficients, which corresponds to coarse scale isotropic atoms. Only the next detail scales, i.e. those corresponding to {M j } 0≤j<J , have to be processed further.

Angular separation Let N a be the number of angles on one edge of one face of the finest scale, for a total of N 2 a angles on each face, and thus six times more bands for the entire considered scale (see Figure 1). This number varies with the scale because the number of angles decreases with the scales becoming coarser to obey the parabolic scaling law of curvelets [START_REF] Candès | Fast discrete curvelet transforms[END_REF], with a parabolic scaling matrix diag(2 -j , 2 -j/2 , 2 -j/2 ) (one short direction and two long ones). This property is essential to ensure that the 3-D curvelets are a basis for sparsely representing smooth trivariate functions with 2-D smooth surface-like singularities.

The vector indexing the angular locations on a face at the j th scale may be expressed with l = (l, l ′ ) ∈ {0, • • • , 2 ⌊-j/2⌋ N a -1} 2 , where ⌊•⌋ is the integer part of its argument 1 . Recall that a wedge is the trapezoidal region sharply localized along a given angle at a given scale, see the dark gray area in Figure 1. The center of the wedge is on the line going from the origin to the point (1, θ l , θ ′ l ′ ), with

θ l = -1 + 2l + 1 2 ⌊-j/2⌋ N a , θ ′ l ′ = -1 + 2l ′ + 1 2 ⌊-j/2⌋ N a . (1)
We can now define the angular separation by multiplying the dyadic annuli corresponding to the wavelet detail subbands by the smooth angular windows V j,l in the Fourier domain. The angular windows are built from a smooth real-valued function V supported on [-1, 1] and satisfying the partition property

∞ l=-∞ V 2 (t -2l) = 1 ∀t ∈ R. (2)
The angular window at scale j and orientation l = (l, l ′ ) is then constructed as

V j,l (ω) = V 2 ⌊-j/2⌋ N a ω 2 -θ l ω 1 ω 1 • V 2 ⌊-j/2⌋ N a ω 3 -θ ′ l ′ ω 1 ω 1 , (3)
where θ l and θ ′ l ′ are defined in [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF]. Note the scaling factor 2 -j/2 as dictated by the parabolic scaling. Owing to (2), the family of angular windows {V j,l } l makes a uniform partition of the dyadic annulus at scale j, i.e.

l V 2 j,l (ω) = 1 . (4)
However, because of the support constraint on V , this relation does not hold for all ω, and a special care should be taken at the corners where only three out of usually four windows overlap. We thus need to redefine them for (4) to hold for any ω. Here is a simple remedy to this problem. Let l a , l b , l c be the indices of the three corner windows, we redefine them on their overlapping domain Ω as

∀ω ∈ Ω, ∀l ∈ {l a , l b , l c } V j,l (ω) ← V j,l (ω) V 2 j,la (ω) + V 2 j,l b (ω) + V 2 j,lc (ω)
.

(

) 5
Piecing all ingredients together, the scale-angular wedge at scale-orientation (j, l) is extracted by the frequency window

W j,l (ω) = M j (ω) • V j,l (ω) , (6)
which is sharply localized near the trapezoid

(ω 1 , ω 2 , ω 3 ) : 2 j+1 < ω 1 < 2 j , ω 2 ω 1 -θ l < 2 ⌊j/2⌋ /N a , ω 3 ω 1 -θ ′ l ′ < 2 ⌊j/2⌋ /N a . (7)
Once a wedge is extracted, an inverse Fourier transform must be applied in order to get the curvelet coefficients at the corresponding scale and orientation. Prior to this, the trapezoidal wedge has to be transformed to a convenient form for which the 3-D FFT algorithm applies. As it can be seen from Figure 1, the wedge can be inscribed inside a 3-D parallelepiped which is ∼ 2 j/2 long on the (ω 2 , ω 3 ) coordinates (i.e. tangentially), and 2 j on ω 1 , i.e. radially. Although this expands the area including the wedge, we can still wrap it inside a rectangular parallelepiped of dimensions ∼ (2 j , 2 j/2 , 2 j/2 ) centered at the origin aligned with the axes of the grid (see Section 2.2.2 for further details about wrapping). With appropriate choice of the size of the rectangular parallelepiped, the data does not overlap with itself after wrapping. With the wrapping trick, an inverse 3-D FFT can be readily applied to the rectangular parallelepiped to obtain the curvelet coefficients at the selected scale and orientation. Algorithm 1 summarizes our implementation of the 3-D FCT and outlines its main steps.

Angular separation:

foreach scale j = 0 to J do foreach orientation l = (l, l ′ ) do Multiply the wavelet cube at scale j with the angular window V j,l in Fourier;

Wrapping: wrap the result in a rectangular parallelepiped centered at the origin of minimal size (2 j 3/8 × 2 j/2+1 /N a × 2 j/2+1 /N a ) ;

Apply a 3-D inverse FFT to the rectangular parallelepiped to collect the curvelet coefficients ;

end Result: 3-D FCT of f .
We are now ready to explain the differences we introduced in the 3-D FCT implementation to reduce its redundancy.

A new implementation

The overall implementation of this new FCT differs from the original one in several points. The main difference lays in the way to apply the Meyer wavelet transform to the data. Other changes have been introduced as described hereafter.

The multiscale separation

The extra redundancy of the curvelets as implemented in Curvelab originates mainly from the way the radial window is implemented, especially at the finest scale. For the sake of clarity, we will first explain the 1-D case for the Meyer wavelets and show the implications in higher dimensions later on.

The Meyer wavelets ψ(ξ) are defined in Fourier domain as follows :

ψ(ξ) =    exp -i2πξ sin( π 2 ν(6|ξ| -1)), if 1/6 < |ξ| ≤ 1/3 exp -i2πξ cos( π 2 ν(3|ξ| -1)), if 1/3 < |ξ| ≤ 2/3 0 elsewhere
, where ν is a smooth function, that goes from 0 to 1 on [0, 1] and satisfies ν(x)

+ ν(1 -x) = 1.
Associated to this wavelet is the Meyer scaling functions defined by

φ(ξ) =    1, if |ξ| ≤ 1/6 cos( π 2 ν(6|ξ| -1)), if 1/6 < |ξ| ≤ 1/3 0 if |ξ| > 1/3
. There is a pair of conjugate mirror filters (h, g) associated to (φ, ψ) whose Fourier transforms ( ĥ, ĝ) can be easily deduced from ( φ, ψ). ĥ and ĝ are thus compactly supported. As a consequence, the Meyer wavelet transform is usually implemented in the Fourier domain by a classical cascade of multiplications by ĥ and ĝ 2 . However, the wavelet at the finest scale is supported

Frequencies -2/3 2/3 1/2 -1/2 -1/4 1/4 0 0 1 | ψ0| | ψ1| | ψ2| | φ2| | ψj| = | ψ(2 j ξ)| | φj| = | φ(2 j ξ)| Frequencies 0 1 -2/3 2/3 1/2 -1/2 -1/4 1/4 0
Figure 2: Meyer scaling and wavelets functions in Fourier domain. In the discrete case, we only have access to the Fourier samples inside the Shannon band [-1/2, 1/2], while the wavelet corresponding to the finest scale (solid red line) exceeds the Shannon frequency band to 2/3. Top: In the Curvelab implementation, the Meyer wavelet basis is periodized in Fourier, so that the exceeding end of the finest scale wavelet is replaced with the mirrored dashed line on the plot. Bottom : In our implementation, the wavelets are shrunk so that they fit in the [-1/2, 1/2] Shannon band, and the decreasing tail of the finest scale wavelet is replaced by a constant (dashed red line) to ensure a uniform partition of the unity. on [-2/3, -1/6[∪]1/6, 2/3], hence exceeding the Shannon band. This necessitates to know signal frequencies that we do not have access to.

As the FCT makes central use of the FFT, it implicitly assumes periodic boundary conditions. Moreover, it is known that computing the wavelet transform of a periodized signal is equivalent to decomposing the signal in a periodic wavelet basis. With this in mind, the exceeding end of the finest scale wavelet is replaced with its mirrored version around the vertical axis at |ξ| = 1/2, as shown in dashed line in the top plot of Figure 2. Consequently, the support of the data to treat is 4/3 larger than the original one, hence boosting the redundancy by a factor (4/3) d in d-D. The left part of Figure 3 shows the implication of periodization in two dimensions for the original curvelets in Fourier domain. The finest scale curvelets have the end of their decreasing window at high frequencies periodized 3 .

In our implementation, we followed another path. First, toward our quest of reducing the redundancy, the supports of the scaling and wavelet functions (hence filters) were shrank by a factor of 4/3. Furthermore, to maintain the uniform partition of unity, which plays an important role for isometry of the transform, following [START_REF] Starck | Image reconstruction by the wavelet transform applied to aperture synthesis[END_REF], the finest scale wavelet is modified by suppressing its decreasing tail so that the wavelet becomes a constant over ] -1/2, -1/4]∪]1/4, 1/2] (see the dashed line at the bottom of Figure 2). The right part of Figure 3 shows the impact of the proposed modifications on the 2-D curvelets in the frequency domain. This strategy and the conclusions carry over to the 3-D case.

This proposed modification of the Meyer wavelets indeed reduces the redundancy of the transform (see Section 2.3). However, this comes at the cost of some changes undergone by the curvelet atoms at the finest scale. First of all, they are less sharply supported in the spatial domain than the original curvelets because of the discontinuity of their Fourier transform, while the decay of the other curvelets remain unchanged. This is one possible reason underlying the slight difference in denoising performance between the original and the modified FCT as shown in Section 2.5. Secondly, they obey a parabolic scaling but with a different constant compared to the curvelets at the other scales.

Ensuring zero-mean subbands

In the original wrapping-based FCT [START_REF] Candès | Fast discrete curvelet transforms[END_REF][START_REF] Demanet | Wave Atoms, and Wave Equations[END_REF], the wedges are wrapped around the origin using a simple modulo operator, which makes every point fit into a well-sized rectangular parallelepiped centered at the origin whose size is designed so as the data does not overlap with itself after wrapping. However, nothing prevents the center of the parallelepiped from receiving a significant non-zero wrapped Fourier coefficient. After an inverse FFT of the wrapped wedge, it is likely to obtain curvelet coefficient subbands with non-zero means. This is obviously unsuitable since curvelet coefficients are expected to represent high frequency content, and typical thresholdingbased processing (e.g. denoising) will be hampered in such a situation. Hopefully, the size and position of the wedges are so that this misleading phenomenon is generally prevented in practice. Nevertheless, this is not guaranteed in general. Therefore, in order to ensure zero-mean curvelet subbands, a straightforward solution is to translate each rectangular parallelepiped where a wedge has to be wrapped in such a way that the center (zero frequency) gets a true zero coefficient, i.e. a point out of the wedge support, and then to wrap the data around the translated box. Doing so, the curvelet subbands are ensured to be zero-mean valued after wrapping. Figure 4 illustrates the difference between the two wrapping strategies in 2-D for the sake of legibility. The technique extends readily to the 3-D case.

Redundancy

In this section, we shall quantify analytically the redundancy reduction induced by our new FCT in any dimension d. Without loss of generality, we assume that the data object is a d-D hypercube of unit side. Let N a be the number of angles along an edge on one face at the finest scale, for a total of N d-1 a orientations on each face. Let N f = 2d be the number of faces of the d-D data hypercube.

The redundancy of the Cartesian coronization or multiscale separation assuming a dyadic frequency tiling is given by

J j=0 1 2 d j , (8)
which is upper-bounded by R w = 2 d 2 d -1 . As explained in Section 2.2.1, an extra-redundancy R add comes into play in the Meyer wavelet transform with the original FCT implementation, while this is avoided by ours:

R add = 4 3 d Original FCT, 1 Proposed FCT. (9)
At the finest scale and on each face, there are N d-1 a wedges, where the size of each of them is

3 8 × 2 N a × ... × 2 N a d-1 times . ( 10
)
The factor 3/8 corresponds to the radial depth of the scale; see Section 2.2.1 and Figure 2. In the other orthogonal directions, a wedge has a size of 1 Na which we double because of overlapping. The redundancy of a face at the finest scale is then

R f = N d-1 a • 3 8 2 N a d-1 = 3 • 2 d-4 . (11)
As it can be seen, the R f redundancy is independent of N a , and is therefore valid at all scales. For a large enough number of scales, it can be reasonably assumed that coarsest (wavelet) scale has the same redundancy as a curvelet subband at the same scale. Consequently, the overall redundancy of the FCT is upper-bounded by (see ( 8))

R = N f • R f • R w • R add = 3d 2 2d-3 2 d -1 • R add . (12)
In the case where wavelets are used instead of curvelets at the finest scale, the redundancy upperbound is changed to

R ′ = (N f • R f • (R w -1) + 1) • R add = 3d 2 d-3 2 d -1 + 1 • R add . ( 13
)
Table 1 compares numerically the redundancy of the original and the proposed FCT in 2-D and 3-D, when wavelets (W) or curvelets (C) are implemented at the finest scale. It may be worth mentioning that for the practitioner, the memory storage requirement of the original FCT (as implemented in Curvelab) is twice larger than the one predicted by the redundancy formula. Indeed, the original curvelets are complex and real curvelets are obtained by hermitian symmetry. This explains the redundancy 404 claimed by [START_REF] Lu | Multidimensional directional filter banks and surfacelets[END_REF] emphasizing the need of lower redundancy curvelets.

Properties

This section enumerates the main properties enjoyed by the FCT according to our new implementation.

• Reduced redundancy: with a reduction factor of (4/3) d compared to the original version. This is one of the distinctive properties of the proposed FCT, and was the main goal underlying our modifications in the first place. For example, our FCT with full curvelets at all scales is (almost) as redundant as the original one with wavelets at the finest scale. In 3-D, redundancy implied by our implementation is 2.5 times lower that the original FCT with curvelets at all scales. In short, our implementation achieves a low redundancy while maintaining the directional selectivity property at the finest scale unlike the original FCT where wavelets are advocated at the finest scale to lower the redundancy [START_REF] Ying | 3D Discrete Curvelet Transform[END_REF][START_REF] Candès | Fast discrete curvelet transforms[END_REF].

• Isometry and fast exact reconstruction: owing to the uniform partition property of the Meyer wavelets, and the coverage property (2) of the angular window, the collection of curvelets in Fourier obtained by multiplication of the scale and angular windows also ensures a uniform partition of the unity. Therefore, with a proper normalization of the FFT (wrapping is a simple reindexing), the proposed FCT corresponds to a Parseval tight frame (PTF), i.e. the frame operator CC * = I, where C * is the FCT analysis operator and C its adjoint. With the PTF property, C turns out to be also the inverse transform operator associated to a fast reconstruction algorithm (each step of the forward transform is easily invertible).

• Parabolic scaling: by construction, the curvelets obey the parabolic scaling law with one short and two long sides ∼ (2 -j , 2 -j/2 , 2 -j/2 ). Although at the finest scale, this property is less faithful to the continuous construction compared to the original FCT (see also the discussion at the end of Section 2.2.1).

• Non-equal ℓ 2 norm atoms: although our FCT implements a PTF, the modified curvelets at the finest scale do not have the same ℓ 2 norm as the curvelets in the other (coarsest) scales. These ℓ 2 norms can nonetheless be calculated analytically so as to normalize the associated curvelets coefficients, which is important for instance in every processing which involves thresholding.

• Guaranteed zero-mean subbands: this is a consequence of the wise translation trick prior to wrapping explained in Section 2.2.2. Of course, this operation preserves isometry and ℓ 2 norms.

Comparison with Curvelab

In the sequel, the new variant of FCT is dubbed LR-FCT, where LR stands low-redundancy. In this section, we compare the 3-D LR-FCT to the original FCT as implemented in the CurveLab toolbox in terms of denoising performance over several datasets, where the noise is additive white Gaussian (AWGN). The goal is to assess experimentally how the LR-FCT version positions itself with respect to the two variants of the CurveLab FCT depending whether wavelets or curvelets are used at the finest scale. For fair comparison, the three transforms were applied in exactly the same conditions. Our goal was not to pick up a very sophisticated thresholding/shrinkage estimator operating in the transform domain (the literature is very rich on the subject), but rather to get a fair picture about differential performance of the transforms themselves with a simple and fast statistical estimator. This is why we opted for hard thresholding with the popular 3σ threshold, where σ is the AWGN standard deviation 5 . The experiments have been carried out on different sorts of data (3-D space simulation, hyperspectral images, still and fast-moving videos), and on ten noise realizations for each triplet transform-dataset-standard deviation. Figure 6 depicts the evolution of the output peak signal-tonoise ratio (PSNR) as a function of the noise standard deviation (translated in PSNR). One can see that the LR-FCT implementation is not as good as the full redundant CurveLab version, but provides a very good tradeoff between redundancy and the denoising performance as measured by PSNR with less than 1dB lower than the CurveLab FCT (with curvelets at all scales) in the least favorable case, and 2-4 dB better than the CurveLab variant using wavelets at the finest scale. As stated in Section 2.2.1, we believe that the difference in performance between the all-curvelets original FCT and the LR-FCT can be explained by the more oscillating behavior of the LR-FCT atoms at the finest scale.

A wide range of applications

This section does not aim to give state-of-the-art results, and instead concentrates on demonstrating how the new transform can be harnessed easily and directly to a variety of problems with reasonable -and beyond -results. Further work will be required to fine-tune these methods and optimize their performance in each specific application. All results reported in this section, especially those involving videos, are available in full resolution online at http://arnaud.woiselle.fr/ awaiting inspection by the interested reader.

LR-FCT denoising : a good tradeoff between efficiency and memory storage

In this experiment, we compare the denoising performance using the LR-FCT to several other 3-D multiscale transforms on various types of datasets (3-D spatial data, hyperspectral images and videos). Again, the noise is AWGN and simple hard thresholding is used. The video datasets included in our experiments were the standard videos used in [] as well as mobile, tempete, and coastguard CIF sequences available at www.cipr.rpi.edu. For hyperspectral data, we used dataset from the OMEGA spectrometer on Mars Express (www.esa.int/marsexpress) with 128 wavelength from 0.93µm to 2.73µm. Beside LR-FCT, the other transforms involved in this comparative study are: the dual-tree complex wavelet transform [START_REF] Kingsbury | Complex wavelets for shift invariant analysis and filtering of signals[END_REF][START_REF] Selesnick | The double-density dual-tree DWT[END_REF], the surfacelet transform [START_REF] Lu | Multidimensional directional filter banks and surfacelets[END_REF] and the orthogonal (decimated) and translation-invariant (undecimated) wavelet transforms. Figures 7 and8 show the output PSNR after denoising as a function of the input PSNR for each transform. Each point on each curve is the average output PSNR on ten noise realizations. The reader may have noticed that the wavelet results are much better here than those tabulated in [START_REF] Lu | Multidimensional directional filter banks and surfacelets[END_REF].

The reason behind this is that unlike those authors, we used the Cohen-Daubechies-Fauveau 7/9 filterbank, which is much better for denoising. Figure8 displays the results for the hyperspectral data from Mars Express (see caption for details). From these experiments, it ca be clearly seen that the proposed LR-FCT compares very favorably to the other multiscale geometrical 3-D transforms, and is particularly better at the low PSNR regime. In a nutshell, it can be safely concluded that LR-FCT provides a very good compromise between denoising performance and memory/CPU requirements.

Dual Tree Wavelets

Inpainting

Inpainting is to restore missing data information based upon the still available (observed) cues from destroyed, occluded or deliberately masked subregions of the data. Inpainting has received considerable interest and excitement and has been attacked using diffusion and transport PDE/Variational principles, non-local exemplar region fill-in and sparsity-based regularization; see e.g. [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] and references therein.

Let f ∈ R N be a vectorized form of the sought after 3-D data cube which is

√ N × √ N × √
N , and M ∈ {0, 1} P ×N , P < N be a binary rectangular matrix where each of its rows is zero except at the entry where a voxel is not missing. The observed (incomplete) data g is then the result of applying the lossy operator M to f : g = Mf + ε . where ε is some noise of finite variance σ 2 that may contaminate the observed values. Restoring f from g is an ill-posed which necessitates some form of regularization to reduce the space of candidate solutions. Here, we promote solutions that are sparse in some prescribed overcomplete dictionary of atoms Φ ∈ R N ×L , L ≥ N , meaning that x := Φα (a synthesis prior) can be sparsely represented to a high accuracy by a few number of atoms in Φ. Put formally, we are seeking to solve the following optimization problem :

min α∈R L α 0 s.t. g -MΦα 2 ≤ ǫ(σ) , (14)
where • 0 is the ℓ 0 pseudo-norm that counts the number of nonzero entries of its argument, and ǫ(σ) is the constraint radius that depends on the noise variance. This is very challenging NP-hard optimization problem, and one has to resort to alternative formulations or greedy algorithms to attempt to solve it. For instance, convex ℓ 1 relaxation could be used instead of the ℓ 0 penalty. Here, we use the algorithm devised in [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis[END_REF] which can be viewed as a stagewise hybridization of matching pursuit with block-coordinate relaxation. The adjective "stagewise" is because their algorithm exploits the fact that the dictionary is structured (union of transforms

Φ = [Φ 1 , • • • , Φ K ])
with associated fast analysis and synthesis operators Φ T k and Φ k ; see [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis[END_REF][START_REF] Starck | Sparse Signal and Image Processing: Wavelets, Curvelets and Morphological Diversity[END_REF] for details. For the reader convenience, Algorithm 2 recalls the main steps of this inpainting algorithm.

In the following experiments, this algorithm was used with a dictionary containing two transforms: the LR-FCT and the undecimated discrete wavelet transform (UDWT), in order to better take into account the morphological diversity of the features contained in the data. Figure 9 Algorithm 2: Inpainting Algorithm.

Data: Observed data g and mask M.

Input: Dictionary Φ = [Φ 1 • • • Φ K ],
number of iterations T iter , final threshold τ (e.g. 3). begin Initial components f

(0) k = 0, k = 1, • • • , K. Initial residual r (0) = g. Initial threshold: let k ⋆ = arg max k Φ T k g ∞ , set λ 0 = max k =k ⋆ Φ T k g ∞ . for t = 1 to T iter do for k = 1 to K do Compute marginal residuals r (t) k = r (t-1) + f (t-1) k .
Update kth component coefficients by thresholding α

(t) k = Thresh λt-1 Φ T k r (t) k . Update kth component f (t) k = Φ k α (t) k . Update the inpainted data f (t) = K k=1 f (t) k .
Update the residuals r (t) = g -Mf (t) . Update the threshold λ t = λ 0 -t (λ 0 -τ σ) /T iter . end Result: The estimate f (Titer) of f .

shows the inpainting result on a synthetic cerebral MRI volume available on BrainWeb [START_REF] Remi | MRI simulation-based evaluation of image-processing and classification methods[END_REF] at http://www.bic.mni.mcgill.ca/brainweb/ with two masks: 80% random missing voxels, and 10% missing z slices. We can see that even with 80% missing voxels, we can still see incredibly faint details in the restored anatomical structures such as in the gyri and the cerebellum.

For video inpainting, we artificially created a mask containing 5% random missing pixels on each frame, a missing line, three clusters of 12 dead pixels each, two consecutive missing frames, and a rotating swirl two-pixel wide. The top row of Figure 10 shows the central frame of the original video, as well as one of the two consecutive missing frames and one frame of the masked video. The inpainting of these missing areas was performed either with a dictionary containing the LR-FCT alone, or another dictionary containing both the LR-FCT and the UDWT. The results are displayed respectively in the middle and bottom rows of Figure 10. While inpainting using the LR-FCT alone or combining it with wavelets give similar PSNR values, visually, inpainting with both transforms yields a smoother result.

Video De-interlacing

Video is often acquired in interlaced format, meaning that every even frame contains only the even lines and similarly for odd frames. The problem of reconstructing the missing lines in the sequence is called de-interlacing and has been addressed with many dedicated methods : using simple spatio-temporal neighborhood linear interpolation [START_REF] Tekalp | Digital Video Processing[END_REF][START_REF] Haan | Deinterlacing: An overview[END_REF], using a directional interpolation to better preserve the edges in the images [START_REF] Doyle | Interlaced to sequential conversion for EDTV applications[END_REF][START_REF] Yoo | Direction-oriented interpolation and its application to de-interlacing[END_REF][START_REF] Oh | Spatio-temporal edge-based median filtering for deinterlacing[END_REF], with more complicated methods relying on motion estimation and compensation [START_REF] Wang | Time-recursive deinterlacing for IDTV and pyramid coding[END_REF][START_REF] Yang | Motion compensation assisted motion adaptive interlaced-to-progressive conversion[END_REF][START_REF] Haan | De-interlacing of video data[END_REF], and recently with variational inpainting [START_REF] Ballester | An inpaintingbased deinterlacing method[END_REF]. In this section, it will be demonstrated that adpating the idea of sparsity-based inpainting as described in Section 3.2, and introducing a proper regularization that takes into account the special structure of the mask, will yield very promising de-interlacing results.

De-interlacing could be naively viewed as a 3-D inpainting problem, where the mask is nonzero at the observed even/odd lines in the consecutive even/odd frames; see Figure 11 for a graphical illustration. However, the major issue that rises with such a mask is that it exhibits a strong coherence with some curvelet atoms 6 . Indeed, curvelets are plate-shaped atoms, and those oriented at 45 degrees are very correlated with the inclined planes formed by the alternating even-odd lines along the time axis. One could think of excluding such atoms from the curvelet dictionary (e.g. by just zeroing their coefficients after the forward FCT in the above inpainting algorithm), or weight DOI [START_REF] Yoo | Direction-oriented interpolation and its application to de-interlacing[END_REF] STELA [START_REF] Oh | Spatio-temporal edge-based median filtering for deinterlacing[END_REF] STIM [START_REF] Ballester | An inpaintingbased deinterlacing method[END_REF] Table 2: PSNR the de-interlacing results on the video sequences stockholm and shields using our method and three other ones from the literature.

the FCT coefficients differently depending on those of the mask. These solutions were tried but the results were not conclusive. Following this discussion, we chose to introduce a regularization term acting directly in the spatial domain. In fact, we implicitly assume that the video to be restored does not contain structures that correlate strongly with the mask, otherwise there is no chance that such structures will be identifiable. Thus, a wise regularization would penalize formation of singularities (edges) in the de-interlaced video that coincide with those of the mask, which is also consistent with the way the correlation between the mask and the curvelets manifests itself. Put formally, the de-interlacing optimization problem to solve writes min

αR T α 0 + µ ∇(Φα) ⊙ ∇f mask 1 s.t. g -MΦα 2 ≤ ǫ(σ), (15)
where f mask is the video representing the mask (Figure 11) vectorized in R N , ∇ is the spatiotemporal discrete gradient operator that produces a vector field in R N ×3 , ⊙ is the Hadamard product, and for u

= (u x , u y , u t ) ∈ R N ×3 , u 1 = i |u[i]| where |u[i]| = u x [i] 2 + u y [i] 2 + u t [i] 2 .
It is not difficult to see that the above penalty can be rewritten as a weighted TV semi-norm:

∇(Φα) ⊙ ∇f mask 1 = W mask (∇(Φα)) 1 = Φα TV,W mask ,
where W mask is a (diagonal) weighting operator whose entries are given by ∇f mask , and which acts on vector fields in R N ×3 by component-wise multiplication by its entries. This yields yet another nice interpretation of our penalty, i.e. promote candidate solutions that are smooth at the locations where the mask changes abruptly, which is exactly the goal we were pursuing.

The resulting optimization problem we end up with is even more challenging that the original inpainting problem. Inspired by the work in [START_REF] Starck | Image decomposition via the combination of sparse representations and a variational approach[END_REF] and [START_REF] Elad | Simultaneous cartoon and texture image inpainting using morphological component analysis[END_REF], we slightly modified Algorithm 2 by incorporating a (sub)gradient-descent correction step with respect to this penalty, namely by adding a step of the form:

f (t) ← f (t) + µγdiv W mask G W mask (∇f (t) )
just after the update of f (t) , where G(.) is a subgradient of the ℓ 1 norm in R N ×3 , and γ is the descent step-size. Although this update is quite heuristic, the obtained algorithm works fairly well in practice. Further work is needed to solve exactly the above problem (or rather its convex relaxation) with theoretical guarantees. We will be left to a future work. Figures 12,[START_REF] Romberg | Multiscale wedgelet image analysis: fast decompositions and modeling[END_REF][START_REF] Chandrasekaran | Surflets : A sparse representation for multidimensional functions containing smooth discontinuities[END_REF][START_REF] Chandrasekaran | Representation and compression of multidimensional piecewise functions using surflets[END_REF] show close-ups of the de-interlacing results on the test sequences stockholm and shields, and compare them the methods proposed in [START_REF] Yoo | Direction-oriented interpolation and its application to de-interlacing[END_REF][START_REF] Ballester | An inpaintingbased deinterlacing method[END_REF][START_REF] Oh | Spatio-temporal edge-based median filtering for deinterlacing[END_REF]. We will refer to them respectively as DOI, STIM, and STELA. The quantitative PSNR results are grouped in Table 2. Although the proposed method is not intended to compete with methods specifically designed for de-interlacing in terms of computational cost, this experiment at least provides a proof of concept that sparsity-based regularization with the LR-FCT transform is efficient for de-interlacing. The results clearly show that the proposed method is at the level of the state-of-the-art both visually and in terms of PSNR. For some videos, it even outperforms the other methods with visually more pleasant quality and less artifacts than dedicated methods (see e.g. Figure 13).

Morphological Component Analysis

The concept of morphological diversity, introduced by Starck et al. [START_REF] Starck | Redundant multiscale transforms and their application for morphological component analysis[END_REF][START_REF] Starck | Image decomposition via the combination of sparse representations and a variational approach[END_REF] is based on the fact that complex signals are linear superposition of more elementary signals (the so-called morphological components). There, such complex signals cannot be efficiently (sparsely) represented with a dictionary built upon a single transform, but rather but rather with an overcomplete dictionary built by amalgamating several sub-dictionaries. Once the sub-dictionaries are identified, solving a sparse decomposition problem in such an overcomplete dictionary yields as a by-product a separation of the structural/morphological components. The Morphological Component Analysis (MCA) algorithm proposed [START_REF] Starck | Redundant multiscale transforms and their application for morphological component analysis[END_REF] was proposed to achieve this goal.

A key to the success of this separation process is a good choice of the sub-dictionaries. Each subdictionary has to be chosen such that its atoms match the shape of the (geometrical) structures they are intended to sparsify, while leading to a non-sparse (or at least not as sparse) representation on the other signal content. For instance, the isotropic wavelet transform [START_REF] Starck | Sparse Signal and Image Processing: Wavelets, Curvelets and Morphological Diversity[END_REF] is a very good candidate for representing isotropic singularities. On the other hand, the curvelet transforms have plate-like atoms well adapted for representing smooth 3-D objects with smooth singularities along smooth surfaces. Using a dictionary built as a union of the 3-D isotropic wavelets and LR-FCT sub-dictionaries, we exemplify the MCA algorithm on a toy 3-D dataset, containing parts of 3-D shells, and 3-D Gaussians. The results in the noiseless and noisy cases are shown in Figures 1617.

Conclusion

We have proposed a new variant of the fast curvelet transform in d-D with a special emphasis on 3-D, lowering its redundancy from 25 to about 11 in 3-D, providing a versatile tool to use in a wide range of inverse problems. We have shown good results in video and spectral data denoising, video and MRI inpainting, sparse component separation, and extended the sparse inpainting approach to the de-interlacing problem. Following the philosophy of reproducible research, a C++/Matlab toolbox implementing the proposed LR-FCT and reproducing the experiments in this paper are made available for download at http://arnaud.woiselle.fr.

Aknowledgements

We wish to thank Michael Elad for constructive discussions. We also wish to thank Luis Garrido for giving us the de-interlaced sequences presented in [START_REF] Ballester | An inpaintingbased deinterlacing method[END_REF] to which we compare in Section 3.3.

2 2 . 5

 25 Low redundancy 3-D fast curvelets 2 2.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2.2 A new implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2.1 The multiscale separation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2.2 Ensuring zero-mean subbands . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.3 Redundancy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.4 Properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 Comparison with Curvelab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3 A wide range of applications 11 3.1 LR-FCT denoising : a good tradeoff between efficiency and memory storage . . . . 11 3.2 Inpainting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.3 Video De-interlacing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.4 Morphological Component Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Algorithm 1 : 1 .

 11 The 3-D Fast Curvelet Transform Data: A 3-D data object f of size N = (N x , N y , N z ). Input: Number of scales J, number of angles N a on each face at the finest scale. begin Multiscale separation: apply the 3-D Meyer wavelet transform in Fourier domain, get cubes of sizes N, N/2, ..., N/2 J .

Figure 1 :

 1 Figure 1: Example in 3-D of how the number of directions N a is defined along an edge, (eight in this example), and showing the overlapping region around one wedge. The dark gray area is the main part of the wedge, and the light one represents the overlapping region.

Figure 2

 2 Figure 2 displays in solid lines the graphs of the Fourier transforms of the Meyer scaling and wavelet functions at three scales.There is a pair of conjugate mirror filters (h, g) associated to (φ, ψ) whose Fourier transforms ( ĥ, ĝ) can be easily deduced from ( φ, ψ). ĥ and ĝ are thus compactly supported. As a consequence, the Meyer wavelet transform is usually implemented in the Fourier domain by a classical cascade of multiplications by ĥ and ĝ 2 . However, the wavelet at the finest scale is supported

Figure 3 :

 3 Figure 3: Top : Examples of 2-D (real) curvelets in Fourier domain at three consecutive scales and different orientations, with the zero frequency at the center. From the outermost edge to the inside, a finest scale curvelet and two lower scale curvelets. (a) According to the Curvelab implementation, and (b) with our modified low-redundancy implementation. Bottom : the corresponding Fourier transforms of 1-D Meyer scaling and wavelet functions.

Figure 4 :

 4 Figure 4: (a) Representation of the influence of a wedge and its overlapping region. The centered dotted rectangle corresponds to the minimal size in which the wedge will be wrapped. (a') The result of the wrapping. (b) The wrapping in a rectangle of the same size, but whose center is chosen such that the zero frequency (big dot) falls outside the support of the wedge. (b') The corresponding wedge after translation and wrapping.

Figure 5 :

 5 Figure 5: Left : volume rendering of a 3-D curvelet atom in the spatial domain corresponding to our implementation, cut by a vertical plane to see its inner structure. Right : the magnitude if its Fourier transform. The colorbar scale is valid only for the left image.

  PSNR of noisy data (dB) PSNR of denoised data (dB)LR-FCT Curvelab (curvelets) Curvelab (wavelets)

Figure 6 :

 6 Figure 6: Mean denoising PSNR versus noise level (as measured by the input PSNR) using different FCT implementations. The denoising PSNR was averaged over ten noise realizations and several datasets. The LR-FCT is in blue. Original FCT implementation of Curvelab using curvelets (red) and wavelets (orange) at the finest scale.

  data (dB) PSNR of noisy data (dB)

Figure 7 :

 7 Figure 7: Output PSNR as a function of the input PSNR for three video sequences. (a) mobile, (b) tempete, and (c) coastguard CIF sequence. The redundancy of each transform is indicated in parentheses on the legend.

  data (dB)PSNR of denoised data (dB)

Figure 8 :

 8 Figure 8: Top row: Mars Express observations at two different wavelengths. Bottom-left: two spectra at two distinct pixels. Bottom-right: output PSNR as a function of the input PSNR for different transforms; see Figure 7 for legend of the curves.

Figure 9 :

 9 Figure 9: Top : a sagittal ((y, z)) slice of the original synthetic MRI volume from BrainWeb[28]. Left column : the data with missing areas, random 80% missing voxels, and 10% missing z slices. Right : inpainting results with a LF-FCT+UDWT dictionary.

Figure 12 :

 12 Figure12: Close-ups of the centeral part of a frame from the shields video sequence de-interlaced using different methods. Top : Original, and de-interlaced result using STELA. Bottom : results using STIM and proposed method.

Figure 13 :

 13 Figure 13: Close-ups of a fast-moving part of the shields video sequence de-interlaced using different methods. Top : Original, and de-interlaced result using STELA. Bottom : results using STIM and proposed method.

Figure 14 :

 14 Figure 14: Close-ups of the stairs in the stockholm video sequence de-interlaced using different methods. Top : Original, and de-interlaced result using STELA. Bottom : results using STIM and proposed method.

Figure 15 :

 15 Figure15: Close-ups of a building in the stockholm video sequence de-interlaced using different methods. Top : Original, and de-interlaced result using STELA. Bottom : results using STIM and proposed method.

Figure 16 :

 16 Figure 16: Morphological component separation results from a data cube containing pieces of 3-D shells and 3-D gaussians. From left to right : original dataset, and morphological components extracted with the LR-FCT and the isotropic wavelet transform. Top : a 2-D slice of each 3-D cube. Bottom: maximum intensity projection view.

Figure 17 :

 17 Figure17: The same as Figure16, with an additive white Gaussian noise. The two components are jointly separated and denoised.

Table 1 :

 1 Redundancy of the original FCT and the proposed one in 2-D and 3-D, when wavelets (W) or curvelets (C) are used at the finest scale.

	Original FCT Proposed FCT
	C	W	C	W
	2-D 7.11	3.56	4.00	2.00
	3-D 24.38	5.42	10.29	2.29

Conclusion 16

On a Cartesian grid, the slopes are equispaced not the angles.

Here, we do not deal with an orthonormal version of the Meyer wavelet transform which can be viewed as a lapped orthogonal transform in the Fourier domain and efficiently implemented with folding.

FFT and periodization have another side effect in the spatial domain this time, that curvelet atoms localized near the boundaries wrap around to the opposite edge. Indeed,In the 2-D case, mirror-extended (reflexive boundary conditions) curvelets have been proposed in[START_REF] Demanet | Curvelets and wave atoms for mirror-extended images[END_REF] to mitigate the curvelet-looking artifacts due to periodization, but their redundancy is still as high as regular curvelets.

In fact, it should be ≈

as can be read from Table1.

The curvelet coefficients are normalized properly by the ℓ 2 norms of the associated curvelets so that this threshold is valid at all subbands. See also previous section.

A similar complicating coherence phenomenon occurs also in super-resolution.

Original central frame