Maureen Clerc

Juliette Leblond

Jean-Paul Marmorat

Théodore Papadopoulo

Domaines

J.-P Marmorat

J.-P Marmorat

Source localization using rational approximation on plane sections

Keywords: inverse EEG problem, cortical mapping, source localization, rational approximation problème inverse en Électroencéphalographie, localisation de sources, approximation rationnelle

Introduction

Electromagnetic Source Mapping aims at localizing active sources within the brain from measurements of the electromagnetic eld they produce, which can be measured passively outside the head. This article deals more specically with the electric potential which is measured using Electroencephalography (EEG).

Estimating neural current sources located within the bra1in from outside measurements falls into a category of inverse source problems, that are severely ill-posed in general, mainly due to the lack of continuity and stability, but also to non-uniqueness [START_REF] Isakov | Inverse problems for partial dierential equations[END_REF][START_REF] Vessella | Locations and strengths of point sources: stability estimates[END_REF].

When a limited number of sources are modeled as pointwise and dipolar, there are in general more measurements than unknowns, and it has been proved that the inverse problem of source estimation has a unique solution [START_REF] Badia | An inverse source problem in potential analysis[END_REF]]. However, even in this pointwise and dipolar case, solutions to the inverse source problem are often unstable, in particular with respect to the number of sources.

Several families of methods exist to solve the inverse source localization problem, when sources can be modeled as the superposition of a small number of dipoles [START_REF] Scherg | Multiple source analysis of interictal spikes: goals, requirements, clinical value[END_REF]. Dipole Fitting methods must minimise a nonconvex goal function, yielding an outcome that is unstable with respect to the number of dipoles in the model [START_REF] Cun | A method for localizing EEG sources in realistic head models[END_REF]. Whenever this number is assumed to be known a priori, an algebraic (scanning) method has been proposed in [START_REF] Badia | An inverse source problem in potential analysis[END_REF], which requires rank computation of related matrices. In practise, one does not know this number in advance, and learning this model order is far from trivial [START_REF] Bénar | Statistical maps for EEG dipolar source localization[END_REF]. If the sources are decorrelated in time, analyzing the covariance matrix of the measured data provides an estimate of the number of active dipoles. The number of sources is indeed equal to the number of singular values that are signicantly dierent from zero.

The MUSIC Method rst applies a principal component analysis (PCA) to the measurements, identies a signal subspace of which the analysis subsequently determines the dipole positions [START_REF] Mosher | Multiple dipole modeling and localization from spatio-temporal meg data[END_REF]. In practise, the dimensionality of the signal subspace is dicult to determine, and the dipoles are extracted one at a time by seeking the global maximum of a contrast function among all possible source positions. MUSIC can thus only be applied if the sources are well modeled by a small number of asynchronous dipoles. With the stronger assumption of decorrelated sources, another method, Beamforming, can also estimate active sources, by scanning a region of interest, and by comparing the covariance of the measurement to that of the baseline, measured in time windows that do not contain the activity of interest [START_REF] Van Veen | Beamforming: a versatile approach to spatial ltering[END_REF]. This article proposes a new approach, which, as MUSIC and Beamforming, requires no prior information on the number of sources. However, unlike MUSIC or Beamforming, which require as input consecutive measurements within a time window, the proposed method works instant by instant, and a fortiori, does not require sources to be decorrelated across time.

Our method belongs to a new category of source estimation algorithms that are grounded in Harmonic Analysis and Best Approximation theory, and oer stability [START_REF] Baratchart | Inverse source problem in a 3D ball from best meromorphic approximation on 2D slices[END_REF][START_REF] Kandaswamy | Analytic sensing: Noniterative retrieval of point sources from boundary measurements[END_REF]. These analytical methods directly localize the sources as the singularities of the potential from boundary measurements.

RR n°7704

The type of geometry and of boundary data which they require do not necessarily coincide with actual measurements. Indeed, these methods usually work in a homogeneous domain, and an explicit parametrization may be needed for its boundary.

We present our constructive approach to this inverse problem, in the framework of a (classical) spherical geometry. Under quasi-static assumptions, Maxwell's equations lead to a formulation of the electric potential u as a solution to Laplace's equation. In the innermost layer (the brain), there may be singularities due to the presence of current sources. These singularities are to be localized from available data on the outer boundary (the scalp). The core of our inverse solution relies on approximation schemes that are meant to operate inside homogeneous domains, and not in such a nested geometry. Hence, a preliminary stage consists of mapping to the cortex the data initially measured on the surface of the scalp. This Cortical Mapping problem is a Cauchy (transmission) problem for the Laplace operator [START_REF] Clerc | Cortical mapping by laplace-cauchy transmission using a boundary element method[END_REF].

Source detection from cortical data is another classical inverse problem for the Laplace operator, that consists of recovering an unknown number of pointwise sources within a homogeneous domain from measurements of the potential and its normal derivative on the boundary.

When the domain is a spherical ball, then the above issue is equivalent to a sequence of 2D inverse problems, each of which consists of recovering the singularities of some function f in a disk from the knowledge of f on the boundary circle [START_REF] Baratchart | Inverse source problem in a 3D ball from best meromorphic approximation on 2D slices[END_REF]. Consequently, we apply to these 2D problems a technique inspired by that described in [START_REF] Baratchart | Recovery of pointwise sources or small inclusions in 2D domains and rational approximation[END_REF] that relies on approximating f on the boundary circle by a rational function with poles in the disk.

Finally, we locate the singularities by analyzing the cluster of these poles. Geometrical restriction to spherical domains allows one to make explicit (and not too complex, in a preliminary feasibility study) the behaviour of the 2D singularities with respect to the 3D sources, which is granted by our recovery scheme.

The outline is as follows: Section 2 introduces the inverse problem, and Section 3 presents the solution proposed in this paper. Section 4 demonstrates the method on numerical examples. The paper also includes a conclusion in Section 5 and technical Appendices that detail some mathematical aspects of the method.

The inverse problem 2.1 Model setting

In a simplied spherical model, the head is assumed to be the union of three disjoint homogeneous spherical layers1 Ω 0 , Ω 1 , Ω 2 , namely the brain, the skull, and the scalp, within a non-conductive medium Ω 3 representing the air. Up to a rescaling, one may assume the ball Ω 0 representing the brain to have radius 1 and to be centered at the origin. The spheres separating the volumes Ω i are denoted S 0 , S 1 and S 2 (see Figure 1). The conductivity in each Ω i is denoted σ i . For simplicity and without loss of generality, we assume that σ 0 = 1. Then, we dene a piecewise constant function σ in R 3 by σ |Ωi = σ i . The current

Ω 0 Ω 1 Ω 2 Ω 3 S 0 S 1 S 2

The inverse problem

The potential created by the dipolar sources (C k , p k ) located inside Ω 0 is a solution to the forward problem

(FP)      ∇ • (σ∇u) = S = n k=1 p k . ∇ δ C k in R 3 σ∂ ν u | S 2 = 0 (current ux) (1)
where ν denotes the outward unit normal vector to the surfaces. The homogeneous Neumann boundary condition is due to the fact that the outer medium (air) is non-conductive. The current owing through the neck is neglected.

Let K denote a set of points on S 2 , representing electrode positions. The inverse source localization problem (IP) associated to the forward problem (FP) is then the following:

(IP) Given measurements of u on K, nd the number of unknown pointwise dipolar sources, their positions C k ∈ Ω 0 and their moments p k ∈ R 3 , such that u satises (FP).

Properties of solutions to (IP)

Mathematical properties of (IP) have been established when the data is known in an open subset K of the boundary. In our case, the data is only known on a discrete set, but it is assumed that the underlying potential is smooth enough so that it is well approximated on an open set K from the data. Uniqueness of solutions to (IP) (identiability from boundary measurements) has been established in [START_REF] Badia | An inverse source problem in potential analysis[END_REF]. If two nite distributions of pointwise dipolar sources generate the same potential on some open subset K of S 2 , then they are identical.

RR n°7704

Stability properties for (IP) may only hold for smooth enough boundary data on K (in Sobolev spaces). This is due to the classical ill-posedness of Cauchy-type issues (Section 3.1.1) [START_REF] Alessandrini | The stability for the cauchy problem for elliptic equations[END_REF][START_REF] Isakov | Inverse problems for partial dierential equations[END_REF][START_REF] Phung | Remarques sur l'observabilité pour l'équation de laplace[END_REF], Zghal, 2010]. Concerning the source problem in homogeneous domain (sections 3.3 and 3.4), stability results are established in [START_REF] Badia | Inverse source problem in an anisotropic medium by boundary measurements[END_REF][START_REF] Vessella | Locations and strengths of point sources: stability estimates[END_REF].

3 Solution to the inverse problem (IP)

The resolution of the above inverse problem (IP) consists of two main steps as represented in the owchart in Figure 2.

Data Transmission from S 2 to S 0 , which involves:

Cortical Mapping (Section 3.1.1): the data is transmitted from the surface of the scalp S 2 where it is measured (on electrodes) onto the surface S 0 of the brain.

Harmonic Projection (Section 3.1.2): ltering out possible outer sources by keeping only the information related to the eective inner sources in

Ω 0 .
Source Recovery in Ω 0 from data on S 0 , which involves:

Plane Sections (Section 3.2): the sphere S 0 is sliced along families of parallel planes, perpendicular to a chosen axis, yielding disks inside which the singularities will be sought.

Planar Singularity Detection (Section 3.3): 2D approximation tech- niques are used to nd the planar singularities on the plane sections of Ω 0 (disks).

3D Source Localization (Section 3.4): for a putative number of sources, the sources are localized in 3D by analyzing the sets of planar singularities.

The Data Transmission step uses the approach proposed in [START_REF] Clerc | Cortical mapping by laplace-cauchy transmission using a boundary element method[END_REF] while the Source Recovery step uses the one described in [START_REF] Baratchart | Inverse source problem in a 3D ball from best meromorphic approximation on 2D slices[END_REF] (see also [START_REF] Ben Abda | Sources recovery from boundary data: a model related to electroencephalography[END_REF], in two dimensions).

Data Transmission

Cortical Mapping

The goal, as recalled in Figure 2, is to estimate the values of the potential and the normal current on the cortex, from the values on electrodes of a potential that satises the Forward Problem (1). This Forward Problem can be decomposed in each of the three layers Ω i , i = 0, 1, 2. By assumption, there are no sources outside the inner volume Ω 0 , hence the potential u satises a homogeneous Laplace equation in the layers Ω 1 and Ω 2 : The continuity of the potential and of the normal current across the interfaces are expressed through the following transmission conditions:

∆ u = 0 in Ω i , i = 1, 2. RR n°7704
u + = u -on S i , i = 0, 1, (2)
and

σ i+1 ∂ ν u + = σ i ∂ ν u -on S i , i = 0, 1, (3)
where superscripts + andindicate the limiting values when approaching S i from Ω i+1 (outside) and Ω i (inside), respectively.

The Data Transmission problem which we aim to solve is a Cauchy Problem (CP), for u harmonic within Ω 1 and Ω 2 , satisfying the transmission conditions

(2-3):

(CP) Given measurements u K of u on K ⊂ S 2 and given that ∂ ν u = 0 on S 2 , nd the values g = u |S0 and φ = ∂ ν u - |S0 on S 0 .

Solving (CP) is non-trivial, because the Cauchy problem for the Laplace equation is the prototype of an ill-posed problem. Indeed, (CP) has similar stability properties to those described in Section 2.3 for K ⊂ S 0 . Regularization schemes have been proposed in [START_REF] Atfeh | Bounded extremal and Cauchy-Laplace problems on the sphere and shell[END_REF][START_REF] Kozlov | An iterative method for solving the Cauchy problem for elliptic equations[END_REF]. We solve the Cortical Mapping with a regularized Tikhonov method deriving from a boundary elements formulation of the problem (this step is thus not limited to spherical interfaces): from its values u K on the measurement set K, u is estimated, along with σ∂ ν u, on the three surfaces S 2 , S 1 , and S 0 . This method, originally presented in [START_REF] Clerc | Cortical mapping by laplace-cauchy transmission using a boundary element method[END_REF], is detailed in A.

3.1.2 Harmonic Projection (in R 3 )

After the Cortical Mapping step has provided the potential g = u and its normal derivative φ = ∂ ν u on the surface of the cortex S 0 , the potential u satises in

RR n°7704

Ω 0 an equation of the form:

     ∆ u = S = n k=1 p k . ∇ δ C k in Ω 0 , ∂ ν u = φ , u = g on S 0 . (4)
From g and φ on S 0 , we rst look for the part u a of the potential u which is harmonic outside the ball Ω 0 , vanishes at ∞, and still contains on S 0 all the information on the distribution of sources.

Knowing that the potential u in Ω 0 is solution of ( 4), let u a be the convolution of S with the Green function for the Laplacian in R 3 :

u a (x) = n k=1 p k , x -C k 4π |x -C k | 3 , x = C k , (5)
where the brackets , denote the scalar product in R 3 × R 3 . Note that:

∆ u a = ∆ u = S in Ω 0 , while ∆ u a = 0 in R 3 \ Ω 0 , and lim |x|→∞ |u a (x)| = 0 . Consequently, for x ∈ Ω 0 \ {C k }, we have: u(x) = h(x) + u a (x) = h(x) + n k=1 p k , x -C k 4π |x -C k | 3 ,
for a harmonic function h in Ω 0 .

In practise, u a is computed from the available boundary data g and φ on S 0 by expanding u there on the basis of spherical harmonics [Baratchart et al., 2006, Dautray and[START_REF] Rr | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF]. Indeed, u being a harmonic function in a neighbourhood of S 0 , the coecients of its expansion on S 0 of negative indices coincide with those of u a . They are given through a linear system, by identication with the coecients of the spherical harmonics expansions of the discretized g and φ on S 0 . Figure 3 shows the singular part u a of the potential, computed from its expression ( 5) on S 0 , from n = 2 sources C 1 = (0.5, 0.5, 0.5), C 2 = (0.5, -0.5, -0.4) with moments p 1 = (1, 1, 1), p 2 = (-1, 1, 1).

From 3D to 2D

Given u a , we can now formulate an inverse source recovery problem in Ω 0 : (SP) Given u a on S 0 , nd the number of unknown pointwise dipolar sources, their locations C k ∈ Ω 0 and their moments p k ∈ R 3 , such that u a satises [START_REF]3 Solution to the inverse problem (IP) 6[END_REF].

Following [START_REF] Baratchart | Inverse source problem in a 3D ball from best meromorphic approximation on 2D slices[END_REF], we rst study the singularities of u 2 a on plane sections of Ω 0 , where the function can be analytically extended to the complex plane.

The ball Ω 0 is sliced along a family of P planes, Π p , p = 1, • • • , P , parallel to some plane Π ⊂ R

Π p = {x 3 = x 3,p }. For some xed x 3,p ∈ (-1, 1), let x ∈ T p : x = (x 1 , x 2 , x 3,p ) , r p = 1 -x 2 3,p , z = (x 1 + i x 2 )/r p , (6)
where r p is the radius of the circle T p , and z ∈ T is the normalised complex ax associated to x ∈ T p .

RR n°7704

From u a on T p , we build the complex variable function f p on the unit circle T ⊂ C (the complex plane) as follows.

f p (z) = u 2 a (x 1 , x 2 , x 3,p ) , (7)
For xed p, the function f p coincides with the trace on T of a function dened on C except at singularities: due to the n sources C k , this extended function (that we still call f p ) has n singularities inside the unit disk D (as well as n related singularities outside the closed disk D).

Indeed, let us denote the source coordinates by C k = (x 1,k , x 2,k , x 3,k ), and their corresponding complex ax by z k = x 1,k + ix 2,k , for k = 1, . . . , n. Assuming z k = 02 , we have from ( 5), ( 7), at x ∈ T p and corresponding complex ax z ∈ T through (6):

f p (z) = n k=1 ϕ k,p (z) (z -s k,p ) 3/2 2 (8)
where s k,p are the singularities induced inside D by the source C k and ϕ k,p are smooth functions in D. The localization of s k,p then leads to that of C k . Indeed, the complex argument of s k,p is independent of p, and equal to the argument of z k , which allows us to determine the number n of sources. Further, for xed k, when p varies, the quantity r p |s k,p | attains its maximal value, equal to |z k |, when x 3,k = x 3,p , in the slice p corresponding to C k . Figure 5 shows the trajectories of the singularities s k,p in a two sources case.

Recovering 2D singularities

This section deals with the computation of the singularities s k,p from the sliced boundary data f p on T , given by (7) at xed p. From formula (8), it can be seen that the function f p has the following properties:

If there is a single source (n = 1), then f p is exactly a rational function with a single triple pole in D at position s 1,p :

f p (z) = ϕ 2 1,p (z) (z -s 1,p ) 3 . (9)
If there are multiple singularities (n ≥ 2), f p is no longer a rational function (because of terms with power 3/2 at the denominator of ( 8)). In this situation, the s k,p are both (triple) poles and branchpoints (of order 3/2)3 . Yet, f p can be well approximated on the boundary T by a rational function with poles in D, see C and [START_REF] Baratchart | Inverse source problem in a 3D ball from best meromorphic approximation on 2D slices[END_REF], [START_REF] Baratchart | Meromorphic approximants to complex Cauchy transforms with polar singularities[END_REF].

This gives rise to the following algorithm that provides estimates s k,p of s k,p from sample values f p built from (7). Finding Planar Singularities s k,p s k,p from Cortical Data f p 1. Choose the number n of sources.

Find initial values s *

k,p of s k,p , to be the poles of a rational approximation f * p of f p with appropriate degree (depending on n).

s k,p = arg min s k,p f p -f p .
Remark 1 Though it should become an output of the proposed method, the number n is a necessary preliminary guess in the present algorithm.

In point (iv),

is the l 2 norm on T . The data f p is assumed to be given either by a number of its pointwise values on T , or by a number of its Fourier coecients, using the spherical harmonics expansion of u a on S 0 from Section 3.1.2.

From this Fourier expansion, we shall keep only the part with negative indices, for it is enough to account for the singularities of f p in D (see C.2), that we still call f p , for simplicity.

Strategies for achieving step (ii) of this algorithm are discussed in C.

The present Planar Singularity Detection step must be performed with several slicing directions Π, in order to get more accuracy on the localization process and to separate sources (see Section 3.4 and[START_REF] Marmorat | Matrix rational H 2 approximation: a stage-space approach using Schur parameters[END_REF][START_REF] Marmorat | Rational approximation: RARL2[END_REF]).

From estimated singularities to sources positions and moments

Given a slicing direction Π, the method described in Section 3.3 provides estimates of the singularities s k,p of f p in each slice. But, we know from Section 3.2 that the points (r p s k,p ) k,p are organized along as many lines l k as there are sources (k = 1 . . . n, see the lines formed by green dots in Figure 5).

Sources from a single slicing direction

Each line l k is (generically, for most slicing directions) associated with one of the sources C k and has the following theoretical properties (from section 3.2):

1. l k lies in a single half-plane H k dened by C k and by the direction of the slicing plane Π: H k contains C k and is orthogonal to Π;

2. l k goes through its associated source C k . At this point, the radial distance of the line l k to the boundary of H k (the diameter of S 0 orthogonal to Π)

reaches its maximum.

The rst property is used to group the points r p s k,p into n estimated lines l k .

To do so, these points are clustered in n classes by applying Matlab's algorithm clusterdata to their polar angles. For each class, the best tting half-plane H k is then estimated using a least-squares algorithm and the points are reprojected on that plane, providing us with an estimation l k of the line l k . The polar angle of H k is an approximation of the complex argument Arg z k of Arg z k . We can now compute |z k | using the second above property as:

|z k | = max p∈1,...,P {|r p s k,p | among r p s k,p ∈ l k ⊂ H k } .
We then get z k which, together with the argument p of the above max, provides an estimate C k of the source positions, for direction Π.

Combining information from multiple slicing directions

The above procedure is repeated for a number n p of dierent slicing directions, which yields a family of n p × n estimations of the n source positions. Matlab's algorithm clusterdata is again used to build n separate clusters from these n p × n points. The nal estimations C k of the sources are obtained as the barycenter of each cluster. As in the previous section, the distance between points is dened so as to ensure that each slicing direction contributes only once to each cluster.

Once the source positions are known, the measured potential is a linear function of the moments. These are thus estimated using a simple linear least squares minimization procedure.

Numerical validation

We now present numerical results obtained with FindSources3D, a Matlab code that implements the above algorithm [START_REF] Bassila | Findsources3d[END_REF].

RR n°7704

We simulated two data sets with OpenMEEG, which implements the symmetric Boundary Element Method [START_REF] Gramfort | OpenMEEG: opensource software for quasistatic bioelectromagnetics[END_REF][START_REF] Kybic | A common formalism for the integral formulations of the forward EEG problem[END_REF]. We considered data at the scalp level (potential measured by 128 electrodes), and also at the cortex level (potential and normal current on a 642 points mesh), in order to test the inuence of the Cortical Mapping step on the quality of the source estimation. The spherical 3-layer head model is the one described in Section 2.1.

In the case of cortical data, that is with potential and current at the cortex level as computed by OpenMEEG, Figure 6 displays for twelve dierent slicing directions, the top views of the planar detected singularities. Figure 7 shows the 3D superposition of all these estimated singularities for all axis directions and the estimated source positions for all the slicing directions.

Figure 8 displays true and estimated sources from the two datasets, while Table 1 displays numerical values of the corresponding positions and moments. For a given slicing direction Π, the singularities associated with a source C k lie in a plane H k containing the source itself and the slicing axis (Figure 6).

As a consequence, singularity lines associated to various slicing directions Π intersect at the sources, which allows to estimate their positions (Fig-

ure 7).

RR n°7704

Once these positions are estimated, the linear problem s easily solved to recover the source moments (Figure 8).

As could be expected, estimation is better when data are directly taken at the cortical level on many points. The treatment of more realistic datasets (from scalp electrodes) needs an additional Cortical Mapping step and thus achieves a less precise estimation. However the full procedure proves to be ecient enough:

when estimating sources from the electrodes dataset, the global position error is less than 10% of the sphere radius (order of the cm).

Conclusion

We presented here some insights concerning the resolution of a source estimation problem. The techniques rely on constructive approximation, they are robust and ecient towards the EEG inverse source problem, as is illustrated by preliminary numerics.

More accuracy on source localization may be achieved by extending the present method. A rst possibility would be to take into account several time samples while constraining the source positions to be xed. Also, the computation steps concerning the singular part of the cortical potential and its 3D to 2D transformation could be made more direct, in order to limit the numerical errors.

The number of unknown sources is not yet identied automatically. However, at many steps, information is available to build good estimates of this unknown number: singular values in the rational 2D approximation scheme, residual boundary approximation error, clustering procedure, . . . Work is in progress to make this number an output of the whole process, using techniques such as the Akaike information criterion to decide when increasing the number of sources is no longer signicative with respect to the data.

Magnetic data from MEG (magnetoencephalography) will be incorporated as well, coupled to the available EEG, that may lead to additional precision in the source localization process.

Geometrically, the approach applies in principle to more general smooth 3D

domains [START_REF] Ebenfelt | An inverse problem for the double layer potential[END_REF]], but we did not carry out such generalizations here and only considered spherical models. One may observe that whenever the complexity of the geometry increases, so does the quantity of planar singularities associated to a source. For an ellipsoidal domain, which has been theoretically

studied in [START_REF] Leblond | Source localisation in ellipsoids by best meromorphic approximation in planar sections[END_REF], we already get 2 planar singularities for each source, in each ellipse.

RR n°7704

A Data transmission by Cortical Mapping

The Cortical Mapping Method, originally presented in [START_REF] Clerc | Cortical mapping by laplace-cauchy transmission using a boundary element method[END_REF], proceeds as follows. With the rationale of the symmetric Boundary Element method [START_REF] Kybic | A common formalism for the integral formulations of the forward EEG problem[END_REF], u (resp. the normal current σ∂ ν u) is approximated with continuous, piecewise-linear elements (resp. discontinuous, piecewise-constant elements). The discretization of these two variables on each of the boundaries S i , i = 0, 1, 2 yields a set of values which are combined in a single vector-valued variable X. The harmonic nature of u in Ω 1 ∪ Ω 2 , along with the fact that ∂ ν u = 0 on S 2 and the transmission conditions (2-3) are all handled by saying that X must belong to the kernel of a specic linear operator. This linear operator is represented by a dense matrix H, whose elements involve boundary integral operators. The knowledge of u on K is handled by a measurement operator M , such that M X represents u K , i.e. the measurements on K. Ideally one would like to nd X such that M X = u K and H X = 0, but the ill-posedness of the Cauchy inverse problem makes it necessary to stabilise the system through a regularization. As a consequence, the method seeks X belonging to the kernel of H solving

arg min X∈KerH M X -u K 2 + λ R X 2 .
The norms above are discrete l 2 norms, λ is a real positive Lagrange parameter to be adjusted, and R is an appropriate regularization operator. Once the minimizer X has been computed, it is immediate to extract from X the desired transmitted data u and σ∂ ν u on S 0 . Results obtained by this method are illustrated in Figure 9, where the 2 sources were taken as in Section 3.1.2, and each sphere was meshed with 642 points. This gure shows the propagation of the potential measured on 128 electrodes onto the outer skull surface S 1 and cortical surface S 0 .

Note that other transmission schemes can be obtained by best approximation with harmonic gradients, as in [START_REF] Atfeh | Bounded extremal and Cauchy-Laplace problems on the sphere and shell[END_REF], and robust interpolation issues can be handled using spherical harmonics [START_REF] Rr | Mathematical Analysis and Numerical Methods for Science and Technology[END_REF].

B Link between 3D sources and 2D singularities

Choose a xed slicing direction Π, as in Section 3.2. For x ∈ T p and z ∈ T given by ( 6), let us establish equation (8) for f p , which is equal to u 2 a according to (7), where u a has been dened in [START_REF]3 Solution to the inverse problem (IP) 6[END_REF]. Indeed, with h k,p = x 3,p - - (

x 3,k , |x -C k | 2 = (x 1 -x 1,k ) 2 + (x 2 -x 2,k ) 2 + (x 3,p -x 3,k ) 2 = |r p z -z k | 2 + h 2 k,p , thus: |x -C k | 2 = (r p z -z k )(r p z -zk ) + h 2 k,p = (r p z -z k ) r p z -zk + h 2 k,p because z = 1/z for z ∈ T (use that |z| 2 = z z = 1). Assuming rst that z k = 0, we get (r p z -z k ) r p z -zk + h 2 k,p = - r p zk z z 2 - h 2 k,p + |z k | 2 + r 2 p r p zk z + z k zk .
r p zk z (z -s k,p )(z -σ k,p ) ,
σ k,p = z k zk s k,p and s k,p = z k 2|z k | 2 r p h 2 k,p + r 2 p + |z k | 2 (h 2 k,p + (r p + |z k |) 2 )(h 2 k,p + (r p -|z k |) 2 ) .
) 10
With the standard convention that the square-root is positive for positive arguments (which is used throughout the paper), s k,p is the root with the smallest modulus and (10) ensures that s k,p ∈ D while σ k,p ∈ D. Recalling ( 7) and ( 5), this leads to the expression (8), or equivalently to (11) which shows that f p admits the singularities s k,p in D. See below for the functions ϕ k,p at the numerators. The above computation exposes two useful properties of s k,p which are used in Section 3.4 towards the localization of C k :

f p (z) = n k=1 ϕ 2 k,p (z) (z -s k,p ) 3 + 2 n k,j=1 k =j ϕ k,p (z)ϕ j,p (z) (z -s k,p ) 3/2 (z -s j,p ) 3/2 ,
1. The argument of the complex number s k,p is independent of p, and equal to the argument of z k . In any slicing direction that separates the sources (that is if z k = z j for k = j, which generically holds), this property allows:

to determine the number of sources n, since the quantity of sources should be equal to the number of values taken by the complex argument of s k,p , as k and p vary; for any xed index k 0 , to track s k0,p among all the s k,p in any slice p (the complex argument of s k0,p does not depend on p);

to determine the argument of z k0 .

2. When p varies (for x 3,p ∈ (-1, 1)), the modulus of r p s k,p increases monotonically for x 3,p < x 3,k (decreases monotonically for x 3,k < x 3,p ), and attains a maximum when x 3,k = x 3,p or h k0,p = 0, in which case one has r p |s k0,p | = |z k0 |. This second property allows us to determine |z k0 |, whence nally z k0 .

Also, if we put p k = (p 1,k , p 2,k , p 3,k ) for the moments and k = p 1,k + ip 2,k , we get from ( 5) that:

ϕ k,p (z) = 1 8π ¯ k (r p zk ) 3/2 √ z (σ k,p -z) 3/2 π 2,k,p (z)
are uniformly bounded in D and π 2,k,p are polynomials of degree 2 (see also [START_REF] Baratchart | Inverse source problem in a 3D ball from best meromorphic approximation on 2D slices[END_REF]):

π 2,k,p (z) = r p z 2 + 2 p 3,k h k,p -Re (z k ¯ k ) ¯ k z + k ¯ k r p .
Finally, whenever z k = 0 (that is when the associated C k lies on the vertical axis), then |x -

C k | 2 = r 2 p + h 2
k,p is a constant and the corresponding term of (8) becomes a rational function of z which assumes the form π 2,k,p (z)/z for the above mentioned polynomials π 2,k,p . In this situation, the function f p to be approximated has a unique (double) pole in D at 0, in every planar section p, which will be revealed by the Rational Approximation step. The sum of the RR n°7704 roots of the polynomial π 2,k,p is equal to -

2p 3,k
¯ k rp h k,p . If p 3,k = 0, the behaviour of these roots still allows to compute the index p such that h k,p = 0, and to nally locate the singularity C k . The situation where z k = 0 and p 3,k = 0 however is degenerated w.r.t. the present choice of Π, which is the reason why several slicing directions should be used.

C Best Rational Approximation Schemes C.1 Best Rational Approximation

From the knowledge of f p on the surrounding circle T , the s k,p in D, k = 1, ..., n are localized using Rational Approximation on T , with poles in D. Indeed, as we will see, the poles s * k,p of such approximants accumulate to the singularities s k,p of the approximated function f p .

Let us rst briey explain the best quadratic Rational Approximation techniques that are used.

As explained in Section 3.3, the singularities s k,p can be described both as branchpoints and as triple poles of f p . For a single source, we noted that f p is exactly a degree 3 rational function and the singularity s k,p is itself a triple pole. For multiple sources, the situation is not so simple, but the property that poles lines pass near the singularities still remains. This makes it interesting to consider two Rational Approximation schemes: one with simple poles and the other with triple poles.

We dene R m to be the set of rational functions R m with less than m poles in D: R m = π m q m , where π m and q m are polynomials such that deg π m < deg q m ≤ m, and where the zeroes of q m belong to D.

A best quadratic rational approximant to

f p in R m is a function R * m ∈ R m , verifying: f p -R * m = min Rm∈Rm f p -R m , (12)
for the L 2 (T ) norm, see [START_REF] Baratchart | On a rational approximation problem in the real Hardy space H 2[END_REF]. Existence and non-uniqueness of R * m are discussed in [Baratchart et al., 1992]. Concerning constructive as- pects, ecient algorithms to generate local minima are obtained using Schur parametrization [START_REF] Marmorat | Matrix rational H 2 approximation: a stage-space approach using Schur parameters[END_REF]. Computation of R * m is made eec- tive through suitable parametrization of rational functions, using gradient algorithms [START_REF] Marmorat | Matrix rational H 2 approximation: a stage-space approach using Schur parameters[END_REF].

Generally and properly speaking, for functions f p ∈ L 2 (T ), the approximation class R m should be the set of meromorphic functions with less than m poles in D. Such a meromorphic function is the sum of a rational function πm qm and of a function h holomorphic (or equivalently analytic) in D (which has no singularities in D but may have poles outside D). Because

π m q m + h = π m + hq m q m ,
we see that a meromorphic function coincides with the quotient of a holomorphic function by a polynomial.

Further, whenever f p has singularities only in D and is analytic in C \ D (and vanishes at ∞), its best meromorphic approximant coincides with its best rational approximant with poles in D. This is the reason why we explain in C.2 how RR n°7704

to get, in general, this part of f p , analytic in C\D, which shares the singularities of f p in D but as no singularity outside D.

C.2 2D Analytic Projection

Observe from B that f p possesses singularities both inside and outside the unit disk, that are linked to each other and to the sources C k , see e.g. ( 10). Indeed, the representation (8) involves the singularities s k,p of f p in D but also the additional reected ones σ k,p (hidden in ϕ k,p ) outside D. First, these are linked with each other by (10). Next, the Rational Approximation algorithms available for data on T require the singularities of the approximated functions to belong to a region of the plane limited by T (see Section 3.3;in particular,this leads to a lower degree of the approximants, that also possess strong robustness properties). We thus choose to keep only those of the singularities s k,p of f p that belong to the disk D, and we need to lter the σ k,p out. It is easily seen from ( 8) that the function f p is continuous on T , because none of its singularities belong to T . Thus it belongs to the Lebesgue space L 2 (T ) of functions with square summable modulus on T . Consequently, it can be uniquely written on T as the sum f p = F + F o , where F is the holomorphic projection of f p in C \ D and vanishes at ∞, while F o is holomorphic in D.

Actually, F and F o respectively belong to the Hardy classes of C \ D and of D, their traces on T belonging to L 2 (T ) [START_REF] Rudin | Real and Complex Analysis[END_REF]. The Hardy-Hilbert spaces H 2 -, resp. H 2 + are the sets of functions analytic in C \ D (vanishing at ∞), resp. analytic in D, and bounded in L 2 (T ) norm (i.e. the space of L 2 (T ) functions with vanishing Fourier coecients of positive indices, resp. of strictly negative indices). We can directly compute F from the Fourier series expansion of f p on T :

f p (e iθ ) = l∈Z F l e ilθ , l∈Z |F l | 2 < ∞ ⇒ F (z) = (P -f p )(z) = l<0 F l z l , |z| ≥ 1 , (13)
if P -denotes the orthogonal anti-analytic projection from L 2 (T ) onto H 2 -.

The important point here is that f p and F share the same singularities inside D while F has no singularities outside D, since F possesses an expression analogous to (8), with identical denominators, but numerators given by smooth functions. This is necessary for the best Rational Approximation problem o be solved among rational (no longer meromorphic) functions with poles in D.

We shall then assume from now on, and already for the computations of Section 3.3, that f p is analytic in C \ D, and vanishes at ∞, without loss of generality.

C.3 Behaviour of Simple Poles with respect to Singularities and Sources

If the degree m is not preliminarily given, observe that its estimation can be obtained by computing the boundary error on T (the value of the criterion in ( 12)) for increasing values of m, until it is small enough, see Remark 1. This will happen in principle for m 3 n, the number of singularities of f p in D according to their multiplicity.

RR n°7704

For such a value of m, compute the best approximant R * m itself (that the ap- proximant should be computed only once, for m large enough for the error to be suciently small, is one of the features that make this scheme ecient). Recall that the singularities s k,p that we aim at recovering appear at triple poles of f p , from (8), which motivates the computation of best rational approximants with triple poles.

RR n°7704

Recall that for a single source (n = 1), then ( 8) is to the eect that f p is a rational function with a triple pole in D, see (9). Hence, its best rational approximant with a single triple pole in D should coincide with f p itself. Whenever n > 1, the situation is of course more complicated since f p admits the s k,p both as triple poles and branchpoints. However, as we shall see, the behavior of poles dominates whence its best rational approximant with triple poles still allows to recover the s k,p . These best rational approximants with triple poles in D are functions R * 3m = π 3m /q 3 m ∈ R 3m that satisfy (12), where π and q are polynomials such that deg q m ≤ m, deg π 3m < 3m. An advantage is that the computations can then be performed with a lower degree than in the simple poles case, since m = n is enough.

Again, even though f p is not a rational function (it admits poles and branchpoints located at the same place), it is (numerically) close to rationals of R 3 with a single triple pole (m = 1), in plane sections p close to the one containing (C k ), even for several sources, when n > 1.

For such data f p , close to a rational function of R 3 with a single triple pole, say t, in D, then the single triple pole of the best rational approximant is close to t, see the robustness property (P'), C.6. Hence the above best approximants R * 3m possess an even stronger property: in general a single (m = 1) triple pole is enough in order to localize several singularities s k,p in D, hence several sources (C k ) (n ≥ 1), by varying p. We give a few additional considerations concerning the asymptotic behaviour of the poles.

In situations where f p is already a rational function R N of R N , then its best rational approximant R * m coincides with f p for m ≥ N . This result is robust in the sense of property (P) below. Observe that whenever n ≥ 2, the function f p has poles and branchpoints in D, and it can be shown that the degree of the denominator q m of R * m is in fact equal to m, for each integer m. Property (P) can be deduced from [Baratchart, 1986, Prop. 5]:

(P) Whenever f p is close (in L 2 (T )) to a rational function R N , then the poles of R * m accumulate to those of R N , as m increases.

Further, deep convergence results from potential theory [START_REF] Baratchart | Meromorphic approximants to complex Cauchy transforms with polar singularities[END_REF] assert that, for a function f p as in the present situation (which admits nitely many poles and branchpoints in D and has a smooth behaviour near T ), the m poles of R * m converge (in some weak sense) to the singularities (s k ) = (s k,p ) of f p as m increases (where, for notational simplicity, the index p is xed and has been omitted).

For n = 2, the sequence of counting probability measures of the poles of R * m will asymptotically charge s 1 and s 2 (the poles will accumulate near s 1 and s 2 ), while only nitely many poles stay away from s 1 , s 2 and from the arc of circle orthogonal to T joining them [START_REF] Baratchart | Asymptotic uniqueness of best rational approximants of given degree to Markov functions in L 2 of the circle[END_REF]. These results are related to the fact that, for n = 2, f p can be represented as:

f p (z) = R 6 (z) + R 2 (z) s2 s1 dt (z -t) (t -s 1 )(s 2 -t) ,
for rationals R 2 ∈ R 2 and R 6 ∈ R 6 with, respectively, two simple poles for R 2 and two triple poles for R 6 , at s 1 and s 2 . This result is used in [START_REF] Baratchart | Meromorphic approximants to complex Cauchy transforms with polar singularities[END_REF] to study the behaviour of the poles of best rational approximants.

C.6 Behaviour of Triple Poles of Rational Approximants

We have the robustness property:

(P') For a function f p , close to a rational of R 3 with a single triple pole t ∈ D, then the triple pole t * of its best rational approximant R * 3 (with a single triple pole in D) and t are close to each other.

Property (P') is a restatement of: Proposition 1 Let R 3 (z) = π(z)/(z -t) 3 be a rational function of R 3 with a triple pole t ∈ D, strictly proper (deg π < 3), and irreducible (π(s) = 0). Then, there exists K > 0 such that for all rational function R3 (z) = π(z)/(z -b) 3 with a triple pole b ∈ D, the following inequality holds: |t -b| ≤ K R 3 -R3 , for the L 2 (T ) norm.

Indeed, as a corollary, if f p -R 3 ≤ (in L 2 (T ) norm), for some > 0 and R 3 ∈ R 3 with a single triple pole t in D, then we have the inequality: |t -t * | < K .

Proof of Proposition 1: Let R = π/q, and R = π/q, be two proper rational (irreducible) functions in L 2 (T ): in particular, π, q, π, q are polynomials, and the roots of q and q lie inside the open unit disk D (q and q are called stable polynomials).

Step 1 We look for a lower bound to d(R, R) = π q -π q (L 2 (T ) norm). If n is the degree of q, put q(z) = z n q(1/z) for its reciprocal polynomial. Then q/ q has modulus 1 on the unit circle, and we also have:

d(R, R) = π q q q -π q .

From the orthogonal decomposition of L 2 (T ) into Hardy spaces of analytic functions [START_REF] Rudin | Real and Complex Analysis[END_REF], and because P -(π/ q) = 0, since the poles of π/ q belong to C \ D.

we get that:

d(R, R) ≥ P - π q q q .

(14

)
Remark 2 The right hand-side of ( 14) vanishes if and only if q is a divisor of q.

Step 2 In the particular case where R ∈ R 3 and R have a single triple pole respectively at t ∈ D, and and at b ∈ D, then q(z) = (z -t) 3 , q(z) = (zb) 3 , q(z) = (1 -bz) 3 , and we can evaluate this right hand-side by a fractional decomposition. Indeed, expand π(z)q(z) in powers of (z -t) in order to obtain: P - π q q q = P - 1 (1 -bz) 3 A 1 z -t + A 2 (z -t) 2 + A 3 (z -t) 3 Expanding now 1/ q in a neighbourhood of z = t in D, we get:

1 q(z) = 1 (1 -bz) 3 = 1 (1 -bs) 3 1 + 3 b(z -t) 1 -bt + 6 b2 (z -t) 2
(1 -bt) 2 + O((z-t) 3 ) .

The P -projections in (15) can then be expressed as: (z -t) 2 + αA 1 + βA 2 + γA 3 z -t .

P - 1 (1 -bz) 3 1 z -t = α z-s , P - 1 (1 -bz) 3 1 (z -t) 2 = α (z-t) 2 + β z-t , P - 1 (1 -bz) 3 1 (z -t) 3 = α (z-t) 3 + β (z-t) 2 + γ z-t ,
Using ( 14), we thus get:

d(R, R) ≥ Φ t,b .
Note that Φ t,b , Ψ t,b , and N t,b have continuous coecients in the variable b ∈ D.

In particular, for b = t ∈ D:

Φ t,t (z) = 0 , Ψ t,t (z) = 3π(t) (1 -|t| 2 ) 3 1 z -t .

Step Then, ∀z ∈ D, qn (z)/ qn (z) = (z -b n ) 3 /(1-bn z) 3 converges to the constant -b 3 , whence (π qn )/(q qn ) → -b 3 π/q in L 2 (T ). Since P -(π/q) = π/q = R (because R m is contained in the above mentioned Hardy class of functions analytic in C \ D and by denition), we deduce that φ continuously extends to T and that φ(b) = R , ∀b ∈ T . Hence, ψ also admits a continuous extension on T and ψ(b) = R /|t -b|, ∀b ∈ T .

Step 4 Finally, the continuous positive function ψ attains its minimal value K ≥ 0 in D at some point b 0 ∈ D. In order to establish by contradiction that K > 0, assume that K = 0; then ψ(b 0 ) = 0 which implies φ(b 0 ) = 0. Then b 0 / ∈ T , since φ | T = R . But if b 0 ∈ D, then necessarily b 0 = t, the unique point in D where φ vanishes, from Remark 2. However, because 1 z-t = (1 -|t| 2 ) -1 2 , we have ψ(t) = 3|π(t)|(1 -|t| 2 ) -7 2 . Because R = π/q is irreducible, then π(t) = 0, so ψ(t) = 0, and this is a contradiction. Hence, K > 0.

∃K 0 > 0, t -b 1 -bt ≤ K 0 φ(b) ≤ K 0 d(R, R) . RR n°7704

Figure 1 :

 1 Figure 1: A nested conductor model for the head (left), in which brain, skull and scalp have a homogeneous conductivity. Our 3D spherical nested model (right, 2D view), where the brain and scalp (Ω 0 , Ω 2 ) conductivity is 0.33, the skull conductivity is 0.0042 and the air Ω 3 is non-conductive.

Figure 2 :

 2 Figure 2: Flowchart of the Source Localization method, consisting of two main steps: Data Transmission and Source Recovery. Source Recovery makes use of a 2D Rational Approximation technique, which is run independently on slices of the domain. The Planar Singularities are then analyzed jointly to yield the positions and moments of the sources.

3 .Figure 3 :

 33 Figure 3: Data Transmission. Cortical Mapping: from electrode data (a), yielding the normal current (b) and the potential (c) on the cortical surface. Harmonic Projection: (d) representation on the cortex S 0 of the singular part u a of the potential, whose singularities are restricted to Ω 0 .

Figure 4 :

 4 Figure 4: Plane Sections. Ω 0 is sliced into disks D p by a series of parallel planes Π p .

Figure 5 :

 5 Figure 5: True sources C 1 , C 2 and theoretical singularities s 1,p , s 2,p , p = 1, • • • , P , which when joining the dots form singular lines l k (see Section 3.4: 3D view (left), top view (right).

3 .

 3 If n = 1 then s k,p = s * k,p . 4. Otherwise, for f p linked to s k,p by (8) and starting with the initial values s * k,p , nd s k,p by minimizing (gradient descent) the criterion:

Figure 6 :

 6 Figure 6: Top views of 2D planar singularities computed from cortical dataset for 12 dierent slicing directions Π.

Figure 8 :

 8 Figure 8: True vs estimated sources: cortical data (left), electrodes data (right).

Figure 9 :

 9 Figure 9: Cortical mapping reconstruction. From the potential measured on 128 electrodes (top left), the Cortical Mapping method reconstructs the normal current and the potential on all surfaces of the model. The normal current is not represented on the scalp because it is simply equal to zero in our model. Note that the spatial distribution of the potential is less sharp on the scalp and skull surfaces (top right and middle right) than on the cortex, due to the high resistivity of the skull.

  where the singularities s k = s k,p ∈ D and σ k = σ k,p ∈ D are linked between each other and with the source parameters (C k being determined by the quantities z k , h k,p and r p ) by the relations:

  Also, for each slice index p = 1, • • • , P , |s k,p | |σ k,p | = 1 and arg s k,p = arg σ k,p = arg z k .

  In sections p close to C k , f p is (numerically) close to a rational function with poles at (s k,p ). Thus, for m large enough, typically m ≥ 3n, the m poles of R * m must then be located near the (s k,p ) (see Property (P), C.5), and they should be packed in a number of clusters coinciding with the number n of sources. This is illustrated in Figures 1012 where m = 1, 2 and 3 simple poles cases are respectively shown for the same situation with n = 2 sources as in Section 3.1.2. For k = 1, 2 and varying p, theoretically known singularities (r p s k,p ) are shown in disks D p as green dots whereas estimated poles are shown with red dots. Black dots represent the sources C 1 , C 2 .

Figure

  Figure 10: 1 simple pole Figure 11: 2 simple poles

  Indeed, in the slice p containing a source C k , the single triple pole s * p of R * 3 is close to the associated singularity s k,p . This situation is illustrated in Figure 13 for the same source conguration as the one used for Figures 1012. As previously, sources are indicated with black dots, green dots show the known theoretical positions of singularities while red dots show the triple poles. One can notice that there is only one pole trajectory for both singularity lines. Thus the triple pole line does not follow the singularity lines as in the single pole case. However, close to sources, the triple poles approximate the singularity lines quite well.

Figure 5

 5 Figure 13: 1 triple pole

  s)(t -b) + 3π (t)(t -b) 2 + π (t)(t -b) 3 /2 , A 2 = 3π(s)(t -b) 2 + π (t)(t -b) 3 , A 3 = π(t)(t -b) 3 .

  bt) 3 , β = 3 b (1 -bt) 4 , γ = 6 b2 (1 -bt) 5 . ,b (z) = (t -b) N t,b (z) (z -t) 3 = (t -b)Ψ t,b (z) ,where Φ t,b , Ψ t,b are rational functions, and N t,b is a polynomial of degree 2, as follows from:Φ t,b (z) = αA 3 (z -t) 3 + αA 2 + βA 3

3

 3 Consider the L 2 (T ) norms of Φ t,b and Ψ t,b as functions of b, and put φ(b) = Φ t,b , ψ(b) = Ψ t,b , whence φ(b) = |t -b|ψ(b) .

From ( 14

 14 ), we now have thatd(R, R) ≥ φ(b) .

  The above expressions are then to the eect that φ and ψ are continuous functions of b ∈ D. Further, φ and ψ admit continuous extensions up to the closed disk D, and ψ does not vanish on D. Indeed, let |b| = 1 and b n ∈ D, b n → b.

Step 5

 5 Thus ψ(b) ≥ K > 0 whence, with K = 1/K , φ(b) ≥ |t -b| and: |t -b| ≤ Kφ(b) ≤ Kd(R, R) ,which achieves the proof of Proposition 1, with R 3 = R.Remark 3 The expressions of α, β, γ lead to a similar result for the hyperbolic distance:

Table 1 :

 1 True vs estimated sources from cortical and electrodes datasets.

			positions			moments
	True sources	0.2000	0.3000	0.4000	0.0000	0.2000	0.6000
		-0.3000	-0.2000	0.4000	0.1000	0.0000	0.8000
	Estimated sources	0.1951	0.3056	0.4260	0.0194	0.2068	0.5874
	(cortical data)	-0.3006	-0.2059	0.4208	0.0808	-0.0052	0.7993
	Estimated sources	0.1917	0.2797	0.4160	0.0112	0.2579	0.5563
	(electrodes data)	-0.2798	-0.1777	0.4015	0.0929	-0.0486	0.8487
	Error	0.0049	-0.0056	-0.0260	-0.0194	-0.0068	0.0126
	(cortical data)	0.0006	0.0059	-0.0208	0.0192	0.0052	0.0007
	Error	0.0083	0.0203	-0.0160	-0.0112	-0.0579	0.0437
	(electrodes data)	-0.0202	-0.0223	-0.0015	0.0071	0.0486	-0.0487

In the remainder of this article, all domains are supposed to be open.RR n°7704

This is generically true with respect to the plane Π. See however B for the degenerated situation where z k = 0.RR n°7704

Recall that a pole (or polar singularity) is the zero of some polynomial at the denominator of the function, while a branchpoint is the singular point of some multivalued complex analytic function, as log or square-root.RR n°7704

RR n°7704

M. Zghal.

Problèmes inverses pour l'équation de Laplace en dimension 3, application à l'éléctroencéphalographie. PhD thesis, ENIT, Univ. Tunis El Manar, 2010.

RR n°7704