Florian Vichot
email: florian.vichot@inria.fr

Hubert Cochet

Benoît Bleuz

Nicolas Toussaint

Pierre Jaïs

Maxime Sermesant

Cardiac Interventional Guidance using Multimodal Data Processing and Visualisation: medInria as an Interoperability Platform

medInria is a free medical imaging software developed at Inria, which aims at providing clinicians with state-of-the-art algorithms dedicated to medical image processing and visualization. Efforts have been made to simplify the user interface, while keeping high-level algorithms. In this particular article, we will concentrate on its use in preoperative preparation for cardiac interventions, and how we handle the difficulties arising from the lack of standard format for data types such as meshes or fibers, the absence of a common programming interface for data processing algorithms, notably registration, and the issues of visualisation where display conventions would be beneficial.

1 Introduction

The ever-increasing presence of technology in the hospital means more and more data has to be visualised, analysed and shared. This situation has raised the need for developing standards between medical actors (equipment manufacturers, software providers, etc). Some success has been achieved with initiatives such as DICOM for medical imaging, or HL7 for medical records, but in other areas such as mesh or tensor data, there is still a lack of accepted standards. Moreover, it is not just data exchange that could benefit from the emergence of standards: the analysis of medical data often involves state-of-the-art algorithms, that are often only available in complex software and that require a rewrite everytime they need to be integrated into a new software, with each time a significant time investment, as well as a risk of errors, or of a loss in performance. Finally, the absence of clear conventions in terms of visualisation is also problematic, as for example there is no standard for the colors scale of lookup tables used to display attributes of mesh files, or the choice of color to represent directions in a tensor image. All of these make preoperative data analysis more complex for clinicians, and represent an increased time cost, as well as an increased risk of errors.

medInria tries to bridge these gaps in different fashions. In terms of data format, our approach is to support a broad range of formats for import and export, and to offer the possibility of exporting into a simple "composite" format that should be easily implemented by other software, while retaining some backward compatibility in software that will not support it. Regarding algorithm re-use, we propose the Registration Programming Interface [START_REF]Registration programming interface[END_REF], which is a generic programming interface for most registration algorithms. We now provide jointly a simple and intuitive graphical interface shared by all registration methods in medInria, that greatly simplifies the user experience while helping the author of a new registration method to focus on the algorithm, not on the user interface. Finally, on the subject of conventions in data visualisation, we discuss potential evolutions of medInria to cope with the varied ways of displaying the same datasets, to ease the work of clinicians in comparing them.

medInria Architecture

Nowadays, software packages for the visualisation and processing of medical images are numerous. So why yet another one? We felt that none addressed both the specific needs of clinicians, by providing an ergonomic, efficient interface, and the needs of researchers, by enabling easy addition of functionalities. While many medical imaging software packages exist, they are all too often designed only for one of these two applications. For example, they often provide a very low-level abstraction to algorithms, which requires some expertise from the users to achieve the desired results. Clinicians needs are simple yet very difficult to satisfy in a computer science point of view: they want ergonomic, reactive and intuitive software that offers real-time interactions with data.

medInria is built in C++ using proven components: for our visualisation needs, we use VTK [START_REF] Kitware | Visualization toolkit[END_REF], the Visualisation Toolkit, an open source, freely available visualisation framework designed for efficient data management and visualisation that is used by many other software, such as 3D Slicer [START_REF]3dslicer[END_REF] or GIMIAS [START_REF] Cistib | Graphical interface for medical image analysis and simulation[END_REF]. The image processing is done using ITK [START_REF] Kitware | Insight toolkit[END_REF], the Insight Segmentation and Registration Toolkit, another open source project. DICOM data is loaded using DCMTK [START_REF]Dicom toolkit[END_REF]. Finally, we use the Qt library for the user interface and the system abstraction, allowing an ergonomic finish to medInria, as well as providing us with cross-platform compatibility.

The software uses a plugin architecture. This allows us to limit the dependencies on the core software, while enabling us to bring in new algorithms and components easily (Fig. 1), as well as shipping version of medInria with just the necessary functionalities to clinicians, to avoid feature bloat.

Latest It is currently at its version 2.0.1, released in April 2012, that builds on the success of the previous versions [START_REF] Toussaint | vtkinria3d: A vtk extension for spatiotemporal data synchronization, visualization and management[END_REF] while also gaining in functionalities. We're currently working on developping an SDK that will allow thirdparties to extend medInria via plugins, to support new data types, new visualisations, new algorithms or even add whole new workspaces. The license of the SDK will be very flexible, allowing either open-source or binary plugins. We plan to make it as simple as possible to write plugins, to encourage the research community to write their own. medInria is freeware for non-commercial use, but we are considering eventually releasing it as open-source.

Data Formats

The first issue we have to tackle when we want to provide a software usable by clinicians during the preoperative phase of an intervention, is to be able to import all the necessary data he will need to establish a diagnosis and identify the pathological area to ablate, for example. We then have to provide means to export the result of that preoperative process in the right format, so it can be loaded into the equipment used during the intervention. Importing medical images is mostly a solved problem: the existence of standards such as DICOM means that there are robust open-source implementations of image importers such as DCMTK that we can use, and we can also use ITK image readers for a large array of image formats.

There are still some rough edges though, for example, it is known that some equipments and software generate the orientation matrix of images in an incorrect manner, thus medInria comes with its own heuristics to try to identify those cases and fix them transparently for the user.

Clinicians also need to import other types of data beyond medical images, such as meshes (see Fig. 2) or tensor images. In the specific case of meshes, depending on their origin, they can come in very different formats. In our example of cardiac electrophysiology, they frequently come from Saint-Jude NavX system, BioSense CARTO, or CardioInsight, with various file formats. Currently, in research, the de facto standard for meshes is the VTK format.

Sadly, the fact that the VTK file format [START_REF] Kitware | Vtk file formats[END_REF] is the de facto standard when it comes to meshes is not a satisfactory solution. This new format is based on the concept of "Composite dataset" that currently exists in medInria: similar to formats like the Open Document Format[7], we wrap a collection of files in a ZIP archive, and expose this to the user. It makes file sharing easier as everything is exposed as a single file, making it simpler to keep it all in one place.

Inside the ZIP archive, we store a sequence of meshes in a simple and backward compatible manner, that is as a series of VTK files, thus giving the possibility to users who wish to load those file into software that does not implement this format to unZIP the archive and load the file separately. We can then use regular text files stored inside the archive to implement functionalities like metadata or annotations. As the format in these text files is kept simple and human readable, their content can easily be imported, either by hand or through an ad-hoc importer, into other software.

Processing Algorithms

When we talk of standards we often think of data formats and protocols, but we should consider also the issue of sharing algorithm implementations amongst different software. As clinicians often require a large amount of data to be able to make a diagnosis, different modalities of observation may be needed in order to gather information on the patient. These different modalities may have different origins, either because they were acquired at different times or in different spatial configurations, which justifies the need for algorithms which registers in the same spatial coordinates different datasets.

However, pairing the numerous image modalities and file formats around with the many registration algorithms is far from being straight forward. Indeed, each registration algorithm has its own input/output formats, set of parameters and visualisation. Thus, its use by a non-expert is challenging and the combination or comparison of such algorithms is difficult.

medInria has tried to address this concern by introducing an open source registration API which goal is to propose a simple and intuitive interface shared by all registration methods. This common interface greatly simplifies the user experience while helping the author of a new registration method to focus on the algorithm, not on the user interface. We call this interface RPI for Registration Programming Interface.

ITK also provides a registration framework and one may ask why RPI is not included into this framework.

The ITK registration framework is too specific and does not allow to easily integrate some registration methods that may not be ITK based, and which do not use its registration loop. As a consequence, it has been decided to not base the registration API on this framework but to design a new framework simpler and more generic. It is still partly based on existing ITK objects though. We made that choice as we felt ITK being a powerful and popular open-source library, it was best not to reinvent already proven components.

RPI defines a C++ programming class which wraps the registration algorithm and provides the generic API.

Latest version available at the Insight Journal [ http://hdl.handle.net/10380/3361] Distributed under Creative Commons Attribution License This class gives a developer a strict but simple programming framework, that exposes only a few methods to override to connect the algorithm to the API. Its purpose is to provide a common interface which has been designed to cover most of the registration methods. It handles many different transformations : linear (both rigid and affine), displacement fields, and stationary velocity fields. The two images to register together are stored as itk::Image objects. This allows a wide range of image types and dimensions : pixel type can be either a scalar type (e.g. short, float, etc.) or can be a itk::Vector for diffusion tensor images. The resulting transformations are also stored as ITK objects.

The registration algorithm can be tuned using parameters that can be exposed by the developer via the API.

For the moment, this mechanism is rather simple, so it does not lend itself to automatic GUI generation for example, as parameters have no semantic typing or metadata associated with them yet. This is an improvement currently being considered, allowing us to associate hints such as bounds, units, acceptable values etc to parameters.

The RPI effort is open-source, and its documentation and source are available online [START_REF]Registration programming interface[END_REF]. format of a given interventional platform, and load it in the cathlab. We present in Fig. 4 the meshes generated by the commercial platform on the left, and the one generated by medInria, once loaded into the platform on the right.

Data Visualisation

In the context of cardiac preoperative diagnosis, medInria is also often used by clinicians to study data such as activation maps or potential maps, to identify reentrant circuits that may cause tachycardia, for instance.

To do so, they load mesh files with their associated attributes, and choose a color scale to map values into colors. The usual technique for applying a color map to a mesh's attribute is to compute its range, and "stretch" the values to cover the whole color scale. The issue with this approach is that to compare two meshes loaded in the same view, you need to have the same value mapped to the same color, so in medInria we make sure that the range is computed on the values of all the displayed meshes combined, to make sure the color maps applies the same to each. This also applies to sequences, where care has to be taken to compute the range over the whole sequence, and not over individual time steps.

The second issue with color scales, is that there is no actual standard set of color scales, nor a standard format to import / export them. In medInria, we have taken the practical approach of providing the set of color scales found in the most used imaging software, such as the one found in Osirix.

A similar problem arises concerning the display of diffusion tensors directions (Fig. 5) using color scales.

Figure 5: Tensors display with colors representing the direction of the diffusion

There is to our knowledge no standard as to which color should stand for which axis, so it's common that different software use different colors, and thus have very varying visualisation of the same dataset. Again, in this situation we try to take a practical approach, and we provide a way to select a color for each axis.

Latest version available at the Insight Journal [ http://hdl.handle.net/10380/3361] Distributed under Creative Commons Attribution License

Conclusion

medInria strives to provide clinicians with a reactive, ergonomic and intuitive software, which finds itself at a crossroad between trying to provide access to state-of-the-art while also staying grounded on practical considerations. Thanks to our partnership with the CHU of Bordeaux, we get direct and concrete feedback on what make sense in a live environment, allowing us to refine our proposed solutions, as well as design new tools that fill an actual need of clinicians.

Given the interoperability issues between different equipments and the difficulty of interaction between competing companies, we believe such academic software is needed in order to achieve the integration of multimodal datasets coming from various equipments for diagnosis, therapy planning, and interventional guidance.

medInria is growing fast, and will remain free. Windows, Linux and Mac OS X versions can be downloaded at http://med.inria.fr.

Figure 1 :

 1 Figure 1: General architecture, here demonstrating an image reader, a view, and a registration plugin working together. (inputs are green, output red)

Figure 2 :

 2 Figure 2: Combining data from different sources: a CT-scan image in volume rendering mode, an Ensite NavX mesh and a BioSense CARTO mesh. Such multimodal integration is not possible in the commercial platforms of each of these companies.

Figure 3 :

 3 Figure 3: medInria v2.0.1: the registration workspace

Figure 4 :

 4 Figure 4: Two heart meshes. Left: a heart mesh from the NavX System, Right: a heart mesh segmented from images, exported to NavX with medInria


Islem Rekik
email: islem.rekik@gmail.com

Stéphanie Allassonnière
email: stephanie.allassonniere@polytechnique.edu

Olivier Clatz
email: olivier.clatz@inria.fr

Ezequiel Geremia
email: ezequiel.geremia@inria.fr

Erin Stretton
email: erin.stretton@gmail.com

Hervé Delingette
email: herve.delingette@inria.fr

Nicholas Ayache
email: nicholas.ayache@inria.fr

Tumor Growth Parameters Estimation and Source Localization From a Unique Time Point: Application to Low-grade Gliomas

Keywords: diffusivity ratio, source estimation, Eikonal Equation, reaction-diffusion glioma growth modeling

Coupling time series of MR Images with reaction-diffusion-based models has provided interesting ways to better understand the proliferative-invasive aspect of glial cells in tumors. In this paper, we address a different formulation of the inverse problem: from a single time point image of a non-swollen brain tumor, estimate the tumor source location and the diffusivity ratio between white and grey matter, while exploring the possibility to predict the further extent of the observed tumor at later time points in low-grade gliomas. The synthetic and clinical results show the stability of the located source and its varying distance from the tumor barycenter and how the estimated ratio controls the spikiness of the tumor.

Introduction

Brain gliomas represent about 50% of all primary brain tumors [START_REF] Tovi | MR imaging in cerebral gliomas analysis of tumour tissue components[END_REF] and can be classified according to their grade of malignancy. Low grade gliomas (LGG) are slow invaders of brain tissue as they keep growing for many years, presenting one of the most controversial decision treatment areas. High grade gliomas (HGG) remain unfortunately incurable with an average life expectancy of one year after its discovery, eventually creating symptoms due to an increase of the intracranial pressure or swelling around the tumor. The diagnosis of brain gliomas includes the analysis of various MRI sequences of the brain which partially reveal the tumor invasion. Based on those images and other clinical information, neurologists try to determine the grade of the gliomas and to estimate their current and further spatial extent and if possible their source location.

For more than a decade, mathematical models of brain tumors have been devised to help clinicians answer these questions. Microscopic models study the cellular mechanisms [START_REF] Bearer | Multiparameter computational modeling of tumor invasion[END_REF][START_REF] Frieboes | Computer simulation of glioma growth and morphology[END_REF][START_REF] Sanga | Predictive oncology: a review of multidisciplinary, multiscale in silico modeling linking phenotype, morphology and growth[END_REF][START_REF] Wang | Prognostic significance of growth kinetics in newly diagnosed glioblastomas revealed by combining serial imaging with a novel biomathematical model[END_REF] that explain the growth dynamics of gliomas at a microscopic scale. On the other hand, macroscopic models pioneered by Murray [6] describe the evolution of tumor cell density. However, those quantities cannot be directly observed in clinical medical images, but it is commonly assumed that visible tumor boundaries correspond to an isovalue of this density. More realistic reaction-diffusion models have been proposed by Swanson et al. [START_REF] Swanson | A quantitative model for differential motility of gliomas in grey and white matter[END_REF] based on the fact that tumor cells migrate faster on white matter fibre tracts myelin sheaths [START_REF] Giese | Migration of human glioma cells on myelin[END_REF]. They have been further refined by Jbabdi et al. [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF] and Clatz et al. [START_REF] Clatz | Realistic simulation of the 3d growth of brain tumors in mr im-ages coupling diffusion with mass effect[END_REF] by considering an anisotropic diffusion in the white matter whose diffusion tensor is estimated from those acquired in DT-MRI.

A key issue for those models to answer clinical questions is their personalization, i.e. the estimation of some patient-specific parameters from medical images. The main parameters to be identified based on reaction-diffusion models were pointed out in [START_REF] Angelini | Glioma dynamics and computational models: a review of segmentation, registration, and in silico growth algorithms and their clinical applications[END_REF] as a combination of tumor diffusion tensors in white and grey matter, its proliferation rate, its initial point and its initial time. Several authors have estimated patient-specific parameters manually [START_REF] Clatz | Realistic simulation of the 3d growth of brain tumors in mr im-ages coupling diffusion with mass effect[END_REF] or through major model simplifications [START_REF] Swanson | Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion[END_REF]. Colin et al. in [START_REF] Colin | System identification in tumor growth modeling using semi-empirical eigenfunctions[END_REF] used a reduced model based on Proper Orthogonal Decomposition (POD) in order to identify growth parameters of pulmonary nodules in CT images. Konukoglu et al. [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF][START_REF] Konukoglu | Towards an identification of tumor growth parameters from time series of images[END_REF][START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] have proposed an approach to automatically and accurately personalize brain tumor models. They first remarked that given time series of brain MR images, only the motion of a tumor front can be observed and therefore only the three following parameters can be recovered : diffusivity in the white d w and grey matter d g as well as the initial time T 0 . Furthermore, since the tumor cell density are only observed in MR images through visible tumor boundaries corresponding to an isovalue, the reaction-diffusion equations can be advantageously replaced by an Anisotropic Eikonal Equation (AEE) [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] which models the time at which a tumor front reaches a given point. By minimizing the distance between the segmented tumor and the simulated one, they were able to estimate uniquely those three parameters and test the prediction of future tumor evolution from at least a pair of images.

In this paper, we tackle a slightly different problem than the one approached by previous authors. Instead of estimating the speed of tumor growth from a time series of images, we aim at characterizing the nature of the glioma, more precisely LGG, from a single MR image. Indeed, we hypothesized that the tumor shape is dependent on the proliferating or infiltrating nature of the tumor. Contrary to HGG where the presence of brain edema is common and usually associated with tumor malignancy, LGG are slowly growing tumors with a minimal surrounding edema [START_REF] Deangelis | Brain tumors[END_REF][START_REF] Kaal | The management of brain edema in brain tumors[END_REF][START_REF] Whittle | The dilemma of low grade glioma[END_REF]. Since our methodology main focus is on LGG, we will not consider the edema-induced mass effect in our further formulation of the tumor growth model. Therefore, the anatomical boundaries such as the ventricles' will remain static as the tumor grows.

Given a segmented brain glioma from an MR image, we solve an inverse problem in order to estimate the diffusivity ratio d w /d g and the tumor source position. By localizing the tumor source and estimating the invaded tissue characteristic using this ratio d w /d g our objective is to provide clinicians with new indices that can be used for diagnosis from the first acquired MR images, combined with a subsequent prediction of tumor invasive margins as it grows from the initially observed boundary. This additional information may help in surgical and/or radio-therapeutic treatment planning especially when it comes to determining the margins for applying the therapy. The problem of tumor seed localization was recently raised in [START_REF] Gooya | Deformable registration of glioma images using em algorithm and diffusion reaction modeling[END_REF], where a reaction-diffusion based joint estimation of tumor evolution parameters was combined with a multimodal deformable registration framework. This approach focuses on the MR image registration with an atlas providing the estimation of the initial seed location as a by-product. Later on, it was extended to a joint segmentation and deformable registration between multi-sequence brain MR image with gliomas and a probabilistic brain atlas of a healthy population [START_REF] Gooya | Joint segmentation and deformable registration of brain scans guided by a tumor growth model[END_REF][START_REF] Gooya | GLISTR: Glioma image segmentation and registration[END_REF]. These recent publications targeting the localization of tumor seeds, the quantification of the mass effect and its white matter diffusion coefficient was based on a normal atlas.

In this paper, we investigate a different problem from registration or segmentation in low-grade gliomas. Our approach addresses a spatio-temporal tumor evolution with the estimation of the diffusivity ratio and the position of its source. Additionally, after solving the inverse problem, we evaluate the stability of the location of tumor sources over time and analyze the relationship between the tumor source and tumor barycenter as both locations have been assumed to match in past studies [START_REF] Laigle-Donadey | Correlations between molecular profile and radiologic pattern in oligodendroglial tumors[END_REF][START_REF] Drabycz | An analysis of image texture, tumor location, and mgmt promoter methylation in glioblastoma using magnetic resonance imaging[END_REF][START_REF] Anderson | A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion[END_REF][START_REF] Sottoriva | Cancer stem cell tumor model reveals invasive morphology and increased phenotypical heterogeneity[END_REF][START_REF] Kansal | Emergence of a subpopulation in a computational model of tumor growth[END_REF]. Finally, knowing the diffusivity ratio d w /d g and the tumor source from a single image, we evaluate whether this information gives insights into predicting further tumor shape evolution in two distinct cases of low-grade gliomas.

Material and Method

MR Glioma images

In this work, we assume that for a given patient, one MR FLAIR image of the brain has been acquired, showing visible boundaries of glioma cells. We also assume that a Diffusion Tensor MR Image (DT-MRI) is also available and acquired at the same time as the FLAIR image. While FLAIR images are acquired in routine on patients with brain tumors, this is unfortunately not the case for DT-MRI. The extent of the tumor has been manually segmented in FLAIR images. Similarly, brain masks have been manually delineated on those images from which brain ventricles have been removed by a simple thresholding of the signal. Also white matter regions have been isolated by thresholding the voxels with the largest eigenvalue in the DT-MRI.

Objective: From the segmented tumor in the FLAIR image, our objective is to provide a quantitative analysis of the tumor shape which is not simply based on geometry (spheroid vs star-shaped) but based on simple biophysics growth principles. Indeed, two quantities are estimated in this analysis: the tumor source position and the diffusivity ratio between white matter and grey matter. This information may provide additional hints about the nature and the future progression of the tumor.

Data issues: Collecting LGG data with DT-MRIs is not a very straightforward task since diffusion MR is a quite recent technology [START_REF] Young | Advanced mri of adult brain tumors[END_REF], rarely used in the common clinical practice and furthermore acquired DT-MRI may have various anomalies like holes, low resolution and signal distortion. This is particularly true around any tumor lesion. This lack of information may be compensated by assuming the symmetry of the brain. For this, we perform a symmetrization process to "reinitialize" the region where the tumor grew and induced a diffusion signal distortion. Hence, for a proper simulation of tumor growth, we have corrected the tensor field by making the hypothesis that the DT MRI was originally symmetric with respect to the mid-sagittal plane. Thus the DT-MRI voxels corresponding to the largest extent of the tumor have been modified by symmetrizing and copying the voxels from the healthy brain hemisphere. However, this symmetrization process is prone to the following difficulties: a-the DT MRI is asymmetric especially in white matter where about 50% of the contralateral tumor volume has zero Fractional Anisotropy (FA) values, while in the corresponding affected region where the tumor grew, the diffusion signal exists with a remarkable distortion (see Fig 1-a where the FA tumor map is darker than the contralateral part with a significant absence of symmetry).

b-even with a reliable DT-MRI symmetry, about 50% of the contralateral non-pathological symmetric volume to the tumor volume has no DT (i.e FA) signal. In Fig 2-a, the presence of large black hole in the FA signal in the contralateral part to the tumor invasion area presents a major barrier to diffusion tensor-guided glioma evolution simulation. To cope with this DT signal-missing problem, an interpolation algorithm based on isotropic diffusion (solving the heat equation) was applied to estimate the missing tensors from neighboring regions. This is done by applying a Gaussian convolution separately on the six components of the diffusion tensors.

Dataset: By excluding LGG cases with completely distorted or missing DTI, we succeeded to include four LGG patients. The first case, a.k.a patient A, has developed a second grade astrocytoma classified as a low-grade glioma. Four successive time points of T2 flair MR images with a resolution of 0.99 × 1 × 2.16 mm 3 were acquired and visible tumor boundaries were manually delineated by an expert. A DT-MRI image was also acquired at the first acquisition time point. As expected, the white matter fiber tracts are perturbed by the tumor growth and the DT-MRI signal near the tumor does not capture the original diffusion tensors of the brain at the onset of the disease. Therefore, we used the symmetrization process to reset the affected area to non-pathologically spatially deviated tensors (see Fig 3 where the white matter diffusion tracts were beautifully recovered).

The second case, a.k.a patient B, suffers from a low-grade glioma and four MR images were acquired at distinct time points with a resolution of 0.89 × 0.97 × 0.97 mm 3 . Only one DT-MRI was acquired during the initial scanning process and the tumor region of the DT-MRI has been symmetrized similarly to patient A. In addition, the DT-MRI includes small holes in the opposite region to the tumor that were successfully interpolated.

The third case, a.k.a patient C, with LGG has three acquisition time points with a resolution of 2 × 2 × 2 mm 3 . In this particular case, displayed in Fig 1 , we encountered difficulty (-a-) where the DT-MRI is not fully symmetric.

The fourth case, a.k.a patient D, with a resolution of 2 × 2 × 2 mm 3 was not symmetrized due to the large holes in the symmetric region to the tumor as show in Fig 2-a. Therefore, we interpolated the holes in the affected tissue without the symmetrization process.

Tumor Growth Modeling : from reaction-diffusion to Eikonal equations

Glial cells dynamics are essentially governed by two biological phenomena : proliferation and invasion. They can be jointly modeled by a reactiondiffusion equation which describes the change over time of the normalized tumor cell density u:

∂u ∂t = ∇. (D(x)∇u) + ρu(1 -u) D ∇u. n ∂Ω = 0 ( 1
)
where ρ is the proliferation rate, D the local diffusion tensor, and n ∂Ω is the normal vector at the domain boundary surface ∂Ω. In the first equation, the proliferation of tumor cells follows a logistic growth parameterized by ρ whereas the tumor infiltration into neighboring neural fibers is captured by an anisotropic diffusion parameterized by D. The second equation indicates that there is no flux of tumor cells outside the domain Ω.

The diffusion tensor is a definite positive and symmetric 3 × 3 matrix whose value may be linked to Diffusion Tensor MRI (DT-MRI) [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF]. Indeed, it characterizes the motility of tumor cells that is considered to be isotropic in grey matter but anisotropic in white matter. More precisely, the tumor diffusion tensor (TDT) may be written as D(x) = d g I 3 in grey matter, where d g is the diffusivity coefficient.

In white matter, there are several approaches to link the TDT D(x) with the DT-MRI signal D water (x). Clatz et al. [START_REF] Clatz | Realistic simulation of the 3d growth of brain tumors in mr im-ages coupling diffusion with mass effect[END_REF] proposed to have D(x) proportional to D water (x) whereas Jbabdi et al. [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF] have introduced a formulation which takes into account the possible equality of the two largest eigenvalues corresponding to a possible fiber crossing. Due to high anisotropies of D water (x) in most parts of the white matter, these two approaches lead however to diffusivities that are much lower than d g in the directions orthogonal to the fibers, which is questionable. Furthermore, the high ratios of anisotropy encountered in those two expressions also lead to large computational times.

In this paper, we propose to use the following white matter tumor diffusion tensor:

D(x) = V(x)[diag(e 1 (x)d w , d g , d g )]V(x) T (2)
where d w is the white matter diffusivity coefficient, V(x) represents the matrix of sorted eigenvectors of D water (x) and e 1 (x), is the normalized largest eigenvalue (between 0 and 1) of D water (x). With this choice, tensors have a non-homogeneous anisotropy ratio which is always less than d w /d g but is maximized at the center of the white matter fibers and continuously decreases towards their boundaries. By simply dividing the duration of tumor evolution simulation of its propagating front using our adopted diffusion tensor 2 by the simulation duration as we used the diffusion tensor formulas presented in [START_REF] Jbabdi | Simulation of anisotropic growth of low-grade gliomas using diffusion tensor imaging[END_REF] (we precisely refer the reader to formulas A11 and A12), we have noticed that our choice speeds-up the computational time by a factor of 200 without any significant differences in performance. In fact, the use of a more nonlinear (more anisotropic) diffusion tensor field increases the computational time of the solution as the characteristic direction of the recursive anisotropic fast marching algorithm used to solve the AEE 5 becomes harder to find. This also may be explained by the fact that below a certain anisotropy ratio, the difference in tumor growth simulation is hardly noticeable. The reaction-diffusion equation ( 1) is not practical when dealing with clinical images. Indeed, in MR images tumor cell density u cannot be observed but only the visible tumor boundary can. Hence, a front motion approximation for the reaction-diffusion equation was introduced by Konukoglu et al. [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] assuming that the visible contour is associated with iso-density contour u = 0.4 [START_REF] Tracqui | A mathematical model of glioma growth: the effect of chemotherapy on spatio-temporal growth[END_REF]. They introduced an Anisotropic Eikonal Equation (AEE) describing the time T (x) at which the evolving tumor front passes through the location x. In its simplest form, the AEE writes as:

F √ ∇T T D∇T = 1 with F = 2 √ ρ (3)
However, such approximation of equation 1 is too simple and Konukoglu et al. then proposed to account for the fact that the tumor front speed increases over time to reach an asymptotic value equal to 2 ρ n T Dn where n is the normal direction of the front. The time-dependent speed of the propagating tumor front was introduced in [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] as:

F = 4ρT -3 2 √ ρT (4)
Furthermore, the front curvature κ ef f (x) also plays a role in the front speed as the front slows at high curvature points. This is especially important at the early stage of the tumor growth when the front is similar to a small sphere. With this additional hypothesis, the speed term becomes :

F = 4ρT -3 2 √ ρT -0.3 √ ρ 1 -e -|κ ef f |/0.3 √ ρ (5)
This last formulation is no longer a Hamilton-Jacobi equation due to the second-order curvature term and therefore cannot be solved by fast sweeping methods such as the Anisotropic Fast Marching (AFM) [START_REF] Konukoglu | A recursive anisotropic fast marching approach to reaction diffusion equation: Application to tumor growth modeling[END_REF]. However, a multi-pass approach was proposed [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] to solve efficiently this equation by applying several times the AFM method while estimating the curvature front from previous iterations. The AFM method is recursive and the larger the tensor anisotropy the more iterations are needed to compute the characteristic direction of equation ( 5). Our white matter TDT of equation ( 2) limits the anisotropy ratio and therefore leads to reasonable computational times (typically few minutes for a tumor growth from a seed point).

Parameter estimation problem from a unique MR image

Based on the previously exposed mathematical model, we can simulate the growth of a glioma given its initial source S(x) for which we assume that T (S(x)) = 0. From this boundary condition and the knowledge of diffusivity d w , d g and proliferation rate ρ we can compute the time T (x) at which the visible tumor front reaches a given point. The isocontours of the field T (x) correspond to the successive shapes of the visible tumor boundary over time as shown in Fig 4 . The speed on the front is not constant but its asymptotic value is 2 ρd g in grey matter and 2 √ ρd w in white matter. In this paper, we are interested in solving the following inverse problem: given a visible tumor boundary S Seg in an MR image, can we extract the growth parameters d w , d g , ρ, T Obs and source location S(x) that best explain the observed tumor boundary ? The duration T Obs between the onset of the tumor and the MR image acquisition is indeed also unknown.

Based on [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF], it has been already established that several combinations of ρ, d w , d g lead to the same front speed and therefore the same tumor growth simulations. Therefore, it is sufficient in this inverse problem to consider a fixed value of the proliferation rate ρ corresponding to the tumor grade and to estimate the remaining parameters. However, unlike [START_REF] Konukoglu | Towards an identification of tumor growth parameters from time series of images[END_REF][START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] this problem can be further simplified by realizing that the front speed cannot be estimated since T Obs is unknown. If one multiplies the diffusivities by a scale factor α then one obtains the same isocontours for a propagation time divided by √ α 1 . This means that the simulated tumor isocontours do not depend on that absolute value of d g and d w but on the diffusivity ratio:

r = d w d g
In the remainder, we will show that this ratio is related to the spikiness of the tumor.

A simple sensitivity analysis has led to conclusion that solving the inverse tumor growth problem only depends on the following 2 parameters: the source location S(x) and the spikiness index r. Here, the "spikiness index" represents a biology-driven estimated measure which quantifies the tortuousness of the boundary of the tumor as displayed on MR axial slices (see Figure 4-A), an index related to the frequency at which the tumor shape bends and twists. We consider that {S(x), r} appropriately characterize well a tumor extent if its visible boundary in MRI, S Seg , is an isocontour of the simulated tumor growth initiated at S(x) with diffusivity ratios equal to r. Therefore, we propose to estimate the patient specific parameters by minimizing the following criterion:

C ρ (S(x, y, z), r) = 1 NT x∈S Seg (T (x) -T ) 2 (6)
with

T = 1 N x∈S Seg T (x) ( 7
)
where N is the number of points belonging to the manually delineated tumor boundary S Seg . In these equations, T and C ρ are respectively the mean time value and time standard deviation computed over the tumor boundary. Our motivation to use this criterion to get good estimates of our unknowns (S(x, y, z), r) derives from the fact that a tumor boundary (propagating front) is simultaneously defined as an isotime and an isosurface. Thus, to quantify how good is the estimation of the parameters guiding the spatiotemporal evolution of the tumor shape, we need to quantify how closely the simulated isosurface matches the observed one (manually delineated boundary). From a time perspective, this also implies that when measuring the time T at every point x of the manually delineated tumor isotime boundary, its value T (x) will be constant in the best case scenario where the simulated tumor front exactly matches the observed tumor boundary. Therefore, minimizing the time standard variation criterion C ρ over the delineated tumor boundary S Seg will guide the algorithm towards a better estimation of the two key parameters driving the invasive tumor front into fitting the successive MR observed boundaries. Note that C ρ is normalized by T because the criterion should be made independent of the tumor front speed and therefore the mean time T .

In order to efficiently minimize the previously outlined criterion, we use the multidimensional unconstrained minimization algorithm without gradient introduced by Powell in [START_REF] Powell | The bobyqa algorithm for bound constrained optimization without derivatives[END_REF]. This algorithm suits our case since our parameters are bounded in both biological and geometrical spaces. Moreover, the derivative of minimization criterion C ρ is not easy to compute. To better study the convergence of this algorithm and evaluate its outcome, 15 tests were performed using synthetic tumors. Further evaluation of this method was then studied using real data: two patients with LGG.

Synthetic Data: synthetic tumor generation process

In order to validate our parameter estimation, we produce synthetic tumor MR images, where the initial tumor location and diffusivity ratio are known. The procedure of generating synthetic tumors relies on choosing a seed point S 0 in either white or grey matter and the proliferation rate ρ to one of the following values 0.008, 0.012, 0.024/day. The tumor front propagation is simulated with a white matter diffusion rate d w and grey matter diffusion rate d g whose ratio varies between 1 and 100 (considered to be a biologically valid range) using equation [START_REF] Wang | Prognostic significance of growth kinetics in newly diagnosed glioblastomas revealed by combining serial imaging with a novel biomathematical model[END_REF]. The growth simulation is stopped at a specified time T Obs thus leading to a time distance map as seen in Fig 4 . Different synthetic tumors were created at different anatomical locations with different sizes and asymptotic speeds of growth in both white and grey matter.

Results

Convexity of the minimization function: a convergence study

To check the convergence of the Powell minimization algorithm to the initially set parameters, we study the convexity of the minimization criterion C (see Equation ( 6)) for the four scalar parameters [S x , S y , S z , r], writing separately the coordinates of the tumor source location. We proceed by alternatively fixing some parameters to their ground truth values and optimizing the remaining ones with the proposed minimization process. In Fig 5, we can clearly see the convexity of the minimization surface plotted after setting the diffusivity ratio r and one of the source coordinates to their right values.

The convexity of the minimization criterion C ρ was successfully checked also when optimizing three parameters and setting the fourth one at its true value. The fact that the minimization criterion appears to be convex at the vicinity of the ground truth parameters is reassuring about the observability of the four parameters. However, the functional may still have local minima and practical optimization results will be discussed in the next sections.

Synthetic data

We evaluated our method on 15 synthetically generated tumors with a diffusivity ratio ranging from 1 to 100. The algorithm succeeded to identify the original tumor source with a mean error of 0.42mm and a standard deviation of 0.36mm. Moreover, it always converges to the real diffusivity ratio d w /d g with a mean error equal to 0.18 and a standard deviation of 0.06. Furthermore, we use this synthetic data to provide a better understanding of the diffusivity ratio r. The first row of Fig 6 shows how the shape of the tumor boundary can switch from "sphere-like"to "star-like"as the diffusivity ratio value jumps from 1 to 50. Although the displayed tumors in the third row are of different sizes, the correlation between tumors' shape spikiness and the value of the diffusivity ratio r ∈ [1 -100] is clear.

Besides, we also notice from the second row of Fig 6 similar tumor contours in terms of spikiness and irregularity when fixing the diffusivity ratio at the same value and varying the location of the initial tumor seed. This confirms that the diffusivity ratio value r controls the spikiness of tumor shape.

Clinical data 3.3.1. Evaluation criteria

In clinical data, no ground truth values of the parameters are available and therefore additional criteria must be introduced in order to assess the quality of the parameter estimation.

After the optimization of criterion 6, we extract the simulated isocontour which is closest from the visible tumor boundary in MRI. Computing symmetric distances between the two surfaces provides a quantitative information about how well the tumor shape can be explained by the proposed tumor growth model. We detail below the proposed approach.

We extract the closest isocontour defined by time T Obs by optimizing C isoT ime (50%)(T ) :

T Obs = min T C isoT ime (50%)(T ) ( 8 )
The criterion C isoT ime (50%)(T ) is defined as the median symmetric distance between the visible tumour boundary S Seg and the isosurface ŜisoTime at time T . More precisely, for each voxel of ŜisoTime (resp. S Seg ) its closest distance from S Seg (resp. ŜisoTime ) is computed through a distance map and added to a list. The median value of that list is then taken as C isoT ime (50%)(T ). The optimization of the functional is done with the Powell algorithm [START_REF] Powell | The bobyqa algorithm for bound constrained optimization without derivatives[END_REF] already used for solving the general inverse problem and requiring few estimations of the functional. The range of time value for the optimization is constrained to be in the range [Tδ; Tδ] where δ is the standard deviation:

δ = x∈S Seg (T (x) -T ) 2

N

We also compute other robust distance criteria C isoT ime (y%)(T ) by taking the y% quantile of the symmetric distances.

Once the closest isochrone surface ŜisoTime is estimated, we define the symmetric distance between (S, ŜisoTime ) as the first evaluation criterion: We can also measure the similarity between the visible tumor volume V seg and the simulated tumor volume V isoT ime :

C symDist (S Seg , ŜisoTime ) =
C overlap (V seg , V isoT ime ) = V seg ∩ V isoT ime 0.5 (V seg + V isoT ime ) (10)

Study of Patient A

For patient studies, we take advantage of the fact that several FLAIR MR images are available for different time points. We first check the stability of source locations for the M time points per patient and correlate the spikiness index with the observed tumor evolution. To estimate the two target parameters [S(x), r], the proliferation rate was fixed to 0.012/day. We can clearly observe in Table 1 the stability of the tumor source location and the fact that for the first 3 time points this location is very close to the tumor barycenter since the Euclidean distance between them d(G, S) ranges from 1.42mm to 2mm. Nonetheless, for the fourth time point, the estimated tumor source distance from its barycenter is 5.38mm. Moreover, the diffusion rate r is mainly constant. In fact, patient A showed a very slow and quasistatic evolution of the tumor boundary between the second and the third acquisition time points. Our simulation model depicts this static evolution behavior as the corresponding estimated spikiness index r shown in Table 1 remains quasi-static (r = 40.47 at t 2 and 40.58 at t 3 ), which is in line with the identical appearance of the estimated tumor red boundaries shown in the first and second rows of Figure 7. However, it is worth noting that the algorithm failed to catch the tiny bump that appeared in the tumor boundary at t 3 (see first column, row 1 and 2 in Figure 7). This can be explained by the fact that the algorithm is based on a mean error minimization process, which therefore accounts for a mean morphological change in tumor shape over all the axial slices included in the tumor volume.

C
Second, we study the possibility to predict further spatial tumor boundaries evolution based on parameter evaluation. More precisely, after estimating patient specific parameters at time t i , we run a forward simulation using those estimated parameters trying to predict the tumor invasion process as time evolves. We then proceed by estimating the closest predicted isocontour from the observed tumor boundary at time t i+k as described in section3.3.1. 2 display the outcomes of the evaluation criteria for every time point.

Fig 8 and Table

We notice a limited discrepancy between the observed contour at time point t i and the predicted one at t i+1 since the symmetric distance between these contours C symDist ranges from 1.74 mm to 1.97 mm. Moreover the median distance error C isoT ime (50%) ranges from 1.44 mm to 1.68 mm and the volume overlap C overlap between the real tumor volume and the extracted one ranges from 85.63% to 88.04%. It is remarkable that the estimated parameters at a specific time point t i captures the tumor spatial evolution at the successive time point t i+1 . As highlighted in Fig 8 ,the estimated parameters based prediction shows promising results.

Study of Patient B

The parameter estimation outcome for Patient B at its four time points is summarized in Table 3. The obtained results confirm the stability of the tumor source location. At every time point, the distance between the tumor barycenter and the located source exceeds 4 mm, meaning that the tumor source may not be close to its barycenter. Unlike Patient A, the spikiness index value r increases as time evolves as it can also assessed visually in Fig 9.

Despite the poor signal quality of the DT-MRI, the evaluation criteria confirm that the estimated parameters were adequate to spatially predict the tumor front propagation. In fact, the symmetric distance C symDist between the real front observed at t i and the extracted one at t i+1 ranges from 1.13 mm to 1.53 mm. The volume overlap C overlap ranges from 82.20% to 88.02% (Table 2). The median distance error C isoT ime (50%) ranges from 0.92 mm to 1.16 mm (Fig 8). Fig 10 shows that the algorithm succeeded to spatially capture the geometric behavior of the evolving tumor. However, the maximum distance error C isoT ime (100%) reaches high values ranging from 7.10 mm to 11.97 mm. These outliers are related to the interpolated holes in the DT-MRI where the simulated propagating front fails to capture the real tumor invasion process.

Study of Patients C and D

In both patient C and D, the estimated source remained stable from the first time point to the third one as the distance between the successively estimated sources d(S t 1 , S t 2 ) and d(S t 1 , S t 3 ) did not reach a value over 1.42mm. The estimated diffusivity ratio for Patient C remained constant as time evolves (r = 7.92), a spikiness value in concordance with the round shape of the tumor boundary and its very slow evolution between t 1 and t 3 as shown in Fig 11. The spatial disparity between the estimated source and the barycenter did not exceed 1.8mm. For Patient D, the estimated diffusivity ratio varied from r t 1 = 20 to r t 2 = 16.2 then increased to r t 3 = 23.8. The successively estimated high diffusivity ratios are in line with the tumor boundary spikiness shown in Fig 2-b. The computed distance between the estimated source and the tumor barycenter varied from 1.42mm at t 1 to 3mm at both t 2 and t 3 .

Discussion

The patient-specific parameter estimation method presented in this paper uses two main ingredients: the traveling front approximation presented by Konukoglu et al. [START_REF] Konukoglu | Image guided personalization of reaction-diffusion type tumor growth models using modified anisotropic eikonal equations[END_REF] transformed into the AEE and the Powell minimization algorithm. We have shown using a limited number of LGG cases that our approach could provide clinicians with two fundamental tumor dynamics characteristics: the spatial position of tumor seed point and its diffusivity ratio. Our optimization approach has been successfully evaluated on synthetic data where ground truth is available.

At this point, we have noticed the quasi-stability of the located source in both clinical cases, which is reassuring. Another key finding is related to the distance between the located source point and tumor barycenter. For Patient B, this distance ranges from 4.24mm to 5.47mm and for the 15 synthetically generated tumors, computed distance ranges from 1.42mm to over 13mm. The disparity between the computed source and its barycenter maybe explained by the spikiness index. Indeed, it is intuitive that "spherelike" tumors have a source which might be very close to its barycenter but "star-like" tumors or spiky tumors are prone to greater disparities. In the latter case, the barycenter may even end-up outside the tumor shape. This outcome contradicts the assumption that the tumor source is its barycenter as assumed in some studies as in [START_REF] Laigle-Donadey | Correlations between molecular profile and radiologic pattern in oligodendroglial tumors[END_REF] (tumor classification process based on the correlation between its located barycenter and molecular profile ), in [START_REF] Drabycz | An analysis of image texture, tumor location, and mgmt promoter methylation in glioblastoma using magnetic resonance imaging[END_REF] (MGMT gene methylation status prediction study in glioblastomas) and in [START_REF] Anderson | A hybrid mathematical model of solid tumour invasion: the importance of cell adhesion[END_REF][START_REF] Sottoriva | Cancer stem cell tumor model reveals invasive morphology and increased phenotypical heterogeneity[END_REF][START_REF] Kansal | Emergence of a subpopulation in a computational model of tumor growth[END_REF] (tumor growth simulation using cancerous cell concentration gradient computed from the center of tumor mass). Therefore, tumor barycenter-derived conclusions may not seem fully reliable and could cloud the true biological phenomena influencing tumor growth patterns. Furthermore, as noticed in Table 1, the source-to-barycenter distance ranging from 1 to 5mm in a tumor of 10mm-size demonstrates that choosing the barycenter as a source estimator leads to a large bias.

In Patient A, there is a jump of the source-barycenter distance value from 1.42mm to 5.38mm. This may be explained by a recent study from Bohman et al. [START_REF] Bohman | Magnetic resonance imaging characteristics of glioblastoma multiforme: implications for understanding glioma ontogeny[END_REF] aiming at a better understanding of glioma ontogeny (development from earlier stage) within the complexity of brain anatomy. In fact, Bohman et al. concluded that as the tumor grows, its boundary approaches the ventricles and the closer its "seed position " is to the ventricles, the larger the tumor would be. This conclusion agrees with our observation in Fig 7 where the tumor boundary becomes more distorted as it seems to be more "pulled" towards the left lateral ventricle. Note that while Bohman et al. use a basic geometric method to estimate the tumor barycenter, it mentions a discrepancy between barycenter and tumor source ranging between 1.25 mm and 15 mm. This is in agreement with our synthetic study and we hope that this work could provide a better understanding of tumor ontogeny.

According to our study, the tumor shape irregularities seems to be well correlated with the diffusivity ratio r. Setting the spikiness index leads to similar irregular tumor contours up to the geometry of their neighboring brain structures (presence of the ventricles for instance). As a matter of fact, the quasi-stability of this ratio for Patient A is justified by the quasistability of its shape over time. Its high values ranging from 40.5 to 42 agrees with its "V-like" irregularity (Fig 7). On the other hand, for Patient B, this ratio raises from 6.45 to 47.70 as its shape evolves from a "sphere-like" to a "star-like" (Fig 9). Interestingly, in a recent paper [START_REF] Stretton | Predicting the location of glioma recurrence after a resection surgery[END_REF], a different ratio (d w /ρ) was used to quantify the diffusivity of the tumor boundary where a low white matter diffusivity rate over proliferation ratio indicates a not very infiltrative glioma (bulky), whereas high ratio defines diffuse tumor. The main objective was to determine whether the tumor recurrence is a bulky or diffuse-type recurring tumor. We believe that a comparison between both of these ratios (d w /dg and d w /ρ) in future research works will be a valuable step in determining the patient-specific factors guiding tumor growth.

For patients C and D, which had DT-MRI anomalies as shown in (Fig 1 and2), our method proved its robustness as the estimated tumor source remained stable as time evolved. The symmetrization process induced a tremendous loss in white matter anisotropy as most of the fibers along which the tumor originally grew disappeared as demonstrated in (Fig 1) through FA histograms comparison. Nevertheless, the algorithm converged to a low diffusivity ratio value (r = 7.92). That was visually proven to be in con-cordance with the tumor "roundness" (Fig 11) and also remained constant throughout time, a fact in correlation with the very slow and slightly unchanged tumor boundary from t 1 to t 3 . The Euclidean distance between the estimated tumor source and its barycenter did not go over 2mm, which was quite expected when noticing the roundness of the tumor shape.

Patient D did not undergo the symmetrization process because of the large holes in the symmetric region to tumor segmented tissue (Fig 2-a) and the original diffusion fibers distorted by tumor growth were used. Despite this, the estimated source remained stable with 1mm spatial variation along z-axis. Its spatial disparity to the barycenter reaching a value of 3mm highlights that the increase of the tortuosity of the tumor quantified by the estimated diffusivity ratio implies an increase in the distance between the tumor barycenter and its original source (as in Patient B where it reached 5.47mm). Since a good estimation of the spikiness index is tied to the good quality of the acquired DT-MRI, the unexpected drop in r value from r t 1 = 20 to r t 2 = 16.2 can be explained by the major tumor-infiltration derived anomalies of the reconstructed diffusion tensors and also the interpolation of the small holes therein. Both cases with various DTI anomalies demonstrated a good performance when estimating a stable tumor source and also capturing the spikiness of its shape through the estimated diffusivity ratio. However, the prediction step was overlooked since it strongly requires a very good quality of the symmetric DTI to the affected invaded tissue.

Relying on key biological tissue characteristic (diffusivity in white and grey matters), our estimation of the tumor source location and diffusivity rate (or spikiness index) also showed a promising prediction of the further spatial propagation of the glioma boundary as time evolves. The good prediction in both LGG cases presented in Table 2 and Fig 7 provides additional confidence in the estimation of tumor source and diffusivity ratio. Note however that this spatial prediction was not based on a simulation from the tumor origin but from the tumor boundary at a previous time point using the estimated spikiness index. This is to avoid error accumulation leading to a weak prediction.

DT-MRI serves as the basis for several reaction-diffusion based models to extrapolate glioma invasion margins for radiation [START_REF] Konukoglu | Extrapolating glioma invasion margin in brain magnetic resonance images: Suggesting new irradiation margins[END_REF] and after resection [START_REF] Stretton | Predicting the location of glioma recurrence after a resection surgery[END_REF] using at least two acquired time points. Those two studies have been successively evaluated on one healthy subject and one LGG patient, and they confirm that the importance of using high resolution and high quality DT-MRI which is obviously a strong constraint. Including DT-MRI acquisitions in a standardized brain tumor imaging protocol would increase the number of recruited patients and could allow a thorough validation of these models.

Conclusion

To go beyond the classical tumor radial expansion model which does not take into account the tumor boundary irregularities ( [START_REF] Bohman | Magnetic resonance imaging characteristics of glioblastoma multiforme: implications for understanding glioma ontogeny[END_REF], [START_REF] Mandonnet | Continuous growth of mean tumor diameter in a subset of grade ii gliomas[END_REF], [START_REF] Harpold | The evolution of mathematical modeling of glioma proliferation and invasion[END_REF]), we have introduced a new method to estimate tumor growth parameter based on a traveling front propagation framework. The described method provides two main patient-specific spatial characteristics of a tumor using one clinical image at a single time point: its source location and an index quantifying the irregularity of its shape. The proposed macro-scale approach was evaluated using synthetic data and four patients with LGG. As a consequence of the estimation of both parameters, a promising prediction of tumor spatial behavior was pointed out. Some aspects of this work may be subject to improvements. For instance, a sensitivity analysis of the proposed method to the choice of tumor diffusion tensor in white matter is needed. Our proposed formulation leads to efficient computations, but a quantitative evaluation from experimental datasets should help refining the approach. Also, the white matter and brain mask segmentations could be improved by using state-of-the-art image classification tools on high resolution T1 MRI for instance.

Future steps include the application of the proposed method on highgrade gliomas and on a large database of glioma patients advancing a better validation of the predictive power previously highlighted. This could also lead to a statistical study of the correlation between the tumor source location, its grade and its diffusivity ratio. Thus, we could check whether lower spikiness index values are obtained for HGG since they appear to have less spiky and irregular shapes than LGG. Including the mass effect when addressing HGG growth would be a necessary step to include in the presented methodology. Hopefully, the coupling of mathematical models with clinical images should lead to a better understanding of tumor growth and to useful computational tools for diagnosis and therapy planning. First row: using the algorithm outcome at t 1 , we spatially determine the tumor evolution at t 2 . The blue contour represents the manually delineated tumor at the second time point and the red one is the extracted iso-time surface ŜisoT ime . Second row: using the algorithm outcome at t 2 , we predict the tumor evolution at t 3 . Tumor boundary at t 3 is colored in blue and the estimated one in red. Third row: using the algorithm outcome at t 3 , we predict the tumor evolution at t 4 . Same colors are used to show the similarity between the estimated and the ground-truth tumor contours.

x∈S

  Seg min d(x, ŜisoTime ) + x∈ ŜisoT ime min d(x, S Seg ) |S Seg | + | ŜisoTime | (9) where d denotes the closest Euclidean distance between a point and a surface.

Figure 1 :

 1 Figure 1: Data encountered problems (Patient C): fractional anisotropy (FA) asymmetry. (a): Non-symmetrized FA color map shows the dramatic diffusion signal distortion as the region where the tumor grew seems darker than the unaffected brain tissue. (b): Nonsymmetrized FA color map overlayed on glioma segmentation. (d) Symmetrized FA color map overlayed on glioma segmentation. A preliminary visual comparison between (b) and (d) clearly shows that when symmetrizing the diffusion tensor image and therefore the driven FA map, there is a tremendous loss in white matter anisotropy as most of the fibers along which the tumor originally grew disappear. The comparison of FA histograms only computed in the segmented tumor region (c) and (d) successively corresponding to (b) and (d) tumor figures demonstrates the anisotropy bias introduced by symmetrization process in this case.

Figure 2 :

 2 Figure 2: Data encountered problems (Patient D): large black holes in the DTI. (a) Tumor segmentation in red is overlayed on the fractional anisotropy map, a measure derived from Diffusion Tensor Imaging. The significantly missing diffusion signal in the contralateral region to tumor segmented tissue prohibited our proceeding to the symmetrization process and therefore we used for our estimation the affected tissue tensors by tumor infiltration. (b) The tumor boundary is drawn on the FLAIR image at the first and second time points in light and dark blue (in this slice they are identical) and at the third time point in red.

Figure 3 :

 3 Figure 3: Two axial slices of the first included patient (Patient A) where the manually segmented glioma (in white) is visualized on the diffusion tensor MR image. The symmetrization process in this case led to reliably recovered infiltration fiber tracts where the simulation of tumor growth starting from a seed point is legitimate.

Figure 4 :

 4 Figure 4: Simulation from a source point of the iso-time contours representing the tumor invasion process with a diffusivity ratio d w /d g equal to 25. The synthetic tumor displayed in red is created by thresholding the time distance map at T = T Obs . (A) represents the axial slice, (B) coronal slice and (C) sagittal slice.

Figure 5 :

 5 Figure 5: The minimization surface after fixing r and S z . The x axis of this 3D plot represents x coordinate of tumor source and the y axis represents its y coordinate.

Figure 6 :

 6 Figure 6: Synthetic tumors generated at different locations with different diffusivity ratios and observation times. First row: At a fixed source location, two iso-time contours with the same observation time (T Obs = 900days) are respectively colored in red for a diffusivity ratio (r = 1) and in light blue (r = 50). The darker blue contour represents the tumor further spatial invasion with (r = 50) and (T f = 1200days). Second row: 3 synthetic tumors are generated using the same observation time value T Obs and the same diffusivity ratio (r = 25) for 3 different source locations. Third row: At the same fixed location, five iso-time contours with the same observation time are colored in red (r = 1), in blue (r = 10), in green (r = 50) and in yellow (r = 100).

Figure 7 :

 7 Figure 7: Different axial slices of Patient A showing the discrepancy between the real boundary and the estimated one at 3 different successive time points.First row: using the algorithm outcome at t 1 , we spatially determine the tumor evolution at t 2 . The blue contour represents the manually delineated tumor at the second time point and the red one is the extracted iso-time surface ŜisoT ime . Second row: using the algorithm outcome at t 2 , we predict the tumor evolution at t 3 . Tumor boundary at t 3 is colored in blue and the estimated one in red. Third row: using the algorithm outcome at t 3 , we predict the tumor evolution at t 4 . Same colors are used to show the similarity between the estimated and the ground-truth tumor contours.

Figure 8 :

 8 Figure 8: Evaluation curves plotted using the distance error measure C isoT ime (y%) for y ∈ (50, 90, 100). For each patient (A or B), we try the predict the evolution at a successive time point (A/B, t i ) based on the estimated parameters at the previous one.

Figure 9 :

 9 Figure 9: Four different axial slices of Patient B showing the remarkable tumor boundary spikiness evolution from the initial time point t 1 to the final one t 4 . Tumor contours at successive time points are respectively colored in: dark blue for t 1 , light blue for t 2 , dark red for t 3 and light red for t 4 .

Figure 10 :

 10 Figure 10: Four different axial slices of Patient B illustrating the prediction of the tumor evolution based on two different time points. Using the algorithm outcome at t 2 , the tumor front is predicted at the next time point t 3 . The red contour represent the predicted tumor boundary and the blue one represents the ground truth or the real tumor boundary.

Figure 11 :

 11 Figure 11: Patient C: three different axial slices of the manually segmented tumor boundary at t 1 in dark blue, at t2 in dark blue and at t 3 in red.

Table 1

 1 represents the outcome of the estimation method for the four time points (t 1 , t 2 , t 3 , t 4 ) of Patient A.

		t 1	t 2	t 3	t 4
	S(x)	(120, 111, 31) (120, 111, 31) (122, 115, 30) (121, 116, 30)
	d(G, S) in mm	2.00	1.73	1.42	5.38
	r	41.98	40.47	40.58	41.38

Table 1 :

 1 Outcome

of the minimization algorithm of Patient A which was separately computed at four time points. The distance d denotes the Euclidean distance between the tumor barycenter G and its estimated source S.

Table 2 :

 2 symDist in mm C overlap in % (P atientA, t 2 ) Prediction evaluation based on the error measures C symDist and C overlap computed for Patient A and B to predict the evolution of the tumor at specific time points t i using the minimization outcome of the previous time t i-1 .

		1.97	85.63
	(P atientA, t 3 )	1.79	87.76
	(P atientA, t 4 )	1.74	88.04
	(P atientB, t 2 )	1.53	82.20
	(P atientB, t 3 )	1.48	84.79
	(P atientB, t 4 )	1.13	88.02

Table 3 :

 3 -Outcome of the minimization algorithm of Patient B separately computed at four time points.

		t 1	t 2	t 3	t 4
	S(x)	(34, 145, 117) (35, 141, 117) (34, 141, 118) (34, 140, 119)
	d(G, S) in mm	5.09	5.47	5.09	4.24
	r	6.45	37.15	30.12	47.70

This is not strictly true if one uses the speed term F taking into account the transient speed as in Equation 4 or 5. However, the difference in simulations due to the absolute value of the diffusivities was found to be negligeable.

Acknowledgments

We thank Ender Konukoglu for his assistance in understanding his research work previously referenced. This work was partially funded by the ERC advanced Grant MedYMA.


Jalal M Fadili

Jean-Luc Starck

M J Fadili

J.-L Starck

MONOTONE OPERATOR SPLITTING FOR OPTIMIZATION PROBLEMS IN SPARSE RECOVERY

Keywords: Convex analysis, Non-smooth optimization, Monotone operator splitting, Sparse recovery

HAL is

INTRODUCTION

The complex structures of natural signals and images require tools in order to make use of their intricate redundancies. To capture this complexity, we have witnessed a flurry of research activities where researchers spanning a wide range of viewpoints have advocated the use of sparsity and overcomplete signal/image representations. It has turned out that exploiting sparsity and overcompleteness offers striking benefits in a wide range of signal and image processing areas. The development of sparse representations in inverse problems is facing major challenges both on theoretical and practical sides. These problems include the design of fast optimization algorithms to handle real world data and application to many ill-posed problems. This work focuses on several optimization problems involved in linear inverse problems where the solution s = Φx, Φ : H → K, is assumed to be sparsely represented in an overcomplete dictionary Φ of waveforms (ϕ i ) i∈I . Φ is a generally a frame of H. Our goal is the generic minimization of functions of the form

f = f 1 + f 2 (P) min x∈H f 1 (x) + f 2 (x) , (1)
where f 1 , f 2 are closed convex functions which are not necessarily differentiable, are not infinite everywhere, and their domains have non-empty intersection. f 1 is typically a convex sparsity-promoting penalty and f 2 is a function measuring the consistency to data. Problem (P) covers the case of the ℓ 1 decoder known also as Basis Pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candès | Near optimal signal recovery from random projections: Universal encoding strategies?[END_REF] when f 1 (x) = x 1 , and f 2 (x) is the indicator function of the affine subspace {x ∈ H s = Fx}, F = HΦ is a linear operator which is typically a composition of a linear measurement or degradation operator H and the dictionary Φ. When the observation s is contaminated by noise, the equality constraint must be relaxed to a noise-aware variant. Problem (P) becomes typically Basis Pursuit DeNoising (BPDN) in its augmented Lagrangian form when f 1 (x) = λ x 1 and

f 2 (x) = 1 2 s -Fx 2 [1]
, and is also known as Lasso [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF].

This augmented Lagrangian form has received a considerable attention, see e.g. [START_REF] Figueiredo | An em algorithm for wavelet-based image restoration[END_REF][START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraints[END_REF][START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF]. Only a few algorithms have been proposed to solve the more difficult to solve constrained form where f 2 is the indicator function of the closed convex set {x ∈ H s -Fx σ}; e.g. [START_REF] Candès | Stable signal recovery from incomplete and inaccurate measurements[END_REF][START_REF] Van Den | Probing the pareto frontier for basis pursuit solutions[END_REF][START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF]. The Dantzig selector [START_REF] Candès | The dantzig selector: Statistical estimation when p is much larger than n[END_REF] is also a special instance of (P)

when f 2 is the indicator function of {x ∈ H F * (s -Fx) ∞ τ }.
In this paper, we formalize all these optimization problems within a unified framework of convex optimization theory, and invoke tools from convex analysis (e.g. duality, proximity operators) and maximal monotone operator splitting. We characterize all these optimization problems, and to solve them, we propose fast iterative and convergent algorithms combining forward-backward and/or Peaceman/Douglas-Rachford splitting iterations. With nondifferentiable sparsity-promoting penalties, the proposed algorithms amount to iterative shrinkage. Furthermore, the computational burden of these algorithms is essentially invested in one application at each iteration of fast implicit operators associated to the involved dictionary Φ and the linear measurement operator H and their respective adjoints. This makes them very competitive for large-scale problems such as those rising in image processing.

NOTATION AND TERMINOLOGY

Let H a Hilbert space (here a real vector space) equipped with the inner product ., . and associated norm . . Let I be the identity operator on H. The operator (induced) norm of A is denoted A = sup x∈H Ax x . Let . p , p ≥ 1 be the ℓ p -norm with the usual adaptation for the case p = +∞. Denote B σ p the closed ℓ p -ball of radius σ > 0. Throughout the paper, we denote by ⇀ weak convergence and → strong convergence in H.

A real-valued function f : H → (-∞, +∞] is coercive, if lim x →+∞ f (x) = +∞. The domain of f is defined by dom f = {x ∈ H : f (x) < +∞} and f is proper if dom f = ∅. We say that a real-valued function f is lower semi-continuous (lsc) if lim inf x→x0 f (x) ≥ f (x 0 ). Γ 0 (H) is the class of all proper lsc convex functions from H to (-∞, +∞]. We note ı C the indicator of the convex set C: ı C (x) = 0 if x ∈ C
, and +∞ otherwise. The interior of C will be denoted int C and its relative interior ri C; see e.g. [START_REF] Rockafellar | Convex analysis[END_REF]. We recall that F is the frame synthesis operator associated to a frame of H if its adjoint F * (analysis operator) satisfies the generalized Parseval relation with bounds c 1 and c

2 c 1 x 2 F * x 2 c 2 x 2 , 0 < c 1 c 2 < +∞ .
The frame is tight when c 1 = c 2 = c and FF * = cI.

PROXIMAL CALCULUS

We first define the notion of a proximity operator, which was introduced in [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF] as a generalization of convex projection operator.

Definition 1 (Proximity operator [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF]) Let f ∈ Γ 0 (H).

Then, for every x ∈ H, the function z

→ 1 2 x -z 2 + f (z)
achieves its infimum at a unique point denoted by prox f x.

The uniquely-valued operator prox f = (I + ∂f ) -1 : H → H thus defined is the proximity operator of f . It will be convenient to introduce the reflection operator rprox f = 2 prox f -I.

The following lemma provides an important result on the proximity operator of the pre-composition of f ∈ Γ 0 (H) with a bounded affine operator A := F • -y. It will be at the heart of many algorithms in the rest of the paper.

Lemma 2 Let F be a bounded linear operator such that int(dom(f ) ∩ Im(A)) = ∅. Then f • A ∈ Γ 0 (H) and (i) F is a tight frame. Then prox f •A (x) = x + c -1 F * prox cf -I (Fx -y) . (2) (ii) F is a general frame. Let µ t ∈ (0, 2/c 2 ). Define u (t+1) =µt " I -prox µ -1 t f " • " µ -1 t u (t) + Ap (t) " , (3)
p (t+1) =x -F * u (t+1) . ( 4
)
Then u (t) → ū and p

(t) → prox f •A = x -F * ū.
More precisely, both u (t) and p (t) converge linearly and the best convergence rate O c2-c1 c2+c1 t is attained for

µ t ≡ 2/(c 1 + c 2 ). (iii) If c 1 = 0, apply (3) with µ t ∈ (0, 2/c 2 ). Then u (t) ⇀ ū and p (t) → prox f •A = x -F * ū at the rate O(1/t).
The proof is based on Fenchel-Rockafellar duality [START_REF] Rockafellar | Convex analysis[END_REF] and is detailed in a forthcoming long version of the paper. In the special case of tight frames, we recover the same result as in [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF]. Note that f • A ∈ Γ 0 (H) is always verified when F corresponds to a frame. The recursion ( 3) is basically a forward-backward splitting [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] applied to the dual problem. Other solvers can be used instead such as Nesterov's scheme [START_REF] Nesterov | Gradient methods for minimizing composite objective function[END_REF].

SPARSE RECOVERY INVERSE PROBLEMS

Splitting schemes

Recall that our goal is the minimization of f = f 1 + f 2 in (1) for f 1 , f 2 ∈ Γ 0 (H) such that their domain has a nonempty intersection. Let T the set of solutions of problem (P), T = {x ∈ H x ∈ (∂f ) -1 (0)}. When T = ∅, the proximaltype recursion constructed as x (t+1) = prox βf (x (t) ), β > 0 is a fundamental algorithm for solving (1), and x (t) ⇀ x ∈ T . The main difficulty with this iteration is that prox βf may be hard to compute in general, depending on the function f . This is for instance the case in most inverse problems arising in image and signal processing. Splitting methods for problem (P) are algorithms that do not attempt to evaluate the proximity mapping prox βf of the combined function f but instead perform a sequence of calculations involving separately the individual proximity (resolvent) operators prox βf1 and prox βf2 .

The Douglas/Peaceman-Rachford (DR) family is the most general pre-existing class of monotone operator splitting methods. Given a fixed scalar β > 0 and a sequence α t ∈ (0, 2), this class of methods may be expressed via the recursion

x (t+1) = 1 - α t 2 I + α t 2 rprox βf1 • rprox βf2 x (t) . ( 5
)
The DR splitting has been brought to light in the image and signal processing literature very recently by [START_REF] Combettes | A Douglas-Rachford splittting approach to nonsmooth convex variational signal recovery[END_REF]. Its most general form is due to [START_REF] Combettes | Solving monotone inclusions via compositions of nonexpansive averaged operators[END_REF] where its convergence is established.

Theorem 3 ([16]

) Suppose that T = ∅. Let β ∈ (0, +∞), (α t ) t∈N be a sequence in (0, 2) such that t∈N α t (2α t ) = +∞, and let (a t ) t∈N and (b t ) t∈N be sequences in H such that t∈N α t ( a t + b t ) < +∞. Fix x (0) ∈ H and define the recursion,

x (t+1/2) = prox βf 2 " x (t) " + bt , (6)
x (t+1) =x (t) + αt " prox βf 1 " 2x (t+1/2) -x (t) " + at -x (t+1/2) " . (7)
Then x (t) ⇀ x and prox β∂f2 (x) ∈ T .

The sequences a t and b t in ( 6)-( 7) play a prominent role as they formally establish the robustness of the DR algorithm to numerical errors when computing the proximity operators of f 1 and f 2 (e.g. via a nested inner iteration).

Sparse recovery with equality-constraint

From now, we assume that H is finite-dimensional. Let us first consider the general equality-constrained problem

(P eq ) : min x∈H Ψ(x) s.t. s = HΦx , (8)
Ψ ∈ Γ 0 (H) is a sparsity-promoting potential function, H is a bounded linear operator, and the equality constraint corresponds to an affine subspace that we denote C eq = {x ∈ H s = HΦx}. We assume the following

Assumption 4 (i) We suppose in the sequel that Ψ is addi- tive; that is Ψ(x) = i∈I ψ i (x i ). (ii) Ψ is coercive. (iii) F = HΦ is surjective. (iv) s ∈ F (ri dom Ψ).
The following proposition gives the algorithm to solve (P eq ) and establishes its convergence. Proposition 5 Suppose that Assumptions 4(i)-(iv) are satisfied. Let β, (µ t ) t∈N , (a t ) t∈N and (b t ) t∈N satisfy the conditions in Theorem 3. Fix x (0) ∈ H and apply the DR recursion (6)- [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] with

prox βf2 (x) =P Ceq (x) = x + F * (FF * ) -1 (s -Fx) , ( 9
)
prox βf1 (x) = prox βΨ (x) = prox βψi (x i ) 1 i |I| . ( 10
)
Then x (t) → x and P Ceq (x) is a solution to (P eq ).

To implement this DR splitting algorithm, we need to express prox βψi . Its closed-form expression is given in [START_REF] Fadili | Inpainting and zooming using sparse representations[END_REF] for a large class of penalty functions ψ i . For instance, with

ψ i (x i ) = λ|x i |, prox βψi = ST λ (x)
the soft-thresholding operator.

Sparse recovery with inequality constraint

We then consider the inequality-constrained problem

(P σ ) : min x∈H Ψ(x) s.t. s -HΦx σ , (11)
which can be viewed as a noise-aware variant of (P eq ). (P σ ) only asks the reconstruction be consistent with the data such that the reconstruction error energy is within the noise level σ. Problem (P σ ) is challenging to solve for large-scale problems. The following proposition is a consequence of Lemma 2 and Theorem 3.

Proposition 6 Under Assumptions 4(i)-(iii) and ∈ int (s + B σ

2 )-F (ri dom Ψ). Define β, (µ t ) t∈N , (a t ) t∈N and (b t ) t∈N as before. Fix x (0) ∈ H and apply the DR recursion (6)-( 7) with prox βf1 as in [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF] and

prox βf2 (x) = P Cσ (x) , (12)
where P Cσ (x) is given by Lemma 2 with y = s, prox f (u) =

P B σ 2 (u) = (I-prox σ • )(u) and prox σ • (u) = u 1 -σ u + .
Then x (t) → x and P Cσ (x) is a solution to (P σ ).

If the data fidelity constraint in (P σ ) is replaced by any ℓ pnorm for p ≥ 1, the algorithm described above can be extended easily by appropriate substitution of P B σ 2 (u) by the projector onto the ℓ p -ball B σ p in Proposition 6. For example, in the case of uniform quantization noise, a good choice would be p = ∞, implying P B σ ∞ (u) = ui max(|ui|/σ,1) i . Another interesting variant is when the noise is impulsive where the value p = 1 is advocated, and the projector onto P B σ 1 can be computed through soft-thresholding.

The Dantzig Selector

Let's now turn to the problem

(P DS ) : min x∈H Ψ(x) s.t. F * (s -Fx) ∞ τ . ( 13
)
To avoid the unique trivial solution x = 0, we assume that τ < F * s ∞ . The compact convex constraint set in ( 13) is C DS . The Dantzig Selector [START_REF] Candès | The dantzig selector: Statistical estimation when p is much larger than n[END_REF] is when Ψ is the ℓ 1 -norm. In this case, (P DS ) can be recast as a linear program [START_REF] Candès | The dantzig selector: Statistical estimation when p is much larger than n[END_REF]. Let G = F * F be the (self-adjoint) Gram operator and = F * s. Beside Assumption 4(i) and (iii), suppose that 0 ∈ int (z + B τ ∞ ) -G (ri dom Ψ). Again, by straightforward application of Lemma 2 and Theorem 3 we arrive at Proposition 7 Fix x (0) and apply the DR recursion (6)- [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] with prox βf1 as in [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF] and

prox βf2 (x) =P CDS (x) , ( 14
)
where P CDS (x) is given by Lemma 2(iii) with G in lieu of F, y = z and (Iprox f )(u) = (I -P B τ ∞ )(u) = ST τ (u). Then x (t) → x and P CDS (x) is a solution to (P DS ).

Analysis prior

In the problems considered above, we seek a sparse set of coefficients x and the solution image is synthesized from these representation coefficients s = Φx. Such a prior is called a synthesis prior. A less common approach seeks an image s whose coefficients x = Φ * s are sparse. Such a prior is called an analysis prior. This is typically the case for the Total Variation regularization. (P eq ), (P σ ), its Lagrangian form or (P DS ) can be modified to incorporate such a prior. For example, (P σ ) becomes

(Q σ ) : min s∈K Ψ(Φ * s) s.t. s -Hx 2 σ , (15)
and similarly for the other problems. Despite its apparent difficulty, it turns out that the analysis-prior problems such as (Q σ ) can still be solved using ( 6)-( 7). The only modification is the computation of the proximity operator of f 1 = Ψ • Φ * . The expression of the latter is easily accessible owing to Lemma 2(iii). Indeed, we have the following corollary.

Corollary 8 Let µ t ∈ (0, 2/ Φ 2 ). Define p (t+1) = s -µtΦ " I -prox µ -1 t Ψ " • " µ -1 t u (t) + Φ * (s -Φu (t) ) " .
Then p (t) → prox Ψ•Φ * at the rate O(1/t).

NUMERICAL EXAMPLES

The proposed framework has been applied successfully to several inverse problems such as denoising, component separation, inpainting and CS. For obvious space limitation, we only report experiments on CS reconstruction.

In the first one, we solve (P eq ) with Ψ(x) = x 1 , H is m × 4m random sensing matrix generated from the Hadamard ensemble, Φ = I (canonical basis), the support of x was selected uniformly at random with sparsity level 0.05m and normally distributed non-zero entries. Four algorithms were compared: our DR scheme ( 9)- [START_REF] Donoho | Stable recovery of sparse overcomplete representations in the presence of noise[END_REF], the LARS/LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the lasso[END_REF], linear-programming with interior point solver [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF][START_REF] Candès | Near optimal signal recovery from random projections: Universal encoding strategies?[END_REF] and StOMP [START_REF] Donoho | Sparse solution of underdetermined linear equations by stagewise orthogonal matching pursuit[END_REF]. The results are shown in Table 1 (similar results were observed for Fourier measurements). As claimed, our DR splitting solver is among the fastest, and becomes faster than the other algorithms as m gets large. Fig 1 depicts the image recovered using ( 10)-( 12) from 17% noisy (SNR=30dB) Fourier measurements of the 256 2 Mondrian image. The dictionary Φ was the wavelet transform.

CONCLUSION

In this paper, we proposed a general framework for solving several optimization problems arising in sparse solutions of linear inverse problems. The algorithms were characterized and proven convergent. Sparse recovery from CS measurements illustrated the applicability of this framework.

Fig. 1 .

 1 Fig. 1. Left: original 256 2 Mondrian image. Right: recovery from 17% noisy (SNR=30dB) Fourier measurements with the wavelet sparsity basis, SNR=22dB.

Table 1 .

 1 Computation times in seconds for exact CS reconstruction from Hadamard measurements.

	m	2 8	2 9	2 10	2 11	2 12	2 13	2 14
	DR	0.58	1.17	2.52	3.80	6.43	15.18	36.20
	LARS	0.28	1.00	3.65	6.60	10.54	24.77	59.38
	LP-IP	2.07	9.85	45.35	163.62	631.57	3191.38	13175.63
	StOMP	0.34	0.84	2.36	3.61	6.38	16.65	43.13