Signal processing methods

Abstract

We describe a method of processing an EEG and/or MEG signal to generate image data representing a 3D current distribution, J, within the brain, the method comprising: capturing a plurality of electric and/or magnetic measurements from the exterior of the head; solving an integral equation for a part of said current distribution to generate said image data representing said 3D current distribution, wherein said integral equation comprises an integral of a first function representing said part of said current distribution and of a second function (.sub.Tv.sub.s(r,)) representing the geometry and conductivity of the head independent of said current distribution; wherein said solving comprises: modelling the head as at least two regions separated by at least one internal boundary, and solving a set of partial differential equations, one for each said internal region, each partial differential equation comprising a geometry-conductivity function (w(r,)) representing the geometry and conductivity of the respective region, wherein said solving is subject to a boundary condition that either i) the gradients of the functions across the or each said internal boundary are smooth when conductivity is taken into account, or ii) a normal component of the electric field of said part of said current distribution is continuous across the or each said internal boundary, and wherein said geometry-conductivity function for an outermost said region of said head defines said second function (.sub.Tv.sub.s(r,)).

Claims

1. A method of processing an EEG and/or MEG signal to generate image data representing a 3D current distribution, J, within the brain, the method comprising: capturing a plurality of electric and/or magnetic measurements from the exterior of the head; solving an integral equation for a part of said current distribution to generate said image data representing said 3D current distribution, wherein said integral equation comprises an integral of a first function representing said part of said current distribution and of a second function (.sub.v.sub.s(r,)) representing the geometry and conductivity of the head independent of said current distribution; wherein said solving comprises: modelling the head as at least two regions separated by at least one internal boundary, and solving a set of partial differential equations, one for each said internal region, each partial differential equation comprising a geometry-conductivity function (w(r,)) representing the geometry and conductivity of the respective region, wherein said solving is subject to a boundary condition that either i) the gradients of the functions across the or each said internal boundary are smooth when conductivity is taken into account, or ii) a normal component of the electric field of said part of said current distribution is continuous across the or each said internal boundary, and wherein said geometry-conductivity function for an outermost said region of said head defines said second function (.sub.v.sub.s (r,)).

2. A method as claimed in claim 1 wherein said set of partial differential equations comprise a set of partial differential equations dependent said geometry and conductivity of the respective regions and independent of said 3D current distribution, each representing the effect of geometry and conductivity on an electric field in the region from the next inner region, and wherein the electric field in an innermost said region comprises a field generated by a monopole source.

3. A method as claimed in claim 2 wherein said innermost region contains said 3D current distribution and wherein a partial differential equation for said innermost said region is non-Laplacian.

4. A method as claimed in claim 1, wherein said set of partial differential equations comprises a partial differential equation for an innermost said region in the form:
.sub.c.sup.2w.sub.c.Math.Q(r) where r is a position vector within said innermost region, is a dummy variable, Q is a charge, is the delta function, w.sub.c is the geometry-conductivity function for said innermost region and .sub.c is an electrical conductivity of said innermost region; and wherein one or more partial differential equations for subsequent successively more outer regions j have the form:
.sup.2w.sub.j=0.

5. A method as claimed in claim 4 wherein said boundary conditions define that ${}_{j-1}\frac{w_{j-1}}{n}={}_j\frac{w_j}{n}$ where j1 denotes a region within region j, and where n denotes a unit vector outward normal to the surface of the boundary between regions j1 and j.

6. A method as claimed in claim 1 wherein said geometry-conductivity function for said outermost region of said head, w.sub.s(r,), defines said second function (.sub.v.sub.s(r,)) through the relationship that w.sub.s(r,) is proportional to said second function (w.sub.s(r,).sub.v.sub.s(r,)).

7. A method as claimed in claim 1 wherein said integral equation comprises one or both of: $u_s\left(r\right)=-\frac{1}{4}{}_{{}_c}\left({}^2\left(\right)\right)v_s\left(r,\right)\mathrm{dV}\left(\right),r{}_s,$ or an equation derivable therefrom, where u.sub.s(r) defines measured electric potential on a boundary of said outermost region of the head [the scalp] at position r, .sub.c defines an innermost said region, .sub.s defines a surface of said outermost region of the head, v.sub.s(r,) is an integral of said second geometry-conductivity function .sub.v.sub.s(r,) for said outermost region of said head, and is a scalar value representing an irrotational part of said 3D current distribution; and $\frac{4}{}B\left(r\right).Math.r=-{}_{{}_c}{}^2\left(.Math.A\left(\right)\right)\frac{\mathrm{dV}\left(\right)}{.Math.r-.Math.}+\frac{1}{4}{}_{{}_c}\left({}^2\left(\right)\right)r.Math.H\left(r,\right)\mathrm{dV}\left(\right)r{}_e,$ or an equation derivable therefrom, where B(r) defines a measured magnetic field at a sensor position r, .sub.c defines said innermost region of the head, and .sub.e defines an external region beyond said outermost said region of said head, H(r,) is a vectorial geometry-conductivity function for said head defined by said geometry-conductivity functions for said regions of said head (w(r,)), is a scalar value representing an irrotational part of said 3D current distribution and .Math.A() defines a radial component of a solenoidal part of said 3D current distribution.

8. A method as claimed in claim 1 further comprising determining a scalar value, , representing an irrotational part of said 3D current distribution, J, from said plurality of electric and/or magnetic measurements.

9. A method as claimed in claim 8 wherein said determining of said scalar value, , uses an expansion of in terms of expansion basis functions having respective unknown coefficients defined by a vector , where b=E.Math., where b comprises a vector representation of said measurements of electric potential, b(r.sub.i) at positions r.sub.i, and where elements E(i,j) of matrix E are defined by a combination of a jth respective basis function and said geometry-conductivity function for said outermost region of said head, w.sub.s(r,) at position r.sub.i, wherein E represents the geometry and conductivity of the head independent of said current distribution, the method comprising determining an estimate of said unknown coefficients, , for a solution of b=E.Math., to determine said representation of .

10. A method as claimed in claim 8 wherein said EEG/MEG signal is an EEG signal, wherein said plurality of electric and/or magnetic measurements comprises a plurality of measurements of electric potential on the scalp (b), and wherein said image data includes a representation of .

11. A method as claimed in claim 1 wherein said EEG/MEG signal comprises an MEG signal, the method comprising determining a radial component of a solenoidal part of said 3D current distribution, and wherein said image data comprises a representation of said radial component of a solenoidal part of said 3D current distribution.

12. A method as claimed in claim 11 wherein said determining of said radial component of said solenoidal part of said 3D current distribution, .Math.A() uses an expansion of .Math.A() in terms of expansion basis functions having respective unknown coefficients defined by a vector , where M=c, where c comprises a vector representation of said measurements of magnetic field, B(r.sub.i) at sensor positions r.sub.i, and where elements M(i,j) of matrix M are defined by a combination of a jth respective basis function and said position r.sub.i, wherein M represents the geometry of the head independent of said current distribution.

13. A method as claimed in claim 12 wherein c further comprises a term comprising representing an irrotational part of said 3D current distribution in combination with a vectorial geometry-conductivity function for said head, H(r,) dependent on said geometry-conductivity functions for said regions of said head (w(r,)).

14. A method as claimed in claim 1 wherein said modelling comprises modelling the head as three or more regions, an innermost region of the head defining a cerebellum region, the next regions defining, in order, two or more of: a cerebrospinal fluid shell region, a skull bone shell region, and a scalp shell region.

15. Apparatus for processing an EEG and/or MEG signal to generate image data representing a 3D current distribution, J, within the brain, the apparatus comprising: an input to receive a plurality of captured electric and/or magnetic measurements from the exterior of the head; and a processor coupled to said input, to an output, to working memory and to program memory storing processor control code, wherein the stored code comprises code to: solve an integral equation for a part of said current distribution to generate said image data representing said 3D current distribution, wherein said integral equation comprises an integral of a first function representing said part of said current distribution and of a second function (.sub.v.sub.s(r,)) representing the geometry and conductivity of the head independent of said current distribution; wherein said solving comprises: modelling the head as at least two regions separated by at least one internal boundary, and solving a set of partial differential equations, one for each said internal region, each partial differential equation comprising a geometry-conductivity function (w(r,)) representing the geometry and conductivity of the respective region, wherein said solving is subject to a boundary condition that either i) the gradients of the functions across the or each said internal boundary are smooth when conductivity is taken into account, or ii) a normal component of the electric field of said part of said current distribution is continuous across the or each said internal boundary, and wherein said geometry-conductivity function for an outermost said region of said head defines said second function (.sub.v.sub.s(r,)); and output said image data, wherein said image data comprises data representing a part of said 3D current distribution.

16. Apparatus for processing an EEG signal to generate image data representing a 3D current distribution, J, within the brain, the apparatus comprising: an input to receive a plurality of measurements of electric potential on the scalp (b); and a processor coupled to said input, to an output, to working memory and to program memory storing processor control code, wherein the stored code comprises code to: use a representation of a part of said current distribution as an expansion having a series of unknown coefficients, , of expansion basis functions, where b is dependent on a product of a matrix E and , and where E represents the geometry and conductivity of the head independent of said current distribution, wherein said geometry includes at least one internal boundary within the head, to solve an equation combining b and E.Math. to determine an estimate of , wherein defines said 3D current distribution; wherein E.Math. represents
.sub..sub.c(.sup.2())v.sub.s(r,)dV() where is a scalar value representing J, where v.sub.s(r,) is a function representing the geometry and conductivity of the head independent of said current distribution, .sub.c denotes the region of the brain over which the current distribution is determined, r lies on the surface of the scalp on which EEG measurements are made, and is a dummy variable; and wherein, in said matrix E, a value of an element of the matrix E is defined by an integral, over a region including said current distribution, of a gradient of a said expansion basis function and a function w.sub.s(r,) representing v.sub.s(r,); where w.sub.s(r,) alternatively represents the geometry and conductivity of the head independent of said current distribution, wherein a function w(r,) is determined for each internal region of the head defined by said at least one internal boundary, wherein w.sub.s(r,) is the outermost said function w(r,), and wherein w.sub.s(r,) is determined subject to the constraint that, when conductivity is taken into account, the gradients of the functions w(r,) on either side of different regions of the head separated by a said internal boundary are smooth across each said internal boundary such that a normal component of electric field is continuous across each said internal boundary; and output image data representing said 3D current distribution determined from said estimate of .

17. Apparatus as claimed in claim 16 wherein said functions w(r,) for each internal region of the head are determined by solving a set of partial differential equations for a set of functions w.sub.j(r,), one for each said region where j labels a region, wherein the partial differential equation for each said region except the innermost comprises Laplace's equation for the respective function w.sub.j(r,), and wherein the partial differential equation for the innermost said region comprises Poisson's equation for the innermost region function w.sub.innermost(r,); and wherein said solving is subject to a boundary condition that across each said internal boundary a normal component of electric field is continuous.

18. Apparatus as claimed in claim 16 wherein said equation combining b and E.Math. is b=E.Math. and wherein said expansion basis functions comprise radial basis functions.

19. A medical imaging device comprising the apparatus of claim 16 and a display to display said image data representing said 3D current distribution within the brain.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) These and other aspects of the invention will now be further described, by way of example only, with reference to the accompanying figures in which:

(2) FIGS. 1a to 1d show, respectively, a flow diagram of a computer procedure for solving an inverse EEG problem according to an embodiment of the invention, a flow diagram of a computer procedure for solving an inverse MEG problem according to an embodiment of the invention, a flow diagram of a computer procedure for solving a set of partial differential equations to determine a matrix E and/or vectorial function H for use in numerically solving an inverse EEG/MEG problem according to embodiments of the invention, and a computer system programmed to implement the methods of FIGS. 1a to 1c;

(3) FIGS. 2a and 2b show a comparison of analytical solution and numerical solutions for u.sub.s (r,) and for v.sub.s (r,) respectively;

(4) FIG. 3 illustrates solutions of a boundary value problem for a multi-layer spherical geometry;

(5) FIG. 4 illustrates elements of a matrix .sub.H(r,) problem for a multi-layer spherical geometry;

(6) FIGS. 5a to 5c show, respectively, the cerebrum, .sub.c, the skull, .sub.b, and the head, .sub.s of a three layer realistic head model;

(7) FIG. 6 shows values of a synthetic function () inside the cerebrum of the three layer head model of FIG. 5;

(8) FIG. 7 shows exact and reconstructed functions () for the three layer head model of FIG. 5;

(9) FIG. 8 shows values of a synthetic function A.sup. () inside the cerebrum of the three layer head model of FIG. 5; and

(10) FIG. 9 shows exact and reconstructed functions A.sup. () for the three layer head model of FIG. 5.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

1. Framework

(11) In the following description we denote a domain and its conductivity by R.sup.3 and respectively. Let the bounded domain .sub.c represent the cerebrum, which has conductivity .sub.c. A shell .sub. with conductivity .sub., representing the cerebrospinal fluid (CSF), surrounds the domain .sub.c. The CSF is surrounded by the skull characterized by the domain .sub.b with conductivity .sub.b. Finally, the skull is surrounded by the scalp, which is modelled as a shell .sub.s with conductivity .sub.s. The domain exterior to the head is denoted by .sub.e which is not conductive. The permeability of all domains is equal to the permeability .sub.0 of empty space. Let J.sup.P(), .sub.c, denote the current which is supported within the cerebrum .sub.c. For the above situation of arbitrary geometry and arbitrary current, it is shown in A S Fokas (J. R. Soc. Interface, 6:479-488, 2009, ibid) that the irrotational part of the current, which is denoted by the scalar function (), contributes to both the electric potential on the scalp and to the magnetic field in .sub.e, On the other hand, the solenoidal part of the current, which is denoted by the vectorial function A(), does not affect the electric potential and furthermore only the radial component of A() affects the magnetic field in .sub.e.

(12) Being more precise, an arbitrary vectorial function with support .sub.c, can be expressed as
J.sup.P()=()+A(), .sub.c,(1)
where A() satisfies the constraint .Math.A()=0. This constraint implies that J.sup.P() involves three arbitrary scalar functions, namely the scalar function () and the two independent scalar functions characterizing A().

(13) It is shown in A S Fokas (J. R. Soc. Interface, 6:479-488, 2009, ibid) that the electric potential on the scalp is given by the expression

(14) $\begin{array}{cc}{u}_s\left(r\right)=-\frac{1}{4}{}_{{}_c}\left({}^2\left(\right)\right)v_s\left(r,\right)\mathrm{dV}\left(\right),r{}_s,& \left(2\right)\end{array}$
where .sub.s denotes the surface of the scalp, i.e the smooth boundary of the domain .sub.s. Here, v.sub.s(r,) denotes a harmonic function which is independent of the current; it depends only on the geometry and on the conductivities. It is clear from equation (2) that the electric potential u.sub.s(r) depends only on the Laplacian of the irrotational part (()) of the current. It was also shown in A S Fokas (J. R. Soc. Interface, 6:479-488, 2009, ibid) that the measured magnetic field in .sub.e is given by the expression

(15) 0$\begin{array}{cc}\frac{4}{}B\left(r\right).Math.r=-{}_{{}_c}{}^2\left(.Math.A\left(\right)\right)\frac{\mathrm{dV}\left(\right)}{.Math.r-.Math.}+\frac{1}{4}{}_{{}_c}\left({}^2\left(\right)\right)r.Math.H\left(r,\right)\mathrm{dV}\left(\right)r{}_e,& \left(3\right)\end{array}$
where H(r,) depends only on the geometry and on the conductivities; it is independent of the current. It is clear from equation (3) that in MEG the radial component of the magnetic field is affected by () as well as by the radial component .Math.A() of the solenoidal part of the current.

(16) We can use Green's first identity to rewrite equations (2) and (3) in the form

(17) $\begin{array}{cc}{u}_s\left(r\right)=\frac{1}{4}{}_{{}_c}\left(\right).Math.{}_{}{v}_s\left(r,\right)\mathrm{dV}\left(\right),r{}_s\mathrm{and}& \left(4\right)\\ \frac{4}{}B\left(r\right).Math.r=-{}_{{}_c}{}^2\left(.Math.A\left(\right)\right)\frac{\mathrm{dV}\left(\right)}{.Math.r-.Math.}-\frac{1}{4}\underset{l\left\{1,2,3\right\}}{.Math.}{}_{{}_c}{r}_l\left[\left(\right).Math.{}_{}{H}_l\left(r,\right)\right]\mathrm{dV}\left(\right),r{}_e.& \left(5\right)\end{array}$

(18) Equations (4) and (5) follow from equations (2) and (3) by employing the following identities
.sub..sub.c[H.sub.i(r,).sub..sup.2()+.sub.H.sub.i(r,).Math..sub.()]dV()=[.sub.().Math.{circumflex over (n)}]H.sub.i(r,)dS()i{1,2,3},
H.sub.i(r,)[.sub.().Math.{circumflex over (n)}]dS()=0, i{1,2,3},
v.sub.s(r,)[.sub.().Math.{circumflex over (n)},]dS()=0,(6)

(19) where {H.sub.1, H.sub.2, H.sub.3} denote the Cartesian coordinates of H(r,).

(20) We here present a general framework for obtaining the maximum possible information for the current, or more precisely for the part of the current that affects EEG and MEG, namely
() and .sup.2[.Math.A()], .sub.c,(7)
from the knowledge of the quantities measured in EEG and in MEG, namely the functions
u.sub.s(r)r.sub.s, and r.Math.B(r) r.sub.e,(8)
respectively. This framework is based on (a) equations (4) and (5), and (b) on the existence of fast and accurate codes for the numerical computation of the auxiliary functions v.sub.s(r,) and H(r,). Furthermore, we implement this general approach to a realistic three shell geometry.

2. Construction of vs(r,), and H(r,)

(21) In order to analyse the basic equations (4) and (5) we first construct the functions
v.sub.s(r,)r.sub.s, .sub.c, and H(r,)r.sub.e, .sub.c.(9)

(22) The function v.sub.s(r,) is a harmonic function which is obtained via the solution of the following system of equations:

(23) $\begin{array}{cc}{}^2{v}_c\left(r,\right)=0,r{}_c,{}_c,\frac{}{n}\left[\frac{1}{.Math.r-.Math.}+v_c\left(r,\right)\right]={}_f\frac{v_f\left(r,\right)}{n},r{}_c;& \left(10\right)\\ {}^2{v}_f\left(r,\right)=0,r{}_f,{}_c,v_f\left(r,\right)=\frac{1}{{}_c}\left[\frac{1}{.Math.r-.Math.}+v_c\left(r,\right)\right],r{}_c,{}_f\frac{v_f\left(r,\right)}{n}={}_b\frac{v_b\left(r,\right)}{n},r{}_f,{}_c;& \left(11\right)\\ {}^2{v}_b\left(r,\right)=0,r{}_b,{}_c,v_b\left(r,\right)=v_f\left(r,\right),r{}_f,{}_b\frac{v_b\left(r,\right)}{n}={}_s\frac{v_s\left(r,\right)}{n},r{}_b,{}_c;& \left(12\right)\\ {}^2{v}_s\left(r,\right)=0,r{}_s,{}_c,v_s\left(r,\right)=v_b\left(r,\right),r{}_b,\frac{v_s\left(r,\right)}{n}=0,r{}_s.& \left(13\right)\end{array}$

(24) In the above equations, .sub.c, .sub.b, .sub. and .sub.s denote the smooth boundaries of the domains .sub.c, .sub.b, .sub. and .sub.s, respectively.

(25) As stated earlier, equations (10)-(13) are independent of the current J.sup.P() and depend only on the geometry and on the conductivities .sub.c, .sub.b, .sub. and .sub.s. However, it turns out that it is convenient to construct the functions v(r,) via first constructing certain functions u.sub.j(r,) which are defined in terms of a monopole source. Here the functions u.sub.j (r,) represent the geometry-conductivity functions denoted w(r,) in the summary of the invention section and the claims (here j{c, , b, s} denoting the cerebrum, cerebrospinal fluid, bone and scalp respectively). The function u.sub.s (r,) should not be confused with the EEG measurements as a function of position, u.sub.s (r.sub.i).

(26) The boundary condition of equation (10) makes the set of partial differential equations extremely difficult to solve. However in the functions below this boundary condition has been simplified by changing, and although this means that the partial differential equation for the cerebrum no longer has a Laplacian form, overall the problem becomes tractable: Let the functions
u.sub.j(r,), r.sub.j, .sub.c, j{c,,b,s},(14)
satisfy the following equations:

(27) $\begin{array}{cc}{}_c{}^2{u}_c\left(r\right)=.Math.Q\left(r-\right),r{}_c,{}_c\frac{u_c}{n}={}_f\frac{u_f\left(r\right)}{n},r{}_c;& \left(15\right)\\ {}^2{u}_f\left(r\right)=0,r{}_f,u_f\left(r\right)=u_c\left(r\right),r{}_c,{}_f\frac{u_f\left(r\right)}{n}={}_b\frac{u_b}{n},r{}_f;& \left(16\right)\\ {}^2{u}_b\left(r\right)=0,r{}_b,u_b\left(r\right)=u_f\left(r\right),r{}_f,{}_b\frac{u_b\left(r\right)}{n}={}_s\frac{u_b\left(r\right)}{n},r{}_b;& \left(17\right)\\ {}^2{u}_s\left(r\right)=0,r{}_s,u_s\left(r\right)=u_b\left(r\right),r{}_b,\frac{u_s\left(r\right)}{n}=0,r{}_s.& \left(18\right)\end{array}$

(28) Then u.sub.j and v.sub.j are related by the equation

(29) $\begin{array}{cc}{u}_j\left(r,\right)=\frac{1}{4}Q\left(\right).Math.{}_{}{v}_j\left(r,\right),j\left\{f,b,s\right\},r{}_j,{}_c.& \left(19\right)\end{array}$

(30) Equation (19) above is consistent with the fact that the potential of a dipole is the directional derivative of the potential of a monopole.

(31) 2.1 Spherical Geometry

(32) The advantage of constructing the function v.sub.s(r,) via u.sub.s(r,) becomes clear by considering the particular case of spherical geometry. In this case, a closed form expression for v.sub.s(r,) can be derived, which involves the inversion of a 77 matrix. It was observed in Fokas et al., Inverse Problems, 28, 2012 (ibid) that this matrix is ill conditioned and requires regularization. On the hand, for the case of N spherical layers, with N arbitrary, a closed form expression for u.sub.s(r,) can be obtained. This approach is straightforward and does not require matrix inversion. Indeed, let r.sub.1<r.sub.2< . . . <r.sub.N-1<r.sub.N, and let .sub.1, . . . ,.sub.N denote the radii and the corresponding conductivities of the domains .sub.1, . . . , .sub.N. We consider a single dipole source characterized by (Q,). The position vector of the observation point is denoted by r. Let us introduce the following notations:

(33) $\begin{array}{cc}\hat{}=\frac{}{.Math..Math.},\hat{r}=\frac{r}{.Math.r.Math.},\cos \left(\right)=\hat{}.Math.\hat{Q},\cos \left(\right)=\hat{}.Math.\hat{r},\cos \left(\right)=\left(\hat{}\hat{Q}\right).Math.\left(\hat{}\hat{r}\right).& \left(20\right)\end{array}$

(34) The observed potential at r is given by the following expression:

(35) $\begin{array}{cc}{u}_s\left(r,\right)=\frac{.Math.Q.Math.}{4{}_N{r}^2}\underset{n=1}{\overset{}{.Math.}}{\left(\frac{}{r}\right)}^{n-1}\left[\frac{2n+1}{{\mathrm{nM}}_{2,2}+\left(1+n\right)M_{2,1}}\right]\left[n\cos \left(\right)P_n\left[\cos \left(\right)\right]+\cos \left(\right)\sin \left(\right)P_n^1\left[\cos \left(\right)\right]\right],& \left(21\right)\end{array}$
where the 22 matrix M is given by

(36) $\begin{array}{cc}\left[\begin{array}{cc}{M}_{1,1}& M_{1,2}\\ M_{2,1}& M_{2,2}\end{array}\right]=\left[\frac{1}{{\left(2n+1\right)}^{N-1}}\right]\underset{k=1}{\overset{N-1}{.Math.}}\left[\begin{array}{cc}n+\left(n+1\right)\frac{{}_k}{{}_{k+10}}& \left(n+1\right)\left(\frac{{}_k}{{}_{k+1}}-1\right){\left(\frac{r}{r_k}\right)}^{2n+1}\\ n\left(\frac{{}_k}{{}_{k+1}}-1\right){\left(\frac{r_k}{r}\right)}^{2n+1}& n+1+\frac{n{}_k}{{}_{k+1}}\end{array}\right].& \left(22\right)\end{array}$

(37) Noting that v.sub.s(r,) can be expressed in the form below (A S Fokas, J. R. Soc. Interface, 6:479-488, 2009, ibid)

(38) $\begin{array}{cc}{v}_s\left(r,\right)=\underset{n=1}{\overset{}{.Math.}}{H}_n\frac{2n+1}{4}\frac{{}^n}{r^{n+1}}{P}_n\left(\hat{}.Math.\hat{r}\right),& \left(23\right)\end{array}$
where the unknown coefficients H.sub.n can be computed by inverting the 77 matrix mentioned earlier. By using equations (19) and (23) we bypass this inversion and can directly obtain v.sub.s(r,) in closed form:

(39) $\begin{array}{cc}{v}_s\left(r,\right)=\underset{n=1}{\overset{}{.Math.}}\left[\frac{2n+1}{{\mathrm{nM}}_{2,2}+\left(1+n\right)M_{2,1}}\right]\frac{1}{{}_N}\frac{{}^n}{r^{n+1}}{P}_n\left(\hat{}.Math.\hat{r}\right).& \left(24\right)\end{array}$
2.2 Arbitrary Geometry

(40) In the case of arbitrary geometry, there exists (open source) code suitable for the numerical solution of equations (15)-(18), for example Zeynep Akalm-Acar and Scott Makeig, Neuroelectromagnetic forward head modeling toolbox, J Neurosci Methods, 190(2):258-270, 2010. The numerical construction of v.sub.s(r,) involves the following steps: 1. Fix a source point :=[.sub.1, .sub.2, .sub.3].sub.c and an observation point r.sub.s. Consider three dipoles positioned at the source .sub.c, with the following orthogonal moments:
2. Q.sub.1=[1,0,0], Q.sub.2=[0,1,0], Q.sub.3=[0,0,1].(25) 3. For each of the above three dipoles, solve the boundary value problem described by equations (15)-(18). We denote the solution for the potential due to the dipole oriented in the direction Q.sub.i by u.sub.s(r,;Q.sub.i). We define the vector
4. u.sub.s(r,):[u.sub.s(r,;Q.sub.1)u.sub.s(r,;Q.sub.2)u.sub.s(r,;Q.sub.3)](26) 5. We obtain the gradient of v.sub.s(r,) from equation
.sub.v.sub.s(r,)=4u.sub.s(r,).(27) 6. Then, by using the gradient theorem, we can compute v.sub.s(r,) from the equation
v.sub.s(r,)v.sub.s(r,0)=.sub.[0,].sub.v.sub.s(r,l).Math.dl.(28) 7. We can approximate the solution of the above equation by the expression

(41) 0$\begin{array}{cc}{v}_s\left(r,\right)4\underset{i=1}{\overset{N}{.Math.}}{l}_{i}u\left(r,{}_i\right).Math.\hat{},& \left(29\right)\end{array}$ 8. where =(l.sub.1+ . . . +l.sub.N){circumflex over ()}. 9. We repeat the above steps .sub.c,r.sub.s. It should be noted that for a given dipole characterized by (Q,), most existing numerical solvers provide the solution for all discrete points r.sub.s in a single run.
2.3 Construction of H(r,)

(42) The vectorial function H(r,) is uniquely defined in terms of the scalar functions v.sub.j(r,):

(43) $\begin{array}{cc}H\left(r;\right)={}_{{}_c}\left[\frac{1}{.Math.r^{}-.Math.}+v_c\left(r^{},\right)-{}_f{v}_f\left(r^{},\right)\right]\left[\hat{n}\left(r^{}\right){}_{r^{}}\frac{1}{.Math.r-r^{}.Math.}\right]\mathrm{dS}\left(r^{}\right)+{}_{{}_f}\left({}_f{v}_f\left(r^{},\right)-{}_b{v}_b\left(r^{},\right)\right)\left[\hat{n}\left(r^{}\right){}_{r^{}}\frac{1}{.Math.r-r^{}.Math.}\right]\mathrm{dS}\left(r^{}\right)+{}_{{}_b}\left({}_b{v}_b\left(r^{},\right)-{}_s{v}_s\left(r^{},\right)\right)\left[\hat{n}\left(r^{}\right){}_{r^{}}\frac{1}{.Math.r-r^{}.Math.}\right]\mathrm{dS}\left(r^{}\right)+{}_s{}_{{}_s}{v}_s\left(r^{},\right)\left[\hat{n}\left(r^{}\right){}_{r^{}}\frac{1}{.Math.r-r^{}.Math.}\right]\mathrm{dS}\left(r^{}\right)r{}_e,{}_c.& \left(30\right)\end{array}$
2.3.1 Construction of .sub.H(r,) in Spherical Geometry

(44) In the particular case of the multi-layer spherical geometry, an analytical expression for the magnetic field due to a dipole characterized by (Q,) is as follows:

(45) $\begin{array}{cc}B\left(r,;Q\right)=\frac{{}_0}{4F^2\left(r,\right)}\left(F\left(r,\right)Q-\left(Q\right).Math.\left[r{}_{r}F\left(r,\right)\right]\right),& \left(31\right)\end{array}$

(46) where F(r,) and .sub.rF(r,) are given by

(47) $\begin{array}{cc}F\left(r,\right)=d\left(rd+r^2-\left(.Math.r\right)\right),{}_{r}F\left(r,\right)=\left(\frac{d^2}{r}+\frac{\left(r-\right).Math.r}{d}+2d+2r\right)r-\left(d+2r+\frac{\left(r-\right).Math.r}{d}\right),d:=.Math.r-.Math.,r=.Math.r.Math..& \left(32\right)\end{array}$

(48) For arbitrary geometry, an expression for the magnetic field due to a dipole characterized by (Q(),) was derived in A S Fokas J. R. Soc. Interface, 6:479-488, 2009, ibid:

(49) $\begin{array}{cc}\frac{4}{\mathrm{.Math.}}B\left(r,;Q\right)=Q\left(\right){}_r\frac{1}{.Math.r-.Math.}-\frac{1}{4}Q\left(\right).Math.{}_{}H\left(r,\right),& \left(33\right)\end{array}$

(50) where H(r,) is the vectorial function defined in equation (30). The gradient .sub.H(r,) is the following 33 matrix:

(51) $\begin{array}{cc}{}_{}H\left(r,\right)=\left(\begin{array}{ccc}\frac{H_1}{{}_1}& \frac{H_2}{{}_1}& \frac{H_3}{{}_1}\\ \frac{H_1}{{}_2}& \frac{H_2}{{}_2}& \frac{H_3}{{}_2}\\ \frac{H_1}{{}_3}& \frac{H_2}{{}_3}& \frac{H_3}{{}_3}\end{array}\right).& \left(34\right)\end{array}$

(52) Using equations (31) and (33), it is straightforward to show that

(53) $\begin{array}{cc}\frac{H_j}{{}_i}=4\left[Q_i{}_{}\frac{1}{.Math.r-.Math.}\right].Math.Q_j-\frac{{\left(4\right)}^2}{}B\left(r,;Q_i\right).Math.Q_j,i,j\left\{1,2,3\right\}.& \left(35\right)\end{array}$
2.3.2 Numerical Construction of H(r,)

(54) There exists several numerical solvers for computing the magnetic field problem due to a dipole or a collection of dipoles, for example Zeynep Akalm-Scott et al., ibid. In what follows, we outline a strategy to compute H(,r) using an existing solver. We denote by B(r,;Q) the solution of the magnetic field due to a dipole characterized by (,Q). 1. Fix a source point :=[.sub.1, .sub.2, .sub.3].sub.c and an observation point r.sub.e. Consider three dipoles positioned at the source .sub.c with the following orthogonal moments:
2. Q.sub.1=[1,0,0], Q.sub.2=[0,1,0], Q.sub.3=[0,0,1].(36) 3. For each of the above three dipoles, solve the magnetic field problem. We denote by B(r,;Q.sub.i) the solution of the magnetic field associated with the dipole oriented in the direction Q.sub.i. 4. We can now compute the elements of the matrix .sub.H(,r) using the expression (16). 5. Repeat the above steps for .sub.c, and for all sensors positions r.sub.i .sub.e, where i denotes the ith sensor. The existing numerical solvers (ibid) compute the magnetic field at all sensors positions in a single run. 6. The gradient theorem,
H.sub.j(,r)H.sub.j(0,r)=.sub.[0,].sub.H.sub.j(l,r).Math.dl j{1,2,3},(37) 7. can be employed to compute the components of the vectorial function H(r,).
2.4 Radial Basis Functions

(55) Consider the Beppo-Levi space of distributions with square integrable second derivatives. This space is equipped with a rotation invariant semi norm: If sBL.sup.(2)(R.sup.3) then this norm is defined by

(56) $\begin{array}{cc}{.Math.s.Math.}^2={}_{R^3}\left[{\left(\frac{{}^{2}s\left(x\right)}{x^2}\right)}^2+{\left(\frac{{}^{2}s\left(x\right)}{y^2}\right)}^2+{\left(\frac{{}^{2}s\left(x\right)}{z^2}\right)}^2\right]\mathrm{dxdydz}+{}_{R^3}\left[2{\left(\frac{{}^{2}s\left(x\right)}{xy}\right)}^2+2{\left(\frac{{}^{2}s\left(x\right)}{xz}\right)}^2+2{\left(\frac{{}^{2}s\left(x\right)}{yz}\right)}^2\right]\mathrm{dxdydz}.& \left(38\right)\end{array}$

(57) Consider the following interpolation problem: Given a set of distinct nodes X={x.sub.i}.sub.i=1.sup.N, and function values {.sub.i}.sub.i=1.sup.N, we construct an interpolation function s(x) such that
s(x.sub.i)=.sub.i, i=1, . . . ,N.(39)

(58) By employing radial basis functions (RBF), an interpolation function can be expressed in the form

(59) $\begin{array}{cc}s\left(x\right)=p\left(x\right)+\underset{i=1}{\overset{N}{.Math.}}{}_i\left(.Math.x-x_i.Math.\right),& \left(40\right)\end{array}$

(60) where p(x) is the polynomial
p(x)=c.sub.1+c.sub.2x+c.sub.3y+c.sub.4z.(41)

(61) In this setting, {.sub.i:1iN} and {c.sub.1, . . . , c.sub.4}, are coefficients to be determined, whereas (r): R.fwdarw.[0,). Particular choices for the functions (r) include (r)=r, (r)=e.sup.ar.sup.2, (r)=(r.sup.2+c.sup.2).sup.1/2, and (r)=r.sup.2 log(r). The smoothest interpolant is given by

(62) $\begin{array}{cc}{s}^{*}\left(x\right)=p\left(x\right)+\underset{i=1}{\overset{N}{.Math.}}{}_i.Math.x-x_i.Math..& \left(42\right)\end{array}$

(63) In certain problems, the following constraint is also imposed

(64) 0$\begin{array}{cc}\underset{i=1}{\overset{N}{.Math.}}{}_{i}p\left(x_i\right)=0.& \left(43\right)\end{array}$

(65) For EEG, =(), and for MEG, =A.sup. (). The unknown coefficients .sub.1, . . . .sub.N, c.sub.1, . . . , c.sub.4, are linearly related to the function (x), thus the interpolation problem reduces to the following least squares problem:

(66) $\begin{array}{cc}\left[\begin{array}{cc}A& P\\ P^T& 0\end{array}\right]\left[\begin{array}{c}\\ c\end{array}\right]=\left[\begin{array}{c}f\\ 0\end{array}\right],& \left(44\right)\end{array}$

(67) where the matrix AR.sup.NN, is defined by
A.sub.i,j=(x.sub.ix.sub.j)(45)

(68) and every row of the matrix PR.sup.N4 is defined as
P(i,.)=[1 x.sub.i y.sub.i z.sub.i] 1iN.(46)

3. The Inverse Problem for EEG

(69) We analyse equation (4). A radial basis function approximation of () can be written in the form

(70) $\begin{array}{cc}\left(\right)=\underset{j=1}{\overset{N}{.Math.}}{}_j(.Math.-{}_j.Math..Math.),& \left(47\right)\end{array}$

(71) where (.sub.j|)=(.sub.j.sup.2+c.sup.2).sup.1/2. Here, {.sub.j:1iN} are the coefficients to be estimated, and c>0 is a constant.

(72) By substituting equation (47) into equation (4) we can formulate the following least squares problem:
b=E,(48)
where
E(i,j)=.sub..sub.c[.sub.(.sub.j|).Math.u.sub.s(r.sub.i,)]dV(),
b(i)=u.sub.s(r.sub.i)(49)

(73) and i denotes the ith sensor index. We can compute the volume integral of equation (49) numerically. The parametrization of the problem acts as an implicit regularization. We assume that the measurements are corrupted by additive noise. In this setting, the maximum likelihood estimate (MLE) yields
{circumflex over ()}=(E.sup.TC.sub.w.sup.1E).sup.1E.sup.TC.sub.w.sup.1b,(50)

(74) where C.sub.w is the covariance matrix of the measurement noise or the sample covariance matrix. The volume integrals of equation (49) are computationally expensive. However, this computation can be parallelized. More precisely, for every matrix entry (i, j), the process of solving the volume integral is independent from every other matrix entry. Furthermore, for a given geometry, the matrix E in equation (49) needs to be computed only once.

(75) Optionally, since head geometries are similar, an approximate result may be obtained using a common head geometry for a group of patients, for example all adult patients.

4. Inverse Problem for MEG

(76) In embodiments in order to solve the MEG inverse problem we first solve the EEG inverse problem. We assume that the magnetic sensors are oriented radially (though the skilled person will appreciate that this is merely a conveniencefor example a radial component of the field may be derived from a sensor at any orientation). We can rewrite equation (5) in the form

(77) $\begin{array}{cc}{}_{{}_c}{}^2\left(.Math.A\left(\right)\right)\frac{\mathrm{dV}\left(\right)}{.Math.r-.Math.}=-\frac{1}{4}\underset{l\left\{1,2,3\right\}}{.Math.}{}_{{}_c}{r}_l\left[\left(\right).Math.{}_{}{H}_l\left(r,\right)\right]\mathrm{dV}\left(\right)-\frac{4}{\mathrm{.Math.}}B\left(r\right).Math.r.& \left(51\right)\end{array}$

(78) We assume that an estimate of the function {circumflex over ()}() has been computed from EEG data, thus all quantities on the right hand side of equation (51) are known. We denote the radial component of the vectorial function A() by A.sup. ():={circumflex over ()}.Math.A(). A radial basis function expansion of the unknown function A.sup. () takes the form

(79) $\begin{array}{cc}{A}^{}\left(\right)=\underset{j=1}{\overset{N}{.Math.}}{}_j(.Math.-{}_j.Math..Math.).& \left(52\right)\end{array}$

(80) In order to formulate a least squares problem, we introduce the following notations:

(81) $\begin{array}{cc}M\left(i,j\right)={}_{{}_c}{}_{}^2(.Math.-{}_j.Math..Math.)\frac{\mathrm{dV}\left(\right)}{.Math.r_i-.Math.};r_i{}_e,c\left(i\right)=-\frac{1}{4}\underset{l\left\{1,2,3\right\}}{.Math.}{}_{{}_c}{r}_{\left(i,l\right)}\left[\left(\right).Math.{}_{}{H}_l\left(r_i,\right)\right]\mathrm{dV}\left(\right)-\frac{4}{}B\left(r_i\right).Math.r_i.& \left(53\right)\end{array}$

(82) where i denotes the ith sensor index. Then equation (51) yields
M=c.(54)

(83) As in the case of EEG, the employed parametrization acts as an implicit regularization. The MLE estimates yields
{circumflex over ()}=(M.sup.TC.sub.w.sup.1M).sup.1M.sup.TC.sub.w.sup.1c,(55)

(84) where C.sub.w is the covariance matrix of the noise associated with the magnetic field measurement system. The volume integrals of equation (53) are computationally expensive. However, these computations can be parallelized. More precisely, for an arbitrary matrix entry (i, j), the process of solving the volume integral is independent from every other matrix entry. Moreover, for a given geometry, the matrix M defined in equation (53) needs to be computed only once.

(85) Note that in the above care should be taken to distinguish between c.sup.+ (a radial basis function parameter) and c.sup.N (processed MEG data, where N denotes the number of magnetometers).

(86) Computer Implementation

(87) Referring now to FIG. 1a, this shows a flow diagram of a procedure to numerically implement the solution of the inverse problem for EEG on a computing system. The computing system may either be a general purpose computing system, or a digital signal processor, and/or may comprise dedicated hardware.

(88) At step S100 the procedure inputs the EEG data, represented as a vector b (using the notation of Section 3 above), and at step S102 calculates an estimate of the vector by solving b=E by any of many numerical methods which will be well known to those skilled in the art. Once an estimate of has been found this can be used to determine an estimate for , from the radial basis functions, which in turn represents the current distribution J. In principle output data from the procedure may be provided in the form of data representing either , or J, and/or the current distribution data may be converted into a 3D representation of the EEG data, for example mapped onto a 3D representation of the head and/or brain (cerebrum)step S114. The skilled person will appreciate that there are many techniques which may be employed to provide a 3D representation of the EEG data, including representing the data as one or more 2D slices.

(89) The procedure of FIG. 1a employs data from the matrix E as previously described, and this may either be calculated or retrieved from storage, as shown in step S101. A procedure for calculating E is described below with reference to FIG. 1c.

(90) FIG. 1b shows a flow diagram of a procedure which may be used to image a current distribution in the brain derived from MEG data. The procedure of FIG. 1b reflects the description in section 4 above. Thus at step S106 the procedure inputs magnetic field data B and, at step S110, determines an estimate of vector from the equation M=c, again using any one of many techniques which will be well known to those skilled in the art. The matrix M may again either be calculated or retrieved from memory (step S111). In embodiments the procedure calculates the parameter vector c, for example using the head (cerebrum) geometry, the previously calculated approximation for iv, and the previously described vectorial function H. Again H is dependent only on the geometry (internal and external), and conductivity, of the head, and may again be either calculated or retrieved from memory (step S109). Then an estimate of vector can be determined from M=c.

(91) The procedure then outputs (S112) data representing the vectorial function A, in particular the radial component of this function. This may be determined from a combination of vector and the radial basis functions (although in principle the output data may comprise rather than data directly representing A). Then, again, the procedure may generate a 3D representation of the solenoidal part of the current distribution determined from the measured magnetic field B.

(92) FIG. 1c shows a flow diagram of a procedure for calculating the elements of matrix E and the elements of vectorial function H, by determining current-independent harmonic functions u.sub.j as previously described in Section 2 above. Thus at step S120 the procedure inputs geometry and conductivity data for the head, in particular defining boundaries of the cerebrum, cerebrospinal fluid, bone and scalp, and their respective electrical conductivities. Standard values may be used for the electrical conductivities. The relevant geometry data may be derived, for example, from a magnetic resonance imaging scan of the patient and/or a more generic model may be employed.

(93) At step S122 the procedure solves a set of partial differential equations, as outlined in Section 2 above, to determine the set of current-independent harmonic functions u.sub.j for each of the cerebrum, CSF, bone and scalp. From these the related functions v.sub.j may then be determined (although this step is not essential). The procedure then computes the elements of matrix E (S124), as described in Section 3 above and/or may compute the elements of the vectorial function H as described in Section 4 above. One or both of matrix E and vector function H may then be stored for later use.

(94) FIG. 1d shows an embodiment of a computer system 200 programmed to implement the methods of FIGS. 1a to 1c. Thus in the illustrated embodiment the computer system comprises a processor 204 coupled to working memory 206 and to non-volatile memory 208 storing processor control code to control the process of 204 to implement the previously described procedures. A graphical use interface 212 is provided for displaying a 3D representation of the determined current distribution(s) to the user, although the skilled person will appreciate that additionally or alternatively other forms of storage/output may be employed. The computer system 200 has an input 202 from EEG and/or MEG apparatus; this may, for example, comprise a connection over a computer network. In a similar manner the computer system 200 has access to a data store 210, again optionally over a computer network. The optional but preferable data store 210 may be employed to store one, some or all of the parameters described with reference to FIG. 1c for use in solving the inverse problem for EEG and/or MEG, in particular brain geometry/conductivity data. Optionally a pre-computed inverse matrix for EEG and/or MEG may be stored or, for example, a pre-computed inverse matrix for each of a number of topologies may be provided as a library, for example indexed by user.

(95) Although an example embodiment has been described which employs a general purpose computing system, the skilled person will appreciate from the foregoing discussion that the final imaging step reduces to a matrix by vector multiplication (the matrix inverse may be pre-computed). Thus it will be appreciated that once the procedure of FIG. 1c has been run, the remaining processing load is small and may be run in, say, a dedicated application specific integrated circuit (ASIC) and/or on a mobile device such as a tablet computer or smart phone. A data store storing the pre-computed inverse matrix may be at a remote location accessible (with suitable security) via a wired or wireless network.

5. A Numerical Estimate of vs(r,), and H(r,)

(96) The preceding approach was validated by performing a forward calculation and checking that the solution to the inverse problem agreed with the initial data. In this section we compare the analytic expression for v.sub.s(r,) obtained from equation (24) with a numerical estimate obtained via the approach outlined in section 2.2. We take the following values for the radii, and the conductivities:

(97) TABLE-US-00001 r.sub.c = 0.071m r.sub.f = 0.074m r.sub.b = 0.079m r.sub.s = 0.085m .sub.c = 0.33 .sub.f = 1 .sub.b = 0.0125 .sub.s = 0.33

(98) For the source, we let =[0.0226, 0.0068, 0.0258]m,Q=[1,1,1]/{square root over (3)}. We use a numerical implementation based on a boundary element method as outlined in Z. Akalin-Akar ibid. We compute u.sub.s(r,) at 1026 discrete points on the outer surface r=0.085 m. The results for a 256 point subset of arbitrary discrete points on the outer surface are presented in FIG. 2a. The circles show the analytic solution, and the solid line with dots is the solution from the boundary element method. The small discrepancy between the numerical and the analytic solution varies between numerical solvers, and it depends on several factors including the mesh, element types, the numerical integration, and specific assumptions used in the implementation. The comparison of the analytical solution for v.sub.s(r,) for a 256 point subset on the outer surface is shown in FIG. 2b. The circles show the analytical solution given by equation (24); the solid line with dots is the numerical solution outlined in section 2.2 above.

(99) The inverse MEG formulation relies on the complete solution of the boundary value problem described by equations (10)-(13). We consider the multi-layer spherical geometry described below:

(100) TABLE-US-00002 r.sub.c = 0.071m r.sub.f = 0.074m r.sub.b = 0.079m r.sub.s = 0.085m .sub.c = 0.33 .sub.f = 1 .sub.b = 0.0125 .sub.s = 0.33

(101) We present a solution for a finite set of observation points, and an arbitrary source point. We introduce the following notation:
.sub.i.sup.:=.sub.i{circumflex over (n)}.sub.i ic,,b

(102) where {circumflex over (n)} is a unit vector normal to the surface .sub.i. In our numerical tests we set =10.sup.3. This corresponds to a small movement towards inside () and outside (+) of the domain .sub.i. We solve v.sub.c(r,) on the surface .sub.c.sup., and v.sub.(r,) on the surface .sub.c.sup.+. For the surface .sub., we solve v.sub.(r,) on the surface .sub..sup., and v.sub.b(r,) on the surface .sub..sup.+. For the surface .sub.b, we solve v.sub.b(r,) on the surface .sub.b.sup., and v.sub.s(r,) on the surface .sub.b.sup.+. For the surface .sub.s, we present the analytical and numerical solution of the function v.sub.s(r,) in FIG. 3. We consider an arbitrary source point =[0.0226, 0.0068, 0.0258]m and 140 arbitrary discrete observation points r on each surface. It is clear that the solution at the observation points on surfaces closer to the source point .sub.c, namely .sub.i.sup., have a slightly larger amplitude than the solution on surfaces .sub.i.sup.+. The solutions converge as .fwdarw.0.

(103) FIG. 2a shows a comparison of the analytical solution and of a numerical solution of u.sub.s(r,) which is described by equations (15)-(18). The radii are [0.071, 0.074, 0.079, 0.085]m, and the corresponding conductivities are [0.33, 1, 0.0125, 0.33]. A dipole described by =[0.0226, 0.0068, 0.0258]m, Q=[1,1,1]/{square root over (3)} is used as a source. The line with circles is the analytical solution and the solid dot line is a numerical solution.

(104) FIG. 2b shows a comparison of the analytical solution and of a numerical solution of v.sub.s(r,) which is described by equations (10)-(13). The line with circles is the analytical solution, equation (24), and the solid line with dots is the solution using the approach outlined in section 2.2 above. The radii are [0.071, 0.074, 0.079, 0.085]m, and the corresponding conductivities are [0.33, 1, 0.0125, 0.33]. A source position has been fixed at =[0.0226, 0.0068, 0.0258]m. A total of 256 arbitrary points have been selected on the outer surface.

(105) FIG. 3 shows the solution of the boundary value problem described by equations (10)-(13) for every surface of the multi-layer spherical geometry. The radii are [0.071, 0.074, 0.079, 0.085]m, and the corresponding conductivities are [0.33, 1, 0.0125, 0.33]. A source position has been fixed at =[0.0226, 0.0068, 0.0258]m. For every surface, 140 arbitrary points have been selected. FIG. 3a shows v.sub.c(r,) and v.sub.(r,) on the surface .sub.c. FIG. 3b shows v.sub.(r,) and v.sub.b(r,) on the surface .sub.. FIG. 3c shows v.sub.b(r,) and v.sub.b(r,) on the surface .sub.b. FIG. 3d shows the analytical and numerical solution of v.sub.b(r,) on the surface .sub.s.

(106) FIG. 4 shows the elements of the matrix .sub.H (r,) given by equation (39) for the multi-layer spherical geometry. The radii are [0.071, 0.074, 0.079, 0.085]m, and the corresponding conductivities are [0.33, 1, 0.0.125, 0.33]. A source position has been fixed at =[0.0397, 0.0284, 0.0057]m. On the subfigures, the x-axis values are sensor indexes and the y axis shows the functions

(107) $\frac{H_j}{{}_i}.$
The solid line with circles is the analytical solution and the solid line with dots is the numerical solution. Each of the nine functions is identified by a numerical header. They are as follows:

(108) $\begin{array}{cc}\frac{H_1}{{}_1},& \left(1\right)\\ \frac{H_1}{{}_2},& \left(2\right)\\ \frac{H_1}{{}_3},& \left(3\right)\\ \frac{H_2}{{}_1},& \left(4\right)\\ \frac{H_2}{{}_2},& \left(5\right)\\ \frac{H_2}{{}_3},& \left(6\right)\\ \frac{H_3}{{}_1},& \left(7\right)\\ \frac{H_3}{{}_2},& \left(8\right)\\ \frac{H_3}{{}_3}.& \left(9\right)\end{array}$
The results of FIGS. 2-4 demonstrate that the previously described approach provides an accurate numerical solution to the inverse problem.

6. Numerical Result: Three Layer Head Model

(109) In the case of the three layer realistic head model, we have 86 electrodes. We consider a radial basis function expansion with 75 coefficients. The head geometry is shown in FIG. 5. Thus FIG. 5 shows the three layer realistic head model used for the results described herein. In more detail, FIGS. 5a to 5c show, respectively, the cerebrum, .sub.c, the skull, .sub.b and the head, .sub.s. The BEM (boundary element method) discretization involves a total of 2476 vertices and 4940 triangles.

(110) 6.1 Numerical Estimate of the Current Via EEG

(111) We consider the synthetic function () defined by

(112) $\begin{array}{cc}\left(\right)=\underset{l=1}{\overset{L}{.Math.}}\underset{i=1}{\overset{2}{.Math.}}{k}_i{h}_i^l\left(\frac{1}{l}-\ln \left[\frac{}{c_1}\right]\right)c_1^{2l+1}\frac{{\left(2l+1\right)}^2}{4}{P}_l\left({}_i.Math.\hat{}\right)& \left(56\right)\end{array}$

(113) We set L=50 and use the following values:

(114) $h_1=0.3,{}_1=\left[\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}}\right],k_1=0.029,h_2=0.4,{}_2=\left[\frac{1}{\sqrt{2}},0,-\frac{1}{\sqrt{2}}\right],k_2=0.015.$

(115) We consider a radial basis function (RBF) expansion of the function (). All information about the function () is encoded into the parametrization {(.sub.j, .sub.j):1jN,(r)=(r.sup.2+c.sup.2).sup.1/2}. The parameters of interest are the components of the vector . In order to avoid the inverse crime, we choose a different discretization of equation (4) to generate the synthetic data than the one used to perform the reconstruction. The RBF expansion of the function (56) is shown in FIG. 6. We have used 75 equally spaced points in the cerebrum, i.e N=75 in equation (47). It is clear from FIG. 6 that this expansion is sufficiently accurate to represent the synthetic function (). We assume a signal to noise ratio (SNR) of 20 dB, and model the noise as white Gaussian noise i.e wN(0,.sup.2I). The reconstructions are shown in FIG. 7.

(116) Thus FIG. 6 shows the RBF (radial basis function) expansion of the function (6.1), showing the synthetic function () and the corresponding radial basis function interpolation of () inside the domain .sub.c of the three layer head model. We used 75 coefficients and sampled the function on a Cartesian grid inside the domain .sub.c. The top figure shows the function () defined by equation (56) as well as the corresponding radial basis function interpolation. The bottom figure shows the absolute error between the exact () and the corresponding interpolation.

(117) FIG. 7 shows the exact and reconstructed functions () for the three Layer head model. We have used 75 coefficients and 86 electrodes. The signal to noise ratio (SNR) is 20 dB. In FIG. 7a, the solid line with circles shows the exact , equation (47), whereas the solid line with dots is the reconstructed {circumflex over ()}, equation (50). In FIG. 7b, the solid line with circles is () as given by equation (56), whereas the solid line with dots is the reconstructed current computed using the estimated {circumflex over ()}. FIG. 7c shows the error e():=|(){circumflex over ()}() of the functions shown in FIG. 7b.

(118) 6.2 Numerical Estimate of the Current Via MEG

(119) We consider the synthetic function used in Fokas, Inverse Problems, ibid, to represent A.sup.(), but with different numerical values. The function A.sup.() is given by

(120) 0$\begin{array}{cc}{A}^{}\left(\right)=\underset{i=1}{\overset{2}{.Math.}}{k}_i\left[\frac{1}{4}\frac{1-h_i^2}{{\left(1+h_i^2-2h_i\left({}_i.Math.\hat{}\right)\right)}^{1.5}}-\frac{1}{4}\right].& \left(57\right)\end{array}$

(121) We use the following values:

(122) $h_1=0.3,{}_1=\left[\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}}\right],k_1=0.029,h_2=0.4,{}_2=\left[\frac{1}{\sqrt{2}},0,-\frac{1}{\sqrt{2}}\right],k_2=0.0047.$

(123) As in the case with EEG, we employ a RBF expansion for the function A.sup.(). All information about the function A() is encoded into the parametrization {(.sub.j, .sub.j): 1jN, (r)=(r.sup.2+c.sup.2).sup.1/2}. The parameters of interest are the components of the vector . In order to avoid the inverse crime we choose a different discretization of equation (5) to generate the synthetic data than the one used to perform the reconstruction. The RBF expansion of the function (57) is shown in FIG. 8. We used 75 equally spaced points in the cerebrum, i.e N=75 in equation (52). It is clear from FIG. 8 that this expansion is sufficiently accurate to represent the synthetic function A.sup.(). We assume a signal to noise ratio (SNR) of 20 dB and model the noise as white Gaussian noise i.e wN(0,.sup.2I). We have 102 real sensor positions to obtain the measurements. The reconstructions are shown in FIG. 9.

(124) Thus FIG. 8 shows the synthetic function A.sup.() and the corresponding radial basis function interpolation of A.sup.() inside the domain .sub.c of the three layer head model. We used 75 coefficients and sampled the function on a Cartesian grid inside the domain .sub.c. The top figure shows the function A.sup.96 () as given by equation (57), as well as the corresponding radial basis function interpolation. The bottom figure shows the absolute error between the exact A.sup.() and the corresponding interpolation.

(125) FIG. 9 shows the exact and reconstructed A.sup.() and the three layer head model. We used 75 coefficients and 102 real sensor positions. The signal to noise ratio (SNR) is 20 dB. In FIG. 9a, the solid line with circles shows the exact , equation (52), whereas the solid line with dots is the reconstructed {circumflex over ()} as given by equation (57). In FIG. 9b, the solid line with circles is A.sup.() as given by equation (57), whereas the solid line with dots is the reconstructed current computed using the estimated {circumflex over ()}. FIG. 9c shows the error e():=|A.sup.().sup. ()| of the functions shown in FIG. 9b.

7. Conclusion

(126) We have introduced an algorithmic approach for extracting the maximum possible information about the neuronal current J.sup.P(), .sub.c, from the measurement of the electric potential u.sub.s(r), r.sub.s, on the scalp, and from the measurement of the radial part of the magnetic field r.Math.B(r) r.sub.e in the exterior of the head, relying on a numerical approach for computing the auxiliary functions .sub.v.sub.s(r,) and .sub.H(r,). In section 2.1 we derived analytic expressions for v.sub.s(r,) and .sub.H(r,) for the particular case of the spherical geometry. This allowed us to verify the effectiveness of the above numerical codes; the agreement with the analytical expressions is excellent. In section 6 we implemented our general approach to a realistic three shell geometry. The above algorithm yields accurate reconstructions both for the component () of the current from the data set {u.sub.s (r.sub.i):1iN.sub.e}, as well as for the radial component of the vectorial function A() (denoted by A.sup. ()) from the data set {r.sub.k.Math.B (r.sub.k):1iN.sub.s, r.sub.k .sub.e} and the estimated function (). Indeed, starting with the particular function () given by equation (56), we generate a set of synthetic data for the electric potential on the scalp; using this set we show that we can reconstruct the function () even in the presence of noise. Similarly, starting with the functions () and A.sup.() given by equations (56) and (57) respectively, we generate a set of synthetic data for the radial part of the magnetic field outside the head; then, using this set as well as the reconstructed function (), we show that we can reconstruct the radial part of A() even in the presence of noise. The assumed noise model in both the EEG and MEG is additive white Gaussian with a signal to noise ratio of 20 dB.

(127) It can thus be seen that embodiments of the above described algorithm are well-suited to address the distributed source problem, in particular in real time, and can handle reasonable noise levels in the measurements.

(128) EEG can yield as much information about the neuronal current as MEG, but at a fraction of the cost since EEG devices are inexpensive and widely available. A further advantage of EEG (and MEG) by comparison with other techniques such as fMRI, PET and the like is the fast time resolutionEEG and MEG can make 10s or 100s of measurements per millisecond. Furthermore, it is straightforward to reconstruct images of the current distributions from EEG data and MRI data (for the head geometry).

(129) Thus a further aspect of the invention contemplates an upgrade for a patient imaging machine, such as a magnetic resonance imaging machine, provided by combining the machine with EEG (and/or MEG) signal measurement apparatus and EEG (and/or MEG) signal processing software or hardware as described above. Such an upgrade may be of substantial benefit in diagnosis and treatment of many neurological diseases including, for example, epilepsy. In another application such an approach may be employed to differentiate mild cognitive impairment (a transitional stage between normal aging and Alzheimer's disease) from Alzheimer's disease proper.

(130) EEG technology is also finding its way into the consumer marketfor example low cost EEG headsets are becoming available for games. Thus a further application of the embodiments of the invention uses processing as described above to provide data for performing an action and/or controlling a consumer electronic device, for example for controlling a computer game or the like. More generally, embodiments of the techniques we describe may be employed to detect a user's thoughts and/or feelings and/or expressions in real time.

(131) The skilled person will appreciate that numerous variations of the above described techniques are possible. For example, although we have described the use of radial basis functions for expanding the neuronal current other basis functions may also be used. Either or both of the inverse EEG and inverse MEG algorithms may be partially, or substantially wholly, parallelised.

(132) Optionally Bayesian filtering and/or dynamic parameter estimation or other time series analysis may be applied to either or both of the raw data and processed data. The skilled person will appreciate that the EEG and/or MEG signals used may be pre-processed in many different ways, for example by filtering to attenuate noise, using any of a range of techniques. Similarly there are many different ways in which a graphical user interface may be provided to operators of the system, for example to facilitate a user in conducting a real-time experiment. Against the skilled person will appreciate that there are many ways in which data from the procedures we have described may be post-processed.

(133) No doubt many other effective alternatives will occur to the skilled person. It will be understood that the invention is not limited to the described embodiments and encompasses modifications apparent to those skilled in the art lying within the spirit and scope of the claims appended hereto.