# Vector tomography for reconstructing electric fields with non-zero
divergence in bounded domains (with application in EEG Source Imaging)

Alexandra Koulouri<EMAIL_ADDRESS>Mike Brookes  Ville Rimpiläinen
Institute for Computational and Applied Mathematics, University of Münster,
Einsteinstrasse 62, 48149 Münster, Germany  Department of Electrical and
Electronic Engineering, Imperial College London, Exhibition Road, London SW7
2BT, United Kingdom  Department of Mathematics, University of Auckland,
Private bag 92019, Auckland 1142, New Zealand

###### Abstract

In vector tomography (VT), the aim is to reconstruct an unknown multi-
dimensional vector field using line integral data. In the case of a
2-dimensional VT, two types of line integral data are usually required. These
data correspond to integration of the parallel and perpendicular projection of
the vector field along integration lines. VT methods are non-invasive, non-
intrusive and offer more information on the field than classical point
measurements; they are typically used to reconstruct divergence-free (or
source-free) velocity and flow fields. In this paper, we show that VT can also
be used for the reconstruction of fields with non-zero divergence. In
particular, we study electric fields generated by dipole sources in bounded
domains which arise, for example, in electroencephalography (EEG) source
imaging. To the best of our knowledge, VT has not previously been used to
reconstruct such fields. We explain in detail the theoretical background, the
derivation of the electric field inverse problem and the numerical
approximation of the line integrals. We show that fields with non-zero
divergence can be reconstructed from the longitudinal measurements with the
help of two sparsity constraints that are constructed from the transverse
measurements and the vector Laplace operator. As a comparison to EEG source
imaging, we note that VT does not require mathematical modelling of the
sources. By numerical simulations, we show that the pattern of the electric
field can be correctly estimated using VT and the location of the source
activity can be determined accurately from the reconstructed magnitudes of the
field.

###### keywords:

Vector tomography , electric field , Radon transform , line integral , inverse
problem , sparsity constraint

††journal: Journal of Computational Physics

##  1 Introduction

Vector fields such as gravitational and electromagnetic fields are fundamental
objects of study in physics. Vector tomography (VT) is a framework that can be
used to reconstruct unknown vector fields using line integral measurements
[50, 53]. The longitudinal line measurements are performed by projecting the
studied field on lines that trace the domain and then integrating the
projected field along the lines. The transverse line measurement are carried
out similarly but now the field components that are perpendicular to these
lines are integrated. VT methods are attractive because they are non-invasive,
non-intrusive and give larger amount of data compared to the one-sensor one-
measurement set-up.

VT studies have been carried out both in theoretical level and applications,
concentrating mainly on the reconstruction of smooth vector fields [53].
Theoretical analyses for the reconstruction of smooth velocity fields have
been reconstructed [40, 41, 52, 32, 8, 38, 48, 55, 36, 28] using such methods
as the inverse Radon transform [24], the inverse Fourier transform with
central slice theorem [39, 13] and back projection (parallel beam tomography)
[37, 51]. The VT framework has been used for the reconstruction of particle
distributions [5], ion fields in plasma [16, 26, 3, 4], velocity fields in
blood veins [32, 29, 53], magnetic field of the corona of the sun [36], Kerr
effect in optical polarization tomography [25] and micro-structures in
oceanographic tomography [48]. Both linear and non-linear iterative algorithms
have been proposed for vector functions with appropriate smoothness [55, 39,
28, 46, 17, 35, 53].

###  1.1 Unbounded domain

The theoretical basis for reconstructing smooth vector fields that decay
sufficiently rapidly to zero was introduced in [40]. Based on Helmholtz’s
theorem [2], vector fields can be decomposed as a sum of irrotational (curl-
free) and solenoidal (divergence-free or source-free) components and it was
first shown that, for a 2-dimensional field, the solenoidal component can be
imaged with the help of longitudinal measurements [40]. It was subsequently
shown that the transverse measurements were required in order to recover the
remainder of the field [8].

The problem was extended to three dimensions in [47] using the formalism of
the 3-dimensional (3D) vector Radon transform. First, a generalization of the
integral measurement was introduced. It was called the probe transform (or
general product measurement) and it was formulated as an inner product between
the Radon transform of a field and a unit-vector in a specific direction. It
was also shown that three types of measurements were required for the recovery
of a 3D field. In [53], the analysis was generalized to multidimensional
cases.

However, in most practical situations, it is difficult or even impossible to
perform the transverse measurements (i.e. the probe transform in the
transverse direction). For example in Doppler techniques [15] or in geophysics
[41], this type of measurement is not physically realizable. In fact, the
transverse measurements can be obtained only in very specific set-ups [8, 33].

###  1.2 Bounded domain

In practical applications, vector fields are investigated in a bounded domain
where the field is non-zero on the boundary, and it is often the boundary that
partially defines the field itself. The VT framework was extended to non-
homogenous boundary conditions in [8, 44]. The field decomposition included an
additional harmonic field component that was used to satisfy the boundary
conditions [8]. Based on the inverse Radon transform, it was show that the
solenoidal component and only one half of the harmonic component can be imaged
with the longitudinal measurement whereas the transverse measurement gives the
irrotational component and the other half of the harmonic component [8]. In
[44], the results were generalized. In particular, it was shown that the
longitudinal measurement can be used to recover both the solenoidal component
and the part of the harmonic component that arises from the field component
that is tangential to the boundary. Additionally, transverse measurement
reconstructs the irrotational component and the harmonic part that arises from
the field component that is normal to the boundary.

###  1.3 Electric field with non-zero divergence

All the previous studies have investigated divergence-free vector fields.
However, also vector fields with non-zero divergence are common in physics.
This kind of fields can be, for example, gravitational or electromagnetic
fields that are generated by unknown sources (and/or sinks) that are located
inside the domain of interest. In this paper, the aim is to use VT to
reconstruct such non-zero divergence vector fields. In particular, we employ
VT for the reconstruction of electric fields in a bounded domain that includes
a dipole source. Strategies to estimate such electric activity is of great
interest especially in medical imaging modalities such as
electroencephalography (EEG) in which the imaging problem is traditionally
parametrisized using source spaces [22, 56]. The proposed VT modelling assumes
the same physical conditions as the dipole source imaging problem i.e. the
underlying electric field is irrotational. The existence of a dipole inside
the domain implies that the field has a singularity, and it has been shown
that VT can be used in such cases [14].

The use of VT rather than traditional inverse source approaches [21, 42]
offers two advantages. First, the continuous VT problem for the recovery of
the electric field using a set of line integrals is only a moderately ill-
posed problem [39] whereas the inverse source problem is a severely ill-posed
problem that cannot be solved from boundary measurements without a priori
knowledge [1]. In practice, however, prior information is also required by the
VT formulation (e.g. introduced as a penalty term) in order to obtain a stable
solution since only a finite number of measurements is available for the
reconstruction (incomplete data problem [39]). Second, the VT approach does
not require an explicit mathematical model of the underlying sources. For
example, in EEG source imaging there is an extensive literature on different
mathematical models of neural sources [31, 30, 9, 27, 49, 42, 6].

In the proposed VT approach, we use a set of line integrals that trace a
conductive 2-dimensional domain and result in a linear system of equations. We
show that the longitudinal measurements are determined by the electric
potentials on the domain boundary and by employing the vectorial Radon
properties that the transverse measurements give information on the underlying
current sources. We describe in detail the theoretical background, the
numerical approximation of the line integrals and finally present the
discretized electric field inverse problem which is solved with the help of
the $L_{1}$-norms of the transverse measurements and the discrete vector
Laplace operator. The resulting non-linear minimization problem is solved
using convex optimization. Finally, we show by numerical simulations that
electric fields with non-zero divergence can be reconstructed in a bounded
domain under a homogeneous Neumann condition.

##  2 Mathematical preliminaries

In this section, we explain the notation used and define the different Radon
measurements. More general information on the Radon transform can be found
from [24, 13].

###  2.1 Radon transform of a scalar function

We denote by $\mathrm{x}\in\mathbb{R}^{2}$ a point, $f(\mathrm{x})$ an
arbitrary scalar function and
$L(l,\mathrm{\hat{s}}^{\perp}):=\\{\mathrm{x}=(x,y)\in\mathbb{R}^{2}:\mathrm{x}\cdot\mathrm{\hat{s}}^{\perp}=l\\}$
a line where ${l}$ is the signed distance of the line from origin and
$\mathrm{\hat{s}}^{\perp}=(\cos\phi,\sin\phi)$ is the unit normal vector of
$L$ (see, Figure 1 for details).

Figure 1: The integration path $L$ is defined with the help of the signed
distance from the origin $l$ and the unit normal vector $\hat{s}^{\perp}$. In
addition, $\hat{s}$ denotes the unit vector along $L$.

The mapping defined by the projection (or the line integral) of
$f(\mathrm{x})$ along a line $L$ is the two-dimensional Radon transform of $f$
and is given by

$\tilde{f}=\mathcal{R}\\{f\\}(l,\hat{s}^{\perp})=\int_{\mathrm{x}\in
L}f(\mathrm{x})~{}d\ell(\mathrm{x}),$ (1)

where $d\ell(\mathrm{x})$ is an increment of length along $L$ [39].

The corresponding scalar inverse Radon transform is

$f=\mathcal{R}^{-1}\\{\tilde{f}\\}=\frac{1}{4\pi}\mathcal{R}^{\\#}\mathcal{H}\frac{\partial}{\partial
l}\tilde{f},$ (2)

where $\mathcal{R}^{\\#}$ is the adjoint operator of $\mathcal{R}$,
$\mathcal{H}\frac{\partial}{\partial l}$ is the filtration part of the inverse
transform and $\mathcal{H}$ denotes the Hilbert transform [13].

###  2.2 $\Omega$-Radon transform of a vector function

In our analysis, we operate in a bounded convex domain
$\Omega\subset\mathbb{R}^{2}$ on a vector function
$\mathrm{f(x)}=(f_{x},f_{y}):\Omega\rightarrow\mathbb{R}^{2}$. We restrict
line $L$ now on the domain $\Omega$. The $\Omega$-Radon transform
$\mathcal{R}_{\Omega}$ of $f$ is the restriction of the Radon transform to the
domain $\Omega$

$\tilde{\mathrm{f}}_{\Omega}=(\tilde{f}_{\Omega x},\tilde{f}_{\Omega
y})=\mathcal{R}_{\Omega}\\{\mathrm{f}\\}(l,\mathrm{\hat{s}}^{\perp}).$ (3)

Using the Dirac-delta function, the Radon integral can be expressed as the
surface integral

$\tilde{\mathrm{f}}_{\Omega}=\mathcal{R}_{\Omega}\\{\mathrm{f}\\}(l,\mathrm{\hat{s}}^{\perp})=\int_{\Omega}\mathrm{f}(\mathrm{x})\delta(l-\mathrm{x}\cdot\mathrm{\mathrm{\hat{s}}}^{\perp})~{}d\mathrm{x}.$
(4)

Additionally, we will find useful the $\Omega$-Radon transform of a
directional derivative along a unit vector $\mathrm{\hat{\omega}}$

$\mathcal{R}_{\Omega}\\{\mathrm{\hat{\omega}}\cdot\nabla\mathrm{f}\\}=\mathrm{\hat{\omega}}\cdot\mathrm{\hat{s}}^{\perp}\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}\\{\mathrm{f}\\}.$ (5)

###  2.3 Line integral data

In 2-dimensional VT, two different types of line integral measurements are
used to reconstruct vector fields. The first is the line integral

$I_{L}^{\parallel}(l,s^{\perp})=\int_{L}\hat{\mathrm{s}}\cdot\mathrm{f}(\mathrm{x})~{}d\ell(\mathrm{x})=\hat{\mathrm{s}}\cdot\mathcal{R}_{\Omega}\\{\mathrm{f}\\},$
(6)

where $\hat{\mathrm{s}}$ is the unit vector in the direction of line $L$. The
second integral is

$I_{L}^{\perp}(l,s^{\perp})=\int_{L}\hat{\mathrm{s}}^{\perp}\cdot\mathrm{f}(\mathrm{x})~{}d\ell(\mathrm{x})=\hat{\mathrm{s}}^{\perp}\cdot\mathcal{R}_{\Omega}\\{\mathrm{f}\\},$
(7)

where $\hat{\mathrm{s}}^{\perp}$ is the unit vector that is orthogonal to the
line $L$. The line integral (6) is also called as longitudinal measurement and
(7) as transverse measurement [50].

##  3 Theory

###  3.1 Electric field with a current source in a bounded domain

Let us consider a bounded domain $\Omega\subset\mathbb{R}^{2}$ with linear
material properties $\sigma(\mathrm{x})$, $\mathrm{x}\in\Omega$ and that the
domain contains an electric source with (primary) current density
$\mathrm{j}^{s}(\mathrm{x}):\Omega\rightarrow\mathbb{R}^{2}$. The electric
source induces an electric field
$\mathrm{e}(\mathrm{x}):\Omega\rightarrow\mathbb{R}^{2}$. The total current
density in the medium can be presented as a sum of the primary and secondary
current [34], i.e.

$\mathrm{j}(\mathrm{x})=\mathrm{j}^{s}(\mathrm{x})+\sigma(\mathrm{x})\mathrm{e}(\mathrm{x}).$
(8)

For current signals with low frequencies, the quasi-static Maxwell equations
can be used

$\displaystyle\nabla\times\mathrm{e}(\mathrm{x})$ $\displaystyle=0$ (9a)
$\displaystyle\nabla\times\mathrm{h}(\mathrm{x})$
$\displaystyle=\mathrm{j}(\mathrm{x}),$ (9b)

where $\mathrm{h}(\mathrm{x})$ is the magnetic field intensity. The divergence
of equation (9b) gives

$\displaystyle\nabla\cdot\nabla\times\mathrm{h}(\mathrm{x})$
$\displaystyle=\nabla\cdot\mathrm{j}(\mathrm{x})$ (10a) $\displaystyle 0$
$\displaystyle=\nabla\cdot(\mathrm{j}^{s}(\mathrm{x})+\sigma(\mathrm{x})\mathrm{e}(\mathrm{x}))$
(10b) $\displaystyle\nabla\cdot\sigma\mathrm{e}(\mathrm{x})$
$\displaystyle=-\nabla\cdot\mathrm{j}^{s}(\mathrm{x}),$ (10c)

which relates the electric field to the current source.

Because the electric field is irrotational, Equation (9a), the field can be
expressed as

$\mathrm{e(x)}=-\nabla u~{}~{}~{}\mbox{in}~{}~{}\Omega,$ (11)

where $u$ is a scalar potential. In this paper, we consider that $u$ is
uniquely defined as the solution of the Poisson equation with the following
boundary conditions

$\displaystyle\nabla\cdot\sigma\nabla u$
$\displaystyle=\mathrm{j}^{s}(\mathrm{x})$ (12a)
$\displaystyle\int_{\partial\Omega}u~{}d\mathrm{x}$
$\displaystyle=0,\,\mathrm{x}\in\partial\Omega$ (12b)
$\displaystyle\frac{\partial u}{\partial\hat{n}}$
$\displaystyle=0,\,\mathrm{x}\in\partial\Omega,$ (12c)

where $\partial\Omega$ is the boundary and $\hat{n}$ is the outward unit
normal vector. The homogeneous Neumann condition (12c) means that the electric
field is tangential at the boundary, that is
$\hat{n}\cdot\mathrm{e}|_{\partial\Omega}=0$.

###  3.2 Line integrals through direct substitution

In this vector tomography framework, we consider two types of line integral
measurements of the electric field. If we directly evaluate the integrals (6)
and (7) by substituting the electric field with the negative gradient of the
scalar potential, first, the longitudinal integral measurements along the line
$L$ are

$I_{L(\mathrm{x}_{a},\mathrm{x}_{b})}^{\parallel}=\int_{L(\mathrm{x}_{a},\mathrm{x}_{b})}\mathrm{e}(\mathrm{x})\cdot\hat{\mathrm{s}}~{}d\ell(\mathrm{x})=\int_{L(\mathrm{x}_{a},\mathrm{x}_{b})}-\nabla
u\cdot\hat{\mathrm{s}}~{}d\ell(\mathrm{x})=u(\mathrm{x}_{a})-u(\mathrm{x}_{b}),$
(13)

where $u(\mathrm{x}_{a})$ and $u(\mathrm{x}_{b})$ are the electric potential
values at the intersections of the line $L$ and the boundary $\partial\Omega$
and $\hat{\mathrm{s}}=(s_{x},s_{y})$ is the unit vector along the line.
Second, the transverse line integral measurements are

$I_{L(\mathrm{x}_{a},\mathrm{x}_{b})}^{\perp}=\int_{L(\mathrm{x}_{a},\mathrm{x}_{b})}\mathrm{e}(\mathrm{x})\cdot\hat{{\mathrm{s}}}_{\perp}~{}d\ell(\mathrm{x})=\int_{L(\mathrm{x}_{a},\mathrm{x}_{b})}\frac{\partial
u}{\partial y}s_{x}-\frac{\partial u}{\partial x}s_{y}~{}d\ell(\mathrm{x}).$
(14)

As it can be seen, the longitudinal integral measurements are directly
determined by the boundary potentials; however, the transverse integral
measurements in practice cannot be measured directly (or only under special
circumstances [8, 33]). This turns out to be a problem because the full
recovery of the electric field requires both types of integral measurements.
In the Appendix, it is shown that only the harmonic component of the electric
field can be reconstructed from the longitudinal line integrals, whereas the
remaining irrotational part requires the transverse measurements. The
transverse integral formulations nevertheless turn out to be useful since, as
will be seen in the following section, they give information about the
underlying current sources that generate the field.

###  3.3 Transverse line integral through $\Omega$-Radon transform

We begin with a similar analysis to that in [44] and we formulate the
transverse line integral measurements using the properties of the Radon
transform. For simplicity, the electric conductivity $\sigma$ of the domain is
considered constant. Let $\mathrm{e}^{c}=(e^{c}_{x},e^{c}_{y})$, defined on
$\mathbb{R}^{2}$, be such a vector-valued function such that
$\mathrm{e}^{c}=\mathrm{e}$ on $\Omega$. Then the $\Omega$-Radon transform of
$\mathrm{e}$ can be expressed as

$\tilde{\mathrm{e}}_{\Omega}=\mathcal{R}_{\Omega}\\{\mathrm{e}\\}(l,\mathrm{\hat{s}}^{\perp})=\mathcal{R}\\{\mathrm{e}^{c}v_{\Omega}\\}(l,\mathrm{\hat{s}}^{\perp})=\int_{\mathbb{R}^{2}}\mathrm{e}^{c}(\mathrm{x})~{}v(h(x))~{}\delta(l-\mathrm{x}\cdot\mathrm{\hat{s}}^{\perp})~{}d\mathrm{x},$
(15)

where $v:\mathbb{R}\rightarrow\mathbb{R}$ is an indicator (Heaviside) function
defined as

$v(s)=\begin{cases}1,&s>0\\\ 0,&s\leq 0\end{cases}$ (16)

and $h(\mathrm{x})$ is a signed distance function satisfying

$h(\mathrm{x})\begin{cases}>0,&\mathrm{x}\in\Omega\\\
=0,&\mathrm{x}\in\partial\Omega\\\
<0,&\mathrm{x}\notin\bar{\Omega}.\end{cases}$ (17)

For simplicity, we define $v_{\Omega}=v(h(\mathrm{(}x)))$. For the signed
distance, we have that $\nabla h(\mathrm{x})=-\hat{n}(\mathrm{x})$ on a
(piecewise) smooth boundary $\partial\Omega$, where $\hat{n}$ is the outward
unit normal vector [43].

We estimate $\nabla\cdot$ and $\nabla\times$ of the field
$\mathrm{e}^{c}v_{\Omega}$

$\displaystyle\nabla\cdot(\mathrm{e}^{c}v_{\Omega})$ $\displaystyle=$
$\displaystyle(\nabla\cdot\mathrm{e}^{c})~{}v_{\Omega}+(\nabla
h)\cdot\mathrm{e}^{c}~{}\delta_{\Omega}$ (18)
$\displaystyle\nabla\times(\mathrm{e}^{c}v_{\Omega})$ $\displaystyle=$
$\displaystyle(\nabla\times\mathrm{e}^{c})~{}v_{\Omega}+(\nabla
h)\times\mathrm{e}^{c}~{}\delta_{\Omega},$ (19)

where $\delta_{\Omega}=\delta(h(\mathrm{x}))$. First, we notice that due to
the irrotational property of the field, $\nabla\times\mathrm{e}^{c}=0$. Then,
from the gradient of the signed distance function and the boundary conditions
(12b) and (12c), we get $\nabla h\cdot
e|_{\partial\Omega}=-\hat{n}\cdot\mathrm{e}|_{\partial\Omega}=\frac{\partial
u}{\partial\hat{n}}|_{\partial\Omega}=0$. Now, we can re-write Equations (18)
and (19) as

$\displaystyle\nabla\cdot(\mathrm{e}^{c}v_{\Omega})$ $\displaystyle=$
$\displaystyle(\nabla\cdot\mathrm{e}^{c})~{}v_{\Omega}$ (20)
$\displaystyle\nabla\times(\mathrm{e}^{c}v_{\Omega})$ $\displaystyle=$
$\displaystyle(\nabla h)\times\mathrm{e}^{c}~{}\delta_{\Omega}.$ (21)

Employing the Radon properties of the derivatives [39], the Radon transforms
of $\nabla\cdot(\mathrm{e}^{c}v_{\Omega})$ and
$\nabla\times(\mathrm{e}^{c}v_{\Omega})$ are

$\displaystyle\mathcal{R}\\{\nabla\cdot(\mathrm{e}^{c}v_{\Omega})\\}$
$\displaystyle=$
$\displaystyle\mathrm{\hat{s}}^{\perp}\cdot\frac{\partial}{\partial
l}\mathcal{R}\\{\mathrm{e}^{c}v_{\Omega}\\}=\mathrm{s}^{\perp}\cdot\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}\\{\mathrm{e}\\}$ (22)
$\displaystyle\mathcal{R}\\{\nabla\times(\mathrm{e}^{c}v_{\Omega})\\}$
$\displaystyle=$
$\displaystyle\mathrm{\hat{s}}^{\perp}\times\frac{\partial}{\partial
l}\mathcal{R}\\{\mathrm{e}^{c}v_{\Omega}\\}=\mathrm{\hat{s}}^{\perp}\times\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}\\{\mathrm{e}\\}.$ (23)

Based on [44, 2], a vector field can be expressed as
$\mathrm{f}=(\mathrm{\hat{s}}^{\perp}\cdot\mathrm{f})\mathrm{\hat{s}}^{\perp}+(\mathrm{\hat{s}}^{\perp}\times\mathrm{f})\times\mathrm{\hat{s}}^{\perp}$.
When we substitute $\mathrm{f}=\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}\\{\mathrm{e}\\}$ and use Equations (20) and (21), we
get

$\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}\\{\mathrm{e}\\}=\mathcal{R}\\{(\nabla\cdot\mathrm{e}^{c})~{}v_{\Omega}\\}\mathrm{\hat{s}}^{\perp}+\mathcal{R}\\{(\nabla
h)\times\mathrm{e}^{c}~{}\delta_{\Omega}\\}\times\mathrm{\hat{s}}^{\perp}.$
(24)

From the dot-product of $\hat{s}^{\perp}$ and Equation (24) we get

$\mathrm{s}^{\perp}\cdot\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}\\{\mathrm{e}\\}=\mathcal{R}\\{\nabla\cdot\mathrm{e}^{c}~{}v_{\Omega}\\}=\mathcal{R}_{\Omega}\\{\nabla\cdot\mathrm{e}\\}.$
(25)

By using the relation
$\nabla\cdot\mathrm{e}=-\frac{1}{\sigma}\nabla\cdot\mathrm{j}^{s}$, Radon
properties [39] and the assumption that the (extended) source function is zero
outside the domain $\Omega$, we obtain

$\mathrm{s}^{\perp}\cdot\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}\\{\mathrm{e}\\}=-\mathrm{\hat{s}}^{\perp}\cdot\frac{1}{\sigma}\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}\\{\mathrm{j}^{s}\\}.$ (26)

Now, the inverse Radon transform (2) gives us

$\displaystyle\mathcal{R}^{-1}\\{\mathrm{\hat{s}}^{\perp}\cdot\mathrm{\tilde{e}_{\Omega}}\\}$
$\displaystyle=\frac{1}{4\pi}\mathcal{R}^{\\#}\mathcal{H}\frac{\partial}{\partial
l}[\mathrm{\hat{s}}^{\perp}\cdot\mathcal{R}_{\Omega}\\{\mathrm{e}\\}]$ (27a)
$\displaystyle=-\frac{1}{4\pi\sigma}\mathcal{R}^{\\#}\mathcal{H}\frac{\partial}{\partial
l}[\mathrm{\hat{s}}^{\perp}\cdot\mathcal{R}_{\Omega}\\{\mathrm{j}^{s}\\}]$
(27b)
$\displaystyle=-\frac{1}{\sigma}\mathcal{R}^{-1}\\{\mathrm{\hat{s}}^{\perp}\cdot\tilde{\mathrm{j}}^{s}\\}.$
(27c)

Therefore, we have that the transverse measurements are

$I_{L}^{\perp}=\mathrm{\hat{s}}^{\perp}\cdot\mathcal{R}_{\Omega}\\{\mathrm{e}\\}=-\mathrm{\hat{s}}^{\perp}\cdot\frac{1}{\sigma}\mathcal{R}_{\Omega}\\{\mathrm{j}^{s}\\}.$
(28)

In other words, we have shown that the transverse integral measurements give
us information about the source activity inside the bounded domain.

###  3.4 Dipole sources

In this paper, we consider focal source activity and dipoles in particular.
Dipole sources can be described with the help of Dirac delta functions as

$\mathrm{j}^{s}(\mathrm{x})=\sum_{i=1}^{N_{s}}\mathrm{q}_{i}\delta(\mathrm{x}-\mathrm{x}_{i}),$
(29)

where $\mathrm{q}_{i}$ is the dipole moment, $\mathrm{x}_{i}$ the dipole
source location and $N_{s}$ the total number of dipole sources. When this
expression is substituted into Equation (28), we obtain

$I_{L}^{\perp}=-\frac{1}{\sigma}\sum_{i=1}^{N_{s}}\int_{\Omega}\mathrm{\hat{s}}^{\perp}\cdot\mathrm{q}_{i}\delta(\mathrm{x}-\mathrm{x}_{i})\delta(l-\mathrm{\hat{s}}^{\perp}\cdot\mathrm{x})~{}d\mathrm{x}=-\frac{1}{\sigma}\sum_{i=1}^{N_{s}}\mathrm{\hat{s}}^{\perp}\cdot\mathrm{q}_{i}\delta(l-\mathrm{\hat{s}}^{\perp}\cdot\mathrm{x}_{i}).$
(30)

Equation (30) implies that the transverse measurement is non-zero only when
the line of integration passes through the source location and the line is not
parallel to the dipole moment, such as the black line in Figure 2. In VT, we
have a set of lines with different directions and only few of them meet these
criteria. This knowledge that only few of the transverse measurements are non-
zero is later used as a sparsity constraint in the electric field inverse
problem.

Figure 2: Electric field (arrows) that is generated by a dipole source
(circle). Of the three shown lines of integration, the gray one is parallel
and the black one perpendicular to the orientation of the dipole. The
transverse line integral along the gray line is zero due to the dot-product in
Equation (30), however the transverse integral along the black line is non-
zero. The transverse line integral along the dashed line is zero because the
line does not pass through the source, and also the longitudinal line integral
is close to zero because the line is far from the dipole source.

###  3.5 Discrete observation model

For the numerical evaluation, the domain in discretized and the electric field
is represented as

$\mathrm{e(x)}\approx\mathrm{e}_{h}(\mathrm{x})=\sum_{i=1}^{N}\mathrm{e}_{i}\phi_{i}(\mathrm{x}),$
(31)

where $h$ refers to the discretization level, $\phi_{i}(\mathrm{x})$ are the
chosen basis functions, $N$ is the number of basis functions and
$\mathrm{e}_{i}=(e_{ix},e_{iy})$ contains the electric field components. In
the following, we denote
${e}=[e_{1x},e_{2x},\cdots,e_{Nx},e_{1y},e_{2y},\cdots,e_{Ny}]^{\mathrm{T}}$
as the vector representation of $\mathrm{e}_{h}(\mathrm{x})$.

The line integrals are evaluated along a set of straight lines which are
formed by connecting pairs of points on the boundary $\partial\Omega$. In
practice, a finite number of points (or electrodes) are used for the
measurements. If the number of electrodes is $n$, then the number of possible
line integral measurements is $m=n(n-1)/2$.

We stack the longitudinal line measurements into a vector
$I^{\parallel}=[I^{\parallel}_{1},\cdots,I^{\parallel}_{m}]^{\mathrm{T}}\in\mathbb{R}^{m}$
and present the observations in a matrix form as

$I^{\parallel}={R}^{\parallel}{e}+\varepsilon,$ (32)

where ${R}^{\parallel}\in\mathbb{R}^{m\times 2N}$ is called the longitudinal
ray matrix and it consists of the integration coefficients of the $m$
longitudinal line integrals, and $\varepsilon\in\mathbb{R}^{m}$ is the
measurement noise that is assumed to be random (Gaussian) white noise.

For the transverse line measurements, a similar matrix can be formulated. This
matrix is called the transverse ray matrix $R^{\perp}\in\mathbb{R}^{m\times
2N}$ and it consists of the integration coefficients of the $m$ transverse
line integrals. The numerical approximation of the ray matrices is described
in detail in Section 4.1.

###  3.6 Discrete electric field inverse problem

For the estimation of the field, the VT methods require the values of both
types of integrals for all possible lines. However, two difficulties arise.
The first one is that the transverse integral measurements cannot be carried
out using physical means. The second one is that in practice we have only a
finite number of lines and measurements available (limited data problem). This
means that there are areas in the domain that are not covered by any of the
line integrals, thus no information can be retrieved from these areas.

Here, we deal with these problems by using two penalty terms. We formulate the
electric field inverse problem as a minimization problem as follows

${\hat{{e}}}=\min_{{e}}\\{\|R^{\parallel}{{e}}-I^{\parallel}\|_{2}^{2}+\alpha\|R^{\perp}\mathrm{{e}}\|_{1}+\beta\|W{\mathrm{e}}\|_{1}\\}.$
(33)

The first term is the data fidelity term, the second term is the $L_{1}$-norm
of the transverse line integral measurements with a tuning parameter
$\alpha>0$, and the third term is the $L_{1}$-norm of the discretized vector
Laplace operator with a tuning parameter $\beta>0$. As discussed in Section
3.4, even though the transverse integrals cannot be measured directly, we can
still say that only a small number of them are non-zero. Therefore, we employ
this knowledge by formulating an $L_{1}$-type sparsity prior that promotes
such behaviour with the help of the transverse ray matrix $R^{\perp}$.

To alleviate the limited data problem, we utilize the weighted vector Laplace
operator. Loosely speaking, the vector Laplace operator [18] relates the local
field values to the average of the surrounding points and thus imposes
“connectivity” between the neighboring points. The vector Laplace is also
related to the current sources as

$\nabla^{2}\mathrm{e}=\nabla(\nabla\cdot\mathrm{e})-\nabla\times(\nabla\times\mathrm{e})=\nabla(\nabla\cdot\mathrm{e})=-\frac{1}{\sigma}\nabla(\nabla\cdot\mathrm{j}^{s}).$
(34)

Because of this and the assumed sparsity of the current sources, we use also
here the $L_{1}$-norm. Furthermore, because it is known that minimizing the
$L_{1}$-norm of the vector Laplace yields harmonic solutions that have their
maxima on the boundary [10], we also use weight factors in (33). The matrix
$W$ has the form $W=\mathrm{diag}(w)D$ where the weights $w\in\mathbb{R}^{2N}$
have been estimated similarly to those in [23] and $D$ is the discretized
vector Laplace operator as in [45].

##  4 Numerical methods

In this section, we describe how the numerical approximations of the line
integrals and ray matrices were carried out and how the approach was tested
with numerical experiments.

###  4.1 Numerical approximation of ray matrices

The domain $\Omega$ is divided into $N_{E}$ disjoint triangular elements,
$\Omega=\cup_{j=1}^{N_{E}}\Omega_{j}$ with $N$ nodes and the electric field is
expressed in a vector form as
$e=[\mathrm{e}_{x},\mathrm{e}_{y}]^{\mathrm{T}}\in\mathbb{R}^{2N}$ (as in
Section 3.5). We use straight lines as the integration paths. The same lines
are used for both the longitudinal and transverse integrals. For the $i$th
longitudinal measurement along the line $L_{i}$, we can write

$I^{\parallel}_{i}=\int_{L_{i}}\mathrm{e}(\mathrm{x})\cdot\mathrm{\hat{s}_{i}}~{}d\ell(\mathrm{x})=\sum_{j=1}^{N_{E}}\int_{\Delta
L_{ij}}\mathrm{e}(\mathrm{x})\cdot\mathrm{\hat{s}}_{i}~{}d\ell(\mathrm{x}),$
(35)

where $L_{i}=\sum_{j=1}^{N_{E}}\Delta L_{ij}$ gives the line segmentation and
$\Delta L_{ij}=L_{i}\bigcap\Omega_{j}$ is the $j$th segment of $L_{i}$ (that
is inside the element $\Omega_{j}$). Figure 3 illustrates these variables.
Note that $\Delta L_{ij}\neq\varnothing$ only if the line $L_{i}$ passes
through the element $\Omega_{j}$.

Figure 3: Left: The computational domain is discretized into triangular
elements, and the domain is traced with lines that connect pairs of points at
the boundary. One of these lines $L_{i}$ passes through the elements that are
colored with gray. Right: The line segment $\Delta L_{ij}$ is the intersection
of the line $L_{i}$ and the triangular element $\Omega_{j}$. The coordinates
of the points on this line segment between $A_{i}$ and $B_{i}$ can be
expressed as a linear combination of the coordinates of the corner points of
the triangle $\mathrm{x}_{j}^{(1)}$, $\mathrm{x}_{j}^{(2)}$ and
$\mathrm{x}_{j}^{(3)}$.

####  4.1.1 Line segment in an element

Let us first examine a line segment $\Delta L_{ij}\neq\varnothing$ that passes
through an element $\Omega_{j}$ that has corner points $\mathrm{x}_{j}^{(1)}$,
$\mathrm{x}_{j}^{(2)}$ and $\mathrm{x}_{j}^{(3)}$ as shown in Figure 3. A
point on a line segment $\Delta L_{ij}$ has the position vector
$\mathrm{p}\in\mathbb{R}^{2}$

$\mathrm{p}=\mathrm{x}^{A_{i}}_{j}+(\mathrm{x}_{j}^{B_{i}}-\mathrm{x}_{j}^{A_{i}})t=\mathrm{x}^{A_{i}}_{j}+\Delta\mathrm{x}_{ij}~{}t,$
(36)

where $\mathrm{x}^{A_{i}}_{j}$ and $\mathrm{x}^{B_{i}}_{j}$ are the
intersecting points between the line $L_{i}$ and the edges of element
$\Omega_{j}$ (see Figure 3),
$\Delta\mathrm{x}_{ij}=\mathrm{x}_{j}^{B_{i}}-\mathrm{x}_{j}^{A_{i}}$ and
$t~{}\in~{}[0,1]$. By changing variables, we obtain
$\mathrm{\hat{s}}_{i}~{}d\ell(\mathrm{x})=\Delta\mathrm{x}_{ij}dt$ and the
line integral in this element becomes

$I^{\parallel}_{ij}=\int_{0}^{1}\Delta\mathrm{x}_{ij}\cdot\mathrm{e}(\mathrm{p}(t))~{}dt=\|\Delta\mathrm{x}_{ij}\|\int_{0}^{1}\mathrm{\hat{s}}_{i}\cdot\mathrm{e}(\mathrm{p}(t))~{}dt,$
(37)

as $\Delta\mathrm{x}_{ij}=\|\Delta\mathrm{x}_{ij}\|\mathrm{\hat{s}}_{i}$ where
$\|.\|$ denotes the length of the vector.

The electric field values at the corner points of the element are
$\mathrm{e}^{(1)}_{j}$, $\mathrm{e}^{(2)}_{j}$ and $\mathrm{e}^{(3)}_{j}$. We
approximate the field inside the element using the linear interpolation

$\mathrm{e(\mathrm{p}}(t))=\mathrm{e}^{(1)}_{j}+\begin{bmatrix}\mathrm{e}^{(2)}_{j}-\mathrm{e}^{(1)}_{j}&\mathrm{e}^{(3)}_{j}-\mathrm{e}^{(1)}_{j}\end{bmatrix}\mathrm{d}=\mathrm{e}^{(1)}_{j}+K_{j}\mathrm{d},$
(38)

where
$K_{j}=[\mathrm{e}^{(2)}_{j}-\mathrm{e}^{(1)}_{j}~{}~{}\mathrm{e}^{(3)}_{j}-\mathrm{e}^{(1)}_{j}]\in\mathbb{R}^{2\times
2}$ and $\mathrm{d}\in\mathbb{R}^{2\times 1}$ are the interpolation
coefficients (or bary-centric coordinates) [12]. For the estimation of
$\mathrm{d}$, we employ an iso-parametric mapping in which the position vector
in the element and the field are represented by the same interpolation
polynomial [54]. In particular, a point $\mathrm{x}$ in $\Omega_{j}$ has
position vector

$\mathrm{x}=\mathrm{x}^{(1)}_{j}+\begin{bmatrix}\mathrm{x}^{(2)}_{j}-\mathrm{x}^{(1)}_{j}&\mathrm{x}^{(3)}_{j}-\mathrm{x}^{(1)}_{j}\end{bmatrix}\mathrm{d}=\mathrm{x}^{(1)}_{j}+J_{j}\mathrm{d}$
(39)

where $\mathrm{d}=[d_{1}\,\ d_{2}]^{\mathrm{T}}$ with $d_{1}\geq 0$,
$d_{2}\geq 0$ and $d_{1}+d_{2}\leq 1$ and
$J_{j}=[\mathrm{x}^{(2)}_{j}-\mathrm{x}^{(1)}_{j}~{}~{}\mathrm{x}^{(3)}_{j}-\mathrm{x}^{(1)}_{j}]\in\mathbb{R}^{2\times
2}$. Now $\mathrm{d}$ can by solved as

$\mathrm{d}={J}_{j}^{-1}(\mathrm{x}-\mathrm{x}^{(1)}_{j}).$ (40)

When we set $\mathrm{x}=\mathrm{p}$ and combine Equations (36) and (40), we
can estimate the interpolation coefficient on the line segment with respect to
$t$

$\mathrm{d}={J}_{j}^{-1}(\mathrm{x}^{A_{i}}_{j}+\Delta\mathrm{x}_{ij}~{}t-\mathrm{x}^{(1)}_{j}).$
(41)

Now, this results in

$\mathrm{e(\mathrm{p}}(t))=\mathrm{e}^{(1)}_{j}+K_{j}{J}_{j}^{-1}(\mathrm{x}^{A_{i}}_{j}+\Delta\mathrm{x}_{ij}~{}t-\mathrm{x}^{(1)}_{j}),$
(42)

and Equation (37) becomes

$\begin{split}I^{\parallel}_{ij}=&\|\Delta\mathrm{x}_{ij}\|\left(\mathrm{\hat{s}}_{i}\cdot\mathrm{e}_{j}^{(1)}+\int_{0}^{1}\mathrm{\hat{s}}_{i}\cdot
K_{j}{J}_{j}^{-1}(\mathrm{x}^{A_{i}}_{j}+\Delta\mathrm{x}_{ij}~{}t-\mathrm{x}^{(1)}_{j})~{}dt\right)\\\
=&\|\Delta\mathrm{x}_{ij}\|~{}\mathrm{\hat{s}}_{i}\cdot\left(\mathrm{e}_{j}^{(1)}+K_{j}J_{j}^{-1}(\frac{1}{2}\Delta\mathrm{x}_{ij}-\mathrm{x}^{(1)}_{j})\right).\\\
\end{split}$ (43)

Finally when we write $K_{j}$ explicitly and denote
$C_{ij}=J_{j}^{-1}(\frac{1}{2}\Delta\mathrm{x}_{ij}-\mathrm{x}^{(1)}_{j})=[c_{1},c_{2}]^{\mathrm{T}}\in\mathbb{R}^{2\times
1}$ we get

$I^{\parallel}_{ij}=\|\Delta\mathrm{x}_{ij}\|~{}\mathrm{\hat{s}}_{i}\cdot\left[(1-c_{1}-c_{2})\mathrm{e}_{j}^{(1)}+c_{1}\mathrm{e}_{j}^{(2)}+c_{2}\mathrm{e}_{j}^{(3)}\right].$
(44)

We use Equation (44) to create a procedure for constructing the ray matrices.

####  4.1.2 Construction of ray matrices

The ray matrices are used to operate to the electric field in order to obtain
the line integral measurements. We denote the relationship between the
longitudinal measurements and the ray matrix operator as follows

$I^{\parallel}=R^{\parallel}e=\left[R_{x}^{\parallel}~{}~{}R_{y}^{\parallel}\right]\left[\begin{array}[]{c}{e}_{x}\\\
{e}_{y}\end{array}\right],$ (45)

where
$I^{\parallel}=\left[I_{1},\cdots,I_{m}\right]^{\mathrm{T}}\in\mathbb{R}^{m}$
contains all the longitudinal line measurements,
$R^{\parallel}\in\mathbb{R}^{m\times 2N}$ is the longitudinal ray matrix
operator that is divided into $R_{x}^{\parallel}\in\mathbb{R}^{m\times N}$ and
$R_{y}^{\parallel}\in\mathbb{R}^{m\times N}$ that contain the contributions of
the integral coefficients of the $x$ and $y$ field components, respectively.
The procedure that was used to construct the longitudinal ray matrix is shown
in Table 1. A similar procedure can be carried out to construct the transverse
ray matrix $R^{\perp}$ by exchanging the vector $\mathrm{\hat{s}}_{i}$ with
$\mathrm{\hat{s}}_{i}^{\perp}$ in the first step.

Table 1: Procedure to construct the longitudinal ray matrix. for i=1:m, go
through all the integration lines
---
take line $L_{i}$ and determine the corresponding unit vector
$\hat{s}_{i}=(s_{x},s_{y})$ along the line.
for j=1:$N_{E}$, go through all the elements
take element $\Omega_{j}$ that has corner points $\mathrm{x}_{j}^{(1)}$,
$\mathrm{x}_{j}^{(2)}$ and $\mathrm{x}_{j}^{(3)}$, where
$j^{(1)}$, $j^{(2)}$ and $j^{(3)}$ are the corresponding node indices.
if $\Delta L_{ij}=L_{i}\cap\Omega_{j}=\varnothing$
go to the next element.
else
calculate $\Delta\mathrm{x}_{ij}$ and $\|\Delta\mathrm{x}_{ij}\|$.
calculate $J_{j}$ and $C_{ij}=[c_{1},c_{2}]^{\mathrm{T}}$.
update the following entries of the $i$th row of the ray matrix.
$R^{\parallel}(i,j^{(1)})=R^{\parallel}(i,j^{(1)})+\|\Delta\mathrm{x}_{ij}\|s_{x}(1-c_{1}-c_{2})$
$R^{\parallel}(i,j^{(2)})=R^{\parallel}(i,j^{(2)})+\|\Delta\mathrm{x}_{ij}\|s_{x}c_{1}$
$R^{\parallel}(i,j^{(3)})=R^{\parallel}(i,j^{(3)})+\|\Delta\mathrm{x}_{ij}\|s_{x}c_{2}$
$R^{\parallel}(i,j^{(1)}+N)=R^{\parallel}(i,j^{(1)}+N)+\|\Delta\mathrm{x}_{ij}\|s_{y}(1-c_{1}-c_{2})$
$R^{\parallel}(i,j^{(2)}+N)=R^{\parallel}(i,j^{(2)}+N)+\|\Delta\mathrm{x}_{ij}\|s_{y}c_{1}$
$R^{\parallel}(i,j^{(3)}+N)=R^{\parallel}(i,j^{(3)}+N)+\|\Delta\mathrm{x}_{ij}\|s_{y}c_{2}$
end
end
end

###  4.2 Numerical experiments

In the experiments, we study electric fields generated by dipole sources in a
bounded 2-dimensional circular domain with homogeneous conductivity. We note
that the same methods can also be applied to any other convex domain. The
domain contained $n=32$ equally spaced electrodes around the boundary. Lines
for the integration were formed by connecting all pairs of electrodes which
resulted in total of $m=496$ lines. Two computational meshes were used: the
finer one consisted of $\bar{N}_{E}=9721$ triangular elements linking
$\bar{N}=3045$ nodes, and the coarser one of $N_{E}=1408$ triangular elements
with $N=760$ nodes.

####  4.2.1 Simulated forward fields and integral data

For the estimation of the longitudinal integral data and the forward electric
field, the finer mesh was used. First, a single dipole source was selected and
the corresponding scalar electric potential distribution
$\bar{u}\in\mathbb{R}^{\bar{N}}$ was computed by solving the Poisson problem
(12a) with boundary conditions (12b) and (12c) using finite element method
(FEM) with linear nodal basis functions [56]. The dipole source function was
numerically approximated using the mathematical dipole model [49]. The
longitudinal integral observations were calculated by taking the differences
of the potential values at the electrodes that locate at the ends of the
integration lines. Gaussian white noise was added to the data using two
signal-to-noise ratios, 40 dB and 20 dB. The signal-to-noise ratio (SNR) is
given by

$\mathrm{SNR}=20\,\mathrm{log}_{10}\frac{\|I^{\parallel}\|_{2}}{\|\varepsilon\|_{2}}$
(46)

The forward field, given by $\mathrm{\bar{e}}=-\nabla\bar{u}$, was estimated
numerically by applying the linear gradient reconstruction approach [11]. The
forward field ${\bar{{e}}}\in\mathbb{R}^{2\bar{N}}$ was then projected to the
inverse mesh ${e}=P\bar{{e}}$ where $P\in\mathbb{R}^{2N\times 2\bar{N}}$ is a
linear reduction mapping operator in order to be able to compare the
reconstructed field with the correct one.

####  4.2.2 Electric field reconstructions

The electric fields were reconstructed using the coarse mesh. First, the
longitudinal and transverse ray matrices were constructed for the mesh by
using the procedure described in Section 4.1. Then, the non-linear
minimization problem (33) was solved using convex optimization techniques [7,
20, 19]. The tuning parameters were kept constant in all experiments
($\alpha=0.06$ and $\beta=0.016$ in our case). The $\alpha$ value was chosen
to be high (compared to $\beta$) in order to ensure the sparsity of the
transverse measurements. The $\beta$ value was used to scale the coefficients
of the discrete weighted Laplace operator (that is used to promote
connectivity between neighbouring nodes) to match the average distance between
the mesh nodes.

####  4.2.3 Reconstruction error metrics

For the evaluation of the reconstructions, we used two different measures.
First, the average magnitude ratio (MR) between the reconstructed and the
actual field was computed as follows

$\mathrm{MR}=\frac{1}{N}\sum_{i=1}^{N}\frac{\sqrt{\hat{e}_{ix}^{2}+\hat{e}_{iy}^{2}}}{\sqrt{{e}_{ix}^{2}+{e}_{iy}^{2}}}.$
(47)

The closer the MR is to one the better the fields match with respect to the
magnitude.

Second, the average cosine similarity (CS) was used to quantify the difference
between the directions of the reconstructed and the actual field

$\mathrm{CS}=\frac{1}{N}\sum_{i=1}^{N}\mathrm{cos}(\mathrm{\hat{e}}_{i},\mathrm{e}_{i})=\frac{1}{N}\sum_{i=1}^{N}\frac{\hat{e}_{ix}e_{ix}+\hat{e}_{iy}e_{iy}}{\sqrt{\hat{e}_{ix}^{2}+\hat{e}_{iy}^{2}}\sqrt{{e}_{ix}^{2}+{e}_{iy}^{2}}}.$
(48)

CS has a value between –1 and +1, and the closer the value is to one the
better the directions of the fields match. Moreover, CS is close to zero when
the reconstructed and the actual field are perpendicular to each other, and
finally CS is close to –1 when the directions of the fields are opposite.

##  5 Results and discussion

Figures 4–6 show the reconstruction results of the electric fields produced by
a single dipole source which has radial orientation (Figure 4), tangential
orientation (Figure 5) and is located in the centre of the domain (Figure 6).
In each figure, the first column shows the true magnitude of the electric
field that is generated by a dipole source and the corresponding unit-length
field lines with the location and orientation of the underlying dipole source
(circle). The last two columns show the average electric field magnitude and
field lines that are calculated over 10 reconstructions. The 10
reconstructions were computed using 10 different realizations of noisy
boundary data with SNR = 40 dB (second column) and SNR = 20 dB (third column).

Figure 4: First column: Top picture shows the magnitude of the electric field
that is generated by a radial dipole source. Bottom picture shows the
corresponding unit length field lines and the location and orientation of the
underlying dipole source. Second column: Top picture shows the average
electric field magnitude that is calculated over 10 reconstructions and the
circle shows the correct location of the underlying dipole source. The 10
reconstructions were computed using 10 different realizations of noisy
boundary data with SNR = 40 dB. Bottom picture shows the corresponding unit-
length electric field lines. Third column: Top and bottom pictures show
similarly the reconstructed electric field when SNR = 20 dB. The average
magnitude ratio (MR) and cosine similarity (CS) values are given under the
reconstructions. Figure 5: First column: Top picture shows the magnitude of
the electric field that is generated by a tangential dipole source. Bottom
picture shows the corresponding unit length electric field lines with the
underlying dipole source. The other columns show the electric field
reconstructions as explained in Figure 4. Figure 6: First column: Top picture
shows the magnitude of the electric field that is generated by a dipole source
in the centre of the domain. Bottom picture shows the corresponding unit
length electric field lines with the underlying dipole source. The other
columns show the electric field reconstructions as explained in Figure 4.

”From the upper row of the figures (showing the field magnitude), we see that
the locations where the reconstructed fields get their maximum value are very
close to the correct locations (of the source) in all the test cases: in fact,
the location is exactly the same in the test cases shown in Figures 4 and 5,
and the location differences in Figure 6 correspond to the distance of a
single node in the mesh.

The dipole may be viewed as a positive and an equivalent negative charge that
are separated by a vanishingly small distance. The field lines that are caused
by a dipole are depicted in the first column of each picture: the arrows point
outwards from the positive and inwards from the negative end of the dipole.
From the reconstructed field lines in the lower row of Figures 4 and 5, it can
be seen that the locations of the positive and the negative charge are found
correctly but that they are separated by a small non-zero distance. When the
dipole source is in the centre of the domain, Figure 6, the locations of the
positive and negative charge are not evident from the field lines.

The maximum values of the reconstructed fields are lower by an order of
magnitude when compared to the actual ones. However, in the low noise cases
(SNR = 40 dB) in Figures 4 and 5, the magnitude ratio values are still very
high which indicates that the magnitude errors are present only near the
dipole sources and elsewhere the magnitudes are reconstructed accordingly.
Similarly in the field orientations, there are differences merely close to the
dipole sources, and especially for the low noise cases the CS values are high
which indicates that the field orientations are correct in most parts of the
domain. These reconstruction errors close to the source were expected because
the dipole source causes a discontinuity in the field.

By comparing the test cases, we can observe that the magnitude and orientation
errors are larger when the dipole is deep in the domain than close to the
boundary. We also see that the reconstruction accuracy decreases with
increasing noise as expected.

##  6 Conclusions and future work

In this paper, we have shown that vector tomography can be utilized for the
reconstruction of electric fields generated by dipole sources in bounded
domains. To the best of our knowledge, this is the first time that VT has been
used for the recovery of a field under non-zero divergence conditions. This
type of problem is of great interest due to its many applications, especially
in the EEG source imaging. As a comparison to the widely used source imaging
methods, we argue that the VT framework could be beneficial because it does
not require an explicit mathematical for the sources.

We first showed that the longitudinal measurements are directly determined by
the electric potentials at the boundary and the transverse measurements are
related to the underlying current sources. We explained that even though the
transverse measurements cannot normally be physically measured, they still can
be utilized in the electric field inverse problem as a penalty term.

Our numerical test cases included reconstructions of non-zero divergence
electric fields generated by a focal source with varying direction and
location. We showed that the pattern of the electric field magnitude could be
reconstructed correctly using VT, even though there were errors near the
dipole source. Nevertheless, for example, the correct location of the source
activity could be determined based on the reconstructed field magnitudes.
Also, the reconstructed field lines follow similar trajectories to the real
ones with some deviations only near the dipole sources. These reconstruction
errors were expected because the field is discontinuous at the dipole source.

In the future, in order to improve the feasibility and to study other kinds of
non-zero divergence fields, we shall extend the proposed approach for
3-dimensional VT problems with non-homogeneous material properties.

## Acknowledgment

The authors would like to thank the John S. Latsis Public Benefit Foundation
and the Finnish Cultural Foundation for financial support.

## Appendix: Electric field decomposition and line integrals

According to Helmholtz decomposition [2], any vector field can be decomposed
into a sum of irrotational (curl-free) and solenoidal (divergence-free)
component. When the field is zero on the domain boundary, the irrotational and
solenoidal components are unique, and they can be recovered from the
transverse and longitudinal measurements, respectively [8, 47]. However, when
non-zero boundary conditions are present the decomposition is not unique
anymore. For this case, however, a unique decomposition can be found by adding
a harmonic component which is both irrotational and solenoidal into the sum
[8, 44].

Let’s consider a similar electric field
$\mathrm{e(x):\Omega\rightarrow\mathbb{R}^{2}}$ in a bounded domain as in
Section 3. We set $\mathrm{e}|_{\partial\Omega}\neq 0$ and
$\mathrm{e}\cdot\hat{n}|_{\partial\Omega}=0$ on the boundary. The non-zero
boundary conditions imply that the decomposition of the field is of the form

$\mathrm{e}=\mathrm{e}_{\mathrm{I}}+\mathrm{e}_{\mathrm{S}}+\mathrm{e}_{\mathrm{H}},$
(49)

where $\mathrm{e}_{\mathrm{I}}$, $\mathrm{e}_{\mathrm{S}}$ and
$\mathrm{e}_{\mathrm{H}}$ are the irrotational, solenoidal and harmonic
component, respectively. The following properties apply for the irrotational
and solenoidal components
$\nabla\cdot\mathrm{e}=\nabla\cdot\mathrm{e}_{\mathrm{I}}$,
$\nabla\times\mathrm{e}=\nabla\times\mathrm{e}_{\mathrm{S}}$ and the harmonic
component satisfies both $\nabla\cdot\mathrm{e}_{\mathrm{H}}=0$,
$\nabla\times\mathrm{e}_{\mathrm{H}}=0$ and in addition also the boundary
conditions.

Under the quasi-static approximation, the electric field is irrotational. It
follows that the irrotational component can be expressed using a scalar
potential $\mathrm{e}_{\mathrm{I}}=-\nabla q$ and the solenoidal component
vanishes $\mathrm{e}_{\mathrm{S}}=0$. Furthermore, we write the harmonic
component as $\mathrm{e}_{\mathrm{H}}=-\nabla r+\nabla\times q$ where
$\nabla\times=(\frac{\partial}{\partial y},-\frac{\partial}{\partial x})$
corresponds to the 2-dimensional curl-operator. Now, $\nabla r=0$ due to the
boundary conditions. Therefore, the decomposition gets the form

$\mathrm{e}(\mathrm{x})=-\nabla q(\mathrm{x})+\nabla\times p(\mathrm{x}).$
(50)

The Radon transform of $\mathrm{e(x)}$ can be expressed using the Radon
property (5) as

$\begin{split}\tilde{\mathrm{e}}=\mathcal{R}_{\Omega}\\{\mathrm{e}\\}(l,\mathrm{s}^{\perp})&=-\mathrm{\hat{s}}^{\perp}~{}\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}(q)+\mathrm{\hat{s}}~{}\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}(p),\end{split}$ (51)

where ${\hat{s}}^{\perp}$ is the unit normal vector perpendicular and
$\hat{s}$ along the line $l$. Now, the longitudinal (6) and transverse (7)
integral have the form

$\displaystyle I^{\parallel}$ $\displaystyle=$
$\displaystyle\mathrm{\hat{s}}\cdot\mathcal{R}_{\Omega}\\{\mathrm{e}\\}(l,\mathrm{s}^{\perp})=\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}(p)$ (52) $\displaystyle I^{\perp}$ $\displaystyle=$
$\displaystyle\mathrm{\hat{s}}^{\perp}\cdot\mathcal{R}_{\Omega}\\{\mathrm{e}\\}(l,\mathrm{s}^{\perp})=\frac{\partial}{\partial
l}\mathcal{R}_{\Omega}(q).$ (53)

From these measurements, the electric field components can be solved as

$\displaystyle p$ $\displaystyle=$
$\displaystyle\mathcal{R}^{-1}\\{\tilde{p}\\}=\frac{1}{4\pi}\mathcal{R}^{\\#}\mathcal{H}\frac{\partial}{\partial
l}\mathcal{R}\\{p\\}=\frac{1}{4\pi}\mathcal{R}^{\\#}\mathcal{H}I^{\parallel}$
(54) $\displaystyle q$ $\displaystyle=$
$\displaystyle\mathcal{R}^{-1}\\{\tilde{q}\\}=\frac{1}{4\pi}\mathcal{R}^{\\#}\mathcal{H}\frac{\partial}{\partial
l}\mathcal{R}\\{q\\}=\frac{1}{4\pi}\mathcal{R}^{\\#}\mathcal{H}I^{\perp}.$
(55)

Thus, both types of line integral measurements are needed for the full
recovery of an electric field in a bounded domain with non-zero boundary
conditions. $I^{\parallel}$ is associated solely with the harmonic component,
thus the boundary conditions, and $I^{\perp}$ is needed to recover the
irrotational component.

## References

  * [1] A. Antoniadis, J. Bigot, Poisson inverse problems, Ann. Stat. 34 (5) (2006) 2132–2158.
  * [2] G. B. Arfken, H. J. Weber, Mathematical Methods for Physicists, 6th ed., Elsevier Academic Press, 2005.
  * [3] A. Balandin, Y. Ono, The method of series expansion for 3-D vector tomography reconstruction, J. Comput. Phys. 202 (2005) 52–64.
  * [4] A. Balandin, Y. Ono, S. You, 3D vector tomography using vector spherical harmonics decomposition, Computers Math. Applic. 63 (2012) 1433–1441.
  * [5] A. Balandin, N. Preobrazhenskii, A. Sedel’nikov, Tomographic determination of particle distribution by velocities, J. Appl. Mech. Tech. Phys. 30 (6) (1989) 862–865.
  * [6] M. Bauer, S. Pursiainen, J. Vorwerk, H. Köstler, C. H. Wolters, Comparison study for Whitney (Raviart-Thomas)-type source models in finite-element-method-based EEG forward modeling, IEEE Trans. Biomed. Eng. 62 (2015) 2648–2656.
  * [7] S. P. Boyd, L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004\.
  * [8] H. Braun, A. Hauck, Tomographic reconstruction of vector fields, IEEE Trans. Signal Process. 39 (1991) 464–471.
  * [9] H. Buchner, G. Knoll, M. Fuchs, A. Rienacker, R. Beckmann, M. Wagner, J. Silny, J. Pesch, Inverse localization of electric dipole current sources in finite element models of the human head, Electroenceph. Clin. Neurophysiol. 102 (1997) 267–278.
  * [10] M. Burger, H. Dirks, J. Müller., Inverse Problems in Imaging, Rad. Ser. Comp. App., 2013.
  * [11] C. D. Correa, K.-L. Ma, A Comparison of Gradient Estimation Methods for Volume Rendering on Unstructured Meshes, IEEE Trans. Vis. Comput. Graph. 17 (2011) 305–319.
  * [12] H. S. M. Coxeter, Introduction to Geometry, 2nd ed., John Wiley and Sons, 1969.
  * [13] S. R. Deans, The Radon Transform and Some of Its Applications, Dover Publications, 2007.
  * [14] E. Derevtsov, V. Pickalov, Reconstruction of vector fields and their singularities from ray transforms, Num. Anal. Applicat. 4 (1) (2011) 21–35.
  * [15] L. Desbat, A. Wernsdorfer, Direct algebraic reconstruction and optimal sampling in vector field tomography, IEEE Trans. Sign. Process. 43 (8) (1995) 1798 –1808.
  * [16] N. Efremov, N. Poluektov, V. Kharchenko, Tomography of ion and atom velocities in plasmas, J. Quant. Spectr. Rad. Transfer 53 (6) (1995) 723–728.
  * [17] A. Giannakidis, L. Kotoulas, M. Petrou, Virtual sensors for 2d vector field tomography, J. Opt. Soc. Am. A 27 (6) (2010) 1331–1341.
  * [18] L. J. Grady, J. Polimeni, Discrete Calculus: Applied Analysis on Graphs for Computational Science, Springer-Verlag London, 2010.
  * [19] M. Grant, S. Boyd, Graph implementations for nonsmooth convex programs, in: V. Blondel, S. Boyd, H. Kimura (eds.), Recent Advances in Learning and Control, Lecture Notes in Control and Information Sciences, Springer-Verlag Limited, 2008, pp. 95–110.
  * [20] M. Grant, S. Boyd, CVX: Matlab software for disciplined convex programming, version 2.1, http://cvxr.com/cvx (Mar. 2014).
  * [21] R. Grech, T. Cassar, J. Muscat, K. Camilleri, S. Fabri, M. Zervakis, P. Xanthopoulos, V. Sakkalis, B. Vanrumste, Review on solving the inverse problem in EEG source analysis, J. Neuroeng. Rehabil. 5 (1) (2008) 25.
  * [22] M. Hämäläinen, R. Hari, R. J. Ilmoniemi, J. Knuutila, O. V. Lounasmaa, Magnetoencephalography – theory, instrumentation, and applications to noninvasive studies on the working human brain, Rev. Mod. Phys. 65 (2) (1993) 413–497.
  * [23] S. Haufe, V. V. Nikulin, A. Ziehe, K.-R. Müller, G. Nolte, Combining sparsity and rotational invariance in EEG/MEG source reconstruction., NeuroImage 42 (2008) 726–738.
  * [24] S. Helgason, The Radon transform (Progress in Mathematics), vol. 5, 2nd ed., Birkhauser, 1999.
  * [25] H. M. Hertz, Kerr effect tomography for nonintrusive spatially resolved measurements of asymmetric electric field distributions, Appl. Opt. 25 (6) (1986) 914–921.
  * [26] J. Howard, Vector tomography applications in plasma diagnostics, Plasma Phys. Controlled Fusion 38 (4) (1996) 489.
  * [27] G. Huiskamp, M. Vroeijenstijn, R. V. Dijk, G. Wieneke, A. C. van Huffelen, The need for correct realistic geometry in the inverse EEG problem, IEEE Transaction on Biomedical Engineering 46 (11) (1999) 1281–1287.
  * [28] I. Javanovic, Inverse Problems in Acoustic Tomography: Theory and Applications, Ph.D. thesis, EPFL (2008).
  * [29] Johnson, Measurements and reconstruction of three-dimensional fluid flow, Patent US4162630 (1979).
  * [30] C. R. Johnson, Computational and numerical methods for bioelectric field problems (1997) 1–81.
  * [31] C. R. Johnson, R. S. MacLeod, M. A. Matheson, Computational medicine: Bioelectric field problems, Computer 26 (10) (1993) 59–67.
  * [32] S. A. Johnson, J. F. Greenleaf, W. Samayoa, F. Duck, J. Sjostrand, Reconstruction of three-dimensional velocity fields and other parameters by acoustic ray tracing, in: 1975 Ultrasonics Symposium, 1975.
  * [33] I. Jovanovia, L. Sbaiz, M. Vetterli, Acoustic tomography for scalar and vector fields: Theory and application to temperature and wind estimation, J. Atmos. Ocean. Technol.
  * [34] A. F. Kip, Fundamentals of electricity and magnetism, McGraw-Hill Book Company, 1969\.
  * [35] A. Koulouri, M. Petrou, Vector field tomography: reconstruction of an irrotational field in the discrete domain, Proc. IASTED Signal Process. Patt. Recogn. Applic.
  * [36] M. I. Kramar, A feasibility study about the use of vector tomography for the reconstruction of the coronal magnetic field., Ph.D. thesis, University of Gottingen (2005).
  * [37] S. J. Lade, D. Paganin, M. J. Morgan, 3-D vector tomography of Doppler-transformed fields by filtered-backprojection, Opt. Commun. 253 (4-6) (2005) 382–391.
  * [38] L. Lovstakken, S. Bjaerum, D. Martens, H. Torp, Blood flow imaging - a new real-time, flow imaging technique, IEEE Trans. Ultrason. Ferroelectri. Freq. Control 53 (2) (2006) 289–299.
  * [39] F. Natterer, The mathematics of computerized tomography, Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2001.
  * [40] S. J. Norton, Tomographic reconstruction of 2D vector fields: Application to flow imaging, Geophysi. J. Int. 97 (1) (1988) 161 – 168.
  * [41] S. J. Norton, Unique tomographic reconstruction of vector fields using boundary data, IEEE Transactions on Image Processing 1 (1992) 406–412.
  * [42] P. L. Nunez, R. Srinivasan, Electric Fields of the Brain: The Neurophysics of EEG, Oxford University Press, 2006.
  * [43] S. Osher, R. Fedkiw, Level Set Methods and Dynamic Implicit Surfaces, No. 153 in Applied Mathematical Sciences, Springer-Verlag New York, 2003.
  * [44] N. F. Osman, J. L. Prince, Reconstruction of vector fields in bounded domain vector tomography, in: Int. Conf. Image Process., 1997.
  * [45] R. D. Pascual-Marqui, C. M. Michel, D. Lehmann, Low resolution electromagnetic tomography: a new method for localizing electrical activity in the brain, Int. J. Psychophysiol. 18 (1) (1994) 49–65.
  * [46] M. Petrou, A. Giannakidis, Full tomographic reconstruction of 2D vector fields using discrete integral data, The Computer J. 53 (9).
  * [47] J. L. Prince, Tomographic reconstruction of 3D vector fields using inner product probes, IEEE Trans. Image Process. 3 (2) (1994) 216–219.
  * [48] D. Rouseff, K. B. Winters, T. E. Ewart, Reconstruction of oceanic microstructure by tomography: a numerical feasibility study, J. Geophys. Res. (96:C5) (1991) 8823–8833.
  * [49] P. H. Schimpf, C. Ramon, J. Haueisen, Dipole models for the EEG and MEG, IEEE Trans. Biomed. Eng. 49 (5) (2002) 409–418.
  * [50] T. Schuster, 20 years of imaging in vector field tomography: A review, in: Y. Censor, M. Jiang, A. K. Louis (eds.), Mathematical Methods in Biomedical Imaging and Intensity-Modulated Radiation Therapy (IMRT), vol. 7 of Publications of the Scuola Normale Superiore, (CRM), Birkhäuser, 2008.
  * [51] T. Schuster, D. Theis, A. K. Louis, A reconstruction approach for imaging in 3D cone beam vector field tomography., J. Biomed. Imaging 2008 (2008) 1–17.
  * [52] A. Schwarz, Three dimensional reconstruction of temperature and velocity fields in a Furnace., Part. Part. Syst. Char. 12 (2) (1995) 75–80.
  * [53] G. Sparr, K. Stråhlen, Vector field tomography, an overview, Center for Mathematical Sciences, Lund Institute of Technology, Sweden, tech. rep. (1998).
  * [54] P. J. Vauhkonen, Second order and infinite elements in three dimensional electrical impedance tomography, Phil. Lic. thesis (1999).
  * [55] A. Wernsdorfer, Complete reconstruction of three–dimensional vector fields, Proc. ECAPT 93 Karlsruhe (1993) 132–135.
  * [56] C. H. Wolters, L. Grasedyck, W. Hackbusch, Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem, Inverse Probl. 20 (4) (2004) 1099–1116.