Djawida Dib
email: djawida.dib@inria.fr

Nikos Parlavantzas
email: nikos.parlavantzas@inria.fr

Christine Morin
email: christine.morin@inria.fr

Towards Multi-Level Adaptation for Distributed Operating Systems and Applications

Keywords: distributed operating system, service-oriented architecture, multi-level adaptation, distributed shared memory

Distributed operating systems simplify building and executing applications on large-scale infrastructures, such as clusters, grids and clouds. These systems operate in a constantly changing environment characterized by varying application needs and varying physical infrastructure capabilities. To handle the diversity and dynamism of both the applications and the underlying infrastructures, the distributed Operating System (OS) should continually adapt to its changing environment. Two challenges arise in this context: how to design the distributed OS in order to facilitate dynamic adaptation, and how to ensure that OS-level adaptation does not conflict with application-level adaptation. This paper proposes to address these challenges by: (1) building the distributed OS as an assembly of adaptable services following the service-oriented architecture; and (2) using a common multi-level adaptation framework to adapt both the OS and the application layers in a coordinated way. Moreover, the paper presents experimental evidence of the usefulness of this approach in adapting the distributed shared memory service of a specific distributed OS.

Introduction

Distributed Operating Systems (OSs) provide common services and abstractions for building applications on large-scale infrastructures such as clusters [START_REF] Morin | Towards an efficient single system image cluster operating system[END_REF][3], grids [START_REF] Morin | Xtreemos: A grid operating system making your computer ready for participating in virtual organizations[END_REF] [START_REF] Grimshaw | The legion vision of a worldwide virtual computer[END_REF] and clouds [START_REF]Openstack[END_REF] [START_REF] Opennebula | [END_REF]. These systems operate in a constantly changing environment characterized by varying application needs and physical infrastructure capabilities. To deal with this dynamicity, it is essential for the distributed OS to support dynamic adaptation. For example, it should support dynamically modifying communication protocols depending on application usage patterns or network conditions.

Service-Oriented Architectures (SOAs) have recently emerged as a popular approach for building flexible software systems [START_REF] Papazoglou | Service oriented architectures: approaches, technologies and research issues[END_REF]. A service-oriented system is an assembly of services, where each service represents a well-defined function.

The research leading to these results is co-funded by the Brittany Region.

The flexibility of SOA results from the dynamic capability to publish, discover and use services as well as modify their implementations without impacting their consumers. While SOA has been widely applied to build distributed applications [START_REF] Tosi | Towards autonomic service-oriented applications[END_REF], it is less explored in building the supporting system software. In this paper, we propose to apply the SOA at the OS level in order to facilitate the OS dynamic adaptation.

Supporting dynamic adaptation of the distributed OS is introducing further challenges, such as detecting when adaptation is needed and deciding what and how to adapt in order to achieve specific objectives. One important challenge, which is not addressed by related research work, is ensuring that OS-level adaptation does not conflict with application-level adaptation. For example, a particular performance problem could be addressed by replacing a service at the application layer or by allocating more resources at the OS layer; performing both adaptations could be inefficient or damaging. To address this challenge, we propose to use a multi-level adaptation framework that coordinates the adaptation of both the application and the OS layers based on configurable policies. Experimental results are provided to demonstrate the usefulness of this approach.

This paper is organized as follows. Section 2 introduces the service-oriented OS. Section 3 discusses the need for a multi-level adaptation framework. Section 4 presents the multi-level adaptation architecture. Section 5 shows the usefulness of the proposed approach with the distributed shared memory example. Section 6 presents some related work. Finally, Section 7 concludes this paper and discusses future work.

Distributed Service-Oriented Operating Systems

Current distributed operating systems provide for distributed infrastructure users what a traditional OS provides for single computer users [START_REF] Barak | The MOSIX Distributed Operating System: Load Balancing for UNIX[END_REF][14] [START_REF] Morin | Xtreemos: A grid operating system making your computer ready for participating in virtual organizations[END_REF]. Specifically, a distributed OS extends in a transparent way traditional OS functionalities. For example, the memory management functionality in a traditional OS is extended to form the Distributed Shared Memory (DSM) functionality in a distributed OS. Previous research works have proposed, compared and classified several algorithms to implement each OS functionality (e.g., DSM [START_REF] Tam | A taxonomy-based comparison of several distributed shared memory systems[END_REF] and resource management [START_REF] Krauter | A taxonomy and survey of grid resource management systems for distributed computing[END_REF]). Depending on the system environment, usually only one algorithm provides optimal performance. For example, a resource management algorithm can be very efficient to manage a small number of resources but it could be less efficient than other algorithms to manage a lot of resources.

Furthermore, distributed OSs operate in a constantly-changing environment, characterized by two levels of dynamicity: the dynamicity of the underlying infrastructure (resources may connect and disconnect at any time leading to variable resource availability), and the dynamicity of the application requirements (such as resource requirements). To accommodate this dynamicity while providing the best quality of service to applications, we believe that the distributed OS should be dynamically adaptable. Specifically, we believe that the OS should support selecting the best algorithm for each functionality at runtime. The se-!" #$%&'(')&*+$,-* .,*+/0(1$23,/ 4),1)&35$'3), #)&&)6('7*(1%, 8*9&3-$'3), "3/+$'3), 47*-:9)3,'3,/ ;*$:( -),131'*,-0 !'+),/( -),131'*,-0 Fig. 1. Service-oriented operating system overview lection could be done either according to the state of the whole system or in a more customized way, for each kind of application. In the latter case, two or more algorithms for the same functionality should be able to operate simultaneously, each one supporting the applications for which it was selected. For example, the fault tolerance functionality may use checkpointing for some applications while using replication for others.

To build such a dynamically-adaptable OS we propose to use SOA, thus enabling its well-known benefits in terms of interoperability and dynamic recomposition to be extended to the OS layer. Then, the OS is built as a composition of services, each service representing a distributed OS functionality. Figure 1 illustrates this idea with examples of some existing algorithms for three distributed OS services. The DSM service may use the weak or the strong consistency models, the fault tolerance service may use replication, migration or checkpointing, and the energy saving service may use a consolidation protocol [START_REF] Ferreto | Server consolidation with migration control for virtualized data centers[END_REF] to minimize the number of used resources or use the "follow the sun" protocol [START_REF] Treinen | Following the sun: case studies in global software development[END_REF] to utilize mostly sun energy.

The main challenge in building such a dynamically-adaptable OS is to determine when and how to change algorithms while taking into account the overhead of the change. Indeed, selecting the appropriate algorithm for each service is a complex task and involves thoroughly evaluating the algorithms as well as the parameters influencing their behaviors. To demonstrate this point, we investigate in Sect.5 two algorithms for the DSM service as well as the parameters influencing their behaviors.

Multi-Level Adaptation Framework

Apart from the distributed OS, it is frequently necessary to dynamically adapt the distributed applications themselves. This may be motivated by changes in user requirements, such as demands for higher precision in calculation results; or in the available computing resources, such as node failures or new nodes joining the system. The application adaptation is required in such cases to enable applications to work correctly despite all possible environment degradations, and to take advantage of new opportunities that can occur. Indeed, modern distributed applications are increasingly built according to SOA in which applications are composed of multiple services and are prepared to handle services appearing and disappearing at any time.

Traditionally, the adaptation of each layer of a distributed system is managed separately based on layer-specific knowledge and objectives. For example, there are research works focusing exclusively on OS layer adaptation [START_REF] Seltzer | Self-monitoring and self-adapting operating systems[END_REF] [START_REF] Teller | Insights into providing dynamic adaptation of operating system policies[END_REF] and other works focusing exclusively on application layer adaptation [START_REF] Tosi | Towards autonomic service-oriented applications[END_REF]. However, the adaptation of one layer clearly impacts the others. As a result, adapting each layer separately can create conflicts, unpredictable results, or redundancies. For example, adapting the OS to accommodate reduced network bandwidth may be unnecessary if the application reacts to the same situation by replacing the remote service with a locally available one. To avoid such interference between adaptation actions, we propose to use a multi-level adaptation framework (such as SAFDIS [START_REF] Gauvrit | SAFDIS: A Framework to Bring Self-Adaptability to Service-Based Distributed Applications[END_REF]), capable of interacting and dynamically adapting both layers in a coordinated way. The framework must maintain knowledge about the whole system and perform adaptation actions according to application-specific and system-wide goals.

Architecture

This section presents a generic architecture to apply the concept of multi-level adaptation. This architecture is composed of two service-oriented layers, application and OS layers, which are on top of a physical infrastructure layer (see Fig. 2). The considered physical infrastructure layer is composed of the hardware and a traditional OS installed on each node. All layers expose rich monitoring mechanisms and are controlled by the multi-level adaptation framework. Service-oriented layers expose also actuation mechanisms enabling the framework to trigger adaptation actions. The adaptation process of the framework follows the four main phases of the MAPE model [START_REF] Kephart | The vision of autonomic computing[END_REF]: Monitoring, Analysis, Planning and Execution. A description of the phases is given below. 1. Monitoring: the framework collaborates with the sensors of each layer in order to provide a dynamic view of the whole system. At the physical infrastructure layer, the sensors inform the framework about the availability and the workload of resources. At the OS layer, the sensors inform the framework about all available algorithms for each OS service, either statically at the start-up of the system or dynamically when a new algorithm component is added to the system. At the application layer, the sensors inform the framework about all available services in the repository as well as applications requirements. The mechanism used to monitor data may be either push or pull, depending on the nature of data and its frequency changes. The gathered information is stored in a knowledge base to be used later in the other phases. 2. Analysis: the framework analyzes the monitored data and uses policies to decide whether an adaptation is required or not. For instance, when the monitored data shows a deterioration of the system performance, the analysis policy determines if the OS, the application, or both should be adapted. The analysis policy triggers next phases only if the estimated overhead of the adaptation is less important than its benefit. 3. Planning: the framework defines a plan to apply the decided adaptation at the analysis phase. The plan consists of set of actions which are scheduled properly. For example, if both application and OS layers should be adapted then the planning strategy will specify which layer will be adapted first. 4. Execution: the framework collaborates with actuators of the service-oriented layers to perform the plan defined at the planning phase. At the OS layer, actuators enable either to switch between algorithms or to ensure that two (or more) algorithms coexist. At the application layer, actuators may modify application parameters or modify the application structure, such as replacing services or changing service interconnections.

Illustrative Examples

In this section we investigate two algorithms of the DSM functionality to illustrate the usefulness of the proposed architecture. The objective of this section is not to improve or to optimize existing DSM algorithms, but to demonstrate that each algorithm can be better than others in particular cases. The DSM service makes the main memory of a cluster of computers look like a single memory, using replicated data consistency models. We have evaluated two algorithms of the sequential consistency model: Broadcast-on-write and Invalidate-on-write. Both algorithms guarantee that any read of a data copy returns the last value assigned to this data. This is done by invalidating all remote data copies when a process access in write mode its local copy. The main difference takes place once the write operation finishes. In the broadcast-on-write algorithm, the owner of the modified copy updates automatically all the other copies. Whereas in the invalidate-on-write algorithm, the invalidated copies will be updated on demand when a page fault exception occurs. If the OS uses the broadcast algorithm, no fault page exception will occur. Nevertheless, some invalidated copies may be unnecessary updated if these copies are never subsequently used by their owners. Therefore, it is difficult to predict which algorithm will perform better. In this perspective we measure, in the next subsection, the execution time of a synthetic application with each algorithm while varying some parameters likely to impact the system performance.

Experiment Environment

To carry out these experiments, we have used and extended the Kerrighed OS [START_REF] Morin | Towards an efficient single system image cluster operating system[END_REF]. The existing implementation of Kerrighed uses the invalidation algorithm for the DSM functionality. We have implemented the boadcast algorithm and an actuator to switch between the two algorithms at runtime. In addition to the Kerrighed modification, we have developed a synthetic distributed application that uses shared memory objects for its inter-process communications. This application consists of a configurable number of processes that share data objects, each process run on a different node. One process creates a configurable number of shared data objects with different sizes and all other processes can access them either in read or in write mode. The application access pattern follows the notion of Write-Run defined in [START_REF] Eggers | A characterization of sharing in parallel programs and its application to coherency protocol evaluation[END_REF]. A write-run is a sequence of write accesses to a shared object by a single process uninterrupted by any access of the other processes. The number of the sequential writes, called the write-run length, represents in these experiments the access pattern parameter. Attention has been paid to these parameters because of their dynamicity during the execution of a distributed shared memory application. The experiments were performed on the paradent cluster of the Grid'5000 site of Rennes that consists of Carri System CS-5393B nodes supplied with 2 Intel Xeon L5420 processors (each with 4 cores at 2.5 GHz), 32 GB of memory, and Gigabit Ethernet network interfaces. The results are discussed in the next subsection.

Results

Figure 4(a) shows the impact of the number of data object copies on the application execution time. The application parameters were set at one data object of 10 bytes size with a write-run length equal to one, only the number of data object copies was adjusted. We notice that the performance with both algorithms is equal when there are two data object copies. The invalidation algorithm performs better when there is only one copy, but the broadcast algorithm performs increasingly better when there are at least three copies.

Figure 4(b) shows the impact of the write-run length on the application execution time. The application parameters were set at one data object of 10 bytes size with 4 data object copies, only the write-run length was adjusted. The application execution times with both algorithms converge when the write-run length is set at 13. If the write-run length is smaller then the convergence point, the broadcast algorithm performs better and conversely if it is greater. $"%%&% $'()*+,-.%/0)1*2034,5"%%%6% %%%%%%7% $'()*+,8.%/0)1*2034,9:"%6% %%%%%%;% $'()*+,9.%/0)1*2034,5"%%%6% %%%%%%<% $'()*+,8.%/0)1*2034,5="> ?"%%&% $'()*+,-.%/0)1*2034,5"%6 %%%%%%7% $'()*+,-.%/0)1*2034,9"%6 %%%%%%;% $'()*+,-.%/0)1*2034,="%6 %%%%%%<% $'()*+,8.%/0)1*2034,9"> Fig. 3. Impact of (a) the number of data object copies, (b) the write-run length, (c) the data object size, and (d) the number of data objects on the execution time of the application.

Figure 4(c), compared to Fig. 4(a) and Fig. 4(b) shows the impact of the data object size on the application execution time. The application parameters were set at one data object of 10 3 bytes size. The number of data object copies and the write-run length were adjusted to define four configurations (A,B,C,D). Each configuration was previously tested with an object size set at 10 bytes (see Fig. 4(a) and Fig. 4(b)). This comparison shows that the object size does not impact the application execution time in our experimental environment.

Figure 4(d) shows the impact of the number of data objects on the application execution time. The application parameters were set at two data objects of 10 bytes size. The number of data object copies and the write-run length was adjusted to define four configurations. The first configuration (A), compared to its analogous one with one shared object in Fig. 4(a), shows an increase of the application execution time with both algorithms. In configurations (B) and (C), we have increased the write-run length and there we notice that the invalidation algorithm, respectively, catches up and overtakes the broadcast algorithm in terms of performance. In the configuration (D), we have decreased the number of data object copies and set the write-run length at 2. In this case the invalidation algorithm outperforms the broadcast algorithm. As a result, we conclude that the increase of the data objects number has an impact on the performance.

The results shown in Fig. 3 demonstrate that the appropriate algorithm depends on application-layer parameters, thus showing the benefit of dynamically adapting the OS based on application-layer knowledge. Based on the data obtained from these experiments and machine learning techniques, we plan to generate a function that predicts from these parameters which algorithm will be better to use. When there is a change in parameters, the adaptation framework will use this function to decide whether to switch algorithms.

Possible Extensions

As a following step, we plan to validate our approach by providing support for adaptation involving multiple layers. We will first take into account the different network requirements of the algorithms. We then plan to extend the analysis phase to take into account not only the application-layer parameters but also the current state of the network.

Another planned extension is investigating coordinated adaptation across both OS and application layers through an application that transfers audio over the network and processes it in a parallel way using distributed shared memory. The audio compression rate of this application can be adapted based on network conditions. The multi-level framework enables this adaptation to be implemented in coordination with OS-level algorithm changes, thus avoiding unpredictable results. Indeed, performing this adaptation in an isolated way, without taking into account OS bandwidth needs, may cause instabilities. For example, if there is no coordinated adaptation and the bandwidth is low, the application increases its compression rate, freeing up available bandwidth. This causes the operating system to use a more bandwidth-demanding algorithm, which reduces the available bandwidth and the cycle continues. A coordinated adaptation, as enabled by our approach, avoids this problem

6R e l a t e d w o r k

The dynamic adaptation of operating systems has been investigated under different angles. Seltzer and Small proposed a design for self monitoring and self adapting an extensible kernel [START_REF] Seltzer | Self-monitoring and self-adapting operating systems[END_REF]. They have defined a mechanism that uses online and offline analysis of the system behavior to adapt its policies. Teller and Seelam investigated the dynamic adaptation of OS policies in order to improve the application and the system performance. They have proposed a multi policy system that aims to satisfy different constraints by dynamically switching between policies. They have demonstrated their concept in the context of I/O scheduling [START_REF] Teller | Insights into providing dynamic adaptation of operating system policies[END_REF].

The above works do not apply the service-oriented approach for structuring their OS. In this context Milanovic and Malek [START_REF] Milanovic | Service-oriented operating system: A key element in improving service availability[END_REF] investigated possible architectures to integrate the OS and the service-oriented computing. Among the investigated architectures, the closest one to our approach is to consider the OS as a set of collaborating services. The main difference with our model is the way of building the kernel. They propose to build it according to available kernel services while we propose to build it by selecting the most appropriate algorithm for each service according to its environment. Another interesting project in this context is SoOS, which tries to extend the limited capabilities of local devices with remote ones [START_REF] Schubert | Service oriented operating systems: Future workspaces[END_REF]. In this project local and remote resources are virtualized on top of the networking and Input/Output interfaces and are then made transparent to the running applications. Thus, the OS is considered as a collection of independent computational entities that appears to its users as a single coherent system. Unlike our proposal, this work focus on the physical infrastructure layer.

Recently, multi-level adaptation frameworks have also been investigated. The authors of [START_REF] Bratanis | Slas for cross-layer adaptation and monitoring of service-based applications: a case study[END_REF] and [START_REF] Kazhamiakin | Cross-layer adaptation and monitoring of service-based applications[END_REF] aim to provide a coherent solution to monitor and adapt service-based applications. The considered service layers are : business service, composition and coordination service, and the infrastructure service. The main difference with our proposal is the considered layers. They focus exclusively on service-oriented application layers while we also include the OS layer.

Conclusion

Modern distributed operating systems must operate reliably in constantly-changing environments. As a result, dynamic adaptability is an essential requirement for such systems. This paper made two contributions in this context. First, it proposed a service-oriented approach for building distributed OSs, which enables adapting the OS algorithms in a dynamic, on-demand fashion. Second, it proposed using a common framework to adapt both the OS and application layer in a coordinated way, thus avoiding any possible conflict or redundancy.

To demonstrate the usefulness of the proposed architecture, the paper investigated the DSM functionality and provided examples and experimental results using the Kerrighed distributed OS. As future work, we plan to implement further adaptation examples involving real applications on Kerrighed and extending the SAFDIS adaptation framework to control and adapt both layers.

Acknowledgment

We cordially thank the late professor Françoise André for her precious ideas and comments which helped to define the objectives of this work. We thank Dr. Louis Rilling for providing documentation support concerning Kerrighed. Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://www.grid5000.fr).

7 <Fig. 2 .

 72 Fig. 2. Multi-level dynamic adaptation architecture -solid arrows represent the communication between sensors of each layer (sensors are not shown in the figure) and the Monitoring component of the framework; dashed arrows represent the communication between actuators of service-oriented layers (actuators are not shown in the figure) and the Execution component of the framework; dotted arrows represent the sequence of phases in an adaptation cycle.


Yves Rozenholc
email: yves.rozenholc@parisdescartes.fr

Markus Reiss
email: mreiss@mathematik.hu-berlin.de

Daniel Balvay
email: daniel.balvay@yahoo.fr

Markus Reiß

Charles A Cuenod
email: charles-andre.cuenod@egp.aphp.fr

Growing time-homogeneous neighborhoods for denoising

research documents, whether they are published or not. The documents may come

Introduction

Dynamic Contrast Enhanced Imaging (DCE-imaging) is beginning to be widely used in medical imaging of brain structures or cancerous tumors [START_REF] Goh | Functional imaging of colorectal cancer angiogenesis[END_REF][2][3][4] [START_REF] Padhani | Angiogenesis imaging in the management of prostate cancer[END_REF] [START_REF] Bisdas | Dynamic contrast-enhanced CT of head and neck tumors: perfusion measurements using a distributed-parameter tracer kinetic model. Initial results and comparison with deconvolution-based analysis[END_REF].

The follow-up of the evolution of a bolus of contrast agent injected during sequential imaging acquisition with Computed Tomography, Magnetic Resonance Imaging or Ultrasound imaging (DCE-CT, DCE-MRI or DCE-US) allows the in vivo assessment of microcirculation parameters such as tissue-blood perfusion, permeability, blood and interstitial volume [START_REF] Miles | Functional CT imaging in oncology[END_REF] [START_REF] Tofts | Modelling tracer kinetics in dynamic Gd-DTPA MR imaging[END_REF][8] [START_REF] De Bazelaire | Comprehensive model for simultaneous MRI determination of perfusion and permeability using a blood-pool agent in rats rhabdomyosarcoma[END_REF]. Taking into account the Arterial Input Function (AIF), such estimates quantify locally the tissue dynamic characteristics [START_REF] Brochot | Evaluation of antiangionenic treatment effects on tumors' microcirculation by Bayesian physiological pharmacokinetic modelling and magnetic resonance imaging[END_REF] [START_REF] Krishnamurthi | Functional imaging in small animals using X-ray computed tomography-: study of physiologic measurement reproducibility[END_REF]. Thus DCE-imaging has great potential for cancer detection and characterization, as well as for monitoring in vivo the effects of treatment [START_REF] Jain | First-pass perfusion computed tomography: initial experience in differentiating recurrent brain tumors from radiation effects and radiation necrosis[END_REF][13] [START_REF] Rosen | Dynamic contrast-enhanced magnetic resonance imaging for assessing tumor vascularity and vascular effects of targeted therapies in renal cell carcinoma[END_REF] [START_REF] Zhu | Early antiangiogenic activity of bevacizumab evaluated by computed tomography perfusion scan in patients with advanced hepatocellular carcinoma[END_REF]. Clinicians are now expecting such information both at high resolution (up to pixel-level) and displayed as parametric maps [START_REF] Sorensen | Contrast agents in functional MR imaging[END_REF]. This is of particular interest in cancerous tumors, known to be heterogeneous with regions varying from non-perfused and necrotic to hypervascular "hot spots" [START_REF] Kiessling | Volumetric computed tomography (VCT): a new technology for noninvasive, highresolution monitoring of tumor angiogenesis[END_REF]. Recently, Cao et al. [START_REF] Cao | Developing DCE-CT to Quantify Intra-Tumor Heterogeneity in Breast Tumors With Differing Angiogenic Phenotype[END_REF] showed the ability of DCE-CT to assess intratumor physiological heterogeneity in tumors. This offers an in vivo tool for the evaluation and optimization of new therapeutic strategies. Thanks to the linear relationship between the concentration of the contrast agent and the attenuation, using DCE-CT instead of DCE-MRI for example, one can estimate the contrast agent concentration directly [START_REF] Axel | Cerebral blood flow determination by rapid-sequence computed tomography: Theoretical analysis[END_REF].

In the last decade, as CT techniques expose to X-ray, for ethical reason, an emphasize has been on reducing the total X-ray dose exposure. So-called "lowdose" CT have been developed to ensure a good (2D-or 3D-) tomographic Radon reconstruction when the Nyquist theorem is strongly violated. Techniques used in low-dose CT techniques use

• iterative reconstruction which rely on the accurate modeling of the distribution of noise in the acquired data [START_REF] Thibault | A Recursive Filter for Noise Reduction in Statistical Iterative Tomographic Imaging[END_REF];

• sparsity inside the tomographic image and compressed sensing [START_REF] Donoho | Compressed sensing[END_REF][START_REF] Candes | Near-optimal signal recovery from random projections: universal encoding strategies?[END_REF] : Prior Image Constrained Compressed Sensing (PICCS) [START_REF] Chen | Prior image constrained compressed sensing (PICCS): a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets[END_REF] or Vastly undersampled imaging with projections (VIPR);

• correlations between voxels in time-dependent images : Highly Constrained Back Projection (HYPR) (see [START_REF] Mistretta | Undersampled Radial MR Acquisition and Highly Constrained Back Projection (HYPR) Reconstruction: Potential Medical Imaging Applications in the Post-Nyquist Era[END_REF] and references within);

• a regularization or an interpolation of the sinogram before reconstruction. [START_REF] Happonen | A Comparative Study of Angular Extrapolation in Sinogram and Stackgram Domains for Limited Angle Tomography[END_REF][START_REF] Wang | Penalized weighted least-squares approach to sinogram noise reduction and image reconstruction for low-dose X-ray computed tomography[END_REF].

Allowing under-sampled measurements during tomographic acquisitions, all these techniques offer the possibility to reduce the X-ray dose in CT by often a factor of 25 to 100. Implemented by the manufacturers as black-boxes inside the scanner, they work image per image or at least in a steady state of the contrast agent concentration.

In our context of DCE-CT and in the typical experiments imagined to track angiogenesis phenomena and micro-vascularization estimation, we want to benefit from a non-steady state of the contrast agent concentration and the sequence is made of about 100 acquisition times. These acquisition times correspond to the same amount of 2D-slices (or 3D-volumes) of CT which are located on the patient at the same position (up to potential movements of the patient). Hence even in the context of low-dose CT described above, decreasing the X-ray dose by acquisition time is of main interest as the patient will receive a total X-ray dose proportional to the number of acquisition times. In our context, it would be highly desirable to decrease the X-dose for each CT-acquisition by the same factor 100 in order to use so-called "hyper-low-dose" CT. As a consequence, such X-ray reduction produces highly noisy raw images : at each location of the (spatially unique) 2D-slice (or 3D-volume), the time series (dynamics) which represents the contrast agent evolution has a so low signal-to-noise ratio that it is hardly usable.

Multi-band filtration [START_REF] Chen | A Multi-Image Restoration Method for Image Reconstruction from Projections[END_REF], which uses multiple acquisitions from the same image in order to increase the signal to noise ratio by using a Maximum A Posteriori (MAP) approach, relies on the assumption that images do not change during acquisition. This hypothesis is clearly violated in DCE-CT due to the variations of the contrast agent concentrations that we aim to exploit.

After a Radon reconstruction performed at each acquisition time, the low signal to noise ratio at the pixel-level resulting from (high-frequency) hyper-lowdose CT acquisition may be improved by

• combining the dynamics over a large number of pixels in a Region Of Interest (ROI);

• in a time by time approach, using spatial averaging [START_REF] Buades | Nonlocal Image and Movie Denoising Int[END_REF];

• using filtering or smoothing techniques such as parametric modeling (usually derived from pharmacokinetic models) [START_REF] Barbier | Methodology of brain perfusion imaging[END_REF] [START_REF] Collins | Dynamic magnetic resonance imaging of tumor perfusion[END_REF].

Unfortunately, these procedures present the following drawbacks:

• mixing the complex dynamics of the contrast agent which may not be homogeneous, will lead to inaccurate parameter estimation;

• information loss due either to excess of smoothing out the dynamics or to misspecification of the model [START_REF] Balvay | Mapping the zonal organization of tumor perfusion and permeability in a rat glioma model by using dynamic contrast-enhanced synchrotron radiation CT[END_REF].

We aim to improve the signal to noise ratio of DCE-CT after the Radon reconstruction of each acquisition time. Because the DCE-CT sequences, that we focus on are made from about 100 CT-slices at the same location, it would be highly desirable to reduce the X-ray dose of each slice by the same 100 factor in order to be equivalent to a single "low dose". Clearly, our aim is not to replace the techniques used to achieved "low dose" CT but to complement such techniques by taking advantage of the dynamical acquisition in order to go further in the reduction of the X-ray dose which remains an overall concern more longitudinal follow-up or pediatric radiology. Here the increase of signal to noise ratio is achieved in two steps using a global procedure, based on statistical multiple testing to compare the dynamics observed at each voxel of the image. Using a test of zero mean for the difference between enhancements or estimated enhancements, neither specific modeling assumption of the dynamics are made nor needed. Being adaptive to the unknown smoothness of the curves, multiple testing have enough power to properly differentiate between dynamics.

The first step is a voxel-wise denoising step : at each voxel, the dynamic is compared to those from the voxels adjacent to the selected location. Homogeneous dynamics are aggregated iteratively in order to build a denoised version of the dynamic coming from the selected voxel. At each iteration, using multiple tests, the dynamics from the next closest voxels around the selected voxel, are compared to previous denoised versions in order to check for homogeneity. As a result, around the selected voxel, a neighborhood of locations is built which aims to satisfy the bias-variance paradigm of non-parametric statistics : being large enough to reduce the variance and improve the signal to noise ratio while not being so large as to endanger homogeneity between selected voxels. This step provides high resolution information with an improvement in signal to noise ratio which allows, for example, the construction of parametric maps.

The second step automatically builds a clustering of all the dynamics, preserving dynamic homogeneity in each cluster. The clustering procedure uses as input the neighborhoods grown in the denoising procedure and their associated denoised enhancements. This clustering depends neither on an a priori knowledge of the number of classes as does the k-means algorithm [START_REF] Jain | Algorithms for Clustering Data[END_REF] nor on any modeling assumption of behavior inside a class as does the EM algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF] in the context of Gaussian mixture. It is based on the same multiple tests as those previously used in the denoising step. The clustering procedure constructs classes made from a set of voxels with enhancements statistically close to an estimated enhancement (built as a class centroid).

The clusters can be viewed as an automatic ROI-partition selection with respect to the typical behavior of the dynamics. This automatic ROI-partition provides a summary of the full dynamic sequence into few typical denoised dynamics which preserves the heterogeneity of the tissues in general and of the tumor in particular. It is similar to a piecewise constant representation of a function having functional values.

From a clinical point of view, this denoising and clustering procedure is a necessary step to allow a relevant evaluation of the microcirculation in vivo by using pharmacokinetic models on a pixel by pixel basis.

Only one hyper-parameter is used in these two stages, namely the level of the multiple test. In the denoising step, it controls how easily the dynamics are aggregated. In the clustering procedure, it plays a role equivalent to a penalization ensuring an adaptive control on the number of classes.

The article is organized as follows: we first introduce the statistical framework and comment on the assumptions in Section 2. In Section 3 we outline our denoising method for DCE-CT. The automatic and unsupervised classification of tissues is tackled in Section 4. Finally in Section 5, we show an application to DCE-CT data on liver metastases and some simulations in order to validate our methods. The two main statistical tools used to construct our denoising procedure are described in the Appendix : multiple testing in 7.1 and neighborhood/ring growth in 7.2.

Statistical framework

For the rest of this paper, we consider a DCE-CT sequence as a finite sequence of noisy images indexed by both time and space:

I = {I x (t), x ∈ X , t ∈ {t 1 , t 2 , . . . , t K }} ,
where I x (t) denotes the noisy enhancement at time t and voxel location x. We will denote by I x the vector (I x (t 1 ), . . . , I x (t K )), of the discretely observed dynamics at location x for all observation times t 1 ,. . . , t K . Here X denotes the finite voxel grid and t j , j = 1 . . . K are the acquisition times.

We assume that the observable gray level I x (t) may be written as

I x (t) = i x (t) + σε x (t), (1)
where i x (t) denotes an unobservable true gray level, ε x (t) denotes a standardized noise and σ the noise level. The noise ε x (t j ) are supposed independent with respect to both space x and time t j . While time independence of the noise may be easily justifiable, the assumption of independence between locations is known not to hold in CT. In particular, it is well known that CT creates radial artifacts. We use this assumption as a simplified model. A more complete model could be written as:

I x (t) = i x (t) + η x (t) + σε x (t), ( 2
)
where η is a time-independent but space-correlated noise which models the spatial artifacts.

Even if these spatial artifacts are clearly visible to the human eye, the noise η may be considered as negligible from a statistical point of view. In the framework of DCE-CT where, because of low X-ray dose, σ is clearly large, our results (see Section 5) indicate that the simplified model (1) suffices.

To simplify our the presentation, we will assume in the following that the parameter σ is known. The knowledge of σ may be ensured by either a proper calibration of the scanner or the construction of an independent estimate. On one hand, assuming that the distribution of the errors ε x (t) is Gaussian, a natural estimate is based on the mean of the residuals. On the other hand, a simple histogram of the DCE-CT noise (see Figure 1) shows that this noise distribution is closer to a Laplace distribution. In this context, median-based estimate is more advisable (see [START_REF] Reiß | Pointwise adaptive estimation for robust and quantile regression[END_REF]) and we have used this in practice (see Section 5.2). Furthermore, from an industrial point of view, it is possible to avoid extra calibration or noise level estimation by using the result in [START_REF] Durot | An adaptive test for zero mean[END_REF]. In such case the only assumption needed is that, for any couple of locations x and y, the distribution of ε x (t)ε y (t) is symmetrical without mass in 0 for all time t.

Denoising DCE-CT

We outline here the statistical procedure used to denoise DCE-CT. This procedure is summarized by the flowchart in Figure 2 and Algorithm 1 in Appendix 7.3. It is based on the two statistical tools introduced in Appendices 7.1 and 7.2.

Our method is based on the statistical comparison between two observed enhancement vectors I x and I y at spatial locations x and y as presented in Appendix 7.1. Two enhancement vectors are considered indistinguishable if their difference does not deviate significantly from the zero vector. This is controlled by a statistical multiple test for whether the difference vector I x -I y has mean zero vector. In this case, we will call I x and I y statistically close or time homogeneous and write I y ≡ α σ 2 I x where α refers to the test level and σ 2 to the noise variance.

Such comparisons of vectors using multiple test procedures have been developed in the Gaussian framework by [START_REF] Spokoiny | Adaptive hypothesis testing using wavelets[END_REF][START_REF] Baraud | Adaptive tests of linear hypotheses by model selection[END_REF] and for heterogenous symmetric noise by [START_REF] Durot | An adaptive test for zero mean[END_REF]. The power of these test procedures is known to be adaptive to the regularity of the underlying signal. Hence, we do not need to specify the behavior of the enhancements in order to be able to deal adaptively and automatically with their particularities. We only require that the enhancements are not too wild in the sense that their differences should have a certain minimal (Hölder ) regularity [START_REF] Wasserman | All of nonparametric statistics[END_REF]. The latter is clearly satisfied in the context of contrast enhancements in DCE-CT.

The use of differences ensures that noise can be assumed to be symmetrically distributed, thus avoiding typical problems that spring from known features in tomographic sequences as described by [START_REF] Gravel | A method for modeling noise in medical images[END_REF].

At each spatial location x ∈ X called center, we aim to construct a spatial x-neighborhood V x made of voxels y ∈ X such that I y is statistically close to I x and such that the statistical error may be controlled. The estimated enhancement vector at location x is then derived from this neighborhood as a centroid given for example by the empirical mean or by the generalized median (see Eq. ( 3) and (4) in Appendix 7.2). Controlling the statistical error means that we aim to include as many voxels as possible in V x to reduce the variance while the bias due to over-smoothing is kept small.

To that end, as presented in detail in Appendix 7.2, at each spatial location (called center), an increasing sequence of "time homogeneous" neighborhoods of voxels is grown according to the following steps: (i) voxels spatially close to the center with enhancement statistically close from the enhancement of the center are selected; (ii) a sequence of estimates -each built on one neighborhoodhaving decreasing variance is built; (iii) from this sequence, the largest statistically acceptable neighborhood is selected. The estimate associated to this last is expected to realize the bias-variance tradeoff.

Step (i) described more in detail in Appendix 7.2 is a pre-selection step which allows the denoising procedure to obtain "neighborhoods" built from different objects of the same type over long distances. This is useful in medical images where the same kind of tissue (e.g. small arteries, see Figure 5-a) may reappear in different areas. At step (ii) for each neighborhood, taking a centroid (mean or median) of the enhancements associated to the voxels included in the spatial neighborhood provides a enhancement estimate. The last step (iii) involves a generalization of the above multiple test in order to decide whether or not two estimated enhancements ÎV and ÎW , constructed on two disjoint sets V and W , are statistically indistinguishable.

Because of the neighborhood grow, the known variances of these estimates decrease and can be used to select an estimate that achieves a good statistical trade-off between bias and variance in step (iii) [START_REF] Reiß | Pointwise adaptive estimation for robust and quantile regression[END_REF]. The precise description of the neighborhood growth is described in Appendix 7.2.

Clustering

Using the test procedure denoted by ≡ α σ 2 introduced in the previous section, we develop an automatic clustering procedure. It does not rely on the knowledge of the number of classes as in the k-means algorithms [START_REF] Hartigan | Clustering algorithms[END_REF][START_REF] Jain | Algorithms for Clustering Data[END_REF]. In addition, it does not need the Gaussian assumptions, often made to describe the behavior inside a class. The classes are built in order to keep statistical homogeneity. Let us denote by C the set of all classes result of the clustering procedure. The estimated enhancement associated to a class c ∈ C, called the center of c, will be denoted Îc . The clustering algorithm is a recursive algorithm with a main loop which can be decomposed in three main steps.: The algorithm stops when all single voxels have been assigned to a class and when classes cannot be merged anymore.

Using the construction described in Appendix 7.2, Step (a) can be done in an efficient way. Suppose that, using the denoising procedure described in Section 3, we have grown at each voxel location x of X a neighborhood V x . Given a list L ⊂ X of voxels of interest (a Region Of Interest (ROI) or the full image), we define the child/ancestor relation by "z y if z ∈ V y " and call "children of x in L" the set

N L x = {z ∈ L, ∃y 1 , ..., y k ∈ X s.t. z y k . . . y 1 x}.

The list N L

x of all children of x is constructed in the function Children, see Algorithm 3 in the Appendix.

While the list L is non empty, the following four steps are made:

Step 0 -Next : Consider one voxel x in the list L with the largest neighborhood size

|V x |. Compute its child list N L x and set c = N L x .
Step a -Class construction : From c build a centroid Ĵ (see Appendix 7.2). A new (possible) class c Ĵ is defined by its center Ĵ and the list c Ĵ of voxels y in c such that Îy ≡ Step c -Update Cluster List : The class c Ĵ is added as a new class to the cluster list.

The clustering procedure is described in Algorithm 2 in the Appendix 7.3.

We now give the details of the robust center construction. Given a set c of voxels, we iteratively define Ĵ = Îc and update c as the set of the y in N L x such that the test Îy ≡ α σ 2 ρ(1,min(|Vy|,|c|)) Ĵ is "strongly" accepted. This procedure, which can be viewed as a robust k-means procedure with k = 1 converges quickly and only requires a few iterations.

Remark 1: In this setting, the larger class attracts the smaller class in its statistical surrounding defined by a "noise level radius" of order σ 2 around its center. Moreover the center of the smaller class is located up to a statistical distance of order σ 2 / min(|c|, |c ′ |). Hence, the two classes are merged if the "statistical distance" between their centers is less than ρ(1, min(|V y |, |c|)) (see Appendix 7.2). Remark 2: Thanks to the decreasing order in the neighborhood sizes, the last unclusterized voxels are mainly due to body movements and are far from all previously constructed centers. Because of this observation one can stop the clustering algorithm when the neighborhood size is less than a prescribed v 0 . The remaining voxels then define an extra class of unselected voxels. This extra class of movements can be used as a prior input for a registration algorithm.

Results

Data material

We use a DCE-CT sequence of 53 images obtained in 90 seconds at the same level of the upper abdomen on a patient who was asked to hold his breath. Acquisition parameters were fixed at 80 Kv and 50 mAs. These images were obtained with an in-plane resolution of 512×512 pixels. The sequence is split into three periods of 30-second breath-hold separated by 8-second pauses to allow the patient to breathe. Each period is characterized by the time delay between two images: one second between images for the first 30 images; two seconds for the next 15 images and three seconds for the last 8 images. After the beginning of the acquisition, at time 3 seconds, an intravenous bolus of 80 ml of Iobitridol (Xenetix, Guerbet, France), an iodinated contrast medium, is injected at the rate of 4-5 ml/s and followed inside the tissues through their enhancements. At time t, the contrast agent concentration is proportional to the difference of gray levels between time t and time 0. On CT images, gray-level ranges of the different tissues, which are measured in Hounsfied Units (HU), are related to their chemical composition (e.g. fat, air, bone, water) and their content in contrast media. On the images, the liver, spleen, aorta, stomach, a vertebra and some blood vessels are visible. However, the images suffer from a poor signal to noise ratio, due to the limited radiation dose used for each acquisition time. Figure 4 shows typical enhancements obtained from voxels in manually selected ROI (of size larger than 100 voxels) within the aorta (top), known to be homogeneous, and within the tumor (bottom), where heterogeneity is expected. In order to provide a visual idea of the noise level and the variability inside a manual ROI, for each selected ROI, we have constructed a mean curve, and the 5% and 95% quantile curves by considering each time separately and we have drawn 5 curves of randomly selected individual voxels.

Using the first four images without enhancement (baseline images), we have derived a sample of the noise distribution as follows: (1) compute a (voxelby-voxel) baseline median image from the first four images; (2) compute the residual images of the differences between each baseline image and the baseline median image; (3) use the values inside the residual images as a sample (of size 4×512×512) of the noise distribution. We expect in step (1), if there are only a few movements that, in each image, most of the voxels represent the same tissue and the variations in their gray levels are just due to noise. The histogram of the noise distribution obtained from this sample is represented in Figure 1 together with the fitted Laplace distribution.

Results for the denoising procedure

In Appendix 7.2, we fix the increasing sizes n i of the neighborhoods to be successively 1, 5, 15, 38, 91, 211, 476. As our procedure involves several multitests, we use a False Discovery Rate (FDR) approach [START_REF] Benjamini | Controlling the False Discovery Rate: a practical and powerful approach to multiple testing[END_REF] calibrated to take into account the multiplicity due to time (see Appendix 7.1) and due to neighborhood growth (see Appendix 7.2). We have tried 3 different values of σ obtained in the following ways: The estimated curve shows a quicker and larger enhancement providing a stronger initial slope that shows a perfusion and a more important blood volume than for the voxel within the tumor itself (left column). See Figure 5-a for more information.

• using the noise distribution sample made from the baseline images (see Section 5.1), a Laplace distribution has been fitted as shown on Figure 1 : its parameter is 1/29 corresponding to a standard deviation σ = 29 √ 2 ≈ 41 HU.

• using the same sample, the classical estimate of the standard deviation gives σ = 60 HU.

• a value σ = 50 HU has been used as an "in-between" mean value.

The difference between these values may be explained by the heavy tails due to tomographic artifacts and movements. Due to the maximal size (476) of the neighborhood, we may expect for the minimal residual noise to be σ/ √ 476 ≈ 2.75. Hence, in the worst case considering that σ = 60HU , our maximal gain in the signal to noise ratio is about 60/2.75 ≈ 22.

Because our procedures rely on tests, a good trade-off between the choice of σ and α should be made. The value σ = 41 HU, which ideally fits our model, does not denoise sufficiently for any α as it leads to neighborhoods which are too small, where growth has been stopped artificially by tomographic artifacts. For the value σ = 50 HU, this phenomenon remains at least in the lower left corner of the sequence. With σ = 60 HU and values of α in the range 0.2 to 0.001 our algorithm behaves stably, leading to similar denoised sequences and similar clustering showing more details as α gets larger. With σ = 60 HU the sequence is properly denoised and most of the tomographic artifacts are removed.

The clustering procedure using σ = 50 HU or σ = 60 HU with the corresponding α value for the tests involved in the clustering procedure leads to similar and stable segmentations. The parameter α is a good tuning parameter which may be used by doctors to tune the clustering in order to focus on the desired level of detail. It is remarkable that, when using this statistical tuning parameter α, the geometry of the clusters is stable and varies slowly while preserving the morphological details.

In what follows, we use the setting σ = 60HU and α = 0.05 to illustrate both the denoising and the clustering procedures.

First in Figures 3, we present the result of our method on two images of our DCE-CT sequence at time 25s and 45s, respectively : (top) the original image; (middle) the denoised image and (bottom) the residuals (i.e. the difference between original and denoised images) . At time 25s, the image shows a clear enhancement in the aorta and in the tumor. At time 45s, the details of the return to the venous system are visible. On both images, the same quality of details can be seen. Morphological information such as shapes and borders of the organs are clearly visible on the denoised image, while most of the noise has been removed. One can notice two types of structures on the residual image. The first comes from movement (see for example on the ribs in the upperleft corner or around the stomach) as we apply our method on the sequence without using any kind of (pre-) registration or motion correction algorithm. Our algorithm suffers from these movements on the borders of moving tissues and should be used with a registration procedure in case of strong movements. The second comes from the CT radial artifacts. In the model given by Eq. ( 2), the noise η is also removed when using our procedure with σ = 60HU . This is remarkable, as nothing, from a theoretical point of view, has been used to identify this radial artifacts. Tomographic artifacts are organized in an image along directions which are not necessarily the same at different times. Thanks to the use of the whole temporal structure in testing, the directions are mixed and the tomographic artifacts do not have a significant effect on the denoising result. The benefit of using an approach which compares not only voxels in one image locally but uses the dynamics of these voxels in the full sequence is clear : one can apply a stronger denoising procedure without losing any details in the 2D structure. In Figure 5-a, the neighborhood V x of the voxel x is not connected. This voxel x has been selected in a small artery covering only a few voxels and most of its neighborhood is included in the aorta and made from voxels disconnected with the original location. This is not a surprise, as physiologically, the flow in the artery system does not allow the enhancements in different arteries to be distinguished at our time resolution.

These figures show the second benefit of our approach : the sequence is not denoised in a "slice-by-slice " approach as a movie but as a single image showing dynamics. Hence, it is not the gray levels of images which are denoised but the full dynamics, over time, at each voxel. This approach clearly improves the signal to noise ratio by reducing the noise without changing the signal and one can then proceed to a proper inspection of kinetic enhancement curves to derive the characteristics of the underlying physiological processes.

It is also worth noting that the maximum of the estimated enhancements in the hot-spot inside the tumor (Figure 5-b) is larger than that associated to the surroundings of the tumor (Figure 5-b). This is in accordance with clinical knowledge as the flow in such voxels is larger. We present here the classifications obtained by our automatic clustering procedure described in Section 4 applied to the denoised sequence obtained with the setting σ = 60 HU and α = 0.05 (see Section 5.2 for discussion about the choice of these parameters). In Section 4, the voxels with neighborhoods of size smaller than 60 voxels, which correspond to movements, are not classified. They appear in Figures 6 and7 in dark blue. The final clustering of this series contains more than 600 clusters. Most of these clusters are small and correspond to movements.

Results for the clustering procedure

Figure 6-(top) shows (with artificial colors) the result of our procedure applied to the enhancements in order to focus on functional phenomena. The enhancements are obtained after removing, at each voxel, a baseline intensity estimated on the first 5 times prior to the contrast agent injection. Only 30 clusters have a size larger than 60 voxels. In this figure, the 30 clusters are shown and nine typical clusters are localized. Air and tissues without enhancement are localized in Cluster 4. Cluster 1 mainly consists of the liver. The aorta is well distinguished (Cluster 3). Hepatic veins correspond to Cluster 5 and smaller vessels to Cluster 6. The tumor is split into four dynamical behaviors (7, 2, 8 and 9) extending from tissues at the periphery, which are compressed by the tumor to hot spots. Due to similar enhancement patterns in reaction to the specific contrast agent used to obtain the DCE-CT sequence, the tissues of the spleen and the stomach are not properly distinguished. While Cluster 5 corresponds to the center of the hepatic veins, Cluster 6 which contains smaller vessels also contains the borders of the hepatic veins. This may be due to a partial volume effect or due to surface tension in large vessels: the flow on their borders is not as fast as in their centers.

The centers (estimated enhancement) of these nine relevant clusters are presented in Figure 6-(bottom) with their associated label number. These curves show a very high signal to noise ratio obtained without any (time) regularization and are a proper summary of the dynamic information that exists in the sequence. From this summary, one can easily derive the characteristics of the underlying physiological processes.

Figure 7 shows (with artificial colors) the result of our procedure directly applied to the original image, without baseline removal. This strategy leads to a differentiation of tissues with a different baseline even if they have the same enhancement. For tissues having sufficient enhancement, the clustering is similar to what is obtained in Fig. 6 (see tumor tissues [START_REF] Goh | Quantitative colorectal cancer perfusion measurement using dynamic contrast-enhanced multidetector-row computed tomography: effect of acquisition time and implications for protocols[END_REF][START_REF] Tofts | Modelling tracer kinetics in dynamic Gd-DTPA MR imaging[END_REF][START_REF] Brix | Microcirculation and microvasculature in breast tumors: pharmacokinetic analysis of dynamic MR image series[END_REF][START_REF] De Bazelaire | Comprehensive model for simultaneous MRI determination of perfusion and permeability using a blood-pool agent in rats rhabdomyosarcoma[END_REF] for example). The other tissues (4) are here distinguished through their baseline, offering different information to the physiologist. In this clustering, 45 clusters have more than 60 voxels. Only these 45 clusters are displayed.

Simulation results

In order to evaluate our method, we have run artificial experiments on simulated data for the ideal case when no movements occur during the sequence. We have built a synthetic dynamical image (120×60×64) without baseline with 64 observation times using the synthetic enhancements shown at the bottom of Figure 8-a. We have chosen time 30 to present our synthetic dynamical image and to illustrate our simulation results. At the top of Figure 8-a, the homogenous regions are surrounded by artificial fine white lines in order to help the visualization of the pattern. These homogenous regions represent air (0), normal tissue [START_REF] Goh | Functional imaging of colorectal cancer angiogenesis[END_REF], two types of venous systems (2, 3) and the aorta [START_REF] Brix | Microcirculation and microvasculature in breast tumors: pharmacokinetic analysis of dynamic MR image series[END_REF]. The large ball in the upper-right corner of the Figure 8-a represents a tumor with (dynamical) properties that vary continuously in a piecewise linear fashion from the center to the border. Inside this synthetic tumor, enhancements may be computed from the radius using the profile given by the thick line in Figure 9 and enhancements 4 to 7 in Figure 8-a. At radius 0.7 for example, the profile yields the value 5.2 which corresponds to a mixture of enhancements 5 and 6 in proportion 20% and 80%, respectively.

From this synthetic dynamical sequence, we have constructed a simulated noisy version by adding (time and space) independent Laplace noise with standard deviation 41HU as shown, at time 30, in Figure 8-b.

Our denoising algorithm (see Section 3) has been applied to this sequence with σ = 60HU and α = 0.05 using the same values as in Section 5.2. The result of this denoising step is presented (at time 30) with the same gray-scale as the observations in Figures 8-c and (to show the gain in contrast) with its full gray-scale range in Figure 8-d. To help the comparison, this last gray-scale is also used in the top Figure 8-a which shows the original data.

Considering the amount of noise in Figure 8-b, any slice-by-slice denoising procedure will fail to show the weak enhancements in the synthetic venous systems (2 and 3) already visible in Figure 8-c and even clearer in Figure 8-d.

Using the same grayscale, the result of our clustering algorithm (see Section 4) applied to this sequence is presented at the top of Figure 8-e with the associated estimated enhancements inside each cluster shown at the bottom of the same Figure . The estimated clusters and their associated estimated enhancements (centroids of each cluster) are shown and linked by capital letters. The enhancements are well recovered as shown by the comparison of the bottom part of Figures 8-a and8-e. The clustering of the synthetic tumor follows the circular geometry used to build the synthetic tumor in this sequence. This clustering has to be understood as the construction of an adaptive piecewise constant map with unknown number of steps. The selected number of steps depends on the underlying enhancement, on the standard deviation of the noise and on the choice of α, the level of the test procedure, and on the unknown localizations. In Figure 9, the fine dotted line shows a possible stepwise profile which could be associated to this clustering: each step represents the 4 estimated enhancements D to G along a diameter of the synthetic tumor.

The choice of a continuous synthetic tumor may seem surprising but it is driven by two ideas: on the one hand continuous changes exist inside tumors, and on the other this condition is more challenging for the algorithm. In addition, we have performed simulations (not presented here) with a synthetic tumor showing either a piecewise constant profile or split into one continuous part and one piecewise constant part. In each case, the estimation and the clustering showed the same quality of behavior. The estimation of each enhancement clearly improved as soon as the synthetic tumor exhibited some constant part in its profile. The abscissa is the radius and the ordinate defines the enhancement mixture: a value 5.2 obtained at radius 0.7 represents a mixture of 20% of enhancement 5 combined with 80% of enhancement 6. The fine dotted line represents a feasible piecewise constant representation of this profile.

Discussion

The goal of DCE-imaging is to create functional images (pictures) at voxellevel and not to compromise by using ROI which can mix several functional behaviors and so result in loss of information. For example, the heterogeneity of tumors is believed to provide important diagnostic information which may be properly evaluated using DCE-CT. As we remarked in the introduction, it is still desirable to reduce X-ray doses even if CT acquisition techniques continue to improve. Therefore images with poor signal to noise ratio, hardly usable to the proper evaluation of micro-circulation parameters, will still be produced. This is why we propose our two step procedure, based on the same multiple tests to compare random vectors, in order to improve the signal to noise ratio of the DCE-CT sequences. The first step denoises the dynamical sequence voxelby-voxel. The second step builds a spatial segmentation of the tissues based on the differentiation of the full dynamics.

The denoising procedure constructs around each voxel a neighborhood of (dynamically) homogeneous voxels with a size related to both the size of the homogenous tissue it belongs to; and the maximal gain in signal to noise ratio provided by the user. Through this step, at each voxel location, an individual denoised enhancement is obtained, providing complete denoised information on the heterogeneity of the dynamics. The size of each neighborhood, which plays the role of a window or bandwidth, is chosen adaptively, thus preventing underor over-smoothing.

In addition to this local approach, the clustering step pulls together the homogeneous voxels (enhancement or dynamics) in order to create a synthetic map (segmentation) of dynamical behaviors that sums up all the information contained in the DCE-CT sequence. Although an over-reduction of the information in this clustering step might be feared, the amount of detail found in the shape of the organs and the good level of heterogeneity in the tumor image show the proper behavior of this technique. The clustering provides an adaptive piecewise constant representation of the dynamical behaviors for which clear steps at the borders of organs and more continuous changes inside the tumor may be seen. This is illustrated by the dynamics 2, 7, 8 and 9 of Figure 6, which sum up in a piecewise way the vascularity of the tumor, ranging from low vascularization to hot-spot.

In an image-by-image (e.g. time-by-time) approach, an adaptive "bandwidth" selection aims at choosing the largest subset of pixels with statistically similar intensities in order to improve the signal to noise ratio. For DCE-CT, since each individual voxel is too noisy, such image-by-image approach fails. Given a voxel x, we propose here to find voxels with statistically similar enhancements: those enhancements together provides a estimated enhancement at voxel x for the all sequence of acquisition times. Clearly, we avoid a (time-bytime) maximal approach and build potentially smaller (at some acquisition time) dynamical homogenous neighborhoods. Hence, around a voxel x, the selected voxels form a large subset with respect to the size of the homogenous tissue they belong to, which is not necessarily large with respect to the gray-level obtained at a specific time.

Going further, it would be of interest to take advantage of the dynamical homogeneity of the tissues before the Radon reconstructions (of all acquisition times) by, for example, either comparing statistically partial sinograms or by looking for a sparse representation of the full dynamics. We have not addressed this issue because usually scanners automatically provide one reconstructed image per acquisition time without giving access to the Radon reconstruction used internally. Our expectation is that the (black-box) Radon reconstruction has already benefited from the larger local homogeneity inside each image. Applying a dynamical approach like ours directly during the Radon reconstruction in order to benefit from dynamical homogeneity remains a challenge.

During the long acquisition period needed in DCE imaging, movements are unavoidable because either patient or peristaltic movements. In Figures 6 and7, we clearly see navy bands at the periphery of the organs (interface regions). They are characteristic, easy to identify and remarkably elongated spatially. Most movements should be removed by prospective or retrospective techniques like respiratory gating or registration. These different techniques do not conflict and can complement each other: registration techniques can be improved by the knowledge of these navy bands and a denoised dynamics, or part of it, could be the target of a registration algorithm. Denoising and clustering based on the comparison of entire dynamics are not strongly affected by small movements. This is a promising avenue for further developments combining a dynamical approach with techniques that deal with movement.

For clarity, we have used a homoscedastic setting with a known constant noise level σ assuming independence in time and in space. This assumption is certainly not perfect. First, the noise level may depend on the intensity which varies not only among the tissues but also with time as the contrast agent goes through the tissue leading to contrast beam hardening artifacts. Second, CTartifacts are not spatially independent. Surprisingly, although our construction ignores CT-artifacts, they are properly removed without affecting the overall quality by using a larger σ than that prescribed directly by the estimation. This is an advantage of our approach that combines spatial locations, preserving time structures thanks to the use of statistically powerful multiple tests. The first issue above is addressed by [START_REF] Durot | An adaptive test for zero mean[END_REF] who propose a test of zero mean in random vectors even when the noise distribution presents strong departures from the simpler Gaussian homoscedastic case: the noise distribution needs only to be symmetric with an unknown noise level which may even vary in time. As our construction is based on differences, this symmetric assumption is fulfilled and ensures that the hardening artifacts will be taken into account. Even if the results we obtain are already satisfying, an industrial application will further benefit from such an implementation.

In the simplified setting we have used, α and σ play the role of tuning parameters. Using the baseline image, parameter σ could be evaluated as we did or by taking a ROI outside the patient. Hence, the true remaining tuning parameter of our algorithm is α. In the denoising step, at each voxel, it controls the level of the test used to choose the proper neighborhood and controls how smooth spacially the denoised dynamical image (meaning the full sequence) is. In the clustering step, it controls the number of clusters used to describe all the tissues. It plays a similar role as a penalty in the description of a function by a piecewise constant function. Decreasing α will produce fewer clusters preserving the easiest locations where the behaviors of the tissues are constant or slowly varying. In the case of DCE-CT, it will preserve large organs and strong details like the venous system, aorta, bones, etc. Inside a heterogenous area like the tumor where the profile varies slowly, the clustering will produce fewer details when α decreases but will keep a realistic piecewise constant description of the profile and hence of the tumor (see Figure 9). Changing the value of α in the clustering will mostly affect only the description of the navy bands.

Our clustering procedure does not require knowledge of the number of classes and a proper initialization as the k-means algorithm does [START_REF] Hartigan | Clustering algorithms[END_REF][START_REF] Jain | Algorithms for Clustering Data[END_REF]. In addition, it does not rely on assumptions about the behavior inside the classes like the Expectation Maximization algorithm does [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF]. Nor does it create a synthetic representation of the enhancements based on linear combinations of a few artificial enhancements, eigenvectors of a functional PCA [START_REF] Ramsay | Functional data analysis[END_REF]. In our clustering, the classes are built only to achieve a statistical homogeneity offering a good piecewise constant description involving a few realistic enhancements.

We have implemented the denoising algorithm (see Section 3) in Matlab c . The processing of a 512×512×53 dynamical image as presented in Section 5 takes three hours on a bi-3 Ghz Quad-core Intel Xeon PowerMac as it uses a loop on the 512×512 voxels. As the processing at one location does not involve processing at other locations, this loop may be highly distributed with a cheap parallelized implementation in open-CL on a Graphic Processing Unit (GPU): after a transfer in the shared memory of the GPU of the global information of the DCE-CT (here the 512×512×53 numbers encoded on 8 bits), an industrial implementation would send to each chip a number of the 512×512 voxels inversely proportional to the number of chips on the GPU.

It is worth noting that even if we have presented our technique on 2D DCE-CT, it could be straightforwardly generalized to 3D DCE-CT sequences using 3D neighborhoods instead of 2D neighborhoods. In such case, the result will be a clustering of the dynamical behaviors existing in the tissues of the 3dimensionnal image.

Conclusion

Using a two-step procedure, each step based on statistical multiple hypothesis testing, we introduce a novel algorithm to denoise and clusterize dynamical image, where each point of the picture codes a complete time series. Based on the comparison of the dynamics, these algorithms preserve their full structure. They do not rely on any a priori knowledge of the features in the image and run unsupervised using only one tuning parameter, which has a clear statistical interpretation as a significance level. The efficiency of our algorithm is shown on DCE-CT data used to follow the vascular and tissular distribution of a contrast agent.

The quality of the denoised dynamical image is shown by the many details which can be found by clinicians and allows a clear delineation of the tumor's heterogeneity. As a byproduct and without special tuning, the tomographic artifacts are removed. The clustering algorithm is based on the result of the denoising procedure. Efficient, it needs neither prior knowledge of the number of classes nor of the distribution inside a class nor of the modelization of typical enhancements. Here the statistical tuning parameter acts as a penalization which control the number of clusters. The result of the clustering provides Regions Of Interest that are meaningful at a physiological level. They automatically sum up the typical dynamics of the tissue behavior for further analysis of the microcirculation. Artificial experiments validate these results on simulated data.

Our new technique for denoising and clusterizing of DCE-CT by working directly with the time series offers the proper tools to allow estimation of microcirculation parameters. Coupled with registration techniques, this promising approach is a necessary step for in vivo evaluation and longitudinal follow-up of cancer tumors. Given j in 0 . . . d -1, let us denote the projection of Z onto the space generated by the vectors with same components on each time index T j l :

Π j Z = (m j 1 , . . . , m j 1 2 d-j times , . . . , m j 2 j , . . . , m j 2 j 2 d-j times ),
where

m j l = 1 2 d-j t∈T j l Z t .
The test is based on the comparison to zero of the squared Euclidean norm Π j Z 2 K equal to

2 d-j 2 j l=1 m j l 2 = 1 2 d-j 2 j l=1   t∈T j l Z t   2 .
Under H 0 , the difference vector Z = I y -I x is a centered Gaussian vector with covariance matrix 2σ 2 Id K where Id K denotes the identity matrix in R K . Hence, under H 0 , Π j Z 2 K /σ 2 follows a χ 2 -distribution with 2 j degrees of freedom. Our test procedure works as follows:

Reject H 0 at level α if for any j = 0 . . . d -1,

Π j Z 2 K /2σ 2 > Ψ -1 2 j (α/d),
where Ψ -1 D denotes the quantile function of a χ 2 (D)-distribution. To simplify the presentation, we have used a Bonferroni correction to ensure that this multiple testing procedure is of level α. Some finer corrections could be applied as proposed by [START_REF] Baraud | Adaptive tests of linear hypotheses by model selection[END_REF] and [START_REF] Durot | An adaptive test for zero mean[END_REF] or, as we do in practice (see Section 5), by using an FDR approach [START_REF] Benjamini | Controlling the False Discovery Rate: a practical and powerful approach to multiple testing[END_REF].

From a clinical point of view, it is clear that the unobservable true enhancements i x and i y are functions of time. So is the function F defined by F (t) := i x (t)i y (t). Moreover F is smooth (regular) in this context. We have f k = F (t k ) and here these multi-tests are adaptive with respect to the unknown Hölder regularity s (see [START_REF] Wasserman | All of nonparametric statistics[END_REF]) of the function F and, for a given fixed power, this test automatically achieves the best rate of testing ρ s (K) for all regularities s > 1/4 (see [33, Thm 1], [34, Thm 4 and 5] and [START_REF] Gayraud | Adaptive minimax testing in the discrete regression scheme[END_REF]).

We now have a tool to compare enhancement sequences at two voxel locations x and y with respect to the known noise level σ.

We write "I y ≡ α σ 2 I x " when the enhancement difference average vector i x -i y is accepted to be the zero vector at level α with respect to a noise level σ 2 following the above construction. This defines precisely the "statistical closeness" used in Section 3.

Growing Time Homogeneous Neighborhood

We now present the construction of the neighborhood V x for a fixed spatial location x.

Given an estimation (or denoising) procedure of the enhancement using a set of locations V , for example the empirical mean defined by

ÎV = 1 |V | y∈V I y (3)
or a generalized median defined, for example, by ÎV = arg min

J∈R K y∈V K i=1 |I y (t i ) -J(t i )|, (4)
given a maximal number of iterations K and a increasing sequence of integers n 1 , . . . , n K with a geometrical growth, we consider in a first step the set

W x = {y ∈ X such that y = x and I y ≡ α σ 2 I x } (5)
of the spatial locations for which enhancements are statistically similar to those of x with respect to the multi-test introduced in Appendix 7.1 and we set V 0 = ∅ and i = 0. Then, setting i = i + 1, we sequentially grow rings, denoted W i , around x and set the neighborhood V i = V i-1 ∪ W i . The ring W i is made of the n i closest points in W x not in V i-1 . At each step, W i is tested for statistical coherence, defined below, with previously built neighborhoods V j , j < i. When the statistical coherence of W i+1 is refused or when i = K, the algorithm stops and returns V x = V i as the selected neighborhood and Îx = ÎVi as the denoised enhancement.

The outline of our construction is given by the flowchart in Figure 2 and Algorithm 1 in Appendix 7.3. This method is illustrated in Figure 10. It shows a zoom on the axial upper abdominal section focused on the right posterior part of the sequence presented in Section 5: a vertebra (1), the aorta (2), the liver (3) and a tumor (4). This figure shows for a specific voxel x (white dot designed by the white arrow) inside the tumor the first four rings W i , i = 1 . . . 4 (in order : red, yellow, green and blue). Due to the pre-selection of voxels in W x , the rings are neither convex nor connected and follow the heterogeneity of the tumor.

During the iterations, the statistical coherence between W i and the previously build neighborhood V j , j < i, is ensured by the test defined by where Ĵi = ÎWi and Îj = ÎVj . This is a generalization of the multi-test introduced in Appendix 7.1 which compares Ĵi , the estimate on the ring W i , with all estimates Î1 , . . . , Îi built on the previous nested neighborhoods V 1 , ..., V i . This multiple testing procedure compares the statistical hypothesis H 0 : "E( Ĵi -Îj ) = 0, for all j = 1 . . . i", against H 1 : "E( Ĵi -Îj ) = 0, for at least one j", where EZ denotes the expectation of a vector Z.

The use of α/i is a classical Bonferroni correction to ensure that the level of this test is α. Other corrections are possible as described in Appendix 7.1. The correction in the noise level, defined by ρ(|V j |, |W i |), aims to take into account the fact that the estimates come from independent samples with respective sizes |V j | and |W i |. This correction depends on the choice of the estimates: for example, if the estimate on a set of locations derives from an empirical mean defined in (3) we set

ρ(|V j |, |W i |) = |V j | -1 + |W i | -1 ,
which derives from the mean equality test for two Gaussian samples with same known variance (see Wald test or Likelihood test in [START_REF] Van Der Vaart | Asymptotic Statistics[END_REF]). The sizes n i of the rings W i are chosen to ensure that the sizes of the V i grow at least geometrically and that the number of tests is of logarithmic order.

Remark 1: The pre-selection step specified by the set W x , which is neither necessarily convex nor connected, allows the denoising construction to obtain "neighborhoods" built from different objects of the same type over long distances. This is useful in medical images where the same kind of tissue may reappear in different areas. Remark 2: This method is closely related to Lepski's method (see e.g. [START_REF] Lepskii 0 | Asymptotically minimax adaptive estimation I: Upper bounds. Optimally adaptive estimates[END_REF] and [START_REF] Lepskii 0 | Optimal pointwise adaptive methods in nonparametric estimation[END_REF]), which method relies on comparing directly Îi+1 with the previously built estimates Î0 . . . Îi . Because they are built on nested subsets, these estimates are not independent. This lack of independence is the drawback of this method from a practical point of view : the growth of neighborhoods often stops too late, leading to over-smoothing [START_REF] Reiß | Pointwise adaptive estimation for robust and quantile regression[END_REF] and in the DCE-CT case to mixing dynamics. This drawback is corrected in our approach by the use of rings. Remark 3: The neighborhood comparison protects from extra bias which could appear in a one step procedure using a subset of W x without more control (see [START_REF] Reiß | Pointwise adaptive estimation for robust and quantile regression[END_REF]). The sizes grow geometrically to ensure that the number of tests remains logarithmic with respect to the number of used voxels, allowing a good control of the power of this multi-test procedure. Remark 4: In the second step of the algorithm, if the number of locations in W x \ V i (set of voxels in W x , but not in V i ) is too small, the algorithm ends and returns the last estimate. Remark 5: The maximal size of V i which controls the gain in the signal to noise ratio is fixed for computational time reasons. Remark 6: In the case of a generalized median defined in (4), to take into account the fact that (i) for Laplace-distributed noise the asymptotic the median differs from that of the mean by a factor 2 in variance (see [START_REF] Ferguson | Asymptotic Joint Distribution of Sample Mean and a Sample Quantile[END_REF] or [42, cor. 21

Algorithms

In this section, we describe synthetically our denoising and clustering algorithms. Given two lists A and B, the notations [A; B], A ∪ B and A \ B denote the concatenation without deletion, the union and the complement of B relative to A, respectively.

Figure 1 :

 1 Figure 1: Histogram of the noise in DCE-CT baseline images with fitted Laplace distribution with parameter λ = 29 (σ = 41HU )

Figure 2 :

 2 Figure 2: Flowchart of the spatially pointwise denoising procedure. See Appendix 7.1 for more informations on the test denoted ≡ α σ 2 and Appendix 7.2 for precisions on ρ. Here, the centroids Îi and Ĵ are the means of the enhancements in their corresponding neighborhoods; other possibilities such as generalized medians are discussed in Appendix 7.2.

Figure 3 :

 3 Figure 3: Axial CT images of the abdomen centered on a liver metastasis. These images belong to a dynamic series. Left column : acquired 25 seconds after the beginning of the acquisition, during the arterial phase. Right column : 45 seconds after the beginning of the acquisition, after the arterial phase. For each column we show the original (top), denoised (middle) and residuals (bottom).

  (a) a class is defined from a single voxel; (b) a new class center is checked to be statistically close from any existing class center; (c) the classes with statistically close centers are gathered.

α σ 2

 2 ρ(1,min(|Vy|,|c|)) Ĵ is accepted. Set L := (L ∪ c) \ c Ĵ to remove selected points from the list. Step b -Class checking : The new class c Ĵ with center Ĵ made at the previous step is checked for closeness with all previously constructed classes. If a class c is found to satisfy Îc ≡ α σ 2 ρ(1,min(|c|,|c Ĵ |)) Ĵ, then c Ĵ and c are merged together to form a new class c and the algorithm returns to Step a/ with this new c.

Figure 4 :

 4 Figure 4: Enhancement curves obtained from manually drawn ROI containing more than 100 voxels within the aorta (top) and the tumor (bottom). For each sub-figure, the mean curve inside the ROI is the thick plain line, the 5% and 95% quantiles of the curves of the selected points are in thick dotted lines. The curves of 5 randomly selected points inside the ROI are drawn as fine lines.

Figure 5

 5 Figure 5-a: Left column : Voxel within the aorta -Right column : Voxel within a small dorsal artery on the right side of the spine. For each column, (Top) the arrow points to the selected voxel x on one axial slice of the sequence; the pink dots show the constructed neighborhood Vx. (Bottom) the thick dotted black line shows the tissue enhancement at voxel x; the thick black line is the denoised dynamic, centroid of the Vx dynamics. In addition, individual curves of 5 randomly selected voxels in Vx are shown in the background as fine line.

Figure 5

 5 Figure 5-b: Left column : Voxel within the tumor -Right column : Voxel within a tumor "hot spot" :The estimated curve shows a quicker and larger enhancement providing a stronger initial slope that shows a perfusion and a more important blood volume than for the voxel within the tumor itself (left column). See Figure5-a for more information.

Figures 5 -

 5 a to 5-b present the result of our method applied to four voxels x in (respectively): (a) the aorta; (b) a small dorsal artery that is located near the spine and crosses our CT-plane on only one or two voxels; (c) the tumor; (d) a tumor hot-spot.Figures 5-a to 5-b show the result of our procedure within a voxel and are divided into two sub-figures.

Figure 6 :

 6 Figure 6: clustering result based on the enhancement curves (variation from the baseline) -(top) clusters in artificial colors ; (bottom) estimated enhancements, "centroids" of the clusters. The numbers link the estimated enhancements to the clusters.

Figure 7 :

 7 Figure 7: clustering result using the entire signal temporal dynamics including the specific baseline level of each tissue -colors are artificial. For tissues showing sufficient enhancement, the clustering is similar to that obtained in Fig.6 (see tumor tissues (2,7,8,9) for example). The other tissues (4) are distinguished here through their baseline.

Figure 8 Figure 8

 88 Figure 8-a: Artificial Simulated Dynamical Image : (bottom) The nine artificial enhancements used in our simulations. The dotted vertical line shows time 30: used for the presentation of the results of denoising and clustering on these artificial data. The numbers link the used enhancements with their locations. (top) The dynamical image shown at time-slice 30.In the upper-right of the sub-figure, the large ball simulates a tumor with behaviors changing continuously from behavior 7 (in the middle) to 4 (on the border). The exact profile along a diameter of this tumor is given in Figure9. The dark homogenous areas are surrounded by a fine white line to emphasize the borders. The gray scale comes from Figure8-c

Figure 8 -Figure 8

 88 Figure 8-d: Denoised version of the artificial data using its full grayscale. Dynamical image shown at time-slice 30

Figure 9 :

 9 Figure9: Representation of the profile inside the simulated tumor (thick line). The abscissa is the radius and the ordinate defines the enhancement mixture: a value 5.2 obtained at radius 0.7 represents a mixture of 20% of enhancement 5 combined with 80% of enhancement 6. The fine dotted line represents a feasible piecewise constant representation of this profile.

  (|Vj |,|Wi|) Îj , j = 1 . . . i -1,

Figure 10 :

 10 Figure 10: First four rings W i (in order : red, yellow, green and blue) selected around the selected voxel x specified by the white dot and the white arrow.

1 + 7 /

 17 .5]); (ii) the median and mean are the same when sets contain only one element, we suggest the approximationρ(|V j |, |W i |) = 1 4 |V j | |V j | + 1 + 7/|W i | |W i | .

Acknowledgment

This work was supported in 2008 by a grant "Bonus Qualité Recherche" from the University Paris Descartes for the project "Cancer Angiogénèse et Outils 30 Mathématiques".

Appendix

"Multiple testing for the comparison of random vectors"

Given two spatial locations x and y, we present the statistical test used to compare to the zero vector the enhancement difference vector with components

For the sake of simplicity, we write Z = f + σε and we introduce these tests in their simplified version: the Gaussian case with σ known. Such a test of comparison to a zero vector is derived from the theoretical work of [START_REF] Spokoiny | Adaptive hypothesis testing using wavelets[END_REF][START_REF] Baraud | Adaptive tests of linear hypotheses by model selection[END_REF] and [START_REF] Durot | An adaptive test for zero mean[END_REF] who consider general frameworks where σ can be unknown and where ε need not necessarily be Gaussian but is at least symmetrical, ensured by the use of differences. We aim to test whether the mean vector f is zero or not and hence consider the hypotheses

Let us emphasize that for a prescribed level, due to their good properties of adaptation, these tests allow good control of the error of second kind (false negative), important in diagnostic settings.

We suppose that K, the time number of the DCE-CT sequence, is of the form K = 2 d . We consider the regular dyadic decomposition of the observation times t 1 . . . t 2 d . For j = 0 . . . d -1, we denote by T j 1 . . . T j 2 j the 2 j intervals with 2 d-j time indices

Denoising algorithm

Input: x ∈ X Output: V x and Îx Initialization : // define first neighborhood i := 0 ; V 0 := ∅ ; W 0 := {x} ; accepted := 1 // Main loop : find new neighbors and test statistical closeness while accepted do // build new neighborhood i

if not( Ĵi ≡ level var Îj ) then accepted := 0 (See Section 7.2 for precisions on ρ and these tests.)

Algorithm 1: Spatially pointwise denoising algorithm

The Algorithm 1 is summarized by the flowchart given by Figure 2.

Clustering algorithm