Tristan Glatard

Carole Lartizien

Bernard Gibaud

Rafael Ferreira Da Silva

Germain Forestier

Frederic Cervenansky

Martino Alessandrini

Hugues Benoit-Cattin

Olivier Bernard

Sorina Camarasu-Pop

Frédéric Cervenansky

Nadia Cerezo

Patrick Clarysse

Alban Gaignard

Patrick Hugonnard

Hervé Liebgott

Simon Marache

Adrien Marion

Johan Montagnat

Joachim Tabary

Denis Friboulet

A virtual

Keywords: Medical image simulation, ontology, workflows, distributed computing infrastructures

I. INTRODUCTION

M Edical images can be simulated from digital models of the human body for a variety of applications in research and industry, including fast prototyping of new devices and the evaluation of image analysis algorithms [START_REF] Gilliam | Echocardiographic simulation for validation of automated segmentation methods[END_REF], [START_REF] Prastawa | Simulation of brain tumors in MR images for evaluation of segmentation efficacy[END_REF], [START_REF] Wagenknecht | Simulation of 3D MRI brain images for quantitative evaluation of image segmentation algorithms[END_REF]. Several image modalities are commonly simulated, among which Magnetic Resonance Imaging (MRI), ultrasound imaging (US), Positron Emission Tomography (PET), and Computed Tomography (CT).

However, image simulation remains mastered only by a few, due to the variety, complexity and heaviness of simulation pipelines. Both the simulation code and the physical model of the imaged object can be very elaborate and specific to the imaging modality. They are usually not designed by the same person or group, and need to be shared. Sharing models requires consistent descriptions so that they can be reused in simulations of different modalities. This is challenging because simulators need different types of information, for instance relaxation times for MR, or radiopharmaceutical activity for PET. Moreover, models contain information not only about anatomy but also about pathologies and external entities present in the body during image acquisition.

Additionally, computing times and the volume of produced data also limit the level of realism and the size of the simulation scene. Although distributed computing infrastructures (DCI) can help in supporting the computing load and storing the generated data, using them should not become an additional burden to simulation users. High-level interfaces should allow their transparent exploitation, together with customized deployments and optimization.

A platform making simulations less ominous, especially for newcomers, is therefore needed. This platform should facilitate access to models and physical parameters, to simulators of different modalities, and to appropriate computing power and storage. This paper describes VIP, an openly-accessible online platform sharing models and simulators, and supporting the execution of heavy simulations. It extends the summary presented in [START_REF] Glatard | Multi-modality image simulation with the virtual imaging platform: illustration on cardiac MRI and echography[END_REF] with a more detailed platform description and additional simulations illustrating different image modalities.

The manuscript is organized as follows. After a review of related work in Section II, we describe in Section III how ontologies were created to describe models in VIP. In Section IV we show how a workflow-based template was used to integrate simulators of 4 different imaging modalities and deploy them on a DCI. Methods and techniques employed to support simulation execution are described in Section V, highlighting computation management, data provenance recording, storage of models and simulated data, and interface. Finally, Section VI exemplifies VIP on PET, US, MR and CT simulations of various models. Usage statistics are also reported.

II. RELATED WORK

Several databases of simulated images are available for whole-body or brain PET imaging [START_REF] Tomei | OncoPET DB: a freely distributed database of realistic simulated whole body [18f]FDG PET images for oncology[END_REF], [START_REF] Reilhac | PET-SORTEO: Validation and Development of Database of Simulated PET Volumes[END_REF], [START_REF] Castiglioni | A publicly accessible Monte-Carlo database for validation purpose in emission tomography[END_REF], [START_REF] Stute | Lucas: Efficient monte carlo simulations of highly realistic pet tumor images[END_REF], and for brain MRI, e.g. Brainweb [START_REF] Kwan | MRI simulation-based evaluation of image-processing and classification methods[END_REF], offering an online simulation service for MRI. The platform in [START_REF] Dikshit | An online interactive simulation system for medical imaging education[END_REF] provides simulationbased interactive tutorials on imaging modalities. Simulators are usually dedicated to a particular modality, for instance SPECT [START_REF] Beekman | Fast SPECT simulation including object shape dependant scatter[END_REF], CT [START_REF] Lazos | An integrated research tool for X-ray imaging simulation[END_REF], mammography [START_REF] Peplow | Digital mammography image simulation using Monte Carlo[END_REF], MRI [START_REF] Drobnjak | Development of a functional magnetic resonance imaging simulator for modeling realistic rigid-body motion artifacts[END_REF], or US [START_REF] Wein | Simulation and fully automatic multimodal registration of medical ultrasound[END_REF], [START_REF] Aiger | Real-time ultrasound imaging simulation[END_REF], [START_REF] Shams | Real-time simulation of medical ultrasound from CT images[END_REF], with the notable exception of GATE [START_REF] Jan | Gate v6: a major enhancement of the gate simulation platform enabling modelling of ct and radiotherapy[END_REF], which was validated for PET, SPECT, CT, and even radiation therapy. VIP targets a multi-modality online simulation service relying on existing codes.

The sharing of object models is important for initiatives related to the Virtual Physiological Human (VPH) program. FieldML [START_REF] Christie | FieldML: concepts and implementation[END_REF] was developed for this program as a standard to represent information contained in models. It can structure repositories and databases such as those available in eu-Heart [START_REF] Smith | euHeart: Personalized and integrated cardiac care using patient-specific cardiovascular modelling[END_REF], focusing on the modeling of the heart. Several ontologies are available to describe human anatomy [START_REF] Rosse | The foundational model of anatomy ontology[END_REF], pathology [START_REF] Schofield | Phenotype ontologies for mouse and man: bridging the semantic gap[END_REF], qualities [START_REF] Gkoutos | Using ontologies to describe mouse phenotypes[END_REF], and radiological terms [START_REF] Langlotz | Radlex: a new method for indexing online educational materials[END_REF], but none can completely describe the exhaustive list of physical and biological parameters targeted for VIP models. A specific development is therefore needed.

One approach to medical image simulation is to simulate images from geometrical models parametrized with distributions of physical properties such as magnetic properties, echogenicity, radioactivity, or chemical composition. This approach is flexible since geometrical models can be deformed to a particular individual [START_REF] Tobon-Gomez | Realistic simulation of cardiac magnetic resonance studies modeling anatomical variability, trabeculae, and papillary muscles[END_REF], [START_REF] Lemaitre | Incorporating patient specific variability in the simulation of realistic whole body 18F-FDG distributions for oncology applications[END_REF]. However, obtaining realistic simulations is challenging because parameter distributions are not easily estimated and human morphology is difficult to model from geometrical objects alone. An alternate [START_REF] Kutter | Visualization and GPU-accelerated simulation of medical ultrasound from CT images[END_REF] is to create a digital object from a real acquisition by extracting object geometry and physical parameters from signal properties. Resulting images usually look very realistic, but this method requires an in-vivo acquisition for each simulation, which constrains investigation. Both model-based and imagebased approaches are targeted in VIP.

Some online platforms, commonly called scientific gateways, integrate medical imaging software with access to computing and storage resources. neuGRID [START_REF] Frisoni | Virtual imaging laboratories for marker discovery in neurodegenerative diseases[END_REF], CBRAIN1 , and LONI pipeline [START_REF] Rex | The LONI pipeline processing environment[END_REF] (Laboratory Of Neuro Imaging) target neuroimaging data sharing and analysis using DCIs. In neuGRID, users can start remote desktop sessions on machines where image analysis tools and clients are pre-installed to access resources of the European Grid Initiative (EGI) 2 . CBRAIN is a web portal where neuroscientists can launch data processings such as segmentation tools on several clusters. The location of files and executions is controlled from the interface. In the LONI pipeline, users describe their own image processing pipelines executed on local resources. Neurolog [START_REF] Montagnat | NeuroLOG: a community-driven middleware design[END_REF] focuses on the sharing and reuse of heterogeneous data and tools produced by distributed sites. The OntoNeuroLOG ontology describing involved datasets and entities is used to harmonize database schemas, and annotate tools. e-bioinfra [START_REF] Shahand | Front-ends to Biomedical Data Analysis on Grids[END_REF] is a web portal for executing image processing tools on distributed computing resources, e.g. Freesurfer and the Functional MRI of the Brain Software Library (FSL3 ). Some scientific gateways such as EUMEDGrid4 , DECIDE 5 , and GISELA [START_REF] Barbera | Supporting e-science applications on e-infrastructures: Some use cases from latin america[END_REF] are built from a set of reusable portlets providing basic functions such as authentication and data transfer.

Compared with the reviewed alternatives, VIP is distinguished by (i) an interface dedicated to medical image simulation, including simulators of 4 image modalities and a repository to store physical and biological models, (ii) a workflow-based methodology that allows fast porting of new simulators with minor adaptations, (iii) a web interface that totally relieves users from resource management decisions.

To objectively compare the activity of the described platforms, the consumed cumulative CPU time metric gives an indication on their usage, therefore usability. In 2011, the ac-counting system of the EGI6 reported that neuGRID consumed 10.3 CPU years, EUMEDGrid consumed 3 CPU years, and ebioinfra consumed 46.3 CPU years. The consumption of VIP is reported in Section VI.

III. MODEL INTEGRATION

We chose ontologies to share model semantics in the platform. Apart from structuring the model repository, two additional roles of ontologies are interesting [START_REF] Gruber | Toward principles for the design of ontologies used for knowledge sharing[END_REF]: (1) they provide a standard vocabulary to refer to conceptual entities, and (2) they model conceptual entities formally, using axioms expressed in a logical language. In VIP, model files describing geometrical objects (represented as voxel maps or meshes) are annotated with concepts defined in our application ontology called OntoVIP (described below). Annotated models are stored in the repository shown on Fig. 1.

A. Ontology design

OntoVIP integrates components describing representational objects (model layers, model layer parts, etc.) and their associated real-world entities (anatomical/pathological structures, foreign bodies, etc.). To increase the potential sharing and interoperability of our semantic model with other systems, we reused existing ontologies or extracted relevant parts from them. Anatomical terms were extracted from the Foundational Model of Anatomy (FMA) [START_REF] Rosse | The foundational model of anatomy ontology[END_REF], pathological structures were extracted from the Mammalian pathology ontology (MPATH [START_REF] Schofield | Phenotype ontologies for mouse and man: bridging the semantic gap[END_REF]), qualities from the phenotype and trait ontology (PATO [START_REF] Gkoutos | Using ontologies to describe mouse phenotypes[END_REF]), and foreign bodies and external agents from RadLex [START_REF] Langlotz | Radlex: a new method for indexing online educational materials[END_REF]. Although extractions are easily performed to extend the list of terms, they have to be performed in moderation to avoid slowing down the reasoning engine by confusing annotation with too many available entities.

The above components were imported as modules into OntoVIP following the integration approach used in [START_REF] Temal | Towards an ontology for sharing medical images and regions of interest in neuroimaging[END_REF]. This integration relies on the foundational ontology of Descriptive Ontology for Linguistic and Cognitive Engineering (DOLCE [START_REF] Masolo | WonderWeb Deliverable D18. The WonderWeb Library of Foundational Ontologies and the DOLCE ontology[END_REF]) and on other core ontologies guaranteeing overall consistency. Models are described using two basic entities: (1) values-layers describe data files associating each point of a 3D map with a particular value, e.g. the activity of a radiopharmaceutical; (2) object-layers describe data files containing labels associated with object-layer-parts linked to anatomical, pathological, foreign-body or externalagent objects. Each object can be associated with a mathematical-distribution, to specify the range of an associated physical parameter, e.g., T1 values. Layers can be associated to different time references in order to model time-varying phenomena, e.g., movement of an organ, metabolism of a substance, or growth of a tumor. Two different time scales are considered: (1) instants denote variations during an image acquisition and (2) time-points denote changes between different imaging procedures. A detailed description of this ontology, specifying the various representational entities and how they refer to real-world entities is provided in [START_REF] Forestier | Sharing object models for multi-modality medical image simulation: a semantic approach[END_REF].

As discussed above, the motivation for ontological modeling is two-fold. First, it provides an organized vocabulary to describe information. This is useful to query the repository using search criteria corresponding to model categories (e.g. static vs dynamic) or model content, such as anatomical structures or contrast agents. Queries may use terms extracted from a relational database, but the ontology brings flexibility with the desired specificity. For example, a lesion may be annotated using the generic term neoplasm or the more specific term glioma. Thanks to relations in the ontology, searching for models containing neoplasm returns models annotated with neoplasm and glioma. Other inferences could be made from part-whole relationships, such as FMA's constitutionalPart / constitutionalPartOf, to enable selecting entities that are part of a structure of interest.

The second motivation to ontological modeling is their reasoning capabilities. For instance, rules were implemented to check the compatibility of models with modalities and simulators (see top bar on Fig. 1). A model is compatible with a simulator if it has all the required physical parameters. These rules are important to prevent users from launching unsuccessful simulations on incomplete models.

B. Semantic model repository

Models can be uploaded to VIP in the form of data files bundled with OntoVIP annotations expressed as Resource Description Framework (RDF7 ) triplets. For instance, in the mymodel rdf:type ontovip:static-model triplet, ontovip:static-model denotes the ontology class modeling static models. During model import, the consistency of assertions is checked before being added to the repository.

A standalone software currently assists users with model annotation. Knowledge stored in OntoVIP is leveraged to enrich annotations with limited manual effort from the user. For example, when a user annotates a data file with the term Brain, the software automatically detects that it is an anatomical-object using subsumption information in OntoVIP. This annotation software is progressively integrated in the online platform.

Models uploaded to the repository can be browsed according to their model parts, time information (static or dynamic), and presence of pathology, foreign body or external agent. Beyond that, rules are applied to check if a model can be used in a simulation of a particular modality, i.e., all the required physical parameters are defined. This is one of the added values of the ontology-based approach for manipulating simulation models. Fig. 1 shows a cardiac model in the model repository. Models stored in the repository can also be visualized in a 3D interface where simulation scenes can be defined, as seen on Fig. 2. VIP model repository showing a cardiac model ready for MRI simulation, where physical parameters of other modalities are missing (for readability purposes, the object list in the anatomical layer was truncated). Fig. 2.

Screenshot of the simulation scene interface. A thorax model was selected from the semantic model repository and an US transducer is positioned to simulate an apical echocardiography (black object on the right with sectorial imaging plane in gray).

IV. SIMULATION WORKFLOWS

VIP integrates simulators based on their workflow description, and without modifying their code so that their validation remains under the responsibility of their developers. The workflow representation was chosen because it provides structured, graph-based representations of applications, which is useful for several reasons. Firstly, graph representations are parallel languages, which facilitates workflow deployment on distributed computing infrastructures. Workflows also facilitate data annotation because dependencies between data files, parameters, programs and results (a.k.a data provenance) are clearly expressed, which structures output data. Finally, workflows foster reusability of software components among simulators. Ultimately, it is envisioned that the platform could facilitate simulator integration by suggesting workflow components based on existing tools (e.g. format converters).

Integrating image simulators in VIP with minimal software refactoring is difficult due to the variety of their interfaces, parameters and characteristics related to (i) the format and in- formation required in the models; (ii) the acquisition protocols and imaging sequence and (iii) the computing model and its parallelization. To facilitate simulator integration, we proposed in [START_REF] Marion | Multi-modality medical image simulation of biological models with the Virtual Imaging Platform (VIP)[END_REF] a workflow template composed of the components presented in Fig. 3. Object preparation converts a repositoryannotated model into the native format suitable for the simulator. It combines the different model layers, defines the scene geometry from the geometrical transformations produced from the 3D interface, splits dynamic models into multi-static models, and adjusts physical parameters if needed. Parameter generation generates simulator parameter files from parameter values. Depending on the type of users, some parameters may be pre-set and hidden. Core simulation invokes the simulation code. It splits the simulation into independent chunks of data on which the code is iterated. Finally, post-processing performs image reconstruction and format conversion.

To facilitate the integration of simulators, the Conceptual Workflows formalism was developed [START_REF] Cerezo | Scientific workflows reuse through conceptual workflows[END_REF]. Conceptual Workflows represent a family of workflows complying with a given template as illustrated on Fig. 3 for VIP. This formalism not only provides a frame for the design of new workflows, but it also assists users in workflow creation and validation. Following this formalism, image simulators of 4 modalities were integrated using the Gwendia (Grid Workflow Efficient Enactment for Data Intensive Applications) workflow language [START_REF] Montagnat | A data-driven workflow language for grids based on array programming principles[END_REF]: SIMRI [START_REF] Benoit-Cattin | The SIMRI project : a versatile and interactive MRI simulator[END_REF] for MRI, FIELD-II [START_REF] Jensen | Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers[END_REF] for US imaging, SORTEO [START_REF] Reilhac | PET-SORTEO: Validation and Development of Database of Simulated PET Volumes[END_REF] for PET, and SINDBAD [START_REF] Tabary | Realistic X-Ray CT simulation of the XCAT phantom with SINDBAD[END_REF] for CT.

Object preparation for SIMRI applies the scene transformation (produced by the interface on Fig. 2), format conversion, layer combination, and time splitting to the model. Core simulation is based on the concurrent computation of the magnetization of spin vectors and uses the Message Passing Interface (MPI 8 ) since the simulator was already parallelized using this framework. Image reconstruction by inverse Fourier transform of the k-space is included in SIMRI.

Object preparation for SORTEO prepares the model for both emission and transmission simulations. Radioactivity is extracted from the model description and added to the description of the acquisition protocol. Photo-electric and Compton cross-sections of the materials are obtained from specific look-up tables. For core simulation, SORTEO uses a two-step Monte-Carlo simulation. First, single events are simulated to evaluate the probability of a photon to be detected 8 http://www.mcs.anl.gov/research/projects/mpi on a detector element when emitted in a region of the model. Second, the simulator generates the true, scattered and random coincidence events using a deadtime model which incorporates these singles rates. As Monte-Carlo simulations, these steps can be split in any number of jobs lower than the number of independent random events. SORTEO produces sinograms or list-mode data. Although data correction and reconstruction are not included in VIP, the documentation refers to the STIR software 9 for both analytic and iterative reconstructions. Reconstruction examples will also be available in the portal.

Object preparation for FIELD-II either samples scatterer positions and amplitudes from distributions of the model, or applies the scene transformation to existing scatterers. Core simulation concurrently simulates the radio-frequency (RF) lines involved in the simulation. It supports both 1D (linear and sectorial) and 2D transducers. Once all lines are simulated, an RF matrix is assembled and the final image is reconstructed using envelope detection and Cartesian reconstruction.

Object preparation for SINDBAD generates cross-sections from material properties (chemical composition), and writes the scene transformation in SINDBAD format. Core simulation was implemented as a two-level splitting: each projection of the 3D scan is computed concurrently, and in addition, Monte-Carlo simulations in each projection are split depending on the number of involved photons. The reconstruction of 3D images from projections is not provided, but the documentation refers to the image reconstruction toolbox of Fessler10 that can be used for that purpose.

Different types of scanners can be used in VIP, depending on the capabilities of the simulators. For SORTEO and SINDBAD, scanners can be defined in parameter files; for FIELD-II, the definition of the transducer is one of the simulation parameters; for SIMRI, different values of B0 and antenna profiles can be set. VIP is also extensible to other simulators, for instance to include another modality, or to compare simulators of a same modality. A new simulator can be integrated as a Gwendia workflow, following the template on Fig. 3. Due to technical constraints, workflow development still requires a developer's assistance. More image processing tools can also be integrated in VIP. Around 10 tools are currently available, including the GATE simulator, and Freesurfer for neuroimaging 11 .

V. SIMULATION EXECUTION

A. Deployment on distributed computing resources

Simulation parallelization relies on data parallelism to avoid refactoring simulator codes. Simulators are iterated on subsets of the simulation scene or of input parameter sets, and the resulting partial results are eventually merged by another process. Simulations are launched on the biomed virtual organization of the European Grid Infrastructure (EGI) 12 . EGI has access to more than 100 computing clusters worldwide, however these are not dedicated and performance can be significantly hampered by the activity of other users, or resource downtimes. To address this, VIP accesses resources by means of the Distributed Infrastructure with Remote Agent Control (DIRAC [START_REF] Tsaregorodtsev | DIRAC3 . The New Generation of the LHCb Grid Software[END_REF]) pilot-job system: computing jobs are not directly queued to the infrastructure, but they are kept on the VIP server until agents actually reach resources and fetch them. This technique is well used in large distributed systems as it improves fault tolerance and overall performance.

As described in [START_REF] Ferreira Da Silva | Multi-Infrastructure Workflow Execution for Medical Simulation in the Virtual Imaging Platform[END_REF], the pilot-job system was extended so that VIP users can also use their own clusters to complement grid resources. Results in [START_REF] Ferreira Da Silva | Multi-Infrastructure Workflow Execution for Medical Simulation in the Virtual Imaging Platform[END_REF] show that introducing even a small fraction of dedicated computing resources has an important impact on the performance of the simulation.

B. Data provenance

Simulations produce large amounts of diverse data, difficult to manipulate by non-expert users. To address this, provenance metadata is recorded during execution.

Provenance is registered on-the-fly in an OPM-compliant RDF repository (Open Provenance Model [START_REF] Moreau | The Open Provenance Model: an overview[END_REF]) tracking simulation component invocations, and the related consumed and produced data. The relations among files and parameters produced during a simulation are recorded as provenance graphs.

Such a provenance model enables inference of new statements about the simulated data. For instance, rules can be applied to provenance paths to produce annotations describing simulation results and build a shared, structured database of simulated data from user activity.

The scalability of semantic repositories is challenging since they contain information about all the data generated by the platform. For instance, a single execution of the SORTEO core simulation workflow produces 15,000 triplets. To tackle this issue, we distinguish provenance metadata from the annotations of simulated data. Provenance metadata is generic, technical, and volatile whereas simulated data annotations represent only a few statements with a long-term added-value. Consequently, a dedicated short-term semantic repository is deployed to store provenance metadata of each simulation execution. Once generated, simulated data annotations are pushed to a long-term semantic repository. If necessary, the short-term provenance repository can then be cleaned-up.

C. Storage

VIP stores files on the EGI data management system. Files are distributed on storage resources at the different sites supporting the biomed virtual organization. These resources expose a homogeneous interface through which files can be listed and transferred. A central logical file catalog provides a common indexing space for these files. 3.5 PB are available for storage, among which 2 PB are currently used.

Although storage availability is not ensured (availability usually ranges from 80% to 95%), replication is possible by linking several physical files to the same logical name. Availability of critical files (e.g. application workflows) is thus ensured. The platform integrates its own local data manager, where critical files are replicated and accessed when EGI storage is down. Files are also cached by jobs on the execution nodes to further decrease transfer error rate.

D. Interface

VIP is available as a web portal where users can access the model repository (browsing and importing), define a simulation scene through the 3D interface, launch new simulations, monitor performance, and transfer input/output files. Documentation and a messaging system to contact the support team are available in the portal 13 .

Performing a simulation consists of (i) selecting a model from the repository (see Fig. 1), (ii) defining the simulation scene (see Fig. 2), (iii) adding simulation parameters, (iv) launching the simulation, monitoring jobs, and downloading results from file transfer interface. In case the user is already familiar with the simulator parameters, then the first three steps are reduced to only parametrizing the core simulator.

Authentication is done with login and password. Users are mapped to a robot certificate used for all grid authentications. For security purposes, the portal keeps track of all user operations, in particular file transfers and simulation execution. The VIP open access policy fosters experimentation: accounts are created based on valid email address only. Different user levels are distinguished: beginners can launch only one simulation at a time, and cannot write in shared folders. Advanced users have extended rights, but they must briefly describe their activities, and have a grid certificate registered in a virtual organization. To avoid technical problems, the certificate is not used to log in to the portal but only for administrators to check user identity. Users are organized in groups defining access rights to applications and data.

A specific interface component enables file transfer between local user machines and grid storage. The upload process consists in (i) uploading the file from the local machine to the portal and (ii) transferring the file to the grid through an asynchronous pool of transfers processed sequentially. Download is performed similarly, in the opposite direction. This two-step process avoids connectivity issues between user machines and distributed grid hosts. The transfer pool manages the load of concurrent transfers performed by the server to avoid network clogging on the portal machine. It also replicates files to ensure availability, and it caches a local copy of critical files such as application workflows. An optimized file browser is available to interact with the transfer service. It caches grid directory content and uses native catalog commands to ensure fast browsing. File permissions are also enforced by the browser: users have a private folder readable only by themselves, and each group has a shared folder.

VI. RESULTS

A. Integrated Models

A heart-thorax model (ADAM) and a full-body model (Zubal) are currently integrated. ADAM [START_REF] Haddad | Un modéle anthropomorphique et dynamique du thorax respirant et du coeur battant[END_REF] contains geometrical descriptions for lungs, thorax, aorta, myocardium, spine, atria, and ventricles, with associated physical parameters for T1, T2, T2*, proton density, echogenicity, and radioactivity. It also has a tumor and a needle. The Zubal phantom (full body with arms) 14 contains 126 labels with an emphasis on brain structures. Simulations can also be conducted with external models (e.g. the 4D NURBS-based Cardiac-Torso, XCAT [START_REF] Segars | 4D XCAT phantom for multimodality imaging research[END_REF]), but the 3D simulation scene interface is then unavailable and the model has to be manually converted to the file formats required by the simulators.

B. Image Simulations

Whole-body and cardiac images simulated with VIP in 4 modalities are described below along with the acquisition parameters, models and resulting images. Reported cumulative CPU times give an estimation of the simulation time on a single CPU core. Simulation times obtained with VIP are only indicative since they depend on the parallelization parameters and uncontrolled grid load during simulation execution.

1) Whole-body imaging: a whole-body PET/CT acquisition was simulated using SORTEO for PET and SINDBAD for Xray CT scan. The anatomical and activity distribution models were constructed from the XCAT phantom which was fitted to a specific patient anatomy. In addition to the main organs, XCAT includes the vessels and the airway tree which are necessary for high resolution simulations such as CT.

The PET acquisition protocol consisted of a static 224 s acquisition described in [START_REF] Marache | Incorporating patient-specific variability in the OncoPET DB database[END_REF]. The scanner geometry was that of the PET/CT Philips GEMINI system (Philips Healthcare, Cleveland, OH). The 18F-FDG (fluorodeoxyglucose) activity distribution was obtained from [START_REF] Tomei | OncoPET DB: a freely distributed database of realistic simulated whole body [18f]FDG PET images for oncology[END_REF]. Three bed fields were used to cover thorax and abdomen. Each bed field had the dimension of the axial field of view (FoV) of the scanner, i.e, 47 cm. The overlap between beds was set to 50% of the axial FoV to compensate for the sensitivity loss on both extremities of the axial FoV. List-mode data were reconstructed with onepass list-mode expectation minimization [START_REF] Reader | One-pass list-mode EM algorithm for high-resolution 3-D PET image reconstruction into large arrays[END_REF] using 5 iterations and 8 subsets, resulting in 128x128x81 images with 4 mm isotropic voxels.

The CT simulation was exclusively analytical although SINDBAD also has a Monte-Carlo mode. The scanner model was a simplified version of a Philips Scanner (Philips MX8000) which is the CT component of the PET/CT GEMINI scanner. We assumed a point source located at 0.75 m from the center of the thorax (also considered as the rotation axis of the CT system) and at 1.2 m from the plane detector. We used a cone-beam geometry and a standard X-ray energy spectrum with a tube voltage of 110 kVp and an aluminium filter of 2.5 mm. The typical values for intensities and scanner rotation length corresponded to a very low-dose acquisition (resp. 1 mA and 0.5 s). According to the magnification factor scale in this scanner geometry (from 1.2 to 2.4 for the phantom), the voxel size of the CT model was chosen to be half the pixel size of the detector. Precisely, the XCAT phantom was converted into a 500x500x689 volume with an isotropic voxel size of 1 mm. A set of 480 projections were reconstructed with a Feldkamp cone-beam reconstruction algorithm resulting in a 450x450x450 matrix with an isotropic pixel size of 1mm. 4 shows the XCAT model and coronal slices extracted from the CT and PET 18F-FDG simulations. Validation of the first and second order simulated PET image statistics against clinical ones has been performed in [START_REF] Tomei | OncoPET DB: a freely distributed database of realistic simulated whole body [18f]FDG PET images for oncology[END_REF]. For CT, the airway tree of the XCAT model was well rendered in the simulation, in spite of the low resolution of the simplified scanner model (2x2 mm detector size) and the simple analytical model. Adding the noise from the scattered radiation and non-ideal detector would allow to reproduce accurate count rates and noise texture in the reconstructed image. Such validation is out of scope here. The CT simulation represented a cumulative CPU time of 12.8 hours, computed in 1.9 hours on VIP.

2) Heart imaging: A 2D MR balanced steady state free precession (bSSFP) sequence at 1.5 T was simulated on a cardiac cycle (14 instants) extracted from the ADAM model. Scan parameters similar to the simulation conducted by Tobon et al. in [START_REF] Tobon-Gomez | Realistic simulation of cardiac magnetic resonance studies modeling anatomical variability, trabeculae, and papillary muscles[END_REF] were chosen: repetition/echo time = 2.9/1.2 ms and flip angle = 45 • , matrix size 256x256, number of signal average: 1. Magnetic parameters for blood, lung and myocardium were also extracted from [START_REF] Tobon-Gomez | Realistic simulation of cardiac magnetic resonance studies modeling anatomical variability, trabeculae, and papillary muscles[END_REF] while those of fat, muscle, spine, and spinal cord were obtained from [START_REF] Schmitt | Inversion recovery TrueFISP: quantification of T1, T2, and spin density[END_REF], [START_REF] Stanisz | T1, T2 relaxation and magnetization transfer in tissue at 3T[END_REF] and [START_REF] Melhem | MR of the spine with a fast T1-weighted fluid-attenuated inversion recovery sequence[END_REF]. Fig. 5 (top-left) shows one instant of the simulated sequence. The whole sequence can be seen from the VIP gallery 15 . Despite the limited number of anatomical structures compared to real images, contrasts in the blood pool and left-ventricle are coherent with the image obtained in [START_REF] Tobon-Gomez | Realistic simulation of cardiac magnetic resonance studies modeling anatomical variability, trabeculae, and papillary muscles[END_REF] (see Fig. 5, bottomleft), with a blood hypersignal in the cavities as opposed to myocardium. More information about the validation of SIMRI is found in [START_REF] Benoit-Cattin | The SIMRI project : a versatile and interactive MRI simulator[END_REF]. Each instant represented a cumulative CPU time of 31±4 min, computed in 5±2 min on VIP.

A healthy and a pathological FDG-PET acquisition were simulated with ADAM using SORTEO. The scanner geometry of ECAT EXACT HR+ (CTI/Siemens Knoxville) was used, and 18F-FDG was the radiotracer to study glucose metabolism. The simulated activities were respectively 45.1 MBq for the healthy case and 43.8 MBq for the pathological case as used in [START_REF] Haddad | Un modéle anthropomorphique et dynamique du thorax respirant et du coeur battant[END_REF]. The raw data was corrected for attenuation and reconstructed using a standard 3D filtered back-projection A 2D+t echocardiographic sequence was also simulated with an image-based approach. It was obtained by deforming a scattering map with a known motion model. The acquisition was simulated with FIELD-II using a 64-element sectorial probe at 3.75 MHz. A sampling frequency of 40 MHz was set. The view angle was 66 • and the pitch was set to half of the wavelength to avoid grating lobes effects. One frame of the obtained simulated sequence is shown on Fig. 5 (top-right) with the corresponding real image (bottom-right). For a better dynamic perception, the full sequence can be found online 16along with the associated benchmark velocity field. As detailed in [START_REF] Alessandrini | Highly realistic simulation of echocardiographic image sequences for groundtruth validation of motion estimation[END_REF], both visual appearance and motion are very realistic due to the use of a real sequence as a template. As the true motion for the synthetic sequences is known, this sequence can be used as benchmark for the evaluation of myocardium motion estimation algorithms. The cumulative CPU time of the simulation was 42 hours, obtained in 4 hours on VIP.

C. Platform usage

To date, 230 users are registered in VIP, among which 64 logged in during the last month. In 2011, the average monthly CPU consumption was 2.75 years and the yearly cumulative was 33 years. This level of activity is comparable to the main grid platforms integrating medical imaging applications (see values reported in Section II). VIP is significantly used, which is an indicator of its availability and usability. VII. CONCLUSION VIP is a versatile, open-access platform for multi-modality medical image simulation. This article provides a complete overview of VIP and describes design issues and solutions selected to enable sharing of object models, the integration of simulators, and the execution of simulations on distributed computing resources. Results show that VIP can run simulations of 4 modalities and different organs. Usage statistics show that it has a good potential to go beyond a proof-ofconcept level. The development of VIP started in 2009 and mobilized 3 developers for 2 years. Operations require about 2 hours daily. Some limitations remain, partly due to the adoption of a web infrastructure. A centralized web portal simplifies usability and user support, but it is a network bottleneck when large files are transferred. To address this, files could be directly transferred from user hosts to storage sites. Performance of small simulations is another issue due to the sharing of computing resources among several users. Although pilot jobs are used to reduce latency, jobs still have to queue from 1 minute to 1 hour. A few dedicated resources could be used to address this problem. The model repository also has to be exemplified on more models. For instance, some terms will probably have to be added to the ontology when specific medical applications are targeted. Finally, simulators and tools still cannot be integrated directly from the web interface because it requires substantial knowledge about the infrastructure.

VIP operates as an image simulation and image processing service for the academic community. Maintenance costs remain reasonable thanks to the use of external, open-source software, and the exploitation of public distributed computing and storage resources.

Fig. 1 .

 1 Fig. 1.VIP model repository showing a cardiac model ready for MRI simulation, where physical parameters of other modalities are missing (for readability purposes, the object list in the anatomical layer was truncated).

Fig. 3 .

 3 Fig. 3. VIP template represented with the Conceptual Workflow formalism.

Fig. 4 .

 4 Fig. 4. Whole-body FDG-PET and CT simulations (coronal slices). From left to right and top to bottom: simulated 0.5 s CT acquisition, simulated static 224 s FDG-PET acquisition, overlaid CT and FDG-PET simulations, XCAT-based voxel model.

Fig.

  Fig.4shows the XCAT model and coronal slices extracted from the CT and PET 18F-FDG simulations. Validation of the first and second order simulated PET image statistics against clinical ones has been performed in[START_REF] Tomei | OncoPET DB: a freely distributed database of realistic simulated whole body [18f]FDG PET images for oncology[END_REF]. For CT, the airway tree of the XCAT model was well rendered in the simulation, in spite of the low resolution of the simplified scanner model (2x2 mm detector size) and the simple analytical model. Adding the noise from the scattered radiation and non-ideal detector would allow to reproduce accurate count rates and noise texture in the reconstructed image. Such validation is out of scope here. The CT simulation represented a cumulative CPU time of 12.8 hours, computed in 1.9 hours on VIP.2) Heart imaging: A 2D MR balanced steady state free precession (bSSFP) sequence at 1.5 T was simulated on a cardiac cycle (14 instants) extracted from the ADAM model. Scan parameters similar to the simulation conducted by Tobon et al. in[START_REF] Tobon-Gomez | Realistic simulation of cardiac magnetic resonance studies modeling anatomical variability, trabeculae, and papillary muscles[END_REF] were chosen: repetition/echo time = 2.9/1.2 ms and flip angle = 45 • , matrix size 256x256, number of signal average: 1. Magnetic parameters for blood, lung and myocardium were also extracted from[START_REF] Tobon-Gomez | Realistic simulation of cardiac magnetic resonance studies modeling anatomical variability, trabeculae, and papillary muscles[END_REF] while those of fat, muscle, spine, and spinal cord were obtained from[START_REF] Schmitt | Inversion recovery TrueFISP: quantification of T1, T2, and spin density[END_REF],[START_REF] Stanisz | T1, T2 relaxation and magnetization transfer in tissue at 3T[END_REF] and[START_REF] Melhem | MR of the spine with a fast T1-weighted fluid-attenuated inversion recovery sequence[END_REF]. Fig.5(top-left) shows one instant of the simulated sequence. The whole sequence can be seen from the VIP gallery15 . Despite the limited number of anatomical structures compared to real images, contrasts in the blood pool and left-ventricle are coherent with the image obtained in[START_REF] Tobon-Gomez | Realistic simulation of cardiac magnetic resonance studies modeling anatomical variability, trabeculae, and papillary muscles[END_REF] (see Fig.5, bottomleft), with a blood hypersignal in the cavities as opposed to myocardium. More information about the validation of SIMRI is found in[START_REF] Benoit-Cattin | The SIMRI project : a versatile and interactive MRI simulator[END_REF]. Each instant represented a cumulative CPU time of 31±4 min, computed in 5±2 min on VIP.A healthy and a pathological FDG-PET acquisition were simulated with ADAM using SORTEO. The scanner geometry of ECAT EXACT HR+ (CTI/Siemens Knoxville) was used, and 18F-FDG was the radiotracer to study glucose metabolism. The simulated activities were respectively 45.1 MBq for the healthy case and 43.8 MBq for the pathological case as used in[START_REF] Haddad | Un modéle anthropomorphique et dynamique du thorax respirant et du coeur battant[END_REF]. The raw data was corrected for attenuation and reconstructed using a standard 3D filtered back-projection

Fig. 5 .

 5 Fig. 5. Left-top: cardiac short-axis MRI simulated with VIP from ADAM. Left-bottom: simulation obtained by Tobon et al. in [25] with XCAT. Center: simulated cardiac FDG-PET at end diastole, transverse view (top: healthy; bottom: pathological). Right: simulated (top) and real (bottom) images from a 2D+t echocardiography.

http://cbrain.mcgill.ca

http://www.egi.eu

http://www.fmrib.ox.ac.uk/fsl

http://applications.eumedgrid.eu

http://applications.eu-decide.eu

http://accounting.egi.eu

http://www.w3.org/RDF

http://stir.sourceforge.net

http://web.eecs.umich.edu/ ∼ fessler/code

http://surfer.nmr.mgh.harvard.edu

http://www.egi.eu

http://vip.creatis.insa-lyon.fr/documentation

http://vip.creatis.insa-lyon.fr/gallery

http://www.creatis.insa-lyon.fr/us-tagging/news

ACKNOWLEDGMENT

This work is funded by the French National Research Agency (ANR) under grant ANR-09-COSI-03. We thank the European Grid Initiative and "France-Grilles".


Frédéric Havet
email: frederic.havet@inria.fr.

A Karolinna Maia

Min-Li Yu

Complexity of greedy edge-colouring

Keywords: Edge colouring, greedy colouring, greedy algorithm, line graphs, caterpillars, NPcomplete Coloration d'arêtes, coloration gloutonne, algorithme glouton, graphe des lignes, chenilles, NP-complet

The Grundy index of a graph G = (V, E) is the greatest number of colours that the greedy edge-colouring algorithm can use on G. We prove that the problem of determining the Grundy index of a graph G = (V, E) is NP-hard for general graphs. We also show that this problem is polynomial-time solvable for caterpillars. More specifically, we prove that the Grundy index of a caterpillar is ∆(G) or ∆(G) + 1 and present a polynomial-time algorithm to determine it exactly.

Introduction

All the graphs considered in this paper are loopless.

A k-colouring of a graph G = (V, E) is a surjective mapping c : V → {1, 2, . . . , k} such that c(u) = c(v) for any edge uv ∈ E. The chromatic number is χ(G) = min{k | G admits a k-colouring}. On the algorithmic point of view, finding the chromatic number of a graph is a hard problem. For all k ≥ 3 it is NP-complete to decide if a graph admits a k-colouring (see [START_REF] Garey | Computers and intractability. A guide to the theory of NP-completeness[END_REF]). Furthermore, it is NP-hard to approximate the chromatic number within |V (G)| ε0 for some positive constant ε 0 as shown by Lund and Yannakakis [START_REF] Lund | On the hardness of approximating minimization problems[END_REF].

Hence lots of heuristics have been developed to colour a graph. The most basic and widespread because it works on-line is the greedy algorithm. Given a vertex ordering σ = v 1 < • • • < v n of V (G), this algorithm colours the vertices in the order v 1 , . . . , v n , assigning to v i the smallest positive integer not used on its lower-indexed neighbours. A colouring resulting of the greedy algorithm is called a greedy colouring. The Grundy number Γ(G) is the largest k such that G has a greedy k-colouring. Easily, χ(G) ≤ Γ(G) ≤ ∆(G) + 1.

Zaker [START_REF] Zaker | Results on the Grundy chromatic number of graphs[END_REF] showed that for any fixed k, one can decide in polynomial time if a given graph has Grundy number at most k. However determining the Grundy number of a graph is NP-hard [START_REF] Zaker | Results on the Grundy chromatic number of graphs[END_REF], and given a graph G, it is even NP-complete to decide if Γ(G) = ∆(G)+1 as shown by Havet and Sampaio [START_REF] Havet | On the Grundy number of a graph[END_REF]. In addition, Asté et al. [START_REF] Asté | Grundy number and products of graphs[END_REF] showed that for any constant c ≥ 1, it is NP-complete to decide if Γ(G) ≤ c • χ(G).

Graph colouring of many graph classes has also been studied. One of the classes is the one of line graphs. The line graph of a graph G, denoted L(G), is the graph whose vertices are the edges of G, with ef ∈ E(L(G)) whenever e and f share an endvertex. Colouring line graphs corresponds to edge-colouring. A k-edge-colouring of a graph G is a surjective mapping φ : E(G) → {1, . . . , k} such that if two edges e and f are adjacent (i.e share an endvertex), then φ(e) = φ(f ). A k-edge colouring may also be seen as a partition of the edge set of G into k disjoint matchings M i = {e | φ(e) = i}, 1 ≤ i ≤ k. By edge-colouring we mean either the mapping φ or the partition.

The chromatic index χ (G) of a graph G is the least k such that G admits a k-edge-colouring. It is easy to see that χ (G) = χ(L(G)). Obviously, ∆(G) ≤ χ (G) and Shannon's and Vizing's theorems state that χ (G) ≤ max{ 3 2 ∆(G); ∆(G) + µ(G)}, where µ(G) is the maximum number of edges between two vertices of G. Holyer [START_REF] Holyer | The NP-completeness of edge-coloring[END_REF] showed that for any k ≥ 3, it is NP-complete to decide if a k-regular graph has chromatic index k.

One can apply the greedy algorithm to colour a line graph. It corresponds to the following greedy algorithm for edge-colouring. Given a graph G = (V, E) and an edge ordering θ = e 1 < • • • < e n , assign to e i the least positive integer that was not already assigned to lower-indexed edges adjacent to it. An edge-colouring obtained by this process is called a greedy edge-colouring and it has the following property: For every j < i, every edge e in M i is adjacent to an edge in M j .

(P) Note that an edge-colouring satisfying (P ) is a greedy edge-colouring relative to any edge ordering in which the edges of M i precede those of M j when i < j.

The Grundy index Γ (G) of a graph G is the largest number of colours of a greedy edgecolouring of G. Notice that Γ (G) = Γ(L(G)). By definition, χ (G) ≤ Γ (G). Furthermore, as an edge is incident to at most 2∆(G) -2 other edges (∆(G) -1 at each endvertex), colouring the edges greedily uses at most 2∆(G) -1 colours. So ∆(G) ≤ Γ (G) ≤ 2∆(G) -1. There are graphs for which the Grundy index equals the maximum degree: stars for example. On the opposite, for any ∆ there is a tree with maximum degree ∆ and Grundy index 2∆(G) -1.

In this paper, we study the complexity of finding the Grundy index of a graph. We prove that it is NP-hard by showing that the following problem is co-NP-complete.

Minimum Greedy Edge-Colouring

Instance: A graph G. Question: Γ (G) = ∆(G)?
The proof is a reduction from 3-Edge-Colourability of Cubic Graphs which was proved to be NP-complete by Holyer [START_REF] Holyer | The NP-completeness of edge-coloring[END_REF]. We recall that a cubic graph is a 3-regular graph. The reduction also proves that it is co-NP-complete to decide if Γ (G) = χ (G).

3-Edge-Colourability of Cubic Graphs

Instance: A cubic graph G. Question: Is G 3-edge colourable?

We then extend the result to a more general problem.

f -Greedy Edge-Colouring Instance: A graph G. Question: Γ (G) ≤ f (∆(G))?
We show that for any function

f such that k ≤ f (k) ≤ 2k -2, the problem f -Greedy Edge-Colouring is co-NP-Complete.
Since determining the Grundy index is NP-hard, a natural question to ask is for which class of graphs it can be done in polynomial time. In Section 3, we consider caterpillars which are trees such that the deletion of all leaves results in a path, called backbone. We show that if T is a caterpillar then Γ (T ) ≤ ∆(T ) + 1 and then give a linear-time algorithm to compute the Grundy index of a caterpillar.

Co-NP-completeness results

The aim of this section is to prove that f -Greedy Edge-Colouring is co-NP-complete for every function f such that k ≤ f (k) ≤ 2k -2 for all k.

For sake of clarity, we first show that Minimum Greedy Edge-Colouring is co-NP-Complete.

Theorem 1. Minimum Greedy Edge-Colouring is co-NP-Complete.

Minimum Greedy Edge-Colouring is clearly in co-NP, because a greedy edge-colouring of a graph G with at least ∆(G) + 1 colours is a certificate that Γ (G) > ∆(G).

We now prove the co-NP-completeness by reduction from 3-Edge-Colourability of Cubic Graphs.

Let H be a cubic graph on n vertices w 1 , . . . w n . Let G be the graph defined by

V (G) = V (H) ∪ {u 1 , . . . , u n } ∪ {v, a, b, c} and E(G) = E(H) ∪ {u i w i | 1 ≤ i ≤ n} ∪ {vu i | 1 ≤ i ≤ n} ∪ {av, bv, cv}. See Figure 1.
In G, d(v) = n + 3, while the degree of all other vertices is at most 4. Thus,

∆(G) = d(v) = n + 3 because n ≥ 4 has H is cubic. Moreover, every edge of G is adjacent to at most n + 3 edges so Γ (G) ≤ n + 4 = ∆(G) + 1. Hence the Grundy index of G is either ∆(G) or ∆(G) + 1. The co- NP-completeness of Minimum Greedy Edge-Colouring follows directly from the following claim. Claim 1. χ (H) = 3 if and only if Γ (G) = ∆(G) + 1.
Proof. (⇒) Suppose that there exists a 3-edge-colouring φ of H. Let us extend φ into a greedy edge-colouring of G with ∆(G)

+ 1 = n + 4 colours. Set φ(av) = 1, φ(bv) = 2, φ(cv) = 3, and for all 1 ≤ i ≤ n, φ(u i w i ) = 4 and φ(u i v) = i + 4. Notice that every vertex w i is incident to an Inria Figure 1: Graph G obtained from a cubic graph H. edge of H of each colour in {1, 2, 3} since H is cubic. Then it is straightforward to check that φ is a greedy (n + 4)-edge-colouring of G.
(⇐) Suppose that there is a greedy (n + 4)-edge-colouring of G. Some edge is coloured n + 4. But such an edge has to be adjacent to at least n + 3 edges and thus to be one of the vu i , say vu n . The edge vu n is adjacent to exactly n + 3 edges. So by Property (P), all edges adjacent to vu n receive distinct colours in {1, . . . , n + 3}.

Let us first prove by induction on 1 ≤ j ≤ n that the edge e j incident to vu n labelled n + 5 -j is one of the vu i , the result holding for j = 1. Suppose now that j ≥ 2. The edge e j must have degree at least n+5-j since it is adjacent to vu n and one edge of each colour in {1, . . . , n+4-j} by Property (P). Hence e j must be incident to v since u n w n is adjacent to four edges. Then e j must have degree at least n + 3 since it is adjacent to the j -1 edges e l for 1 ≤ l < j and one edge of each colour in {1, . . . , n + 4 -j}. Hence e j is one of the vu i .

Hence, without loss of generality, we may assume that φ(vu i ) = i + 4 for all 1 ≤ i ≤ n. The edge vu i is adjacent to an edge coloured 4. This edge must be u i w i since the edges av, bv and cv are adjacent to at most 2 edges coloured in {1, 2, 3}. Thus φ(u i w i ) = 4 for all 1 ≤ i ≤ n. Now every edge u i w i is adjacent to three edges, one of each colour in {1, 2, 3}. Since φ(vu i ) ≥ 5, these three edges must be the three edges incident to w i in H. Thus all the edges of H are coloured in {1, 2, 3}. Hence the restriction to φ to H is a 3-edge-colouring.

Remark 1. Observe that the graph G has chromatic index ∆(G). Indeed colour the edges adjacent to v with the colours 1, . . . , ∆(G) and then extend greedily this colouring to the other edges. Since all the remaining edges are adjacent to at most 6 edges they will all get a colour less or equal to 7. Since ∆(G) ≥ 7, we obtain a ∆(G)-edge colouring. Hence the above reduction shows that it is co-NP-complete to decide if Γ (G) = χ (G).

Theorem 1 may be generalized as follows.

Theorem 2. Let f be a function such that k ≤ f (k) ≤ 2k -2 for all k ∈ N. f -Greedy Edge-Colouring is co-NP-Complete.
Proof. f -Greedy Edge-Colouring is clearly in co-NP, because a greedy edge-colouring of a graph G with more than f (∆(G)) colours is a certificate that Γ (G) > f (∆(G)).

We now prove the co-NP-completeness by reduction from 3-Edge-Colourability of Cubic Graphs.

Let H be a cubic graph on n vertices w 1 , . . . , w n and let G be the graph defined as in the proof of Theorem 1. Set p = f (n + 3) -(n + 3). Then 0 ≤ p ≤ n + 1. For 1 ≤ i ≤ p, let T i be the tree with vertex set {a i , b i , c i , t i } ∪ {a i,j , b i,j , c i,j , s i,j , t i,j | 1 ≤ j ≤ n -1} and edge set

{a i t i , b i t i , c i t i } ∪ n-1
j=1 {a i,j t i,j , b i,j t i,j , c i,j t ij , t i,j s i,j , s i,j t i }. Let G be a graph obtained from the disjoint union of G and the T i by adding the edge u n t i for all 1 ≤ i ≤ p. See Figure 2. (⇒) Suppose that there exists a 3-edge-colouring φ of H. Let us extend φ into a greedy edge-colouring of G with f (∆(G )) + 1 = n + p + 4 colours. We first extend it into a greedy (n + 4)-colouring of G as we did in the proof of Theorem 1. In particular, we have φ(u n w n ) = 4 and φ(u n v) = n + 4. For all 1 ≤ i ≤ p and all 1 ≤ j ≤ n -1, we set φ(t i a

i ) = 1, φ(t i b i ) = 2, φ(t i c i ) = 3, φ(t i,j a i,j ) = 1, φ(t i,j b i,j ) = 2, φ(t i,j c i,j ) = 3, φ(t i,j s i,j ) = j+3, and φ(t i u n ) = n+4+i.
Then it is straightforward to check that φ is a greedy (n + p + 4)-edge-colouring of G .

(⇐) Suppose that G admits a greedy (n + p + 4)-edge-colouring φ. For all 1 ≤ i ≤ p, there is an edge e i coloured n + 4 + i. This edge must has to be adjacent to at least n + 3 + i edges by Property (P). So all the e i must be in F = {vu n } ∪ {u n t i | 1 ≤ i ≤ p}. Now the edge e p is adjacent to an edge e 0 coloured n + 4. This edge is adjacent to at least n + 4-edges, one of each colour in {1, . . . , n + 3} and e p . Hence it also has to be in F . Since |F |, all the edges in F are coloured with distinct labels in {n + 4, . . . , n + p + 4}. Now applying the same reasoning as in the proof of Theorem 1, we derive that the restriction to φ to H is a 3-edge-colouring.

Greedy edge-colouring of caterpillars

In this section, we show a polynomial-time algorithm solving Greedy Edge-Colouring for caterpillars. Recall that caterpillars are trees such that the deletion of all leaves results in a path, called backbone.

Inria

We first show some properties of a greedy edge-colouring of a caterpillar where we can see that the Grundy index of a caterpillar T is at most ∆(T ) + 1. We then give a polynomial-time algorithm that computes its Grundy index.

Grundy index of a caterpillar

Lemma 3. Let T be a caterpillar and v a vertex in its backbone. In every greedy edge-colouring of T , the colours 1, . . . , d(v) -2 appear on the edges incident to v.

Proof. By the contrapositive. Let c be a greedy edge-colouring of T . Suppose that a colour α ∈ {1, . . . , d(v) -2} is not assigned to any edge incident to v. Then, since all the edges incident to v have different colours, at least three colours strictly greater than d(v) -2 appear on three edges incidents to v. One of these colours, say β must appear on an edge e incident with a leaf. But e is uniquely adjacent to edges incident to v. So e is adjacent to no edge coloured α. Since α ≤ d(v) -2 < β, the edge-colouring c is not greedy. Proof. Suppose by way of contradiction that one of these two edges, say e 1 , is incident to a leaf. Then e 1 is adjacent to d(v) -1 other edges, and one of them, namely e 2 , is assigned a colour greater than d(v) -1. Thus e 1 is adjacent to at most d(v) -2 edges whose colour is less or equal to d(v) -1. So, there is a colour α in {1, . . . , d(v) -1} such that no edge incident to e 1 is coloured α. This contradicts the fact that c is greedy. Hence Proof. Set ∆(T ) = ∆. Suppose by way of contradiction that it is possible to greedily colour T with ∆ + 2 colours. Let e be an edge coloured ∆ + 2. It must be adjacent to at least ∆ + 1 edges, one of each colour 1, . . . , ∆ + 1. Thus, the edge e is in the backbone. According to Lemma 4, the edges e 1 and e 2 adjacent to e with colours ∆ and ∆ + 1 are in the backbone. Furthermore all the edges adjacent to e which are neither e 1 nor e 2 are coloured in {1, . . . , ∆ -2}. Hence e is adjacent to no edge coloured ∆ -1, a contradiction. Theorem 5 is tight. For instance, consider the caterpillar T k with backbone (t, u, v, w) for which the vertex t has degree k -1, u has degree 2 and v and w degree k. An edge-colouring in which the k -2 edges incident to t and a leaf are coloured with 1, . . . , k -2, the k -1 edges incident to w and a leaf with 1, . . . , k -1, the k -2 edges incident to v and a leaf with 1, . . . , k -2, the edge tu with k -1, the edge vw with k and the edge uv with k + 1 is greedy.

Finding the Grundy index of a caterpillar

Let us now examine in more details when a caterpillar T has Grundy index ∆(T ) + 1.

Let T be a caterpillar with backbone P = (v 1 , v 2 , . . . , v n ). The first edge of P is v 1 v 2 . For any edge e = v i v i+1 ∈ P , removing e from T gives two caterpillars T - e and T + e , the first one containing v i and the second one containing v i+1 . For convenience, the backbone of T - e is

P - e = (v i , v i-1 , . . . , v 1 ) and the backbone of

T + e is P + e = (v i+1 , . . . , v n ). Hence the first edge of T - e is (v i , v i-1
) and the first edge of T + e is (v i+1 , v i+2 ).

Lemma 6. Let T be a caterpillar of maximum degree ∆ with backbone P = (v 1 , . . . , v n ). Then Γ (T ) = ∆ + 1 if and only if there is an edge e ∈ E(P ) \ {v 1 v 2 , v n-1 v n } such that 1. one endvertex of e has degree ∆, and 2. one of the two caterpillars T - e and T + e has a greedy edge-colouring such that the first edge of its backbone is coloured ∆ and the other has a greedy edge-colouring such that its first edge of its backbone is coloured ∆ -1.

Proof. Assume that T has a greedy (∆ + 1)-edge-colouring. Let e be an edge coloured ∆ + 1. By Lemma 4, e is in the backbone and incident to a vertex of degree ∆, proving (i). Moreover, the edge e is adjacent to an edge coloured ∆ and another one labelled ∆-1. Again by Lemma 4 these two edges must also be in the backbone. In particular, e is not v 1 v 2 nor v n-1 v n because these two edges are adjacent to a unique edge of the backbone. Moreover the greedy edge-colourings induced on T - e and T + e clearly satisfy (ii). Conversely, assume that there is an edge e ∈ E(P ) \ {v 1 v 2 , v n-1 v n } satisfying (i) and (ii). Let φ -and φ + be the greedy edge-colourings of T - e and T + e respectively as in (ii). Let φ be the edge-colouring of T defined by φ(e) = ∆ + 1, φ(f ) = φ -(f ) for all f ∈ T - e and φ(f ) = φ + (f ) for all f ∈ T + e . We claim that φ is a greedy edge-colouring. Clearly, since φ -and φ + are greedy, it suffices to prove that e is adjacent to an edge of every colour i in {1, . . . , ∆}. Since φ + and φ - satisfy (ii), then e is adjacent to an edge labelled ∆ and an edge labelled ∆-1. Now, e is incident to a vertex v of degree ∆. This vertex is incident to e and an edge f in the backbone. The edge f is the first edge of a tree T f in {T + e , T - e }. In the greedy edge-colouring of T f , the edge f has a colour greater than ∆ -2, so the ∆ -2 edges incident to v which are not e nor f have all one colour in 1, . . . , ∆ -2. Hence e is adjacent to an edge of every colour in {1, . . . , ∆}.

In view of this lemma, it is useful to decide when a caterpillar has a greedy edge-colouring with a prescribed colour on the first edge of its backbone.

Lemma 7. Let T be a caterpillar with backbone P with first edge is e = uv. Then T has a greedy edge-colouring such that e is coloured k if and only if one of the following holds:

1. d(u) ≥ k or d(v) ≥ k; 2. d(u) = k -1 and T +
e admits a greedy edge-colouring such that the first edge of P + e is coloured k -1.

Proof. Let e = uv with u the first vertex of P . Assume first that T has a greedy edge-colouring such that e is coloured k and that e is incident to no vertex of degree k. Then the edges incident to u must be coloured by 1, . . . , d(u) -1 and the edges incident to v and a leaf are coloured by 1, . . . d(v) -2. Hence the edge incident to e and coloured k -1 must be the first edge of P + e is coloured k -1 by Property (P). So the edge incident to e and coloured k -2 must be incident to u, and thus d

(u) -1 ≥ k -2, that is d(u) ≥ k -1.
Assume now that (i) holds. Let x be a vertex in {u, v} with degree at least k. One can colour all the edges incident to x with 1, . . . , d(v) such that e is coloured k and then extend this edge-colouring greedily to obtain the desired greedy edge-colouring of T .

Finally assume that (ii) holds. Let φ be a greedy-edge colouring of T + e such that the first edge of P + e is coloured k -1. One can extend it by assigning k to e, 1, . . . , k -2 to the k -2 edges incident to u and leaves . It is routine to check that this is a greedy-edge colouring of T .

Inria

Lemma 7 translates easily in a linear-time recursive algorithm. Using this, and Lemma 6, one derive a linear-time algorithm for computing the Grundy index of a given caterpillar.

Theorem 8. Determining the Grundy index of a caterpillar T can be done in O(|V (T )|).

Proof. Theorem 5 and Lemma 6 imply that Algorithm 1 return the Grundy index of T provided that we have a subroutine FirstEdge(T ,P ,k) that returns 'yes' if a caterpillar T with backbone P admits a greedy-edge colouring such that the first edge of P is coloured k. (1) Let P = (v 1 , v 2 , ..., v n ) be the backbone of T . Compute d(v i ) for all i ≤ v n and compute ∆ = ∆(T ).

(2) for i = 2 to n -2 do Let us now examine the complexity of Algorithm 1. Let us first observe that FirstEdge(T , P , k) makes a constant number of operations before calling FirstEdge(T -u, P -u, k -1). Hence an easy induction show that it makes O(k) operations in total.

Figure 2 :

 2 Figure 2: The graph G obtained from a cubic graph H.

Claim 2 .

 2 χ (H) = 3 if and only if Γ (G ) = f (∆(G )) + 1.

Lemma 4 .

 4 Let c be a greedy edge-colouring of a caterpillar T and v a vertex in the backbone of T . If two edges e 1 and e 2 incident to v receive colours greater than d(v) -1, then e 1 and e 2 are two edges of the backbone and the edges incident to v and leaves are coloured 1, . . . , d(v) -2.

e 1 and e 2

 2 are edges of the backbone. Now by Lemma 3, there must be edges incident to v of each colour in {1, . . . , d(v) -2}. So the d(v) -2 edges distinct form e 1 and e 2 , which are the edges linking v and leaves are coloured in {1, . . . , d(v) -2}. We use Lemma 4 to prove the following theorem. Theorem 5. If T is a caterpillar, then Γ (T ) ≤ ∆(T ) + 1.

Algorithm 1 :

 1 GrundyIndex(T)Input: A caterpillar T . Output: Γ (T ).

( 3 )( 2 )( 4 )

 324 e := v i v i+1 ; (4) if d(v i ) = ∆ or d(v i+1 ) = ∆ then(5) if FirstEdge(T + e , P + e , ∆)=TRUE and FirstEdge(T - e , P - e , ∆ -1)=TRUE then (6) return ∆ + 1; (7) if FirstEdge(T + e , P + e , ∆ -1)=TRUE and FirstEdge(T - e , T - e , ∆)=TRUE then (8) return ∆+1; (9) Return ∆; Such a subroutine FirstEdge may be obtained by Algorithm 2 according to Lemma 7. Algorithm 2: FirstEdge Input: A caterpillar T with backbone P and an integer k. Output: T RU E if there is a greedy k-edge-colouring of T with first edge of P coloured k, and F ALSE otherwise. (1) Let u be the first vertex of P and v its second. (So uv is the first edge.) If d(u) ≥ k or d(v) ≥ k then (3) return T RU E; If d(u) ≥ k -1 then (5) return FirstEdge(T -u, P -u, k -1); (6) return F ALSE;

RR n°8171

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Algorithm 1 first compute (line 1) the degrees of all the v i , which can be done in time O(|V (T )|) and then takes the maximum of all this values which can also be done in time O(|V (T )|.

In a second phase (line 2 to 8), for each edge e ∈ P which is incident to a vertex of degree ∆, Algorithm 1 makes at most four calls of FirstEdge with last parameter ∆ -1 or ∆. Hence for each e ∈ P it makes O(∆) operations. Let S be the set of vertices of degree ∆, The number of edges of P incident to a vertex of degree ∆ is at most 2|S|. But every vertex in S is adjacent to at least ∆ -2 leaves. Hence |V (T )| ≥ |S| + (∆ -2)|S|, so |S| ≤ |V (T )|/(∆ -1). Hence, in this second phase, the algorithm makes at most O 2 × |V (T )| ∆-1 ∆ = O(|V (T )|) operations. Thus, in total, Algorithm 1 makes O(|V (T )|) operations.