Jean-Marc Jézéquel

David Méndez-Acuña

Thomas Degueule

Benoit Combemale

Olivier Barais

When Systems Engineering Meets Software Language Engineering

The engineering of systems involves many different stakeholders, each with their own domain of expertise. Hence more and more organizations are adopting Domain Specific Languages (DSLs) to allow domain experts to express solutions directly in terms of relevant domain concepts. This new trend raises new challenges about designing DSLs, evolving a set of DSLs and coordinating the use of multiple DSLs for both DSL designers and DSL users. This paper explores various dimensions of these challenges, and outlines a possible research roadmap for addressing them. The message of this paper is also to claim that if language engineering techniques to design any single (disposable) language are mature, the language engineering community needs to fundamentally change its view on software language design. We need to take the next step and adopt the perspective that a software language is, fundamentally, software too and thus the result of a composition of design decisions. These design decisions should be represented as first-class entities in the software languages workbench and it should be possible, during the language lifecycle, to add, remove and change language design decisions with limited effort to go from continuous design to continuous meta-design.

A Language-Oriented Vision for Systems Engineering

The engineering of complex software intensive systems involves many different stakeholders, each with their own domain of expertise. It is particularly true in the context of systems engineering in which rather than having everybody working with code/model defined in general-purpose (modeling/programming) languages, more and more organizations are turning to the use of Domain Specific Languages (DSLs). DSLs allow domain experts to express solutions directly in terms of relevant domain concepts, and use generative mechanisms to transform DSL specifications into software artifacts (e.g., code, configuration files or documentation), thus abstracting away from the complexity of the rest of the system and the intricacies of its implementation.

The adoption of DSLs has major consequences on the industrial development processes. This approach, a.k.a. Language-Oriented Programming [START_REF] Ward | Language-oriented programming[END_REF], breakdowns the development process into two complementary stages (see Figure 1): the development, adaptation or evolution by language designers of one or several DSLs, each capitalizing the knowledge of a given domain, and the use of such DSLs by language users to develop the different system concerns. Each stage has specific objectives and requires special skills. Figure 1 depicts the two interdependent processes that continuously drive each other's. The main objective of the language engineering process is to produce a DSL which tackles a specific concern encountered by engineers in the development of a complex system, together with its tooling. Once an appropriate DSL is made available to systems engineers, it is used to express the solution to this specific concern in the final system. However, by definition, DSLs are bounded to evolve with the domain they abstract. Consequently, systems engineers need to be well aware of end users' expectations in order to report their new requirements to the language designers. A new evolved DSL is then produced by the language designers, which is in turn used by systems engineers and so on and so forth. It is worthwhile to note that, although this is unlikely in large companies, these roles can be alternatively played by the same people in smaller organizations.

As a matter of fact, while DSLs have been found useful for structuring development processes and providing abstractions to stakeholders [START_REF] Hutchinson | Empirical assessment of mde in industry[END_REF], their ultimate value has been severely limited by their user-understanding ambiguity, the cost of tooling and the tendency to create rigidity, immobility and paralysis (the evolution of such languages is costly and error-prone). The development of software languages is a challenging task also due to the specialized knowledge it requires. A language designer must own not only quite solid modeling skills but also the technical expertise for conducting the definition of specific artifacts such as grammars, metamodels, compilers, and interpreters. "Software languages are software too" [START_REF] Favre | Empirical language analysis in software linguistics[END_REF] and, consequently, languages development inherits all the complexity of general software development; concerns such as maintainability, reusability, evolution, user experience are recurring requirements in the daily work of software language engineers. As a result, there is room for application of software engineering techniques that facilitate the DSL construction process. This fact permitted the emergence of what we know as Software Language Engineering that is defined as the application of systematic, disciplined, and measurable approaches to the development, use, deployment, and maintenance of software languages [START_REF] Kleppe | Software Language Engineering: Creating Domain-Specific Languages Using Metamodels[END_REF].

The message of this paper is twofold. First, we claim that language engineering techniques for designing disposable DSLs are close to maturity. However, as we will see, some challenges such as composition, modularity or evolution still need to be addressed. Hopefully, decades of research in software engineering already paved the way and software language engineering should leverage these facilities in order to tackle these challenges. Second, we claim that the common view on software language design should fundamentally evolve. Rather than abstract syntax trees, metamodels, type checkers, parsers, code generators, compilers, etc., we need to model and represent a software language as the composition of a set of language design decisions, concerning, among others, the existing language-units solutions, variation points, features and usage scenarios that are needed to satisfy the requirements. Once we are able to represent software languages, in several phases of the lifecycle, in terms of the aforementioned concepts, changing and evolving software languages will be considerably simplified.

The remainder of this paper is organized as follows. We investigate the underlying challenges of the adoption of DSLs in the development of complex software intensive systems firstly from the points of view of the language designer (Section 2) and secondly from the language user point of view (Section 3). Finally, Section 4 draws some perspectives and concludes.

Challenges for SLE from Language Designers' Point of View

From a language designer point of view, the development of new DSLs as well as the evolution of existing ones becomes daily activities. Evolving a DSL usually requires the co-evolution of all its tooling (parsers, textual syntax and graphical syntax editors, compilers, code generators, . . . ). Besides, language users generally require backward compatibility or tooling for supporting the migration. Consequently, each evolution is costly and error-prone and the software language engineering community still needs to come up with new solutions. To enable this vision that the language design decisions should be represented as first-class entities in the software languages workbench and the it should, during the language lifecycle, be possible to add, remove and change language design decisions against limited effort, his section explores some required software engineering techniques that have been used in the context of software languages engineering for improving i) the reuse, ii) the variability management and iii) the verification and validation. Specifically, we highlight the main challenges that remain to be addressed in each case.

Reuse

Reusability of software artifacts is a central notion that has been thoroughly studied and used by both academics and industrials since the early days of soft-ware construction. Essentially, designing reusable artifacts allows the construction of large systems from smaller parts that have been separately developed and validated, thus reducing the development costs by capitalizing on previous engineering efforts. It is however still hardly possible for language designers to design typical language artifacts (e.g. language constructs, grammars, editors or compilers) in a reusable way. The current state of the practice most of the time prevents the reusability of language artifacts from one language to another, or from one system to another, consequently hindering the emergence of real engineering techniques around software languages.

Conversely, concepts and mechanisms that enable artifacts reusability abound in the software engineering community. In this section, we present the timehonored concepts of substitutability, inheritance and components, show their relevance for language designers and draw some perspectives and challenges for their inclusion in software language engineering. Substitutability In its broadest sense, substitutability is the mechanism that allows the replacement of one software artifact (e.g. code, object, module) with another one under certain conditions. In the context of software language engineering, the considered artifacts (languages, models, abstractions, tools, etc.) are all candidates for substitutability mechanisms, allowing reusing them in different contexts. We propose the notion of types as interfaces that express the constraints that different artifacts must verify in order to be substituted one another.

The substitution principle has been thoroughly investigated in the objectoriented programming community. It states whether given objects in a program can safely be substituted to other objects based on the subtyping relation that stands between their types [START_REF] Liskov | A behavioral notion of subtyping[END_REF]. Most object-oriented programming languages feature a mechanism of subtype polymorphism that allows considering the same object through different interfaces (i.e. types), provided they are subtypes one another, enabling facilities such as code reuse or dynamic binding.

In the context of software language engineering, the definition of such types and subtyping relations enables model (i.e. graph of objects) polymorphism, namely the possibility to manipulate a model or program created using a specific language through different tools initially designed for similar yet different languages [START_REF] Steel | On model typing[END_REF]. Model polymorphism allows tackling a wide range of scenarios that are commonly faced by system engineers. As a concrete example, consider the management of evolution on complex languages such as UML. It is difficult for engineers to deal with this evolution, as all efforts concentrated around a language are lost with subsequent versions; e.g. a transformation defined for UML2.1 cannot be reused for models created using UML2.2 since these are, although semantically close, different languages. Specifying the parameter of such a transformation in terms of model interface allows reusing it for any model that matches this interface: if a subtyping relation can be established between the two versions, model polymorphism allows the reuse of all the tooling across them, even benefiting from dynamic binding for prospective specialization.

However, the currently prevalent modeling frameworks do not provide such type of substitutability mechanisms, and only a few recent research works address them (e.g. [START_REF] Lara | Generic meta-modelling with concepts, templates and mixin layers[END_REF][START_REF] Guy | On model subtyping[END_REF][START_REF] Sun | Using model types to support contract-aware model substitutability[END_REF]). Challenges such as complete language semantics substitutability or concrete syntax replacement still need to be addressed. More generally, using interfaces for specifying the expected features and properties of a language paves the way for language-agnostic, generic manipulation of models and programs. This is particularly relevant in system engineering as engineers need to deal with many different domains and stakeholders, each using his own domain-specific, independently-evolving language.

Extension The need for language extension arises in many scenarios. DSLs are initially designed and implemented for a restricted set of users with a finite set of features that support their requirements but, most of the time, new requirements will emerge once the language gets effectively used: they tend to grow with the users' needs. Moreover, DSLs are now more and more scattered among different set of users that tackles the same domain, but with their own specificities. In this case, language designers should be able to reuse an existing DSL that contains the basic constructs and features for a particular domain, and extend it with their own business distinctive features.

To support such scenarios where extensions are most of the time unforeseen, DSLs must be designed in a way that facilitates their reuse and extensibility. Conversely, a language designer that extends an existing language should be able to concentrate on its business specificities, seamlessly reusing the base language along with all its tooling. In this regard, real extensibility mechanisms should support the introduction of new constructs, abstractions, or tools without having to understand or recompile the base language source code.

Modularity and composability Modularization of software (i.e., componentsbased software development) is considered an effective mechanism for achieving software reuse and, consequently, reducing costs and time to market. The main principle is to structure software applications as sets of interconnected building blocks that, in turn, are designed in such a way that allows later re-use in other applications. In this context, three of the most important challenges are (1) design and implement components with real potential re-use; (2) design the interfaces that enable crosscutting collaboration among components in a system; and (3) provide components models [START_REF] Lau | Software component models[END_REF] that offer composition mechanisms for integrating a set of components.

All the aforementioned ideas apply also when the software under construction is a software language; there are benefits in terms of the reduction of construction effort [START_REF] Cleenewerck | Component-based dsl development[END_REF]. Nevertheless, those general challenges gain some special connotations analyzed below: i) Components design (how to breakdown a language?): Decomposing a language in several language modules (a.k.a., language units) is not only about offering a languages benchmark that enables modular definition of metamodels, grammars and semantics; it is just the technical part. The other important challenge is to understand how the language units should be defined so they can be reused in other contexts. What is the correct level of granularity? What are the "services" that a language unit should offer for being considered reusable? What is the meaning of a "service" in the context of software languages? What is the meaning of a "services composition" in the context of software languages? ii)Languages interfaces (how language units are specified?): The construction of a language unit is not only about implementing a subset of the language but also about specifying its boundary (i.e., the set of services it offers to other language units and the set of services it requires from other language units). This fact refers to the classical idea of required and provided interfaces introduced by components-based software engineering approaches. But... what is the meaning of "provided and required services" in the context of software languages? We argue that the answer to that question must consider at least two facts: composability and (one more time) substitutability. In the case of composability, required and provided interfaces should provide a mechanism for exposing providing services so they can be consumed by the required services of other language units. In other words, interfaces are a mechanism for interaction between language units. By the other hand, substitutability refers to the possibility of implementing provided services in different language units so the provided interface becomes also a set of constraints that ensure the safe replacement of the given implementations. iii) Language units composition (how two language units do interact?): The nature of the interaction between two language units might be different depending on some architectural decisions; extending one language unit with another one is a different situation from a required service of a language unit consuming one or several provided services from other language units. In the case of extension, the base language unit is usually independent of the extensions whereas the extensions have little sense without the base language unit [START_REF] Erdweg | Language composition untangled[END_REF]. In the case of required and provided interactions, the requiring language unit usually cannot work without the provided one. We argue that there is a need of composition operators that explicitly define the role of each language unit in a composition.

Variability management and language families

One of the main limitations of components-based software development -and of course it is also true in the case of software languages modularization-is the difficulty of designing components that can be actually re-used in other contexts. Despite the design principles and patterns, reusing of software components is not guaranteed. In fact, in many cases the effort of building modularized software is not compensated with the re-usability opportunities.

One of the answers that the software community has found is the idea of variability management. Variability management is a mechanism for explicitly representing the commonalities and differences among a family of software products. A family of products is defined as a set of software applications that have similar purposes and that share some functionality but that is specialized in a particular type of users or situation. The idea is to effectively reuse the implementation of such common functionality and having a repository of "common assets" that implement product features. The process of creating a product by using the family of products is called product derivation. To do so, it is necessary to select the desired product features and to offer a mechanism of composition for integrating the assets corresponding to each feature. This is the main principle of what we know as Software Product Line Engineering (SPLE) that is a software engineering approach that has demonstrated important benefits in the general case of software development.

As demonstrated in [START_REF] Vacchi | Variability support in domainspecific language development[END_REF], variability management -and the ideas behind SPLC in general-can be applied in the context of software languages for increasing the re-usability and then increasing the productivity of software language engineers. In this context, a family of products actually is a family of languages where there are some commonalities and some differences (consider as an example the family of OCL variants presented in [START_REF] Wende | A role-based approach towards modular language engineering[END_REF]).

Some of the challenges that should be considered are:

Alignment with the modularization approach: It is worth noting that modularization is a prerequisite for addressing variability management. In fact, at the implementation level software modularization and variability management are strongly linked. Each concrete feature expressed in the variability model must correspond to a software component in the architecture so a given configuration can be derived in a concrete functional product. In the case of software languages each feature should be mapped to one (or more) language units that offers the corresponding services. Moreover, in [START_REF] Van Der Linden | Software Product Lines in Action: The Best Industrial Practice in Product Line Engineering[END_REF] van der Linden et. al. present a set of three variability realization techniques at the level of the software modularization schema. Those techniques can be viewed as a set of requirements in terms of modularization and composition of the architecture and they are quite related with the concepts of extension, substitutability and adaptation, some of them discussed in the previous section. How to conjugate all those concepts for effectively define an approach that allows the construction of families of software languages?

Multi-stage orthogonal variability modeling: Typically, a software language specification is intended to define the abstract syntax, the concrete syntax and the semantics of a language. As a result, language units have to contribute to each of those dimensions. In other words, each language unit specification includes a partial definition of the abstract syntax, the concrete syntax, and the semantics. The whole language specification is obtained by putting all the language units together. In [START_REF] Grönniger | Modeling language variability[END_REF] the authors observed that there exists some variability between each of those dimensions. Thereby, one language construct (i.e., a concept in the abstract syntax) may be represented in several ways (i.e., several possible concrete syntaxes) and/or may have different meanings (several possible semantics). This analysis remains the same for both the whole language specification and each segment defined in language units. Consequently, we have at least three different dimensions of variability each of them regarding one field of the tuple:

-Abstract syntax variability or "functional variability": This variability refers to the capability of selecting the desired language constructs for a particular product as long as the dependencies are respected. Consider for example a family of languages for state machines where concepts such as timed transitions, composite states, or history pseudo-states are optional and are only included if the user of the language needs them. This variability dimension is quite similar to the classical concept of functional variability of SPLC where each feature represents a piece of functionality that may be or not included depending on the specific requirements of a user. -Concrete syntax variability or "representation variability": This variability refers to the capability of offering different representations for the same concept. Consider for example a language for state machines that can have textual or graphical representations. Note that this type of variability is especially relevant in the context of metamorphics DSLs that we will explain later in this paper. -Semantics variability or "interpretation variability": This variability refers to capability of offering the different interpretations to the same concept. Consider for example the semantics differences that exist between state machines languages explored in [START_REF] Crane | Uml vs. classical vs. rhapsody statecharts: not all models are created equal[END_REF]. In that work, we can see how, for example, the priorities between conflicting transitions in a state machine are resolved with different criteria. If we are able to manage such variability, the reuse opportunities are drastically increased since we are able to reuse the entire language infrastructure (e.g., editors, abstract syntax trees) for the implementation of different languages that are interpreted according to the needs of specific users.

Note that both representation variability and interpretation variability depend on the functional variability. It makes no sense to select a representation (or interpretation) for a language construct that has not been included as part of the language product. In other words, the configuration of representation and interpretation must be performed only for the construct selected in the functional variability resolution.

Verification & Validation

Just as any complex software artifact, software languages need to be thoroughly verified and validated. Their complex nature, the different aspects that compose them makes it particularly difficult: is a language really suited for the problems it tries to tackle? Can all programs relevant for a specific domain be expressed in a precise and concise manner? Are all valid programs correctly handled by the interpreter? Does the compiler always generate valid code?

Different techniques have been developed for the V&V of traditional software and are good candidates for adaptation to software languages: among them, we focus on design-by-contract and software testing, and the challenges they need to address for the engineering of software languages.

Design-by-contract [START_REF] Meyer | Applying 'design by contract[END_REF] advocates the definition of precise interfaces for software components, e.g. using preconditions, postconditions, invariants or types. Contracts can then be checked at different levels to assess the correct interaction of components. In the context of software languages, contracts may be defined on the abstract syntax (e.g. using invariants), on the semantics (e.g. using preconditions and postconditions), etc [START_REF] Sun | Using model types to support contract-aware model substitutability[END_REF]. Design-by-contract is especially relevant for system engineering as it raises the level of abstraction in which the interaction between the different domains and languages is considered, and makes explicit some of the original requirements on the language. An integrated design-bycontract process for software languages engineering is expected to bring the same benefits as in traditional software development: precise -structural and behavioral -interfaces, improved error handling, specification-driven definition of artifacts, etc.

Software testing, on the other side, is the most prevalent V&V technique in software engineering. Testing software languages is a challenging activity since all their aspects must be checked: abstract syntax, grammar, semantics, tooling, etc. Furthermore, in this context, test data (i.e. models or programs) are themselves complex artifacts, thus complicating the coverage of representative inputs and the definition of oracles functions [START_REF] Baudry | Barriers to systematic model transformation testing[END_REF]. The extensive use of generative programming techniques also raises additional problems due to the gap between generation time and testing time; i.e. the to-be-tested generated artifacts are not known yet when the tests are written. Workarounds on this issue include the automatic generation of test cases together with generated artifacts, which in turns increases the testing activity complexity. Finally, the inherent nature of multi-languages engineering requires not only the different languages to be tested, but also their combination and interaction. Such integration tests should be dedicated to the verification and validation of the composition, reusing the testing effort spent on each of its part.

Challenges for SLE From the Language Users' Point of View

From the perspective of the users the emergence of several software languages is also challenging. Despite the overall purpose of constructing DSLs is to facilitate the daily work of systems engineers, dealing with several languages implies not only learning new syntaxes but also interacting with an increasing number of tools: editors, compilers, and code generators among others. The remainder of this section is dedicated to explore some of those challenges that must be addressed to serve this vision that the language design decisions must be represented as first-class entities in the software languages workbench.

Language viewpoint

Domain-Specific Languages (DSLs) are plain languages, in the sense that many difficult design decisions must be taken during their development and mainte-nance, and that they can take different shapes: plain-old to more fluent APIs; internal or embedded DSLs written inside an existing host language; external DSLs with their own syntax and domain-specific tooling. All forms of DSLs have strengths and weaknesses -whether you are a developer or a user of a DSL. The basic trade-offs between internal and external DSLs have already been identified and are subject to extensive discussions and research for several years. A new trend though is observed. DSLs are now so widespread that very different users with separate roles and varied objectives use them. Depending on the kinds of users, roles or objectives, the same form of DSL (external or internal) might not be the best for everybody. Beyond the unification of the different approaches, it is worthwhile for DSLs to support the ability to be self-adaptable to the most appropriate form (including the corresponding IDE) according to a particular usage or task: we call such a DSL a metamorphic DSL.

From the same language description and different interface specifications, we envision the ability to derive various IDEs that can be used accordingly. This vision raises many challenges: systematic methods to evaluate when a form of a DSL meets the expected properties (e.g., learnability); artefact modularization; information sharing, while being able to visualize and manipulate an artefact in a particular representation and in a particular IDE; global mechanism to ensure consistency of the artefacts between these heterogeneous IDEs.

Language evolution

By definition, DSLs are bbounded to evolve with the domain they abstract. Consequently, DSLs users need to learn and understand the newly-created abstractions, syntaxes and tools. This raises new challenges in terms of change management and learnability of languages. In order to facilitate the transition, the migration of models from one version of a language to another, aka co-evolution, must be fully supported by the workbench: automatic migration when possible, diff computation with explicit user refinement when required, etc. The same situation arises when an older language is replaced with a completely new one, defined independently, but abstracting the same domain.

Language integration

The development of modern complex software-intensive systems often involves the use of multiple DSLs that capture different system aspects. In addition, models of the system aspects are seldom manipulated independently of each other. System engineers are thus faced with the difficult task of relating information presented in different models. For example, a system engineer may need to analyze a system property that requires information scattered in models expressed in different DSLs. Current DSL development workbenches provide good support for developing independent DSMLs, but provide little or no support for integrated use of multiple DSLs. The lack of support for explicitly relating concepts expressed in different DSMLs makes it very difficult for developers to reason about information spread across different models.

Supporting coordinated use of DSLs leads to what we call the globalization of modeling languages [START_REF] Combemale | Globalizing modeling languages[END_REF], that is, the use of multiple modeling languages to support coordinated development of diverse aspects of a system. The term "globalization" is used to highlight the desire that DSLs developed in an independent manner to meet the specific needs of domain experts, should also have an associated framework that regulates interactions needed to support collaboration and work coordination across different system domains.

Globalized DSLs aim to support the following critical aspects of developing complex systems: communication across teams working on different aspects, coordination of work across the teams, and control of the teams to ensure product quality.

Conclusion and Perspectives

This paper claims that research conducted in SLE for systems engineering should consider that:

-The first phase of research and development in SLE has matured the technology to a level where industry adoption is wide-spread and few fundamental issues remain for efficiently designing any single (disposable) DSL. -The traditional view on SLE suffers from a number of key problems that cannot be solved without changing our perspective on the notion of language, and especially of DSL. These problems include i) the lack of first-class representation of design decisions in DSL: since design decisions are cross-cutting and intertwined, they are easy to forget and hard to change, leading to high maintenance costs; ii) the lack of support for explicitly relating different DSLs that makes it very difficult for systems engineers to use multiple DSLs while enabling a coordinated development of the diverse system aspects, and to reason about information spread across artifacts built with different DSLs. -As a community, we need to take the next step and adopt the perspective that a software language is, fundamentally, software too, that is, the result of a composition of design decisions. These design decisions should be represented as first-class entities in the software language workbench and it should, during the language lifecycle, be possible to add, remove and change language design decisions with limited effort to go from continuous design to continuous meta-design.

Fig. 1 :

 1 Fig. 1: Language engineering stakeholders

5 Acknowledgments This work is partially supported by the ANR INS Project GEMOC (ANR-12-INSE-0011), the bilateral collaboration between INRIA and Thales Research & Technology VaryMDE, and the ITEA2 European project MERgE.


José Ignacio Orlando

Matthew Blaschko

Learning fully-connected CRFs for blood vessel segmentation in retinal images

Keywords: Blood vessel segmentation, Fundus imaging, Conditional Random Fields, Structured Output SVM

In this work, we present a novel method for blood vessel segmentation in fundus images based on a discriminatively trained, fully connected conditional random field model. Retinal image analysis is greatly aided by blood vessel segmentation as the vessel structure may be considered both a key source of signal, e.g. in the diagnosis of diabetic retinopathy, or a nuisance, e.g. in the analysis of pigment epithelium or choroid related abnormalities. Blood vessel segmentation in fundus images has been considered extensively in the literature, but remains a challenge largely due to the desired structures being thin and elongated, a setting that performs particularly poorly using standard segmentation priors such as a Potts model or total variation. In this work, we overcome this difficulty using a discriminatively trained conditional random field model with more expressive potentials. In particular, we employ recent results enabling extremely fast inference in a fully connected model. We find that this rich but computationally efficient model family, combined with principled discriminative training based on a structured output support vector machine yields a fully automated system that achieves results statistically indistinguishable from an expert human annotator.

Introduction

Retinal images-also known as fundus images or retinographies-are projective color images of the inner surface of the human eye (Figure 1). They allow physicians to observe the retina and its internal parts, including the vascular tree, the optic disc, the fovea, etc. [START_REF] Fraz | Blood vessel segmentation methodologies in retinal images-a survey[END_REF]. Retinal blood vessels provide useful information for several applications, including medical diagnosis and screening of ophthalmological and cardiovascular diseases. The vessels tortuosity, for instance, can be used to characterize hypertensive retinopathy, and the measurement of vessels diameter is utilized to diagnose hypertension [START_REF] Miri | Retinal image analysis using curvelet transform and multistructure elements morphology by reconstruction[END_REF]. Changes in the vasculature distribution are also interpreted as a possible symptom of diabetic retinopathy [START_REF] Fraz | Blood vessel segmentation methodologies in retinal images-a survey[END_REF]. Moreover, vessels are used to assist laser surgeries [START_REF] Kanski | Synopsis of Clinical Ophthalmology[END_REF] and as landmarks for image registration [START_REF] Li | Registration of OCT fundus images with color fundus photographs based on blood vessel ridges[END_REF]. All these applications require the segmentation of the retinal vasculature. In current best practice, trained specialists delineate the vessels manually, although this is a particularly tedious and time-consuming task. Difficulties in the imaging process-such as the inadequate contrast between vessels and background, and uneven background illumination-increase the variability among segmentations performed by different human observers. These facts motivate the development of automatic strategies for blood vessel segmentation without human intervention.

In the last two decades, several approaches have been introduced to solve this problem, but none of them have yet proved to be accurate enough to be assumed as a standard by the medical community [START_REF] Fraz | Blood vessel segmentation methodologies in retinal images-a survey[END_REF]. These segmentation algorithms can be classified into two general categories, supervised and unsupervised. Supervised methods require a set of labeled training samples to learn a model. These samples are composed of pixels with known annotations and their features. Most of the effort in supervised segmentation involves finding new features for training, or better classifiers to perform the pixel classification task. For instance, Ricci and Perfetti [START_REF] Ricci | Retinal blood vessel segmentation using line operators and support vector classification[END_REF] combine a linear support vector machine (SVM) with a line detector feature. Becker et al. [START_REF] Becker | Supervised feature learning for curvilinear structure segmentation[END_REF] learn the features in a supervised way, based on a gradient boosting framework. Unsupervised methods usually involve systems based on clustering, region-oriented approaches [START_REF] Mendonca | Segmentation of retinal blood vessels by combining the detection of centerlines and morphological reconstruction[END_REF][START_REF] Martinez-Perez | Segmentation of blood vessels from red-free and fluorescein retinal images[END_REF] or thresholding after vessel enhancement with filters [START_REF] Miri | Retinal image analysis using curvelet transform and multistructure elements morphology by reconstruction[END_REF][START_REF] Zhang | Retinal vessel extraction by matched filter with first-order derivative of Gaussian[END_REF], line detectors [START_REF] Ricci | Retinal blood vessel segmentation using line operators and support vector classification[END_REF][START_REF] Nguyen | An effective retinal blood vessel segmentation method using multi-scale line detection[END_REF] or morphological operators [START_REF] Fraz | Ensemble classification system applied for retinal vessel segmentation on child images containing various vessel profiles[END_REF]. More complex image processing operations have been combined with vessel centerline detection at different scales to obtain the final segmentation [START_REF] Fraz | Ensemble classification system applied for retinal vessel segmentation on child images containing various vessel profiles[END_REF][START_REF] Vlachos | Multi-scale retinal vessel segmentation using line tracking[END_REF]. The main advantage of unsupervised methods is that they do not need to be trained using manual annotations. However, they have reported worse results than supervised approaches [START_REF] Nguyen | An effective retinal blood vessel segmentation method using multi-scale line detection[END_REF].

Conditional Random Fields (CRFs) have been widely utilized to solve segmentation problems in several applications [START_REF] Li | Markov Random Field Modeling in Image Analysis[END_REF]. To the best of our knowledge, however, they have not been yet applied to retinal blood vessel segmentation in fundus images. This is likely due to that standard pairwise potentials such as in a Potts model assign a low prior to the elongated structures that comprise a vessel segmentation. We overcome this by using a much richer class of potentials.

The main contribution of our paper consists of an automatic method for blood vessels segmentation in retinal images based on fully-connected CRFs. We follow the efficient inference approach proposed in [START_REF] Krähenbühl | Efficient inference in fully connected CRFs with Gaussian potentials[END_REF], and we learn its configuration in a supervised way, using a Structured Output SVM (SOSVM) [START_REF] Joachims | Cutting-plane training of structural SVMs[END_REF]. We validate the approach on the main benchmark data set for vessel segmentation in fundus images, the DRIVE data set [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF]. In contrast with several published works, no test data was utilized to adjust the parameters of the method. Instead, we randomly split the original training data into two subsets, one for training and the other for validation. Once we found the best configuration using the validation set, we applied it over the images in the test set. In this statistically principled setting, we report a fully-automated system that achieves a performance statistically indistinguishable from an expert human annotator.

The remainder of this paper is organized as follows: the formulation of the fullyconnected CRF is described in Section 2; in Section 3 we explain how we learn the CRF model using a SOSVM. In Section 4 we include our results and a comparison against other state-of-the-art works. Finally, Section 5 concludes the paper.

Fully-Connected CRF segmentation

We pose the segmentation task as an energy minimization problem in a fully-connected conditional random field (CRF). In the original definition of CRFs, images are mapped to graphs, where each pixel represents a node, and every node is connected with an edge to its immediate neighbors [START_REF] Li | Markov Random Field Modeling in Image Analysis[END_REF]. In the fully-connected version, each node is assumed to be a neighbor of every other. Following this approach the method is able to take into account not only neighboring information but also long-range interactions between pixels. This property improves the segmentation accuracy, but makes the inference process computationally expensive. Recently, however, Krähenbühl and Koltun [START_REF] Krähenbühl | Efficient inference in fully connected CRFs with Gaussian potentials[END_REF] have introduced an efficient inference approach under the restriction that the pairwise potentials are a linear combination of Gaussian kernels over a Euclidean feature space. This approach, which is based on taking a mean field approximation of the original CRF, is able to provide accurate segmentations in less than a second.

We denote y = {y i } as a labeling over all pixels of the image in the label space L = {-1, 1}, where 1 is associated to blood vessels and -1 to any other class. Its corresponding energy function is given by the sum of its unary energy ψ u and its pairwise energy ψ p :

E(y) = i ψ u (y i , x i ) + i<j ψ p (y i , y j , f i , f j ) (1)
where x i and f i are the unary and pairwise features, respectively. Unary potentials define a log-likelihood over the label assignment y, and they are traditionally computed by a classifier [START_REF] Krähenbühl | Efficient inference in fully connected CRFs with Gaussian potentials[END_REF]. Pairwise potentials define a similar distribution but considering only the interactions between pixels features and their labels. Parameters for both potentials are learned using a Structured Output SVM, as we explain in detail in Section 3.

Unary potentials: We obtain the unary potentials according to the following expression:

ψ u (y i , x i ) = -w uy i , x i -β yi (2)
where w uy i is a weight vector and β yi is a bias term, respectively, both associated to the label y i .

Pairwise potentials: Following the restriction imposed by the inference approach we use, our pairwise potentials are obtained as follows:

ψ p (y i , y j , f i , f j ) = µ(y i , y j ) M m=1 w p (m) k (m) f (m) i , f (m) j (3)
where each k (m) is a Gaussian kernel over an arbitrary feature f (m) , w p (m) is a linear combination weight, and µ(y i , y j ) represents a label compatibility function. In general, our pairwise kernels have the following form:

k (m) f (m) i , f (m) j = exp   - |p i -p j | 2 2θ 2 p - |f (m) i -f (m) j | 2 2θ 2 (m)   (4)
where p i and p j are the coordinate vectors of pixels i and j. We include positions in order to increase the effect of close pixels over distant ones. Parameters θ p and θ (m) control the degree of relevance of the two parts of the kernels into the expression. The scale value θ (m) is estimated as the median of the distances over pairs (f

(m) i , f (m) j
) [START_REF] Schölkopf | Support Vector Learning[END_REF]. In the same way, θ p is obtained using f (m) = p. The compatibility function µ is given by the Potts model, µ(y i , y j ) = [y i = y j ]. It penalizes nearby similar pixels that are assigned to different labels.

3 Learning CRF's parameters using Structured Output SVM Our goal is to learn a vector w = (w u , w β , w p ), where w u , w β and w p are the weights for the unary features, for the bias term and for the pairwise kernels, respectively. We obtain it in a supervised way, using the 1-slack formulation of the Structured Output SVM with margin-rescaling presented in [START_REF] Joachims | Cutting-plane training of structural SVMs[END_REF].

Let the training set S = { s (1) , y (1) , ..., s (n) , y (n) }, where n is the number of training images. Each y (i) corresponds to the ground truth of the i-th image in the training set. Each set s (i) = {x (i) , B, f (i) } contains the set x (i) of unary feature vectors, a bias constant B, and the set f (i) of pairwise features for every pixel in the image.

The weights w are obtained by solving:

min w,ξ≥0 1 2 w 2 + Cξ (5)
s.t.∀ ȳ(1) , ..., ȳ(n) :

n i=1 w, ϕ(s (i) , y (i) ) -ϕ(s (i) , ȳ(i) ) ≥ n i=1 ∆(y (i) , ȳ(i) ) -ξ
where C is a regularization constant; ξ is a slack variable shared across all the constraints ȳ(i) ; ϕ(s, y) is a feature map function that relates a given set s with a given labeling y; and ∆(y, ȳ) is a loss function that evaluates the difference between a ground truth y and a constraint ȳ. In this work, we define ∆ as the typical Hamming loss:

∆(y, ȳ) = i [y i = ȳi ] (6)
Our feature map is defined as follows:

ϕ(s, y) =   k ϕ u (x k , y k ), k ϕ β (B, y k ), k j<k ϕ p (y k , y j , f k , f j )   (7)
where the components represent the sum of the unary feature map, the bias feature map and the pairwise feature map, respectively, for all the pixels in the image. We make precise the definitions of ϕ u , ϕ β , and ϕ p in the sequel. We define a binary vector ϕ y (y i ) ∈ {0, 1} |L| such that:

ϕ y (y i ) (1, 0) if y i = -1 (0, 1) if y i = 1 (8)
The individual feature maps are obtained as follows:

ϕ u (x k , y k ) = x k ⊗ ϕ y (y k ) (9)
ϕ β (B, y i ) = Bϕ y (y i ) (10)
∀m

: [ϕ p (y k , y j , f k , f j )] m = µ(y i , y j )k (m) (f (m) i , f (m) j ) ( 11
)
where ⊗ is the Kronecker product. Eq. ( 5) is solved using a cutting-plane approach [START_REF] Joachims | Cutting-plane training of structural SVMs[END_REF].

Validation & Results

We validate our method using the publicly available data set DRIVE [START_REF] Staal | Ridge based vessel segmentation in color images of the retina[END_REF]. DRIVE is widely utilized in the state-of-the-art to quantify blood vessel segmentation performance. The data set is divided into two sets: one for training and one for test, each of them containing 20 images. The test set provides two manual segmentations generated by two different experts for each image. selection of the first observer is accepted as ground truth and used for performance evaluation in literature. The training set includes a set of manual segmentations made by the first observer. We performed the evaluation of our method in terms of sensitivity (Se) and specificity (Sp):

Se = T P T P + F N , Sp = T N T N + F P , (12)
where T P , T N , F P and F N represent the number of true positives, true negatives, false positives and false negatives, respectively. Sensitivity and specificity are indicators of the capability of the algorithm to identify blood vessels and to exclude all other classes, respectively. The value of C and the combination of features utilized to compute the unary and pairwise energies of the CRF were selected without using test data. We have divided the original training set into two new subsets, training* and validation, containing 75% and 25% of the images, respectively. We utilized training* to train the model, and validation to optimize the combination of parameters. We then used the selected parameters to retrain on the training* set, and we evaluated the final model on the test set. During model selection, we used a criterion that minimized the distance with respect to the Se and Sp of the second human observer.

Several state-of-the-art features were considered, including wavelets [START_REF] Soares | Retinal vessel segmentation using the 2-d Gabor wavelet and supervised classification[END_REF], line [START_REF] Ricci | Retinal blood vessel segmentation using line operators and support vector classification[END_REF][START_REF] Nguyen | An effective retinal blood vessel segmentation method using multi-scale line detection[END_REF] and ridge detectors [START_REF] Martinez-Perez | Segmentation of blood vessels from red-free and fluorescein retinal images[END_REF], matched filter responses [START_REF] Al-Rawi | An improved matched filter for blood vessel detection of digital retinal images[END_REF], vessel enhancers [START_REF] Frangi | Multiscale vessel enhancement filtering[END_REF][START_REF] Marín | A new supervised method for blood vessel segmentation in retinal images by using gray-level and moment invariants-based features[END_REF][START_REF] Sinthanayothin | Automated localisation of the optic disc, fovea, and retinal blood vessels from digital colour fundus images[END_REF][START_REF] Saleh | An efficient algorithm for retinal blood vessel segmentation using h-maxima transform and multilevel thresholding[END_REF][START_REF] Zana | Segmentation of vessel-like patterns using mathematical morphology and curvature evaluation[END_REF], gradient module, Huffman's entropy and image intensities. A forward selection process was followed to select the best combination. All features are obtained over preprocessed images. We first select the green band of the original color image. Then we extend the borders of the FOV [START_REF] Soares | Retinal vessel segmentation using the 2-d Gabor wavelet and supervised classification[END_REF], and finally we subtract an estimated background for bias correction using a median filter with square windows of length 25. For each combination of features, values of C = 10 i , with i = {0, 1, ..., 5}, were explored. The best configuration we found for the unary features utilize 2-D Gabor wavelets [START_REF] Soares | Retinal vessel segmentation using the 2-d Gabor wavelet and supervised classification[END_REF], line detectors [START_REF] Ricci | Retinal blood vessel segmentation using line operators and support vector classification[END_REF] and the vessel enhancement process proposed in [START_REF] Sinthanayothin | Automated localisation of the optic disc, fovea, and retinal blood vessels from digital colour fundus images[END_REF]. The vessel enhanced image obtained after applying the process proposed in [START_REF] Saleh | An efficient algorithm for retinal blood vessel segmentation using h-maxima transform and multilevel thresholding[END_REF] is utilized as a pairwise feature. Two segmentation examples, obtained over a healthy and a pathological image, are included in Figure 1. We also include the ground truth labeling and the labeling performed by the second observer for qualitative comparison.

We have performed an extensive comparison of our method with state-of-the-art automated methods that report their results in terms of Se and Sp, and that are not trained on the test set. Figure 2 shows a scatter plot of the competing methods as well as the performance of the human observer, in terms of Se and Sp. We also include results obtained following the same principled discriminative training but using a localneighborhood based CRF. The numerical values are listed in Table 1.

Conclusions

In this work, we have presented a discriminatively trained segmentation model based on a fully connected CRF for the purpose of blood vessel segmentation in fundus images. Conditional random fields have not been previously applied in this setting, likely due to standard pairwise potentials (such as in a Potts model) leading to poor performance in the presence of thin, branching structures. In contrast, our approach is able to achieve high accuracy by considering a much more expressive, fully connected graph structure. The benefits of this approach are demonstrated empirically.

Extensive comparison with state-of-the-art methods has shown that our approach performed best and is a fully-automated segmentation algorithm that achieves results on par with a human annotator as measured by sensitivity and specificity. Our method is statistically tied with the expert annotator for both measures. This is in part due to previous studies focusing on raw pixel accuracy, which ignores the fact that the number of pixels occupied by blood vessels is a relatively small fraction of the image. As a result, competing methods suffer as measured by sensitivity. In contrast, our method closely matches the performance of an expert human both in sensitivity and specificity, and achieves a substantially higher sensitivity than all other competing Fig. 2. Scatter-plot of Se vs. Sp, comparing existing methods, local-neighborhood based CRF, the proposed method, and the human annotator, based on DRIVE data set. Our method is statistically tied with the performance of the expert annotator and achieves a much higher sensitivity than all other segmentation systems.

Method

Sensitivity Specificity Fraz et al., 2012 [START_REF] Fraz | Ensemble classification system applied for retinal vessel segmentation on child images containing various vessel profiles[END_REF] 0.7152 0.9769 You et al., 2011 [25] 0.741 0.9751 Miri and Mahloojifar, 2011 [START_REF] Miri | Retinal image analysis using curvelet transform and multistructure elements morphology by reconstruction[END_REF] 0.7352 0.9795 Zhang et al., 2010 [9] 0.712 0.9724 Palomera-Pérez et al., 2010 [START_REF] Palomera-Pérez | Parallel multiscale feature extraction and region growing: application in retinal blood vessel detection[END_REF] 0.64 0.967 Vlachos and Dermatas, 2009 [START_REF] Vlachos | Multi-scale retinal vessel segmentation using line tracking[END_REF] 0.747 0.955 Al-Diri et al., 2009 [START_REF] Al-Diri | An active contour model for segmenting and measuring retinal vessels[END_REF] 0.7282 0.9551 Espona et al., 2008 [28] 0.7436 0.9615 Espona et al., 2007 [29] 0 methods. Sensitivity is particularly important as it reflects an accurate estimation of the vessel pixels, the primary goal in vessel segmentation for fundus image analysis. Additional implementation details are available at http://pages.saclay.inria. fr/matthew.blaschko/projects/retina/.

Fig. 1 .

 1 Fig. 1. Examples of results obtained on healthy (top) and pathological (bottom) images from DRIVE. From left to right: original images, ground truth labelings, 2nd human observer labelings, our segmentations.

Table 1 .

 1 Comparison of Se and Sp of our method with respect to the existing blood vessel segmentation algorithms and the human observer, based on DRIVE data set.

		.6634	0.9682
	Martínez-Pérez et al., 2007 [8]	0.7246	0.9655
	Mendonca and Campilho, 2006 [7]	0.7344	0.9764
	Local-neighborhood based CRF 0.7042±0.049 0.9815±0.005
	Our method	0.785±0.045 0.967±0.007
	2nd human observer	0.776±0.059 0.973±0.008

Acknowledgements

This work is partially funded by ERC Grant 259112, and FP7-MC-CIG 334380.


Hanna Becker

Laurent Albera

Pierre Comon

Rémi Gribonval

Isabelle Merlet

Fast, variation-based methods for the analysis of extended brain sources

Keywords: EEG, extended source localization, ADMM, sparsity

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

FAST, VARIATION-BASED METHODS FOR THE ANALYSIS OF EXTENDED BRAIN SOURCES

Hanna Becker * , † , Laurent Albera ‡, §, ¶ , Pierre Comon † , Rémi Gribonval ¶ , Isabelle Merlet ‡, §

INTRODUCTION

Over the last decades, a large number of methods have been proposed to reconstruct the electric activity everywhere within the brain based on surface electroencephalographic (EEG) measurements (see e.g. [START_REF] Becker | A performance study of various brain source imaging approaches[END_REF] and references therein). This process is generally referred to as source imaging and permits to identify the brain regions which are involved in generating characteristic activity patterns such as epileptic spikes, which is of clinical importance. However, in order to delineate the seizure onset zone in epileptic patients, it is not only important to localize the foci of the epileptogenic sources, but also to identify their spatial extent. This is particularly challenging in the context of several highly correlated, simultaneously active source regions corresponding to epileptic activity that spreads from one region to another.

In a recent comparison of different source imaging algorithms [START_REF] Becker | A performance study of various brain source imaging approaches[END_REF], the VB-SCCD algorithm [START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF] showed a good performance for the localization of extended sources. Moreover, it permits to simultaneously localize several highly correlated active source regions, which is problematic with other extended source localization methods such as STWV-DA [START_REF] Becker | A performance study of various brain source imaging approaches[END_REF] and 4-ExSo-MUSIC [START_REF] Birot | Localisation of extended brain sources from EEG/MEG: the ExSo-MUSIC approach[END_REF]. Therefore, VB-SCCD is one of the most promising approaches for the identification of multiple active brain regions in the context of propagation phenomena. However, the algorithm shows some difficulties in separating close sources and tends to combine them into one large source. Furthermore, the implementation of VB-SCCD using Second Order Cone Programming (SOCP) [START_REF] Alizadeh | Second-order cone programming[END_REF] as proposed in [START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF] leads to a high computational complexity, which practically forbids the application of the method for large numbers of time samples.

In this paper, we improve on these points by proposing a new source imaging algorithm, subsequently referred to as sparse VB-SCCD (SVB-SCCD), which includes an additional L 1 -norm regularization term. Such an approach, also known as sparse Total Variation (sparse TV) regularization [START_REF] Baldassarre | Structured sparsity models for brain decoding from fMRI data[END_REF], TV-L 1 regularization [START_REF] Gramfort | Identifying predictive regions from fMRI with TV-l1 prior[END_REF] or fused LASSO [START_REF] Tibshirani | Sparsity and smoothness via the fused LASSO[END_REF], has previously been used in image processing [START_REF] Ma | An efficient algorithm for compressed mr imaging using total variation and wavelets[END_REF] and fMRI prediction [START_REF] Baldassarre | Structured sparsity models for brain decoding from fMRI data[END_REF][START_REF] Gramfort | Identifying predictive regions from fMRI with TV-l1 prior[END_REF], where it has been shown to lead to robust solutions, but is new in the field of brain source imaging. Note though that the combination of sparsity in the original source domain and in a transformed domain that is different from the total variation has been explored in [START_REF] Chang | Spatially sparse source cluster modeling by compressive neuromagnetic tomography[END_REF] for MEG source imaging. As shown in this paper, the SVB-SCCD method permits to obtain more focal source estimates than VB-SCCD and achieves the separation of even close sources. Furthermore, we demonstrate the use of a different optimization technique, called Alternating Direction Method of Multipliers (ADMM) [START_REF] Boyd | Distributed optimization and statistical learning via alternating direction method of multipliers[END_REF], which is much faster than SOCP. This gain on computational complexity enables us to apply the algorithm to large time intervals and to reconstruct the source signals. It also makes it possible to take into account the temporal structure of the data by employing an L 1,2 -norm regularization as first suggested in [START_REF] Ou | A distributed spatiotemporal EEG/MEG inverse solver[END_REF], leading to more robust source estimates. The superior performance of the resulting L 1,2 -SVB-SCCD algorithm in comparison to VB-SCCD, STWV-DA, 4-ExSo-MUSIC, and cLORETA [START_REF] Wagner | Smooth reconstruction of cortical sources from EEG and MEG recordings[END_REF] is demonstrated by means of computer simulations, which are conducted in the context of a realistic head model and highly correlated extended sources emitting epileptic spike activity.

DATA MODEL AND PROBLEM FORMULATION

The electric potential at the surface of the scalp is characterized by the superposition of signals originating from all over the brain. As most of these signals are generated by pyramidal cells located in the gray matter, for modeling purposes, we define a source space that consists of D dipoles located on the cortical surface with a fixed orientation perpendicular to this surface (see also [START_REF] Dale | Improved localization of cortical activity by combining EEG and MEG with MRI cortical surface reconstruction: a linear approach[END_REF]). More particularly, the dipoles are positioned at the centroids of the triangles of a cortical surface mesh. The EEG measurements X ∈ R N ×T recorded with N sensors for T time samples can then be modeled as a weighted sum of the dipole signals. The weights depend on the propagation effects within the head volume conductor and are summarized in the so-called lead field matrix G ∈ R N ×D , which can be computed numerically using boundary element methods (BEM) [START_REF] Gramfort | Mapping, timing and tracking cortical activations with MEG and EEG: Methods and application to human vision[END_REF]. Distinguishing two types of cerebral activity, the signals of interest (e.g., epileptic activity) and the background activity of the brain, characterized by the matrices S ∈ R D×T and S b ∈ R D×T , respectively, this leads to the following data model:

X = GS + GS b = GS + N. (1)
The objective of source imaging algorithms consists in estimating the signal matrix S from the measurements X. As the number of dipoles (D ∼ 10000) is generally much larger than the number of sensors (N ∼ 100), this is an ill-posed problem, which requires additional assumptions to regularize the solution. In this paper, we consider several approaches based on sparsity, which are described in the next section.

SOURCE LOCALIZATION AND EXTRACTION

VB-SCCD

The VB-SCCD algorithm [START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF] assumes a piece-wise constant spatial source distribution, which is achieved by imposing sparsity on the variational map of the sources. The latter can be computed by applying a linear transform, characterized by the matrix V, to the source distribution, which is equivalent to computing the total variation on the discretized cortical surface. The elements V p,d of V, p = 1, . . . , P , d = 1, . . . , D, where P is the number of edges of the triangular grid, are given by:

V p,d =    1 if d = d p,1 -1 if d = d p,2 0 otherwise (2)
where d p,1 and d p,2 are the indices of the dipoles sharing the p-th edge. Thus, VS describes the differences in amplitude between adjacent dipoles. This gives rise to the VB-SCCD cost function, which is minimized with respect to the signal matrix S1 :

min S 1 2 ||X -GS|| 2 F + λ||VS|| 1 . (3)
The regularization parameter λ balances between the reconstruction error and the constraint, corresponding to the first and second term in (3), respectively. As suggested in [START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF], this parameter may be adjusted according to the acceptable upper limit for the reconstruction error, which depends on the noise level.

SVB-SCCD

In practice, it is reasonable to assume that only a small number of the source dipoles contribute to the signals of interest. Hence, we introduce an additional regularization term in (3) that imposes sparsity in the original source domain, leading to the following optimization problem which is equivalent to the sparse TV [START_REF] Baldassarre | Structured sparsity models for brain decoding from fMRI data[END_REF] or fused LASSO [START_REF] Tibshirani | Sparsity and smoothness via the fused LASSO[END_REF] approach:

min S 1 2 ||X -GS|| 2 F + λ(||VS|| 1 + α||S|| 1 ). (4)
This approach does not only permit us to adjust the size of the reconstructed source regions by varying the new regularization parameter α, but also prevents the estimated signal vector from being amplitude-biased, which is a problem that frequently arises using the VB-SCCD algorithm. Setting α = 1 leads to very focal source estimates, whereas small α only avoid the amplitude bias, but do not influence the size of the reconstructed source regions. In our experience, reasonable results can be achieved for 0.01 ≤ α ≤ 1. As for VB-SCCD, the regularization parameter λ regulates the impact of the source priors and is tuned depending on the noise level.

Exploitation of temporal structure

The VB-SCCD and SVB-SCCD algorithms as described in the previous sections consider each time sample indepen-dently and thus do not take into account the temporal structure of the data. However, it can be expected that in the considered time interval, the active source regions stay the same. This hypothesis can be enforced by replacing the L 1 -norm in the above equations by the L 1,2 -norm, which is defined as follows:

||S|| 1,2 = D d=1 T t=1 S 2 d,t
. This permits to obtain more robust source estimates. The resulting source localization approach is subsequently called L 1,2 -SVB-SCCD.

Optimization using ADMM

The optimization problems of the three different algorithms, VB-SCCD, SVB-SCCD, and L 1,2 -SVB-SCCD, can be rewritten in a generalized, constrained optimization framework with latent variables Y and Z:

min S 1 2 ||X -GS|| 2 F + λ(f (Y) + αf (Z)) s. t. Y = VS, Z = S. (5)
Here, f represents the regularizing function that is either the L 1 norm (for SVB-SCCD) or the L 1,2 -norm (for L 1,2 -SVB-SCCD). Note that the VB-SCCD optimization problem can be regarded as a special case of ( 5) where α = 0. Problem ( 5) can be solved using ADMM [START_REF] Boyd | Distributed optimization and statistical learning via alternating direction method of multipliers[END_REF], which is a simple and efficient algorithm for constrained convex optimization. It is based on the idea of alternatingly updating the variables S ∈ R D×T , Y ∈ R P ×T , and Z ∈ R D×T in the augmented Lagrangian of ( 5), as well as computing alternating updates of the scaled Lagrangian multipliers U ∈ R P ×T and W ∈ R D×T . After initialization (for example, by setting all variables to zero), at the k-th iteration, the following update rules can be derived:

S (k+1) = P -1 G T X + ρV T (Y (k) -U (k) ) + ρ(Z (k) -W (k) ) with P = G T G + ρ(V T V + I) Y (k+1) = prox f,λ/ρ VS (k+1) + U (k) Z (k+1) = prox f,λα/ρ S (k+1) + W (k) U (k+1) = U (k) + VS (k+1) -Y (k+1) W (k+1) = W (k) + S (k+1) -Z (k+1)
where ρ > 0 denotes the penalty parameter introduced in the augmented Lagrangian (see [START_REF] Boyd | Distributed optimization and statistical learning via alternating direction method of multipliers[END_REF]). Please note that in practice, the computation of the inverse of the large matrix P ∈ R D×D should be avoided, for example, by resorting to inversion lemma and matrix decompositions (such as the QRdecomposition) which can be computed efficiently. The updates of Y and Z are formulated using the proximity operator, which is given by:

prox f,β (Y) = arg min X 1 2 ||Y -X|| 2 F + βf (X). (6)
Solutions to [START_REF] Gramfort | Identifying predictive regions from fMRI with TV-l1 prior[END_REF] for f corresponding to the L 1 -norm or the L 1,2 -norm of X can be found in [START_REF] Gramfort | Mixednorm estimates for the M/EEG inverse problem using accelerated gradient methods[END_REF]. The algorithm is stopped after convergence or a maximal number of iterations is reached.

SIMULATIONS

In this section, we compare the performance of SVB-SCCD, VB-SCCD, L 1,2 -SVB-SCCD, L 1,2 -VB-SCCD, STWV-DA, 4-ExSo-MUSIC, and cLORETA based on computer simulations. To this end, EEG data is generated for N = 91 electrodes using a realistic head model with three compartments that represent the brain, the skull, and the scalp. The source space consists of D = 19626 dipoles corresponding to the triangles of the cortical surface mesh with orientations perpendicular to the cortical surface. A BEM method2 is used to compute the lead field matrix. We considered three source regions, in the following referred to as patches, each of which is composed of 100 adjacent dipoles. The first patch was attributed an epileptic spike signal comprising T = 200 time samples (at 256 Hz sampling frequency) that was segmented from stereotactic EEG (SEEG) recordings of a patient suffering from epilepsy. We then generated 100 different realizations of this signal, one for each patch dipole, by introducing small variations in amplitude and delay. Assuming that the other patches were activated due to a propagation of the epileptic activity of the first patch, we used the same signals for the dipoles of the second and third patch, but introduced a delay of 4 to 24 ms depending on the distance to the first patch. All source dipoles that do not belong to a patch were attributed Gaussian background activity with an amplitude that was adjusted to the amplitude of the SEEG signals between epileptic spikes, thus leading to realistic Signal to Noise Ratio (SNR) such that ||GS|| 2 F /||N|| 2 F ≈ 1. The EEG data were spatially prewhitened before applying the source localization algorithms. For both VB-SCCD and SVB-SCCD, the regularization parameter λ was adjusted such that the reconstruction error lies within a confidence interval of 95 to 99 % of the noise power. In the case of SVB-SCCD, we used a fixed parameter α = 0.67 because we found that this leads to reasonable results for the considered scenarios. For VB-SCCD, SVB-SCCD, and cLORETA, which provide one source estimate per time sample, we determined the active patches by thresholding the source estimates at the data sample of maximal power, corresponding to the maximum of the epileptic spike. For each identified source region, comprised of adjacent dipoles, we then computed the average of the time signals of all involved source dipoles in order to obtain one estimated time signal per patch. For 4-ExSo-MUSIC, an estimate of the patch signals was computed as Ŝp = H + X. Here, H + denotes the pseudoinverse of the spatial mixing matrix H ∈ R N ×R whose r-th column corresponds to the sum of the lead field vectors associated to the dipoles belonging to the r-th estimated patch. No further processing was necessary in case of STWV-DA, as this algorithm already provides a time signal for each estimated extended source at its output.

The performance of the source localization was assessed ROC curves for (left) scenario 1 (SupFr, InfFr, SupOcc), (middle) scenario 2 (SupOcc, InfPa, InfFr), and (right) scenario 3 (SupOcc, MidTe, OccTe). using the Dipole Localization Error (DLE) and the Receiver Operating Characteristic (ROC) curves, which evaluate the True Positive Fraction (TPF) as a function of the False Positive Fraction (FPF). If I and Î denote the original and estimated sets of indices of all dipoles belonging to an active patch, J represents the set of all dipoles belonging to the source space, Q and Q are the numbers of original and estimated active dipoles, and r k denotes the position of the k-th source dipole, then the DLE, TPF and FPF are defined as:

DLE = 1 2Q k∈I min ℓ∈ Î ||r k -r ℓ || + 1 2 Q ℓ∈ Î min k∈I ||r k -r ℓ || TPF = #(I ∩ Î) #I ; FPF = # Î -#(I ∩ Î) #J -#I .
We use the notation #I for the cardinality of the set I.

The quality of the extracted signals is evaluated by calculating the correlation coefficients between the estimated patch signal and the averaged signal of all dipoles belonging to a patch. We then computed the mean of the correlation coefficients for all patches.

We simulated three different scenarios with patches of different distances. The first scenario comprised three patches of medium to large distance, located in the superior frontal gyrus (patch SupFr), the inferior frontal gyrus (patch InfFr), and the superior occipital gyrus (SupOcc). The second scenario included two close patches, positioned in the superior occipital and the inferior parietal gyri (patches SupOcc and InfPa), with the third patch in the inferior frontal gyrus (patch InfFr). In the third scenario, we considered three close patches, located in the inferior parietal gyrus (patch InfPa), the mid temporal gyrus (patch MidTe), and the occipital temporal gyrus (patch OccTe). The performance achieved with the different source imaging algorithms in terms of DLE and signal correlation coefficient for the three scenarios, averaged over 30 realizations with different patch signals and background activity, is summarized in localization approaches for the considered multi-patch scenarios. The use of the additional L 1 -norm regularization term in the SVB-SCCD approach turns out to be insignificant in the case of three patches with medium distance (scenario SupFr, InfFr, SupOcc), as SVB-SCCD and VB-SCCD exhibited a comparable performance in this case. However, for two close patches (scenario SupOcc, InfPa, InfFr), one can observe a slight improvement of the DLE obtained with SVB-SCCD compared to VB-SCCD, and for three close patches, the SVB-SCCD approach clearly leads to better results than VB-SCCD. This can also be seen in Figure 2, where we illustrate an example of the source imaging results obtained with the different methods for the scenario with the three close patches InfPa, MidTe, and OccTe. Obviously, the SVB-SCCD approach provides a better separation of the sources than the VB-SCCD approach. Furthermore, Figure 2 visualizes the difficulties encountered with the other extended source localization methods. While the cLORETA solution exhibits high dipole amplitudes at the three patch locations, it does not correctly recover the patches' extents. STWV-DA finds only two patches, one of which is dislocated and deformed compared to the original patch, whereas 4-ExSo-MUSIC identifies only one patch of overestimated extent. The exploitation of the temporal structure of the data in the VB-SCCD and SVB-SCCD algorithms hardly has an impact on the source localization results of scenarios 1 and 2, but for SVB-SCCD, it yields more robust solutions in the case of three close patches (cf. DLE for scenario 3). Furthermore, it leads to a better performance in terms of source extraction as demonstrated by the obtained signal correlation coefficients. For fixed parameters, the CPU runtimes of the different algorithms are roughly comparable, except for cLORETA, which is faster 3 .

CONCLUSIONS

In this paper, we have analyzed two extensions of the VB-SCCD algorithm. Following the fused LASSO approach, we have included an additional, sparsity-inducing regularization term, which permits to obtain a better separation of close sources. Furthermore, we have taken into account the temporal structure of the data, which leads to an increased performance in terms of signal extraction. Finally, we have illustrated the use of an efficient algorithm, ADMM, to solve the L 1,2 -SVB-SCCD optimization problem in a much faster way 3 Please note that the CPU times required for the constructions of the Lapacian matrix for cLORETA and of a dictionary of potential sources for STWV-DA and ExSo-MUSIC have not been considered in this analysis. than the previously employed SOCP algorithm. The superior performance of the proposed approach in comparison to the classic VB-SCCD algorithm as well as other state-of-the-art methods for extended source localization has been demonstrated by means of realistic computer simulations.

H.

  Becker was supported by Conseil Régional PACA and by CNRS. The work of P. Comon was funded by the FP7 European Research Council Programme, DECODA project, under grant ERC-AdG-2013-320594. The work of R. Gribonval was funded by the FP7 European Research Council Programme, PLEASE project, under grant ERC-StG-2011-277906. We also acknowledge the support of Programme ANR 2010 BLAN 0309 01 (project MULTIMODEL).

  Fig.1. ROC curves for (left) scenario 1 (SupFr, InfFr, SupOcc), (middle) scenario 2 (SupOcc, InfPa, InfFr), and (right) scenario 3 (SupOcc, MidTe, OccTe).

Fig. 2 .

 2 Fig. 2. Source imaging results obtained with the different tested algorithms.

Table 1 .

 1 The corresponding ROC curves are shown in Figure 1. Both the ROC curves and the DLE values show that the VB-SCCD and SVB-SCCD type algorithms clearly outperform the other extended source

			DLE in cm		corr. coeff. in %
	scenario	1	2	3	1	2	3
	VB-SCCD	0.99 2.57 9.73 94.9 92.5 78.5
	SVB-SCCD	0.94 1.05 3.81 95.5 94.9 89.3
	L1,2-VB-SCCD	0.97 1.09 10.8 97.9 97.5 77.9
	L1,2-SVB-SCCD 1.03 1.06 2.24 98.5 98.3 96.6
	STWV-DA	10.4 12.3 18.7 65.1 80.1 83.6
	4-ExSo-MUSIC	35.3 9.78 9.74 81.0 80.4 74.8
	cLORETA	8.69 4.35 13.1 82.5 90.4 45.7

Table 1 .

 1 Performance

of source imaging algorithms in terms of DLE and signal correlation for scenario 1 (SupFr, InfFr, SupOcc), scenario 2 (SupOcc, InfPa, InfFr), and scenario 3 (SupOcc, MidTe, OccTe).

Please note that we use a different formulation of the optimization problem than employed in[START_REF] Ding | Reconstructing cortical current density by exploring sparseness in the transform domain[END_REF], but which leads to equivalent results for an appropriate value of the regularization parameter λ. Furthermore, the optimization problem was originally stated for only one time sample.

ASA, ANT, Enschede, Netherlands