Fayrouz Soualah-Alila
email: fayrouz.soualah-alila@univ-lr.fr

Cyril Faucher
email: cyril.faucher@univ-lr.fr

Frédéric Bertrand
email: fbertran@univ-lr.fr

Mickaël Coustaty
email: mickael.coustaty@univ-lr.fr

Antoine Doucet
email: antoine.doucet@univ-lr.fr

Applying Semantic Web Technologies for Improving the Visibility of Tourism Data

Keywords: E.2 [Data]: Data Storage Representations; D.2.12 [Software]: Software Engineering-Interoperability Tourism, Semantic Web Tourinflux, TourInFrance, Ontology, Schema.org, Mapping

Tourism industry is an extremely information-intensive, complex and dynamic activity. It can benefit from semantic Web technologies, due to the significant heterogeneity of information sources and the high volume of on-line data. The management of semantically diverse annotated tourism data is facilitated by ontologies that provide methods and standards, which allow flexibility and more intelligent access to on-line data. This paper provides a description of some of the early results of the Tourinflux project which aims to apply semantic Web technologies to support tourist actors in effectively finding and publishing information on the Web.

INTRODUCTION 1.1 Context

The tourism is viewed as information intensive industry and a highly dynamic and changing domain where information plays an important role for decision and action making [START_REF] Inkpen | Information technology for travel and tourism[END_REF]. Tourism represents a significant part of the global economy, more than 9% of the global GDP (US $7 trillion) [START_REF] Wttc | World travel and tourism council, travel and tourism, economic impact 2014 world[END_REF]. It becomes increasingly difficult for policy makers, on the one hand to guide choices in order to boost and make their territory attractive, and on the other hand to analyze the comments and the vision of the territory by users. As the Permission to make digital or hard copies of all or part of this work for personal or classroom use is granted without fee provided that copies are not made or distributed for profit or commercial advantage and that copies bear this notice and the full citation on the first page. Copyrights for components of this work owned by others than ACM must be honored. Abstracting with credit is permitted. To copy otherwise, or republish, to post on servers or to redistribute to lists, requires prior specific permission and/or a fee. Request permissions from Permissions@acm.org. World Wide Web has changed people's daily life, it has significantly influenced the way of information gathering and exchanging in the tourism area, with the intensive use of social networks and Web sites specialized in e-tourism (Tri-pAdvisor, Booking.com, etc.). Along with the Web technologies, tourism actors such as tour operators need knowledge about travel, hotels, markets, etc. Indeed, information like hotel prices, hotel localization, offers, market tendencies are used to improve their service quality by enhancing employee's knowledge about customer's preferences. Other actors such as experts from tourism industry are also making use of the Web to promote their territory, and so to attract more visitors. However, the gaps of the existing Web technologies arise two main topics: First, tourism domain is characterized by a significant heterogeneity of information sources and by a high volume of on-line data. Data related to tourism is produced by different experts (travel agents, tourist offices, etc.) and by visitors, constituting semantically heterogeneous data, often incomplete and inconsistent. This data could be composed of information related to tourism objects (hotels, concert, restaurant, etc.), temporal information or opinions. There already exist different taxonomies and catalogues which are designed and used internally by tourism actors to help them to manage heterogeneous tourism data. Efforts are now made to generate standards to facilitate inter and intra tourism data exchange. Second, the information contained on Web pages are originally designed to be human-readable, and so, most of information currently available on the Web are kept in large collections of textual documents. As the Web grows in size and complexity, there is an increasing need for automating some of the time consuming tasks such as information searching, extraction and interpretation. Semantic Web technologies are an answer to this issue by proposing some novel and sophisticated question answering systems. In this way, ontologies provides a formal framework to organize data and to browse, search or access this information [START_REF] Bruijn | Modeling semantic web services[END_REF].

Tourinflux

This work takes place under the Tourinflux 1 project which aims at providing the tourism industry with a set of tools (1) allowing them to handle both their internal data, and the information available on the Web, and (2) allowing to improve the display information available about their terri-tory on the Web. This should allow them to improve their decision process and subsequently, their effectiveness. This paper presents some preliminary results of Tourinflux project. The paper aim to describe the development of an ontology based model for tourism domain. It is aimed primarily to present some aspects of knowledge management: Section 2 reviews the related ontologies developed in the previous projects. In particular, we describe how an ontology can be used to improve the quality of users Web search and data publication. Section 2 also provides description of current tourism standards and a brief background on a particular ontology: the Schema.org 2 model. In an early stage of our project, a partial ontology called TIFSem (Semantic TourInFrance) for the tourism was created, using Protégé 3 and the Web Ontology Language 4 (OWL). A partial view of TIFSem ontology developed will be presented in section 3 with some example of queries. We notice that this work is in progress; we are still gathering new concepts for its taxonomy and new axioms in order to validate them with experts of the domain. Finally, conclusions and outlook are provided in section 4.

STATE OF ART

Tourism ontologies

In the domain of tourism, several researches have focused on the design of ontologies. Several available tourism ontologies show the current status of the efforts: The OTA (OpenTravelAlliance) specification defines XML Schema documents corresponding to events and activities in various travel sectors [START_REF] Ota | opentravel alliance message specifications[END_REF]. The Harmonise ontology is specialised to address interoperability problems focusing on data exchange [START_REF] Dell'erba | Simplified composition overview of the granule Geolocation with Schema.org mapping Figure 5: Simplified composition overview of the granule Geolocation with Schema[END_REF]. Harmonise is based on mapping different tourism ontologies by using a mediating ontology. This central Harmonise ontology is represented in RDF and focuses mainly on tourism events and accommodation types. The Hi-Touch ontology models tourism destinations and their associated documentations [START_REF] Legrand | Semantic web methodologies and tools for intra-european sustainable tourism[END_REF]. The ontology classifies tourist objects, which are linked together in a network by semantic relationships. The top-level classes of the Hi-Touch ontology are documents (any kind of documentation about a tourism product), objects (the tourism objects themselves) and publication (a document created from the results of a query). The QALL-ME ontology provides a model to describe tourism destinations, sites, events and transportation [START_REF] Ou | Development and alignment of a domain-specific ontology for question answering[END_REF] and aims at establishing a shared infrastructure for multilingual and multimodal question answering. The Tourpedia catalogue describes places (points, restaurants, accommodations and interests) related to different locations (Amsterdam, London, Paris, Rome, etc.) [START_REF] Cresci | Towards a dbpedia of tourism: the case of tourpedia[END_REF]. These models focus on different areas of tourism domain but there does not exist one single ontology which matches all the needs of different tourism related applications. We then propose in section 3 an ontology to globally describe tourism information mixing heterogeneous content. Concepts included in the defined ontology will allow us to annotate information sources on tourism. Then the enriched information can be used: (1) from the user side, to match 2 http://schema.org/ 3 http://protege.stanford.edu/ 4 http://www.w3.org/2001/sw/wiki/OWL tailored package holidays to client preferences and (2) from tourism experts' side, to analyze and better manage on-line data about their territory.

Tourism standards

The interoperability of Tourism Information Systems (TIS) is a major challenge for the development of tourism domain. Several national, European and international institutional initiatives have proposed different standards to meet the specific needs of tourism professionals. In France, professionals have adopted the TourInFrance standard (TIF). Since its creation in 1999 by the tourism ministry, TourInFrance is used today by more than 3000 tourist offices in France, by Departmental Tourist Committees (DTC) and by different tour operators, to facilitate data exchange between these different actors. In 2004, the TourInFrance Technical Group (TIFTG) approved the new version of the standard, TIF V3. In this version, the standard has evolved towards XML technologies to facilitate publishing information on the Web and information exchange between systems, and is accompanied by several thesaurus. The last version of TIF can be downloaded from the following Web site: http://tourinfrance.free.fr/. Since 2005, this standard stopped evolving. As a result, tourism professionals have adapted the standard to their own needs and proposed their own evolution in an unorganized way, producing TIS that are not really interoperable among themselves and that cannot directly be shared using international standard. Accordingly, the exploitation of tourism information is trapped in its own territory, and so it is impossible to aggregate these information.

To illustrate this evolution of TIF we show in the following example how the standard was transformed. The figures 1 and 2 show different syntaxes of TIF to describe an address; data represented in the figure 1 respect the standard TIF V3, while data represented in the figure 2 is structured under a derived version of TIF (different tags syntaxes, new tags added, etc.).

To overcome these limitations we want to evolve the TIF standard to share the knowledge it represents and to ensure data interoperability, by applying the concept of ontology to represent the standard terminology [START_REF] Bittner | Ontology and semantic interoperability[END_REF].

Semantic Web standards

TIF standard is unable to easily share and interoperate with global Web standards. Since few years, an initiative led by Bing, Google and Yahoo tends to become a de-facto norm to easily share semantic content. Launched in 2011, Schema.org aims to create and to support a common set of schemas for structured data markup Web pages. The purpose of Schema.org is to provide a collection of schemas (HTML tags) that webmasters can use to markup HTML pages in ways recognized by major search providers, and that can also be used for structured data interoperability (RDFa, JSON-LD, etc.). Upon using these tags in a website, search engines can understand the intention of using resource (text, image, video) in that website better; which in turn would return better search results when a user is looking for a specific resource through a search engine [START_REF] Toma | Improving the online visibility of touristic service providers by using semantic annotations[END_REF]. Schema.org is a very large vocabulary counting hundreds of terms from multiple domains (the full specification of Schema.org is available at https://schema.org/docs/full.html).

In our approach we want to spread enriched semantic tourist data that can be easily indexed by search engines. One solution would be to use directly the Schema.org ontology which provides good indexing by search engines and which is easy to implement. The disadvantages of this solution is that we lose precision in the tourist data of the TIF model, in addition, using directly Schema.org can cause economic problems because a large number of tourist offices already use TIF or derived versions. A second solution would be to realize an ontology by matching terms of TIF with terms of Schema.org by using OWL relations, and work with Schema.org community to extend the schema, either formally by adding new terms or informally by defining how Schema.org can be combined with some additional vocabulary terms. We can thus produce more accurate information than the terms of Schema.org, in addition with a TIF/Schema.org correspondence, producing an interoperable model becomes possible from XML TIF sources. This method is more complex to implement because it requires to connect each term of TIF with a terms of Schema.org. The difficulty is not the transformation of XML TIF to OWL, but searching correspondence between the terms of TIF and those of Schema.org, and seeking what is lacking in each model. We suggest to consider the second solution, which presents a better compromise. We consider the Schema.org as the extension core of the TIF standard to improve the online visibility of touristic service providers. Schema.org includes many of the TIF terms. Of course not all of the Schema.org terms are relevant for the tourism domain. The relevant terms are those that belong to the categories Hotels, Food Establishments, Events, Trips, Place of Interest and News. These tourism terms are distributed in different trees of the schema. For example:

• Thing → Place → LocalBusiness → LodgingBusiness:

-Hostel, Hotel, Motel, etc. • Thing → Event:

-MusicEvent, SocialEvent, SportsEvent, etc. Schema.org offers other interesting terms such as customer reviews and reservation:

•

Thing → Action → AssessAction → ReviewAction • Thing → Intangible → Reservation:
-EventReservation, FoodEstablishmentReservation, LodgingReservation, etc. In addition, several Linked data initiatives have enabled the development of mappings between Schema.org and unavoidable semantic resources such as DBPedia, Dublin Core, FOAF, GoodRelations, SIOC and WordNet (full mappings of Web data vocabularies to Schema.org terms are available at http://schema.rdfs.org/mappings.html). We are interested in establishing the mappings at two levels. One will be a mapping between concepts which have the same name or approximate name, and between concepts which have the same role. Then taking in to account the mapped concepts (equivalentClass, subClassOf), we need to match one property from the Schema.org term with another property from a concept in TIF (equivalentProperty, subPropertyOf).

PROPOSED APPROACH

The tourist information can be used in different areas and different contexts. It is important that this information have to be modular, so that the elements describing the same entity can be used separately, hence the notion of a modular information represented by the IO. An InformationObject (IO) is a modular entity, that can be used and re-used to support tourism activities. An IO can be a hotel, a restaurant, an event, etc., or any other resource that can be attached to a tourism activity. The idea of IO is to create an entity that is:

• interoperable: can plug-and-play with any TIS,

• reusable: can be used or adapted for different use cases,

• accessible: can be stored a way that allows for easy searchability, Table 1: Core concepts of the TIFSem ontology CORE CONCEPTS Description

DUBLIN CORE

The Dublin Core is a set of vocabulary terms that can be used to describe Web resources (video, images, Web pages, etc.).

UPDATE

Information about the updates of a resource, and information used to represent aspects of the evolution of the resource.

MULTIMEDIA

Information about multimedia materials associated to the described resource.

CONTACTS

Information about ways to be in relation with a resource, or a person responsible for this resource.

LEGAL INFORMATION

Information to formally identify an entity and its activities.

CLASSIFICATIONS

Information used to qualify the resource described and identify the quality of services associated.

RELATED SERVICES

Information to reference other resources.

GEOLOCATIONS

Information to locate and describe the environment a resource.

PERIODS

Information about an opening period, a closing period, a reservation period, etc.

CUSTOMERS Customer information.

LANGUAGES Languages spoken in the described resource.

RESERVATION MODES Information about the reservation for the described resource, namely, who to contact customers and if the reservation is required.

PRICES

Services' price of the resource described, as well as the different means of payment.

CAPACITY

The capacity of the resource.

OFFERS SERVICES Information about the services available within the described resource or nearby.

ADDITIONAL DESCRIPTION Additional information about the described resource.

ITINERARIES

Activities associated to the described resource such as hiking.

SCHEDULES

Information about the availability status of services, over defined periods.

• manageable: can be tracked and updated over time.

In TIF, a typical IO is composed by 18 parts called granules (or semantic units). These granules form the main core of the TIFSem ontology. We introduce then the TIFSem ontology for the tourist domain, which contains concepts and relations for tourist resources. It consists of 19 main concepts of tourist information described by the table 1. The main concepts of that ontology have sub-concepts pertinent to the local area, and for all sub-concepts they have is a relationships. Those sub-concepts have properties, restrictions and semantic relations among other concepts.

The division of these main concepts into sub-concepts, was also guided by different data sources. Indeed, in order to elaborate the TIFSem ontology, we have consulted different kinds of sources to enable the understanding and the collecting of concepts related to the specialized domain of tourism, and the corresponding vocabulary. Sources coming from Departmental Tourism Committee of the Charente-Maritime5 (CDT17) and Departmental Tourism Committee of the Aube6 (CDT10) were consulted.

After collecting the concepts and corresponding lexical items from the sources, we started structuring the first version of the ontology, selecting classes and subclasses (step (1), figure 4). The ontology is currently available at this link: http://flux.univ-lr.fr/public/DataTourism.owl. The mapping with Schema.org (step (2), figure 4) complements the ontology to make it more complete, up to date and coherent. Table 2 shows the mappings between the main concepts coming from TIFSem to those in Schema.org. 6 show a simplified composition overview of the two granules Geolocation and ReservationModes, with Schema.org mapping.

The meaning of concepts is mainly expressed through the definitions in natural language. So we used OWL language which is equipped with a rich expressive power, also to enrich Schema.org which is specified in RDF. The use of an upper ontology, especially when there is extensive use of multiple inheritance, has the advantage of providing a manual navigation of the ontology and therefore simplify the recovery of concepts that can serve as a description. The last step in the development process of TIFSem is the integration of instances into the ontology coming from CDT17 and CDT10 (step (3) (4), figure 4). To illustrate our approach, we consider a first example that ranks different accommodation offers based on their proximity to restaurants, bars and events in La Rochelle. A couple wants to spend a week on vacation may choose between hotels not far from downtown. The couple likes dining out, going to bars and enjoying events. Therefore, an additional characteristic that quantifies the aptness of each hotel for those that like to go out, would be of great help. However, such qualitative information is rarely available on tourism on-line Web sites. One possibility would be to compute such semantic annotations based on the geolocalisation of the hotels and those objects related to restaurants, bars and events. Figure 7 shows a simplified SPARQL query reflecting the above example. Distance function filters any potential results whose corresponding return value was not less than a constant value. As second example, the tourist office of La Rochelle wants to have a better visibility of its territory by extracting information about rural events occurring on their territory. They need to know the events audience and information about the public attending these events. Figure 8 shows a simplified SPARQL query reflecting the second example.

CONCLUSIONS

In this paper, we presented some early stage work in the Tourinflux project on identifying tourism domain ontology. Semantic technologies provide methods and concepts facilitating effective integration of tourism information originating from various sources, representing basic notions and conceptual relations in tourism domain. The objective is to improve the display of search results, making it easier for tourism industrial to find data in the right Web pages. In particular having annotations on websites that can be understood by search engines busts the on-line visibility and increases the chances that Web sites are in the search engines results to a relevant query. As part of our current and future work we are developing a Schema.org mapping to facilitate Web searching of tourism data. We are also in the process of extending the TIFSem ontology by collecting contents about more touristic service providers. We also plan to expand the TIFSem model with other types of data (about opinions, temporal data), to provide semantic and contextual answers to queries.

ESAIR' 15 ,

 15 October 23 2015, Melbourne, VIC, Australia c 2015 ACM. ISBN 978-1-4503-3790-8/15/10 ...$15.00 DOI: http://dx.doi.org/10.1145/2810133.2810137.

Figure 1 :

 1 Figure 1: Data respecting the standard TIF V3

Figure 2 :

 2 Figure 2: Data structured under a derived version of TIF

Figure 3 :

 3 Figure 3: Schema.org/TIF mapping

Table 2 :

 2 TIFSem main concepts/Schemas.

		org
	CONCEPT	Schemas.org
	MULTIMEDIA	schema.org:MediaObject
	CLASSIFICATIONS	schema.org:Rating
	CONTACTS	schema.org:ContactPoint
	LEGAL INFORMATION schema.org:Organization
	LANGUAGES	schema.org:Language
	GEOLOCATION	schema.org:Place
	RESERVATION MODES schema.org:Reservation,
		schema.org:LodgingReservation
	PRICES	schema.org:Offer,
		schema.org:PriceSpecification
	Figures 5 and figure

http://tourinflux.univ-lr.fr/

http://www.charente-maritime.org/

http://www.aube-champagne.com/


Eric De Bodt

Marie Cottrell

Patrick Letremy

Michel Verleysen

On the use of self-organizing maps to accelerate vector quantization

Self-organizing maps (SOM) are widely used for their topology preservation property: neighboring input vectors are quantified (or classified) either on the same location or on neighbor ones on a predefined grid. SOM are also widely used for their more classical vector quantization property. We show in this paper that using SOM instead of the more classical Simple Competitive Learning (SCL) algorithm drastically increases the speed of convergence of the vector quantization process. This fact is demonstrated through extensive simulations on artificial and real examples, with specific SOM (fixed and decreasing neighborhoods) and SCL algorithms.

Motivation

Vector quantization (VQ) is a widely used tool in many data analysis' fields. It consists in replacing a continuous distribution by a finite set of quantizers, while minimizing a predefined distortion criterion. Vector quantization may be used in clustering or classification tasks, where the aim is to determine groups (clusters) of data sharing common properties. It can also be used in data compression, where the aim is to replace the initial data by a finite set of quantified ones; labeling the quantified set and using the labels rather than the data themselves makes compression possible. Vector quantization is basically an unsupervised process. Supervised variants exist (LVQ1, LVQ2, in Kohonen [START_REF] Kohonen | Comparison of SOM Point Densities Based on Different Criteria[END_REF]); in this last case, the distortion criterion takes class labels into account.

While most recent learning algorithms in various domains do not consider anymore the limitations due to computational load because of the increasing power of computers, this question is still important in VQ. The reason is that VQ is often used on huge databases in high-dimensional spaces; therefore the learning process of such VQ task may take hours or days to converge, especially if the number of clusters is large too.

In this paper, we study how, and up to which point, the Kohonen Self-Organizing Maps (SOMs) may be used to accelerate the VQ learning process. It is well known that SOMs accelerate vector quantization, compared to other, more traditional VQ algorithms. However, to our knowledge, it exists no quantitative experimental study on this topic. To fill this gap is the main purpose of this paper. After reminding formal definitions of the conventional VQ and SOM algorithms, and in particular of their distortion measures (section 2), we study the rate of convergence of conventional VQ and SOM algorithms, both on artificial databases where the exact solution of VQ can be computed, and on real ones (sections 3 and 5). Section 4 uses a hybrid method with a fixed-neighborhood SOM as initialization to a conventional VQ algorithm, in order to benefit both from the accelerated convergence of SOM and from the lower distortion error after convergence of conventional VQ; this hybrid algorithm is similar to the conventional SOM with decreasing neighborhood, and is aimed to better quantify the gain in convergence speed.

Vector quantization

Principle and distortion

Vector quantization consists in replacing a continuous distribution, in some cases known only through a finite number of samples, by a finite number of quantizers; the number of quantizers must be (much) smaller than the number of known samples. Each quantizer defines a cluster in the space; the principle of vector quantization is to project all samples in a cluster on the corresponding quantizer.

Most of the methods used to perform VQ necessitate setting a priori the number of clusters or quantizers. The choice of this number results from a trade-off between the precision (distortion) of the quantization and the necessity of an efficient description of the resulting clusters (quantity of information kept after quantization).

Once the number of quantizers is predefined, a good criterion of the quality of the clustering is the distortion, which measures the deviation between the data and the corresponding quantizers. Let us recall the main definitions and introduce our notations.

We consider a continuous data space Ω, of dimension d, endowed by a continuous probability density function (pdf) f(x), where the cumulated density (or repartition function) is F(x)=P(X<x) (where P is the probability law, and where the inequality is verified in each dimension).

A vector quantization Φ is an application from the continuous space Ω to a finite subset F (the codebook) formed by n code-vectors or centroids or quantizers q 1 , q 2 ,…,q n of Ω. The positions of the code-vectors are supposed to be computed as a result of a quantization algorithm or learning algorithm.

The aim of a vector quantization (VQ) is to compress the information by replacing all elements x of a cluster C i (subset of Ω) by a unique quantizer (or code-vector, or centroid) q i . For a given number n of code-vectors, vector quantization tries to minimize the loss of information or distortion, measured by the mean quadratic error

∑ ∫ = - = ξ = Φ ξ n i C i n i dx ) x ( f q x ) q , , q , q , f ( ) , f ( 1 2 2 1 Κ . ( 1
)
If a N-sample x 1 , x 2 ,…,x N is available (randomly chosen according to f(x)), this distortion is estimated by the intra-class sum of squares

∑ ∑ = ∈ - = ξ = Φ ξ n i C x i j n i j q x N ) q , , q , q , f ( ) , f ( ˆ1 2 2 1 1 Κ . (2)
All classical VQ algorithms (LBG, k-means,…) minimize this distortion function by choosing appropriate centroid locations. See for example Anderberg [START_REF] Anderberg | Cluster Analysis for Applications[END_REF] or Bishop et al. [START_REF] Bishop | Discrete Multivariate Analysis: Theory and Practice[END_REF] for proofs. There is no unique minimum of the distortion function, and the result strongly depends on the initialization.

Simple Competitive Learning and batch VQ algorithms

There exist many algorithms that deal with the VQ problem. Most of them are very slow in terms of convergence speed. The most popular one is the so-called Simple Competitive Learning algorithm (SCL) that can be defined as follows (see for example [START_REF] Hertz | Introduction to the Theory of Neural Computation[END_REF]):

let Ω be the data space (with dimension d), endowed with a density probability function f(x). The data are randomly drawn according to the density f(x) and are denoted by x 1 , x 2 ,…,x N . The number of desired clusters is a priori fixed to be n. The quantizers q 1 , q 2 , …, q n are randomly initialized. At each step t,

• a data x t+1 is randomly drawn according to the density f(x) ;

• the winning quantizer q win(t) is determined by minimizing the classical Euclidean norm || x t+1q win(t) || = min i || x t+1 -q i || ;

• the quantizer q win(t) is updated by q win(t+1) = q win(t) + ε(t) (x t+1q win(t) ).

where ε(t) is an adaptation parameter which satisfies the classical Robbins-Monro conditions (Σ ε(t) = ∞ and Σ ε 2 (t) < ∞).

The SCL algorithm is in fact the stochastic or on-line version of the Forgy algorithm (also called moving centers algorithm, LBG or Lloyd's algorithm, see for example [4], [8], [9]). In that version of the algorithm, the quantizers are randomly initialized. At each step t, the clusters C 1 , C 2 , …, C n are determined by putting in class C j the data which are closer to q j than to any other quantizer q i . Then the mean values in each cluster are simultaneously computed and taken as new quantizers. The process is then repeated. The Forgy algorithm works off-line as a batch algorithm. It will be referred to as BVQ (for Batch VQ) in the following. It also exits an intermediate version of the algorithm, frequently named the K-means method [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF]. In that case, at each step, only one data is randomly chosen, and only the winning quantizer is updated as the mean value of its class.

It can be proven ([1], [START_REF] Bishop | Discrete Multivariate Analysis: Theory and Practice[END_REF]) and it is well-known that BVQ minimizes the distortion (1) or, more exactly, the estimated one [START_REF] Bishop | Discrete Multivariate Analysis: Theory and Practice[END_REF]. Note that the stochastic SCL algorithm also minimizes this distortion, but only in mean value: at each step, there is a positive probability to increase the distortion, as for any stochastic algorithm.

Let us denote by q 1 *, q 2 *, …, q n * one set of quantizers which (locally) minimizes the distortion. This minimum needs not to be unique and generally depends on the initial values. At a local minimum of the distortion, each q i * is the center of gravity of its class C i , with respect to the density f, and the quantizers are fixed points of the BVQ algorithm. In an exact form,

∫ ∫ = i i C C i dx ) x ( f dx ) x ( f x * q , (3)
estimated by

∑ ∑ ∈ ∈ = i j i j C x C x j i x * q ˆ1 . ( 4
)
Note that the equations are implicit ones, since the C i are defined according to the positions of the q* i .

For example, in the one-dimensional case, the classes

C i (1 ≤ i ≤ n) are intervals defined by C i = [a i , b i ],
with a i = ½ (q* i-1 + q* i ) and b i = ½ (q* i+1 + q* i ), for 1 < i <n, and

a 1 = inf (Ω), b n = sup (Ω).
The main goal of this paper (after a preliminary work presented to ESANN'99, see [START_REF] De Bodt | Using the Kohonen Algorithm for Quick Initialization of Simple Competitive Learning Algorithm[END_REF]) is to evaluate the speed of convergence of VQ algorithms. In situations where the solution is unique and where it is possible to compute the exact values q* i , the performances will be evaluated through the rate at which the values q i converge to q* i (see section 3).

Self-Organizing Maps (SOM)

Let us consider now the SOM algorithm (as defined by T. Kohonen in [START_REF] Kohonen | Self-organizing maps[END_REF]). As we will show it, it can be seen as an extension of the Simple Competitive Learning Algorithm (SCL) in its classical stochastic form, and of the Forgy Algorithm (BVQ) in its batch form. We will consider here the SOM algorithm with a fixed number of neighbors (although the number of neighbors uses to decrease with time in practical implementations).

For a given neighborhood structure, where V(i) denotes the neighborhood of unit i, the SOM algorithm is defined as follows. The quantizers q 1 , q 2 , …, q n are randomly initialized. At each step t,

• a data x t+1 is randomly drawn according to the density f(x) ;

• the winning quantizer q win(t) is determined by minimizing the classical Euclidean norm || x t+1q win(t) || = min i || x t+1 -q i || ;

• the quantizer q win(t) and its neighbors q k (t) for k in V(j) are updated by

q win,k(t+1) = q win,k(t) + ε(t) (x t+1 -q win,k(t)
). where ε(t) is an adaptation parameter which satisfies the classical Robbins-Monro conditions (Σ

ε(t) = ∞ and Σ ε 2 (t) < ∞).
We see that the SCL algorithm is a particular case of the SOM algorithm, when the neighborhood is reduced to zero. Sometimes SCL is called 0-neighbor Kohonen algorithm. There also exists a batch SOM algorithm, similar to the Forgy (BVQ) algorithm. The only difference is that at each step, for a given set of classes C 1 , C 2 , …, C n , the quantizer q j is set to the mean value of the union of the class C j and its neighbors. See [START_REF] Kohonen | Comparison of SOM Point Densities Based on Different Criteria[END_REF] for example.

It appears clearly that the SOM algorithm is different from the SCL algorithm only because a neighborhood structure is defined between the n quantizers. The neighborhood structure of the SOM algorithm is most frequently used for visualization and data interpretation properties. We however only consider it here as a way to accelerate the convergence of the SOM algorithm. In other words, we are only interested here in the VQ property of self-organizing maps, and not in its topological properties.

In the one-dimensional case, and for a one-dimensional structure of neighborhood, if the neighborhood V(i) contains i-1, i, i+1 (two-neighbor case), the limit points q* i of the batch SOM algorithm verify equation (3) or equation [START_REF] Forgy | Cluster analysis of multivariate data: efficiency versus interpretability of classifications[END_REF], where C i is replaced by

C i 2 = C i-1 ∪C i ∪C i+1 = [a i , b i ],
with a i = ½ (q* i-2 + q* i-1 ) and b i = ½ (q* i+1 + q* i+2 ), for 2 < i <n-1, and a 1 = a 2 = inf (Ω), and b n-1 = b n = sup (Ω).

Here again, the batch SOM algorithm is nothing else than the iterative computation of the solutions of equations ( 3) or (4), when C i is replaced by C i 2 .

The batch SOM algorithm and the classical stochastic SOM algorithm do not decrease anymore the distortion (1), but a generalized distortion that is the distortion extended to the neighbor classes (as long as the number of neighbors ν remains fixed) [START_REF] Ritter | On the Stationary State of Kohonen's Self-Organizing Sensory Mapping[END_REF]. This generalized distortion is given by

∑ ∫ = ν ∈ - = ξ = Φ ξ n i C i n v ) i ( V k k dx ) x ( f q x ) q , , q , q , f ( ) , f ( 1 2 2 1 Υ Κ , (5)
where V(i) is the set of indexes in the neighborhood of i, including i. This generalized distortion function can also be estimated through a finite set of samples x 1 , x 2 ,…,x N , similarly to (2).

Experimental results: convergence to the exact solution of the VQ problem

The SOM algorithm is not equivalent to the SCL algorithm: it is deemed to minimize the generalized distortion (5), and not the classical distortion of VQ problems [START_REF] Anderberg | Cluster Analysis for Applications[END_REF]. Despite this fact, we will show that the SOM algorithms perform better than the classical SCL algorithm, i.e. converge faster towards the solution of (1), at least during the first iterations of the algorithm.

We study this phenomenon from two points of view. First, in some cases where it is possible to exactly compute the solutions of equations ( 3) or (4), we evaluate the error between the current values and the optimal values as a function of the number of iterations, for both the SCL and SOM algorithms. This is the topic of this section.

Secondly, for more realistic data, we compare the decreasing slope of the true distortion (1) as a function of the number of iterations, also for both algorithms. This is the topic of section 5

In some one-dimensional cases (d = 1), if the set Ω is a real interval, and if the density f is known and wellbehaved, it is possible to directly compute the solutions q* i , starting from a given set of increasing initial values, by an iterative numerical procedure.

If the initial values are ordered, the current values q 1 , q 2 , …, q n will remain ordered at each iteration of the SCL algorithm. As mentioned in the previous section, the classes C i (1 ≤ i ≤ n) are therefore intervals defined by

C i = [a i , b i ],
with a i = ½ (q i-1 + q i ) and b i = ½ (q i+1 + q i ), for 1 < i <n, and a 1 = inf (Ω), b n = sup (Ω). This constitutes the first set of equations (C i as a function of q i ) used in this iterative procedure.

Equations (3) or (4) have no explicit solutions. However, it is possible to compute analytically the solutions q i , as a function of the limits a i and b i of the intervals C i , for some "easy" densities f(x). This will constitute the second set of equations (equivalent to (3)) used in the iterative procedure. Table 1 presents these recurrent explicit equations for the densities f(x) = 2x, 3x 2 , e -x .

Density f Distribution Function Using alternatively the two sets of equations leads to the convergence to the optimal values q* i of the quantizers. This iterative procedure is similar to the BVQ algorithm. Formulas in Table 1 are the analytical solutions of equation ( 3), while the BVQ algorithm usually involves in practical experiments the use of approximation (4).

a 0 b 0 q i 2x on [0,1] x 2 0 1 2 2 3 3 3 2 i i i i i a b a b q - - = 3x 2 on [0,1] x 3 0 1 3 3 4 4 4 3 i i i i i a b a b q - - = e -x on [0,+∞] 1-e -x 0 +∞ i i i i i i
Knowing the optimal values of the quantizers, it is possible to study the speed of convergence of any Vector Quantization algorithm. In that purpose, we will study the Euclidean distance between the current values of (q i (t)) resulting from some VQ algorithm and the solutions (q* i ), as a function of the numbers of iterations. We define the mean quadratic error

D 2 (t) = D 2 (q(t),q*) = (1/n) Σ 1≤ i ≤ n (q i (t) -q* i )² (6)
which will be the error measure of the Vector Quantization algorithm into consideration.

In practical situations, one can observe that the error measure D 2 (t) decreases to 0 very slowly when using the SCL algorithm. Note that in all simulations with several algorithms we carefully started from the same initial increasing values q i (0), including for the exact computation of the (q* i ), in order to avoid any effect due to the initial conditions.

In Figures 1, 2 and 3, we represent the variations of the error measure D 2 (t), for the SCL algorithm and for the SOM with 2, 4 and 8 neighbors; 100 quantizers were used for all simulations on artificial data. Figures 1, 2 and 3 respectively concern densities f(x) = 2x, 3x 2 and e -x . We can see for example that the SOM with neighbors decreases to the optimal values (q* i ) much faster than the SCL algorithm, even if it finally converges to its own optimal points. These optimal points minimize the generalized distortion extended to the neighbors (5), and are different from the (q* i ) (see section 2). We also measured the evolution of D 2 (t) as a function of the number of iterations, for a Gaussian density N(0,1).

In this case, the exact values q * i have been obtained through equation ( 4), by using very large samples to compute at each step the C i . Figure 4 shows this evolution of D 2 (t) as a function of the number of iterations, respectively for the SCL algorithm and for the SOM with 2, 4, 8 and 16 neighbors. This section has shown that the use of SOM can greatly increase the speed of convergence towards the exact solutions of the VQ problem. Nevertheless, it must not be forgotten that the SOM algorithm will not finally converge to these solutions but rather to a minimum of (5). In the next section, we therefore use a mixed algorithm, beginning by some iterations of the SOM algorithm and ending with a classical SCL procedure, in order to benefit both from the accelerated convergence and from the convergence towards optimal states.

Hybrid algorithm SOM/SCL

Based on the results of the previous section, we propose to use a hybrid VQ algorithm (denoted by KSCL for Kohonen SCL), which consists of an initial phase (a SOM algorithm with ν neighbors), followed by the classical SCL. We compare the value of the error D 2 (t) after the same number of iterations for KSCL and SCL. Note that KSCL is not far from the classical SOM algorithm with decreasing neighborhood; the reason to use KSCL is not to suggest another algorithm, but to better quantify the gain in convergence speed compared to classical SCL, regardless of the specific decreasing function used in a classical SOM.

For example, let us fix a total number of iterations T, the initial ordered points q 1 (0), q 2 (0), …, q n (0), a constant ε and several probability functions :

f(x) = 2x on [0,1], f(x) = 3x 2 on [0,1], f(x) = e -
x on [0, +∞ [, and the standard Gaussian N(0,1). Let us also consider the 2-neighbors SOM algorithm (ν = 2).

In Figures 5,6, 7 and 8, we represent the evolution of the error measure for different KSCL algorithms and for the four probability densities that we took as examples. We consider four KSCL variants where the 2-neighbors SOM algorithm is used respectively during 0%, 30%, 60%, 90% of the total number of iterations T. We can observe in all simulations that the 2-neighbors algorithm greatly accelerates the decrease of the error measure in the beginning of the curves. In all cases, using too early the SCL algorithm slows down the decrease.

Moreover, the performances remain better than those of the SCL algorithm, whatever is the choice of the KSCL variant. Nevertheless, it is also clear that determining the optimal iteration for substituting SOM by SCL strongly depends on the probability density. An optimal choice of this parameter would thus require extensive simulations, which is not the goal searched for here.

We may conclude this section by claiming that, in any case, the SOM algorithm with a fixed number of neighbors can work as an efficient initialization of the SCL algorithm to accelerate the convergence and improve the performances. We verified this statement on many other probability densities, real data and for several values of the number of neighbors' ν.

In practical implementations of the SOM algorithm, the number of neighbors is made decreasing. Our observations confirm that this widely used strategy is very efficient to improve the decreasing of the error measure (6).

It would be interesting to consider this so-called error measure in multidimensional settings. Nevertheless, this concept is not well suited to dimensions greater than 1. The lack of ordering concept in dimension greater than 1 does not facilitate the problem and the correspondence between the current quantizers at a given iteration and their optimal values looses its clear meaning.

We will thus replace the concept of error measure by the concept of distortion as defined by ( 1) or (2). In the next section, we study how the distortion is decreasing along the quantization process, in both cases (SCL without neighbors, or SOM with a decreasing number of neighbors).

Experimental results: comparative evolution of the distortion

In this part, we study the vector quantization speed of SCL and of SOM, by computing the distortion defined in equation ( 2) in the case of real data.

As we mentioned previously, the SCL is supposed to minimize this distortion, while the SOM (with fixed or decreasing number of neighbors) is not. However we can observe that in any case, the SOM algorithm accelerates the decrease of the distortion, at least during a large part of the simulation.

We represent the distortion as a function of the number of iterations, for 5 different quantization algorithms: SCL and 4 variants of the SOM algorithm which differ by the number of neighbors. For a two-dimensional neighborhood structure, we consider successively 3 SOM algorithms with a fixed number of neighbors (SOM5, SOM9 and SOM25, the suffix being the number of neighbors) and then the classical SOM algorithm with a decreasing number of neighbors (from 25 to 1, the last part of the SOM iterations being equivalent to SCL).

We illustrate these simulations on two data sets. The table SAVING contains 5 variables measuring economic ratios for 42 countries between 1960 and 1970; the table TOP500 contains 6 variables relative to 77 American companies in 1986. 1Figure 9 represents the distortion for the SAVING data, with 25 quantizers, square grid (5 by 5) for the SOM algorithms and 500 steps of iterations. Figure 10 represents the distortion for the data SAVING, with 100 quantizers, square grid (10 by 10) for SOM algorithms and 1000 iterations. Figure 11 represents the distortion for the data TOP500, with 100 quantizers, square grid (10 by 10) for SOM algorithms and 1000 iterations. In each simulation, we can see that the SOM algorithm performs as the best quantizer (it leads to a lower minimum of the distortion function). The SCL algorithm is very slow, and the SOM with non-decreasing number of neighbors is powerful at the beginning of the iterations, but allows at some iteration the distortion to increase (or would allow, after a larger number of iterations). In fact the classical SOM algorithm ending with no neighbor appears to be an excellent VQ algorithm. When the quality of the result is the ultimate goal (regardless of the computation time), one can use the SCL algorithm after the SOM one to refine the solution (or in other words, one can increase the number of iterations performed without neighbors in the SOM procedure). Indeed even if a classical SOM usually ends without neighbor, the number of iterations performed without neighbors could be not sufficient to reach an optimal solution. Performing a sufficient number of iterations both with and (then) without neighbors is thus important for the quality of the solution.

Conclusion

The experiments illustrated in this paper, as well as many other ones performed on other data sets, indicate that the quality of the SOM algorithm resides not only in its topology preservation property, but also its vector quantization one. The SOM algorithm may be recommended compared to other VQ techniques like SCL, in order to reach a better minimum of the distortion error with a fixed number of iterations, or to reach faster a similar value of the distortion.

The better convergence properties cannot be proven theoretically. Nevertheless, we can make the analogy with simulated annealing techniques: the use of a neighborhood in the SOM algorithm introduces apparent disorder, making it possible to escape from a local minimum of the objective function and to increase the slope of convergence. Ending the VQ procedure with the SCL algorithm may be compared to ending a simulated annealing technique with a "temperature" parameter equal to zero.

The fact that the SOM method converges faster than classical VQ algorithms is well known; however, few studies bring systematic evidences supporting this claim. This paper shows quantitatively, through extended simulations, that the speed of convergence is increased by the use of the neighborhood property in the SOM. On artificial databases, where the optimal VQ solution may be computed analytically, we measured the VQ speed by the difference between the values of the quantizers during the convergence and their optimal stable solutions. On real databases where the optimal VQ solutions cannot be computed, we measured the VQ speed by the decrease of the distortion error with respect to the number of iterations. In both cases, we showed quantitatively how the use of the SOM algorithm, with fixed or decreasing number of neighbors, accelerates the convergence of the VQ learning. This produces clear indications that SOMs should be preferred to conventional VQ methods, even if topology preservation is not looked for.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Evolution of D 2 (t) as a function of the number of iterations, for the density f(x) = 2x.

Figure 4 :

 4 Figure4: Evolution of D 2 (t) as a function of the number of iterations, for the standard Gaussian N(0,1) density One could argue that the comparisons are made on different algorithms, where the processing time per iteration is different. The comparison in terms of number of iterations is thus not fair if the total processing time is searched for. Nevertheless, as an example, the difference of the processing time in our simulations of a single iteration when using the 2-neighbors SOM algorithm instead of the SCL algorithm is significantly less than 1%. The differences shown in Figures 1 to 4 thus remain striking.

Figure 7 :Figure 8 :

 78 Figure 7: Evolution of D 2 (t) as a function of the number of iterations, for 4 variants of the KSCL algorithm, on the e -x density.

Figure 9 :

 9 Figure 9: Evolution of the distortion as a function of the number of iterations, on the Saving data set; 25 quantizers are used (see text for details on the algorithms). The distortion is illustrated after each iteration.

Figure 10 :

 10 Figure 10: Evolution of the distortion as a function of the number of iterations, on the Saving data set; 100 quantizers are used (see text for details on the algorithms).

Figure 11 :

 11 Figure 11: Evolution of the distortion as a function of the number of iterations, on the Top500 data set; 100 quantizers are used (see text for details on the algorithms).

Table 1 :

 1 Exact computation of the quantizers as a function of the limits of the clusters, for some "easy" examples of densities.

  Evolution of D 2 (t) as a function of the number of iterations, for 4 variants of the KSCL algorithm, on the 2x density.

	0.007 0.02
	0.018
	0.006
	0.016
	0.005 0.014
	0.004 0.012
	0.01
	0.003 0.008
	0.002 0.006
	0.004
	0.001
	0.002
	0 0
	0 0	20000 20000	40000 40000	60000 60000	80000 80000	100000 100000
		KSCL 00 % KSCL 00 %	KSCL 30 % KSCL 30 %	KSCL 60 % KSCL 60 %	KSCL 90 % KSCL 90 %
	0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 Figure 5: 0 0.02
	0	20000	40000	60000	80000	100000
		KSCL 00 %	KSCL 30 %	KSCL 60 %	KSCL 90 %
	Figure 6: Evolution of D 2 (t) as a function of the number of iterations, for 4 variants of the KSCL algorithm, on
			the 3x 2 density.

The data are available on the http://www.dice.ucl.ac.be/neural-nets/Research/Projects/projects.htm WEB page.

Acknowledgements

Michel Verleysen is a Senior Research Fellow of the Belgian National Fund for Scientific Research (FNRS). The authors thank the reviewers for their comments and suggestions.


S Balac

H Benoit-Catton

T Lamotte

C Odet

Analytic Solution to Boundary Integral Computation of Susceptibility Induced Magnetic Field Inhomogeneities

Keywords: Boundary integrals, Susceptibility artifacts, Magnetostatics

A method to compute the magnetic field induced by susceptibility inhomogeneities in magnetic resonance imaging is presented. It is based on a boundary integral representation formulae. The integral is set over the surfaces between media of different magnetic susceptibilities. The computational procedure consists of approximating these surfaces with triangular mesh elements and using analytical expressions to compute the integral over each triangle. The proposed method supplies high accuracy and is easily paralleled. A detailed analysis for the convergence rate of the method is performed. Numerical results obtained for several samples, including a human head, are presented.

INTRODUCTION

When a sample is magnetized by a static magnetic field, as during a magnetic resonance imaging (MRI) experiment, the magnetic field around the sample is distorted. The magnetic susceptibility distribution in the sample determines the local perturbation of the magnetic field. In MRI, differences in susceptibility between a metallic implant and the surrounding anatomic tissues, or between air and tissues, are one of the greatest sources of geometric distortions in the images, see [START_REF] Ludeke | Susceptibility artifacts in NMR imaging[END_REF][START_REF] Bakker | Analysis of machine-dependent and objectinduced geometric distorsion in 2DFT MR imaging[END_REF] for instance. These distortions are harmful for applications that require accurate quantitative measurements from the images, see [START_REF] Chang | A technique for accurate magnetic resonance imaging in the presence of field inhomogeneities[END_REF][START_REF] Weis | Magnetic resonance spectroscopic imaging for visualization and correction of distorsions in MRI: High precision applications in neurosurgery[END_REF].

Since the knowledge of the static magnetic field perturbation is necessary to rectify distorted MR images, many studies deal with its quantification, see [START_REF] Bui | A quantitative study of the pixel-shifting, blurring and nonlinear distortions in MRI images caused by the presence of metal implants[END_REF][START_REF] Sumanaweera | Quantifying MRI geometric distortion in tissue[END_REF], and its computation. Analytic methods have been proposed to study simple geometries objects such as spheres or cylinders, see [START_REF] Ludeke | Susceptibility artifacts in NMR imaging[END_REF]. For more complex geometries, a precise computation of the magnetic field perturbation involves partial-differential equations derived from Maxwell's equations and requires the use of numerical methods. Among the classical numerical methods are the finite-difference method (FDM), the finite-element method (FEM), and the boundary elements method (BEM) (we refer to [START_REF] Sadiku | Numerical Techniques in Electromagnetism[END_REF] for a comprehensive presentation of these numerical methods in electromagnetism). These methods have been widely used to compute the magnetic field perturbation generated by susceptibility inhomogeneities during an MRI experiment. For instance, Bhagwadien et al. [START_REF] Bhagwandien | Numerical analysis of the magnetic field for arbitrary magnetic distribution in 3D[END_REF] developed an iterative finite-difference method for computing the static field for arbitrary magnetic susceptibility distributions. Li et al. [START_REF] Li | A computer simulation of the static magnetic field distribution in the human head[END_REF] used commercial finite-element software to compute the three-dimensional magnetic field in a human head. Based on the previous finite-element approach, Truong et al. [START_REF] Truong | Three dimensional numerical simulations of susceptibility induced magnetic field inhomogeneities in the human head[END_REF] proposed a multiresolution computation of the magnetic field and its application to several anatomical models. Jenkinson et al. [START_REF] Jenkinson | Perturbation calculation of B0 field for non conducting material[END_REF] proposed an approach based on a perturbation method which is adapted to a nonconductive object represented by rectangular voxels restricted to air and singular tissue. The boundary element method has been used by Hwang et al. [START_REF] Hwang | The calculation of the susceptibility induced magnetic field from 3D NMR images with applications to trabecular bone[END_REF] to simulate the magnetic field perturbation on trabecular bone MR images. The BEM is used as well by De Munck et al. [START_REF] De Munck | The computation of MR image distortions caused by tissue susceptibility using the boundary element method[END_REF] to compute distortions on MRI images from a human head. Yoder et al. [START_REF] Yoder | MRI simulator with static field inhomogeneity[END_REF] simplified the BEM approach in the case of homogeneous voxels perpendicular to the static field and integrated the computed field perturbation in a MRI simulator. The BEM relies on a computation over the interfaces between domains of different susceptibilities. The magnetic field perturbation is obtained through the resolution of a linear system whose coefficients are determined by evaluating surface integrals. In the quoted works, the surface integrals were approximated considering only small susceptibility values corresponding to anatomical tissues values.

In [START_REF] Balac | Induced magnetic field computations using a boundary integral formulation[END_REF], Balac and Caloz present a method to compute the magnetic field induced by metallic bodies with high constant magnetic susceptibility (such as medical implants) based on an integral representation formula over the boundary of the body. The method uses analytic expressions to compute the integral over the flat panels of the boundary and a piecewise quadratic interpolation of the surface for the curved panels. It supplies both high accuracy and low computation time. In this paper, we generalize the method presented in [START_REF] Balac | Induced magnetic field computations using a boundary integral formulation[END_REF] to deal with bodies composed of several domains of arbitrary magnetic susceptibilities. The magnetic field is computed through an explicit integral representation formula over the boundary interfaces between the interior domains. The advantages of our method compared to the FEM or to the FDM are numerous. First, the computation depends only on the interface boundaries between media of different susceptibilities so that the computational problem is reduced from a 3D problem to a 2D problem. Furthermore, in the FDM/FEM the exterior domain must be truncated and an approximation of the behavior of the field at infinity must be introduced on an artificial boundary. In our approach, as well as with the BEM, the behavior at infinity is always exactly satisfied. Finally, the FEM as the FDM or the BEM lead to large linear systems to be solved, whereas in the proposed method the solution is obtained pointwise by evaluating a surface integral. Our computational method assumes the interface boundaries between media of different susceptibilities are meshed in a collection of flat triangles. The advantage is that the surface integral can be evaluated exactly with no limitation on the susceptibility values or on the static field direction. The only approximation in the method stands in the way the true interfaces are approached by the union of flats triangles and the approximation error is controlled by the size of the triangles. Thus, the computation of susceptibility induced magnetic field inhomogeneities is done by a unique and accurate method with a low computation time for both soft tissues and metallic objects. Moreover, our method can be parallelized in a very simple and efficient way, the computation time being divided by a factor equal to the number of computation nodes.

Our paper is organized as follows. In Section 2, we formulate the magnetostatic problem and give the integral formula used to compute the magnetic field. Section 3 is devoted to the presentation of the computational method and to a detailed analysis of its accuracy. In Section 4, we present results obtained with our computational method. The validation is done on a diamagnetic and a paramagnetic ball. The convergence rate and the computation time are discussed. Finally, results obtained on a hollow ball and on a human head are presented.

FORMULATION OF THE PROBLEM

The Magnetostatic Problem

We denote by Ω the domain filled by the whole sample. Ω is an open bounded domain in R 3 , whose boundary is denoted by Σ. We assume that Ω is the union of N domains Ω i with distinct magnetic properties. For convenience, we denote by Ω N +1 the exterior domain. We denote by Σ ij the common boundary to Ω i and Ω j , 1 ≤ i < j ≤ N + 1 (possibly we have Σ ij = ∅) and set

I = ? (i, j) ∈ N 2 | 1 ≤ i < j ≤ N + 1 and Σ ij ? = ∅ ? .
We denote by ∂Ω k the whole boundary of Ω k , k ∈ {1, . . . , N }.

The sample is assumed to be embedded in the magnetic field B 0 of the MRI device. We consider B 0 constant in intensity and direction. In this modeling, we neglect the effects of the RF field B 1 of the MRI device on the static field inhomogeneities.

Let B be the total magnetic flux density and H be the total magnetic field intensity. The basic equations of magnetostatics are, see [START_REF] Jackson | Classical Electrodynamics, Third Edition[END_REF],

div B = 0, in R 3 , rot H = 0, in R 3 . (1)
At a distance far from Ω the magnetic field tends to become homogeneous; we have lim

P →∞ B(P ) = B 0 . (2)
In the exterior region Ω N +1 , H and B are related through the relation

H = 1 µ 0 B, (3)
where µ 0 is the magnetic permeability of vacuum, while in Ω they are connected to the magnetization M by the relation

H = 1 µ 0 B -M. (4)
For the magnetic flux induction B ? = B -B 0 the problem reads: find

B ? in the space L 2 (R 3 ) 3 such that div B ? = µ 0 rot M, in R 3 , rot B ? = 0, in each Ω i , i ∈ {1, . . . , N + 1}, [B ? ∧ n] = µ 0 [M ∧ n] ,
at the interfaces Σ ij , (i, j) ∈ I,

where n is the unit normal to Σ ij from Ω i to Ω j and [v] denotes the jump of v across the boundary Σ ij : [v] = v| Ωi -v| Ωj . To deal with a well-posed problem, we still need to add some relation between the magnetization M and the field B ? . The domains Ω i are assumed to be isotropic linear paramagnetic or diamagnetic materials. In each domain Ω i , the magnetization M is then related to the applied field B 0 through the linear relation

M = χ (i) m µ 0 B 0 , (6)
where

χ (i)
m is the magnetic susceptibility in Ω i . The first consequence of this assumption is that on each Ω i we have rot M = 0. The second one is that the interface conditions in (5) take the form

[B ? ∧ n] = j Σij , (7)
where

j Σij = [χ m ] Σij (B 0 ∧ n) = ? χ (i) m -χ (j) m ? (B 0 ∧ n). (8)
We can check that div j Σij = 0, so that j Σij can be interpreted as a fictitious surface current density over the surface Σ ij . Finally, we look for B ? such that div

B ? = 0, in R 3 , rot B ? = 0, in each Ω i , i ∈ {1, . . . , N + 1}, [B ? ∧ n] = j Σij , at the interfaces Σ ij , (i, j) ∈ I. (9)
2.2. The Integral Formula for the Magnetic Field

Our method originates from the Biot and Savart law that states the solution to problem (9) can be expressed as, see [START_REF] Jackson | Classical Electrodynamics, Third Edition[END_REF],

B ? (P ) = 1 4π ? (i,j)∈I ? ? Σij j Σij (Q) ∧ r(P, Q) r 3 (P, Q) dσ Q , (10)
where r(P, Q) = PQ. This relation is valid in any point P that does not belong to any of the interfaces Σ ij . However, it is not convenient for numerical purposes in its present form and we use the following vector identity to transform relation [START_REF] Truong | Three dimensional numerical simulations of susceptibility induced magnetic field inhomogeneities in the human head[END_REF] in a more suitable form: if u, v, w denote three vectors in R 3 , we have

(v ∧ w) ∧ u = (u • v) w -(u • w) v.
It follows that

(B 0 ∧ n) ∧ r(P, Q) r 3 (P, Q) = ? r(P, Q) r 3 (P, Q) • B 0 ? n - ? r(P, Q) r 3 (P, Q) • n ? B 0
and the solution to problem [START_REF] Li | A computer simulation of the static magnetic field distribution in the human head[END_REF] given by ( 10) is

B ? (P ) = 1 4π ? (i,j)∈I [χ m ] Σij ? - ? ? Σij ? r(P, Q) r 3 (P, Q) •n ? B 0 dσ Q + ? ? Σij ? r(P, Q) r 3 (P, Q) •B 0 ? n dσ Q ? . (11)
As the source magnetic field B 0 is constant, the first integral term in [START_REF] Jenkinson | Perturbation calculation of B0 field for non conducting material[END_REF] reads

- 1 4π
? (i,j)∈I [χ m ] Σij ? ? Σij ? r(P, Q) r 3 (P, Q) • n ? B 0 dσ Q = - 1 4π   ? (i,j)∈I [χ m ] Σij ? ? Σij r(P, Q) r 3 (P, Q) • n dσ Q   B 0 = ? χ m (P )B 0 , if P ∈ Ω, 0, if P / ∈ Ω. (12)
Indeed, if S is a close surface then the integral

A(P, S) = 1 4π ? ? S r(P, Q) r 3 (P, Q) • n dσ Q
is the so-called "solid angle" whose value is

A(P, S) = ? -1, if P ∈ V, 0, if P / ∈ V,
where V denotes the domain enclosed in S.

Finally, the Biot and Savart integral [START_REF] Truong | Three dimensional numerical simulations of susceptibility induced magnetic field inhomogeneities in the human head[END_REF] had been transformed into the following expressions. If P belongs to Ω but at none of the interfaces Σ ij , then

B ? (P ) = χ m (P ) B 0 + 1 4π ? (i,j)∈I [χ m ] Σij ? ? Σij ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q = χ m (P ) B 0 + 1 4π N ? k=1 χ (k) m ? ? ∂Ω k ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q . ( 13
)
If P belongs to one of the interface Σ ij , then B(P ) is undefined. The reason being that the magnetic induction is not continuous across boundaries of medium of different susceptibilities, see equations [START_REF] Bui | A quantitative study of the pixel-shifting, blurring and nonlinear distortions in MRI images caused by the presence of metal implants[END_REF]. If P belongs to the exterior domain Ω N +1 , then

B ? (P ) = 1 4π ? (i,j)∈I [χ m ] Σij ? ? Σij ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q = 1 4π N ? k=1 χ (k) m ? ? ∂Ω k ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q . (14)
Thus, in both cases, the computation of the magnetic field perturbation B ? turned to the calculation of the following surface integrals over the boundaries ∂Ω k , k = 1, . . . , N ,

I k = ? ? ∂Ω k ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q . ( 15
)
3. THE COMPUTATIONAL METHOD

Overview of the Method

We assume that each boundary ∂Ω k , k = 1, . . . , N , has a parameterization F k : D -→ ∂Ω k , where D is a polygonal domain in the plane which is two times continuously differentiable (or that it can be decomposed in a union of smooth closed surfaces having this property). To the polygonal domain D, we associate a triangulation Th , Th = { Ki , 1 ≤ i ≤ N } where Ki are triangles in the plane. The mesh size h of the triangulation is defined by

h = max 1≤i≤N diam ? Ki ? , where
v i = F (v i ), i = 1, 2, 3
). Thus, we can write I k as As the elements K in T h are curved triangles in the space, it is not in general possible to proceed to an exact calculation of I k . We then compute an approximation of I k by using a piecewise linear interpolation of the surface ∂Ω k . The approximated surface ? ∂Ω k will be composed of flat triangles K in the space that are obtained as the linear interpolants of the curved triangular elements K in T h . We denote by Th the triangulation of ?

I k (P ) = ? K∈T h ? ? K ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q . ( 16
)
∂Ω k obtained in this way. We have

∂Ω k = ? K∈T h K ≈ ? K∈ Th K = ? ∂Ω k . ( 17
)
If K ∈ T h is the image of the triangle K ∈ Th through the parameterization F , the corresponding triangle K ∈ Th is the image of K through the linear interpolant FK of F over K. It follows that the integral I k over the boundary ∂Ω k can be approximated by the following integral over the approximate boundary ? ∂Ω k :

Ĩk (P ) = ? ? ? ∂Ω k ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q = ? K∈ Th ? ? K ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q . ( 18
)
Since the triangles K ∈ Th are flat, the normal vector n is constant on each triangle K. Moreover, as the magnetic field B 0 is constant in strength and direction, Ĩk may be written as

Ĩk (P ) = ? K∈ Th ? B 0 • ?? ? K r(P, Q) r 3 (P, Q) dσ Q ? n ? . (19)
Thus, the problem now reduces to the computation of the following integral over each triangle K ∈ Th :

JK (P ) = ? ? K r(P, Q) r 3 (P, Q) dσ Q . ( 20
)
There are several ways to compute analytically this integral, see [START_REF] Balac | Magnetic susceptibility artifacts in magnetic resonance imaging, Calculation of the magnetic field disturbances[END_REF][START_REF] Okon | The potential integral for a linear distribution over a triangular domain[END_REF] for instance. Of course, in the implementation of the method, we only use the triangulation Th of ?

∂Ω k and never refer to the exact boundary ∂Ω k , nor to the triangulations T h and Th . The triangulation Th is generated using standard commercial meshing softwares, see Figure 2.

The only approximation in the method stands in the way the surface ∂Ω k is approached by the surface ?

∂Ω k . It is obvious that the approximation of the surface ∂Ω k becomes more and more accurate as the number of triangles increases, and therefore, that I k tends to Ĩk when h tends to zero. In the next section, we investigate the way I k converges to Ĩk , that is to say how the approximation error is controlled by the size of the triangles (or equivalently by the number of triangles).

Error Estimate for the Method

In order to study the approximation error in computing the magnetic field perturbation with the method presented in the previous section, we follow the general method presented in [START_REF] Chien | Numerical evaluation of surface integrals in three dimensions[END_REF]. A similar error analysis is presented in [START_REF] Balac | Induced magnetic field computations using a boundary integral formulation[END_REF] when the surface is approached using piecewise quadratic interpolation. The approximation error related to the computation of the integral over the interface ∂Ω k is given by

E k (P ) = ? ? ∂Ω k Φ(P, Q) n dσ Q - ? ? ? ∂Ω k Φ(P, Q) n dσ Q , (21)
where we have set

Φ : (P, Q) ∈ R 3 ? -→ r(P, Q) r 3 (P, Q) • B 0 . ( 22
)
With the notations of Section 3.1, it is clear that the error estimate consists of studying E k (P ) = ?

K∈ Th ? ? ?

F ( K) Φ(P, Q) n dσ Q - ? ? F ( K) Φ(P, Q) n dσ Q ? = ? K∈ Th e K (P ). ( 23
)
Our goal is now to carefully overestimate the local error e K (P ). We introduce the reference element σ defined by σ = {(s, t) ∈ [0, 1] 2 | 0 ≤ s + t ≤ 1} and denote by ρ 1 , ρ 2 , and ρ 3 its vertices, see Figure 3. For all K ∈ T h with K = F ( K), we define the mapping m K by where L j , j ∈ {1, 2, 3}, denote the j th basis function for the linear Lagrange interpolation over the unit simplex σ. Namely, for all (s, t) ∈ σ, we have

∀ (s, t) ∈ σ, m K (s, t) = F (L 1 (s, t)v 1 + L 2 (s, t)v 2 + L 3 (s, t)v 3 ), (24)
L 1 (s, t) = 1 -(s + t), L 2 (s, t) = s, L 3 (s, t) = t. (25)
We denote by x i (s, t) the i th component of m K (s, t) for i ∈ {1, 2, 3}. A straightforward calculation shows that K is the image of σ by mK where the mapping mK is defined for all (s, t) ∈ σ by

mK (s, t) = F (L 1 (s, t)v 1 + L 2 (s, t)v 2 + L 3 (s, t)v 3 ) = 3 ? j=1 m K (ρ j ) L j (s, t). ( 26
)
If we denote by xi (s, t) the i th component of mK (s, t), we have

xi (s, t) = 3 ? j=1 x i (ρ j ) L j (s, t) = x i (0, 0)L 1 (s, t) + x i (1, 0)L 2 (s, t) + x i (0, 1)L 3 (s, t). (27)
In the following, we refer to the derivation with respect to the variable s by ∂ s or by the subscript s and set Φ P (.) = Φ(P, .). Introducing the change in variables defined by m K and mK , the local error e K (P ) may be written as e K (P ) = ? ?

F ( K) Φ(P, Q) n dσ Q - ? ? F ( K) Φ(P, Q) n dσ Q = ? ? σ (Φ P (m K (s, t)) (∂ s m K (s, t) ∧ ∂ t m K (s, t)) -Φ P ( mK (s, t)) (∂ s mK (s, t) ∧ ∂ t mK (s, t))) ds dt. ( 28
)
The local error is decomposed into two terms e K (P ) = e 1, K (P ) + e 2, K (P ),

where e 1, K (P ) = ? ? (31)

σ Φ P (m K (s, t)) × (∂ s m K (s, t) ∧ ∂ t m K (s, t) -∂ s mK (s, t) ∧ ∂ t mK (s, t)) ds dt, (30)
The first term gives the error when approximating the normal to the surface and the second one the error when computing the integrand on the approximated surface. We start with an estimate on each triangle between the normal to the surface and the normal to the approximated surface.

Lemma 1. For each triangle K ∈ T h and for all (s, t) ∈ σ, we have

∂ s m K (s, t) ∧ ∂ t m K (s, t) -∂ s mK (s, t) ∧ ∂ t mK (s, t) = O ? h 3 ? . Proof. We set e n (s, t) = ∂ s m K (s, t) ∧ ∂ t m K (s, t) -∂ s mK (s, t) ∧ ∂ t mK (s, t). The component i of e n reads e i n (s, t) = x i+1 s (s, t)x i+2 t (s, t)-x i+2 s (s, t)x i+1 t (s, t)-xi+1 s (s, t)x i+2 t (s, t)+ xi+2 s (s, t)x i+1 t (s, t), ( 32
) for i ∈ {1, 2, 3} with the convention i + 2 = 1 if i = 2 and i + 1 = 1, i + 2 = 2 if i = 3.
To get estimates for e n (s, t), we need to develop the functions x i (s, t) and xi (s, t) according to Taylor formula in a neighborhood of (0, 0). It is important to note that the first derivatives of x i behave in O(h), the second derivatives of x i behave in O(h 2 ) and so on. Indeed, the derivatives of x i with respect to (s, t) give rise to formula involving v2 -v1 and v3 -v1 . For example,

∂ s x i (s, t) = ∇F i ((1 -s -t)v 1 + sv 2 + tv 3 ) • (v 2 -v1 ) = O(h), ∂ t x i (s, t) = ∇F i ((1 -s -t)v 1 + sv 2 + tv 3 ) • (v 3 -v1 ) = O(h).
The Taylor expansion of x i (s, t) in a neighborhood of (0, 0) gives x i (s, t) = x i (0, 0) + sx i s (0, 0) + tx i t (0, 0) + 1 2 s 2 x i ss (0, 0) + stx i st (0, 0) +

1 2 t 2 x i tt (0, 0) + O ? h 3 ? . (33)
Then the function xi (s, t) in (27) admits the development xi (s, t) = x i (0, 0)L 1 (s, t) + x i (1, 0)L 2 (s, t) + x i (0, 1)L 3 (s, t)

= x i (0, 0)L 1 (s, t) + ? x i (0, 0) + x i s (0, 0) + 1 2 x i ss (0, 0) + O ? h 3 ? ? L 2 (s, t) + ? x i (0, 0) + x i t (0, 0) + 1 2 x i tt (0, 0) + O ? h 3 ? ? L 3 (s, t) = x i (0, 0) + sx i s (0, 0) + tx i t (0, 0) + 1 2 sx i ss (0, 0) + 1 2 tx i tt (0, 0) + O ? h 3 ? . (34)
We make the difference x i (s, t) -xi (s, t) in ( 33) and (34) to obtain

x i (s, t) = xi (s, t) + H i (s, t) + O ? h 3 ? , (35)
where

H i (s, t) = 1 2
?? s 2s ?

x i ss (0, 0) + 2stx i st (0, 0) + ? t 2t ?

x i tt (0, 0)

? ( 36
)
behaves in O(h 2 ). Equation (35) permits us to express e i n (s, t) as

e i n (s, t) = E i 3 (s, t) + O ? h 4 ? , (37)
where

E i 3 (s, t) = x i+1 s (0, 0)H i+2 t (s, t) + x i+2 t (0, 0)H i+1 s (s, t) -x i+2 s (0, 0)H i+1 t (s, t) -x i+1 t (0, 0)H i+2 s (s, t) (38)
is a collection of term which are of order 3 in h.

The next step in our error analysis uses the notion of pair of symmetric triangles introduced by Chien in [START_REF] Chien | Numerical evaluation of surface integrals in three dimensions[END_REF]. We call a pair of symmetric triangles, two triangles K1 , K2 ∈ Th with vertices v1 , v2 , v3 and v1 , v4 , v5 , respectively, having the following property, see Figure 3,

v1 -v2 = -(v 1 -v4 ), v1 -v3 = -(v 1 -v5 ). ( 39
)
If the triangulation Th of the polygonal domain D (assumed to be a right-angled isosceles triangle) is composed of N = η 2 right-angled, isosceles triangles, then it is shown in [START_REF] Balac | Induced magnetic field computations using a boundary integral formulation[END_REF] that there exists (1/2)(η 2η) pairs of symmetric triangles and η unmatched triangles if η is odd and there exists (1/2)(η 2 -3η + 2) pairs of symmetric triangles and 3η -2 unmatched triangles if η is even. The key result in our error analysis is given in Lemma 2. It states that a cancellation happens among the local error terms from two symmetric triangles.

Lemma 2. Let K1 ∈ Th and K2 ∈ Th be a pair of symmetric triangles, then

E i 3 (s, t) ? ? K1 = -E i 3 (s, t) ? ? K2 , ∀ (s, t) ∈ σ. (40)
Proof. Let v1 , v2 , v3 and v1 , v4 , v5 be the vertices of the two triangles K1 and K2 as shown in Figure 3. Let (a i , b i ) denote the coordinates of vi . We have

a 1 -a 2 = a 4 -a 1 , b 1 -b 2 = b 4 -b 1 , a 1 -a 3 = a 5 -a 1 , b 1 -b 3 = b 5 -b 1 . (41)
Then for all i ∈ {1, 2, 3},

∂ s x i (0, 0) ? ? K1 = (a 3 -a 1 ) ∂ 1 F i (a 1 , b 1 ) + (b 3 -b 1 ) ∂ 2 F i (a 1 , b 1 ) = -(a 5 -a 1 ) ∂ 1 F i (a 1 , b 1 ) -(b 5 -b 1 ) ∂ 2 F i (a 1 , b 1 ) = -∂ s x i (0, 0) ? ? K2 . (42)
A similar development leads to ∂ t x i (0, 0)| K1 = -∂ t x i (0, 0)| K2 and

∂ 2 ss x i (0, 0) ? ? K1 = ∂ 2 ss x i (0, 0) ? ? K2 , ∂ 2 st x i (0, 0) ? ? K1 = ∂ 2 st x i (0, 0) ? ? K2 , (43)
∂ 2 tt x i (0, 0)

? ? K1 = ∂ 2 tt x i (0, 0) ? ? K2 . ( 44
)
This last result combined with (36) shows that

H i 3 (s, t)| K1 = H i 3 (s, t)| K2 .
The result given in the lemma is then obvious from (38) and ( 43),(44).

Lemma 3. The global error E 1 (P ) =

? K∈ Th e 1, K (P ) behaves in O(h 2 ). Proof. The local error e 1, K on the element K is given by (30). A Taylor expansion of Φ P (m K (s, t)) at (s, t) = (0, 0) leads to

Φ P (m K (s, t)) = Φ P (m K (0, 0)) + O(h). ( 45
)
Using Lemma 1 and equation (37), we get the following expression for the i th component of e 1, K :

e i 1, K (P ) = ? σ Φ P (m K (0, 0)) E i 3 (s, t) ds dt + O ? h 4 ? . ( 46
)
Due to pairs of symmetric triangles, cancellation will occur among the terms contributing to the error. Let us consider a pair of symmetric triangles K1 and K2 and their corresponding images K 1 and K 2 in T h . We have Φ P (m K1 (0, 0)) = Φ P (m K2 (0, 0)) and from Lemma 2,

E i 3 (s, t)| K1 = -E i 3 (s,
+ e i 1, K2 = O ? h 4 ? . ( 48
)
To conclude, we use the assumption made on the structured triangulation of D. We split the global error E 1 into two sums. One corresponding to the pairs of symmetric triangles for which the local error behaves in O(h 4 ). The number of pairs of symmetric triangles is proportional to K∈ Th e 2, K (P ) behaves in O(h 2 ). Proof. The local error e 2, K on the element K is given by (30). A Taylor expansion of Φ P (m K (s, t)) at (s, t) = (0, 0) leads to

η 2 = N
Φ P (m K (s, t)) = Φ P (m K (0, 0)) + s ∇Φ P (m K (0, 0)) •∂ s m K (0, 0) + t ∇Φ P (m K (0, 0)) • ∂ t m K (0, 0) + O ? h 2 ? , (50)
whereas a Taylor expansion of Φ P ( mK (s, t)) at (s, t) = (0, 0) leads to Φ P ( mK (s, t)) = Φ P ( mK (0, 0)) + s ∇Φ P ( mK (0, 0))

•∂ s mK (0, 0) + t∇Φ P ( mK (0, 0))

• ∂ t mK (0, 0) + O ? h 2 ? . ( 51
)
We have mK (0, 0) = m K (0, 0) and

∂ s m K (0, 0) = ∂ s mK (0, 0) + O ? h 2 ? , ∂ t m K (0, 0) = ∂ t mK (0, 0) + O ? h 2 ? . It follows that Φ P (m K (s, t)) -Φ P ( mK (s, t)) = O(h 2 ). As well, ∂ s mK (s, t) ∧ ∂ t mK (s, t) = ∂ s mK (0, 0) ∧ ∂ t mK (0, 0) + O ? h 3 ? = O ? h 2 ? . (52)
Thus, the local error e 2, K behaves in O(h 4 ) and the sum over the N = O(h 2 ) triangles of the triangulation, gives

E 2 (P ) = O(h 2 ).
Collecting the results of both Lemmas 3 and 4, we easily deduce Proposition 1.

Proposition 1. The error E k (P ) in approximating

I k (P ) = ? ? ∂Ω k ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q , by Ĩk (P ) = ? ? ? ∂Ω k ? r(P, Q) r 3 (P, Q) • B 0 ? n dσ Q as presented in Section 3.1 behaves in O(h 2 )
where h is the mesh size of the triangulation of ? ∂Ω k . So far we have fixed the point P and got estimates for the error E k (P ). We know that the integrand Φ P behaves like 1/d 2 where d is the distance from P to the interface ∂Ω k . Thus, if the point P is at a distance O(h) from the interface ∂Ω k , then E k (P ) = O(1). That means that for a given triangulation of the interface ∂Ω k , the error is not proportional to h 2 anymore when P is very close to the interface. This lack of accuracy of the method for points very close to the interface is not so harmful since these points are in an area where the approximation of the interface ∂Ω k by ?

∂Ω k leads to uncertainty on their belonging to the domain Ω k .

RESULTS

In this section, we first validate the accuracy of the method and illustrate some of its properties in the case where the sample is a ball made of diamagnetic or paramagnetic material. For this simple geometry, an analytical expression for the induced magnetic field is available, see [START_REF] Jackson | Classical Electrodynamics, Third Edition[END_REF] for instance. As well, we illustrate the convergence rate of the method by varying the mesh resolution of the ball boundary and give results on the computation time obtained on a single computer and on a computer cluster. Finally, we present results obtained on a hollow ball and on a human head.

Validation of the Method

The proposed method is validated with the case of a ball for which the exact values of the induced magnetic field is known, see [START_REF] Jackson | Classical Electrodynamics, Third Edition[END_REF]. The ball has a radius of 3 cm and is embedded in a magnetic field B 0 of 1 Tesla. The ball boundary is meshed using 10080 triangles. We have tested two different materials for the ball: one with diamagnetic properties (χ m = -9.05 10 -6 ) and the other one with paramagnetic properties (χ m = 10 -4 ). We have computed the magnetic field perturbation along a scan line oriented perpendicularly to the static field B 0 and passing through the center of the ball. We observe in the two cases that the numerical solutions and the exact one coincide. For instance, Figure 4 shows the component of the induced magnetic field along B 0 in the case of the paramagnetic ball. The error between the numerically computed field and the exact one rapidly decreases as the distance from the boundary increases. The maximal error is located on the vicinity of the object boundary but is lower than 0.6 percent of the computed field for both materials: 0.018 µT for the diamagnetic material and 0.19 µT for the paramagnetic one (the error is depicted in Figure 5). First, we can notice that the relative error does not depend on the magnetic material involved since our method does not used any approximation nor assumption on the susceptibility range of values. Furthermore, the only approximation of the method lies in the way the interfaces between domains of different susceptibilities is approached using the mesh. This explains why the error is larger close to the boundary than far away from it. The defect in the surface approximation is greater for a point close to the boundary than for a point at a large distance of it. This phenomenon is in complete accordance with the analysis of the error behavior for our method presented in Section 3.2.

Convergence Rate

To illustrate the behavior of the method as the number of triangles increases, we consider a point P at different distances from the center of the ball. The magnetic field perturbation B ? is computed with meshes having an increasing number of triangles from 48 to 10080. We calculate in each case the error between the exact and the numerically computed magnetic field perturbation B ? . In Figures 6 and7, we present the logarithm of the error E versus the logarithm of the number N of triangles for points chosen outside the ball (see Figure 6) and points chosen inside the ball (see Figure 7). In the two cases, the error is proportional to 1 / N or equivalently to h 2 . This is in complete agreement with the results of the error analysis given in Section 3.2. One can see from Figure 6 that the constant of proportionality between the error and 1/N varies as a function of the distance to the sphere. This is in connection with the fact that the intensity of the magnetic field perturbation B ? rapidly decreases as the distance to the sphere increases. On the contrary, the constants of proportionality are the same for the three points chosen inside the ball since the magnetic field perturbation B ? is constant inside the ball. As mentioned in the error analysis, the integrand Φ P behaves like 1/d 2 where d is the distance from P to the interface ∂Ω k . The consequence is a lesser accuracy of the method for points at a distance of the boundary comparable to h. This explains in Figure 7 the behavior of the curve for the point the closest to the boundary for small values of N . Such points inside the ball and closed to the boundary do not "view" the boundary approximation with the same accuracy than the outside points or or than the inside points far from the boundary. As a consequence, they are more sensible to the surface approximation. However, even if the error is more important for small values of N , it lessens quickly as N increases. Finally, these results illustrate the fact that the error in our method depends only on the surface approximation and decreases linearly with the number of triangles. The numerical precision of the method can be determined in advance when choosing the number of triangles describing the boundaries.

Computation Time and Parallelization

Equation [START_REF] Okon | The potential integral for a linear distribution over a triangular domain[END_REF] clearly shows that the computation of the magnetic field perturbation B ? in a point P is obtained by a summation over the triangles of the triangulation. As a consequence, the computation time τ p needed to evaluate B ? in one point depends linearly of the time τ K required to compute the integral over one triangle in the following way:

τ p = N × τ K , ( 53
)
where N is the number of triangles in the triangulation of the boundaries. The total time T needed to compute the magnetic field perturbation B ? in a volume of interest involving (X × Y × Z) voxels is given by This linear dependence of the computation time on the number of triangles is illustrated in Figure 8, where we give the time required to compute the magnetic field perturbation in 256 points for different meshes with an increasing number of triangles. Since the magnetic field perturbation is computed by summing the contribution of all the triangles of the triangulation, the computation algorithm is simple to parallelize with a "divide and conquer" parallelization scheme. We have implemented a parallel version of the computation method on an 8 PC cluster using MPI (www-unix.mcs.anl.gov/mpi). Each node in the cluster receives one eighth of the triangles and computes their contributions (i.e., the related integrals) at each point of the volume of interest. At the end of the process, one node gets all the contributions and sums them. Using this simple parallelization process divides the computation time by a factor equal to the number of computation nodes.

T = X × Y × Z × N × τ K . (54)

Test Example of a Hollow Ball

In [START_REF] Chauvel | Cancellation of metal-induced MRI artifacts with dual component paramagnetic and diamagnetic material: Mathematical modelization and experimental verification[END_REF], the authors show that a judicious association of selected paramagnetic and diamagnetic material in a hollow ball can completely cancel the induced magnetic field outside the ball. They propose this dual component approach as a way of reducing artifacts generated by implanted metallic objects. We have tested our method in the case of a hollow ball with the values given in [START_REF] Chauvel | Cancellation of metal-induced MRI artifacts with dual component paramagnetic and diamagnetic material: Mathematical modelization and experimental verification[END_REF]. A spherical shell of titanium (with external radius of 3.6 cm and magnetic susceptibility χ m = 23 10 -4 ) encapsulates a smaller ball of bismuth (with radius of 3.1 cm and susceptibility χ m = -11.9 10 -4 ). This object is magnetized by a static magnetic field B 0 of 1 Tesla. Figure 9 shows the component of the magnetic field perturbation along B 0 on a central slice parallel to B 0 obtained by our computational method when each of the spherical surfaces is meshed with 966 triangles. One can note that outside the ball, the field cancels despite the high magnetic field observed inside the ball. These results are in concordance with those given in [START_REF] Chauvel | Cancellation of metal-induced MRI artifacts with dual component paramagnetic and diamagnetic material: Mathematical modelization and experimental verification[END_REF] and contribute to validate our computation algorithm with two different surfaces corresponding to an interface air/titanium and an interface titanium/bismuth and involving paramagnetic and diamagnetic materials.

Test on a Human Head

In a last numerical experiment, we have applied our method to compute the field induced by the susceptibility variations between the air and a human head. We use the Montreal brain phantom [START_REF] Collins | Design and construction of a realistic digital brain phantom[END_REF] and the Marching cube algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF] to get a mesh description by 31064 triangles of the interface head/air including nasal fossa and ears, see Figure 2. Using a single Pentium IV 2.4 Ghz, the computation of the magnetic field perturbation on a slice of 256× 256 points requires 20 hours. We present the magnetic field perturbation obtained on a sagittal slice (Figure 10), on a frontal slice (Figure 11), and on an axial slice (Figure 12). These figures clearly show the variation of the magnetic field inside the head due to the high magnetic susceptibility variation between the air in the nasal fossa and the surrounding tissue.

CONCLUSION

In this paper, we have proposed an analytic method to compute MRI susceptibility induced magnetic field inhomogeneities. Our approach is based on a description of the susceptibility interfaces. It does not make any assumption on the range of the susceptibility values or on the direction of the static field. The easy parallelization of the method makes the proposed approach usable on a large data set. Furthermore, since the surface integral involved in the perturbed field computation is computed exactly, the perturbed field could be estimated with a given precision by adjusting only the mesh of the involved surfaces. Such a property has to be explored in further works to define a multiresolution method based on multiresolution meshes, see [START_REF] Valette | A multiresolution wavelet scheme for irregularly subdivided 3D triangle mesh[END_REF][START_REF] Lounsbery | Multiresolution analysis for surfaces of arbitrary topological type[END_REF]. For instance, to compute the perturbed field for points far from the interfaces, a coarse mesh could be used and inversely a fine version of the mesh for points close to the interface.

  Ki |P -Q|, with | | being the usual Euclidian norm. The parameters h and N are linked by the relation N = O(1/h 2 ). Without loss of generality, we can assume that the polygonal domain D is a right-angled, isosceles triangle and that the triangulation Th of D is composed of N rightangled, isosceles triangles. A triangulation T h of the surface ∂Ω k can be obtained from the triangulation Th of the polygonal domain D through the parameterization F , each element K in T h (with vertices v 1 , v 2 , v 3 ) being the image of a triangle K in Th (with vertices v1 , v2 , v3 such that

Figure 1 .

 1 Figure 1. The situation under consideration.

Figure 2 .

 2 Figure 2. Interface between a head and air defined by a mesh of 31064 triangles.

Figure 3 .

 3 Figure 3. Notation used in the error analysis.

  (m K (s, t)) -Φ P ( mK (s, t))) × (∂ s mK (s, t) ∧ ∂ t mK (s, t)) ds dt.

Figure 4 .

 4 Figure 4. Magnetic field perturbation along a scan line oriented in a direction orthogonal to the magnetic field B 0 (1 T.) and passing through the center of a paramagnetic ball (χm = 4 10 -3 ) of radius 3 cm. The numerical computation is done using 10080 triangles for the mesh of the ball surface.

Figure 5 .

 5 Figure 5. Difference between the exact and the numerically computed magnetic field perturbation in the situation of Figure4.

Figure 6 .

 6 Figure 6. Logarithm of error between the exact and the numerically computed magnetic field perturbation as a function of the logarithm of the number of triangles for points outside the ball and located at different distances from the center of the ball of radius R.

Figure 7 .

 7 Figure 7. Logarithm of error between the exact and the numerically computed magnetic field perturbation as a function of the logarithm of the number of triangles for points inside the sphere located at different distances from the center of the sphere of radius R.

Figure 8 .

 8 Figure 8. Computation time of the magnetic field perturbation in 256 points on Pentium IV 2.4 Ghz computer as a function of the number of triangles.

Figure 9 .Figure 10 .

 910 Figure 9. Induced magnetic field (in µT ) on a central slice parallel to B 0 for a hollow ball associating titanium (χm = 23 10 -4 ) and bismuth (χm = -11.9 10 -4 .)

Figure 11 .

 11 Figure 11. Induced magnetic field in µT on an axial slice of the brain crossing the nasal fossa.

Figure 12 .

 12 Figure 12. Induced magnetic field in µT on a frontal slice of the brain crossing the internal ears.

  t)| K2 . It follows that the quantity ?

	σ	Φ P (m K (0, 0)) E i 3 (s, t) ds dt	(47)
	from two symmetric triangles cancels giving
		e i 1, K1

  and behaves in O(1/h 2 ). Thus, this sum over the pairs of symmetric triangles behaves in O(h 2 ). The second sum corresponds to the triangles that cannot be matched for which the local error behaves in O(h 3 ). The numbers of such triangles is proportional to η = √ N and behaves in O(1/h). Thus, this sum over the unmatched triangles behaves in O(h 2 ) as well. In other terms,

	we write	E 1 (P ) =	?	e 1, K (P ) =	?	e 1, K (P ) +	?	e 1, K (P )
		K∈ Th = O ? η 2 ? = O ? h -2 ? O O ? ? h 4 ? + O(η) O matched triangles h 4 ? + O ? h -1 ? ?	O h 3 ? ?	unmatched triangles h 3 ? = O ? . h 2 ?	(49)

Lemma 4. The global error E 2 (P ) = ?

This work is in the scope of the scientific topics of the PRC-GdR ISIS Research Group of the French National Center for Scientific Research CNRS. This work has been also funded by the INSA Lyon. We want to thank B. Belaroussi for his contribution on the hollow ball results and F. Bellet for his help on the parallel implementation of the method.