A Lebret

E Petit

B Durning

J Hodel

A Rahmouni

P Decq

Volume Assessment of the Cerebrospinal Fluid Spaces for Computer Aided Diagnosis

Keywords: cerebrospinal fluid, hydrocephalus, 3D segmentation, topology preservation, computer-aided diagnosis

The present work aims to provide support in the diagnosis of hydrocephalus, which requires an assessment to the volumes of the cerebrospinal fluid (CSF) within its total, ventricular and subarachnoid spaces. In this paper we describe a fully automatic method to estimate the CSF volumes from a new 3D whole body MR imaging sequence. The method was developed using image properties as well as anatomical and geometrical features, completed with a topological assumption on the CSF shape. Experiments on phantoms and clinical data were performed and evaluated by comparing our assessments of volumes with those derived from a segmentation controlled by expert physicians. Then we show that a robust distinction between pathological cases and healthy adult people can be achieved by a linear discriminant analysis on volumes of the ventricular and intracranial subarachnoid spaces. We also find that healthy adults maintain a proportional relationship between these volumes.

Introduction

Hydrocephalus is caused by an abnormal accumulation of the cerebrospinal fluid (CSF) into cerebral ventricles and/or subarachnoid cisterns of the brain (see [START_REF] Harold | A contemporary definition and classification of hydrocephalus[END_REF] for more details). Obviously the measurement of mechanical features such as CSF pressures or volumes may be of great interest for assisted diagnosis and therapeutic [START_REF] Peña | Communicating hydrocephalus: The biomechanics of progressive ventricular enlargement revisited[END_REF][START_REF] Linninger | Cerebrospinal fluid flow in the normal and hydrocephalic human brain[END_REF].

As MRI is a reference technique for the CSF imaging, a new sequence that significantly highlights the CSF, has recently been developed for intracranial images [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF] and is also tested for whole body images. Due to the poor dynamic of the luminance perception of the human visual system, such images appear to be binary but in fact they are not (see Fig. 3,4), and furthermore they are disturbed by noise and artifacts.

For now, the volume assessment of the CSF spaces is achieved by expert physicians by extracting the total volume through thresholding and interactive removal of noise and unwanted structures. Next, this volume is separated into its ventricular, spinal and intracranial subarachnoid spaces. This leads to a functional parameter calculated by the ratio between the ventricular CSF space volume and the intracranial subarachnoid one. This ratio could be efficient to discriminate healthy adults from patients affected by an hydrocephalus. However, it should be pointed out that most of the previous steps used to be manually performed and that the overlapping of ratio ranges between healthy people and patients may result to misclassification.

In order to solve these problems, we first propose to fully automate the above segmentation and separation steps, and then to refine the classification by applying a linear discriminant 1 analysis to the ventricular and intracranial subarachnoid CSF volumes. We also show that healthy adults maintain a proportional relationship between these volumes.

To solve the segmentation and separation problems, we do not only use image properties, but also anatomical and geometrical features as well as a topological assumption about the CSF shape. The topological assumption becomes necessary because of both, the complex geometry of the CSF shape as shown in Fig. 1 (left), and the numerous image artifacts (especially motion artifacts), which can lead to changes in its topology. Benefits of associating topological constraints to medical image segmentation methods are shown in [START_REF] Ségonne | Segmentation of Medical Images under Topological Constraints[END_REF]. These constraints may prevent from any change in the topology of the anatomical structures to be extracted, such as unwanted created connections and/or holes. From a topological point of view, each anatomical structure including the CSF may be considered as a simple object such as a filled or a hollow sphere [START_REF] Bazin | Topology-preserving tissue classification of magnetic resonance brain images[END_REF][START_REF] Miri | Topology-preserving discrete deformable model: Application to multi-segmentation of brain MRI[END_REF][START_REF] Rueda | Topology-corrected segmentation and local intensity estimates for improved partial volume classification of brain cortex in MRI[END_REF]. Several works have successfully added such constraints to their approach in order to preserve the topology of anatomical objects throughout the segmentation process (deformable surface with an energy minimization goal [START_REF] Mangin | From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations[END_REF][START_REF] Poupon | Multi-object deformable templates dedicated to the segmentation of brain deep structures[END_REF], level set [START_REF] Han | Topology preserving level set method for geometric deformable models[END_REF], fast marching [START_REF] Bazin | Topology-preserving tissue classification of magnetic resonance brain images[END_REF], discrete deformable model [START_REF] Miri | Topology-preserving discrete deformable model: Application to multi-segmentation of brain MRI[END_REF], morphological operators [START_REF] Dokládal | Topologically controlled segmentation of 3D magnetic resonance images of the head by using morphological operators[END_REF][START_REF] Rueda | Topology-corrected segmentation and local intensity estimates for improved partial volume classification of brain cortex in MRI[END_REF], statistical an topological atlases registration approach [START_REF] Bazin | Homeomorphic brain image segmentation with topological and statistical atlases[END_REF]).

Here, from a topological point of view, the CSF space is considered as a filled sphere as shown in Fig. 1 (right). This property, used as a refinement constraint, brings both accuracy and robustness to our segmentation method.

The paper is organized as follows. Section 2 describes the segmentation methodology for the whole CSF. The separation of the CSF volume into its ventricular, intracranial and spinal subarachnoid sub-volumes is presented in section 3. In section 4, we report experiments on phantoms and clinical dataset. Concluding remarks are given in section 5.

Segmentation of the whole CSF

The total CSF volume is segmented by using image properties with the hypothesis that its shape is topologically equivalent to a filled sphere as shown in Fig. 1 (right). The segmentation is performed as follows:

1. Extract an inner object (C inner ) that is within the whole CSF as close as possible to its actual borders and has the same topology.

2. Extract an outer object of the CSF (C lim ) that includes the whole CSF but do not necessarily have the same topology.

3. Homotopically thicken C inner preserving its topology as long as it remains within C lim . The transformation is guided by a priority function based on a distance criterion.

Note that obtaining both, inner and outer objects, as close as possible to the object of interest, allows to speed up significantly the processing time.

Extraction of an inner object of the CSF

An inner object is extracted by applying the moment-preserving thresholding described in [START_REF] Tsai | Moment-preserving thresholding: a new approach[END_REF] on the initial image followed by a larger connected component extraction (Fig. 4-b).

Moment-preserving thresholding

This method is suited to "binary nature" images with a low contrast. When the histogram presents an excessive overlapping between classes, the moment-preserving method becomes more efficient than conventional clustering-based or entropy-based methods. In an axial cross-section of the spinal chord, the CSF has a disk topology (left).

The moment preservation assumes that statistical moments on original and thresholded images are identical. Let F = {f ijk } be a 3D image of L intensity levels. For each intensity level z, its associated normalized histogram is given by p(z) and the statistical moment M of order m for m = 1, 2, 3 is defined as:

M m = L-1 z=0 p(z) z m .
The resulting thresholded image G = {g ijk } may be considered to be composed of two gray levels z 0 and z 1 (z 0 < z 1 ). The optimal threshold T inner is obtained such that first three moments for F , the mean M 1 , the variance M 2 and the skewness M 3 , are equals to those for G when g ijk = z 0 (resp.

g ijk = z 1 ) if f ijk ≤ T inner (resp. f ijk > T inner ).
Let P 0 and P 1 be the percentages of voxels of G on both sides of T inner , the statistical moment B m of order m is: B m = 1 l=0 P l z m l . Moments preservation leads to solve the system:

1 l=0 P l z m l = L-1 z=0 p(z) z m , (1)
with the constraint P 0 + P 1 = 1 for m = 1, 2, 3. Equation (1) has four unknowns z 0 , z 1 , P 0 and P 1 , whose solution can be found in [START_REF] Tsai | Moment-preserving thresholding: a new approach[END_REF] and leads to a value for T inner that is given by:

P 0 = 1 L Tinner z=1 p(z). (2)
Figure 2 shows an histogram with the threshold T inner (= 38) obtained from the momentpreserving method, which allows to retrieve the inner CSF. T lim (= 13) is the threshold using the triangle thresholding, which is presented in Section 2.2 and allows to retrieve the outer object of the CSF. For comparison, the Otsu threshold T Otsu (= 53) is also shown.

In Tab. 1, the results obtained by the other thresholding methods reviewed in [START_REF] Sezgin | Survey over image thresholding techniques and quantitative performance evaluation[END_REF], for the histogram on Fig. 2, are also presented for comparison to that of the moment-preserving and result in more distant threshold values.

Retrieval of the inner object

Since G contains several thresholded objects, a connected component labeling is performed to extract the three largest connected components (the inner object C inner and both eyeballs E Table 1: Threshold values obtained on the histogram from Fig. 2 using various methods also reviewed in [START_REF] Sezgin | Survey over image thresholding techniques and quantitative performance evaluation[END_REF].

Method Threshold

Attribute_Huang [10] 1 Cluster_Kittler [START_REF] Kittler | Minimum error thresholding[END_REF] 5 Cluster_Ridler [START_REF] Ridler | Picture thresholding using an iterative selection method[END_REF] 52 Entropy_Kapur [START_REF] Kapur | A new method for gray-level picture thresholding using the entropy of the histogram[END_REF] 43 Entropy_Li [START_REF] Li | An iterative algorithm for minimum cross entropy thresholding[END_REF] 4 Entropy_Sahoo [START_REF] Sahoo | Threshold selection using renyi's entropy[END_REF] 43 Entropy_Shanbag [START_REF] Shanbag | Utilization of information measure as a means of image thresholding[END_REF] 241 Entropy_Yen [START_REF] Yen | A new criterion for automatic multilevel thresholding[END_REF] 53

E 2
). E 1 and E 2 are retained to determine a region of interest in Section 3.1. If C inner has holes, then T inner is automatically adjusted by increasing its value until the topological assumption is satisfied (E 1 and E 2 are not affected by this operation due to their intensity). By extracting the largest components (Fig. 4-b), residual noise as well as other small objects such as salivary ducts are also removed.

Thresholding of an outer object of the CSF

To extract an outer object of the CSF, another histogram-based thresholding method is carried out on the initial image as described in [START_REF] Zack | Automatic measurement of sister chromatid exchange frequency[END_REF]. It is followed by a largest connected component extraction to retrieve an outer object of the CSF (C lim ).

The triangle thresholding algorithm in [START_REF] Zack | Automatic measurement of sister chromatid exchange frequency[END_REF] is suitable for bimodal images for which the object to extract has a small amplitude and a large variance relatively to the background. Let p(z) be the histogram function for gray levels z. Maximum and minimum values, p max and p min , are determined from p(z) with their z values, z min and z max , and the points (z min , p min ) and (z max , p max ) are connected by a line segment δ. Let d e (a, b) be the Euclidean distance between two points a and b. For each point a = (z, p(z)), we can define the minimum distance from the line segment δ such that: d min (z, p(z)) = min b∈δ {d e (a, b)}. The optimal threshold T lim is obtained such that:

T lim = argmax z {d min (z, p(z))}. (3)
C lim is the object obtained as the largest connected component of the thresholded image.

Homotopic thickening of the inner object

An homotopic thickening is carried out on the inner object C inner . The thickening is guided by a priority function Ψ that is related to C lim and based on a distance criterion.

Priority function

The priority function Ψ = {ψ ijk } controls the order in which voxels are processed and eventually added to C inner . This function uses an Euclidean distance map denoted by D = {d ijk }, which is calculated by the linear algorithm proposed in [START_REF] Saito | New algorithms for euclidean distance transformation on an N-dimensional digitized picture with applications[END_REF]. The value ψ ijk of each voxel (i, j, k) in the outer object C lim is its distance to the nearest non-zero voxel of C inner such that:

ψ ijk = d ijk if (i, j, k) ∈ C lim +∞ otherwise. (4)

Guided homotopic thickening

Guided by the previous priority function Ψ, an homotopic thickening is applied to the inner object C inner by adding simple points.

A point is simple if its addition to or its removal from a binary object C does not change the topology of the object and of the background. Therefore, removing simple points does not change the number of connected components and holes (cavities and/or tunnels) of C and its complement C. It is shown in [START_REF] Bertrand | A new characterization of 3-dimensional simple points[END_REF] that a simple point x is locally characterized by two topological numbers denoted by T n (x, C) and T n(x, C) for (n, n) ∈ {(6, 26), [START_REF] Sahoo | Threshold selection using renyi's entropy[END_REF][START_REF] Duda | Pattern Classification[END_REF])}, such that T n (x, C) = 1 and T n(x, C) = 1. An efficient computation of T n (x, C) and T n(x, C) that only involves the 26-neighborhood is also described in [START_REF] Bertrand | A new characterization of 3-dimensional simple points[END_REF].

Let ∂C inner denotes the boundary of C inner such that

∂C inner = {x ∈ (C inner ∩ C lim ) | N n (x) ∩ C inner = ∅},
where C inner is the complement of C inner and N n (x) is the n-neighborhood of x with n ∈ {6, 26}. The thickening process consists to add to C inner simple points from ∂C inner as presented in Algorithm 1. It shows that boundary points (i.e. voxels) of ∂C inner are chosen iteratively, according to their priorities, and are incorporated into C inner if they are simple. The resulting thickened object is illustrated in Fig. 3-d and4-c, d.

Separation into CSF volumes

The separation of the CSF spaces is required to analyze the relationship between the ventricular and intracranial subarachnoid volumes. First, the spinal and intracranial regions are separated using geometrical and topological properties of the CSF at their interface. Second, to retrieve the ventricular space volume from that of the intracranial subarachnoid one, the cerebral aqueduct must be detected because we use it as a morphological reconstruction marker.

Algorithm 1: Guided homotopic thickening

Input: C inner , Ψ Output: C inner Make a priority queue Q from ∂C inner using ψ ijk in the ascending order

while Q = ∅ do Poll x = (i, j, k) from Q ; if x is simple for C inner then C inner ← C inner ∪ {x} ; Add new points of ∂C inner around x to Q ; end if end while (a) (b) (c) (d)

Separation of spinal and intracranial spaces

Expert physicians usually perform the disconnection of the spinal and intracranial spaces by tracking the foramen magnum (the hole in the bottom of the skull through which the spinal cord passes in order to be connected to the brain), which is even for an expert, difficult to detect on this type of images. The separation into spinal and intracranial spaces requires to locate the closest axial cross-section plane to that of the foramen magnum. Thus, the disconnection is achieved by browsing the 2D connected component on each axial plane from the spinal part, and by detecting any change in its topology (see the axial cross-section in Fig. 5). A region of interest (ROI) is first determined to reduce the processing time as well as to improve the robustness.

Determination of the ROI

An ROI R is used on one hand to locate the cutting axial cross-section, and on the other hand to retrieve the cerebral aqueduct in order to separate the intracranial space into its ventricular and subarachnoid spaces (Section 3.2). Apart from pathological cases, the eyeballs typically face the cerebral aqueduct in the axial plane. We consider the ROI forming a rectangular cuboid, each of whose edge is parallel to one of the i-, j-, k-axes of the image and such that:

• The ROI is located in the i-axis direction between the coronal plane that passes through the maximum posterior boundary of E 1 and E 2 , and the last coronal plane that still contains non-zero voxels.

• Let D be the larger diameter of E 1 and E 2 . The height of the ROI is set empirically to be 8D, such that the top face of the ROI is located with the distance 2D above the eyeballs and the bottom face is located with the distance 5D below the eyeballs. The height may be interactively resized to fit human variations.

• The depth is determined by the two sagittal planes between the medial boundaries of E 1 and E 2 .

Determination of the cut-off plane

2D connected components are counted on each axial plane from the bottom to the top of the ROI R. The first axial plane that shows a change in the topology of the connected components is the plane just above that of the searched cut-off (Fig. 5).

Intracranial volumes separation

The process to extract the ventricular space may fail because of some subarachnoid regions that are very close to the third and fourth ventricles. Figure 6 shows, for instance, the proximity between the interpeduncular cistern and the infundibular recess of the third ventricle.

To overcome this problem we first detect a specific tubular object. The cerebral aqueduct is a thin tubular structure (diameter: 1 to 3 mm; length: ∼14 mm) located in the median sagittal plane, which connects the third and fourth ventricles. It is also the longest tubular structure with the highest intensity range in the most median sagittal planes of the ROI R. Its enhancement is performed by a vessel segmentation method (see [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF]). Finally, ventricles are recovered by a morphological reconstruction [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF] using the cerebral aqueduct as a marker.

Detection of the cerebral aqueduct

Several methods to enhance vessel-like structures in 2D and 3D grayscale images are reviewed in [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF]. Some efficient methods are based on properties of eigenvalues of the Hessian matrix H [START_REF] Lorenz | Multi-scale line segmentation with automatic estimation of width, contrast and tangential direction in 2D and 3D medical images[END_REF][START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF][START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF][START_REF] Krissian | Model based detection of tubular structures in 3D images[END_REF][START_REF] Manniesing | Vessel enhancing diffusion: A scale space representation of vessel structures[END_REF]. These multi-scale methods operate in a Gaussian scale space on which they calculate second-order derivatives, build the Hessian matrix H, and decompose it depending on its eigenvalues λ 1 , λ 2 and λ 3 such that |λ 1 | < |λ 2 | < |λ 3 |. The eigenvalues are analyzed to determine the likelihood for each voxel (i, j, k) to belong to a curvilinear structure. This analysis is based on the following assumptions: 1) λ 1 ≈ 0; 2) λ 2 ≈ λ 3 < 0; and 3)

|λ 1 | ≪ |λ 2 |.
The tubular structure filter implemented in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] is used to enhance the cerebral aqueduct. Let F σ = {f σ ijk } the filter in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] such that:

f σ ijk = 0 if λ 2 > 0 or λ 3 > 0, (1 -e -R 2 A 2α 2 )e - R 2 B 2β 2 (1 -e -S 2 2c 2 ) otherwise, (5)
where

R A = |λ 2 | / |λ 3 |, R B = |λ 1 | / |λ 2 λ 3 |, S = λ 2 1 + λ 2 2 + λ 2 3
, with some weights α ∈ [0, 1], β ∈]0, 1] and c ∈ [0, 500]. Figure 7 compares extracted curvilinear structures in the ROI R using methods from [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF] (Fig. 7-b) and [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] (Fig. 7-c). It shows that unlike the filter in [START_REF] Sato | Three-dimensional multi-scale line filter for segmentation and visualization of curvilinear structures in medical images[END_REF], that in [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF] allows to dissociate strict tubular structures from others. The cerebral aqueduct part has more highlighted values in Fig. 7-c by adjusting the filter with a higher weight to detect lines (α = 0.7). The aqueduct is the longest connected component found in the most median sagittal planes.

(a) (b) (c)

Ventricles reconstruction

Considering the cerebral aqueduct part as the marker, we would like to reconstruct all the ventricles. For this problem, methods such as fast marching [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF], level set [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] or active contour models [START_REF] Laurent | On active contours models and balloons[END_REF] are well known. However, the following two problems may hinder the reconstruction by the above methods: 1) the cerebral aqueduct part is too small as a marker to initialize a model; 2) motion artifacts may generate significant intensity inhomogeneities even in the ventricular spaces.

Grayscale reconstruction [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF] is a very useful operator from mathematical morphology that provides satisfactory results here. Let δ

f,B d (g) be the elementary geodesic dilation of a grayscale image g inside f (g is called the marker image and f is the mask) such that δ

f,B d (g) = (g⊕B d )∧f , where ∧ stands for the point-wise minimum and g⊕B d is the geodesic dilation of g by an isotropic structuring element B d chosen as the unit ball. The geodesic dilation of size n ≥ 0 is obtained by: δ

(n) f,B d (g) = δ (1) f,B d • δ (1) f,B d • • • • • δ (1) f,B d (g) n times . (6)
The grayscale reconstruction by dilation ρ f,B d (g) of f from g is calculated by iterating geodesic dilations of g inside f until idempotence such that:

ρ f,B d (g) = ∨ n≥1 δ (n) f,B d (g) , (7)
where ∨ stands for the point-wise maximum.

The ventricular space volume is retrieved by dilation of the previously detected cerebral aqueduct used as the marker image, using the parallel reconstruction algorithm proposed in [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF]. Figure 8 shows the result of such a reconstruction using a 6-connectivity.

CSF volumes analysis

As the ground truth is not available for clinical data, the segmentation method for the total volume assessment was first validated on phantom MR images. Experiments were carried out on four phantoms of known volume (image size: 250 × 250 × 128 voxels) and 17 clinical MR images (image size: ∼ 260 × 790 × 160 voxels).

Images acquisition

MR images were acquired in the sagittal plane on a 1.5 T system (Magnetom Avanto; Siemens Medical Solutions, Erlangen, Germany). The MR sequence in [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF], called SPACE (Sampling Perfection with Application optimized Contrast using different flip-angle Evolution), is a variant of the T2-weighted turbo spin echo sequence with variable flip-angles. The sequence was as follows: repetition time TR (ms)/echo time TE (ms), 2400/762; turbo factor of 141; 250 × 250 mm field of view; 256×256 acquisition matrix; 1 mm isotropic resolution; number of excitations 1.4; 160 slices; acquisition time, 3 min [9].

Complexity and computation time

The complexity in time is dominated by the homotopic thickening step and is O(n log n) where n is the number of points in C inner ∩ C lim . It may be reduced to be linear by parallelizing the homotopic thickening process. The total computation time for segmentation and separation was performed from 3 to 6 minutes on an Intel Pentium Dual Core 1.60 GHz / 2 GB for each clinical MR image.

Validation on phantoms

Images were acquired from blood bags and synthetic resin phantoms with various shapes and volumes (see an example in Fig. 9). Phantoms 1 and 2 have regular shape while phantoms 3 and 4 have more tortuous one.

It has been observed (Fig. 10) that both segmentations, the semi-manual by an expert as explained in [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF] and the grayscale reconstruction [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF], have resulted in an overestimation. Instead, the methods of homotopic thickening and level set initialized by a fast marching have produced underestimated but close results. The guided homotopic thickening on phantoms has resulted in an accuracy of 98.5% ± 1.8 against 94.6% ± 4.0 for the semi-manual segmentation by an expert.

Statistical analysis on clinical data

Clinical images were acquired from different subjects (ages between 25 and 84 years): seven healthy volunteers, two patients with a communicating hydrocephalus just two days after their surgery, and eight patients among which seven communicating hydrocephalus (CH) and one non-communicating hydrocephalus (NCH). Due to the limited number of clinical images, the CH and NCH patients have not been distinguished. In addition, the two patients who just have undergone their surgery were removed from the analysis because according to physicians, the CSF has not had enough time to reach the whole subarachnoidian space.

Classification of clinical data

Following the proposal in [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF] to use ratios to classify data, a linear discriminant analysis [START_REF] Duda | Pattern Classification[END_REF] was first performed on the clinical dataset using volumes of the ventricular and intracranial subarachnoid spaces as the input variables. Let TP, FP, TN and FN the respective numbers of true positives, false positives, true negatives and false negatives. The sensitivity S e and the specificity (or recall) S p are respectively defined as S e = TP/(TP+FN) and S p = TN/(TN+FP). The additional F 1 -score metric is also used to validate the classifier. It is defined as F 1 = 2pS e /(p + S e ) where p stands for the precision, defined by p = TP/(TP + FP). The predicted classification on the clinical data results in a sensitivity/specificity/F 1 -score of 100%, which is a suitable result to discriminate between healthy people and pathological cases (Fig. 11). A similar observation was made on data from [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF].

This first result appears to confirm such as in [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF] that the relationship between the ventricular CSF space volume and the intracranial subarachnoid could be an efficient functional parameter to distinguish healthy people from pathological cases.

Linear regression model for healthy adults

Observing the previous result has provided us the idea of carrying out a linear regression on the healthy people subset that results from the linear discriminant analysis step. Let V V and V S the respective volumes of the ventricular and intracranial subarachnoid spaces. We use the linear model: V S = β 0 + β 1 V V + ǫ, and optimize the intercept β 0 and the slope β 1 by minimizing the residual error ǫ. We obtain β 0 = 0 and β 1 = 0.1071, namely:

V V = 0.1071V S (8)
Table 2 includes some statistics about the model: 1) the coefficient of determination R 2 (= 98.43%) indicates that the predictor explains rather well the answer; 2) the F statistic (= 312.8) and the p-value (= 1.06e-05 ≪ 0.05) show further that the null hypothesis (β 1 = 0) can be rejected with a strong presumption. In addition, the graphical residuals analysis reinforces the linear regression assumptions (Fig. 12). The left figure of Fig. 13 shows the distribution of ventricular and intracranial subarachnoid volumes with the model obtained from [START_REF] Han | Topology preserving level set method for geometric deformable models[END_REF]. The right figure illustrates the same model with the distribution of the data from [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF].

Conclusion

This work uses a new sequence of whole body MRI and succeeds in providing good estimations of the volumes of the CSF. The volumes are automatically retrieved through segmentation and separation steps, which use image properties, anatomical and geometrical features, as well as topological assumptions on the CSF. The proposed method seems to confirm the assumption that the proportional relationship between the ventricular and intracranial subarachnoid volumes allows to discriminate healthy and pathological cases. As the database is currently in expansion, our tool facilitates the refinement of these results in future work.

This new ability to quantify CSF volumes, which is a main parameter for the diagnosis, for the treatment decisions, for the monitoring of patients, and finally for clinical research, may represent an important breakthrough in the field of computer-aided neuro-imaging. Guided Homotopic Tickening Geodesic Reconstruction [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF] Fast Marching + Level Set [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] Figure 10: Comparisons of the total volume assessment of phantoms for different methods. Top: total volume measurement for objects that result from: the semi-manual segmentation by an expert as described in [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF]; the initial inner object extraction (Section 2.1); the homotopic thickening process (Section 2.3); a grayscale reconstruction method [START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF] and a fast marching followed by a level set method [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF], both using the inner object as a marker. Bottom: errors on the total volume measurement for the previous methods.

Figure 1 :

 1 Figure1: The CSF around the brain is considered as topologically equivalent to a filled sphere (right). In an axial cross-section of the spinal chord, the CSF has a disk topology (left).

Figure 2 :

 2 Figure 2: Image histogram sample from the clinical dataset and its zooming (vertical log scale) at the upper-right corner.

Figure 3 :Figure 4 :RFigure 5 :

 345 Figure 3: Segmentation process (axial cross-section), which is easily observed in the red square region. (a) original image; (b) contrast-enhanced original image; (c) inner object image; (d) homotopic thickened image.

1 infundibular recess 2 Figure 6 :

 126 Figure 6: Proximity of the interpeduncular cistern (1) to the infundibular recess of the third ventricle (2). The small rectangular part in the left figure is zoomed in the right figure.

Figure 7 :

 7 Figure 7: Detection of the cerebral aqueduct in the ROI R (median plane cross-section). From left to right: (a) original thickened image; (b) filtered image by [28] (σ = 0.5); (c) filtered image by [7] using 5 scale spaces σ from 0.4 to 0.9 and with α = 0.7, β = 0.2 and c = 200.

Figure 8 :

 8 Figure 8: Ventricular space retrieving using a grayscale reconstruction by dilation of the cerebral aqueduct inside the original image: (a) the original image (the mask) and the cerebral aqueduct (the marker) colored in red; (b) result of the grayscale reconstruction.

Figure 9 :

 9 Figure 9: Example of a synthetic resin phantom MR image. Cross-section of the phantom (left). The 3D reconstruction of the rectangular region is shown in the right figure.

  Semi-manual[START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF] Moments preserving + Connected Component Extr.Guided Homotopic Tickening Geodesic Reconstruction[START_REF] Vincent | Morphological grayscale reconstruction in image analysis: Applications and efficient algorithms[END_REF] Fast Marching + Level Set[START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] [START_REF] Hodel | Intracranial cerebrospinal fluid spaces imaging using a pulse-triggered three-dimensional turbo spin echo MR sequence with variable flip-angle distribution[END_REF] Moments preserving + Connected Component Extr.

Table 2 :

 2 Optimal linear regression model from the clinical data set.

	Model	R 2	σ 2	Standard Error	F	p-value
	VV = 0.1071VS	0.9843	12.25	0.0061	312.8	1.06e-05


Tieyong Zeng
email: zeng@hkbu.edu.hk

Malgouyres

Matching Pursuit Shrinkage in Hilbert Spaces

Keywords: Dictionary, Matching Pursuit, shrinkage, sparse representation. EDICS Category: SMR-REP, TEC-RST

In this paper, we study a variant of the Matching Pursuit named Matching Pursuit Shrinkage. Similarly to the Matching Pursuit it seeks for an approximation of a datum living in a Hilbert space by a sparse linear expansion in an enumerable set of atoms. The difference with the usual Matching Pursuit is that, once an atom has been selected, we do not erase all the information along the direction of this atom.

Doing so, we can evolve slowly along that direction. The goal is to attenuate the negative impact of bad atom selections.

We analyse the link between the shrinkage function used by the algorithm and the fact that the result belongs to an l p space.

I. INTRODUCTION

A. Recollection on sparse approximation

Finding a sparse approximation of a data in a Hilbert space is a reccurent problem in applied science.

The problem is to approximate a datum v ∈ H (H is a Hilbert space of finite or infinite dimension) by a linear expansion in a dictionary of known atoms (ψ i ) i∈I :

v ∼ i∈I λ i ψ i ,
where (λ i ) i∈I ∈ R I . The approximation is needed because v is usually corrupted by noise. Also, it is sometimes preferable to search for an approximation which is coarser than the noise requires. Doing so we favors desired/expected properties of the coordinates (λ i ) i∈I .

Moreover, the dictionary is usually overcomplete. This offers the freedom to select among all the possible sets of coordinates one of those agreeing with some prior knowledge or desired property of the coordinates. The property receiving most of the attention is sparsity. Heuristically, we select the set of coordinates offering the "simplest" explanation of the datum. Rigorously, for a given accuracy after reconstruction, we want l 0 ((λ i ) i∈I ) def = #{i ∈ I, λ i = 0}, to be as small as possible, where # denotes the cardinality of a set.

Unfortunately, problems similar to    minimize l 0 ((λ i ) i∈I )

under the constraint i∈I λ i ψ i -v ≤ τ (1)
where τ > 0 and . is the norm associated with the scalar product of the considered Hilbert space, are known to be NP-Hard in general (see [START_REF] Davis | Adaptive greedy approximations[END_REF]).

As a conclusion, solving [START_REF] Davis | Adaptive greedy approximations[END_REF] is both an open and interesting problem. It receives a lot of attention and it is impossible to list all the contributions to its resolution. Before describing the most popular technics, we give in the next section the algorithm studied in this paper. It will then be simpler to motivate our proposal.

B. The Matching Pursuit Shrinkage

The Matching Pursuit Shrinkage (MPS) is very similar to the usual Matching Pursuit (MP) algorithm (see [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]). The main difference is that it uses a shrinkage 1 function θ : R → R . We describe the algorithm in Table I.

• Input : A datum v, a dictionary (ψi)i∈I , a shrinkage function θ and α ∈ [0, 1]

• Output : Coordinates (sn, γn) n∈N

• The algorithm

-Initialize R 0 v = v -Repeat until convergence (loop in n)
1) Select a well correlated atom ψγ n such that

| ψγ n , R n v | ≥ α sup i∈I | R n v, ψi |; (2)
2) Evolve along ψγ n

R n v = snψγ n + R n+1 v, (3)
where Several convergence criterion might be considered but, for simplicity, we always assume that the algorithm stops whenever s n = 0.

sn = θ(Mn) with Mn = R n v, ψγ n . (4)
Whenever they exist, we can construct coordinates

λ i = n∈N:γn=i s n , ∀i ∈ I (5)
from the result of the MPS. We also consider (when they exist)

u = i∈I λ i ψ i = +∞ n=0 s n ψ γn .
Notice that if we sum (3) for n = 0 . . . N -1, we obtain

v = N -1 n=0 s n ψ γn + R N v. (6)
This explains the name "residual error" for R N v.

C. Other algorithms promoting sparsity

One of the oldest and simplest algorithm for building a sparse approximation is the Matching Pursuit (MP) [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF] or Projection Pursuit [START_REF] Friedman | Projection pursuit regression[END_REF]. It corresponds to the algorithm of Table I when θ is the identity (i.e.

s n = M n ).
July 8, 2009 DRAFT

In finite dimension (see [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]) and in infinite dimension but under restrictive conditions on the dictionary and the signal (see [START_REF] Gribonval | On the exponential convergence of matching pursuits in quasi-incoherent dictionaries[END_REF]), the MP is known to converge exponentially. When no hypotheses are made on the dictionary, we only know that the MP converges (see [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]). Some examples show that we cannot expect a "good" converge rate in the most general setting (see [START_REF] Devore | Some remarks on greedy algorithms[END_REF]). Though the MP and the best k-term approximation have a similar convergence, when the dictionary is "quasi-orthogonal" (see [START_REF] Temlyakov | Greedy algorithms and m-term approximation with regard to redundant dictionaries[END_REF]).

There exists "fast" variants of the MP (see [START_REF] Gribonval | Harmonic decomposition of audio signals with matching pursuit[END_REF]). Also, a real-time implementation of the MP is available for audio signal processing (see [START_REF] Krstulovic | MPTK: Matching Pursuit made tractable[END_REF]). The improved performance are obtained by carefully optimizing the structures, algorithms and their implementation. In particular, the update of ( R n v, ψ i ) i∈I and the computation of γ n satisfying (2) (in Table I) are implemented in a very efficient way. Each iteration of the MP is typically of complexity O(log(#I)). These optimization are possible because one coordinate only is updated. If K coordinates are modified at each iteration, we obtain a complexity O(K +log(#I)).

This might be less favorable when K is large. Althought its approximation performances are not as good as most modern models/algorithms, these acceleration make the MP a usefull algorithm.

The accelerations decribed in [START_REF] Krstulovic | MPTK: Matching Pursuit made tractable[END_REF] can be applied to the MPS, as described in Table I. The potential advantage of introducing a shrinkage function θ is to attenuate the mistakes in the selection of a coordinate γ n . Let us underline that avoiding wrong selection of coordinates is one of the key ingredient of modern variants of the MP such as CoSaMP [START_REF] Needell | Cosamp: Iterative recovery from incomplete and inaccurate samples[END_REF], Subspace Pursuit [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] and Iterative Hard Thresholding [START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF].

However, especially when the solution we are looking for is moderately sparse, those algorithms are more computationaly intensive.

Let us go back in time. The most famous variant of the MP is the Orthogonal Matching Pursuit (OMP)

(see [START_REF] Pati | Orthogonal matching pursuit : Recursive function approximation with applications to wavelet decomposition[END_REF]). In Table I, it replaces the update rule (3) by an orthogonal projection onto the subspace generated by the selected atoms. It is known to provide sparser solutions than the regular MP. From the computational point of view, it has two drawbacks. Firstly, although several attemps have been made to optimize it (see [START_REF] Mallat | Adaptive time frequency decompositions[END_REF], [START_REF] Cotter | Forward sequential algorithms for best basis selection[END_REF]), the orthogonal projection is computationaly expensive and often requires too much memory. Secondly, every selected coordinate is modified. As a consequence, the adaptation of the optimization performed in [START_REF] Krstulovic | MPTK: Matching Pursuit made tractable[END_REF] would only be efficient when the result is very sparse. Algorithms such as the Gradient Pursuit (see [START_REF] Blumensath | Gradient pursuits[END_REF]) approximately solve the OMP at a cost more similar to the cost of the MP. However, at each iteration, they typically update all the selected coordinates. The computational cost of the Gradient Pursuit is therefore more important than the cost of a fast implementation of the MP, when the solution is moderetely sparse.

Finally, the l 1 regularization (also named Basis Pursuit and Basis Pursuit Denoising, see [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]

λ i ψ i 2 + β i∈I |λ i |
and it is very efficient for providing sparse approximations of v ∈ H. However, its resolution remains (and will probably remain in a near futur) a challenge for large scale problems. A famous (and representative) solver of the l 1 regularization problem is the Iterative Soft Thresholding (see [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]). It updates all the coordinates at each iteration and often requires many iterations before it reaches a suitable convergence level. It is interesting to notice that, in this context, the impact of the choice of the shrinkage function is well understood (see [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF]): Every proximal threholding function corresponds to a different objective function.

Inspired by the l 1 regularization problem, a "coordinatewise optimization algorithms" has been proposed in [START_REF] Friedman | Pathwise coordinate optimization[END_REF]. It performs a soft thresholding, sequencially on each coordinate. The "greedy coordinate descent" proposed in [START_REF] Wu | Coordinate descent algorithms for lasso penalized regression[END_REF] is similar but selects the coordinates according to a criteria similar to the MP. Because they only update one coordinate at each iteration, these algorithms can benefit from the optimization proposed in [START_REF] Krstulovic | MPTK: Matching Pursuit made tractable[END_REF].

D. Notations

The following notations and hypotheses hold all along the paper.

The datum v belongs to a Hilbert space H. The space H might be of finite or infinite dimension. For any two elements u and w in H, their scalar product is denoted by u, w . As usual, the norm of u ∈ H is defined by

u def = u, u . The dictionnary (ψ i ) i∈I is made of atoms ψ i ∈ H, such that ψ i = 1, for all i ∈ I.
We sometimes denote the dictionary by D. For simplicity, we assume that I is enumerable. In particular, the supremum in (2) may not be reached. In such a case, the MPS is only defined for α < 1.

For any u ∈ H, we denote u

D def = sup i∈I | u, ψ i |. We denote def = Span{D} (7)
the closed linear span of the elements of D. We denote V ⊥ the orthogonal complement of V in H. We denote the orthogonal projection onto V and V ⊥ by P V and P V ⊥ .

The sequences (s n ) n∈N , (γ n ) n∈N , (R n v) n∈N are always defined according to Table I. The coordinates (λ i ) i∈I are according to [START_REF] Devore | Some remarks on greedy algorithms[END_REF].

We also use the standard notations : sgn(t) = 1, if t ≥ 0 and -1, if t < 0; # denotes the cardinal of a set; ⌊.⌋ is the floor function.

July 8, 2009 DRAFT

E. Overview

In Section II, we define shrinkage, thresholding and gap functions. We also illustrate these definitions by several examples. In section III, we prove that as soon as θ is a shrinkage function: (R n v) n∈N converges and n∈N s n ψ γn exists. We also prove that (s n ) n∈N is square summable. In Section IV, we prove that when θ is a thresholding function, (s n ) n∈N is absolutely summable. This implies in particular that (λ i ) i∈I exists and is absolutely summable. In Section V, we prove that when θ is a gap function, the sequence (s n ) n∈N is finite. Again, this implies that (λ i ) i∈I exists and is finite. Finally, in Section VI, we evaluate n∈N s n ψ γn -P V v D , when θ is a shrinkage function.

II. GENERAL SHRINKAGE FUNCTIONS

A. Definitions Definition 1: A function θ(•) : R → R is called a shrinkage function if and only if it satisfies:

1) θ(•) is nondecreasing, i.e, ∀t, t ′ ∈ R, t ≤ t ′ =⇒ θ(t) ≤ θ(t ′ ); 2) θ(•) shrinks the amplitude, i.e, ∀t ∈ R, |θ(t)| ≤ |t|.
Notice that this implies

θ(0) = 0, and
θ(-t) ≤ 0 ≤ θ(t), ∀t ≥ 0. (8)
Therefore, for any shrinkage function θ(•) and any t ∈ R, we know that:

if t ≥ 0, 0 ≤ θ(t) ≤ t and 0 ≤ θ(t)(t -θ(t)), if t ≤ 0, 0 ≥ θ(t) ≥ t and 0 ≤ θ(t)(t -θ(t)).
As a conclusion,

∀t ∈ R, θ(t)(t -θ(t)) ≥ 0. (9)
The inequality (8) also garantees that

∀t ∈ R, |t| |θ(t)| = tθ(t). (10
τ -> 0, i.e. ∃τ > 0, ∀x ∈ R, |x| ≤ τ ⇒ θ(x) = 0. (11)
If θ(•) is a thresholding function, we trivially have

0 < τ -≤ τ + .
The internal and external thresholds are illustrated on Figure 1.

t 0 y y = θ(t) y = t -τ + τ -
gap(θ) def = inf t:θ(t) =0 θ 2 (t) + 2θ(t)(t -θ(t)). (12)
If gap(θ) > 0, we call θ a gap shrinkage function and, if gap(θ) = 0, the function is called a non-gap shrinkage function.

July 8, 2009 DRAFT

The following relation exists between the gap and the internal threshold of a shrinkage function. It proves in particular that any gap shrinkage function is a thresholding function.

Proposition 1: For any gap function θ(•), we have

gap(θ) ≤ τ -
where τ -is the internal threshold of θ(•).

Proof: The proof is given in Appendix.

B. Examples

Let us illustrate the above definitions through some examples.

1) For τ > 0, the soft thresholding function ρ τ (•) defined by

ρ τ (t) = sgn(t) • max(|t| -τ, 0).
is a thresholding function and it is a non-gap shrinkage function, i.e., gap(ρ τ ) = 0.

2) For τ > 0, the hard thresholding function defined by

h τ (t) =    t , if |t| > τ, 0 , otherwise.
is a thresholding function and it is a gap shrinkage function with gap τ .

3) The identity function defined as:

i(t) = t, ∀t ∈ R, (13)
is not a thresholding function and it is a non-gap shrinkage function.

4) For τ > 0, the Non-Negative Garrote threshold function (see [START_REF] Gao | Wavelet shrinkage denoising using the non-negative garrote[END_REF]) defined as:

δ G τ (t) = t max 0, 1 - τ 2 t 2 , ∀t ∈ R, (14)
is a thresholding function and it is non-gap.

5) For 0 < τ 1 < τ 2 , the firm shrinkage function (see [START_REF] Gao | Waveshrink with firm shrinkage[END_REF]) defined as:

δ τ1,τ2 (t) =          0, if |t| ≤ τ 1 ; sgn(t) τ2(|t|-τ1) τ2-τ1 if τ 1 < |t| < τ 2 ; t, if |t| ≥ τ 2 , (15)
is a thresholding function and it is non-gap.

July 8, 2009 DRAFT 6) For p ∈ N, τ > 0, the generalized threshold function (see [START_REF] Zhao | Wavelet shrinkage denoising by generalized threshold function[END_REF]) defined as:

δ p τ (t) =    t, if |t| ≤ τ ; t -τ p t p-1 (sgn(t) p ), if |t| > τ, (16)
a thresholding function and it is non-gap.

III. CONVERGENCE OF THE MP SHRINKAGE FOR A SHRINKAGE FUNCTION

This section is devoted to prove that under mild condition, the MP shrinkage algorithm converges.

Proposition 2: Let (ψ i ) i∈I be a normed dictionary, v ∈ H and θ(•) be a shrinkage function. For any M > 0 and any v ∈ H, the quantities defined in Table I satisfy:

v 2 = M -1 n=0 s 2 n + 2s n (M n -s n ) + R M v 2 . ( 17
)
As a consequence, we have

v 2 ≥ M -1 n=0 s 2 n + R M v 2 , ( 18
) +∞ n=0 s 2 n < +∞, (19)
+∞ n=0 |s n | |M n | < +∞, ( 20
) ( R n ) n∈N is nonincreasing. ( 21
)
Proof: We can deduce from

R n+1 v = R n v -s n ψ γn ,
and ψ γn , ψ γn = 1 that

R n+1 v 2 = R n v 2 -2s n R n v, ψ γn + s 2 n = R n v 2 -2s n (M n -s n ) -s 2 n .
Summing these equalities for all n = 0, . . . , M -1, we obtain after simplification

R M v 2 = R 0 v 2 - M -1 n=0 (s 2 n + 2s n (M n -s n )).
We then obtain (17

) from R 0 v = v.
Using (9), we know that

s n (M n -s n ) = θ(M n )(M n -θ(M n )) ≥ 0.
July 8, 2009 DRAFT

Together with [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] this leads to [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF].

Notice that this also provides [START_REF] Gao | Wavelet shrinkage denoising using the non-negative garrote[END_REF]. Moreover, [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF] garantees that ( M n=0 s 2 n ) M ∈N is a bounded increasing sequence. It converges and ( 19) holds. We also have

2 M -1 n=0 |s n | |M n | = 2 M -1 n=0 s n M n from (10) = v 2 -R M v 2 + M -1 n=0 s 2 n from (17) ≤ v 2 + +∞ n=0 s 2 n .
This ensures that (20) holds.

Now we can prove the convergence of the MP algorithm.

Theorem 1: Let (ψ i ) i∈I be a normed dictionary, v ∈ H and θ(•) be a shrinkage function. The sequences defined in (4) satisfy:

(R n v) n∈N converges.
As a consequence, +∞ n=0

s n ψ γn exists.

We denote the limit of (R n v) n∈N by R +∞ v and we trivially have

v = +∞ n=0 s n ψ γn + R +∞ v.
Proof: The proof is based on Jones' proof for the convergence of projection pursuit regressions (see [START_REF] Huber | Projection pursuit[END_REF]) and the proof of Theorem 1 in [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF].

First notice that the statement of the proposition is trivial for v = 0. We further assume that v = 0.

In order to prove the theorem, we prove that the sequence (R n v) n∈N is a Cauchy sequence. Before doing so, let us start with some preliminaries.

Notice first that for all w 1 , w 2 ∈ H, we have:

w 1 -w 2 2 = w 1 2 -w 2 2 -2 w 2 , w 1 -w 2 ≤ ≤ w 1 2 -w 2 2 + 2| w 2 , w 1 -w 2 |. (22)
Moreover, for N 2 > N 1 ≥ 0, from (6) we have

R N1 v -R N2 v = N2-1 n=N1 s n ψ γn . ( 23
)
July 8, 2009 DRAFT Finally, for any n ≥ 0 and any m ≥ 0,

| R m v, s n ψ γn | = |s n | | ψ γn , R m v | ≤ |s n | sup i∈I | ψ i , R m v | ≤ 1 α |s n | |M m |. ( 24
)
Let us now consider N 2 > N 1 ≥ 0. Using ( 22), ( 23) and ( 24), we obtain

R N1 v -R N2 v 2 ≤ R N1 v 2 -R N2 v 2 + 2| R N2 v, N2-1 n=N1 s n ψ γn | ≤ R N1 v 2 -R N2 v 2 + 2 α |M N2 | N2-1 n=N1 |s n |. ( 25
)
Using ( 21) of Proposition 2, we know that the sequence ( R n v ) n∈N is non-negative and nonincreasing. Therefore, it converges to some value R ∞ and for any ǫ > 0, there exits

K > 0 such that for all m > K, R 2 ∞ ≤ R m v 2 ≤ R 2 ∞ + ǫ 2 .
As a consequence, for any

N 2 > N 1 ≥ K, R N1 v -R N2 v 2 ≤ ǫ 2 + 2 α |M N2 | N2 n=N1 |s n | (26)
Using [START_REF] Wu | Coordinate descent algorithms for lasso penalized regression[END_REF] • The first one is that: +∞ n=0 |s n | < +∞. In this case, we know that there is K ′ > 0 such that for any

N 2 > N 1 ≥ K ′ N2 n=N1 |s n | ≤ α 2 v ǫ 2 .
Moreover, from [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] we know that

|M N2 | = | R N2 v, ψ γN 2 | ≤ R N2 v ≤ v .
So [START_REF] Candes | The dantzig selector: Statistical estimation when p is much larger than n[END_REF] becomes : for any ǫ > 0 there are K and K ′ > 0 such that for any

N 2 > N 1 ≥ max(K, K ′ ) R N1 v -R N2 v 2 ≤ ǫ 2 + ǫ 2 .
As a conclusion (R n v) n∈N is a Cauchy sequence.

July 8, 2009 DRAFT

• The second one is that: lim inf q→+∞ |M q | q n=0 |s n | = 0. In this case, let ǫ > 0 and let p > 0 be an integer. We are going to estimate R m v -R m+p v , for m > K (K is such that (26) holds).

First, there is q > m + p such that

|M q | q n=0 |s n | ≤ α 2 ǫ 2 . ( 27
)
Moreover, we can decompose

R m v -R m+p v ≤ R m v -R q v + R m+p v -R q v .
Applying [START_REF] Candes | The dantzig selector: Statistical estimation when p is much larger than n[END_REF] with N 1 = m and N 2 = q and using ( 27) we obtain

R m v -R q v 2 ≤ ǫ 2 + ǫ 2 .
Similarly, applying [START_REF] Candes | The dantzig selector: Statistical estimation when p is much larger than n[END_REF] for N 1 = m + p and N 2 = q and using ( 27) we obtain

R m+p v -R q v 2 ≤ ǫ 2 + ǫ 2 .
Hence, we finally obtain

R m v -R m+p v ≤ 2 √ 2ǫ, which proves that (R n v) n∈N is a Cauchy sequence.
As a conclusion, (R n v) n∈N converges. The second statement directly follows from (6).

Proposition 2 ensures that

+∞ n=0 |s n | 2 (28) exists.
IV. l 1 NORM BOUNDS In general, when H is an infinite dimensional space, we have no guarantee that

+∞ n=0 |s n | (29)
exists. A simple counter example consists in considering (ψ i ) i∈I a Riesz basis (for definition, see [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF])

of H, v = i∈I s i ψ i ∈ H such that i∈I |s i | diverges and θ(t) ≡ t.
Below, we prove that (29) exists, whatever v ∈ H and whatever the dictionary, as soon as θ is a thresholding function.

July 8, 2009 DRAFT Proposition 3: Let (ψ i ) i∈I be a normed dictionary, v ∈ H and θ(•) be a thresholding function. The quantities defined in Table I satisfy:

+∞ n=0 |s n | ≤ v 2 -R +∞ v 2 τ - ≤ v 2 τ -, (30)
where τ -> 0 denotes the internal threshold as defined in the Definition 2.

Proof: Let M ∈ N fixed. Using [START_REF] Combettes | Proximal thresholding algorithm for minimization over orthonormal bases[END_REF], we know that

M -1 n=0 s 2 n ≤ v 2 -R M v 2 .
Together with [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], this leads to

M -1 n=0 s n M n = 1 2 v 2 + M -1 n=0 s 2 n -R M v 2 ≤ v 2 -R M v 2 .
Using [START_REF] Dai | Subspace pursuit for compressive sensing signal reconstruction[END_REF] and the fact that θ(•) is a thresholding function, for any n ∈ N, we have:

s n M n = |s n ||M n | ≥ τ -|s n |,
where the last inequality is obtained via the discussing on two cases: s n = 0 or s n = 0.

As a conclusion for all M ∈ N we have

M -1 n=0 |s n | ≤ v 2 -R M v 2 τ - . (31)
Letting M go to infinity, we obtain (30).

Remark 1:

The above upper bound does not depend on the dictionary (ψ i ) i∈I . It holds for any v ∈ H.

We therefore do not expect this bound to be tight in any dedicated or applicative context.

Remark 2: As a side effect, the above proposition garantees that the coordinates λ i exist for all i ∈ I (see [START_REF] Devore | Some remarks on greedy algorithms[END_REF]). We even know that

i∈I |λ i | < +∞. V. l 0 BOUNDS If θ(•)
is a gap shrinkage function the MP shrinkage stops automatically after a finite number of iterations.

Proposition 4: Let (ψ i ) i∈I be a normed dictionary, v ∈ H and θ(•) be a gap shrinkage function (i.e.

gap(θ) > 0). The sequence (s n ) n∈N defined in Table I satisfies:

#{n|s n = 0} ≤ ⌊ v 2 gap(θ) 2 ⌋.
Proof: Suppose that the sequence (s n ) n∈N contains M non-zero terms. Observing Definition 3, for each s n = 0, we have:

s 2 n + 2s n (M n -s n ) ≥ gap(θ) 2 ,
where we recall that M n = R n v, ψ γn , s n = θ(M n ).

From [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF], we know that:

v 2 ≥ n∈N:sn =0 s 2 n + 2s n (M n -s n ) ≥ M • gap(θ) 2 .
Noting that M is integer, we have:

M ≤ ⌊ v 2 gap(θ) 2 ⌋.
Remark 3: An interesting consequence of the proposition is that

#{i ∈ I, λ i = 0} ≤ ⌊ v 2 gap(θ) 2 ⌋.
In words, v is approximated with less than ⌊ v 2 gap(θ) 2 ⌋ non-zero coordinates.

VI. BOUND ON THE RESIDUAL ERROR

In this section, we are interested in the residual error norm. The result concerns shrinkage functions.

Before stating the result, let us give the following lemma:

Lemma 1: Let (ψ i ) i∈I be a normed dictionary, v ∈ H and θ(•) be a shrinkage function. The sequence (M n ) n∈N defined in Eq.( 4) satisfies: The latter statement is impossible since, from [START_REF] Friedman | Pathwise coordinate optimization[END_REF], we know that lim n→+∞ s n = 0. This proves (32).

lim sup n→+∞ M n ≤ sup t:θ(t)=0 t, (32)
The proof of (33) is similar.

In particular, if the external threshold of θ(•) is zero (i.e. τ + = 0),

lim n→+∞ M n = 0, since sup t:θ(t)=0 t = inf t:θ(t)=0 t = 0.
Recall that we have defined the semi-norm on H as

u D def = sup i∈I | u, ψ i |, ∀u ∈ H.
Notice that • D is a norm as soon as D generates H. Geometrically,

{u ∈ H, u D ≤ τ } is a polyhedron, for τ ≥ 0.
Recall that in [START_REF] Gribonval | Harmonic decomposition of audio signals with matching pursuit[END_REF] we denote V def = Span((ψ i ) i∈I ), the closure of vector space spanned by the dictionary (ψ i ) i∈I , V ⊥ its orthogonal complement and we denote the orthogonal projection onto V and V ⊥ by P V and P V ⊥ respectively.

Proposition 5: Let (ψ i ) i∈I be a normed dictionary, v ∈ H and θ(•) be a shrinkage function. The limits defined in Theorem 1 satisfy

+∞ n=0 s n ψ γn -P V v D = R +∞ v -P V ⊥ v D ≤ τ + α ,
where τ + is the external threshold of θ(•), as defined in Definition 2.

Proof: Let ǫ > 0, from Lemma 1, we know that for any k ≥ 0 there is

n k ≥ k inf t:θ(t)=0 t -ǫ ≤ M nk ≤ sup t:θ(t)=0 t + ǫ.
Given the definition of τ + , we therefore know that

-τ + -ǫ ≤ M nk ≤ τ + + ǫ.
We rewrite

|M nk | ≤ τ + + ǫ.
Moreover, since P V is contractive and given the construction of M nk , we know that

|M nk | ≥ α sup i∈I | R nk v, ψ i | ≥ α sup i∈I | P V (R nk v), ψ i |.
Therefore, for all i ∈ I,

| P V (R nk v), ψ i | ≤ τ + α + ǫ α .
Since (R nk v) k∈N converges to R +∞ v (see Theorem 1), we finally have

| P V (R +∞ v), ψ i | ≤ τ + α + ǫ α ,
for all i ∈ I. Since the above inequalities hold for any ǫ > 0, we obtain

P V (R +∞ v) D ≤ τ + α .
Moreover, using Theorem 1, we know that

P V ⊥ R +∞ v = P V ⊥ (v) -P V ⊥ +∞ n=0 s n ψ γn = P V ⊥ (v) .
We therefore obtain

R +∞ v -P V ⊥ v D = P V (R +∞ v) D ≤ τ + α .
Using Theorem 1 (again), we also know that

+∞ n=0 s n ψ γn = P V +∞ n=0 s n ψ γn = P V (v) -P V (R +∞ v) Therefore, +∞ n=0 s n ψ γn -P V (v) D = P V (R +∞ v) D ≤ τ + α .
This finishes the proof of the theorem.

Remark 4: A consequence of the above proposition is that when the MPS is used with a thresholding function, it provides a feasible point for the "Dantzig selector" (see [START_REF] Candes | The dantzig selector: Statistical estimation when p is much larger than n[END_REF]). The "Dantzig selector" consists in the optimization problem:

min (λi)i∈I |λ i | subject to i∈I λ i ψ i -P V v D ≤ τ + α .
From Proposition 3, we know that the MPS provides a set of coordinates (λ i ) i∈I (see [START_REF] Devore | Some remarks on greedy algorithms[END_REF]) such that min Xiaolong (Peking Univ.) for the helpful discussions. We would also like to thank all the anonymous reviewers as their comments and suggestions contributed greatly to the quality of the paper.

APPENDIX

Proof of Proposition 1

Proof of gap(θ) ≤ inf t:θ(t) =0 |t|. Let t 0 ∈ R be such that t 0 > inf t:θ(t) =0 |t|. We cannot simultaneously have θ(t 0 ) = 0 and θ(-t 0 ) = 0, since θ(•) is nondecreasing. Let us denote

t =    t 0 , if θ(t 0 ) = 0 -t 0 , if θ(t 0 ) = 0
We have θ(t) = 0 and given the definition of the gap, we know that gap(θ) |t|.

Lemma used in the proof of Theorem 1

This lemma is a variation on the Lemma used for the proof of Theorem 1 in [START_REF] Mallat | Matching pursuits with time-frequency dictionaries[END_REF]. x k y k < +∞.

One of the following alternatives holds : x k = 0.

Proof: First, since (y k ) k∈N is a sequence of nonnegative real numbers, its inferior limit always exists.

We

• either have lim inf k→+∞ y k > 0,

• or lim inf k→+∞ y k = 0

Let us first assume that lim inf k→+∞ y k > 0.

There exists ǫ > 0 and n > 0 such that for any k ≥ n, y k ≥ ǫ. Therefore, we have As a conclusion, for any ǫ > 0 and any m ≥ 0, there is j ≥ m such that y j j k=0

x k < ǫ.

This means that the second alternative holds.

Fig. 1 .Definition 3 :

 13 Fig. 1. Example of a thresholding function θ. It is non-gap. Its internal and external thresholds are not equal.

(

  λi)i∈I |λ i | is finite. Proposition 5 garantees that the constraint is satisfied. July 8, 2009 DRAFT ACKNOWLEDGEMENT We would like to thank Prof. Alain Trouvé (ENS-Cachan), Prof. Michael Ng (HKBU) and Dr. Li

Lemma 2 :

 2 Let (x k ) k∈N and (y k ) k∈N be two sequences such that∀k ∈ N, 0 ≤ x k ≤ y k(34)and +∞ k=0

  consider ǫ > 0 and m ≥ 0.Since +∞ k=0 x k y k < +∞, there is n ≥ m such that +∞ k=n x k y k < ǫ 2 . (35)Since lim inf k→+∞ y k = 0, there is p ≥ 0 such that j ∈ {n, . . . n + p} be such that y j ≤ y k , ∀k ∈ {n, . . . n + p}. k from (36) and (34)< ǫ from (35).

TABLE I THE

 I

MATCHING PURSUIT SHRINKAGE (MPS).

  , we know that +∞ n=0 |M n ||s n | < +∞. Moreover, 0 ≤ |s n | ≤ |M n | for all n ∈ N. So Lemma 2 (see Appendix) can be applied with x n ≡ |s n | and y n ≡ |M n |. Two situations might occur :

  Let us prove the first statement. If sup t:θ(t)=0 t = +∞ the statement is trivial. We therefore focus on the case sup t:θ(t)=0 t < +∞. Let us assume that (32) does not hold. Then there exists ǫ > 0 and an increasing sequence (k n ) n∈N ∈ N N such that So there exists an increasing sequence (k n ) n∈N ∈ N N such that s kn = θ(M kn ) ≥ θ( sup

	This means that
	lim sup n→+∞	s n > 0.
	and
	inf t:θ(t)=0	t ≤ lim inf n→+∞	M n .	(33)
	Proof: M kn ≥ sup t:θ(t)=0	t + ǫ, ∀n ∈ N.
					t + ǫ) > 0
				t:θ(t)=0
	July 8, 2009				DRAFT

The rigorous definition of shrinkage functions is given in Section II.July 8,

DRAFT

[START_REF] Blumensath | Iterative hard thresholding for compressed sensing[END_REF] DRAFT