Computational framework for applying electrical impedance tomography to head imaging

TL;DR

This paper presents a computational framework for head imaging using electrical impedance tomography that accurately estimates head shape and electrode positions without prior geometric information, enabling detection of internal conductivity variations.

Contribution

It introduces a principal component model for head shapes and integrates shape and electrode estimation into a regularized Newton reconstruction algorithm.

Findings

Can detect internal conductivity variations despite incomplete geometric info

Uses a library of fifty head shapes for modeling

Demonstrates robustness in numerical experiments

Abstract

This work introduces a computational framework for applying absolute electrical impedance tomography to head imaging without accurate information on the head shape or the electrode positions. A library of fifty heads is employed to build a principal component model for the typical variations in the shape of the human head, which leads to a relatively accurate parametrization for head shapes with only a few free parameters. The estimation of these shape parameters and the electrode positions is incorporated in a regularized Newton-type output least squares reconstruction algorithm. The presented numerical experiments demonstrate that strong enough variations in the internal conductivity of a human head can be detected by absolute electrical impedance tomography even if the geometric information on the measurement configuration is incomplete to an extent that is to be expected in practice.

Figures

Equations

$$\sigma\in L_{+}^{\infty}(\Omega):=\big{\{}\varsigma\in L^{\infty}(\Omega)\ |\ {\rm Re}(\varsigma)>c\ \,{\rm a.e.}\ {\rm in}\ \Omega\ {\rm for}\ {\rm some}\ c>0\big{\}}.$$

$${\rm Re}(\zeta)\geq 0,\qquad\zeta_{\partial\Omega\setminus{\overline{E}}}\equiv 0,\qquad{\rm Re}\big{(}\zeta|_{E_{m}}\big{)}\not\equiv 0$$

$$\mathbb{C}^{M}_{\diamond}\,:=\,\Big{\{}J\in\mathbb{C}^{M}\,\Big{|}\,\sum_{m=1}^{M}J_{m}=0\Big{\}}$$

$$U=\sum _{m=1}^M{U}_m{\chi }_m\in L^2(\partial \Omega ),$$

$$\begin{array}[]{ll}\displaystyle{\nabla\cdot(\sigma\nabla u)=0}&{\rm in}\;\;\Omega,\\[6.0pt] {\displaystyle{\nu\cdot\sigma\nabla u}=\zeta(U-u)}&{\rm on}\;\;\partial\Omega,\\[2.0pt] {\displaystyle\int_{E_{m}}\nu\cdot\sigma\nabla u\,{\rm d}S}=I_{m},&m=1,\ldots,M,\\[4.0pt] \end{array}$$

$${\mathcal{H}}^s:=H^s(\Omega )\oplus {\mathbb{C}}_{\diamond }^M,s\in \mathbb{R}.$$

$$\Vert (u,U){\Vert }_{{\mathcal{H}}^1}\le C\Vert I{\Vert }_2,$$

$$\Vert (u,U){\Vert }_{{\mathcal{H}}^s}\le C\Vert I{\Vert }_2\text{for}\text{any}s\in \mathbb{R},$$

$$S_{j}:\left\{\begin{array}[]{l}\mathbb{S}_{+}\to\mathbb{R}^{3},\\[2.84526pt] \hat{x}\mapsto r_{j}(\hat{x})\,\hat{x},\end{array}\right.$$

$$\mathbb{S}_{+}=\big{\{}x\in\mathbb{R}^{3}\;|\;\|x\|_{2}=1\ {\rm and}\ x_{3}>0\big{\}},$$

$$\bar{r}=\frac{1}{n}\sum _{j=1}^n{r}_j\mathrm{and}{\rho }_j=r_j-\bar{r},j=1,\dots ,n,$$

$$\sum _{j=1}^n\underset{\eta \in V_{\tilde{n}}}{\min }\Vert {\rho }_j-\eta {\Vert }_{H^1({\mathbb{S}}_{+})}^2\le \sum _{j=1}^n\underset{\eta \in W}{\min }\Vert {\rho }_j-\eta {\Vert }_{H^1({\mathbb{S}}_{+})}^2,$$

$$R_{ij}=({\rho }_i,{\rho }_j{)}_{H^1({\mathbb{S}}_{+})},i,j=1,\dots ,n,$$

$$V_{\tilde{n}}=\mathrm{span}\{{\hat{\rho }}_1,\dots ,{\hat{\rho }}_{\tilde{n}}\},$$

$$w_k=\frac{1}{\sqrt{{\lambda }_k}}{v}_k,k=1,\dots ,\tilde{n},$$

$$\sum _{j=1}^n\underset{\eta \in V_{\tilde{n}}}{\min }\Vert {\rho }_j-\eta {\Vert }_{H^1({\mathbb{S}}_{+})}^2$$

$$=\sum _{j=1}^n\Vert {\rho }_j{\Vert }_{H^1({\mathbb{S}}_{+})}^2-\sum _{j=1}^n\sum _{k=1}^{\tilde{n}}({\rho }_j,{\hat{\rho }}_k{)}_{H^1({\mathbb{S}}_{+})}^2.$$

$$\sum _{j=1}^n\sum _{k=1}^{\tilde{n}}({\rho }_j,{\hat{\rho }}_k{)}_{H^1({\mathbb{S}}_{+})}^2$$

$${\hat{\rho }}_k=w_k^{\mathrm{T}}\rho ,k=1,\dots ,\tilde{n},$$

$$({\hat{\rho }}_l,{\hat{\rho }}_k{)}_{H^1({\mathbb{S}}_{+})}=w_l^{\mathrm{T}}R{w}_k={\delta }_{lk},l,k=1,\dots ,\tilde{n},$$

$$\sum _{j=1}^n\sum _{k=1}^{\tilde{n}}({\rho }_j,{\hat{\rho }}_k{)}_{H^1({\mathbb{S}}_{+})}^2$$

$$=\sum _{k=1}^{\tilde{n}}\Vert R{w}_k{\Vert }_2^2=\sum _{k=1}^{\tilde{n}}{w}_k^{\mathrm{T}}{R}^2{w}_k$$

$$\sum _{j=1}^n\Vert {\rho }_j{\Vert }_{H^1({\mathbb{S}}_{+})}^2=\mathrm{tr}(R)=\sum _{j=1}^n{\lambda }_j.$$

$$\sum _{j=1}^n\Vert {\rho }_j{\Vert }_{H^1({\mathbb{S}}_{+})}^2-\sum _{j=1}^n\sum _{k=1}^{\tilde{n}}({\rho }_j,{\hat{\rho }}_k{)}_{H^1({\mathbb{S}}_{+})}^2=\sum _{j=1}^n{\lambda }_j-\sum _{k=1}^{\tilde{n}}{w}_k^{\mathrm{T}}{R}^2{w}_k=\sum _{j=\tilde{n}+1}^n{\lambda }_j$$

$$e(\tilde{n}):=\frac{{\sum }_{j=1}^n{\min }_{\rho \in V_{\tilde{n}}}\Vert {\rho }_j-\rho {\Vert }_{H^1({\mathbb{S}}_{+})}^2}{{\sum }_{j=1}^n\Vert {\rho }_j{\Vert }_{H^1({\mathbb{S}}_{+})}^2}=\frac{{\sum }_{j=\tilde{n}+1}^n{\lambda }_j}{{\sum }_{j=1}^n{\lambda }_j}.$$

$$S(\hat{x};\alpha)=\Big{(}\bar{r}(\hat{x})+\sum_{k=1}^{\tilde{n}}\alpha_{k}\hat{\rho}_{k}(\hat{x})\Big{)}\hat{x},\qquad\hat{x}\in\mathbb{S}_{+},$$

$$\alpha^{(j)}:=\big{[}(\rho_{j},\hat{\rho}_{1})_{H^{1}(\mathbb{S}_{+})},\dots,(\rho_{j},\hat{\rho}_{n})_{H^{1}(\mathbb{S}_{+})}\big{]}^{\rm T},\qquad j=1,\dots,n,$$

$$({\Gamma }_{\alpha }{)}_{kl}=\frac{1}{n-1}\sum _{j=1}^n({\rho }_j,{\hat{\rho }}_k{)}_{H^1({\mathbb{S}}_{+})}({\rho }_j,{\hat{\rho }}_l{)}_{H^1({\mathbb{S}}_{+})}=\frac{1}{n-1}{w}_k^{\mathrm{T}}{R}^2{w}_l=\frac{{\lambda }_k}{n-1}{\delta }_{kl},$$

$$E_{m}:=\big{\{}x\in\partial\Omega\cap B_{R_{0}}(x_{m})\ |\ \|Q_{x_{m}}x\|_{2}<R_{m}\big{\}},\qquad m=1,\dots,M,$$

$$\hat{E}_{m}:=\big{\{}x\in\partial\Omega\ |\ \|Q_{x_{m}}x\|_{2}<\mu R_{m}\big{\}}$$

$$U:(\sigma ,\zeta ,\alpha ,\beta ;I)\mapsto U(\sigma ,\zeta ,\alpha ,\beta ;I),$$

Full text

1 Introduction

Electrical impedance tomography (EIT) is a noninvasive imaging method that is based on current and voltage measurements on the boundary of the examined physical body. To be more precise, a set of contact electrodes is employed to drive current patterns into the object and the resulting electric potential is measured at (some of) the electrodes. Such measurements obviously depend on the interior conductivity of the object, and the aim of any algorithm designed for EIT is to (partially) invert this dependence and reconstruct (useful information about) the conductivity. The inverse problem of EIT is nonlinear and highly illposed; we refer to [3, 5, 32] for general information about EIT and its applications. The aim of this work is to develop a computational framework for applying EIT to head imaging without accurate information on the electrode positions or the shape of the studied patient’s head.

In medical applications of EIT, the exact shape of the imaged object is usually unknown. Even if one has some generic information about the average shape of the imaged part of the human body, the natural variations between different subjects are often significant: consider, e.g., the head shapes and sizes over the whole human population. Unfortunately, EIT is known to be extremely sensitive to geometric mismodeling of the imaged object [2, 4, 23]. The most straightforward remedy to inaccurate geometric modeling is difference imaging [1]: if electrode measurements are performed at two separate times and the aim is to reconstruct the corresponding change in the conductivity, the modeling errors partially cancel out and it is possible to get reasonable reconstructions. However, waiting for a change in the conductivity distribution is not a plausible option in many applications such as detection of stroke.

There exist a few previous algorithms for handling geometric uncertainties in absolute EIT imaging. The first one was introduced in [21, 22] for two-dimensional EIT and, intuitively speaking, it is based on compensating for the mismodeled geometry by reconstructing a mildly anisotropic conductivity. In [27, 28], the so-called approximation error approach [20] was applied to EIT with geometric uncertainties. To put it short, the error caused by the mismodeling was included as an extra additive noise process in the measurement model, its statistics were estimated in advance based on heavy simulations and prior knowledge on all unknown parameters, and finally the actual inversion was performed within the Bayesian paradigm. Recently, [14] built a polynomial surrogate for the dependence of the boundary measurements of EIT on all unknowns, including the parametrized measurement geometry, and employed this surrogate in straightforward Tikhonov regularization in two spatial dimensions. Of the aforementioned methods, the ones introduced in [21, 22] and [14] have only been implemented in two spatial dimensions, and the approximation error approach employed in [27, 28] requires a vast and expressive teaching sample describing the possible types of geometric mismodeling.

In this work we extend the conceptually straightforward way of handling imprecisely known body shape and electrode positions introduced in [8, 9, 10]. The basic idea is to utilize the derivatives of the electrode potentials with respect to the exterior boundary shape and the electrode locations in a regularized Newton-type output least squares algorithm that simultaneously reconstructs all unknowns in the measurement setup. The main weakness of the algorithm in [9, 10] is arguably the need to employ relatively dense finite element (FE) meshes in order to overcome an instability in the computation of the shape derivatives. We facilitate this stability issue by adopting the so-called smoothened complete electrode model (CEM) from [16] as our forward model; consult [6, 31] for information about the original CEM. Moreover, unlike in [9, 10] where the imaged three-dimensional objects were homogeneous in one spatial dimension, our parametrization of the imaged body is truly three-dimensional — as is the human head. See also [18] where the derivatives introduced in [8] were employed in simultaneous reconstruction of conductivity and electrode positions in a realistic but accurately modeled head geometry.

To be able to implement the algorithm of [9, 10] for head imaging by EIT, one needs (i) an appropriate parametrization for the natural variations in head shapes and sizes over the human population and (ii) a robust and efficient way of meshing a given (parametrized) human head model to enable forward solution by a finite element method (FEM). To build the former, we employ a library of fifty human heads from [26] and form a principal component model for the associated geometric variations. It turns out that a few first principal components already give a reasonably accurate geometric representation. In particular, this means that we only need to estimate a few geometric parameters as a part of the reconstruction algorithm. On the other hand, the latter requirement is satisfied by introducing a customized workflow for mesh generation, compatible with attaching electrode patches anywhere on the surface of a parametrized head.

Our numerical experiments demonstrate that strong enough alterations in the internal conductivity of a human head can be reconstructed by EIT even if the electrodes are misplaced and the geometric model for the head is inaccurate to the extent that is to be expected taking into account the natural variations. In fact, it seems that the misplacement of electrodes is a worse hindrance for the application of absolute EIT to head imaging than a moderately mismodeled head shape (cf. [15]): in many cases, altogether ignoring the inaccuracies in the head model and only including the estimation of the electrode locations and contact conductances in the algorithm already leads to a tolerable reconstruction of the conductivity.

This text is organized as follows. Section 2 reviews the smoothened CEM. The parametrized head model and the associated algorithm for forming head meshes are introduced in Section 3. Section 4 considers the computation of the (shape) derivatives needed in our Newton-type reconstruction algorithm. The numerical experiments are documented in Section 5 and the conclusions drawn in Section 6.

2 Forward model

We employ as our forward model the smoothened CEM [16] that has favorable regularity properties compared to the traditional CEM [31]. Preliminary tests with water tank data have also indicated that the smoothened version models real-world measurements approximately as accurately as the traditional CEM [16]

Consider a bounded Lipschitz domain $\Omega\subset\mathbb{R}^{3}$ and assume there are $M\in\mathbb{N}\setminus\{1\}$ electrodes $E_{1},\dots,E_{M}$ attached to its boundary. The electrode patches are identified with the nonempty, connected and open surface patches that they cover. We assume the electrodes are well separated, i.e., $\overline{E}_{m}\cap\overline{E}_{l}=\emptyset$ if $m\not=l$, and denote $E=\cup E_{m}$. The admittivity $\sigma:\Omega\to\mathbb{C}$ describing the electric properties of $\Omega$ is assumed to be bounded and strictly positive in the real part, that is,

[TABLE]

In other words, it is assumed that there are no ideally resistive or conductive regions of nonzero measure inside $\Omega$. As a specialty of the smoothened CEM, the contact admittance between the electrodes and the imaged object $\Omega$ is interpreted as a single function $\zeta\in L^{\infty}(\partial\Omega)$ such that

[TABLE]

for all $m=1,\dots,M$ in the topology of $L^{\infty}(\partial\Omega)$. These conditions can be interpreted as follows: the contact conductance cannot be negative, the contact admittance vanishes away from the electrodes, and on each electrode there must be a region of nonzero measure where the contact conductance is positive, i.e., there must be a possibility for feeding current through any of the electrodes.

When a single measurement by an EIT device is performed, net currents $I_{m}\in\mathbb{C}$, $m=1,\dots,M$, are driven through the corresponding electrodes and (noisy versions of) the resulting constant electrode potentials $U_{m}\in\mathbb{C}$, $m=1,\dots,M$, are measured. As it is natural to assume there are no sinks or sources inside $\Omega$, the current pattern $I=[I_{1},\dots,I_{M}]^{\rm T}$ belongs to the mean-free subspace

[TABLE]

due to conservation of electric charge. The vector containing the electrode potentials, $U=[U_{1},\dots,U_{M}]^{\rm T}\in\mathbb{C}^{M}$, is isomorphically identified with

[TABLE]

where $\chi_{m}$ is the characteristic function of the respective electrode $E_{m}$ on $\partial\Omega$. In the following, an electrode potential pattern $U$ may be treated either as a piecewise constant function supported on $\overline{E}$ or as an element of $\mathbb{C}^{M}$; the correct interpretation should be clear from the context.

The forward problem corresponding to the smoothened CEM is as follows [16]: the electromagnetic potential $u$ inside $\Omega$ and the piecewise constant electrode potential $U$ weakly satisfy

[TABLE]

where $\nu\in L^{\infty}(\partial\Omega,\mathbb{R}^{3})$ denotes the exterior unit normal of $\partial\Omega$. If one chooses $\zeta$ to be constant on each electrode, defines $z_{m}:=(1/\zeta)|_{E_{m}}$, $m=1,\dots,M$, and explicitly accounts for the second condition of (2), the forward problem (4) reduces to that associated to the standard CEM (cf., e.g., [31]). Hence, the CEM can actually be considered a special case of its smoothened counterpart.

We look for the solution $(u,U)$ of (4) in $\mathcal{H}^{1}$ defined via

[TABLE]

Note that the use of $\mathbb{C}_{\diamond}^{M}$ in (5) corresponds to systematically choosing the ground level of potential so that the mean of the electrode potentials is zero. The following proposition summarizes the unique solvability of (4) and the regularity properties of the associate forward solution.

Proposition 1**.**

Under the assumptions (1) and (2), the problem (4) has a unique solution $(u,U)\in\mathcal{H}^{1}$ that depends boundedly on $I\in\mathbb{C}_{\diamond}^{M}$, that is,

[TABLE]

where $C=C(\Omega,E,\sigma,\zeta)>0$ is independent of $I$. If in addition $\sigma\in\mathcal{C}^{\infty}(\overline{\Omega},\mathbb{C}^{n\times n})$, $\zeta\in\mathcal{C}^{\infty}(\partial\Omega)$ and $\partial\Omega\in\mathcal{C}^{\infty}$, then the solution to (4) satisfies

[TABLE]

where $C=C(\Omega,E,\sigma,\zeta,s)>0$ is still independent of $I\in\mathbb{C}_{\diamond}^{M}$.

Proof.

The assertions follow by slightly modifying the proofs of [16, Theorem 2.2 & Theorem 2.5] that are based on the use of quotient spaces instead of forcing the electrode potential patterns to be mean-free. ∎

The latter part of Proposition 1 could actually be made more explicit: the solution to (4) depends boundedly on the current pattern in the topology of $\mathcal{H}^{s}$ for a given $s\geq 1$ assuming that $\sigma$, $\zeta$ and $\partial\Omega$ are in certain regularity classes dictated by $s$; see [16, Theorem 2.5] for more information. This is the main reason why we resort to the smoothened CEM. By modeling the contact conductance as a smooth function, one achieves higher Sobolev regularity for the forward solution and thereby faster convergence of the numerical shape derivatives needed in our Newton type reconstruction algorithm; see [16] for more information. For comparison, the highest attainable regularity for the forward solution of the traditional CEM is $(u,U)\in\mathcal{H}^{2-\epsilon}$, $\epsilon>0$, even if $\partial\Omega$ and $\sigma$ are infinitely smooth.

3 Parametrized head model

To build a realistic model for the variations in the shape and size of the human head, we exploit the library of $n=50$ heads from [26]; see the top row of Figure 1 for two examples. Each of these heads allows a star-shaped parametrization if the origin is chosen to be (approximately) the center of mass of its bottom face. To be more precise, we represent the crown of the $j$th head as the graph of a function

[TABLE]

where $\mathbb{S}_{+}$ is the upper unit hemisphere, i.e.,

[TABLE]

and $r_{j}:\mathbb{S}_{+}\to\mathbb{R}_{+}$ gives the distance from the origin to the surface of the $j$th head as a function of the direction $\hat{x}\in\mathbb{S}_{+}$. Because $r_{j}$ completely determines $S_{j}$, we somewhat misleadingly call $r_{1},\dots,r_{n}$ the functional representations for our library of heads. The bottom row of Figure 1 shows such representations for the two example heads on the top row as functions of the polar $\theta\in(0,\pi/2)$ and azimuthal $\phi\in[0,2\pi)$ angles.

3.1 Optimal parametrization

Our aim is to form an accurate parametrization for the variations in the family of head representations $r_{1},\dots,r_{n}$ using as few degrees of freedom as possible. To this end, we set

[TABLE]

i.e., $\bar{r}$ describes the average head and $\rho_{j},\dots,\rho_{n}$ are the corresponding perturbations that define the employed library of heads. We then look for an $\tilde{n}$-dimensional subspace $V_{\tilde{n}}\subset H^{1}(\mathbb{S}_{+})$, $1\leq\tilde{n}\leq n$, that satisfies

[TABLE]

for all $\tilde{n}$-dimensional subspaces $W$ of the Sobolev space $H^{1}(\mathbb{S}_{+})$. In other words, the aim is to find an $\tilde{n}$-dimensional subspace that on average contains the best approximations for the perturbations $\rho_{j},\dots,\rho_{n}$ if the quality of fit is measured by the squared norm of $H^{1}(\mathbb{S}_{+})$. Here, it is implicitly assumed that $\rho_{j},\dots,\rho_{n}\in H^{1}(\mathbb{S}_{+})$.

For simplicity, suppose the functions $\rho_{1},\dots,\rho_{n}$ are linearly independent. In the following lemma, we utilize the matrix $R=\mathbb{R}^{n\times n}$ defined componentwise as

[TABLE]

where $(\cdot,\cdot)_{H^{1}(\mathbb{S}_{+})}:H^{1}(\mathbb{S}_{+})\times H^{1}(\mathbb{S}_{+})\to\mathbb{R}$ is the real $H^{1}$ inner product on $\mathbb{S}_{+}$. It is straightforward to check that $R$ is symmetric and positive definite. The eigenvalues of $R$ are repeated according to their multiplicity and denoted by $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq\lambda_{n}>0$. A corresponding set of orthonormal eigenvectors is $v_{1},\dots,v_{n}\in\mathbb{R}^{n}$. We also adopt the notation $\rho=[\rho_{1},\dots,\rho_{n}]^{\rm T}:\mathbb{S}_{+}\to\mathbb{R}^{n}$.

Lemma 2**.**

Assume that $\rho_{1},\dots,\rho_{n}\in H^{1}(\mathbb{S}_{+})$ are linearly independent and $1\leq\tilde{n}\leq n$. An $\tilde{n}$-dimensional subspace $V_{\tilde{n}}\subset L^{2}(\mathbb{S}_{+})$ that satisfies (9) is given by

[TABLE]

where $\hat{\rho}_{k}:=w_{k}^{\rm T}\rho$, $k=1,\dots,\tilde{n}$, are mutually orthonormal functions in $H^{1}(\mathbb{S}_{+})$ with

[TABLE]

being suitably scaled versions of the first $\tilde{n}$ eigenvectors for $R\in\mathbb{R}^{n\times n}$.

Proof.

First of all, any $\tilde{n}$-dimensional subspace $V_{\tilde{n}}\subset H^{1}(\mathbb{S}_{+})$ can of course be given as a linear span of $\tilde{n}$ functions, say, $\hat{\rho}_{1},\dots,\hat{\rho}_{\tilde{n}}\in H^{1}(\mathbb{S}_{+})$ that are orthonormal in the inner product of $H^{1}(\mathbb{S}_{+})$. Hence, by the Hilbert projection theorem,

[TABLE]

Finding $V_{\tilde{n}}$ that satisfies $\eqref{eq:minim}$ is thus equivalent to introducing orthonormal functions $\hat{\rho}_{1},\dots,\hat{\rho}_{\tilde{n}}\in H^{1}(\mathbb{S}_{+})$ that maximize the sum of squared projections

[TABLE]

appearing as the second term on the right-hand side of (3.1).

It follows easily from (9) that an optimal $V_{\tilde{n}}$ must, in fact, be a subspace of ${\rm span}\{\rho_{1},\dots,\rho_{n}\}$: otherwise one could reduce the dimension of $V_{\tilde{n}}$ by cutting away its component orthogonal to ${\rm span}\{\rho_{1},\dots,\rho_{n}\}$ without affecting the value on the left-hand side of (9). Hence, the maximizing orthonormal functions for (12), if they exist, can be written as

[TABLE]

for some $w_{k}\in\mathbb{R}^{n}$, $k=1,\dots,\tilde{n}$. The previously imposed condition that $\hat{\rho}_{1},\dots,\hat{\rho}_{\tilde{n}}$ are orthonormal is equivalent to asking that

[TABLE]

i.e., the coefficient vectors $w_{1},\dots,w_{\tilde{n}}\in\mathbb{R}^{n}$ are required to be orthonormal in the inner product induced by the symmetric and positive definite matrix $R\in\mathbb{R}^{n\times n}$.

To finalize the proof, we algebraically manipulate the expression (12) into a matrix form:

[TABLE]

as $R$ is symmetric. It is straightforward to show that this expression is maximized under the constraint (13) if the vectors $w_{1},\dots,w_{\tilde{n}}\in\mathbb{R}^{n}$ are chosen according to (10): one can, e.g., set $z_{k}=R^{1/2}w_{k}$ in order to rewrite the maximization problem in a more standard form and subsequently deduce that $z_{k}=v_{k}$, $k=1,\dots,\tilde{n}$. This completes the proof. ∎

Notice that

[TABLE]

In consequence, it follows from (3.1) and (3.1) that the left-hand side of (9) takes the value

[TABLE]

for an optimal $\tilde{n}$-dimensional subspace $V_{\tilde{n}}$. In particular, one has a simple formula for the relative squared error in the optimal representations for the considered library of heads:

[TABLE]

Figure 2 shows $e(\tilde{n})$ as a function of $\tilde{n}$. In particular, the optimal subspace $V_{\tilde{n}}$ is able to explain $58$% of the (squared) geometric variations in the head library already for $\tilde{n}=5$.

Encouraged by Lemma 2 and Figure 2, we introduce a parametrization for the crown of the human head as

[TABLE]

where $\bar{r}$ is the mean representation from (8), $\hat{\rho}_{1},\dots,\hat{\rho}_{\tilde{n}}$ are the optimal orthonormal basis functions introduced in Lemma 2, and $\alpha\in\mathbb{R}^{\tilde{n}}$ contains the free shape parameters. In the numerical experiments of Section 5, the number of free parameters appearing in (16) is chosen to be $\tilde{n}=5$; the corresponding basis functions $\hat{\rho}_{k}$, $k=1,\dots,5$, are visualized together with the average head in Figure 3.

Remark 3**.**

It is easy to see that the minimizing $\tilde{n}$-dimensional subspace $V_{\tilde{n}}$ for (9) is unique if $\lambda_{\tilde{n}}>\lambda_{\tilde{n}+1}$. Moreover, Lemma 2 holds in exactly the same form even if $\rho_{1},\dots,\rho_{n}$ are linearly dependent as long as there are at least $\tilde{n}$ linearly independent functions among them. On the other hand, if $\tilde{n}$ is larger or equal to the number of linearly independent functions among $\rho_{1},\dots,\rho_{n}$, then the left-hand side of (9) obviously vanishes for any optimal $\tilde{n}$-dimensional subspace $V_{\tilde{n}}$. Finally, observe that one could as well have used any other Sobolev space $H^{s}(\mathbb{S}_{+})$, $s\in\mathbb{R}$, in the above analysis, which would have led to $s$-dependent optimal representations.

3.2 Interpretation as principal components

In this section, we continue to assume that $\rho_{1},\dots,\rho_{n}\in H^{1}(\mathbb{S}_{+})$ are linearly independent and consider the limiting case $\tilde{n}=n$. We interpret $\alpha\in\mathbb{R}^{n}$ appearing in the parametrization (16), with $\tilde{n}=n$, as a random variable and the projection vectors

[TABLE]

onto the optimal basis $\hat{\rho}_{1},\dots,\hat{\rho}_{n}$ as its realizations over our library of heads. This is a natural interpretation because substituting $\alpha=\alpha^{(j)}$ in (16) leads to an exact representation for the $j$th head shape, i.e., for the $j$th realization of the human head in our sample (cf. (7) and Remark 3).

It follows trivially from (8) that the sample mean of $\{\alpha^{(1)},\dots,\alpha^{(1)}\}\subset\mathbb{R}^{n}$ is zero. Moreover, the associated sample covariance matrix $\Gamma_{\alpha}\in\mathbb{R}^{n\times n}$ can be given componentwise as

[TABLE]

where we used the same logic as in (3.1) together with the orthogonality of the scaled eigenvectors of $R$ defined by (10). In light of the available sample, the random variables $\alpha_{1},\dots,\alpha_{n}$ are thus uncorrelated and the corresponding variances form a non-increasing sequence as do the eigenvalues of $R$. Since $\alpha_{1},\dots,\alpha_{n}$ appear as the multipliers of the basis functions $\hat{\rho}_{1},\dots,\hat{\rho}_{n}$ in the parametrization (16) that is exact for the $j$th head if $\alpha$ takes its $j$th realization, the functions $\hat{\rho}_{1},\dots,\hat{\rho}_{n}$ can be interpreted as the principal components of the sample $\{\rho_{1},\dots,\rho_{n}\}$; see, e.g., [19] for more information on principal component analysis.

In the numerical experiments of Section 5, the shape coefficient vector $\alpha\in\mathbb{R}^{n}$ is assumed to be a priori distributed as $\alpha\sim\mathcal{N}(0,\Gamma_{\alpha})$ with the diagonal covariance matrix $\Gamma_{\alpha}\in\mathbb{R}^{n\times n}$ defined by (18). To visually examine the validity of such a hypothesis, Figure 4 depicts the histograms of the first six components of $\alpha\in\mathbb{R}^{n}$ over the sample (17). Due to the relatively small sample size, Figure 4 provides no conclusive evidence for the normality assumption, but it cannot be excluded either.

3.3 Building a FEM mesh with electrodes

We start by introducing an exact definition for a set of circular electrodes $E_{1},\dots,E_{M}\subset\partial\Omega$ with preassigned center points $x_{1},\dots,x_{M}\in\partial\Omega$ and radii $R_{1},\dots,R_{M}$. To this end, assume that the portion of $\partial\Omega$ to which the electrodes are attached is of class $\mathcal{C}^{1}$. Let $Q_{x}:\mathbb{R}^{3}\to\mathbb{R}^{2}$ be an orthogonal projection onto the tangent plane of $\partial\Omega$ at $x\in\partial\Omega$, with the tangent plane identified with $\mathbb{R}^{2}$ via choosing on it a local orthonormal coordinate system for which $x$ acts as the origin.111It is easy to check that the nonuniqueness in the choice of such a coordinate system does not affect the definition of a set of circular electrodes by (19). Let $B_{R}(x)\subset\mathbb{R}^{3}$ denote an open ball of radius $R>0$ around $x$. Using a compactness argument, it can be shown that the restriction $Q_{x}:B_{R}(x)\cap\partial\Omega\to\mathbb{R}^{2}$ is injective for all $x\in\partial\Omega$ if $R\leq R_{0}$ for some $R_{0}=R_{0}(\Omega)>0$; see [8] for closely related results. We then define

[TABLE]

which results in electrodes that have circular projections onto the tangent planes at their respective centers if $R_{1},\dots,R_{M}$ are small enough compared to $R_{0}$. In other words, the electrodes are modeled as elastic: although their projections are disks, the electrodes themselves elastically take the shape of the boundary they are attached to. This can be considered a reasonable approximation as long as the electrodes are small compared to the characteristic variations in the shape $\partial\Omega$.

Our workflow for generating a tetrahedral mesh for the head model consists of three steps: generation of an initial surface mesh, insertion of electrodes, and tetrahedral mesh generation. The initial surface mesh is constructed by subdividing $k$ times a coarse surface partition consisting of four triangles, where $k\in\mathbb{N}$ can be chosen by the operator of the algorithm. An electrode $E_{m}$, $m=1,\dots,M$, is inserted onto the surface mesh by finding the triangles that intersect the corresponding ‘extended electrode’

[TABLE]

with a suitably chosen $\mu>1$. The union of these (closed) triangles is denoted by $S_{m}$. The nodes of the surface mesh on the boundary $\partial S_{m}$ are then projected onto the tangent plane at $x_{m}$ by $Q_{x_{m}}$. A (dense) mesh $\mathcal{T}_{m}$ for the resulting polygonal domain is generated using the Triangle software [29]. To avoid hanging nodes, no new vertices are introduced on the boundary of the polygon at this stage. Subsequently, the surface mesh is modified by replacing the original triangles on $S_{m}$ with the (denser) triangulation $Q_{x_{m}}^{-1}\mathcal{T}_{m}$. After inserting all $M$ electrodes, the process is completed by generating a tetrahedral partition for the whole volume by TetGen [30] starting from the formed surface mesh.

4 Computation of derivatives

The electrode potentials $U\in\mathbb{R}_{\diamond}^{M}$ corresponding to a given electrode current pattern $I\in\mathbb{R}_{\diamond}^{M}$ can be interpreted as a function of five variables:

[TABLE]

where $\sigma\in L^{\infty}_{+}(\Omega)$ is the conductivity, $\zeta:\partial\Omega\to\mathbb{C}$ is the contact admittance satisfying (2), $\alpha\in\mathbb{R}^{\tilde{n}}$ is the parameter vector determining the shape of the computational head model in (16), and $\beta=(\beta_{1},\dots,\beta_{M})\in[\mathbb{S}_{+}]^{M}$ define the directions of the electrode centers from the origin of the parametrization (16). The aim of this section is to summarize how the partial derivatives of $U$ can be efficiently computed. To maintain readability, we usually suppress the dependence of $U$ on the variables that are not in the focus of our attention at a particular point in the text.

4.1 Derivative with respect to $\sigma$

It is well known that the electrode potential pattern $U$ is Fréchet differentiable with respect to $\sigma$; see, e.g., [12, 25]. To be more precise, the Fréchet derivative of $U$ at $\sigma\in L^{\infty}_{+}(\Omega)$ in the direction $\eta\in L^{\infty}(\Omega)$, i.e. $D_{\sigma}U(\sigma;\eta)\in\mathbb{C}_{\diamond}^{M}$, can be assembled using the formula (see, e.g., [17])

[TABLE]

where $(u,U)\in\mathcal{H}^{1}$ is the solution to (4) for the applied current pattern $I\in\mathbb{C}_{\diamond}^{M}$ and $(\tilde{u},\tilde{U})\in\mathcal{H}^{1}$ is the forward solution for an auxiliary current pattern $\tilde{I}$. We frequently employ similar ‘tilde notation’ in what follows.

Observe that knowing the projection $D_{\sigma}U(\sigma;\eta)\cdot\tilde{I}$ for all $\tilde{I}$ in a basis for $\mathbb{C}_{\diamond}^{M}$ uniquely determines $D_{\sigma}U(\sigma;\eta)\in\mathbb{C}_{\diamond}^{M}$, and since $M$ is typically well below hundred, this poses no computational challenges. On the other hand, if $M-1$ linearly independent current patterns are injected through the electrodes when taking the EIT measurements (as is the typical procedure), one must in any case solve the forward problem (4) with these same current patterns at each iteration of a Newton-type reconstruction algorithm. In consequence, one need not solve any extra forward problems to be able to perform the necessary evaluations on the right-hand side of (21). This means that the extra cost for computing the needed derivatives with respect to the conductivity corresponds to evaluating the integral on the right-hand side of (21) when $\eta$ runs through the degrees of freedom in the conductivity parametrization and $(u,U)$, $(\tilde{u},\tilde{U})$ over all pairs of the already computed forward solutions. The same conclusion also applies to the formulas (23) and (26) below.

4.2 Derivative with respect to $\zeta$

In order to be able to discuss differentiation with respect to $\zeta$, it is convenient to first describe a suitable parametrization for the considered contact admittances. To this end, we introduce a contact admittance (half) shape

[TABLE]

We then parametrize $\zeta:\partial\Omega\to\mathbb{C}$ as

[TABLE]

where the orthogonal projection $Q_{x_{m}}:\mathbb{R}^{3}\to\mathbb{R}^{2}$ is as defined in Section 3.3 and $\zeta_{m}\in\mathbb{C}$, with ${\rm Re}(\zeta_{m})>0$, defines the relative admittance on the $m$th electrode. Based on the material in [34], it is straightforward to show that the derivative of $U$ with respect to $\zeta_{m}$, $m=1,\dots,M$, can be assembled via (cf. [17, (26)])

[TABLE]

where we have used the same notation as in (21).

4.3 Derivative with respect to $\beta$

To begin with, let us consider a general way of perturbing the shapes and sizes of the electrodes; see [8, 14] for more information. Assume the portion of $\partial\Omega$ to which the electrodes are attached is of class $\mathcal{C}^{2}$. Let $a\in\mathcal{C}^{1}(E,\mathbb{R}^{3})$ and define a perturbed set of electrodes via

[TABLE]

where $P_{x}:\mathbb{R}^{3}\supset B_{R}(x)\to\partial\Omega$, for some $R>0$, is the projection in the direction of $\nu(x)$ onto $\partial\Omega$. The contact admittance $\zeta:\partial\Omega\to\mathbb{C}$ is assumed to stretch accordingly, that is, the admittance on the perturbed electrodes is defined via

[TABLE]

and it vanishes on $\partial\Omega\setminus\overline{E}^{a}$. It can be shown that there exists a constant $c=c(E,\Omega)>0$ such that the definitions (24) and (25) are unambiguous if $\|a\|_{\mathcal{C}^{1}(E,\mathbb{R}^{3})}<c$; see [8] for the proof of a closely related result.

The electrode potential pattern $U=U(a)\in\mathbb{C}_{\diamond}^{M}$ can obviously be interpreted as a function of (a small enough) $a\in\mathcal{C}^{1}(E,\mathbb{R}^{3})$. Indeed, for a given $a$ one simply considers the perturbed electrodes $E_{1}^{a},\dots,E_{M}^{a}$ and contact admittance $\zeta^{a}$ in (4) in place of the original unperturbed ones, which naturally leads to an $a$-dependent forward solution $(u(a),U(a))\in\mathcal{H}^{1}$.

The mapping $\mathcal{C}^{1}(E,\mathbb{R}^{m})\ni a\mapsto U(a)\in\mathbb{R}_{\diamond}^{M}$ is Fréchet differentiable at the origin [8, 16]. The associated Fréchet derivative in the direction $a\in\mathcal{C}^{1}(E,\mathbb{R}^{m})$ can be assembled using

[TABLE]

where $a_{\tau}:E\to\mathbb{R}^{3}$ is the tangential component of $a\in\mathcal{C}^{1}(E,\mathbb{R}^{3})$ and ${\rm Grad}$ denotes the surface gradient [7]. For a piecewise constant $\zeta$, as in the traditional CEM, the right-hand side of (26) reduces to an integral over $\partial E$ [8].

In order to utilize (26) in a reconstruction algorithm, we still need to define vector fields that (approximately) correspond to changing the direction of, say, the $m$th electrode, i.e., the parameter $\beta_{m}\in\mathbb{S}_{+}$ in (20). To this end, we represent $\beta_{m}=\beta_{m}(\theta_{m},\phi_{m})$ as a function of its polar $\theta_{m}\in(0,\pi/2)$ and azimuthal $\phi_{m}\in[0,2\pi)$ angles over the upper hemisphere $\mathbb{S}_{+}$. We denote by $\hat{\theta}_{m}\in\mathbb{R}^{3}$ and $\hat{\phi}_{m}\in\mathbb{R}^{3}$ the tangent vectors for $\partial\Omega$ at the center point $x_{m}$ of $E_{m}$ obtained by (numerically) differentiating the parametrization (16) with respect to $\theta$ and $\phi$, respectively.

We extend $\hat{\theta}_{m}$ as a constant vector field over the whole electrode $E_{m}$ and define

[TABLE]

to be its tangential component. (Note that $\tilde{a}^{\theta}_{m}$ depends on the location on the electrode as does the tangent plane.) To ensure that all points on $E_{m}$ move as much in the tangential direction, i.e., that the size of the electrode is approximately maintained, we finally normalize

[TABLE]

and extend $a^{\theta}_{m}$ as zero to the rest of $E$. Although the perturbation of $E_{m}$ by such a vector field $a^{\theta}_{m}$ in the spirit of (24) does not exactly maintain the circular shape of the elastic electrode $E_{m}$, it seems feasible that the variations in its shape and size are typically of higher order compared to the moved distance — especially if the electrodes are small compared to the characteristic variations in the shape of $\partial\Omega$. A perturbation field for $E_{m}$ in the direction of $\hat{\phi}_{m}$ is introduced in the same manner.

To summarize, we slightly abuse the notation by identifying $U(\beta)=U(\theta,\phi)$, where $\theta=[\theta_{1},\dots,\theta_{M}]^{\rm T}$ and $\phi=[\phi_{1},\dots,\phi_{M}]^{\rm T}$, and approximate the corresponding derivatives of $U$ based on the formulas

[TABLE]

and

[TABLE]

for $m=1,\dots,M$.

Remark 4**.**

It would definitely be worthwhile to analyze how the vector fields $a_{m}^{\theta}:E\to\mathbb{R}^{3}$ and $a_{m}^{\phi}:E\to\mathbb{R}^{3}$ should be defined so that they would exactly correspond to differentiation with respect to the polar and azimuthal angles for the center point of the $m$th electrode. However, as the above introduced crude approximations seem to function satisfactorily in the reconstruction algorithm introduced in Section 5, we postpone such analysis to future studies. In any case, it should be emphasized that utilizing the formulas (27) and (28) is computationally far more attractive than resorting to difference approximations: Using, e.g., central differences to approximate the derivatives on the left-hand sides of (27) and (28) would require remeshing and computing forward solutions of (4) for $4M$ extra measurement geometries. In our numerical studies with $M=32$ electrodes, this would correspond to forming $128$ extra FEM meshes and solving the associated forward problems with $M-1=31$ different current patterns at each step of the reconstruction algorithm.

4.4 Derivative with respect to $\alpha$

In principle, one could try to employ the sampling formulas introduced in [9, 16] to compute the Fréchet derivative of $U=U(\alpha)$ with respect to the shape parameter vector $\alpha\in\mathbb{R}^{\tilde{n}}$ in (16). However, the analysis in [9, 16] is based on the assumption that the electrodes stretch accordingly when the shape of $\Omega$ changes. For cylindrical domains, as are the ones considered in [9, 10], compensating for such stretching is straightforward because an electrode is essentially determined by its two end points; see [9]. However, in our inherently three-dimensional setting, analyzing how the changes in the sizes and shapes of the electrodes affect the shape derivatives becomes more complicated, and we choose to leave such consideration for future studies.

As discussed in Section 3.1, we use $\tilde{n}=5$ in our numerical experiments, and thus approximating the derivatives with respect to such a low number of shape parameters can be done by resorting to, e.g., central differences without increasing the computational load out of proportion. That is, we numerically approximate

[TABLE]

where ${\rm e}_{j}$ is the $j$th Cartesian basis vector and $\epsilon>0$ is chosen appropriately.

Observe that (29) is in any case computationally significantly less attractive than (21), (23) and (26) since it requires extra solutions of the forward problem (4) as well as remeshing of the domain. To be more precise, the use of (29) demands solving the forward problem with all employed current patterns for $2\tilde{n}$ extra domain shapes. Hence, it is imperative to keep $\tilde{n}$ as low as possible in the reconstruction algorithm (cf. Remark 4).

5 Numerical experiments

This section presents our numerical examples. We start by briefly reviewing the reconstruction algorithm that is based on a simple Gauss–Newton iteration; for more information consult, e.g., [10, 15]. Next, we explain how realistic measurement data are simulated. Finally, some reconstructions are presented with emphasis on analyzing how much ignoring different kinds of geometric inaccuracies in the forward model deteriorate their quality.

Although all theory presented above applies to complex-valued admittivities and contact admittances as well as to the resulting complex-valued potential measurements, we assume from this point on that the frequency of the alternating current fed through the electrodes is so low that the phase information of the measurements can be ignored and the amplitudes of the net currents and electrode potentials can be treated as if they resulted from measurements with direct current; see [33] for detailed analysis on the validity of such an approximation. In particular, we assume that $\sigma$ and $\zeta$ are real-valued and dub them the conductivity and the contact conductance, respectively, in what follows.

5.1 The algorithm

We assume to be able to drive $M-1$ linearly independent current patterns $I^{(1)},\dots,I^{(M-1)}\in\mathbb{R}_{\diamond}^{M}$ through the electrodes and measure the corresponding noisy electrode potentials denoted by $V^{(1)},\dots,V^{(M-1)}\in\mathbb{R}^{M}$. The potential measurements are packed into a single vector

[TABLE]

and we analogously introduce the stacked forward map

[TABLE]

via

[TABLE]

Here, the conductivity $\sigma\in\mathbb{R}_{+}^{N}$ is identified with the $N\in\mathbb{N}$ degrees of freedom used to parametrize it, the contact conductance $\zeta\in\mathbb{R}^{M}$ is identified with the coefficient vector in (22), $\alpha\in\mathbb{R}^{\tilde{n}}$ is the parameter vector in (16) determining the shape of the computational head model, and $\theta\in(0,\pi/2)^{M}$ and $\phi\in[0,2\pi)^{M}$ define the polar and azimuthal angles of the electrode center points, respectively. As the shape of the computational domain changes during the iteration, the degrees of freedom $\sigma\in\mathbb{R}_{+}^{N}$ are chosen to be the coefficients of a piecewise linear parametrization for the conductivity with respect to a FE mesh in a large storage head that encloses all incarnations of the computational head model encountered by our iterative reconstruction algorithm. The conductivity parametrization is transferred between the current computational head and the storage head via linear projections between the corresponding FE meshes.

The reconstruction algorithm aims at computing a minimizer for the Tikhonov functional

[TABLE]

where the weighted norms are defined by $\|x\|_{A}^{2}:=x^{\rm T}\!Ax$. A minimizer of (5.1) can be interpreted as a maximum a posteriori (MAP) estimate for the to-be-reconstructed parameters, assuming the measurements are corrupted by additive Gaussian noise with zero mean and covariance $\Gamma_{\eta}\in\mathbb{R}^{M(M-1)\times M(M-1)}$ and the prior distributions for $\sigma,\zeta,\alpha,\theta$ and $\phi$ are also Gaussian with the mean-covariance pairs $(\bar{\sigma},\Gamma_{\sigma})$, $(\bar{\zeta},\Gamma_{\zeta})$, $(\bar{\alpha},\Gamma_{\alpha})$, $(\bar{\theta},\Gamma_{\theta})$ and $(\bar{\phi},\Gamma_{\phi})$, respectively. The Cholesky factors of the above introduced positive definite covariance matrices are denoted by $L_{*}$, where $*$ stands for the variable in question. Moreover,

[TABLE]

denotes the Jacobian matrix of $\mathcal{U}(\sigma,\zeta,\alpha,\theta,\phi)$ that can be numerically approximated by following the guidelines in Section 4.

Algorithm 1**.**

Assume the covariance matrices $\Gamma_{\eta},\Gamma_{\sigma},\Gamma_{\zeta},\Gamma_{\alpha}$, $\Gamma_{\theta}$, $\Gamma_{\phi}$ and the expected values $\bar{\alpha}$, $\bar{\theta}$, $\bar{\phi}$ are given. Compute Cholesky factorizations for the positive definite covariance matrices and form a block-diagonal matrix

[TABLE]

Choose the prior means for the conductivity and contact conductances to be the homogeneous estimates $\bar{\sigma}=\tau_{\sigma}\mathbf{1}$ and $\bar{\zeta}=\tau_{\zeta}\mathbf{1}$ where

[TABLE]

and ${\bf 1}=[1,\ldots,1]^{\rm T}$ is at each occurrence a constant vector of the appropriate length.

Initialize

[TABLE]

and iterate as follows starting from $j=0$:

Form

[TABLE] 2. 2.

Solve the direction $\Delta b$ from

[TABLE] 3. 3.

Set $\tilde{b}^{(j+1)}=b^{(j)}-q\Delta b$ where $0<q\leq 1$ is fixed by the operator of the algorithm. 4. 4.

1. (a)

If $F(\tilde{b}^{(j+1)})<F(b^{(j)})$, set $b^{(j+1)}=\tilde{b}^{(j+1)}$, update $j=j+1$ and return to step 1. 2. (b)

If $F(\tilde{b}^{(j+1)})\geq F(b^{(j)})$, perform a brief line search on $[b^{(j)},\tilde{b}^{(j+1)}]$. If decrease in the value of $F$ is observed, dub the corresponding minimizer $b^{(j+1)}$ and return to step 1; otherwise, terminate the iteration.

Define $[\sigma_{*}^{\rm T},\zeta_{*}^{\rm T},\alpha_{*}^{\rm T},\theta_{*}^{\rm T},\phi_{*}^{\rm T}]^{\rm T}:=b^{(j)}$ to be the reconstruction.

We choose $q=0.5$ in all numerical examples documented in Section 5.3, which seems to provide a suitable trade-off between robustness and speed of convergence for the considered setups.

5.2 Simulation of measurement data

We consider a couple of different target heads. In each case, the parameters $\alpha_{\rm trgt}\in\mathbb{R}^{n_{0}}$ defining the shape of the target are drawn from (a variant of) the distribution $\mathcal{N}(0,\Gamma_{\alpha})$, where, allowing a slight abuse of notation, $\Gamma_{\alpha}\in\mathbb{R}^{n_{0}\times n_{0}}$ is the diagonal covariance matrix defined componentwise by (18); for a motivation of this choice, see the principal component construction in Section 3.2. The number of components in $\alpha_{\rm trgt}$ is chosen considerably higher than the number of shape parameters reconstructed by Algorithm 1, i.e., $n_{0}=10$ vs. $\tilde{n}=5$. In other words, the heads used for simulating data are more detailed than the reconstruction algorithm is capable to accurately represent.

The expected values for the polar and azimuthal angles of the electrode centers, i.e. $\bar{\theta}$ and $\bar{\phi}$, are chosen to correspond to the angular positions to which one originally aimed to attach the electrodes. The actual central angles of the target electrodes, i.e. $\theta_{\rm trgt}$ and $\phi_{\rm trgt}$, are then drawn from the distributions $\mathcal{N}(\bar{\theta},\Gamma_{\theta})$ and $\mathcal{N}(\bar{\phi},\Gamma_{\phi})$, where

[TABLE]

Here, $\mathbb{I}\in\mathbb{R}^{M\times M}$ is the identity matrix and $\varsigma_{\theta},\varsigma_{\phi}>0$ determine the standard deviations in the two angular directions. Notice that $\varsigma_{\theta}$ and $\varsigma_{\phi}$ must be chosen so small that the electrodes are not at a risk to overlap or move outside the crown of the computational head. The radii of all electrodes are chosen to be the same, say, $R>0$, which is the most common setup encountered in practice.

The relative contact conductances $(\zeta_{\rm trgt})_{m}\in\mathbb{R}_{+}$, $m=1,\dots,M$, appearing in (22) are independently drawn from $\mathcal{N}(\tilde{\tau}_{\zeta},\tilde{\varsigma}_{\zeta}^{2})$, where $\tilde{\tau}_{\zeta}>0$ is chosen so much larger than $\tilde{\varsigma}_{\zeta}>0$ that negative contact conductances are practically impossible. The corresponding conductance shape is chosen to be $\hat{\zeta}_{\rm trgt}\equiv 1$ in (22), that is, we employ the traditional CEM with constant contact conductances when simulating data. Finally, $\sigma_{\rm trgt}\in\mathbb{R}_{+}^{N_{0}}$, $N_{0}\in\mathbb{N}$, corresponds to a piecewise linear parametrization for the studied target conductivity on a dense FE mesh associated to the target head and the electrode positions defined by $\alpha_{\rm trgt}$, $\theta_{\rm trgt}$ and $\phi_{\rm trgt}$.

We approximate the ideal data $\mathcal{U}(\sigma_{\rm trgt},\zeta_{\rm trgt},\alpha_{\rm trgt},\theta_{\rm trgt},\phi_{\rm trgt})$ by FEM with piecewise linear basis functions and denote the resulting (almost) noiseless data by $\mathcal{U}_{\rm trgt}\in\mathbb{R}^{M(M-1)}$. The actual noisy data is then formed as

[TABLE]

where $\eta\in\mathbb{R}^{M(M-1)}$ is a realization of a zero-mean Gaussian with the diagonal covariance matrix

[TABLE]

with $\mathbb{I}\in\mathbb{R}^{M(M-1)\times M(M-1)}$. The free parameter $\varsigma_{\eta}>0$ can be tuned to set the relative noise level. Such a noise model has been used with real-world data, e.g., in [10].

5.3 Reconstructions

Let us then focus on the parameter choices for Algorithm 1. In (5.1) we employ $\Gamma_{\eta}\in\mathbb{R}^{M(M-1)\times M(M-1)}$, $\Gamma_{\theta}\in\mathbb{R}^{M\times M}$ and $\Gamma_{\phi}\in\mathbb{R}^{M\times M}$ as defined in (33) and (32), respectively. In other words, we assume to know an accurate model for the measurement noise as well as to have a good understanding of the ways the electrodes are typically misplaced. Moreover, it is not far-fetched to assume the intended angular positions for the electrodes, i.e. $\bar{\theta}$ and $\bar{\phi}$, are also known. The matrix $\Gamma_{\alpha}\in\mathbb{R}^{\tilde{n}\times\tilde{n}}$, with $\tilde{n}=5$, in (5.1) is defined by (18) and the corresponding expected value $\bar{\alpha}$ is set to zero; see Section 3.2.

In the reconstruction algorithm, we choose the form

[TABLE]

for the contact shape $\hat{\zeta}:[0,1)\to\mathbb{R}_{+}$ appearing in (22). This leads to a smooth contact conductance $\zeta:\partial\Omega\to\mathbb{R}_{+}\cup\{0\}$ with $\zeta(x_{m})=1$ for each electrode center $x_{m}$, $m=1,\dots,M$.222Replacing $2$ by any other positive real number in (34) would also lead to a smooth contact conductance, but this particular choice seems to function well in our numerical examples. The ‘expected value’ $\bar{\zeta}$ in (5.1) is selected automatically by Algorithm 1. We choose the corresponding prior covariance to be $\Gamma_{\zeta}=(2\tilde{\varsigma}_{\zeta}\bar{\zeta}/\tilde{\tau}_{\zeta})^{2}\mathbb{I}$, where $\tilde{\tau}_{\zeta}$ and $\tilde{\varsigma}$ are the expected value and standard deviation employed when drawing the target contact conductance values. The reason for this ad hoc choice is two-fold: First of all, as the contact conductance models for the data simulation and the reconstruction process are different, there is no reason to expect that the distributions of the associated parameters are the same; for further analysis, see [16], where it is noted that the smoothened model often leads to significantly higher peak values for the contact conductance function. Secondly, the estimation of the contact conductances is usually a rather stable process which allows one to employ only mild regularization, i.e., twice as large ratio between the prior standard deviation and expected value than used in the data simulation step.

The prior covariance matrix for the conductivity is defined componentwise as

[TABLE]

where $\varsigma_{\sigma}>0$ is the pointwise standard deviation, $l>0$ is the correlation length, and $x_{i},x_{j}\subset\mathbb{R}^{2}$ are coordinates of FE nodes in the storage head. Notice that there is no need to choose the ‘expected value’ $\bar{\sigma}$ in (5.1) as it is automatically selected by Algorithm 1.

The electrode potentials and associated Jacobian matrices needed in Algorithm 1 were approximated by FEM with piecewise linear basis functions on meshes with about 20 000 nodes and 85 000 tetrahedra; see Section 3.3 and, e.g., the right-hand image of Figure 5. The number of degrees of freedom for the conductivity parametrization on the storage head was about $N=9000$. All computations were performed with a nonoptimized MATLAB implementation on a laptop with 16 GB RAM and an Intel CPU having clock speed 3.5 GHz.

5.3.1 Test 1: one insulating and one conductive inclusion

The left-hand image of Figure 5 illustrates the first target head with the misplaced electrodes on its boundary. The azimuthal and polar angles of the electrodes were drawn as explained in Section 5.2 with $\varsigma_{\theta}=\varsigma_{\phi}=0.03$ radians in (32). The right-hand image of Figure 5 shows the geometric initial guess for Algorithm 1, i.e., the average head shape together with electrodes at their intended positions organized into three belts with constant polar angles. Counting upwards from the bottom belt, there are altogether $M=16+10+6=32$ electrodes and their common radius is $R=0.75$ cm. The mean and standard deviation used for drawing the relative contact conductances $\zeta_{\rm trgt}\in\mathbb{R}^{M}$ were $\tilde{\tau}_{\zeta}=100\,{\rm S}/{\rm m}^{2}$ and $\tilde{\varsigma}_{\zeta}=10\,{\rm S}/{\rm m}^{2}$, respectively. Observe that such contacts are comparable to those documented in [6, 10, 13, 14], if the quality of contact is measured by the ratio between the contact conductance and the conductivity level of the imaged body. The relative noise level in (33) was set to $\varsigma_{\eta}=10^{-3}$ (cf. [10]), and the pointwise standard deviation and correlation length for the prior in (35) were chosen to be $\varsigma_{\sigma}=0.1\,{\rm S}/{\rm m}$ and $l=3.3\,{\rm cm}$, respectively.

The target conductivity consists of two balls with constant conductivities $\sigma_{1}=2\,{\rm S}/{\rm m}$ and $\sigma_{2}=0.02\,{\rm S}/{\rm m}$ embedded in a homogeneous background characterized by $\sigma_{0}=0.2\,{\rm S}/{\rm m}$. Take note that such a background level is in line with conductivities reported for white and gray matter in the human brain [11, 24]. The center and radius of the conductive inclusion are $y_{1}=(-2,-2,6)$ cm and $r_{1}=1.5$ cm, respectively, whereas the insulating inclusion has radius $r_{2}=2$ cm and it is centered around $y_{2}=(3,3,3)$ cm.

Figure 6 compares the target conductivity with three reconstructions. The top right reconstruction corresponds to the complete Algorithm 1; the Gauss–Newton iteration converged in six rounds and about twenty minutes. The bottom left image resulted from not taking the mismodeling of the head shape into account, that is, the reconstruction was formed by setting $\alpha=0$ in Algorithm 1 and only estimating the other four unknowns. As the computation of the derivatives with respect to the shape parameters is the most time consuming part of Algorithm 1, forming this second reconstruction took only a few minutes, although the Gauss–Newton iteration lasted for one extra round. Finally, the bottom right reconstruction corresponds to altogether ignoring geometric inaccuracies and fixing $\alpha=0$, $\theta=\bar{\theta}$ and $\phi=\bar{\phi}$ in Algorithm 1. Computing this third reconstruction in seven Gauss–Newton iterations was somewhat faster than the second one, although the difference was not significant.

The reconstruction corresponding to the complete Algorithm 1 in the top right image of Figure 6 and the one ignoring (only) the mismodeled head shape in the bottom left image are comparable in quality. One can approximately locate the two inhomogeneities in each of these images, although the reconstructions are significantly blurred and unable to capture the full dynamic range of the target conductivity (due to the employed smoothness prior). The reconstruction corresponding to the complete algorithm carries, arguably, less significant artifacts close to the object boundaries, but it also overestimates the size of the imaged head. This is a recurring phenomenon: quite often the complete algorithm is unable to estimate the overall size of the head, even if it is otherwise able to reproduce the head shape. Our hypothesis is that this phenomenon can be explained by the scaling of the object and the average level of the contact conductances having a similar effect on the measurements. Although the reconstruction ignoring all uncertainties in the measurement setup in the bottom right image of Figure 6 can be considered tolerable, it suffers from severe artifacts and is certainly inferior to the other two.

Figure 7 compares the target parameter values $\alpha_{\rm trgt}$, $\zeta_{\rm trgt}$, $\theta_{\rm trgt}$ and $\phi_{\rm trgt}$ with those reconstructed by the complete Algorithm 1. Recall that the algorithm only reconstructs five shape parameters, i.e. $\alpha_{*}\in\mathbb{R}^{5}$, while the target head shape is described by $\alpha_{\rm trgt}\in\mathbb{R}^{10}$. Moreover, the comparison between the reconstructed $\zeta_{*}$ and the corresponding target $\zeta_{\rm trgt}$ is not straightforward as they correspond to different contact shapes in (22); see [16] for further information on this matter. Apart from $\alpha_{1}$ that mainly controls the overall size of the object (cf. Figure 3), the other four estimated shape parameters are reasonable, yet by no means completely accurate. Moreover, the algorithm seems to be able to detect general patterns in the polar and azimuthal angles of the electrodes relative to each other, although many of the electrode positions themselves are not reconstructed accurately. The interplay between the effects the electrode positions, the object shape and the contact conductances have on the measurements is altogether nontrivial, which is definitely reflected in the estimated parameter values illustrated in Figure 7. Be that as it may, the main objective of Algorithm 1 is achieved because the conductivity reconstruction in the top right image of Figure 6 can be considered informative.

5.3.2 Test 2: one inclusion with more uncertain electrode positions

The second test considers the same intended electrode configuration as the first one, but this time around we increased the uncertainty about the electrode positions and the contact conductances by choosing $\varsigma_{\theta}=\varsigma_{\phi}=0.04$, $\tilde{\tau}_{\zeta}=100\,{\rm S}/{\rm m}^{2}$ and $\tilde{\varsigma}_{\zeta}=20\,{\rm S}/{\rm m}^{2}$ when drawing the target values $\theta_{\rm trgt}$, $\phi_{\rm trgt}$ and $\zeta_{\rm trgt}$. The target shape parameters $\alpha_{\rm trgt}$ were redrawn to obtain a new head shape. The to-be-reconstructed conductivity still has a constant background described by $\sigma_{0}=0.2\,{\rm S}/{\rm m}^{2}$ with an embedded, insulating cylindrical inclusion with constant conductivity level $\sigma=0.02\,{\rm S}/{\rm m}$, height $H=7$ cm, and radius $r=1.5$ cm centered at $(1,3,3)$ cm. The central axis of the cylinder is oriented orthogonal to the coronal plane. The prior for the conductivity and the noise model are the same as in the first test.

The target geometry is compared with the one reconstructed by the complete Algorithm 1 in Figure 8; this time the Gauss–Newton iteration converged in eight rounds and about 25 minutes. The shape parameter vector converged to $\alpha_{*}=[2.47,-5.72,1.82,1.76,0.32]^{\rm T}$, whereas the first five components of $\alpha_{\rm trgt}$ are $1.62$, $-5.50$, $0.24$, $2.57$ and $-3.17$. Even though the estimates for the shape parameters are not very accurate, the reconstructed head shape anyway resembles the target head to a certain degree. Even more obviously, the electrode pattern has converged toward the target setup; see the right-hand image of Figure 5 for the geometric initial guess.

Figure 9 presents the target conductivity together with three reconstructions. As in Figure 6 of the first numerical test, the reconstructions correspond to the complete Algorithm 1, to ignoring the mismodeled head shape, and to paying no attention to the inaccuracies in either the electrode positions or the head model. This time around, the reconstruction ignoring all inaccuracies in the measurement setting did not converge in twenty Gauss–Newton iterations, leading to the intolerably bad reconstruction of the conductivity in the bottom right image of Figure 9. The other two reconstructions in Figure 9 contain useful information about the conductivity inside the target head, although their quality is worse than in test 1, most probably due to the initial higher degree of uncertainty about the electrode positions and the contact conductances. The reconstruction produced by the complete algorithm in the top right image of Figure 9 is arguably the best one, separating itself from the one ignoring the inaccuracy in the head model in the bottom left image by exhibiting a more constant background conductivity level. It is also worth noticing the shapes of the cross-sections in the reconstruction corresponding to the complete algorithm clearly mimic those of the target head.

5.3.3 Test 3: two conductive inclusion and larger electrodes

In our final experiment, we consider an alternative measurement configuration with just two electrode belts and a lower number of $M=14+8=22$ electrodes. The electrode radius is also larger than before, namely $R=1\,{\rm cm}$. The probability distributions from which the parameters $\zeta_{\rm trgt}$, $\theta_{\rm trgt}$ and $\phi_{\rm trgt}$ were drawn are the same as in the first experiment, apart from the means $\bar{\theta}$ and $\bar{\phi}$ that are determined by the intended electrode positions and a slightly lower level of uncertainty about the electrode position, i.e. $\varsigma_{\theta}=\varsigma_{\phi}=0.02$. As another novelty, the covariance for the head shapes is multiplied by two: $\alpha_{\rm trgt}\in\mathbb{R}^{10}$ was drawn from $\mathcal{N}(0,2\Gamma_{\alpha})$ where $\Gamma_{\alpha}\in\mathbb{R}^{10\times 10}$ is defined by (18). This corresponds to allowing larger geometric variations in the head shape compared to what the library of [26] predicts. We assume to be aware of this fact and also replace $\Gamma_{\alpha}\in\mathbb{R}^{5\times 5}$ by $2\Gamma_{\alpha}$ in the Tikhonov functional (5.1) and Algorithm 1. The noise model and the prior for the conductivity remain the same as in the preceding two tests. The target conductivity is characterized by two balls with a common constant conductivity level $\sigma_{1}=\sigma_{2}=2\,\,{\rm S}/{\rm m}$, a common radius $r_{1}=r_{2}=1.5$ cm and center points $y_{1}=(-5,0,5)$ cm and $y_{2}=(4,-2,3)$ cm, respectively. The background conductivity level is still $\sigma_{0}\equiv 0.2\,{\rm S}/{\rm m}$.

The target head with the misplaced electrodes and the initial geometric guess for Algorithm 1 are visualized in Figure 10. Figure 11, which is organized in the same way as Figures 6 and 9 above, visualizes the target conductivity and three reconstructions. The complete algorithm converged after eight rounds of Gauss–Newton iteration in about 25 minutes; the other two reconstructions were once again obtained considerably faster. The main conclusions are the same as in the previous two tests: The complete Algorithm 1 and the one ignoring (only) the uncertainties in the head shape produce comparable reconstructions, with the former being somewhat more accurate and also being able to mimic the main characteristics of the shape of the target head. On the other hand, the reconstruction that also ignores the inaccuracies in the electrode positions is contaminated by more severe artifacts.

The geometric target parameters $\alpha_{\rm trgt}$, $\theta_{\rm trgt}$ and $\phi_{\rm trgt}$ are compared with the ones reconstructed by the complete Algorithm 1 in Figure 12. The reconstructed parameters once again exhibit similar patterns as the target values, but the correspondence is not perfect. It is worth mentioning that we also ran the complete algorithm with $\tilde{n}=10$, i.e., with the number of the reconstructed and target shape parameters being the same, but this had essentially no effect on the reconstructions presented in Figures 11 and 12.

6 Conclusion

We have introduced a computational framework for applying absolute EIT to practical head imaging, with the leading idea being to reconstruct the head shape and the electrode positions as a part of a Newton-type output least squares algorithm. Our numerical experiments demonstrated that significant enough variations in the conductivity of an imaged head can be detected by EIT even if the uncertainties about the head shape and electrode positions are on the level that is to be expected in practice. Moreover, it seems that the natural variations in the head shape can be handled to a certain extent by (only) including the estimation of the electrode positions in a reconstruction algorithm implemented in a fixed average head geometry. That is, mismodeling the head shape can be partially compensated by letting the electrodes move freely during the reconstruction process (cf. [15]).

Our investigations did not touch the most severe hindrance for introducing EIT as a head imaging modality, namely the highly resistive skull that makes it difficult to get reliable information about the conductivity in the brain. The most natural way to include the effect of skull in our framework is to employ a head library to form a joint principal component model for the head and skull shapes that are presumably strongly correlated. The joint head and skull geometry could then potentially be parametrized by a low number of shape parameters whose estimation could be included in a Newton-type reconstruction algorithm. Such considerations and testing our algorithm with real-world data are left for future studies.